4.2 Accretion onto black holes4 Hydrodynamical Simulations in Relativistic 4 Hydrodynamical Simulations in Relativistic

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 [172, 173]. Such codes are designed to integrate the coupled system of Einstein field equations and relativistic hydrodynamics, sidestepping Newtonian approaches.

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 [191Jump To The Next Citation Point In The Article, 134Jump 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 non-spherical collapse (see [206] 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., [262Jump To The Next Citation Point In The Article, 259Jump To The Next Citation Point In The Article, 89Jump To The Next Citation Point In The Article, 260Jump To The Next Citation Point In The Article] and references therein).

4.1.1 Core collapse supernovae 

At the end of their thermonuclear evolution, massive stars in the (Main Sequence) mass range from tex2html_wrap_inline5180 to tex2html_wrap_inline5182 develop a core composed of iron group nuclei, which becomes dynamically unstable against gravitational collapse. The iron core collapses to a neutron star or a black hole releasing gravitational binding energy of the order tex2html_wrap_inline5184, which is sufficient to power a supernova explosion. The standard model of type II/Ib/Ic 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 process is depicted in Figure  3 . The shock is produced after the bounce of the inner core when its density exceeds nuclear density. Numerical simulations have tried to elucidate whether this shock is strong enough to propagate outward 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 through the outer layers of the star. The accepted scenario, which has emerged from numerical simulations, contains the following two necessary ingredients for any plausible explosion mechanism (see [134] for a review): (i) the existence of the shock wave together with neutrino 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 post-shock pressure due to neutrino energy deposition [298, 30]); and (ii) the convective overturn which rapidly transports energy into the shocked region [59, 29] (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.

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

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 [299] 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 [293Jump 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 [49Jump 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. [24Jump To The Next Citation Point In The Article] 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, tex2html_wrap_inline5186 (its precise value, nowadays preferred around tex2html_wrap_inline5188  [107], is still a matter of debate). Mayle et al. [174] computed neutrino spectra from stellar collapse to stable neutron stars for different collapse models using, as in [49Jump To The Next Citation Point In The Article], multigroup flux-limited diffusion for the neutrino transport. Van Riper [294Jump To The Next Citation Point In The Article] and Bruenn [50Jump 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 [249Jump To The Next Citation Point In The Article, 248Jump To The Next Citation Point In The Article] and [178Jump 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 [249, 248] the neutrino transport is described by moments of the general relativistic Boltzmann equation in polar slicing coordinates [22Jump To The Next Citation Point In The Article], [178] used maximal slicing coordinates. Miralles et al. [184Jump To The Next Citation Point In The Article], 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. [24Jump To The Next Citation Point In The Article], who used a phenomenological EOS. Swesty et al. [280] 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., [293Jump To The Next Citation Point In The Article, 24, 294, 50, 184]). Swesty and coworkers used a finite temperature compressible liquid drop model EOS for the nuclear matter component [148]. 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_inline5190 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. [103], which used an Eulerian second-order PPM scheme (see [191] for a review of the present status). There are only a few relativistic simulations, so far restricted to spherical symmetry [129, 244Jump To The Next Citation Point In The Article, 211].


Click on thumbnail to view movieClick on thumbnail to view movie

Click on thumbnail to view movieClick on thumbnail to view movie

Figure 4: Movie showing the animations of a relativistic adiabatic core collapse using HRSC schemes ( snapshots of the radial profiles of various variables are shown at different times). The simulations are taken from [244Jump To The Next Citation Point In The Article]: velocity (top-left), logarithm of the rest-mass density (top-right), gravitational mass (bottom-left), and lapse function squared (bottom-right). See text for details of the initial model. Visualization by José V. Romero.

Martí et al. [162Jump To The Next Citation Point In The Article] implemented Glaister's approximate Riemann solver [106] 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. [136Jump 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. [162]. The state-of-the-art implementation of the tensorial artificial viscosity terms in [136], 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 [163Jump 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 [38Jump 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 [39], 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 [244Jump To The Next Citation Point In The Article] (see [309] 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 tex2html_wrap_inline5192, where tex2html_wrap_inline5194 is the bounce density, and tex2html_wrap_inline5196 if tex2html_wrap_inline5198 and tex2html_wrap_inline5200 otherwise [293Jump To The Next Citation Point In The Article]. A set of animations of the simulations presented in [244Jump To The Next Citation Point In The Article] is included in the four movies in Figure  4 . They correspond to the rather stiff Model B of [244Jump To The Next Citation Point In The Article]: tex2html_wrap_inline5202, tex2html_wrap_inline5204, and tex2html_wrap_inline5206 . The initial model is a white dwarf having a gravitational mass of tex2html_wrap_inline5208 . 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 tex2html_wrap_inline5212 . The gravitational mass of the inner core at the time of maximum compression is tex2html_wrap_inline5214 . The minimum value that the quantity tex2html_wrap_inline5216 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 tex2html_wrap_inline5218).

Axisymmetric Newtonian simulations.    Beyond spherical symmetry most of the existing simulations of core collapse supernova are Newtonian. Axisymmetric investigations have been performed by [190Jump To The Next Citation Point In The Article, 37, 42Jump To The Next Citation Point In The Article] 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, 319Jump To The Next Citation Point In The Article] have performed Newtonian parameter studies of the axisymmetric collapse of rotating polytropes.


Click on thumbnail to view movie

Figure 5: Movie showing the animation of the time evolution of the entropy in a core collapse supernova explosion [138Jump To The Next Citation Point In The Article]. The movie shows the evolution within the innermost tex2html_wrap_inline4636 of the star and up to tex2html_wrap_inline4638 after core bounce. See text for explanation. Visualization by Konstantinos Kifonidis.

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 [189Jump To The Next Citation Point In The Article]). 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 tex2html_wrap_inline4638 after core bounce and shock formation (Figure  5 depicts just one intermediate snapshot at tex2html_wrap_inline5226). The movie shows the evolution of the entropy within the innermost tex2html_wrap_inline4636 of the star.

The initial data used in these calculations is taken from the tex2html_wrap_inline5230 pre-collapse model of a blue supergiant star in [308]. The computations begin tex2html_wrap_inline5232 after core bounce from a one-dimensional model of [51]. This model is obtained with the one-dimensional general relativistic code mentioned previously [49], 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 [138] 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 [301] first computed neutron star bounces of tex2html_wrap_inline5234 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, 66Jump To The Next Citation Point In The Article, 67Jump To The Next Citation Point In The Article], who computed the gravitational radiation emitted in such events. The Einstein equations were formulated using the so-called conformally flat metric approximation [304Jump To The Next Citation Point In The Article]. 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 [66Jump To The Next Citation Point In The Article].

Dimmelmeier et al. [67Jump To The Next Citation Point In The Article] 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 tex2html_wrap_inline5238 . Collapse is initiated by decreasing the adiabatic index to some prescribed fixed value tex2html_wrap_inline5240 with tex2html_wrap_inline5242 . 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 [319Jump To The Next Citation Point In The Article] - 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 [67Jump To The Next Citation Point In The Article], 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.


Click on thumbnail to view movie
Reduced size (4.6 MB)
    Click on thumbnail to view movie
Original size (12.6 MB)

Figure 6: Animation showing the time evolution of a relativistic core collapse simulation (model A2B4G1 of [67Jump To The Next Citation Point In The Article]). Left: velocity field and isocontours of the density. Right: gravitational waveform (top) and central density evolution (bottom). The model exhibits a multiple bounce collapse (fizzler) with a type II signal. The camera follows the multiple bounces. Visualization by Harald Dimmelmeier.

The movie presented in Figure  6 shows the time evolution of a multiple bounce model (model A2B4G1 in the notation of [66]), 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 [67Jump To The Next Citation Point In The Article] can be found at the MPA's website [189Jump To The Next Citation Point In The Article].

The relativistic models analyzed by Dimmelmeier et al. [67Jump To The Next Citation Point In The Article] cover almost the same range of gravitational wave amplitudes (tex2html_wrap_inline5244 for a source at a distance of tex2html_wrap_inline5246) and frequencies (tex2html_wrap_inline5248) as the corresponding Newtonian ones [319Jump To The Next Citation Point In The Article]. Averaged over all models, the total energy radiated in the form of gravitational waves is tex2html_wrap_inline5250 in the relativistic case, and tex2html_wrap_inline5252 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. [67] is available at the MPA's website [189]. 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 [260Jump To The Next Citation Point In The Article].

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 [262Jump To The Next Citation Point In The Article] and supermassive stars [265Jump To The Next Citation Point In The Article] (see Section  4.1.2), for the head-on collision of two neutron stars [183Jump To The Next Citation Point In The Article], and for the coalescence of neutron star binaries [258Jump To The Next Citation Point In The Article, 266Jump To The Next Citation Point In The Article, 267Jump To The Next Citation Point In The Article] (see Section  4.3).

Concerning Newtonian studies, Bonazzola and Marck [42] 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 [190]). 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 [319] 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 [261Jump To The Next Citation Point In The Article] (see Section  4.3.1).

