By browsing through the literature one immediately realizes that the numerical study of gravitational collapse has been neatly split in two main directions since the early investigations. First, the numerical astrophysical community has traditionally focused on the physics of the (type II) Supernova mechanism, i.e., collapse of unstable 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, using 3D Cartesian grids with advanced HRSC schemes. Progress in this direction has been achieved, however, at the expense of necessary approximations 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 a recent review of the current status in this direction see  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 a non-spherical collapse. Much of the effort has focused on developing reliable numerical schemes to evolve the gravitational field and much less emphasis, if any, has been given to the details of the microphysics of core collapse. This approach has also 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. We turn next to their description.
The path to reach such conclusions has been delineated from numerical simulations in one spatial dimension. Fully relativistic simulations of microphysically detailed core collapse off spherical symmetry are still absent and they could introduce some modifications. The overall picture described above has been demonstrated in multidimensions using sophisticated Newtonian models, as is documented extensively in .
May and White's formulation and their corresponding one-dimensional code formed the basis of the majority of subsequent spherically symmetric codes to study the collapse problem. Wilson  investigated the effect of neutrino transport on the stellar collapse concluding that heat conduction by neutrinos does not produce the ejection of material, in contrast to previous results [54, 13]. 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 with a particularly small value of the incompressibility modulus of symmetric nuclear matter at zero temperature, (whose 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 , 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 [189, 188] and  the neutrino transport was first handled without approximation by using a general relativistic Boltzmann equation. Whereas in [189, 188] 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 for densities above nuclear density, obtained results qualitatively similar to those of Baron et al. , who used a phenomenological EOS, i.e., the explosion was favoured by soft 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 are almost independent of the nuclear incompressibility, once the EOS is not unphysically softened as in earlier simulations (e.g. [221, 20, 222, 44, 141]). Swesty and coworkers used a finite temperature compressible liquid drop EOS . 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.
All the aforementioned investigations used artificial viscosity techniques to handle shock waves. Together with a detailed description of the microphysics, the correct numerical modeling of the shock is the major issue in stellar collapse. In this context, the use of HRSC schemes, specifically designed to capture discontinuities, is much more recent. It started in the late eighties with the Newtonian simulations performed by Fryxell et al.  using 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. Most of them have been performed by the Valencia group, first employing Newtonian gravity , and later relativistic gravity [101, 33, 98, 184]. 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 of the discrepancies.
As mentioned in Section 2.2, 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 could be applied, hence, to both systems simultaneously (only coupled through the source terms of the equations). Results for simple models of adiabatic collapse following this approach can be found in [123, 33, 98].
A comprehensive study of adiabatic one-dimensional core collapse using explicit HRSC schemes was presented in  (a similar computation using implicit schemes was later performed by ). 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 varies with 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 MPEG movies in Fig. 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 show the time evolution of the radial profiles of the following fields: velocity (movie 1), logarithm of the rest-mass density (movie 2), gravitational mass (movie 3) and the square of the lapse function (movie 4).
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 the quantity , the lapse function squared, reaches is 0.14, which indicates the highly relativistic character of these simulations (at the surface of a typical neutron star ).
Studies of black hole formation in strong field hydrodynamical simulations of gravitational collapse are less frequent [221, 196, 24, 156]. In , using the (Lagrangian) May and White code described above, the mass division between neutron star and black hole formation was considered by means of general relativistic simulations. For the typical (cold) EOS, the critical state was found to lie between the collapses of 1.95 and 1.96 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 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) this code was able to follow relativistic configurations for many collapse time scales 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, Shapiro and Teukolsky . 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. These 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 . In  the collapse of an unstable polytrope, followed to black hole formation, was first achieved using null slicing. In  the collapse of a 1.4 polytrope with to a black hole was compared to , using the Misner-Sharp equations, finding a 10% agreement. However, the evolution stopped after the formation of the apparent horizon. This work demonstrated numerically that the use of a retarded time coordinate allowed for much longer evolutions after black hole formation. The Lagrangian character of both codes has prevented their multidimensional extension.
Recently, Neilsen and Choptuik  have studied critical collapse of a spherically symmetric perfect fluid with the EOS using a conservative form of the hydrodynamic equations and HRSC schemes. Critical phenomena in gravitational collapse were first discovered numerically in simulations of the massless Klein-Gordon field minimally coupled to gravity . Since then critical phenomena arising at the threshold of black hole formation, in parametrized families of collapse, have been found in a variety of physical systems, including the perfect fluid model. The Einstein equations coupled to the hydrodynamic equations were integrated in  in spherical symmetry adopting the ADM formalism. The use of a conservative formulation and numerical schemes well adapted to describe ultrarelativistic flows  made it possible to find evidence for the existence of critical solutions for large values of the adiabatic index (), extending the previous upper limit.
One-dimensional numerical simulations of self-gravitating matter accreting onto a black hole, formed after the gravitational collapse of a massive star, have been performed in . Using the formulation outlined in Section 2.2.2 and HRSC schemes, Papadopoulos and Font  performed the simulations adopting an ingoing null foliation of a spherically symmetric black hole spacetime. Such a foliation has the interesting property of penetrating the black hole horizon, allowing for an unambiguous numerical treatment of the inner boundary. The simulations reported in  describe the capture of an initial Gaussian shell distribution of perfect fluid matter by a non-rotating black hole. The results show that the accretion process initiates, as expected, a rapid increase of the mass of the apparent horizon. This is depicted in Fig. 5 . The horizon almost doubles its size during the first 10 M -15 M (this is enlarged in the insert of Fig. 5). Once the main accretion process has finished, the mass of the horizon slowly increases in a quasi-steady manner, whose rate depends on the mass accretion rate imposed at the world-tube of the integration domain.
All numerical studies of axisymmetric gravitational collapse in general relativity, which will be briefly discussed next, have used Wilson's formulation of the hydrodynamic equations and finite difference schemes with artificial viscosity. The problem has not yet been considered in two (and three) dimensions using conservative formulations and HRSC schemes.
The investigations of the collapse problem started with the work of Nakamura , who was the first to calculate a 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 collapsing (perfect fluid) matter and the formation of a black hole. However, the scheme was not accurate enough to compute the emitted gravitational radiation. The reason is that the energy emitted in this process as gravitational radiation is very small compared to the total rest mass energy, making its accurate numerical computation very challenging.
In a series of papers [18, 211, 177, 212], Bardeen, Stark and Piran revisited this problem studying the collapse of rotating relativistic polytropic stars to black holes, succeeding in computing the gravitational radiation emission. 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. Their parameter space survey showed how black hole formation (of the Kerr class) occurs only for angular momenta less than a critical value. Their numerical results also demonstrated that the gravitational wave emission from axisymmetric rotating collapse to a black hole was very low, , and that the general waveform shape was nearly independent of the details of the collapse. Concerning this we note that, in a previous investigation, Müller  obtained the first numerical evidence of the low gravitational efficiency of the core-collapse scenario. By means of axisymmetric Newtonian simulations he found that was radiated as gravitational waves.
Evans , building on previous work by Dykema , also studied the gravitational collapse problem for non-rotating matter configurations. His numerical scheme to integrate the matter fields was more sophisticated than in the previous approaches as it included monotonic upwind reconstruction procedures and flux limiters. Discontinuities appearing in the flow solutions were stabilized by adding artificial viscosity terms in the equations, following Wilson's approach. Evans  was able to show that Newtonian gravity and the quadrupole formula for gravitational radiation were inadequate to study the problem.
Most of the aforementioned axisymmetric codes adopted spherical polar coordinates. Numerical instabilities are encountered at the origin (r =0) and, mainly, at the polar axis () where some fields diverge due to the coordinate singularities. Evans did important contributions towards 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 3+1 code using spherical coordinates to three dimensions has been accomplished by Stark , although no applications in relativistic astrophysics have yet been presented.
Alternatively, existing axisymmetric codes employing the characteristic formulation of the Einstein equations , such as the (vacuum) one presented in , do not share the axis instability problems of 3+1 codes. Hence, axisymmetric characteristic codes, once conveniently coupled to matter fields, could be a promising way of studying the axisymmetric collapse problem.
In a much simpler approach some investigations have considered dynamical evolutions of a collisionless gas of particles in full general relativity [197, 198, 199, 200, 4]. Such pressureless particles move along geodesics of the spacetime. The geodesic equations are, in comparison to the hydrodynamic equations, trivial to integrate. All efforts can focus then on solving the gravitational field equations. Shapiro and Teukolsky [197, 198] built a code to solve the Vlasov equation in general relativity by N -body particle simulation, computing the gravitational field using the ADM formalism. They studied the origin of quasars and active galactic nuclei via the collapse of relativistic star clusters to supermassive black holes. They also considered the collapse of prolate spheroids finding that if those are sufficiently compact, the final singularities are hidden inside black holes. On the controversial side, however, they found that when the spheroids are large enough there are no apparent horizons, in disagreement with the cosmic censorship hypothesis , because a naked singularity may form in non-spherical relativistic collapse . Collapses with small enough angular momentum also showed the same behavior . However, such computations have not yet been attempted using a fully hydrodynamical code, with the inclusion of pressure effects. Hence, the possibility of naked singularity formation in more realistic numerical simulations still remains to be analyzed.
Semi-analytical studies of finite-sized collections of dust, shaped in the form of stars or shells, falling isotropically onto a black hole are available in the literature [151, 91, 201, 163, 176]. These investigations approximate gravitational collapse by a dust shell of mass m falling into a Schwarzschild black hole of mass . These studies have shown that for a fixed amount of infalling mass the gravitational radiation efficiency is reduced compared to the point particle limit, never exceeding that of a particle with the same mass, the reason being cancellations of the emission from distinct parts of the extended object.
In  such conclusions were corroborated with a numerical simulation of the gravitational radiation emitted during the accretion process of an extended object onto a black hole. The first order deviations from the exact black hole geometry were approximated using curvature perturbations induced by matter sources whose non-linear evolution was integrated using a (non-linear) hydrodynamical code (adopting the Valencia conservative formulation and HRSC schemes). All possible types of curvature perturbations are captured in the real and imaginary parts of the Weyl tensor scalar (see, e.g., ). In the framework of the Newman-Penrose formalism the equations for the perturbed Weyl tensor components decouple, and when written in the frequency domain even separate . Papadopoulos and Font  used the limiting case for Schwarzschild black holes, i.e. the inhomogeneous Bardeen-Press equation . The simulations showed the gradual excitation of the black hole quasi-normal mode frequency by sufficiently compact shells.
An example of these simulations appears in the MPEG movie of Fig. 6 . This movie shows the time evolution of the shell density (left panel) and the associated gravitational waveform during a complete accretion/collapse event. The (quadrupolar) shell is parametrized according to with and . Additionally denotes a logarithmic radial (Schwarzschild) coordinate. The animation shows the gradual collapse of the shell onto the black hole. This process triggers the emission of gravitational radiation. In the movie it can be clearly seen how the burst of the emission coincides with the more dynamical accretion phase, when the shell crosses the peak of the potential and is subsequently captured by the hole. The gravitational wave signal coincides with the quasinormal ringing frequency of the black hole, 17 M . The existence of an initial burst, causally separated from the physical burst, is also noticeable in the movie. It just reflects the gravitational radiation content of the initial data (see  for a detailed explanation).
We end this section by pointing out that, 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 been done employing Newtonian physics. We point out, however, that this situation may 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 performed for rapidly-rotating supramassive neutron stars , for the head-on collision of two neutron stars , or for the coalescence of neutron star binaries [202, 205] (see Section 4.3).
Bonazzola and Marck  performed the first three-dimensional simulations, using pseudo-spectral methods, 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. Their results are only applicable to the pre-bounce phase of the supernova collapse as the simulations do not consider shock propagation after bounce.
The more detailed multi-dimensional non-relativistic simulations are those performed by the MPA/Garching group (an on-line sample can be found at their website ). As mentioned before these include sophisticated microphysics and employ accurate HRSC integration schemes. To illustrate the degree of sophistication achieved in Newtonian simulations we present in the MPEG movie in Fig. 7 an animation of the early evolution of a core collapse supernova explosion up to 220 ms after core bounce and shock formation (only an intermediate snapshot at 130 ms is depicted in Fig. 7). The movie shows the evolution of the entropy within the innermost 3000 km of the star.
The initial data used in these calculations are taken from the pre-collapse model of a blue supergiant star of . The computations start 20 ms 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.
|Numerical Hydrodynamics in General Relativity
José A. Font
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to email@example.com