4.2 Accretion onto black holes4 Simulations of Astrophysical Phenomena4 Simulations of Astrophysical Phenomena

4.1 Gravitational collapse 

The study of gravitational collapse of massive stars has been largely pursued numerically over the years. This problem was already the main motivation when May and White built the first one-dimensional code [132, 133]. It is worth noticing that this code was conceived to integrate the coupled system of Einstein field equations and relativistic hydrodynamics, sidestepping Newtonian approaches.

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 [146Jump To The Next Citation Point In The Article] 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.

4.1.1 One dimensional simulations 

The standard model of type II Supernovae involves the presence of a strong shock wave formed at the edge between the homologous inner core and the outer core which falls at roughly free-fall speed. A schematic illustration of the dynamics of this system is depicted in Fig.  3 . The shock is produced after the bounce of the inner core when its density exceeds nuclear density. Numerical simulations have tried to elucidate if this shock would be strong enough to propagate outwards through the star's mantle and envelope given certain initial conditions for the material in the core, an issue subject to important uncertainties in the nuclear EOS, as well as in the outer layers of the star. In the accepted scenario which has emerged from the numerical simulations, the following features have been found to be necessary for all plausible explosion mechanisms: the existence of the shock wave together with neutrino re-heating that re-energizes it (in the so-called delayed mechanism by which the stalled prompt supernova shock is reactivated by an increase in the pressure behind it due to neutrino energy deposition [228, 27]), and convective overturn which rapidly transports energy into the shocked region [53, 26] (and which can lead to large-scale deviations from spherically symmetric explosions).


Click on thumbnail to view image

Figure 3: Schematic profiles of the velocity versus radius at three different times during core collapse: at the point of ``last good homology'', at bounce and at the time when the shock wave has just detached from the inner core.

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 [146Jump To The Next Citation Point In The Article].

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 [221Jump To The Next Citation Point In The Article] 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 [43Jump To The Next Citation Point In The Article] developed a code in which the neutrino transport was handled with an approximation, the multigroup flux-limited diffusion. Baron et al. [20Jump To The Next Citation Point In The Article] 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, tex2html_wrap_inline3937 (whose precise value, nowadays preferred around tex2html_wrap_inline3939  [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 [43Jump To The Next Citation Point In The Article], multigroup flux-limited diffusion for the neutrino transport. Van Riper [222Jump To The Next Citation Point In The Article] and Bruenn [44Jump To The Next Citation Point In The Article] verified that a softer supranuclear EOS, combined with general relativistic hydrodynamics, increases the chances for a prompt explosion. In [189Jump To The Next Citation Point In The Article, 188Jump To The Next Citation Point In The Article] and [136Jump To The Next Citation Point In The Article] 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 [18Jump To The Next Citation Point In The Article], [136] used maximal slicing coordinates. Miralles et al. [141Jump To The Next Citation Point In The Article], 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. [20Jump To The Next Citation Point In The Article], 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. [221Jump To The Next Citation Point In The Article, 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 tex2html_wrap_inline3941 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 [124Jump To The Next Citation Point In The Article], and later relativistic gravity [101, 33Jump To The Next Citation Point In The Article, 98Jump To The Next Citation Point In The Article, 184Jump To The Next Citation Point In The Article]. Martí et al. [124Jump To The Next Citation Point In The Article] 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. [105Jump To The Next Citation Point In The Article] 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 [125Jump To The Next Citation Point In The Article], where the characteristic fields of the one-dimensional (spherically symmetric) system were derived. The first astrophysical application was stellar collapse. In [33Jump To The Next Citation Point In The Article] 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].


Click on thumbnail to view movieClick on thumbnail to view movie

Click on thumbnail to view movieClick on thumbnail to view movie

Figure 4: Animations of a relativistic adiabatic core collapse using HRSC schemes. The simulations are taken from Ref. [184Jump To The Next Citation Point In The Article]: Velocity (movie 1, top-left), logarithm of the rest-mass density (movie 2, top-right), gravitational mass (movie 3, bottom-left) and lapse function squared (movie 4, bottom-right). See text for details of the initial model.

A comprehensive study of adiabatic one-dimensional core collapse using explicit HRSC schemes was presented in [184Jump To The Next Citation Point In The Article] (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 tex2html_wrap_inline3943, where tex2html_wrap_inline3517 is the bounce density and tex2html_wrap_inline3947 if tex2html_wrap_inline3949 and tex2html_wrap_inline3951 otherwise [221Jump To The Next Citation Point In The Article]. A set of animations of the simulations presented in [184Jump To The Next Citation Point In The Article] is included in the four MPEG movies in Fig.  4 . They correspond to the rather stiff Model B of [184Jump To The Next Citation Point In The Article]: tex2html_wrap_inline3953 and tex2html_wrap_inline3955 . The initial model is a white dwarf having a gravitational mass of tex2html_wrap_inline3957 . 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 tex2html_wrap_inline3961 . The gravitational mass of the inner core at the time of maximum compression is tex2html_wrap_inline3963 . The minimum value the quantity tex2html_wrap_inline3965, 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 tex2html_wrap_inline3969).

Studies of black hole formation in strong field hydrodynamical simulations of gravitational collapse are less frequent [221Jump To The Next Citation Point In The Article, 196Jump To The Next Citation Point In The Article, 24Jump To The Next Citation Point In The Article, 156Jump To The Next Citation Point In The Article]. In [221Jump To The Next Citation Point In The Article], 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 tex2html_wrap_inline3971 and 1.96 tex2html_wrap_inline3971 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 tex2html_wrap_inline3975 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 [221Jump To The Next Citation Point In The Article]) was developed by Miller and Motta [139Jump To The Next Citation Point In The Article] and later by Baumgarte, Shapiro and Teukolsky [24Jump To The Next Citation Point In The Article]. The novelty of these codes was the use of an outgoing null coordinate u (an ``observer-time'' coordinate, as suggested previously by [96Jump To The Next Citation Point In The Article]) 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 [221Jump To The Next Citation Point In The Article]. 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 tex2html_wrap_inline3971 polytrope with tex2html_wrap_inline3491 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.


Click on thumbnail to view image

Figure 5: Evolution of the black hole apparent horizon mass during the spherical accretion of a self-gravitating perfect fluid whose density is parametrized according to tex2html_wrap_inline3985 with tex2html_wrap_inline3511, tex2html_wrap_inline3513 and tex2html_wrap_inline3515 . The background density tex2html_wrap_inline3517 is given by the spherical Michel solution [137Jump To The Next Citation Point In The Article]. The grid extends from tex2html_wrap_inline3519 to tex2html_wrap_inline3521 . The mass of the apparent horizon shows a rapid increase in the first 15 M (enlarged in the insert), in coincidence with the most dynamical accretion phase. The slow, quasi-steady growth at later times is the quiet response of the black hole to the low mass accretion rate imposed at the world tube.

Recently, Neilsen and Choptuik [156Jump To The Next Citation Point In The Article] have studied critical collapse of a spherically symmetric perfect fluid with the EOS tex2html_wrap_inline4001 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 tex2html_wrap_inline3731 (tex2html_wrap_inline4005), 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 [171Jump To The Next Citation Point In The Article]. Using the formulation outlined in Section  2.2.2 and HRSC schemes, Papadopoulos and Font [171Jump To The Next Citation Point In The Article] 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.

4.1.2 Multidimensional studies 

The use of general relativistic axisymmetric codes in numerical astrophysics has been largely devoted to the study of the gravitational collapse and bounce of a rotating stellar core and the subsequent emission of gravitational radiation. These studies of axisymmetric scenarios started in the early eighties when the available computer resources were still inadequate to attempt three dimensional simulations. The collapse scenario was by then the main motivation for code builders in numerical relativity, as it was considered one of the most important gravitational wave sources, in the same way as the coalescence and merging of neutron star (and/or black hole) binaries is driving the development of fully three dimensional codes in the nineties (and possibly will continue during the forthcoming years) in correlation with the increasing availability of computer power.

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, 211Jump To The Next Citation Point In The Article, 177, 212Jump To The Next Citation Point In The Article], 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, tex2html_wrap_inline4011, 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 tex2html_wrap_inline4013 was radiated as gravitational waves.

