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 [146] 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 [146].

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 [225] 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 [221] 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 [43] developed a code in which the neutrino transport was handled with an approximation, the multigroup flux-limited diffusion. Baron et al. [20] 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 [83], is still a matter of debate). Mayle et al. [134] computed neutrino spectra from stellar collapse to stable neutron stars for different collapse models using, as [43], multigroup flux-limited diffusion for the neutrino transport. Van Riper [222] and Bruenn [44] verified that a softer supranuclear EOS, combined with general relativistic hydrodynamics, increases the chances for a prompt explosion. In [189, 188] and [136] 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 [18], [136] used maximal slicing coordinates. Miralles et al. [141], 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. [20], who used a phenomenological EOS, i.e., the explosion was favoured by soft EOS. Swesty et al. [214] 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 [114]. 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. [79] using an Eulerian second order PPM scheme (see [146] 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 [124], and later relativistic gravity [101, 33, 98, 184]. Martí et al. [124] implemented Glaister's approximate Riemann solver [82] 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. [105] repeated this computation with a different EOS using a PPM second order Godunov-type scheme, disagreeing with Martí et al. [124]. The state-of-the-art implementation of the tensorial artificial viscosity terms in [105] 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 [125], where the characteristic fields of the one-dimensional (spherically symmetric) system were derived. The first astrophysical application was stellar collapse. In [33] 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 [34], 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 [184] (a similar computation using implicit schemes was later performed by [236]). 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 [221]. A set of animations of the simulations presented in [184] is included in the four MPEG movies in Fig. 4 . They correspond to the rather stiff Model B of [184]: 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 [221], 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 [196] 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 [142] formulation for spherically symmetric gravitational collapse
(as in [221]) was developed by Miller and Motta [139] and later by Baumgarte, Shapiro and Teukolsky [24]. The novelty of these codes was the use of an outgoing null
coordinate
*u*
(an ``observer-time'' coordinate, as suggested previously
by [96]) 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 [96] (or, alternatively, the Misner-Sharp equations) were solved by
an explicit finite difference scheme similar to [221]. In [139] the collapse of an unstable polytrope, followed to black hole
formation, was first achieved using null slicing. In [24] the collapse of a 1.4
polytrope with
to a black hole was compared to [221], 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 [156] 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 [50]. 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 [156] in spherical symmetry adopting the ADM formalism. The use of a
conservative formulation and numerical schemes well adapted to
describe ultrarelativistic flows [157] 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 [171]. Using the formulation outlined in Section
2.2.2
and HRSC schemes, Papadopoulos and Font [171] 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 [171] 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 [148], 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) [149]. 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 [145] 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 [65], building on previous work by Dykema [59], 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 [65] 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 [65]. 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 [210], although no applications in relativistic astrophysics have yet
been presented.

Alternatively, existing axisymmetric codes employing the characteristic formulation of the Einstein equations [232], such as the (vacuum) one presented in [85], 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 [174], because a naked singularity may form in non-spherical
relativistic collapse [199]. Collapses with small enough angular momentum also showed the
same behavior [200]. 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 [170] 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., [49]). 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 [216]. Papadopoulos and Font [170] used the limiting case for Schwarzschild black holes, i.e. the inhomogeneous Bardeen-Press equation [19]. 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 [170] 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 [203], for the head-on collision of two neutron stars [140], or for the coalescence of neutron star binaries [202, 205] (see Section 4.3).

Bonazzola and Marck [37] 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 [3]). 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 [235]. The computations start 20 ms after core bounce from a one-dimensional model of [45]. This model is obtained with the one-dimensional general relativistic code mentioned previously [43] 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 [107] for a more detailed description of a more energetic initial model.

Numerical Hydrodynamics in General Relativity
José A. Font
http://www.livingreviews.org/lrr-2000-2
© Max-Planck-Gesellschaft. ISSN 1433-8351 Problems/Comments to livrev@aei-potsdam.mpg.de |