Since early investigations the numerical study of gravitational collapse has been neatly divided between two main communities. First, the computational astrophysics community has traditionally focused on the physics of the (type II/Ib/Ic) Supernova mechanism, i.e., the 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), along with the incorporation of methodology for the magnetic field. These studies are currently performed routinely in multiple 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, and only recently this situation has started to change (see below). In this approach the emission of gravitational radiation is computed through the quadrapole formula. For reviews of the current status in this direction see [271, 183, 422] 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 nonspherical collapse (see  for a 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 has been given to the details of the microphysics of the collapse (yet again this situation is starting to change). In addition, this approach has mostly considered gravitational collapse leading to black-hole formation, employing relativistic hydrodynamics, MHD, and gravity. Both approaches, which had barely interacted over the years except for a handful of simulations in spherical symmetry, have begun a much closer interaction very recently. This has been mainly the result of the progress achieved in numerical relativity regarding long-term stable formulations of the Einstein equations, as we will illustrate in the following Section 5.1.1.
However, before proceeding any further, we must point out that the current version of this article does not cover, as earlier versions did, 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 . 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. The interested reader is directed to  and references therein for details.
At the end of their thermonuclear evolution, massive stars in the (Main Sequence) mass range from to 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 , 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. In the accepted scenario, which has emerged from numerical simulations, the existence of the shock wave together with neutrino heating that re-energizes it (in the delayed mechanism by which the stalled-prompt supernova shock is reactivated by an increase in the post-shock pressure due to neutrino energy deposition [414, 44]), and convective overturn, which rapidly transports energy into the Ledoux-unstable shocked region [81, 43] (and which can lead to large-scale deviations from spherically-symmetric explosions), are all necessary ingredients for any plausible explosion mechanism. The most recent simulations have also highlighted the importance of a low-mode instability in the standing accretion shock for large-scale explosion asymmetries, pulsar kicks, and for the development of neutrino-driven explosions (see [422, 184] for details on the degree of sophistication achieved in present-day supernova modelling).
However, the path to reach such conclusions has mostly been delineated from numerical simulations in one spatial dimension. Fully relativistic simulations of microphysically-detailed core collapse off spherical symmetry are only beginning to become available, and they might well introduce some modifications. The broad picture described above has been demonstrated in multiple dimensions using sophisticated Newtonian models. Nevertheless, while studies based upon Newtonian physics are highly developed nowadays, state-of-the-art simulations still fail, broadly speaking, to generate successful supernova explosions under generic conditions.
There are only a few fully relativistic simulations with HRSC schemes available in the literature [178, 339, 294], which, however, do not account for neutrino transport. On the other hand, there also exist state-of-the-art neutrino-radiation–hydrodynamics codes in spherical symmetry, which incorporate the relativistic effects of the gravitational field in an exact or approximate manner. The code of the Oak Ridge-Basel group (AGILE-BOLTZTRAN) belongs to the former class and consists of a general-relativistic time-implicit discrete-angle Boltzmann solver, coupled to a conservative, general relativistic, time-implicit, hydrodynamics TVD solver. On the other hand, the VERTEX code of the Garching group calculates neutrino transport by a variable Eddington-factor closure of neutrino energy, number, and momentum equations, and the hydrodynamics step is based on conservative high-order Riemann solvers, as those used in its precursor code PROMETHEUS. However, the relativistic effects of the gravitational field are only incorporated in VERTEX through an effective gravitational potential. Comparisons of these two supernova codes are reported in .
Despite the incomplete physics employed in adiabatic approximations of stellar core collapse it is still useful to present them here in order to get acquainted with the complex hydrodynamics involved in the problem. A comprehensive study of adiabatic, one-dimensional, general relativistic core collapse using explicit upwind HRSC schemes was presented in  (see  for a similar computation using implicit schemes). In this investigation the equations for the hydrodynamics and the geometry are written using radial gauge polar slicing (Schwarzschild-type) coordinates. The collapse is modeled with an ideal gas EOS with a nonconstant adiabatic index, which depends on the density as , where is the bounce density and if and otherwise . A set of animations of the simulations presented in  is included in the four movies in Figure 4. They correspond to the rather stiff Model B of : and . The initial model is a white dwarf with a gravitational mass of . The animations show the time evolution of the radial profiles of the following fields: velocity (Movie 4a), logarithm of the rest-mass density (Movie 4b), gravitational mass (Movie 4c), and the square of the lapse function (Movie 4d).
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 (), is . The gravitational mass of the inner core at the time of maximum compression is . The minimum value that the quantity reaches is , 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 ).
Among the more detailed multi-dimensional nonrelativistic hydrodynamic simulations are those performed by the MPA/Garching group (an online sample can be found at their website ). As mentioned earlier, these include advanced microphysical EOS, state-of-the-art treatments of neutrino physics, and they 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 220 ms after core bounce and shock formation (Figure 5 depicts just one intermediate snapshot at 130 ms). The movie shows the evolution of the entropy within the innermost 3000 km of the star.
The initial data used in these calculations is taken from the pre-collapse model of a blue super-giant star in . The computations begin 20 ms after core bounce from a one-dimensional model of . This model is obtained with the one-dimensional general-relativistic code mentioned previously , which includes a detailed treatment of the neutrino physics and a realistic EOS, and which is able to follow core collapse, bounce and the associated formation of the supernova shock. Because of neutrino cooling and energy losses due to the dissociation of iron nuclei, the shock initially stalls inside the iron core.
The movie shows how the stalling shock is “revived” by neutrinos streaming from the outer layers of the hot, nascent neutron star in the center. A negative entropy gradient builds up between the so called “gain-radius”, which is the position where the cooling of hot gas by neutrino emission switches into net energy gain by neutrino absorption, and the shock further out. Convective instabilities, which are characterized by large rising blobs of neutrino-heated matter and cool, narrow downflows, therefore, set in. 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 directed to  for a more detailed description of a more energetic initial model.
Multidimensional models, such as those developed by the MPA/Garching group or the Oak Ridge–Basel group, have shown that such strong radial mixing, as in the previous movie, requires the development of hydrodynamic instabilities near the core and even in an early stage of the explosion. In recent years axisymmetric numerical models have highlighted the fact that the standing accretion shock is generically unstable to nonradial deformations even when convection is absent or highly damped, to the point that the term “standing accretion shock (advective/acoustic) instability” (SASI) has been coined in the literature (see  and references therein). The SASI shows a preferential growth of low ( and ) shock-deformation modes at the time of shock revival, leading to large-amplitude bipolar oscillations, which have important implications for large-scale explosion asymmetries, pulsar kicks, and the development of neutrino driven explosions.
A new line of investigation led by Burrows et al.  has recently claimed successful supernova explosions by resorting to an alternative mechanism based on acoustic power generated in the inner core. The simulations have been performed with the Newtonian, two-dimensional radiation/hydrodynamic code VULCAN/2D with multi-group, flux-limited diffusion. Such acoustic power is generated in two phases. It starts in the first 100 ms following core bounce in the inner turbulent region and, most importantly, it becomes more significant at later times, by the excitation and sonic damping of core -mode oscillations. The simulations of  for an progenitor show that an mode grows at late times to be prominent around roughly 500 ms after bounce to drive a successful explosion. This mechanism has been met with skepticism by .
In addition to understanding the explosion mechanism, one of the prime motivations for performing multidimensional core-collapse simulations is to compute theoretical gravitational-radiation waveforms, which may be confronted by experimental data currently being collected by a number of detectors. As early as 1982, Müller  obtained the first numerical evidence of the low gravitational wave efficiency of the core collapse scenario, with an energy smaller than being radiated as gravitational waves during core bounce. In the 1990s, the first Newtonian parameter studies of the axisymmetric collapse of rotating polytropes were performed by [123, 425, 439], which provided the first investigations of the gravitational wave signals from core bounce and early post-bounce phase. In core collapse events, gravitational waves are initially dominated by a burst associated with the hydrodynamic bounce. If rotation is present, the post-bounce wave signal shows large amplitude oscillations associated with pulsations in the collapsed core [439, 328] and (possibly) low rotational dynamic instabilities [306, 305] ( and being the kinetic rotational energy and the potential energy, respectively). On the other hand, it has recently been shown that gravitational waves from convection have preeminence over the purely burst signal of the core bounce on longer timescales . In this case, the inherent long emission timescale can yield high energy in a continuous signal.
It is well known that neutron stars have intense magnetic fields () or even larger at birth (). In extreme cases such as magnetars, the magnetic field can be so strong as to affect the internal structure of the star . The presence of such intense magnetic fields render magneto-rotational core collapse simulations mandatory. We note that as early as 1979 magneto-rotational core collapse was already proposed by  as a plausible supernova explosion mechanism.
In recent years, an increasing number of authors have performed axisymmetric magneto-rotational core collapse simulations (within the ideal MHD limit) employing Newtonian treatments of magneto-hydrodynamics, gravity, and microphysics (see, e.g., [296, 69] and references therein). Relativistic magneto-rotational core collapse simulations have only very recently become available (see below). The weakest point of all existing magneto-rotational core collapse simulations to date is the fact that the strength and distribution of the initial magnetic field in the core are unknown. If the magnetic field is initially weak, there exist several mechanisms which may amplify it to values, which can have an impact on the dynamics, among them differential rotation (-dynamo) and the magneto-rotational instability (MRI hereafter). The former transforms rotational energy into magnetic energy, winding up any seed poloidal field into a toroidal field, while the latter, which is present as long as the radial gradient of the angular velocity of the fluid is negative, leads to exponential growth of the field strength.
Specific magneto-rotational effects on the gravitational wave signature were first studied in detail by  and , who found differences with purely hydrodynamic models only for very strong initial fields (). The most exhaustive parameter study of magneto-rotational core collapse to date has been carried out very recently by [296, 295]. The axisymmetric simulations of [296, 295] have employed rotating polytropes, Newtonian (and modified Newtonian) gravity and hydrodynamics, and ad-hoc initial poloidal magnetic field distributions (as no self-consistent solution is yet known). These authors have found that for weak initial fields (, which is the most relevant case, astrophysics-wise) there are no differences in the collapse dynamics nor in the resulting gravitational wave signal, when comparing with purely hydrodynamic simulations. However, strong initial fields () manage to slow down the core efficiently (leading even to retrograde rotation in the proto-neutron star), which causes qualitatively different dynamics and gravitational wave signals. For some models,  even find highly bipolar, jet-like outflows. The creation and propagation of MHD jets is also studied in the recent multigroup, radiation MHD simulations of supernova core collapse in the context of rapid rotation performed by . These simulations suggest that for rotating cores a supernova explosion might be followed by a secondary, weak MHD jet explosion, which might be generic in the collapsar model of GRBs.
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 was performed by Dimmelmeier et al. [93, 94, 95, 96], who computed the gravitational radiation emitted in such events. The Einstein equations were formulated using the conformally-flat metric approximation  (CFC hereafter). 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 numerical code are given in  (see also Section 4.4). These simulations employed simplified models to account for the thermodynamics of the process, in the form of a polytropic EOS conveniently modified to account for the stiffening of the matter once nuclear matter density is reached . The inclusion of relativistic effects results primarily in the possibility of achieving deeper effective potentials. Moreover, higher densities than in Newtonian models are reached during bounce, and the resulting proto-neutron star is more compact.
The simulations show that the three different types of rotational supernova core collapse (and their associated gravitational waveforms) identified in previous Newtonian simulations by Zwerger and Müller  – regular collapse, multiple bounce collapse, and rapid collapse (Types I, II, and III, respectively) – are also present in relativistic gravity, at least when the initial models are simple rotating polytropes. As mentioned before, the gravitational wave signal is characterized by a distinctive burst associated with the hydrodynamic bounce followed by a strongly-damped ring-down phase of the newly-born proto-neutron star. While indeed the maximum density reaches higher values in relativistic gravity than in Newtonian gravity, the gravitational wave signals were found to be of comparable amplitudes. In fact, the simulations performed by  reveal that this is a general trend, the gravitational radiation emitted in the relativistic models being strikingly similar to those previously obtained from Newtonian simulations . On average, thus, the CFC simulations of  show that the gravitational wave signals of relativistic models have similar amplitude but somewhat higher frequency than their Newtonian counterparts .
Among the most interesting results of these investigations was the fact that in rotational core collapse, the collapse type can change from multiple centrifugal bounce (Type II) to standard single-bounce collapse (Type I). Centrifugal bounce was highly suppressed in relativistic gravity, yet it was still possible for simple EOS (but see below for the situation regarding improved microphysics). The main conclusion of  is that only the gravitational signal of a Galactic supernova (i.e., within a distance of about 10 kpc) might be unambiguously detectable by first generation detectors. Signal recycling techniques in next generation detectors may be needed, however, for successful detection of more distant core collapse events (with an increased event rate), up to distances of the Virgo cluster. An online waveform catalogue for all models analyzed by Dimmelmeier et al.  is available at the MPA’s website .
The movie presented in Figure 6 shows the time evolution of a multiple bounce model (model A2B4G1 in the notation of ), with a type II gravitational wave signal. The left panel shows isocontours of the logarithm of the density together with the corresponding velocity field distribution. The two panels on the right show the time evolution of the gravitational wave signal (top panel) and of the central rest-mass density. In the animation the “camera” follows the multiple bounces by zooming in and out. The panels on the right show how each burst of gravitational radiation coincides with the time of bounce of the collapsing core. The oscillations of the nascent proto-neutron star are therefore imprinted on the gravitational waveform. Additional animations of the simulations performed by  can be found at the MPA’s website .
Relativistic simulations with improved dynamics and gravitational waveforms were reported by , who used an approximation for the field equations which they called CFC+. This approximation incorporates additional degrees of freedom to CFC such that the spacetime metric is exact up to the second post-Newtonian order. As in the CFC simulations of  a simple (modified polytrope) EOS was used, and the initial models were unmagnetized. The CFC simulations of  cover the basic morphology and dynamics of core collapse types studied by , including the extreme case of a core with strong differential rotation and torus-like structure. In either case no significant differences were found in all models investigated by , both regarding the dynamics and gravitational wave signals. Therefore, second post-Newtonian corrections to CFC do not significantly improve the results when simulating the dynamics of core collapse to a neutron star. (It remains to be seen whether this conclusion changes in the strongest gravitational fields as, e.g., in the collapse to a black hole.)
Comparisons of the CFC approach with fully general-relativistic simulations (employing the BSSN formulation) have been reported by  in the context of axisymmetric core collapse simulations. As in the case of CFC+, the differences found in both, the collapse dynamics and the gravitational waveforms are minute, which highlights the suitability of CFC (and CFC+) for performing accurate simulations of these scenarios without the need for solving the full system of Einstein’s equations. It is also worth pointing out the detailed comparison of a comprehensive set of models (some including collapse to black holes) reported recently by , to which the interested reader is directed for further information. Owing to the excellent approximation offered by CFC for studying core collapse, extensions to improve the microphysics of the numerical modelling, through the incorporation of microphysical EOS and simplified neutrino treatments, as well as extensions to account for three spatial dimensions, have been reported in the last few months.
The most realistic simulations of rotational core collapse in general relativity so far (both in axisymmetry and in three dimensions, but without incorporating magnetic fields) are those performed by [305, 98]. These simulations employ both the CFC and BSSN formulations of Einstein’s equations, and have been carried out using state-of-the-art numerical codes, both in spherical polar coordinates (the CoCoNut code of ) and in Cartesian coordinates (the CACTUS/WHISKY codes; see [244, 27, 305] and references therein). As in the earlier works of [96, 365, 74] the gravitational wave information from these collapse simulations is extracted using (variations of) the Newtonian standard quadrapole formula.
The initial models of [305, 98] use nonrotating pre-supernova models to which an artificially-parameterized rotation is added. In addition, the relevant microphysics of core collapse is accounted for by employing a finite-temperature nuclear EOS  together with an approximate (parameterized) neutrino treatment with deleptonization and neutrino pressure effects included, as recently suggested by . In the simulations of [305, 98] the electron fraction profiles are obtained from general-relativistic spherical models with Boltzmann neutrino transport, and are used as input data for the time-dependent electron fraction prescription suggested by  (with no analytic fit). The physics included in the modelling provide a reasonably good approximation up until shortly after core bounce. Therefore, the setup is valid for obtaining reliable gravitational wave signals from core bounce.
The simulations reported by [305, 98] allow for a straight comparison with models that have simple (polytropic) EOS and the same rotation parameters as those of [96, 365, 74]. This sheds light on the influence of rotation and of the EOS on the gravitational waveforms. This comparison shows that for microphysical EOS the influence of rotation on the collapse dynamics and waveforms is somewhat less pronounced than in the case of simple EOS. In particular, the most important result of these investigations is the suppression of core collapse with multiple centrifugal bounces and its associated Type II gravitational waveforms, as that shown in Figure 6, which is explained as mainly due to the influence of deleptonization on the collapse dynamics, which removes energy from the system. As a result, the improved simulations of [98, 305] reveal an enhanced uniformity of the gravitational wave signal from core bounce as, in particular, the low frequency signals are suppressed (as they were associated with the multiple centrifugal bounce scenario). In addition, signals for slow and moderate rotation have almost identical frequencies, that is, the signal peaked at a frequency of 670 Hz for the models considered. These important results have implications for signal analysis and signal inversion.
On the other hand, to further improve the realism of core collapse simulations in general relativity, the incorporation of magnetic fields in numerical codes via solving the MHD equations is also currently being undertaken [362, 75, 76]. The work of , which is based on the conservative formulation of the GRMHD equations discussed in Section 3.1.2, is focused on the collapse of initially strongly magnetized cores (). Although these values are probably astrophysically not relevant, because the stellar evolution models of  predict a poloidal field strength of in the progenitor, they permit one to resolve the scales needed for the MRI to develop, since the MRI length scale grows with the magnetic field. We note that, in MRI unstable regions, the instability grows exponentially for all length scales larger than a critical length scale , where is the Alfvén speed and is the angular velocity. The fastest growing MRI mode develops at length scales near on a typical time scale of , where is the distance to the rotation axis. Therefore, in order to numerically capture the MRI, one has to resolve length scales of about the size of , which can be prohibitively small to capture for initially small magnetic fields (say, on the order of meters for a progenitor), especially for three-dimensional computations, where the affordable grid resolution is still limited. The results of  show that the magnetic field is mainly amplified by the wind-up of the magnetic field lines by differential rotation. Consequently, the magnetic field is accumulated around the inner region of the proto-neutron star, and MHD outflow forms along the rotation axis removing angular momentum from the star.
A different approach is followed by [75, 76]. Their progenitors are chosen to be weakly magnetized (), which is in much better agreement with predictions from stellar evolution. Under these conditions it is appropriate to use the “passive” magnetic field approximation, in which the magnetic field entering into the energy-momentum tensor is negligible when compared with the fluid part, and the components of the anisotropic term of satisfy . With such simplification, the evolution of the magnetic field, governed by the induction equation, does not affect the dynamics of the fluid, which is governed solely by the hydrodynamics equations. However, the magnetic field evolution does depend on the fluid evolution, due to the presence of velocity components in the induction equation. This approximation has interesting numerical advantages as the seven eigenvalues of the GRMHD Riemann problem (entropy, Alfvén, and fast and slow magnetosonic waves) reduce to three. Using this approach, the work of  presents the first relativistic simulations of magneto-rotational core collapse, which take into account the effects of a microphysical EOS  and a simplified neutrino treatment. These effects have been incorporated into the code following the procedure employed by  and  in their unmagnetized simulations. The results show that for the core collapse models with microphysics the saturation of the magnetic field cannot be reached within dynamic time scales by winding up the poloidal magnetic field into a toroidal one. Cerdá-Durán et al.  have also estimated the effect of other amplification mechanisms including the MRI and several types of dynamos, concluding that for progenitors with astrophysically expected (i.e., weak) magnetic fields, the MRI is the only mechanism that could amplify the magnetic field on dynamic time scales. In addition, all microphysical models exhibit post-bounce convective overturn in regions surrounding the inner part of the proto-neutron star. Since this has a potential impact on enhancing the MRI, it deserves further investigation with more accurate neutrino treatment or alternative microphysical equations of state.
The restriction to the passive magnetic-field approximation employed by  in studying magneto-rotational core collapse of weakly magnetized progenitors can be justified if the MRI is indeed inefficient. Otherwise, an “active” magnetic field approach becomes necessary (as in the work of ). However, the use of active magnetic fields alone for core collapse simulations will probably not be sufficient to model all the effects amplifying the magnetic field, since the numerical resolution needed to correctly describe them (less than 10 m) is not affordable in current numerical simulations. In addition, most of the relevant effects should be investigated in three dimensions, turning the computational work that lies ahead into a daunting task.
Three-dimensional and fully general-relativistic simulations of core collapse have been performed by  and . The former study focuses on the likelihood of nonaxisymmetric dynamic instabilities in the proto-neutron stars that form following core collapse, for a comprehensive set of initial rotating polytropes with a hybrid (polytrope and thermal) EOS. In this work, the early stage of the collapse is performed in axisymmetry and only after the stellar core becomes compact enough the three-dimensional simulation is started, adding to the collapsed core a bar-mode nonaxisymmetric density perturbation. The simulations show that the dynamic bar-mode instability sets in if the progenitor already is rapidly rotating () and has a high degree of differential rotation, which seems unlikely according to stellar evolution models.
On the other hand,  have performed 3D simulations of general-relativistic core collapse within the BSSN framework and a finite-temperature nuclear EOS  and an approximate (parameterized) neutrino treatment . These simulations provide yet another confirmation of the satisfactory accuracy of CFC for the core-collapse supernova problem. All models investigated in 3D are identical to those evolved in axisymmetry with the CoCoNut code in . Typical simulation grids with nine fixed–mesh-refinement levels and extending up to 3000 km in length were used. All models were followed up to 20 ms after core bounce, which showed that they all stay essentially axisymmetric through bounce, in good agreement with the models of . One of the models, labelled E20A in , was further evolved until 70 ms after core bounce, as it showed the largest gravitational wave amplitude of the sample and a value of the ratio in the newly-formed proto-neutron star high enough to eventually develop low nonaxisymmetric (bar-mode-like) deformations. This had already been found in a similar model using Newtonian physics in .
The left panel of Figure 7 displays the time evolution of the normalized mode amplitudes of the four lowest azimuthal density modes () in the equatorial plane, and shows that the one-armed mode (red curve) becomes dominant from ms after core bounce. Correspondingly, the right panel of Figure 7 plots the gravitational wave strains and along the polar axis, where the late-time increase in amplitude during the growth of the instability is visible. An animation of the entire evolution of this collapse model is shown in Figure 8, where the dynamics and growth of various modes are manifest. The movie begins at 1 ms before core bounce and ends at 71 ms at which point the simulation was terminated. Initially, Cartesian-grid–induced features are prevalent in the post-bounce dynamics of the proto-neutron star. These become gradually washed out as the star goes nonaxisymmetrically unstable and develops , , and features. The mode structure is nonstationary and exhibits strong dynamic variation with radius and time.
Ott et al.  find that the characteristic wave strain of this model peaks at a frequency of 930 Hz, which is the pattern frequency associated with its nonaxisymmetric unstable one-armed structure. Such a value is slightly larger than the typical frequencies of the burst signals from core collapse bounces ( 300–800 Hz). If confirmed when incorporating additional physics in the modelling, namely magnetic fields, and such dynamic, low nonaxisymmetric instabilities turn out to be a generic feature of the early phase of post-bounce dynamics, the associated gravitational waves might be detectable out to much larger distances than the Milky Way.
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. The inclusion of magnetic fields in those simulations has barely been initiated. In the following we discuss only multidimensional work, directing the interested reader to earlier versions of this review that include basic coverage of one-dimensional simulations.
Evans , building on previous work by Dykema , also studied the axisymmetric gravitational collapse problem for nonrotating 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 studies discussed so far have adopted spherical polar coordinates. Numerical instabilities are encountered at the origin () and at the polar axis (), where some fields diverge due to the coordinate singularities. Evans made important contributions toward regularizing the gravitational field equations in such situations . These coordinate problems have deterred researchers from building three-dimensional codes in spherical coordinates. In this case, Cartesian coordinates are adopted which, despite the desired property of being everywhere regular, present the important drawback of not being adapted to the topology of astrophysical systems. As a result, this has important consequences on the grid resolution requirements. To the best of our knowledge the only extensions of axisymmetric 2+1 codes in spherical coordinates to three dimensions has been accomplished by Stark  and by Dimmelmeier et al. . In the former case, no applications in relativistic astrophysics have yet been presented. For the latter, the CoCoNuT code is being applied to a number of problems such as core collapse (discussed in the preceding section), black-hole formation, and neutron star instabilities (see below).
Recently, Alcubierre et al.  proposed a method (cartoon) that does not suffer from stability problems at coordinate singularities and in which, in essence, Cartesian coordinates are used even for axisymmetric systems. Using this method, Shibata  investigated the effects of rotation on the criterion for prompt adiabatic collapse of rigidly and differentially rotating () polytropes to a black hole. Collapse of the initial approximate equilibrium models (computed by assuming a conformally-flat spatial metric) was induced by pressure reduction. In  it was shown that the criterion for black-hole formation depends strongly on the amount of angular momentum, but only weakly on its (initial) distribution. Shibata also studied the effects of shock heating using a gamma-law EOS, concluding that it is important in preventing prompt collapse to black holes in the case of large rotation rates. Using the same numerical approach, Shibata and Shapiro  have also 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 0.75. 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 black hole.
Also using the cartoon approach,  performed long-term simulations of the gravitational collapse of rapidly rotating neutron stars to black holes in axisymmetry, employing the code described in  (see Section 4.4). Together with the cartoon method the key novel ingredient added to the original code was an algorithm to excise the black-hole singularity when the black hole first forms, which allowed one to extend the evolutions far beyond what could be achieved without it. In agreement with the early works of [277, 386] it was found that unstable rotating stars collapse promptly to black holes only if they are “sub-Kerr”, namely if , being the angular momentum and the mass of the system. This result was also independently obtained by  in a parameter study of rotational collapses of supramassive neutron stars, employing the cartoon technique as well (but no black-hole excision). Otherwise, for “supra-Kerr” models,  found that the stars collapse to tori, which later fragment. These types of systems have also been studied in three dimensions as we discuss below.
Alternatively, axisymmetric codes employing the characteristic formulation of the Einstein equations  do not share the axis instability problems of most 2+1 codes. Such axisymmetric characteristic codes, conveniently coupled to hydrodynamics codes, are a promising way of studying the axisymmetric collapse problem and, in particular, of extracting accurate information regarding gravitational radiation. First steps in this direction have been taken by , where an axisymmetric Einstein-perfect fluid code is presented, and in  where this code was applied to study the core collapse of perturbed spherically-symmetric polytropes. 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.
Further recent improvements in the numerical modelling have accounted for physical effects, which may alter the structure of differentially rotating stars on secular timescales, such as the viscosity and magnetic fields. The former effect involves the solution of the relativistic Navier–Stokes equation, which has been attempted by , both in axisymmetry and in full 3D. It has been found that the role of viscosity in a hypermassive star generically leads to the formation of a compact, uniformly-rotating core, often unstable to gravitational collapse, surrounded by a low-density disk. In some models viscous heating becomes important enough to prevent the prompt collapse to a black hole, which is nevertheless the ultimate outcome once the star cools.
On the other hand, the collapse of magnetized neutron stars in full general relativity has also been numerically undertaken recently by [360, 105]. Both works study the fate of a hypermassive neutron star as it collapses to form a rotating black hole surrounded by a massive torus, a process which does not happen in prompt timescales due to the transport of angular momentum via magnetic braking and the MRI. Analysis of the lifetime and energetics of the black-hole–torus system formed in these computations suggests that the collapse of hypermassive neutron stars is a possible candidate for the central engine of short gamma-ray bursts.
In  the gravitational collapse of uniformly-rotating neutron stars to Kerr black holes was studied in three dimensions, also using the cactus/whisky code for numerical relativity. Long-term simulations were possible by excising the region of the computational domain that includes the curvature singularity when the black hole forms and lies within the apparent horizon. It was found that for sufficiently high angular velocities, the collapse can lead to the formation of a differentially-rotating unstable disk. Gravitational waveforms from such simulations were reported in . Despite good qualitative agreement being found with waveforms obtained in previous axisymmetric simulations , the newly computed amplitudes are about an order of magnitude smaller due to the more realistic rotation rates of the collapse progenitors.
More recently,  have succeeded in computing the gravitational waveforms from the gravitational collapse of a neutron star to a rotating black hole well beyond what had so far been possible, even providing information on the ring-down signal through which the formed black hole relaxes to its asymptotic, stationary state. These new three-dimensional computations, contrary to most approaches, did not excise the region inside the apparent horizon once the black hole formed in order to improve, in principle, the numerical stability of the code by removing the curvature singularity. Instead, they were performed without excision, which, combined with suitable singularity-avoiding gauge conditions and the use of minute (artificial) Kreiss–Oliger type numerical dissipation in the spacetime evolution and gauge equations, improved dramatically their long-term stability. An example of this remarkable result is provided by Figure 9, which shows a snapshot during the evolution of one of the collapse models of 1. The figure displays the gravitational waveform train extracted by matching the numerical spacetime with nonspherical perturbations of a Schwarzschild black hole described in terms of gauge-invariant odd and even-parity metric perturbations, namely the lowest-order odd-parity multipole . Moreover, in the central parts of the figure the collapse of the lapse function and the formation of a rotating black hole are visible.  shows that (stellar mass) black-hole formation events would produce gravitational radiation potentially detectable by the present generation of laser-interferometer detectors if the events took place at Galactic distances.
This work is licensed under a Creative Commons Attribution-Noncommercial-No Derivative Works 2.0 Germany License.