Evans [65Jump To The Next Citation Point In The Article], 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 [65Jump To The Next Citation Point In The Article] 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 (tex2html_wrap_inline4017) 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 [197Jump To The Next Citation Point In The Article, 198Jump To The Next Citation Point In The Article, 199Jump To The Next Citation Point In The Article, 200Jump To The Next Citation Point In The Article, 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 tex2html_wrap_inline4023 . 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 [170Jump To The Next Citation Point In The Article] 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 [170Jump To The Next Citation Point In The Article] 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.


Click on thumbnail to view movie

Figure 6: Time evolution of the accretion/collapse of a quadrupolar shell onto a Schwarzschild black hole. The left panel shows isodensity contours and the right panel the associated gravitational waveform. The shell, initially centered at tex2html_wrap_inline3525, is gradually accreted by the black hole, a process which perturbs the black hole and triggers the emission of gravitational radiation. After the burst, the remaining evolution shows the decay of the black hole quasi-normal mode ringing. By the end of the simulation a spherical accretion solution is reached. Further details are given in [170Jump To The Next Citation Point In The Article]. (Editorial note: this movie has been reduced in size substantially compared to the original version provided by the author. To download the author's original version please go to this article's download page.)

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 tex2html_wrap_inline4027 with tex2html_wrap_inline4029 and tex2html_wrap_inline4031 . Additionally tex2html_wrap_inline4033 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 [203Jump To The Next Citation Point In The Article], for the head-on collision of two neutron stars [140Jump To The Next Citation Point In The Article], or for the coalescence of neutron star binaries [202Jump To The Next Citation Point In The Article, 205Jump To The Next Citation Point In The Article] (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 tex2html_wrap_inline4037 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 [107Jump To The Next Citation Point In The Article] for a more detailed description of a more energetic initial model.


Click on thumbnail to view image

Figure 7: Animation of the time evolution of the entropy in a core collapse supernova explosion [107]. The movie shows the evolution within the innermost 3000 km of the star and up to 220 ms after core bounce. See text for explanation.

4.2 Accretion onto black holes4 Simulations of Astrophysical Phenomena4 Simulations of Astrophysical Phenomena

image Numerical Hydrodynamics in General Relativity
José A. Font
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to livrev@aei-potsdam.mpg.de