By browsing through the literature one realizes that the numerical study of gravitational collapse has been neatly divided between two main communities since the early investigations. First, the computational astrophysics community has traditionally focused on the physics of the (type II/Ib/Ic) supernova mechanism, i.e., collapse of unstable stellar iron cores followed by a bounce and explosion. Nowadays, numerical simulations are highly developed, as they include rotation and detailed state-of-the-art microphysics (e.g., EOS and sophisticated treatments for neutrino transport). These studies are currently performed, routinely, in multi-dimensions with advanced HRSC schemes. Progress in this direction has been achieved, however, at the expense of accuracy in the treatment of the hydrodynamics and the gravitational field, by using Newtonian physics. In this approach the emission of gravitational radiation is computed through the quadrupole formula. For reviews of the current status in this direction see [191, 134] and references therein.
On the other hand, the numerical relativity community has been more interested, since the beginning, in highlighting those aspects of the collapse problem having to do with relativistic gravity, in particular, the emission of gravitational radiation from non-spherical collapse (see  for a recent Living Reviews article on gravitational radiation from gravitational collapse). Much of the effort has focused on developing reliable numerical schemes to solve the gravitational field equations and much less emphasis, if any, has been given to the details of the microphysics of the collapse. In addition, this approach has mostly considered gravitational collapse leading to black hole formation, employing relativistic hydrodynamics and gravity. Quite surprisingly, both approaches have barely interacted over the years, except for a handful of simulations in spherical symmetry. Nevertheless, numerical relativity is steadily reaching a state in which it is not unrealistic to expect a much closer interaction in the near future (see, e.g., [262, 259, 89, 260] and references therein).
However, the path to reach such conclusions has been mostly delineated from numerical simulations in one spatial dimension. Fully relativistic simulations of microphysically detailed core collapse off spherical symmetry are still absent and they might well introduce some modifications. The broad picture described above has been demonstrated in multi-dimensions using sophisticated Newtonian models, as is documented in .
Spherically symmetric simulations. May and White's formulation and their corresponding one-dimensional code formed the basis of most spherically symmetric codes built to study the collapse problem. Wilson  investigated the effect of neutrino transport on stellar collapse, concluding that, in contrast to previous results, heat conduction by neutrinos does not produce the ejection of material [60, 14]. The code solved the coupled set of hydrodynamic and Einstein equations, supplemented with a Boltzmann transport equation to describe the neutrino flow. Van Riper  used a spherically symmetric general relativistic code to study adiabatic collapse of stellar cores, considering the purely hydrodynamical bounce as the preferred explosion mechanism. The important role of general relativistic effects to produce collapses otherwise absent in Newtonian simulations was emphasized. Bruenn  developed a code in which the neutrino transport was handled with an approximation, the multigroup flux-limited diffusion. Baron et al.  obtained a successful and very energetic explosion for a model in which the value of the incompressibility modulus of symmetric nuclear matter at zero temperature was particularly small, (its precise value, nowadays preferred around , is still a matter of debate). Mayle et al.  computed neutrino spectra from stellar collapse to stable neutron stars for different collapse models using, as in , multigroup flux-limited diffusion for the neutrino transport. Van Riper  and Bruenn  verified that a softer supranuclear EOS, combined with general relativistic hydrodynamics, increases the chances for a prompt explosion. In [249, 248] and  the neutrino transport was first handled without approximation by using a general relativistic Boltzmann equation. Whereas in [249, 248] the neutrino transport is described by moments of the general relativistic Boltzmann equation in polar slicing coordinates ,  used maximal slicing coordinates. Miralles et al. , employing a realistic EOS based upon a Hartree-Fock approximation for Skyrme-type nuclear forces at densities above nuclear density, found that the explosion was favoured by soft EOS, a result qualitatively similar to that of Baron et al. , who used a phenomenological EOS. Swesty et al.  also focused on the role of the nuclear EOS in stellar collapse on prompt timescales. Contrary to previous results they found that the dynamics of the shock is almost independent of the nuclear incompressibility once the EOS is not unphysically softened as in earlier simulations (e.g., [293, 24, 294, 50, 184]). Swesty and coworkers used a finite temperature compressible liquid drop model EOS for the nuclear matter component . For the shock to propagate promptly to a large radius they found that the EOS must be very soft at densities just above nuclear densities, which is inconsistent with the neutron star mass constraint imposed by observations of the binary pulsar PSR 1913+16.
From the above discussion it is clear that numerical simulations demonstrate a strong sensitivity of the explosion mechanism to the details of the post-bounce evolution: general relativity, the nuclear EOS and the properties of the nascent neutron star, the treatment of the neutrino transport, and the neutrino-matter interaction. Recently, state-of-the-art treatments of the neutrino transport have been achieved, even in the general relativistic case, by solving the Boltzmann equation in connection with the hydrodynamics, instead of using multigroup flux-limited diffusion [177, 232, 233, 154]. The results of the few spherically symmetric simulations available show, however, that the models do not explode, neither in the Newtonian or in the general relativistic case. Therefore, computationally expensive multi-dimensional simulations with Boltzmann neutrino transport, able to account for convective effects, are needed to draw further conclusions on the viability of the neutrino-driven explosion mechanism.
From the numerical point of view, many of the above investigations used artificial viscosity techniques to handle shock waves. Together with a detailed description of the microphysics, the correct numerical modeling of the shock wave is the major issue in stellar collapse. In this context, the use of HRSC schemes, specifically designed to capture discontinuities, is much more recent, starting in the late 1980s with the Newtonian simulations of Fryxell et al. , which used an Eulerian second-order PPM scheme (see  for a review of the present status). There are only a few relativistic simulations, so far restricted to spherical symmetry [129, 244, 211].
Martí et al.  implemented Glaister's approximate Riemann solver  in a Lagrangian Newtonian hydrodynamical code. They performed a comparison of the energetics of a stellar collapse simulation using this HRSC scheme and a standard May and White code with artificial viscosity, for the same initial model, grid size and EOS. They found that the simulation employing a Godunov-type scheme produced larger kinetic energies than that corresponding to the artificial viscosity scheme, with a factor of two difference in the maximum of the infalling velocity. Motivated by this important difference, Janka et al.  repeated this computation with a different EOS, using a PPM second-order Godunov-type scheme, disagreeing with Martí et al. . The state-of-the-art implementation of the tensorial artificial viscosity terms in , together with the very fine numerical grids employed (unrealistic for three-dimensional simulations), could be the reason for the discrepancies.
As mentioned in Section 2.1.3, Godunov-type methods were first used to solve the equations of general relativistic hydrodynamics in , where the characteristic fields of the one-dimensional (spherically symmetric) system were derived. The first astrophysical application was stellar collapse. In  the hydrodynamic equations were coupled to the Einstein equations in spherical symmetry. The field equations were formulated as a first-order flux-conservative hyperbolic system for a harmonic gauge , somehow ``resembling'' the hydrodynamic equations. HRSC schemes were applied to both systems simultaneously (only coupled through the source terms of the equations). Results for simple models of adiabatic collapse can be found in [161, 38, 126].
A comprehensive study of adiabatic, one-dimensional core collapse using explicit upwind HRSC schemes was presented in  (see  for a similar computation using implicit schemes). In this investigation the equations for the hydrodynamics and the geometry are written using radial gauge polar slicing (Schwarzschild-type) coordinates. The collapse is modeled with an ideal gas EOS with a non-constant adiabatic index, which depends on the density as , where is the bounce density, and if and otherwise . A set of animations of the simulations presented in  is included in the four movies in Figure 4 . They correspond to the rather stiff Model B of : , , and . The initial model is a white dwarf having a gravitational mass of . The animations in Figure 4 show the time evolution of the radial profiles of the following fields: velocity (top-left), logarithm of the rest-mass density (top-right), gravitational mass (bottom-left), and the square of the lapse function (bottom-right).
The movies show that the shock is sharply resolved and free of spurious oscillations. The radius of the inner core at the time of maximum compression, at which the infall velocity is maximum (v =-0.62 c), is . The gravitational mass of the inner core at the time of maximum compression is . The minimum value that the quantity reaches is 0.14, which indicates the highly relativistic character of these simulations (at the surface of a typical neutron star the value of the lapse function squared is ).
Axisymmetric Newtonian simulations. Beyond spherical symmetry most of the existing simulations of core collapse supernova are Newtonian. Axisymmetric investigations have been performed by [190, 37, 42] using a realistic EOS and some treatment of weak interaction processes, but neglecting neutrino transport, and by [135, 188, 132, 101, 102] employing some approximate description of neutrino transport. In addition, [84, 310, 319] have performed Newtonian parameter studies of the axisymmetric collapse of rotating polytropes.
Among the more detailed multi-dimensional non-relativistic hydrodynamical simulations are those performed by the MPA/Garching group (an on-line sample can be found at their website ). As mentioned earlier, these include advanced microphysics and employ accurate HRSC integration schemes. To illustrate the degree of sophistication achieved in Newtonian simulations, we present in the movie in Figure 5 an animation of the early evolution of a core collapse supernova explosion up to after core bounce and shock formation (Figure 5 depicts just one intermediate snapshot at ). The movie shows the evolution of the entropy within the innermost of the star.
The initial data used in these calculations is taken from the pre-collapse model of a blue supergiant star in . The computations begin after core bounce from a one-dimensional model of . This model is obtained with the one-dimensional general relativistic code mentioned previously , which includes a detailed treatment of the neutrino physics and a realistic EOS, and which is able to follow core collapse, bounce, and the associated formation of the supernova shock. Because of neutrino cooling and energy losses due to the dissociation of iron nuclei, the shock initially stalls inside the iron core.
The movie shows how the stalling shock is ``revived'' by neutrinos streaming from the outer layers of the hot, nascent neutron star in the center. A negative entropy gradient builds up between the so called ``gain-radius'', which is the position where cooling of hot gas by neutrino emission switches into net energy gain by neutrino absorption, and the shock further out. Convective instabilities therefore set in, which are characterized by large rising blobs of neutrino-heated matter and cool, narrow downflows. The convective energy transport increases the efficiency of energy deposition in the post-shock region by transporting heated material near the shock and cooler matter near the gain radius where the heating is strongest. The freshly heated matter then rises again while the shock is distorted by the upward streaming bubbles. The reader is addressed to  for a more detailed description of a more energetic initial model.
Axisymmetric relativistic simulations. Previous investigations in general relativistic rotational core collapse were mainly concerned with the question of black hole formation under idealized conditions (see Section 4.1.2), but none of those studies has really addressed the problem of supernova core collapse, which proceeds from white dwarf densities to neutron star densities and involves core bounce, shock formation, and shock propagation.
Wilson  first computed neutron star bounces of polytropes, and measured the gravitational wave emission. The initial configurations were either prolate or slightly aspherical due to numerical errors of an otherwise spherical collapse.
More than twenty years later, and with no other simulations in between, the first comprehensive numerical study of relativistic rotational supernova core collapse in axisymmetry has been performed by Dimmelmeier et al. [64, 65, 66, 67], who computed the gravitational radiation emitted in such events. The Einstein equations were formulated using the so-called conformally flat metric approximation . Correspondingly, the hydrodynamic equations were cast as the first-order flux-conservative hyperbolic system described in Section 2.1.3 . Details of the methodology and of the numerical code are given in .
Dimmelmeier et al.  have simulated the evolution of 26 models in both Newtonian and relativistic gravity. The initial configurations are differentially rotating relativistic 4/3-polytropes in equilibrium, which have a central density of . Collapse is initiated by decreasing the adiabatic index to some prescribed fixed value with . The EOS consists of a polytropic and a thermal part for a more realistic treatment of shock waves. Any microphysics, such as electron capture and neutrino transport, is neglected. The simulations show that the three different types of rotational supernova core collapse (and their associated gravitational waveforms) identified in previous Newtonian simulations by Zwerger and Müller  - regular collapse, multiple bounce collapse, and rapid collapse - are also present in relativistic gravity. However, rotational core collapse with multiple bounces is only possible in a much narrower parameter range in general relativity. Relativistic gravity has, furthermore, a qualitative impact on the dynamics: If the density increase due to the deeper relativistic potential is sufficiently large, a collapse that is stopped by centrifugal forces at subnuclear densities (and thus undergoes multiple bounces) in a Newtonian simulation becomes a regular, single bounce collapse in relativistic gravity. Such collapse type transitions have important consequences for the maximum gravitational wave signal amplitude. Moreover, in several of the relativistic models discussed in , the rotation rate of the compact remnant exceeds the critical value at which Maclaurin spheroids become secularly, and in some cases even dynamically, unstable against triaxial perturbations.
The movie presented in Figure 6 shows the time evolution of a multiple bounce model (model A2B4G1 in the notation of ), with a type II gravitational wave signal. The left panel shows isocontours of the logarithm of the density together with the corresponding meridional velocity field distribution. The two panels on the right show the time evolution of the gravitational wave signal (top panel) and of the central rest-mass density. In the animation the ``camera'' follows the multiple bounces by zooming in and out. The panels on the right show how each burst of gravitational radiation coincides with the time of bounce of the collapsing core. The oscillations of the nascent proto-neutron star are therefore imprinted on the gravitational waveform. Additional animations of the simulations performed by  can be found at the MPA's website .
The relativistic models analyzed by Dimmelmeier et al.  cover almost the same range of gravitational wave amplitudes ( for a source at a distance of ) and frequencies () as the corresponding Newtonian ones . Averaged over all models, the total energy radiated in the form of gravitational waves is in the relativistic case, and in the Newtonian case. For all collapse models that are of the same type in both Newtonian and relativistic gravity, the gravitational wave signal in relativistic gravity is of lower amplitude. If the collapse type changes, either weaker or stronger signals are found in the relativistic case. Therefore, the prospects for detection of gravitational wave signals from axisymmetric supernova rotational core collapse do not improve when taking into account relativistic gravity. The signals are within the sensitivity range of the first generation laser interferometer detectors if the source is located within the Local Group. An on-line waveform catalogue for all models analyzed by Dimmelmeier et al.  is available at the MPA's website . Finally, we note that a program to extend these simulations to full general relativity (relaxing the conformally flat metric assumption) has been very recently started by Shibata .
Three-dimensional simulations. To date, there are no three-dimensional relativistic simulations of gravitational collapse in the context of supernova core collapse yet available. All the existing computations have employed Newtonian physics. This situation, however, might change in the near future, as very recently the first fully relativistic three-dimensional simulations of gravitational collapse leading to black hole formation have been accomplished, for rapidly rotating, supermassive neutron stars  and supermassive stars  (see Section 4.1.2), for the head-on collision of two neutron stars , and for the coalescence of neutron star binaries [258, 266, 267] (see Section 4.3).
Concerning Newtonian studies, Bonazzola and Marck  performed pioneering three-dimensional simulations, using pseudo-spectral methods that are very accurate and free of numerical or intrinsic viscosity. They confirmed the low amount of energy radiated in gravitational waves regardless of the initial conditions of the collapse: axisymmetric, rotating or tidally deformed (see also ). This result only applies to the pre-bounce phase of the supernova collapse as the simulations did not consider shock propagation after bounce. More recently, [234, 48] have extended the work of  studying, with two different PPM hydrodynamics codes, the dynamical growth of non-axisymmetric (bar mode) instabilities appearing in rotational post-collapse cores. A relativistic extension has been performed by  (see Section 4.3.1).
Spherically symmetric simulations. Van Riper , using a (Lagrangian) May and White code, analyzed the mass division between the formation of a neutron star or a black hole after gravitational collapse. For the typical (cold) EOS used, the critical state was found to lie between the collapses of and cores.
In , a one-dimensional code based on Wilson's hydrodynamical formulation was used to simulate a general relativistic (adiabatic) collapse to a black hole. The fluid equations were finite-differenced in a mixed Eulerian-Lagrangian form with the introduction of an arbitrary ``grid velocity'' to ensure sufficient resolution throughout the entire collapse. The Einstein equations were formulated following the ADM equations. Isotropic coordinates and a maximal time-slicing condition were used to avoid the physical singularity once the black hole forms. Due to the non-dynamical character of the gravitational field in the case of spherical symmetry (i.e., the metric variables can be computed at every time step solving elliptic equations), the code developed by  could follow relativistic configurations for many collapse timescales after the initial appearance of an event horizon.
A Lagrangian scheme based on the Misner and Sharp  formulation for spherically symmetric gravitational collapse (as in ) was developed by Miller and Motta  and later by Baumgarte et al. . The novelty of these codes was the use of an outgoing null coordinate u (an ``observer-time'' coordinate, as suggested previously by ) instead of the usual Schwarzschild time t appearing in the Misner and Sharp equations. Outgoing null coordinates are ideally suited to study black hole formation as they never cross the black hole event horizon. In these codes, the Hernández-Misner equations  (or, alternatively, the Misner-Sharp equations) were solved by an explicit finite difference scheme similar to the one used by . In , the collapse of an unstable polytrope to a black hole was first achieved using null slicing. In , the collapse of a polytrope with was compared to the result of , using the Misner-Sharp equations and finding a 10% agreement. This work showed numerically that the use of a retarded time coordinate allows for stable evolutions after the black hole has formed. The Lagrangian character of both codes has prevented their multi-dimensional extension.
Linke et al.  analyzed the gravitational collapse of supermassive stars in the range . In the same spirit as in Refs. [181, 26], the coupled system of Einstein and fluid equations was solved employing coordinates adapted to a foliation of the spacetime with outgoing null hypersurfaces. From the computed neutrino luminosities, estimates of the energy deposition by -annihilation were obtained. Only a small fraction of this energy is deposited near the surface of the star, where it could cause the ultrarelativistic flow believed to be responsible for GRBs. However, the simulations show that for collapsing supermassive stars with masses larger than , the energy deposition is at least two orders of magnitude too small to explain the energetics of observed long-duration bursts at cosmological redshifts.
The interaction of massless scalar fields with self-gravitating relativistic stars was analyzed in  by means of fully dynamic numerical simulations of the Einstein-Klein-Gordon perfect fluid system. A sequence of stable, self-gravitating K =100, N =1 relativistic polytropes, increasing the central density from to () was built. Using a compactified spacetime foliation with outgoing null cones, Siebel et al.  studied the fate of the relativistic stars when they are hit by a strong scalar field pulse with a mass of the order of , finding that the star is either forced to oscillate in its radial modes of pulsation or to collapse to a black hole on a dynamical timescale. The radiative signals, read off at future null infinity, consist of several quasi-normal oscillations and a late time power-law tail, in agreement with the results predicted by (linear) perturbation analysis of wave propagation in an exterior Schwarzschild geometry.
Axisymmetric simulations. Beyond spherical symmetry the investigations of black hole formation started with the work of Nakamura , who first simulated general relativistic rotating stellar collapse. He adopted the (2+1)+1 formulation of the Einstein equations in cylindrical coordinates and introduced regularity conditions to avoid divergence behavior at coordinate singularities (the plane z =0) . The equations were finite-differenced using the donor cell scheme plus Friedrichs-Lax type viscosity terms. Nakamura used a ``hypergeometric'' slicing condition (which reduces to maximal slicing in spherical symmetry), which prevents the grid points from hitting the singularity if a black hole forms. The simulations could track the evolution of the collapse of a ``core'' of a massive star with different amounts of rotational energy and an initial central density of , up to the formation of a rotating black hole. However, the resolution affordable due to limitations in computer resources ( grid points) was not high enough to compute the emitted gravitational radiation. Note that the energy emitted in gravitational waves is very small compared to the total rest mass energy, which makes its accurate numerical computation very challenging. In subsequent works, Nakamura  (see also ) considered a configuration consisting of a neutron star (, ) with an accreted envelope of , which was thought to mimic mass fall-back in a supernova explosion. Rotation and infall velocity were added to such configurations, simulating the evolution depending on the prescribed rotation rates and rotation laws.
Bardeen, Stark, and Piran, in a series of papers [22, 276, 228, 277], studied the collapse of rotating relativistic () polytropic stars to black holes, succeeding in computing the associated gravitational radiation. The field and hydrodynamic equations were formulated using the 3+1 formalism, with radial gauge and a mixture of (singularity avoiding) polar and maximal slicing. The initial model was a spherically symmetric relativistic polytrope in equilibrium of mass M, central density , and radius . Rotational collapse was induced by lowering the pressure in the initial model by a prescribed fraction, and by simultaneously adding an angular momentum distribution that approximates rigid-body rotation. Their parameter space survey showed how black hole formation (of the Kerr class) occurs only for angular momenta less than a critical value. The numerical results also demonstrated that the gravitational wave emission from axisymmetric rotating collapse to a black hole was , and that the general waveform shape was nearly independent of the details of the collapse.
Evans , building on previous work by Dykema , also studied the axisymmetric gravitational collapse problem for non-rotating matter configurations. His numerical scheme to integrate the matter fields was more sophisticated than in previous approaches, as it included monotonic upwind reconstruction procedures and flux limiters. Discontinuities appearing in the flow solution were stabilized by adding artificial viscosity terms in the equations, following Wilson's approach.
Most of the axisymmetric codes discussed so far adopted spherical polar coordinates. Numerical instabilities are encountered at the origin (r =0) and at the polar axis (), where some fields diverge due to the coordinate singularities. Evans made important contributions toward regularizing the gravitational field equations in such situations . These coordinate problems have deterred researchers from building three-dimensional codes in spherical coordinates. In this case, Cartesian coordinates are adopted which, despite the desired property of being everywhere regular, present the important drawback of not being adapted to the topology of astrophysical systems. As a result, this has important consequences on the grid resolution requirements. The only extension of an axisymmetric 2+1 code in spherical coordinates to three dimensions has been accomplished by Stark , although no applications in relativistic astrophysics have yet been presented.
Recently, Alcubierre et al.  proposed a method (``cartoon'') that does not suffer from stability problems at coordinate singularities and in which, in essence, Cartesian coordinates are used even for axisymmetric systems. Using this method, Shibata  investigated the effects of rotation on the criterion for prompt adiabatic collapse of rigidly and differentially rotating () polytropes to a black hole. Collapse of the initial approximate equilibrium models (computed by assuming a conformally flat spatial metric) was induced by a pressure reduction. In  it was shown that the criterion for black hole formation depends strongly on the amount of angular momentum, but only weakly on its (initial) distribution. Shibata also studied the effects of shock heating using a gamma-law EOS, concluding that it is important in preventing prompt collapse to black holes in the case of large rotation rates. Using the same numerical approach, Shibata and Shapiro  have recently studied the collapse of a uniformly rotating supermassive star in general relativity. The simulations show that the star, initially rotating at the mass-shedding limit, collapses to form a supermassive Kerr black hole with a spin parameter of . Roughly 90% of the mass of the system is contained in the final black hole, while the remaining matter forms a disk orbiting around the hole.
Alternatively, existing axisymmetric codes employing the characteristic formulation of the Einstein equations , such as the (vacuum) code presented in , do not share the axis instability problems of most 2+1 codes. Hence, axisymmetric characteristic codes, once conveniently coupled to hydrodynamics codes, are a promising way to study the axisymmetric collapse problem. First steps in this direction are reported in , where an axisymmetric Einstein-perfect fluid code is presented. This code achieves global energy conservation for a strongly perturbed neutron star spacetime, for which the total energy radiated away by gravitational waves corresponds to a significant fraction of the Bondi mass.
Three-dimensional simulations. Hydrodynamical simulations of the collapse of supermassive () uniformly rotating neutron stars to rotating black holes, using the code discussed in Section 3.3.1, are presented in . The simulations show no evidence for massive disk formation or outflows, which can be related to the moderate initial compactness of the stellar models (). A proof-of-principle of the capabilities of the code discussed in Section 3.3.2 to study black hole formation is presented in , where the gravitational collapse of a spherical, unstable, relativistic polytrope is discussed. Similar tests with differentially rotating polytropes are given in  for a recent artificial viscosity-based, three-dimensional general relativistic hydrodynamics code.
Evans and Coleman  first observed critical phenomena in the spherical collapse of radiation fluid, i.e., a fluid obeying an EOS with and . The threshold of black hole formation was found for a critical exponent , in close agreement with that obtained in scalar field collapse . Their study used Schwarzschild coordinates in radial gauge and polar slicing, and the hydrodynamic equations followed Wilson's formulation. Subsequently, Maison , using a linear stability analysis of the critical solution, showed that the critical exponent varies strongly with the parameter of the EOS. More recently, Neilsen and Choptuik , using a conservative form of the hydrodynamic equations and HRSC schemes, revisited this problem. The use of a conservative formulation and numerical schemes well adapted to describe ultrarelativistic flows  made it possible to find evidence of the existence of critical solutions for large values of the adiabatic index (), extending the previous upper limit.
|Numerical Hydrodynamics in General Relativity
José A. Font
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to firstname.lastname@example.org