4.1.2 Black hole formation 

Apart from a few one-dimensional simulations, most numerical studies of general relativistic gravitational collapse leading to black hole formation have used Wilson's formulation of the hydrodynamics equations and finite difference schemes with artificial viscosity. The use of conservative formulations and HRSC schemes to study black hole formation, both in two and three dimensions, began rather recently.

Spherically symmetric simulations.    Van Riper [293Jump To The Next Citation Point In The Article], 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 tex2html_wrap_inline5254 and tex2html_wrap_inline5256 cores.

In [255Jump To The Next Citation Point In The Article], 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 [255] could follow relativistic configurations for many collapse timescales tex2html_wrap_inline5258 after the initial appearance of an event horizon.

A Lagrangian scheme based on the Misner and Sharp [185] formulation for spherically symmetric gravitational collapse (as in [293Jump To The Next Citation Point In The Article]) was developed by Miller and Motta [181Jump To The Next Citation Point In The Article] and later by Baumgarte et al. [26Jump 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 [124Jump To The Next Citation Point In The Article]) 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 [124] (or, alternatively, the Misner-Sharp equations) were solved by an explicit finite difference scheme similar to the one used by [293Jump To The Next Citation Point In The Article]. In [181Jump To The Next Citation Point In The Article], the collapse of an unstable polytrope to a black hole was first achieved using null slicing. In [26Jump To The Next Citation Point In The Article], the collapse of a tex2html_wrap_inline4658 polytrope with tex2html_wrap_inline4616 was compared to the result of [293], 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. [156] analyzed the gravitational collapse of supermassive stars in the range tex2html_wrap_inline5268 . 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 tex2html_wrap_inline5270 -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 tex2html_wrap_inline5272, 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 [270Jump To The Next Citation Point In The Article] 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 tex2html_wrap_inline5278 to tex2html_wrap_inline5280 (tex2html_wrap_inline5282) was built. Using a compactified spacetime foliation with outgoing null cones, Siebel et al. [270] studied the fate of the relativistic stars when they are hit by a strong scalar field pulse with a mass of the order of tex2html_wrap_inline5284, 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 [193], 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) [195]. 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 tex2html_wrap_inline5288 ``core'' of a massive star with different amounts of rotational energy and an initial central density of tex2html_wrap_inline5290, up to the formation of a rotating black hole. However, the resolution affordable due to limitations in computer resources (tex2html_wrap_inline5292 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 [194] (see also [197]) considered a configuration consisting of a neutron star (tex2html_wrap_inline5294, tex2html_wrap_inline5296) with an accreted envelope of tex2html_wrap_inline5298, 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, 276Jump To The Next Citation Point In The Article, 228, 277Jump To The Next Citation Point In The Article], studied the collapse of rotating relativistic (tex2html_wrap_inline5234) 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 tex2html_wrap_inline5304, and radius tex2html_wrap_inline5306 . 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 tex2html_wrap_inline5308, and that the general waveform shape was nearly independent of the details of the collapse.

Evans [79Jump To The Next Citation Point In The Article], building on previous work by Dykema [74], 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 (tex2html_wrap_inline5312), where some fields diverge due to the coordinate singularities. Evans made important contributions toward regularizing the gravitational field equations in such situations [79]. 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 [275], although no applications in relativistic astrophysics have yet been presented.

Recently, Alcubierre et al. [6] 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 [259Jump To The Next Citation Point In The Article] investigated the effects of rotation on the criterion for prompt adiabatic collapse of rigidly and differentially rotating (tex2html_wrap_inline5234) 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 [259] 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 [265] 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 tex2html_wrap_inline5316 . 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 [305], such as the (vacuum) code presented in [109], 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 [269Jump To The Next Citation Point In The Article], 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 (tex2html_wrap_inline5234) uniformly rotating neutron stars to rotating black holes, using the code discussed in Section  3.3.1, are presented in [262Jump To The Next Citation Point In The Article]. The simulations show no evidence for massive disk formation or outflows, which can be related to the moderate initial compactness of the stellar models (tex2html_wrap_inline5320). A proof-of-principle of the capabilities of the code discussed in Section  3.3.2 to study black hole formation is presented in [89Jump To The Next Citation Point In The Article], where the gravitational collapse of a spherical, unstable, relativistic polytrope is discussed. Similar tests with differentially rotating polytropes are given in [73] for a recent artificial viscosity-based, three-dimensional general relativistic hydrodynamics code.

4.1.3 Critical collapse 

Critical phenomena in gravitational collapse were first discovered numerically by Choptuik in spherically symmetric simulations of the massless Klein-Gordon field minimally coupled to gravity [56Jump To The Next Citation Point In The Article]. Since then, critical phenomena arising at the threshold of black hole formation have been found in a variety of physical systems, including the perfect fluid model (see [114] for a review).

Evans and Coleman [81] first observed critical phenomena in the spherical collapse of radiation fluid, i.e., a fluid obeying an EOS tex2html_wrap_inline5322 with tex2html_wrap_inline4616 and tex2html_wrap_inline5326 . The threshold of black hole formation was found for a critical exponent tex2html_wrap_inline5328, in close agreement with that obtained in scalar field collapse [56]. Their study used Schwarzschild coordinates in radial gauge and polar slicing, and the hydrodynamic equations followed Wilson's formulation. Subsequently, Maison [159], using a linear stability analysis of the critical solution, showed that the critical exponent varies strongly with the parameter tex2html_wrap_inline5330 of the EOS. More recently, Neilsen and Choptuik [203], 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 [204] made it possible to find evidence of the existence of critical solutions for large values of the adiabatic index tex2html_wrap_inline4726 (tex2html_wrap_inline5334), extending the previous upper limit.

4.2 Accretion onto black holes4 Hydrodynamical Simulations in Relativistic 4 Hydrodynamical Simulations in Relativistic

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