Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity

This article presents a comprehensive overview of numerical hydrodynamics and magneto-hydrodynamics (MHD) in general relativity. Some significant additions have been incorporated with respect to the previous two versions of this review (2000, 2003), most notably the coverage of general-relativistic MHD, a field in which remarkable activity and progress has occurred in the last few years. Correspondingly, the discussion of astrophysical simulations in general-relativistic hydrodynamics is enlarged to account for recent relevant advances, while those dealing with general-relativistic MHD are amply covered in this review for the first time. The basic outline of this article is nevertheless similar to its earlier versions, save for the addition of MHD-related issues throughout. Hence, different formulations of both the hydrodynamics and MHD equations are presented, with special mention of conservative and hyperbolic formulations well adapted to advanced numerical methods. A large sample of numerical approaches for solving such hyperbolic systems of equations is discussed, paying particular attention to solution procedures based on schemes exploiting the characteristic structure of the equations through linearized Riemann solvers. As previously stated, a comprehensive summary of astrophysical simulations in strong gravitational fields is also presented. These are detailed in three basic sections, namely gravitational collapse, black-hole accretion, and neutron-star evolutions; despite the boundaries, these sections may (and in fact do) overlap throughout the discussion. The material contained in these sections highlights the numerical challenges of various representative simulations. It also follows, to some extent, the chronological development of the field, concerning advances in the formulation of the gravitational field, hydrodynamics and MHD equations and the numerical methodology designed to solve them. To keep the length of this article reasonable, an effort has been made to focus on multidimensional studies, directing the interested reader to earlier versions of the review for discussions on one-dimensional works. Electronic Supplementary Material Supplementary material is available for this article at 10.12942/lrr-2008-7.


Introduction
The description of important areas of modern astronomy, such as high-energy astrophysics or gravitational wave astronomy, requires general relativity. Einstein's theory of gravitation plays a major role in astrophysical scenarios involving compact objects such as neutron stars and black holes. High-energy radiation is often emitted in regions of strong gravitational fields near such compact objects. The production of relativistic radio jets in active galactic nuclei, explained by either hydrodynamic or electromagnetic mechanisms, involves rotating supermassive black holes. The discovery of kHz quasi-periodic oscillations in low-mass X-ray binaries extended the frequency range over which these oscillations occur into timescales associated with the relativistic, innermost regions of accretion disks. A relativistic description is also necessary in scenarios involving explosive collapse of very massive stars (∼ 30 M ) to black holes (in the collapsar and hypernova models), or during the last phases of the coalescence and merge of neutron star binaries and neutron-star-black-hole binaries. These catastrophic events are believed to occur at the central engine of the most highly energetic events in nature, gamma-ray bursts (GRBs). Astronomers have long been scrutinizing these systems using the complete frequency range of the electromagnetic spectrum. Nowadays, they are the main targets of ground-based laser interferometers of gravitational radiation. The direct detection of these elusive ripples in the curvature of spacetime, and the wealth of new information that could be extracted from them, is one of the driving motivations of present-day research in relativistic astrophysics.
This field is currently in a state of very rapid development. On the one hand, steady progress is being driven by observational facilities, either ground based or in space, provided by a large number of high-energy X-ray and γ-ray telescopes such as MAGIC and HESS, and satellites such as ASCA (which stopped operating in 2001), RXTE, Chandra, XMM-Newton, INTEGRAL, Suzaku, and GLAST (whose launch is scheduled for the first half of 2008), along with the first generation of laser-interferometer gravitational-wave detectors, which have just become operational (LIGO, VIRGO, GEO, and TAMA). On the other hand, a second factor responsible for such development is provided by the corresponding theoretical studies, which struggle to come up with models to explain those observations. Modern theoretical astrophysics has always relied on numerical simulations as a formidable way to improve our understanding of the dynamics of astrophysical systems. Presently, the rapid increase in computing power through parallel supercomputers and the associated advance in software technologies is making possible large-scale numerical simulations in the framework of general relativity. However, computational astrophysicists and numerical relativists still face a daunting task. In a fairly general case, the equations governing the dynamics of relativistic astrophysical systems are an intricate, coupled system of time-dependent partial differential equations, comprising the general-relativistic hydrodynamic and magneto-hydrodynamic (MHD) equations and the Einstein gravitational field equations. Furthermore, the number of equations must be augmented in situations, which require one to account for nonadiabatic processes, such as viscosity, resistivity, radiative transfer, or the description of the microphysics of neutron-star interiors (realistic equations of state for nuclear matter, nuclear physics, and so on).
However, in some astrophysical situations of interest, such as in studies of accretion of matter onto compact objects, the propagation of relativistic outflows and jets, or oscillations of relativistic stars, the "test-fluid" approximation is good enough to provide an accurate description of the underlying dynamics. In this approximation, the self-gravity of the fluid is neglected in comparison to the background gravitational field. This is best exemplified in accretion problems, where the mass of the accreting fluid is usually much smaller than the mass of the compact object. Additionally, a description employing ideal hydrodynamics and MHD (i.e., with the stress-energy tensor being that of a perfect fluid, which is furthermore a perfect conductor), is also a fairly standard choice in numerical astrophysics.
The scope of this review is limited to discussing the first set of equations, namely hydrodynamics

Equations of General-Relativistic Hydrodynamics
The general-relativistic hydrodynamics equations consist of the local conservation laws of the stressenergy tensor, T µν (the Bianchi identities) and of the matter current density, J µ (the continuity equation): As usual, ∇ µ stands for the covariant derivative associated with the four-dimensional spacetime metric g µν . The density current is given by J µ = ρu µ , u µ representing the fluid 4-velocity and ρ the proper rest-mass density. The stress-energy tensor for a nonperfect (unmagnetized) fluid is defined as (see, e.g., [264]) where ε is the specific energy density of the fluid in its rest frame, p is the pressure, and h µν is the spatial projection tensor h µν = u µ u ν + g µν . In addition, η and ζ are the shear and bulk viscosities. The expansion θ, describing the divergence or convergence of the fluid world lines, is defined as θ = ∇ µ u µ . The symmetric, trace-free, spatial shear tensor σ µν is defined by and, finally, q µ is the energy flux vector. In the following, we will neglect nonadiabatic effects, such as viscosity and heat transfer, assuming the stress-energy tensor to be that of a perfect fluid T µν = ρhu µ u ν + pg µν , where we have introduced the relativistic specific enthalpy, h, defined by Introducing an explicit coordinate chart (x 0 , x i ), the previous conservation equations read where the scalar x 0 represents a foliation of the spacetime with hypersurfaces (coordinatized by x i ). Additionally, √ −g is the volume element associated with the 4-metric, with g = det(g µν ), and the objects Γ ν µλ are the 4-dimensional Christoffel symbols. In order to close the system, the equations of motion (1) and the continuity equation (2) must be supplemented with an equation of state (EOS) relating some fundamental thermodynamic quantities. In general, the EOS takes the form p = p(ρ, ε).
The available EOSs have become sophisticated enough to take into account physical and chemical processes such as molecular interactions, quantization, relativistic effects, nuclear physics, etc. Nevertheless, due to their simplicity, the most widely employed EOSs in numerical simulations Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 9 in astrophysics are the ideal fluid EOS, p = (Γ − 1)ρε, where Γ is the adiabatic index, and the polytropic EOS (e.g., to build equilibrium stellar models), p = Kρ Γ ≡ Kρ 1+1/N , K being the polytropic constant and N the polytropic index. State-of-the-art, microphysical EOSs that describe the interior of compact stars at nuclear matter densities have also been developed. While they are being increasingly used in the numerical modelling of relativistic stars, the true EOS of neutron-star interiors remains still largely unknown, as reproducing the associated densities and temperatures is not amenable to laboratory experimentation.
In the "test-fluid" approximation, where the fluid self-gravity is neglected, the dynamics of the system is completely governed by Equations (1) and (2), together with the EOS (9). In those situations in which such an approximation does not hold, the previous equations must be solved in conjunction with the Einstein gravitational-field equations, which describe the evolution of the geometry in a dynamic spacetime. A detailed description of the various numerical approaches to solving the Einstein equations is beyond the scope of the present article (see, e.g., [217,41] for recent reviews). We briefly mention that the Einstein equations, both in vacuum as well as in the presence of matter fields, can be formulated as an initial value (Cauchy) problem, using the 3+1 decomposition of spacetime [25]. More details can be found in, e.g., [431,158]. Given a choice of gauge, the Einstein equations in the 3+1 formalism split into evolution equations for the 3-metric γ ij and the extrinsic curvature K ij (the second fundamental form), and constraint equations (the Hamiltonian and momentum constraints), which must be satisfied on every time slice. The Arnowitt-Deser-Misner (ADM) 3+1 metric equations have been shown over the years to be intrinsically unstable for long-term numerical simulations, especially for those dealing with blackhole spacetimes. Recently, there have been diverse attempts to reformulate those equations into forms better suited for numerical investigations (see [363,39,5] and references therein). Among the various approaches proposed, the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) reformulation of the ADM system [363,39] is very appropriate for long-term stable numerical work. BSSN makes use of a conformal decomposition of the 3-metric,γ ij = e −4φ γ ij and the trace-free part of the extrinsic curvature, A ij = K ij − γ ij K/3, with the conformal factor φ chosen to satisfy e 4φ = γ 1/3 ≡ det(γ ij ) 1/3 . In this formulation, as shown by [39], in addition to the evolution equations for the conformal three-metricγ ij and the conformal-traceless extrinsic-curvature variablesÃ ij , there are evolution equations for the conformal factor φ, the trace of the extrinsic curvature K and the "conformal connection functions"Γ i =γ jkΓi jk . Indeed, BSSN (or slight modifications thereof) is currently the standard 3+1 formulation for most numerical relativity groups [5,158]. Long-term stable applications include strongly-gravitating systems such as neutron stars (both, isolated and in binary systems) and single and binary black holes. Such binary-black-hole evolutions, possibly the grandest challenge of numerical relativity ever, since the beginning of the field, have only been possible in the last few years (see, e.g., [323] and references therein).
Alternatively, a characteristic initial-value-problem formulation of the Einstein equations was developed in the 1960s by Bondi, van der Burg, and Metzner [59], and by Sachs [344]. This approach has gradually advanced to a state where long-term stable evolutions of caustic-free spacetimes in multiple dimensions are possible, even including matter fields (see [217] and references therein). A comprehensive review of the characteristic formulation is presented in a Living Reviews article by Winicour [419]. Examples of this formulation in general-relativistic hydrodynamics are discussed in various sections of the current article.
Traditionally, most of the approaches for numerical integrations of the general-relativistic hydrodynamics equations have adopted spacelike foliations of the spacetime, within the 3+1 formulation. More recently, however, covariant forms of these equations, well suited for advanced numerical methods, have also been developed. This is reviewed next (Section 2.1) in a chronological way.

Spacelike (3+1) approaches
In the 3+1 (ADM) formulation [25,158], spacetime is foliated into a set of nonintersecting spacelike hypersurfaces. There are two kinematic variables describing the evolution between these surfaces: the lapse function α, which describes the rate of advance of time along a timelike unit vector n µ normal to a surface, and the spacelike shift vector β i that describes the motion of coordinates within a surface.
The line element is written as where γ ij is the 3-metric induced on each spacelike slice.

1+1 Lagrangian formulation (May and White)
For historical reasons it is worth beginning the overview of the formulations with the pioneering numerical work in general relativistic hydrodynamics, namely the one-dimensional gravitational collapse code of May and White [247,248]. Building on theoretical work by Misner and Sharp [263], May and White developed a time-dependent numerical code to solve the evolution equations describing adiabatic spherical collapse in general relativity. This code was based on a Lagrangian finite-difference scheme (see Section 4.1), in which the coordinates are co-moving with the fluid. Artificial viscosity terms were included in the equations to damp the spurious numerical oscillations caused by the presence of shock waves in the flow solution. May and White's formulation became the starting point of a large number of numerical investigations in subsequent years and, hence, it is useful to describe its main features in some detail. For a spherically-symmetric spacetime, the line element can be written as ds 2 = −a 2 (m, t)dt 2 + b 2 (m, t)dm 2 + R 2 (m, t)(dθ 2 + sin 2 θdφ 2 ), (12) m being a radial (Lagrangian) coordinate, indicating the total rest mass enclosed inside the sphere 4 3 πR(m, t) 3 . The co-moving character of the coordinates leads, for a perfect fluid, to a stress-energy tensor whose nonvanishing components are T 1 1 = T 2 2 = T 3 3 = −p, T 0 0 = (1 + ε)ρ. In these coordinates the local conservation equation for the baryonic mass, Equation (2), can be easily integrated to yield the metric potential b: The gravitational field equations, Equation (10), and the equations of motion, Equation (1), reduce to the following quasi-linear system of partial differential equations (see also [263]): ∂ε ∂t + p ∂ ∂t 1 ρR 2 represents the total mass interior to radius m at time t. The final system, Equations (14) - (17), is closed with an EOS of the form given by Equation (9). Hydrodynamics codes based on the original formulation of May and White and on later versions (e.g., [406]) have been used in many nonlinear simulations of supernova and neutron-star collapse (see, e.g., [262,391] and references therein), as well as in perturbative computations of sphericallysymmetric gravitational collapse within the framework of the linearized Einstein equations [346,347]. However, the Lagrangian character of May and White's formulation, together with other theoretical considerations concerning the particular coordinate gauge, has prevented its extension to multiple-dimensional calculations. In spite of this, for one-dimensional problems, the Lagrangian approach adopted by May and White has considerable advantages with respect to an Eulerian approach with spatially-fixed coordinates, most notably the reduction of numerical diffusion.

3+1 Eulerian formulation (Wilson)
The use of Eulerian coordinates in multiple-dimensional numerical-relativistic hydrodynamics started with the pioneering work of Wilson [411,417]. Introducing the basic dynamic variables D, S µ , and E, representing the relativistic density, momenta, and generalized internal energy, respectively, and defined as the equations of motion in Wilson's formulation [411,413] are: with the "transport velocity" given by V µ = u µ /u 0 . We note that, in the original formulation [413], the momentum density equation, Equation (21), is solved only for the three spatial components, S i , and S 0 is obtained through the 4-velocity normalization condition u µ u µ = −1.
A direct inspection of the system shows that the equations are written as a coupled set of advection equations. In doing so, the terms containing derivatives (in space or time) of the pressure are treated as source terms. This approach, hence, sidesteps an important guideline for the formulation of nonlinear hyperbolic systems of equations, namely the preservation of their conservation form. This is a necessary condition to guarantee correct evolution in regions of sharp entropy generation (i.e., shocks). Furthermore, some amount of numerical dissipation must be used to stabilize the solution across discontinuities. In this spirit, the first attempt to solve the equations of general-relativistic hydrodynamics in the original Wilson's scheme [411] used a combination of finite-difference upwind techniques with artificial viscosity terms. Such terms adapted the classic treatment of shock waves introduced by von Neumann and Richtmyer [407] to the relativistic regime (see Section 4.1.1).
Wilson's formulation has been (and still is) widely used in hydrodynamic codes developed by a variety of research groups. (The latest extensions made to incorporate magnetic fields are discussed in Section 3.1.1). Many different astrophysical scenarios were first investigated with these codes, including axisymmetric stellar core collapse [279,277,283,38,386,319,118], accretion onto compact objects [170,317], numerical cosmology [72,73,17] and, more recently, the coalescence of neutron-star binaries [416,418,242] (see [417] for an up-to-date summary of such studies). This formalism has also been employed, in the special-relativistic limit, in numerical studies of heavy-ion collisions [415,249]. We note that in most of these investigations, the original formulation of the hydrodynamic equations was slightly modified by redefining the dynamic variables, Equation (19), with the addition of a multiplicative α factor (the lapse function) and the introduction of the Lorentz factor, W ≡ αu 0 : As mentioned before, the description of the evolution of self-gravitating matter fields in general relativity requires a joint integration of the hydrodynamic equations and the gravitational field equations (the Einstein equations). Using Wilson's formulation for the fluid dynamics, such coupled simulations were first considered in [413], building on a vacuum numerical-relativity code specifically developed to investigate the headon collision of two black holes [382]. The resulting code was axially symmetric and aimed to integrate the coupled set of equations in the context of stellar core collapse [120].
More recently, Wilson's formulation has been applied to the numerical study of the coalescence of neutron-star binaries in general relativity [416,418,242] (see Section 5.3.2). These studies adopt an approximation scheme for the gravitational field, by imposing the simplifying condition that the three-geometry (the 3-metric γ ij ) is conformally flat. The curvature of the three metric is then described by a position-dependent conformal factor φ 4 times a flat-space Kronecker delta, and the line element, Equation (11), reads Therefore, in this approximation scheme all radiation degrees of freedom are removed. Moreover, under the maximal-slicing condition (K = 0), the field equations reduce to a set of five Poisson-like elliptic equations in flat spacetime for the lapse, the shift vector, and the conformal factor. While in spherical symmetry this approach is no longer an approximation, being identical to Einstein's theory, beyond spherical symmetry its quality degrades. In [193] it was shown by means of numerical simulations of extremely relativistic disks of dust that it has the same accuracy as the first post-Newtonian approximation. In less extreme situations, however, as in rotational stellarcore collapse, the approximation yields results comparable to those obtained using full general relativity (see [74,92,365]).
Wilson's formulation showed some limitations in handling situations involving ultrarelativistic flows (W ≥ 2), as first pointed out by Centrella and Wilson [73]. Norman and Winkler [291] performed a comprehensive numerical assessment of such formulation by means of special-relativistic hydrodynamics simulations. Figure 1 reproduces a plot from [291] in which the relative error of the density compression ratio in the relativistic shock reflection problem -the heating of a cold gas, which impacts at relativistic speeds with a solid wall and bounces back -is displayed as a function of the Lorentz factor W of the incoming gas. The source of the data is [73]. This figure shows that for Lorentz factors of about 2 (v ≈ 0.86 c), the threshold of the ultrarelativistic limit, the relative errors are between 5% and 7% (depending on the adiabatic exponent of the gas), showing a linear growth with W .
Norman and Winkler [291] conclude that those large errors were mainly due to the way in which the artificial viscosity terms are included in the numerical scheme in Wilson's formulation. These terms, commonly called Q in the literature (see  Figure 1: Results for the shock heating test of a cold, relativistically inflowing gas against a wall using the explicit Eulerian techniques of Centrella and Wilson [73]. The plot shows the dependence of the relative errors of the density compression ratio versus the Lorentz factor W for two different values of the adiabatic index of the flow, Γ = 4/3 (triangles) and Γ = 5/3 (circles) gases. The relative error is measured with respect to the average value of the density over a region in the shocked material. The data are from Centrella and Wilson [73] and the plot reproduces a similar one from Norman and Winkler [291].
in some cases, namely at the pressure gradient in the source of the momentum equation, Equation (21), and at the divergence of the velocity in the source of the energy equation, Equation (22). However, [291] proposes that one add the Q terms in a consistent way, in order to consider the artificial viscosity as a real viscosity. Hence, the hydrodynamic equations should be rewritten for a modified stress-energy tensor of the following form: T µν = ρ(1 + ε + (p + Q)/ρ)u µ u ν + (p + Q)g µν . (25) In this way, for instance, the momentum equation takes the following form (in flat spacetime): In Wilson's original formulation, Q is omitted in the two terms containing the quantity ρh. In general, Q is a nonlinear function of the velocity and, hence, the quantity QW 2 V in the momentum density of Equation (26) is a highly nonlinear function of the velocity and its derivatives. This fact, together with the explicit presence of the Lorentz factor in the convective terms of the hydrodynamic equations, as well as the pressure in the specific enthalpy, make the relativistic equations more strongly coupled than their Newtonian counterparts. As a result, Norman and Winkler proposed the use of implicit schemes as a way to describe more accurately such coupling. Their code, which in addition incorporates an adaptive grid, reproduces very accurate results even for ultrarelativistic flows with Lorentz factors of about ten in one-dimensional, flat spacetime simulations.
Recently, Anninos and Fragile [19] have compared state-of-the-art artificial viscosity schemes and high-order nonoscillatory central schemes (see Section 4.1.3) using Wilson's formulation for the former class of schemes and a conservative formulation for the latter (similar to the one considered in [311,309]; see Section 2.2.2). Using the three-dimensional Cartesian code cosmos, these authors found that earlier results for artificial viscosity schemes in shock tube tests or shock reflection tests are not improved, i.e. the numerical solution becomes increasingly unstable for shock velocities greater than ∼ 0.95 c. On the other hand, results for the shock-reflection problem with a secondorder finite-difference central scheme show the suitability of such a scheme to handle ultrarelativistic flows (see Figure 8 of [19]), the underlying reason being the use of a conservative formulation of the hydrodynamic equations rather than the particular scheme employed (see Section 4.1.3).

3+1 conservative Eulerian formulation (Ibáñez and colleagues)
In 1991, Martí, Ibáñez, and Miralles [237] presented a new formulation of the (Eulerian) generalrelativistic hydrodynamics equations. This formulation was designed to take fundamental advantage of the hyperbolic and conservative character of the equations, contrary to the formulation discussed in the previous Section 2. 1.2. From the numerical point of view, the hyperbolic and conservative nature of the relativistic Euler equations allows for the use of schemes based on the characteristic fields of the system, bringing to relativistic hydrodynamics existing tools of computational classical fluid dynamics. This procedure departs from earlier approaches, most notably in avoiding the need for artificial dissipation terms to handle discontinuous solutions [411,413], as well as implicit schemes as proposed in [291].
If a numerical scheme written in conservation form converges, it automatically guarantees the correct Rankine-Hugoniot (jump) conditions across discontinuities -the shock-capturing property (see, e.g., [218]). Writing the relativistic hydrodynamic equations as a system of conservation laws, identifying the suitable vector of unknowns, and building up an approximate Riemann solver permitted the extension of state-of-the-art high-resolution shock-capturing schemes (HRSC in the following) from classical fluid dynamics into the realm of relativity [237].
Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Theoretical advances in the mathematical character of the relativistic hydrodynamic equations were first achieved studying the special relativistic limit. In Minkowski spacetime, the hyperbolic character of relativistic hydrodynamics and MHD is exhaustively studied by Anile and collaborators (see [15] and references therein) by applying Friedrichs' definition of hyperbolicity [143] to a quasilinear form of the system of hydrodynamic and MHD equations, where A µ are the Jacobian matrices of the system and w is a suitable set of primitive (physical) variables (see below). The system (27) is hyperbolic in the time direction defined by the vector field ξ with ξ µ ξ µ = −1, if the following two conditions hold: (i) det(A µ ξ µ ) = 0 and (ii) for any ζ such that ζ µ ξ µ = 0, ζ µ ζ µ = 1, the eigenvalue problem A µ (ζ µ − λξ µ )r = 0 has only real eigenvalues {λ n } n=1,···,5 and a complete set of right-eigenvectors {r n } n=1,···,5 . Besides verifying the hyperbolic character of the relativistic hydrodynamic equations, Anile and collaborators [15] obtained the explicit expressions for the eigenvalues and eigenvectors in the local rest frame, characterized by u µ = δ µ 0 . In Font et al. [134] those calculations were extended to an arbitrary reference frame in which the motion of the fluid is described by the 4-velocity u µ = W (1, v i ).
The approach followed in [134] for the equations of special-relativistic hydrodynamics was extended to general relativity in [36]. The choice of state vector (conserved quantities) in the 3+1 Eulerian formulation developed by Banyuls et al. [36] differs slightly from that of Wilson's formulation [411]. It comprises the rest-mass density (D), the momentum density in the j-direction (S j ), and the total energy density (E), measured by a family of observers, which is the natural extension (for a generic spacetime) of the Eulerian observers in classical fluid dynamics. Interested readers are directed to [36] for more complete definitions and geometrical foundations.
In terms of the primitive variables w = (ρ, v i , ε), the conserved quantities are written as: where the contravariant components v i = γ ij v j of the three-velocity are defined as and W is the relativistic Lorentz factor With this choice of variables the equations can be written in conservation form. Strict conservation is only possible in flat spacetime. For curved spacetimes there exist source terms, arising from the spacetime geometry. However, these terms do not contain derivatives of stress-energytensor components. More precisely, the resulting first-order flux-conservative hyperbolic system, well suited for numerical applications, reads: with g ≡ det(g µν ) satisfying √ −g = α √ γ with γ ≡ det(γ ij ). The state vector is given by Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 and the corresponding sources S(w) are The local characteristic structure of the previous system of equations was presented in [36]. The eigenvalues (characteristic speeds) of the corresponding Jacobian matrices are all real (but not distinct, one showing a threefold degeneracy as a result of the assumed directional splitting approach) and a complete set of right eigenvectors exists. More precisely, for a fluid moving along the x-direction, the eigenvalues read: where c s is the speed of sound. We note that the Minkowskian limit of these expressions is recovered properly (see [101]) as well as the Newtonian one (λ 0 = v x , λ ± = v x ± c s ). Hence, system (32) satisfies the definition of hyperbolicity. As will become apparent in Section 4.1.2, the knowledge of the spectral information of the flux-vector Jacobian matrices is essential in order to construct (upwind) HRSC schemes based on Riemann solvers. This information can be found in [36] (see also [137]). Three-dimensional, Eulerian codes, which evolve the coupled system of Einstein and hydrodynamics equations following the conservative Eulerian approach discussed in this section have been developed by Font et al. [137] and by Baiotti et al. [27] (see Section 4.4 for further details). In particular, in [137] the spectral decomposition (eigenvalues and right-eigenvectors) of the generalrelativistic hydrodynamics equations, valid for general spatial metrics, was derived, extending earlier results of [36] for nondiagonal metrics. A complete set of left-eigenvectors was presented by Ibáñez et al. [176]. Due to the paramount importance of the characteristic structure of the equations in the design of upwind HRSC schemes based upon Riemann solvers, we summarize all necessary information in Section 6.2 of the article.
The range of applications considered so far in relativistic astrophysics employing the above formulation of the hydrodynamics equations, Equations (32) - (35), includes the study of gravitational collapse (both stellar-core collapse to neutron stars and black-hole formation), accretion flows on to black holes, as well as neutron-star pulsations and neutron-star-binary mergers (see Section 5). The first applications in general relativity were performed, in one spatial dimension, in [237], using a slightly different form of the equations. Preliminary investigations of gravitational stellar collapse were attempted in [235,53] by coupling the above formulation of the hydrodynamic equations to a hyperbolic formulation of the Einstein equations developed by [54]. Evolutions of fully-dynamic spacetimes in the context of stellar-core collapse, both in spherical symmetry and in axisymmetry, have also been done [178,339,96,97,74]. These investigations are discussed in Section 5.1.1. In the special-relativistic limit this formulation has successfully been applied to simulate the evolution of (ultra-) relativistic extragalactic jets, using numerical models of increasing complexity (see, e.g., [241,7]).

Covariant approaches
General (covariant) conservative formulations of the general-relativistic hydrodynamics equations for ideal fluids, i.e., not restricted to spacelike foliations, have been reported in [117,311,309]. The form invariance of these approaches with respect to the nature of the spacetime foliation implies that existing work on highly-specialized techniques for fluid dynamics (i.e. HRSC schemes) can be adopted straightforwardly. In the next two Sections 2.2.1 and 2.2.2 we describe the existing covariant formulations in some detail.

Eulderink and Mellema
Eulderink and Mellema [115,117] were the first to derive a covariant formulation of the GRHD equations. As in the formulation discussed in Section 2.1.3, these authors took special care to use the conservative form of the system, with no derivatives of the dependent fluid variables appearing in the source terms. Additionally, this formulation is strongly adapted to a particular numerical method based upon a generalization of Roe's approximate Riemann solver. Such a solver was first applied to the nonrelativistic Euler equations in [337] and has been widely employed since in simulating compressible flows in computational fluid dynamics. Furthermore, their procedure is specialized for a perfect-fluid EOS, p = (Γ − 1)ρε, Γ being the (constant) adiabatic index of the fluid.
After the appropriate choice of the state-vector variables, the conservation laws, Equations (7) and (8), are rewritten in flux-conservative form. The flow variables are then expressed in terms of a parameter vector ω as where ω α ≡ Ku α , ω 4 ≡ K p ρh and K 2 ≡ √ −gρh = −g αβ ω α ω β . The vector F 0 represents the state vector (the unknowns), and each vector F i is the corresponding flux in the coordinate direction x i .
Eulderink and Mellema computed the exact "Roe matrix" [337] for the vector (38) and obtained the corresponding spectral decomposition. The characteristic information is used to solve the system numerically using Roe's generalized approximate Riemann solver. Roe's linearization can be expressed in terms of the average statẽ where indices L and R denote the left and right states in a Riemann problem (see Section 4.1.2). Further technical details can be found in [115,117]. The performance of this general-relativistic Roe solver was tested in a number of one-dimensional problems for which exact solutions are known, including nonrelativistic shock tubes, specialrelativistic shock tubes, and the spherical accretion of dust and a perfect fluid onto a (static) Schwarzschild black hole. In its special-relativistic version it has been used in the study of the confinement properties of relativistic jets [116]. However, no astrophysical applications in strong-field general-relativistic flows have yet been attempted with this formulation.

Papadopoulos and Font
In this formulation [311], the spatial velocity components of the 4-velocity, u i , together with the rest-frame density and internal energy, ρ and ε, provide a unique description of the state of the fluid at a given time and are taken as the primitive variables. They constitute a vector in Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 a five-dimensional space w = (ρ, u i , ε). The initial-value problem for Equations (7) and (8) is defined in terms of another vector in the same fluid-state space, namely the conserved variables, U, individually denoted (D, S i , E): Note that the state vector variables differ slightly from previous choices (see, e.g., Equations (19), (28), (29), (30), and (38)). With those definitions the GRHD equations take also the standard conservation law form, The flux vectors F j and the source terms S (which depend only on the metric, its derivatives, and the undifferentiated stress-energy tensor), are given by and The state of the fluid is uniquely described using either vector of variables, i.e., either U or w, and each one can be obtained from the other via the definitions (40) - (42) and the use of the normalization condition for the 4-velocity, g µν u µ u ν = −1. The local characteristic structure of the above system of equations was presented in [311], where the formulation proved well suited for the numerical implementation of HRSC schemes. The formulation of [311] was developed for a perfect fluid EOS. Extensions to account for generic EOS are given in [309], where a comprehensive analysis of general-relativistic hydrodynamics in conservation form is also provided.
A technical remark must be included here. In all conservative formulations discussed in Sections 2.1.3, 2.2.1, and 2.2.2, the time update of a given numerical algorithm is applied to the conserved quantities U. After this update the vector of primitive quantities w must be reevaluated, as those are needed in the Riemann solver (see Section 4.1.2). The relation between the two sets of variables is, in general, not in closed form and, hence, the recovery of the primitive variables is done using a root-finding procedure, typically a Newton-Raphson algorithm. This feature, distinctive of the equations of (special and) general-relativistic hydrodynamics -it does not exist in the Newtonian limit -may lead, in some cases, to accuracy losses in regions of low density and small speeds, apart from being computationally inefficient. The specific details of this issue for each formulation of the equations can be found in [36,117,311]. In particular, for the covariant formulation discussed in Section 2.2.1, there exists an analytic method to determine the primitive variables, which is, however, computationally very expensive since it involves many extra variables and the solving of a quartic polynomial. Therefore, iterative methods are still preferred [117]. On the other hand, we note that the covariant formulation discussed in this section, when applied to null-spacetime foliations, allows for a simple and explicit recovery of the primitive variables, as a consequence of the particular form of the Bondi-Sachs metric.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 2.2.2.1 Lightcone hydrodynamics A comprehensive numerical study of the formulation of the hydrodynamic equations discussed in this section was presented in [311], where it was applied to simulate one-dimensional relativistic flows on null (lightlike) spacetime foliations. The various demonstrations performed include standard shock-tube tests in Minkowski spacetime, perfect-fluid accretion onto a Schwarzschild black hole using ingoing null Eddington-Finkelstein coordinates, dynamic spacetime evolutions of relativistic polytropes (i.e., stellar models satisfying the Tolman-Oppenheimer-Volkoff equations of hydrostatic equilibrium) sliced along the radial null cones, and accretion of self-gravitating matter onto a central black hole.
Procedures for integrating various forms of the hydrodynamic equations on null hypersurfaces are much less common than on spacelike (3+1) hypersurfaces. They are presented in [182] (see [45,46] for more recent implementations). A Lagrangian method, limited to spherical symmetry, is developed in [258]. Recent work in [103] includes an Eulerian, nonconservative, formulation for general fluids in null hypersurfaces and spherical symmetry, including their matching to a spacelike section.
The general formalism laid out in [311,309] has been applied to astrophysical problems of increasing complexity. Applications in spherical symmetry include the investigation of accreting dynamic black holes, which can be found in [311,312]. Studies of the gravitational collapse of supermassive stars are discussed in [225] and studies of the interaction of scalar fields with relativistic stars are presented in [380]. Axisymmetric neutron-star spacetimes are considered in [378], and results for axisymmetric gravitational-core collapse using characteristic numerical relativity can be found in [379]. A proof-of-principle demonstration of the inclusion of matter fields in three dimensions was first given in [45] and short-time evolutions of a self-gravitating star in close orbit around a black hole have only been accomplished recently [46].

Going further: beyond ideal hydrodynamics
Formulations of the equations of nonideal hydrodynamics in general relativity are also available in the literature. The term "nonideal" is loosely employed here to refer to additional physical effects such as dissipative processes like viscosity and thermal conduction, and radiative transfer. These nonadiabatic effects can play a major role in some astrophysical systems, such as accretion disks or relativistic stars. Ideal magneto-hydrodynamics in general relativity is discussed in detail in Section 3.
In the discussion of the preceding section viscous effects have been explicitly neglected. Otherwise, Euler's equation would have to be replaced by the Navier-Stokes' equation (no longer hyperbolic) and the energy equation would contain additional terms accounting for heat transfer and the dissipation of kinetic energy. While their numerical description still represents a challenging task, viscous effects due to the transfer of momentum along velocity gradients by turbulence and thermal motions should, in general, be considered in some representative astrophysical scenarios, namely accretion disks for which shearing motions or steep velocity gradients may appear. In addition, viscosity is also believed to play a significant role in the stability and dynamics of compact stars, destroying differential rotation in rapidly-rotating neutron stars or suppressing the bar-mode instability as well as gravitational-radiation driven instabilities. The equations of viscous hydrodynamics, the Navier-Stokes-Fourier equations, have been formulated in relativity in terms of causal dissipative relativistic fluids (see the Living Reviews article by Müller [274] and references therein). These extended fluid theories, however, remain largely unexplored, numerically, in astrophysical systems. The reason may be the lack of an appropriate formulation that is well suited to numerical studies. Work in this direction was done by Peitz and Appl [315] who provided a 3+1 coordinate-free representation of different types of dissipative relativistic-fluid theories, which possess, in principle, the potential of being adapted well to numerical applications.
The interaction between matter and radiation fields, present in different levels of complexity Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 in all astrophysical systems, is described by the equations of radiation hydrodynamics. The Newtonian framework is highly developed (see, e.g., [257]; the special-relativistic transfer equation is also considered in that reference). Pons et al. [321] discuss a hyperbolic formulation of the radiative transfer equations, paying particular attention to the closure relations and to extending HRSC schemes to those equations. General-relativistic formulations of radiative transfer in curved spacetimes are considered in, e.g., [331], [432], and [314] (see also references therein). In [314] the equations of general-relativistic radiation hydrodynamics are derived from a tensor formalism and, therefore, can be applied to any spacetime or coordinate system. An explicit derivation of the radiation-moment equations is given for the case of the Schwarzschild spacetime.

Equations of General-Relativistic Ideal Magnetohydrodynamics
General-relativistic MHD is concerned with the dynamics of relativistic, electrically-conducting fluids (plasma) in the presence of magnetic fields. Here, we concentrate on purely ideal GRMHD, neglecting the presence of viscosity and heat conduction in the limit of infinite conductivity, i.e., the fluid is assumed to be a perfect conductor (see Section 3.3 for a brief discussion of the validity of this approximation). Like the GRHD equations discussed in the preceding Section 2, the GRMHD equations can also be cast in first-order, flux-conservative, hyperbolic form that is well suited to numerical work. Concerning this issue, the contribution of Anile [15] is remarkable (but see also [220] for an in-depth analysis of these equations), as he carried out a comprehensive study of the mathematical structure of the GRMHD equations. In order to analyze the hyperbolicity of the equations Anile found it convenient to write those equations using the following set of covariant variables V = (u µ , b µ , p, s), where b µ is the magnetic-field four-vector in the fluid rest frame and s is the specific entropy. With respect to these Anile variables, the system of GRMHD equations can be written as a quasi-linear system of the form A µA B ∇ µ V B = 0, where the indices A and B run from zero to nine, as the number of variables. The particular form of the 10 × 10 matrices A µ can be found in [15]. The eigenvalue structure of the system can be obtained by considering a generic characteristic hypersurface of the previous quasi-linear system, φ(x µ ) = 0, such that it defines a characteristic matrix A µ φ µ , whose determinant must vanish, with φ µ ≡ ∂φ/∂x µ . If we now consider a wave propagating in an arbitrary direction x with a speed λ, the normal to the characteristic hypersurface is given by the (spacelike) 4-vector φ µ = (−λ, 1, 0, 0), which can be substituted in the determinant of the characteristic matrix to obtain the corresponding characteristic polynomial, whose zeroes give the characteristic speed of the waves propagating in the x-direction. Three different kinds of waves can be obtained according to which factor in the equation resulting from the condition det(A µ φ µ ) = 0 vanishes, either entropic waves, Alfvén waves, or magnetosonic waves. It is worth noting that in Anile's study [15] both the entropy and Alfvén waves appear as double roots of the characteristic polynomial, a direct result of working with an augmented system of equations. Therefore, those nonphysical waves must be removed from the wave decomposition when building up a numerical scheme based upon such wave structure to solve the GRMHD equations.
In recent years there has been intense research on formulating and solving numerically the MHD equations in general-relativistic spacetimes, either background or dynamic [196,86,40,149,20,201,108,366,24,286,265,91]. Both, artificial viscosity and HRSC schemes have been developed and most of the astrophysical applications previously attempted with no magnetic fields included are currently being revisited, namely for studies of gravitational collapse, neutron-star dynamics, black-hole accretion, and jet formation. We note, however, that all HRSC numerical approaches employed, with the exception of [201,24], are based on central schemes (or incomplete Riemann solvers), for which the use of the full wave decomposition is not needed (see below).
In terms of the (Faraday) electromagnetic tensor F µν , Maxwell's equations read where

Numerical approaches
Maxwell's equations can be simplified if the fluid is a perfect conductor. In this case σ is infinite and, to keep the current finite, the term F µν u ν must become zero, which results in a vanishing electricfield four-vector in the fluid rest frame, e µ = 0. This case corresponds to the ideal MHD condition. Under this assumption, the electric field measured by the Eulerian observer has components and Maxwell's equations ∇ ν * F µν = 0 reduce to the divergence-free condition plus the induction equation for the evolution of the magnetic field: In the induction equation we use the notationṽ i = v i − β i /α. For a fluid endowed with a magnetic field, the stress-energy tensor is the sum of that of the fluid and that of the electromagnetic field, T µν = T µν Fluid + T µν EM , where T µν Fluid is given by Equation (5) for a perfect fluid. On the other hand T µν EM can be obtained from the Faraday tensor as follows: which, in ideal MHD, can be rewritten as where b 2 = b ν b ν . The total stress-energy tensor is thus given by with the definitions p * = p + b 2 /2 and h * = h + b 2 /ρ.

Nonconservative formulations
As happens for the equations of general-relativistic hydrodynamics discussed previously, available approaches to cast the corresponding magnetohydrodynamics equations in forms suitable for numerical work also follow to a great extent the basic distinction between conservative and nonconservative formulations. The latter, which are covered in this section, were first laid out by Wilson in the mid 1970s, who pioneered a new emerging field in numerical astrophysics, just as he had done a few years earlier for numerical Eulerian hydrodynamics. The original approach was followed by [429] and more recently by [86] and [20], who adopted an internal energy formulation with artificial viscosity for shock capturing and a method of characteristics for handling the magnetic fields. The dynamic variables used in the formulation of [86], D, S µ and E, coincide with those reported in Equation (23), save for the definition of the total four-momentum S µ , which includes the norm of the magnetic field four-vector b µ , namely S µ = (ρh + b 2 )W u µ , with W = αu 0 . With this choice, the "conservation" equations of mass density, momentum and energy are again written as a set of advection equations similar to Equations (20) - (22), except for the inclusion of the magnetic field terms, namely: where V i is again the transport velocity defined in Section 2.1.2. As in the purely hydrodynamic case, the time component of the momentum is obtained from the corresponding normalization condition together with the algebraic relations between the magnetic-field components. Similarly, an artificial viscosity must be added for the numerical solution of the last two equations to increase the entropy of the fluid at shocks. In [86] such a viscosity prescription coincides with that used in the unmagnetized case [171], except for the modification of the definition of the inertia to include the magnetic energy.
In addition, a supplementary evolution equation is needed for the constrained transport mag- , the induction equation, which reads: where B i must be subject to the divergence-free constraint, ∂B j ∂x j = 0. A similar internal-energy formulation was presented in [20]. There are, however, important subtle differences with respect to the approach of [86] in the particular form of the evolution equations, due to the different definition of the Lorentz (relativistic boost) factor employed by [20], namely W = √ −gu 0 = α √ γu 0 . In addition, there are a couple of points particular to the formulation of [20], which are worth mentioning. First, in the momentum-evolution equation, the term ∂ ∂x i ( which has significant numerical advantages for achieving stable formulations of sheared Alfvén waves. Second, the induction equation is rewritten in the form where η is a coefficient related to the largest characteristic speed in the flow multiplying the divergence cleanser function. Such a coefficient is included for purely numerical reasons in order to reduce ∇ · B discretization errors in time, either propagating any nonzero divergence off the computational domain or damping it (see [35] for tests of divergence-cleaning methods in classical MHD). This approach is adopted by [20] in order to avoid the implementation of the constrained transport method in their code (see Section 4.3) which, due to the need of an additional grid staggered with respect to the base grid, would much complicate the adaptive-mesh refinement (AMR) module of the code. Similar strategies are adopted in the AMR GRMHD code of [286].

Conservative formulations
We turn now to the discussion of conservative formulations for the GRMHD equations. Very recently, a number of groups have proposed various such conservative approaches to formulating these Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 equations, all aimed at solving them using methods specifically designed for hyperbolic systems of conservation (or balance) laws [149,20,201,108,366,24,286,265,91]. The existing approaches all share the key ingredient of casting the GRMHD equations in first-order flux-conservative form, yet the choice of variables (state-vector, fluxes, etc) differs (albeit slightly) among most of them. Therefore, for the sake of clarity, we start here by discussing one such conservative formulation first, namely that of [24] (which coincides with the one obtained by [286] and is also employed in the codes of [265,151]), leaving to the second part of this section the analysis of similarities and differences present in the additional formulations.
Following [24] the conservation equations for the energy-momentum tensor, together with the continuity equation and the equation for the evolution of the magnetic field, can be written as a first-order, flux-conservative, hyperbolic system equivalent to Equation (32). The state vector and the vector of fluxes of the GRMHD system of equations read where, again,ṽ i = v i − β i /α. Here, the conserved quantities for the continuity, Euler, and energy equations are now respectively defined as D = ρW , The corresponding vector of sources coincides with the one given by Equation (33) save for the use of the complete (fluid plus electromagnetic field) stress-energy tensor (the magnetic-field evolution equation is source-free). We note that expressions (57) and (58) contain components of the magnetic-field four-vector in the co-moving frame and of the magnetic-field three-vector measured by the Eulerian observer. The two are related by The hyperbolic structure of the flux-vector Jacobian matrices corresponding to Equations (57)-(58), the knowledge of which is essential for building up a numerical code based on upwind HRSC schemes, is analyzed in [24] (see also [23]). In the classical MHD case, the wave structure has been analyzed by [61], which shows that there are seven physical waves: two Alfvén waves (with eigenvalues λ a± = v x ± v a , v x and v a being the fluid and Alfvén speeds, respectively), two fast and two slow magnetosonic waves (λ f± = v x ± v f , λ s± = v x ± v s ), and one entropy wave (λ e = v x ). The expressions for the Alfvén and magnetosonic speeds read We note in passing the numerical analysis of the hyperbolic structure of the classical ideal MHD system carried out by [400], where nonuniform convergence to the true solution was found for a number of finite volume methods in coarse grid simulations.
The wave structure of the above GRMHD hyperbolic system is, understandably, more involved than that corresponding to the GRHD case. As mentioned before, Anile [15] (see also [16,404]) studied the characteristic structure of these equations (eigenvalues, right and left-eigenvectors) in the (ten-dimensional) space of covariant variables (u µ , b µ , p, s), s being the specific entropy. It was Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 25 found that only the entropic waves and the Alfvén waves can be explicitly written in closed form. For the formulation of [24] discussed in this section, and along each coordinate direction, they are given by On the other hand, the (fast and slow) magnetosonic waves are given by the numerical solution (through a root-finding procedure) of a quartic characteristic equation. Among the magnetosonic waves, the two solutions with maximum and minimum speeds are called fast magnetosonic waves λ f ± , and the two solutions in between the slow magnetosonic waves λ s± . The seven waves can be ordered as follows which defines the domain of dependence of the GRMHD system. We note that, as a result of working with the unconstrained, 10 × 10 system of equations, there appear superfluous eigenvalues associated with unphysical waves, which need to be removed (the entropy waves as well as the Alfvén waves appear as double roots). Any attempt to develop a numerical procedure based on the wave structure of the RMHD equations must, hence, remove these nonphysical waves (and the corresponding eigenvectors) from the wave decomposition. In the particular case of special-relativistic MHD, [199] and [198] eliminate the unphysical eigenvectors by enforcing the waves to preserve the values of the invariants u µ u µ = −1 and u µ b µ = 0, as suggested by [15]. Correspondingly, in [34] the physical eigenvectors are selected by comparing with the equivalent expressions in the nonrelativistic limit.
We also note that, as in the classical case, the relativistic MHD equations have degenerate states in which two or more wavespeeds coincide. This breaks the strict hyperbolicity of the system and may affect numerical schemes for its solution. These degeneracies are analyzed in [199], which finds that, with respect to the fluid rest frame, the degeneracies in both classical and relativistic MHD are the same; either the slow and Alfvén waves have the same speed as the entropy wave when propagating perpendicularly to the magnetic field (Degeneracy I), or the slow or the fast wave (or both) have the same speed as the Alfvén wave when propagating in a direction aligned with the magnetic field (Degeneracy II). On the other hand, [24] derives a covariant characterization of such degenerate states, which can be checked in any reference frame. In addition, [23] also works out a single set of right and left eigenvectors, which are regular and span a complete basis in any physical state, including degenerate states. Such renormalization procedure can be seen as a relativistic generalization of the work performed in [61] in classical MHD.
Finally, it is worth pointing out that major advances in the physical understanding of the wave structure of the GRMHD equations have been possible in recent years owing to the remarkable derivation of the exact solution of the Riemann problem in special relativistic MHD [340,150].
3.1.2.1 Additional conservative approaches As mentioned previously, a number of conservative formulations for the GRMHD equations have been proposed recently. Most of them show some differences with respect to the formulation just discussed, and some similarities too. Let us now consider these issues.
First, in the formulations of Gammie et al. [149] and Komissarov [201], the choice of conserved variables corresponding to the momentum and energy equations follows directly from the conservation of energy momentum, which is written with the free index down, i.e., ∇ µ T µ ν = 0. The Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 continuity equation, however, remains identical to the expression used in [24]. Therefore, the vector of conserved variables in [201,149] reads The choice of T 0 µ for the energy-momentum equations is not capricious, but rather motivated by the reduction on the number of source terms in the corresponding evolution equations for the case of spacetime metrics with symmetries (ignorable coordinates x µ ; see also [309] for a similar approach in the purely hydrodynamic case). This is not the situation in the formulation of [24] in which the energy-momentum equations are obtained, in the spirit of the 3+1 splitting, from projections of the conservation law, namely γ ν i ∇ µ T µ ν = 0 and n ν ∇ µ T µ ν = 0. Regarding the induction equation, both [149] and [201] use the following expression which differs from the expression used by [24] both from the use ofṽ i and the presence of √ γ and √ −g accompanying the time and spatial derivatives, respectively, in the latter case.
On the other hand, in the approach followed by Shibata and Sekiguchi [366], the evolution equations for the energy-momentum equations are also derived as in [24], from projections of the conservation law. However, the choice of conserved variables is slightly different, mainly due to the presence in all equations of a common factor e 6φ , φ being a spacetime conformal factor, motivated by the formulation used in the numerical code for solving Einstein's equations, namely BSSN [363,39]. (We note that such factors are also present in the formulation used by Shibata in the purely hydrodynamic case; see, e.g., [355].) Another important difference with [24] is the choice of the transport velocity v i = u i /u 0 throughout (instead ofṽ i ) and the fact that in the final form of the energy equation the continuity equation has not been subtracted. Additionally, the induction equation also adopts a different form as it is written for a three-magnetic field given by Correspondingly, in the conservative formulation of Anninos et al. [20] the main differences appear in the energy-momentum equations, leaving the continuity equation unaltered with respect to the nonconservative formulation discussed in Section 3.1.1. Regarding the momentum and energy-conservation equations the choice of variables in [20] is different to those of [24] and [366], as they are components of the energy-momentum tensor and not projections of the conservation law. But the choice is also different to that in [149] and [201], as the indices of the energymomentum tensor are both chosen to be contravariant (instead of T 0 0 and T 0 i ). More precisely, the corresponding conserved variables are S j = √ −gT 0j and ε = √ −gT 00 . Another important point to mention regarding the set of equations of this conservative formulation (and also in their nonconservative internal-energy formulation) is the presence of factors of 4π in some terms, which have not been absorbed by the definition of the electromagnetic tensor (as in most of the other approaches). Finally, regarding the induction equation, we note that it is also cast in fully conservative form and η the same coefficient as in Equation (56).
Next, the formulation of Duez et al. [108] also shares some similarities and has some differences with all previous approaches discussed so far. The transport velocity is used throughout.

27
The choice of conserved variable for the continuity equation, which the authors name ρ , is the same as in [20], being the equation also cast as a source-free advection equation. As for the momentum equation the choice of conserved variable coincides with that of [149,201], i.e., direct components of the energy-momentum tensor, namely √ −gT 0 i . However, for the remaining part of the energy-momentum conservation law, which leads to the energy equation, Duez et al. [108] do not follow [149,201] and, instead of using √ −gT 0 0 , they choose τ , the same variable as [24], also subtracting the continuity equation, namely τ = √ γn µ n ν T µν − ρ . Finally, the induction equation is written as in Equation (67), but without the divergence cleanser coefficient (η = 0). Finally, in the formulation of Del Zanna et al. [91] an effort is made in the choice of the state vector so that the GRMHD equations keep a somewhat close relation to classical MHD and to simplify the expression of the source terms. As in Antón et al. [24] the same definition for the three-velocity is used, along with the same variable for the relativistic density D (densitized with a factor γ 1/2 ). For the Euler equation, components of the stress-energy tensor T 0 j are used, as in the approaches of [149,201,108], which reduces the complexity of the corresponding source term. The most important difference with respect to other formulations is in the choice of conserved variable for the energy equation, for which [91] take √ −gT µν n ν . The continuity equation is not subtracted from the energy equation. As far as the induction equation is concerned, the same conservative expression used by [149,201] is employed.

Recovery of primitive variables
On the other hand, as for the case of the GRHD equations discussed before, iterative (root-finding) algorithms are also required for the GRMHD equations to recover the primitive variables w from the state vector U. Not surprisingly, the recovery procedure is in the GRMHD case more involved than for unmagnetized flows. A number of approaches are available in the literature. In particular, we direct the interested reader to [289] for a comparison of six different methods for performing such inversions from conserved to primitives variables for the case in which the matter thermodynamics is described by an ideal gas EOS. In [289] not only the accuracy of each of the methods is assessed but also its run-time.
For illustrative purposes, we discuss next with some detail the specific approach followed by [24], which is an extension to full general relativity of that developed by [199] in the special relativistic case. The procedure relies on the fact that it is not necessary to solve the system given by the definition of the conserved variables in terms of the primitives (see Equation (57)) for the three components of the momentum, but instead for its modulus S 2 = S i S i . By eliminating the components of b α through Equations (59) and after some algebra, it is possible to write S 2 as where Z ≡ ρhW 2 . In addition, the equation for the total energy can be worked out in a similar way Equations (68) and (69), together with the definitions of D and Z, form a system for the unknowns ρ, p and W , assuming the function h = h(ρ, p) is provided. The roots of this system can be found using a two-dimensional scheme (e.g. bisection or Newton-Raphson).

Going further
Like any approximation, ideal MHD has its own range of validity. This is mainly set by the relevant length and time scales of the system under study. The former needs to be typically much Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 larger than the Larmor radius of the plasma particles, while the latter needs to be much shorter than the time scale of the collisions among the plasma particles. Furthermore, it is assumed that such collisions yield no resistivity, which permits one to neglect any diffusion of the magnetic field lines. There are nevertheless examples in astrophysics in which resistivity effects may become important and ideal MHD may not be valid to describe the dynamics of the plasma, typically in situations involving dilute accretion flows. For example, the Hall effect in protostellar disks is one of those effects, with significant implications regarding the magneto-rotational instability, the mechanism believed to regulate angular momentum transport in disks through sustained accretion. In addition, reconnection of magnetic field lines (in, e.g., solar flares) is another important nonideal MHD effect, which may appear in turbulent plasmas near the scale of energy dissipation. When the fluid is no longer a perfect conductor, kinetic effects become important and the particle-in-cell (or Monte Carlo) methods are usually employed instead of finite-volume methods.
Magnetic reconnection and dissipation processes are in particular recognized to be very important in high-energy astrophysics. Such processes not only affect the global dynamics of the plasma in scenarios involving compact objects, but are also believed to account for the production of the high-energy emission. Magnetic-field-line reconnection has actually been already observed in the existing numerical simulations based upon ideal GRMHD. However, this can only be explained with purely numerical reasons due to the nonvanishing amount of numerical artificial resistivity inherent to every computer code. Only very recently have a few approaches been put forward to develop numerical schemes to handle relativistic plasma with physical resistivity, involving the formulation and solution of the resistive relativistic MHD equations [409,204].
Another interesting area in which new results are likely to be found may come from the coupling of GRMHD codes with particle codes. While the former type of codes give a good description of the large-scale dynamic processes, the latter are particularly suited to capturing the radiative processes operative at a microphysical scale in highly-relativistic objects, whose correct modelling may help extract and interpret observational diagnostics. Monte Carlo methods are extremely powerful in such aspects, as recently shown in the context of GRBs by [141], who confirmed the generation of strong electromagnetic fields by a Weibel-like instability in collisionless shocks and demonstrated their macroscopic propagation. When suitably coupled to GRMHD codes, particle codes could use the motivated initial and boundary conditions supplied by the magneto-fluid model to perform local particle-in-cell simulations.
Finally, it is worth noting that in some situations it may suffice to solve the equations of force-free electrodynamics as an alternative to solving the full set of GRMHD equations. In an astrophysical context this approach, valid in the low inertia limit of MHD, has been considered by, e.g., [203,250] to model neutron star and black-hole magnetospheres. On the one hand, the force-free equations of motion can also be cast as a set of conservation laws (amenable to be solved using the same techniques as for the full GRMHD system) and may be easier to integrate in regions where the rest-mass energy is much smaller than the magnetic energy density. On the other hand, if the inertia of the plasma particles is not negligible, the force-free approximation is inaccurate as it cannot capture the thermal energy evolution of the particles as well as their motion along field lines. Finally, the question of the likelihood of this approach to hold within current sheets, due to the tearing mode instability of magnetic field lines reconnection, should also be investigated.

Numerical Schemes
We turn now to describing the numerical schemes, mainly those based on finite differences, specifically designed to solve nonlinear hyperbolic systems of conservation laws. As discussed in the previous two sections, both the equations of general-relativistic hydrodynamics and magnetohydrodynamics fall in this category, irrespective of the formulation employed. Even though we also consider schemes based on artificial viscosity techniques, the emphasis is on the high-resolution shock-capturing (HRSC) schemes (or Godunov-type methods), based on (either exact or approximate) solutions of local Riemann problems using the characteristic structure of the equations. Such finite-difference schemes (or, in general, finite-volume schemes) have been the subject of diverse review articles and textbooks (see, e.g., [218,219,398,177]). For this reason only the most relevant features will be covered here, directing the reader to the appropriate literature for further details. In particular, an excellent introduction to the implementation of HRSC schemes in specialrelativistic hydrodynamics is presented in the Living Reviews article by Martí and Müller [240]. Alternative techniques to finite differences, such as smoothed particle hydrodynamics, (pseudo-) spectral methods and others, are briefly considered last.

Finite difference schemes
Any system of equations presented in the previous sections can be solved numerically by replacing the partial derivatives by finite differences on a discrete numerical grid and then advancing the solution in time via some time-marching algorithm. Hence, specification of the state vector U on an initial hypersurface, together with a suitable choice of EOS, followed by a recovery of the primitive variables, leads to the computation of the fluxes and source terms. Through this procedure the first time derivative of the data is obtained, which then leads to the formal propagation of the solution forward in time, with a timestep constrained by the Courant-Friedrichs-Lewy (CFL) condition.
The hydrodynamic and MHD equations (either in Newtonian physics or in general relativity) constitute nonlinear hyperbolic systems and, hence, smooth initial data can transform into discontinuous data (cross of characteristics in the case of shocks) in a finite time during the evolution. As a consequence, conventional finite-difference schemes (see, e.g., [218,219,398]) present important deficiencies when dealing with such systems. Typically, first-order accurate schemes are much too dissipative across discontinuities (excessive smearing) and second-order (or higher) schemes produce spurious oscillations near discontinuities, which do not disappear under grid refinement. To avoid these effects, standard finite-difference schemes have been conveniently modified in various ways to ensure high-order, oscillation-free accurate representations of discontinuous solutions, as we discuss next.

Artificial viscosity approach
The idea of modifying the hydrodynamic equations by introducing artificial viscosity terms to damp the amplitude of spurious oscillations near discontinuities was originally proposed by von Neumann and Richtmyer [407] in the context of the (classical) Euler equations. The basic idea is to introduce a purely artificial dissipative mechanism whose form and strength are such that the shock transition becomes smooth, extending over a small number of intervals ∆x of the space variable. In their original work von Neumann and Richtmyer proposed the following expression for the viscosity term: ∂x , v being the fluid velocity, ρ the density, ∆x the spatial interval, and k Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 a constant parameter whose value needs to be adjusted in every numerical experiment. This parameter controls the number of zones in which shock waves are spread.
This type of generic recipe, with minor modifications, has been used in conjunction with standard finite-difference schemes in all numerical simulations employing May and White's formulation, mostly in the context of gravitational collapse, as well as Wilson's formulation. So, for example, in May and White's original code [247] the artificial viscosity term, obtained in analogy with the one proposed by von Neumann and Richtmyer [407], is introduced in the equations accompanying the pressure, in the form: Further examples of similar expressions for the artificial viscosity terms, in the context of Wilson formulation, can be found in, e.g., [411,171,417]. In particular, a state-of-the-art formulation of the artificial viscosity approach is reported in [19]. Correspondingly, the interested reader is also directed to the recent work by [86,20,108] for details on diverse modern implementations of artificial viscosity terms in nonconservative formulations of the GRMHD equations.
The main advantage of the artificial viscosity approach is its simplicity, which results in high computational efficiency. Experience has shown, however, that this procedure is both problem dependent and inaccurate for ultrarelativistic flows [291,19]. Furthermore, the artificial viscosity approach has the inherent ambiguity of finding the appropriate form for Q that introduces the necessary amount of dissipation to reduce the spurious oscillations and, at the same time, avoids introducing excessive smearing at discontinuities. In many instances both properties are difficult to achieve simultaneously. A comprehensive numerical study of artificial-viscosity-induced errors in strong shock calculations in Newtonian hydrodynamics (including also proposed improvements) was presented by Noh [290].

High-resolution shock-capturing (HRSC) upwind schemes
In finite-difference schemes, convergence properties under grid refinement must be enforced to ensure that the numerical results are correct (i.e., if a scheme with an order of accuracy α is used, the global error of the numerical solution has to tend to zero as O(∆x) α as the cell width ∆x tends to zero). For hyperbolic systems of conservation laws, schemes written in conservation form are preferred since, according to the Lax-Wendroff theorem [216], they guarantee that the convergence, if it exists, is to one of the weak solutions of the original system of equations. Such weak solutions are generalized solutions that satisfy the integral form of the conservation system. They are C 1 classical solutions (continuous and differentiable) away from discontinuities and have a finite number of discontinuities.
For the sake of simplicity, let us consider in the following an initial value problem for a onedimensional scalar hyperbolic conservation law, and let us introduce a discrete numerical grid of spacetime points (x j , t n ). An explicit algorithm written in conservation form updates the solution from time t n to the next time level t n+1 as: wheref is a consistent numerical flux function (i.e.,f (u, u, · · · , u) = f (u)) of p + q + 1 arguments and ∆t and ∆x are the timestep and cell width respectively. Furthermore, u n j is an approximation of the average of u(x, t) within the numerical cell Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 The class of all weak solutions is too wide in the sense that there is no uniqueness for the initial value problem. The numerical method should, hence, guarantee convergence to the physically admissible solution. This is the vanishing-viscosity-limit solution, that is, the solution when η → 0, of the "viscous version" of the scalar conservation law, Equation (70): Mathematically, the solution one needs to search for is characterized by the entropy condition (in the language of fluids, the condition that the entropy of any fluid element should increase when crossing a discontinuity). The characterization of these entropy-satisfying solutions for scalar equations was given by Oleinik [301]. For hyperbolic systems of conservation laws it was developed by Lax [215].
The Lax-Wendroff theorem [216] cited above does not establish whether the method converges. To guarantee convergence, some form of stability is required, as Lax first proposed for linear problems (Lax equivalence theorem; see, e.g., [336]). Along this direction, the notion of totalvariation stability has proven very successful, although powerful results have only been obtained for scalar conservation laws. The total variation of a solution at time t = t n , TV(u n ), is defined as A numerical scheme is said to be TV-stable if TV(u n ) is bounded for all ∆t at any time for each initial data. In the case of nonlinear, scalar conservation laws it can be proven that TV-stability is a sufficient condition for convergence [218], as long as the numerical schemes are written in conservation form and have consistent numerical flux functions. Current research has focused on the development of high-resolution numerical schemes in conservation form satisfying the condition of TV-stability, such as the total variation diminishing (TVD) schemes [162] (see below).
Let us now consider the specific system of hydrodynamic equations as formulated in Equation (32) and let us consider a single computational cell of our discrete spacetime. Let Ω be a region (simply connected) of a given four-dimensional manifold M, bounded by a closed three-dimensional surface ∂Ω. We further take the three-surface ∂Ω as the standard-oriented hyper-parallelepiped made up of two spacelike surfaces {Σ x 0 , Σ x 0 +∆x 0 } plus timelike surfaces {Σ x i , Σ x i +∆x i } that join the two temporal slices together. By integrating system (32) over a domain Ω of a given spacetime, the variation in time of the state vector U within Ω is given -keeping apart the source terms -by the fluxes F i through the boundary ∂Ω. The integral form of system (32) is which can be written in the following conservation form, well adapted to numerical applications: Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 where ∆V = Besides its convergence properties, a numerical scheme written in conservation form ensures that, in the absence of sources, the (physically) conserved quantities, according to the partial differential equations, are numerically conserved by the finite difference equations.
The computation of the time integrals of the interface fluxes appearing in Equation (76) is one of the distinctive features of upwind HRSC schemes. One needs first to define the numerical fluxes, which are recognized as approximations to the time-averaged fluxes across the cell interfaces, which depend on the solution at those interfaces, U(x j + ∆x j /2, x 0 ), during a timestep, At the cell interfaces the flow can be discontinuous and, following the seminal idea of Godunov [155], the numerical fluxes can be obtained by solving a collection of local Riemann problems, as is depicted in Figure 2. This is the approach followed by the Godunov-type methods [164,112]. Figure 2 shows how the continuous solution is locally averaged on the numerical grid, a process that leads to the appearance of discontinuities at the cell interfaces. Physically, every discontinuity decays into three elementary waves of any of the following type: shock waves, rarefaction waves, and contact discontinuities (see, e.g., [398]). The complete structure of the Riemann problem can be solved analytically (see [155] for the solution in Newtonian hydrodynamics, [238,322] in specialrelativistic hydrodynamics, and [340,150] in special-relativistic MHD) and, accordingly, used to update the solution forward in time.
For reasons of numerical efficiency and, particularly in multiple dimensions, the exact solution of the Riemann problem is frequently avoided and linearized (approximate) Riemann solvers are preferred. These solvers are based on the exact solution of Riemann problems corresponding to a linearized version of the original system of equations. The spectral decomposition of the fluxvector Jacobian matrices is on the basis of all solvers (extending ideas used for linear hyperbolic systems). After extensive experimentation, the results achieved with approximate Riemann solvers are comparable to those obtained with the exact solver (see [398] for a comprehensive overview of Godunov-type methods, and [240] for an excellent summary of approximate Riemann solvers in special-relativistic hydrodynamics).
Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 continuous solution discrete solution contact discontinuity shock rarefaction , a local Riemann problem is set up as a result of the discretization process (bottom panel), when approximating the numerical solution by piecewise constant data. At time t n these discontinuities decay into three elementary waves, which propagate the solution forward to the next time level t n+1 (top panel). The timestep of the numerical scheme must satisfy the Courant-Friedrichs-Lewy condition, being small enough to prevent the waves from advancing more than ∆x/2 in ∆t.
In the frame of the local characteristic approach, the numerical fluxes appearing in Equation (76) are computed according to some generic flux formula that makes use of the characteristic information of the system. For example, in Roe's approximate Riemann solver [337] it adopts the following functional form: where w L and w R are the values of the primitive variables at the left and right sides, respectively, of a given cell interface. They are obtained from the cell-centered quantities after a suitable monotone reconstruction procedure. For a purely linear system, Equation (80) provides the exact expression for the flux in terms of the conserved variables. Therefore, the above expression for the Roe flux of conserved variables is the natural extension (after linearization) of the upwind flux for characteristic variables (see, e.g., [218]), once the quasilinear system is written in diagonal form, namely with where A is the Jacobian matrix of the quasilinear system, λ i are its eigenvalues, and R is the right-eigenvectors matrix. The upwind flux of the characteristic variables is given by , which yields the Roe flux given by Equation (80) when written in terms of the conserved variables.
The last term in the flux formula, Equation (80), represents the numerical viscosity of the scheme, and it makes explicit use of the characteristic information of the Jacobian matrices of the system. This information is used to provide the appropriate amount of numerical dissipation to obtain accurate representations of discontinuous solutions without excessive smearing, avoiding, at the same time, the growth of spurious numerical oscillations associated with the Gibbs phenomenon. In Equation (80), { λ n , r n } n=1...5 are, respectively, the eigenvalues and right eigenvectors of the Jacobian matrix of the flux vector. Correspondingly, the quantities {∆ ω n } n=1...5 are the jumps of the characteristic variables across each characteristic field. They are obtained by projecting the jumps of the state-vector variables with the left-eigenvector matrix: The "tilde" in Equations (80) and (82) indicates that the corresponding fields are averaged (local linearization) at the cell interfaces from the left and right (reconstructed) values. The way the cell-reconstructed variables are computed determines the spatial order of accuracy of the numerical algorithm and controls the amplitude of the local jumps at every cell interface. If these jumps are monotonically reduced, the scheme provides more accurate initial guesses for the solution of the local Riemann problems (the average values give only first-order accuracy). A wide variety of cell-reconstruction procedures is available in the literature. Among the slope-limiter procedures (see, e.g., [398,219]) most commonly used for TVD schemes [162] are the secondorder, piecewise-linear reconstruction, introduced by van Leer [403] in the design of the MUSCL scheme (Monotonic Upstream Scheme for Conservation Laws), and the third-order, piecewiseparabolic reconstruction developed by Colella and Woodward [80] in their Piecewise Parabolic Method (PPM). Since TVD schemes are only first-order accurate at local extrema, alternative reconstruction procedures for which some growth of the total variation is allowed have also been Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 35 developed. Among those, we mention the total variation bounded (TVB) schemes [377] and the essentially nonoscillatory (ENO) schemes [163] and extensions thereof.
Alternatively, high-order methods for nonlinear hyperbolic systems have also been designed using flux limiters rather than slope limiters, as in the Flux-Corrected Transport (FCT) scheme of Boris and Book [60]. In this approach, the numerical flux consists of two pieces, a high-order flux (e.g., the Lax-Wendroff flux) for smooth regions of the flow, and a low-order flux (e.g., the flux from some monotone method) near discontinuities, [398,219] for further details).
During the last few years most of the standard Riemann solvers developed in Newtonian fluid dynamics have been extended to relativistic hydrodynamics: Eulderink [117], as discussed in Section 2.2.1, explicitly derived a relativistic Roe's Riemann solver [337]; Schneider et al. [345] carried out the extension of Harten, Lax, van Leer, and Einfeldt's (HLLE) method [164,112] [102]. Recently, much effort has been spent concerning the exact special relativistic Riemann solver and its extension to multiple dimensions [238,322,333,334]. The interested reader is directed to [240] and references therein for a comprehensive description of such solvers in special-relativistic hydrodynamics.
Upwind HRSC schemes are general enough to be applicable to any hyperbolic system as long as the wave structure of the equations is known. As discussed before, the relativistic (and classical) MHD system of equations show degeneracies in the eigenvectors of the flux-vector Jacobian matrices. This fact makes it hazardous to solve them with linearized Riemann solvers based on the full spectral decomposition of the flux-vector Jacobians. For the case of special-relativistic MHD, approaches based on such full-wave decomposition have been put forward by [199,34,198]. Advances in this front have been significant, as even the exact solution of the Riemann problem in special-relativistic MHD has been obtained recently [340,150]. In addition, there has been a successful attempt by [24] to develop a GRMHD code in which a full-wave decomposition (Roe-type) Riemann solver based on a single, renormalized set of right and left eigenvectors has been implemented. This set of eigenvectors is regular for any physical state, including degeneracies (see [23] for details). Such a Riemann solver is invoked in the code of [24] after a (local) linear coordinate transformation based on the procedure developed by [320] that allows one to use special-relativistic Riemann solvers in general relativity, and which has been properly extended to include magnetic fields (see [24] for details on such extension and Section 6.1 for details of the procedure in the purely hydrodynamic case). A similar approach is followed in the GRMHD code of Komissarov [200].
Due to the numerical degeneracies of the GRMHD system most approaches to solve those equations use, however, incomplete Riemann solvers, i.e., simpler alternative approaches to compute the numerical fluxes such as the Harten-Lax-van Leer (HLL) single-state Riemann solver of [164] (see also [174] for an improved HLL scheme at contact discontinuities (HLLC) for ideal specialrelativistic MHD), or high-order central (symmetric) schemes, which entirely sidestep the use of the Jacobian matrix eigenvector structure. Such symmetric schemes, which have proven recently to yield results with an accuracy comparable to those provided by full-wave decomposition Riemann solvers in simulations involving both hydrodynamic and magnetohydrodynamic flows, are briefly discussed next.

High-order central schemes
The use of high-order nonoscillatory symmetric (central) TVD schemes for solving hyperbolic systems of conservation laws emerged in the mid 1980s [83,338,427,287] (see also [428] and [398] and references therein) as an alternative approach to upwind HRSC schemes. Symmetric Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 schemes are based on either high-order schemes (e.g., Lax-Wendroff) with additional dissipative terms [83,338,427], or on nonoscillatory high-order extensions of first-order central schemes (e.g., Lax-Friedrichs) [287,226]. One of the nicest properties of central schemes is that they exploit the conservation form of the Lax-Wendroff or Lax-Friedrichs schemes. Therefore, they yield the correct propagation speeds of all nonlinear waves appearing in the solution. Furthermore, central schemes sidestep the use of Riemann solvers, which results in enhanced computational efficiency, especially in multidimensional problems. Its use is, thus, not limited to those systems of equations in which the characteristic information is explicitly known or in which the Riemann problem is prohibitively expensive to compute.
For illustrative purposes let us write the numerical flux function resulting in the fully-discrete central scheme of Kurganov and Tadmor [211], to use in the conservation form algorithm given by Equation (76) This numerical flux depends only on the local propagation speeds a i+1/2 , given by the maximum of the eigenvalues of the Jacobian matrix at the cell interface i+1/2. Due to its simple form, it can be implemented and extended to any spatial order straightforwardly through the appropriate choice of monotone cell-reconstruction algorithms to compute the left and right states at cell interfaces. Conservative central schemes have been gradually developed during the last few years to reach a mature status where a number of characteristic-information-free central schemes of high order can be applied to any nonlinear hyperbolic system of conservation laws. The typical results obtained for the Euler equations show a quality comparable to that of upwind HRSC schemes, at the expense of a small loss of sharpness of the solution at discontinuities [398]. An up-to-date summary of the status and applications of this approach is discussed in [398,211,392].
In the context of special and general-relativistic MHD, Koide et al. [196,197] applied a secondorder central scheme with nonlinear dissipation developed by [83] to the study of black-hole accretion and formation of relativistic jets, investigating issues such as the magnetic extraction of rotational energy of the black hole [194]. More recently [89] assessed a state-of-the-art third-order, convex, essentially nonoscillatory, central scheme [226] in multidimensional special-relativistic hydrodynamics, later extended to relativistic MHD by [90]. These authors obtained results as accurate as those of upwind HRSC schemes in standard tests (shock tubes, shock reflection test). Yet another central scheme has been considered by [19,20] in one-dimensional special and generalrelativistic hydrodynamics and MHD, where results similar to those reported by [89,90] are discussed. The scheme of Kurganov-Tadmor (see Equation (83)) was assessed by [230] in the context of special-relativistic hydrodynamics, using standard numerical experiments such as shock tubes, the shock reflection test, and the relativistic version of the flat-faced step test. As for the other central schemes analyzed by [89,19] the results were comparable to those obtained by HRSC schemes based on Riemann solvers, even well inside the ultrarelativistic regime. Lucas-Serrano et al. [230] used high-order reconstruction procedures such as those provided by the PPM scheme [80] and the piecewise hyperbolic method (PHM) [102], which proved essential for keeping the inherent diffusion of central schemes at discontinuities at reasonably low levels. The scheme also produced accurate results in the case of full general-relativistic hydrodynamic simulations involving dynamic Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 spacetimes, as shown by [361] for oscillations of rapidly rotating neutron stars and the merger of neutron-star binaries. The scheme of Kurganov-Tadmor is also the adopted choice in the GRMHD simulations in dynamic spacetimes of [108,366]. Finally, a MUSCL-type scheme with HLL numerical fluxes is used in the code of [149], which allows a systematic investigation of GRMHD processes in accretion tori around black holes. Such HLL fluxes are also the choice used in the recent approaches of [396,91], which also implement a weighted, essentially nonoscillatory (WENO) method to build up fifth-order convergent numerical codes. Such high-order schemes may be suitable to solve the total (kinetic, thermal, and magnetic) energy equation in GRMHD codes, which deal with the challenging regime posed by high-Mach number flows.
It is worth emphasizing that early pioneer approaches in the field of relativistic hydrodynamics [291,73] used standard finite-difference schemes with artificial viscosity terms to stabilize the solution at discontinuities. However, as discussed in Section 2.1.2, those approaches only succeeded in obtaining accurate results for moderate values of the Lorentz factor (W ∼ 2). A key feature lacking in those investigations was the use of a conservative approach for both the system of equations and, accordingly, for the numerical schemes. In light of the results reported in recent investigations employing central schemes [89,19,230], it appears that this was the ultimate reason preventing the extension of early computations to the ultrarelativistic regime.
The alternative of using high-order component-wise central schemes instead of upwind HRSC schemes becomes apparent when the spectral decomposition of the hyperbolic system under study is unknown. The straightforwardness of a central scheme makes its use very appealing, especially in multiple dimensions where computational efficiency is an issue. Perhaps the most important example in relativistic astrophysics is the system of general-relativistic MHD equations for which, despite some progress, has been achieved in recent years (see, e.g., [34,199,24,23]), much more work is needed concerning their solution with full-wave-decomposition HRSC schemes based upon Riemann solvers. Meanwhile, an obvious choice is the use of central schemes.

Source terms
Most "conservation laws" involve source terms on the right-hand side (RHS) of the equations. In hydrodynamics, for instance, those terms arise when considering external forces such as gravity, which make the RHS of the momentum and energy equations no longer zero (see Section 2). Other effects leading to the appearance of source terms are radiative heat transfer (accounted for in the energy equation) and ionization (resulting in a collection of nonhomogeneous continuity equations for the mass of each species, which is not conserved separately). The incorporation of the source terms in the solution procedure is a common issue to all numerical schemes considered so far. Since a detailed discussion on the numerical treatment of source terms is beyond the scope of this article, we simply provide some basic information in this section, directing the interested reader to [398,219] and references therein for further details.
There are, essentially, two ways of handling source terms. The first approach is based on unsplit methods by which a single finite-difference formula advances the entire equation over one time step (as in Equation (76)): The temporal order of this basic algorithm can be improved by introducing successive substeps to perform the time update (e.g., predictor-corrector, Shu and Osher's conservative high-order Runge-Kutta schemes, etc.) Correspondingly, the second approach is based on fractional-step or splitting methods. The basic idea is to split the equation into different pieces (transport + sources) and apply appropriate methods for each piece independently. In the first-order Godunov splitting, Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 the operator L ∆t f solves the homogeneous PDE in the first step to yield the intermediate value U * . Then, in the second step, the operator L ∆t s solves the ordinary differential equation dU * /dt = S(U * ) to yield the final value U n+1 . In order to achieve second-order accuracy (assuming each independent method is second order) a common fractional-step method is the Strang splitting, This method advances by a half timestep the solution for the ODE containing the source terms, and by a full timestep the conservation law (the PDE) in between each source step.
We note that in some cases the source terms may become stiff, as in phenomena that either occur on much faster time scales than the hydrodynamic timestep ∆t, or act over much smaller spatial scales than the grid resolution ∆x. Stiff source terms may easily lead to numerical difficulties. The interested reader is directed to [219] (and references therein) for further information on various approaches to overcome the problems of stiff source terms.

Other techniques
Two of the most frequently employed alternatives to finite difference schemes in numerical hydrodynamics are smoothed particle hydrodynamics (SPH) and, to a lesser extent, (pseudo-) spectral methods. This section contains a brief description of both approaches. In addition, we also mention the field-dependent variation method and the finite element method. We note, however, that both of these approaches have barely been used so far in relativistic hydrodynamics.

Smoothed particle hydrodynamics
The Lagrangian particle method SPH, derived independently by Lucy [231] and Gingold and Monaghan [152], has shown successful performance in modelling fluid flows in astrophysics and cosmology. Most studies to date consider Newtonian flows and gravity, enhanced with the inclusion of the fluid self-gravity.
In the SPH method a finite set of extended Lagrangian particles replaces the continuum of hydrodynamic variables, the finite extent of the particles being determined by a smoothing function (the kernel) containing a characteristic length scale h. SPH-discretized representations of the equations of motion can be obtained after the introduction of mean-smoothed values for functions and for the gradient and divergence operators, as is briefly outlined below. Reviews of the classical SPH equations are abundant in the literature (see, e.g., [267,271] and references therein).
The main advantage of the SPH method is that it does not require a computational grid, avoiding mesh tangling and distortion. Hence, compared to grid-based finite-volume methods, SPH avoids wasting computational power in multidimensional applications, when, e.g., modelling regions containing large voids. Among the limitations of SPH we note the difficulties in modelling systems with extremely different characteristic lengths and the fact that boundary conditions usually require a more involved treatment than in finite-volume schemes. Comparisons between SPH and HRSC methods have been reported by various authors (see, e.g., [271,4] and references therein). In particular, the recent comparison study of state-of-the-art hydrodynamic codes carried out by Agertz et al. [4] has revealed fundamental differences between SPH and grid methods in modelling interacting multiphase fluids. More precisely, dynamic instabilities such as Kelvin-Helmholtz or Rayleigh-Taylor are much better resolved by Eulerian grid-based methods than by SPH techniques, the reason being the appearance of spurious pressure forces on particles in regions where shocks or steep gradients appear.
In the past few years, implementations of SPH to handle (special) relativistic (and even ultrarelativistic) flows have been developed (see, e.g., [79] and references therein). However, SPH has been so far scarcely applied in simulating relativistic flows in curved spacetimes. Relevant references include [189,212,213,381,300,122,299].
Given a function f (x), its mean smoothed value f (x) , (x = (x, y, z)) can be obtained from where W is the smoothing kernel, h the smoothing length and √ g d 3 x the volume element. The kernel must be differentiable at least once, and the derivative should be continuous to avoid large fluctuations in the force felt by a particle. Additional considerations for an appropriate election of the smoothing kernel can be found in [153]. The kernel is required to satisfy a normalization condition which is assured by choosing a normalization function and Ω(v), a standard spherical kernel. The smooth approximation of gradients of scalar functions can be written as and the approximation of the divergence of a vector reads Discrete representations of these procedures are obtained after introducing the number density distribution of particles n( ..,N the collection of Nparticles where the functions are known. The previous representations then read: ∇f with Ω ab ≡ Ω(x a , x b ; h). These approximations can then be used to derive the SPH version of the general-relativistic hydrodynamic equations. Explicit formulae are reported in [212]. The time evolution of the final system of ODEs is performed in [212] using a second-order Runge-Kutta time integrator with adaptive timestep. As in non-Riemann solver based finite-volume schemes, in SPH simulations involving the presence of shock waves, artificial viscosity terms must be introduced as a viscous pressure term [232].
Recently, Siegler and Riffert [381] have developed a Lagrangian conservation form of the generalrelativistic hydrodynamic equations for perfect fluids (with artificial viscosity) in arbitrary background spacetimes. Within that formulation these authors [381] have built a general-relativistic SPH code using the standard SPH formalism as known from Newtonian fluid dynamics (in contrast to previous approaches, e.g., [232,189,212]). The conservative character of their scheme has allowed the modelling of ultrarelativistic flows including shocks with Lorentz factors as large as 1000.

Spectral methods
The basic principle underlying spectral methods consists of transforming the partial differential equations into a system of ordinary differential equations by means of expanding the solution in a series on a complete basis. The mathematical theory of these schemes is presented in [156,71] (see also [159] for a recent Living Review article on spectral methods for numerical relativity). Spectral methods are particularly well suited to the solution of elliptic and parabolic equations. Good results can also be obtained for hyperbolic equations as long as no discontinuities appear in the solution. When a discontinuity is present, some amount of artificial viscosity must be added to avoid spurious oscillations. In the specific case of relativistic problems, where coupled systems of elliptic equations (i.e., the Einstein constraint equations) and hyperbolic equations (i.e., hydrodynamics) must be solved, an interesting strategy is to use spectral methods to solve the elliptic equations and HRSC schemes for the hyperbolic ones. Using such combined methods, first results were obtained in one-dimensional supernova collapse simulations, both in the framework of a tensor-scalar theory of gravitation [292,294] and in general relativity [293]. Multidimensional approaches for core collapse and studies of neutron star dynamics are available in [97,99].
Following [55] we illustrate the main ideas of spectral methods considering the quasi-linear one-dimensional scalar equation: with u = u(t, x) and λ a constant parameter. In the linear case (λ = 0), and assuming the function u to be periodic, spectral methods expand the function in a Fourier series: From the numerical point of view, the series is truncated for a suitable value of k. Hence, Equation (92), with λ = 0, can be rewritten as Finding a solution of the original equation is then equivalent to solving an "infinite" system of ordinary differential equations, where the initial values of the coefficients a k and b k are given by the Fourier expansion of u(x, 0).
In the nonlinear case, λ = 0, spectral methods proceed in a more convoluted way: first, the derivative of u is computed in the Fourier space. Then, a step back to the configuration space is taken through an inverse Fourier transform. Finally, after multiplying ∂u/∂x by u in the configuration space, the scheme returns again to the Fourier space.
The particular set of trigonometric functions used for the expansion in Equation (93) is chosen because it automatically fulfills the boundary conditions, and because a fast transform algorithm is available (the latter is no longer an issue for today's computers). If the initial or boundary conditions are not periodic, Fourier expansion is no longer useful because of the presence of a Gibbs phenomenon at the boundaries of the interval. Legendre or Chebyshev polynomials are, in this case, the most common base of functions used in the expansions (see [156,71] for a discussion on the different expansions).
Extensive numerical applications using (pseudo) spectral methods have been undertaken by the LUTH Relativity Group at the Observatoire de Paris in Meudon. The group has focused on Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 the study of compact objects, providing stationary initial data for single and binary compact stars and black holes, as well as performing dynamic studies of gravitational collapse to neutron stars and black holes. The Meudon group has developed a fully object-oriented library (based on the C++ computer language) called LORENE [227] where (multiple domain) spectral methods have been implemented within spherical coordinates. This library is currently being used by a growing number of numerical relativity groups worldwide, as can be inferred from the comprehensive summary of applications in general-relativistic astrophysics presented in [159] (see also [55] for an earlier review). Additional numerical relativity codes based on spectral methods have been developed by Ansorg and collaborators, who have built remarkably accurate stationary data for relativistic stars and black holes [21,22], by the Cornell/Caltech group, who are focused on the black-hole-binary problem, and by Faber et al. [122], who have investigated the black-hole-neutron-star binary system for a conformally-flat metric employing LORENE and an SPH code to solve for the hydrodynamics (see below).

Going further
The finite-element method is the preferred choice for solving multidimensional structural-engineering problems since the late 1960s [436]. First steps in the development of the finite-element method for modeling astrophysical flows in general relativity are discussed by Meier [253]. The method, designed in a fully covariant manner, is valid not only for the hydrodynamic equations but also for the Einstein and electromagnetic field equations. The most unique aspect of the approach presented in [253] is that the grid can be arbitrarily extended into the time dimension. Therefore, the standard finite-element method generalizes to a four-dimensional covariant technique on a spacetime grid, with the engineer's "isoparametric transformation" becoming the generalized Lorentz transform. At present, the method has shown its suitability for accurately computing the equilibrium stellar structure of Newtonian polytropes, either spherical or rotating. The main limitation of the finite-element method, as Meier points out [253], is that it is generally fully implicit, which results in severe computer time and memory limitations for three and four-dimensional problems.
Richardson and Chung [335] proposed recently the flow field-dependent variation (FDV) method as a new approach for general relativistic (non ideal) hydrodynamics computations. In the FDV method, parameters characteristic of the flow field are computed in order to guide the numerical scheme toward a solution. The basic idea is to expand the conservation flow variables into a Taylor series in terms of the FDV parameters, which are related to variations of physical parameters such as the Lorentz factor, the relativistic Reynolds number and the relativistic Froude number. The main drawback of the FDV method is the dependence of the solution procedure on a large number of problem-dependent parameters. Richardson and Chung [335] implemented the FDV method in a C++ code called GRAFSS (General-Relativistic Astrophysical Flow and Shock Solver). The test problems reported are the special relativistic shock tube (problem 1 in the notation of [239]) and the Bondi accretion on to a Schwarzschild black hole. While their method yields the correct wave propagation, the numerical solution near discontinuities has considerably more diffusion than with upwind HRSC schemes. Nevertheless, the generality of the approach is worth considering.

The magnetic-field divergence-free constraint
The numerical advantage of using the integral form of the conservation law, Equation (76), for the state vector variables does not apply in a straightforward way in the case of the GRMHD equations for the magnetic field components. While the induction equation is part of the hyperbolic system, there remains an additional constraint equation to be fulfilled, namely the magnetic field divergence-free constraint. Indeed, there is no guarantee that such divergence is conserved numerically when updating the magnetic field using the numerical procedure of Equation (76) for the rest of the components of the state vector (continuity, momentum, and energy equations).
Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Among the methods designed to preserve the divergence of the magnetic field (see [401] for a review), the constrained transport method (CT hereafter) designed by [119] and first extended to HRSC methods by [343] (see also [229] for a recent discussion), stands as the preferred choice for most groups. We note, however, that for AMR codes based on unstructured grids and multiple coordinate systems (such as those of [20,286]) there are other schemes, such as projection methods or hyperbolic divergence cleaning, which appear more suitable than the CT method to enforce the magnetic field divergence-free constraint. The projection method involves solving an elliptic equation for a corrected magnetic field projected onto the subspace of zero divergence solutions by a linear operator. More precisely the magnetic field is decomposed into a curl and a gradient as whose divergence leads to an elliptic (Poisson) equation ∆φ = ∇ B * , which can be solved for the scalar quantity φ. The magnetic field is next corrected according to B = B * − ∇φ. We note that an alternative projection scheme can be implemented by taking the curl of Equation (95) and solving for the vector potential A. Correspondingly, the basic idea of the hyperbolic divergence cleaning approach is to introduce an additional scalar field ψ, which is coupled to the magnetic field by a gradient term ∇ψ in the induction equation. The scalar field ψ is calculated by adding an additional constraint hyperbolic equation given by Therefore, divergence errors are propagated off the grid in a wave-like manner with characteristic speed c h . Further details on these two methods are given, e.g., in [401]. The approach followed by the CT scheme to ensure the solenoidal condition of the magnetic field is based on the use of Stokes theorem after the integration of the induction equation on surfaces of constant t and x i , Σ t,x i . Let us write the induction equation as where we have defined the densitized magnetic field vector B = √ γ B and Ω = (α v − β) × B.
Following [24] we can obtain a discretized version of Equation (97) as follows. At a given time, each numerical cell is bounded by six two-surfaces. Let us consider the two-surface Σ t,x 3 , defined by t = const. and x 3 = const., and the remaining two coordinates spanning the intervals from x 1 to x 1 + ∆x 1 , and from x 2 to x 2 + ∆x 2 . The magnetic flux through this two-surface is given by Integrating Equation (97) on the two-surface Σ t,x 3 , and applying Stokes theorem to the right hand side leads to the following equation with ∂(Σ t,x 3 ) being the boundary of the corresponding two-surface. This equation can be integrated to give Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 where the caret denotes that quantitiesΩ i are calculated at the edges of the numerical cells, where they can be discontinuous. At each edge, these quantities can be calculated using the solution of four Riemann problems between the corresponding faces, whose intersection defines the edge. It is important to emphasize, however, that irrespective of the expression used to computeΩ i , the method to advance the magnetic fluxes at the faces of the numerical cells satisfies, by construction, the divergence constraint (see, e.g., [24] for details). Indeed, the total magnetic flux through the closed surface S formed by the boundary of the numerical cells and enclosing a volume V can be calculated as after applying the Gauss theorem and divergence constraint. This equation implies that no magnetic flux source exists inside the volume V and, therefore, the magnetic flux is a conserved quantity as dΦ T /dt = 0. As a result, every numerical scheme constructed on the basis of the CT scheme will conserve magnetic flux up to machine accuracy. If one is able to generate initial conditions, which satisfy the divergence constraint, i.e., with Φ = 0 across the boundary of each numerical cell, then such a constraint will be preserved during the numerical evolution. Any of the Riemann solvers and flux formulae discussed in Section 4.1.2 can be thus used to calculate the quantitiesΩ i needed to advance in time the magnetic fluxes following Equation (100). At each edge of the numerical cell,Ω i is written as an average of the numerical fluxes calculated at the interfaces between the faces whose intersection define the edge. Let us consider, for illustrative purposes,Ω x . If the indices (j, k, l) denote the center of a numerical cell, an x−edge is defined by the indices (j, k + 1/2, l + 1/2). By definition, Ω x = α(ṽ y B z −ṽ z B y ). Since F y (B z ) =ṽ y B z −ṽ z B y and F z (B y ) =ṽ z B y −ṽ y B z , we can expressΩ x in terms of the fluxes as followŝ whereF y (F z ) refers to the numerical flux in the y (z) direction corresponding to the equation for B z (B y ) and multiplied by α to account for the correct definition of Ω. Analogous expressions for Ω y j+1/2,k,l+1/2 andΩ z j+1/2,k+1/2,l , can be derived in a similar way. We refer the interested reader to [401] for additional details on the constrained transport method.

State-of-the-art three-dimensional codes
We instead direct the interested reader to [323] and references therein for an up-to-date overview of such approaches. We also consider only multidimensional (3D) codes.

Hydrodynamics
We start with those purely hydrodynamic codes which also consider dynamic spacetimes. Hydrodynamic codes, which only solve for the matter dynamics on fixed background metrics (quite more numerous), are listed in Table 1, where we give their main basic features and provide the relevant references for further reading. Correspondingly, Table 2 lists test-fluid GRMHD codes, which could also, in principle, deal with unmagnetized flows by simply setting to zero the magnetic field contribution. In addition, the reader is directed to the Living Review article of [240] where similar tables can be found for the particular case of special relativistic hydrodynamics. CACTUS/WHISKY: The whisky code was the result of a collaboration among several European institutions [246,27]. This code solves the general relativistic hydrodynamics equations formulated as in Section 2.1.3 on a 3D numerical grid with Cartesian coordinates. The code has been constructed within the framework of the Cactus Computational Toolkit (see [244] for details), developed at the Albert Einstein Institute (Potsdam-Golm, Germany) and at Louisiana State University (Baton Rouge, Louisiana). This public-domain code provides high-level facilities such as parallelization, input/output, portability on different platforms and several evolution schemes to solve general systems of partial differential equations. Special attention is dedicated to the solution of the Einstein equations, whose matter terms in nonvacuum spacetimes are handled by the whisky code. The initial development of whisky benefitted in part from the release of a public version of the general-relativistic hydrodynamics code gr astro described in [137,131], and developed mostly by the group at Washington University (St. Louis, Missouri). The main features of the whisky code are (a) cell-reconstruction procedures with minmod and monotonized central (MC) slope limiters, along with PPM and ENO methods; (b) full-wave decomposition and incomplete Riemann solvers to compute the numerical fluxes such as HLLE, Roe-like, and Marquina flux formula; (c) analytic expression for the right and left-eigenvectors of flux-vector Jacobian matrices [176] and compact flux formulae [10] for the Roe-Riemann solver and Marquina's flux formula; (d) a "method of lines" (MoL) approach for the implementation of conservative, high-order time evolution schemes; (e) the possibility to couple the GRHD equations with a conformally-decomposed three-metric in the BSSN formulation [363,39]. It is worth noting that whisky also incorporates excision methods adapted to HRSC schemes for studying black-hole formation, as described in [166]. In particular, a modified PPM reconstruction scheme was developed to handle the boundary of the excised region inside apparent horizons. The accuracy of the simulations performed with the whisky code can also benefit from the mesh refinement driver implemented in the cactus code, called carpet, as recently demonstrated in the gravitational collapse studies of [28,305].
CoCoNuT: This code, described in detail in [97], was designed with the aim of studying generalrelativistic astrophysical scenarios such as rotational core collapse to neutron stars and black holes, as well as pulsations and instabilities of the formed compact objects. Contrary to most 3D codes, which are based on Cartesian coordinates, the CoCoNuT code employs spherical coordinates. Upwind HRSC methods are used for the hydrodynamic equations in the formulation of Section 2.1.3 (Marquina's flux formula and PPM cell-reconstruction procedures) and spectral methods (through the LORENE library) for solving the metric equations which are formulated in the conformal flatness (CF) approximation (an approach, which [97] coined Mariage des Maillages, i.e., French for grid wedding). In this approximation, and under the maximal slicing condition, the 3+1 ADM equations reduce to a set of five coupled elliptic (Poisson-like) nonlinear equations for the metric components (lapse function, shift vector and conformal factor), which are optimally suited to be solved using spectral methods (see Section 4.2.2). The CF approximation has proved to agree remarkably well with BSSN in simulations performed with the whisky code for rotating neutron stars and stellar Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 core collapse [98,305]. Extensions of the CF approximation also implemented in the CoCoNuT code are considered by [74].
Shibata's GRHD code: In this code the hydrodynamics equations are formulated both in a nonconservative way (following Wilson's approach of Section 2.1.2; see [353]) and in a conservative way similar to the approach of Section 2.1.3 (see [355]). In the former case, an important difference with respect to the original system given by Equations (20) - (22), is that an equation for the entropy is solved instead of the energy equation. The hydrodynamic equations are integrated using van Leer's [403] second-order finite difference scheme with artificial viscosity, following the approach of a precursor code developed by [303]. For the conservative system, both upwind HRSC schemes and high-order central schemes are implemented in the code, along with third-order parabolic cell-reconstruction procedures. As mentioned before, the choice of conserved variables in the conservative formulation is, however, slightly different to that of Section 2.1.3, mainly due to the presence of a common factor e 6φ in all quantities, φ being a spacetime conformal factor, motivated by the formulation used in the code for solving Einstein's equations. Those are written using the BSSN formulation and solved with finite differences (see [355] for details particular to the formulation of the spacetime variables and their solution). In addition, [355] uses the transport velocity throughout (instead ofṽ i ) and in the energy equation the continuity equation is not subtracted (contrary to the approach of [36]). The choice of coordinates in Shibata's code depends on the dimensionality of the simulation under study. In the full 3D case, Cartesian coordinates (x, y, z) are employed. While this is also the case in axisymmetry there is an important subtle difference worth mentioning; the hydrodynamics equations are first written using cylindrical coordinates ( , φ, z), with = x 2 + y 2 and φ = arctan(y/x). However, the Einstein's equations are solved in the y = 0 plane using Cartesian coordinates. To be able to compute y derivatives in this plane the code implements the "cartoon" method [6], which makes possible axisymmetric computations with a Cartesian grid. Next, the hydrodynamics equations are rewritten in Cartesian coordinates using appropriate relations at the y = 0 plane (see [355] for details). Comprehensive tests of the code are described in [353,355,361]. Shibata's code has allowed important breakthroughs in the study of the morphology and dynamics of various general-relativistic astrophysical problems, such as black-hole formation resulting from both the gravitational collapse of rotating neutron stars and rotating supermassive stars, dynamic instabilities of rotating neutron stars, and, perhaps most importantly, the coalescence of neutron-star binaries, a long-standing problem in numerical relativistic hydrodynamics. These applications are discussed in Section 5.
Duez et al. [109]: The code of [109] shares many features with that of Shibata just discussed. As in the case of [355] the Einstein equations are formulated in the BSSN approach and solved using an iterative Crank-Nicholson scheme for the time update and second-order centered differencing for the spatial derivatives. Correspondingly, the hydrodynamics equations are formulated in the same conservative approach followed by [355], yet they are not solved using HRSC schemes, but through an artificial viscosity (either quadratic or linear), which allows it to handle shock waves. A noteworthy feature of this code is the incorporation of an algorithm, which makes unnecessary the addition of a low-density atmosphere required in most Eulerian hydrodynamic codes, either Newtonian or relativist, when studying the dynamics of matter sources of compact support (as, e.g., the pulsations of stars). Technical details of boundary conditions and gauge choices, along with a comprehensive list of tests passed by the code, are available in [109]. Extensions of this code to account for viscous fluids through the solution of the relativistic Navier-Stokes equations are reported in [107].
Oechslin et al. [299]: As in the CoCoNuT code described above, in the code of [299] the relativistic hydrodynamics equations are solved together with the Einstein equations in the conformally-flat approximation. The code has evolved from an early Newtonian version, which was designed with the definite aim of studying the merger of neutron star binaries. It has gradually been improved to provide ever more accurate descriptions of such mergers. In its latest version the code incorporates Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 47 microphysical treatment of neutron-star matter through two finite-temperature nuclear EOSs along with a modern treatment to extract gravitational waveform information. Compared to the other codes mentioned before, the most noticeable feature of the code of [299] is the use of SPH techniques to solve for the relativistic hydrodynamics equations, which are formulated in conservation form.
The code implements the new SPH artificial viscosity formalism of [79] in which the artificial viscosity interaction is determined by approximately solving a Riemann problem between particle neighbors. Binary merger simulations performed with this code are discussed in Section 5. Characteristic numerical relativity codes: Although most numerical relativity codes are based on the Cauchy problem, there also exist a couple of characteristic numerical relativity codes, which can handle hydrodynamics. On the one hand, there is the axisymmetric code of [378,379], which has been used to perform dynamic evolutions of neutron stars, both pulsating and those formed after a core collapse. In this code the hydrodynamics equations are implemented following the conservative approach outlined in Section 2.2.2 and solved using upwind HRSC schemes, high-order cell-reconstruction procedures, and conservative Runge-Kutta schemes for the time update. Regarding full 3D, the Pittsburgh characteristic vacuum code has recently been upgraded to account for perfect fluid matter sources [46], using the same conservative formulation of Section 2.2.2. However, instead of relying on Riemann-solver-based methods for their solution, the hydrodynamics equations are solved in the Pitt code with a dissipative central scheme designed by [83]. Using this code, short time evolutions of a self-gravitating star in close orbit around a black hole are discussed in [46].
As an important note we point out that all codes based on the Cauchy approach place the outer boundary of the grid at a finite distance, which may potentially lead to numerical problems caused by unwanted reflections of outgoing waves back into the computational domain. Further work on the development of sophisticated boundary conditions is needed to solve (or at least ameliorate) this type of problem. Alternative solutions, which allow for compactified spacetimes that include future null infinity on the computational grid, are provided by the light-cone approach developed by Winicour et al. [419] or the conformal formalism of Friedrich [142].

Magnetohydrodynamics
As in the previous Section 4.4.1, we focus on magnetohydrodynamic codes, which also consider dynamic spacetimes. GRMHD codes which "only" solve for the matter dynamics on fixed background metrics are listed in Table 2, where we give their main basic features and provide references for further reading. As mentioned before, the development of such test-fluid GRMHD codes has witnessed spectacular progress in the last few years. This has been possible not only thanks to the various conservative formulations, which have been put forward but also, and perhaps most importantly, because of the use of accurate, characteristic-information-free HRSC schemes. Such schemes have rendered possible systematic investigations of strongly magnetized scenarios of relativistic astrophysics, which had thus far remained out of numerical reach. We note that in the characterization followed in Table 2 to classify the various existing codes we have made the explicit distinction between full-wave-decomposition approximate Riemann solvers (such as those used in the codes of [201,24]) and incomplete Riemann solvers in which only a subset of wave speeds are used. In this regard, we place the widely used HLL (and HLLE) scheme in the latter group, together with symmetric schemes, despite the fact that HLL was originally developed as an approximate Riemann solver scheme [164] with a particularly simple numerical flux.
Shibata's GRMHD code [366]: As described in Section 3.1.2 the GRMHD equations are formulated in the code of [366] in a conservative way. For their solution the semi-discrete highresolution central scheme of Kurganov-Tadmor [211] is used (see also [361]), together with a piecewise parabolic interpolation of the primitive variables for the cell-reconstruction scheme. In order to obtain the maximum characteristic speeds needed for the central scheme, instead of solving the Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 quartic equation in the characteristic polynomial, Shibata and Sekiguchi [366] follow the convenient prescription proposed by [149] and approximate the fourth-order equation by a quadratic equation for which there exists an analytic solution. The magnetic field divergence-free constraint is enforced with a constrained transport scheme. Correspondingly, and as in the purely hydrodynamic code, Einstein's equations are also written in the GRMHD code using the BSSN formulation and solved with finite differences (see [355] for further details). The code is axisymmetric and, as in the hydrodynamic case, the "cartoon" method is employed [6] for evolving the BSSN equations, which makes possible the use of Cartesian coordinates in axisymmetric simulations, and a cylindrical grid is used for the GRMHD equations. It has recently been applied in the study of the collapse of magnetized hypermassive neutron stars to black holes in the context of gamma-ray burst mechanisms, as we discuss in Section 5.
Duez et al. [108]: The development of this code has proceeded simultaneously to that of [366], to the point that they share some important features and both have been used to study similar astrophysical problems, whose results have been presented in joint publications [105,106]. This code solves the GRMHD equations in both two dimensions (axisymmetry) and three dimensions, using a Cartesian grid for the latter case and the "cartoon" method in the former (although a cylindrical grid is used for the induction and MHD evolution equations). As in the purely hydrodynamic code of the same group [109] the Einstein equations are formulated in the BSSN approach and solved using a three-step iterative Crank-Nicholson scheme for the time update, which yields second-order accuracy in time. The main modification of the GRMHD code with respect to its predecessor is the complete reformulation of the hydrodynamic sector, whose accuracy has been improved through the use of shock-capturing capabilities. More precisely, the GRMHD evolution equations are solved with the HLL approximate Riemann solver, together with different possibilities for the cell-reconstruction step (either MC, CENO, or PPM). It is worth noting that, contrary to the no-atmosphere approach used in the purely hydrodynamic code [109] to handle vacuum regions outside stars, the GRMHD code of [108] is not suitable for such an approach and a very small positive density must be maintained outside the stars.
WHISKY-MHD: As the name indicates, this code is the extension of the hydrodynamic code whisky discussed above to solve for the full set of GRMHD equations in a dynamic spacetime. The code is described in [151]. It is a fully three-dimensional code in Cartesian coordinates, based on HRSC techniques on domains with adaptive mesh refinement, accomplished as in the whisky code through the AMR driver implemented in the cactus code, called carpet. The whisky-MHD code uses, hence, the cactus computational toolkit which provides the infrastructure for the parallelization and the I/O of the code, along with methods for the solution of the Einstein equations formulated in the BSSN approach [363,39]. The GRMHD equations are implemented using the conservative formulation of [24] discussed in Section 3.1.2 and integrated in time employing the method of lines. The code uses an HLLE approximate Riemann solver, with second-order TVD slope limiters for the cell-reconstruction procedure, and the magnetic-field divergence-free constraint is enforced using the flux-CT scheme briefly discussed in Section 4.3.
Anderson et al. [12]: The most distinctive feature of this code, which is mainly described in [12] (see also [286,14] for further details and [13] for the assessment of the purely hydrodynamic module of the code), is the fact that it uses a sophisticated computational infrastructure, which provides distributed AMR. This has allowed the investigation of both unmagnetized and magnetized neutron-star binary mergers and the computation of the gravitational radiation in the wave zone far from the merger [13,12]. As mentioned in Section 3.1.2 the code uses a conservative formulation of the GRMHD equations. These are solved by means of the HLLE approximate Riemann solver with PPM cell reconstruction, in conjunction with divergence cleaning to control the divergence-free constraint. Correspondingly, the Einstein equations, cast in first-order symmetric hyperbolic form, are solved in the generalized harmonic decomposition.
Cerdá-Durán et al. [75]: The GRMHD code developed by [75] has been designed to study Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 49 Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 magneto-rotational, relativistic, stellar core collapse [76]. It is an extension of the axisymmetric hydrodynamics code developed by [95] (whose 3D extension constitutes the CoCoNuT code described above), in which magnetic fields are included following the approach laid out in [24]. Einstein's equations are formulated using the conformally-flat condition, which has proved very accurate for studying rotational core collapse [96,98] and are solved using spectral methods. In the code, both the HLLE approximate Riemann solver and the Kurganov-Tadmor central scheme are implemented to solve for the MHD evolution equations, along with the flux-CT scheme for the magnetic-field divergence-free constraint. As a first step towards magneto-rotational core collapse simulations, an early version of the code assumed a passive (test) magnetic field, a justifiable assumption since weakly-magnetized fluids are present in such astrophysical scenarios. An extension of the code to relax this assumption has recently been finished.

Hydrodynamic and Magnetohydrodynamic Simulations in Relativistic Astrophysics
With the exception of the vacuum two-body problem (i.e., the coalescence of two black holes), all astrophysical systems and sources of gravitational radiation involve matter. Therefore, although the black-hole-binary problem has always been considered the paradigmatic, and more demanding, problem to solve numerically, the joint integration of the equations of motion for matter and geometry was in the minds of theorists from the very beginning of numerical relativity. Nowadays there is a large (and increasing) body of numerical investigations in the literature dealing with hydrodynamic integrations in static background spacetimes. While such a statement also holds true in the case of the MHD equations, particularly from the last few years onwards, the development still awaits thorough numerical exploration. In the purely hydrodynamic case most of the investigations are based on Wilson's Eulerian formulation of the hydrodynamic equations and use schemes based on finite differences with some amount of artificial viscosity. The use of conservative formulations of the equations, and the incorporation of the characteristic information in the design of numerical schemes, only began in more recent years. However, conservative approaches are turning into the most common choice for the recent code developments, not only for hydrodynamics but particularly for GRMHD (as can be inferred from Table 2).
On the other hand, time-dependent simulations of self-gravitating flows in general relativity (evolving the spacetime dynamically with the Einstein equations coupled to a hydrodynamic source) still constitute a much smaller sample. Although there is much interest in this direction, only the spherically symmetric case (1D) has been extensively studied. In axisymmetry (2D) fewer attempts have been made, with most of them devoted to the study of the gravitational collapse of rotating stellar cores, black-hole formation, and the subsequent emission of gravitational radiation. Threedimensional simulations have only started more recently but are already becoming mainstream as improved theoretical and numerical approaches, along with increasing computational power, become available. Much effort is focused on the study of neutron star instabilities, gravitational collapse, and on the coalescence of compact neutron star binaries, all scenarios being investigated with the incorporation of microphysical models for neutron star matter and magnetic fields. These three-dimensional efforts also have a vacuum counterpart, the black-hole-binary problem, which has witnessed spectacular progress in the last few years [323]. A prime driver for these investigations is the emerging possibility of soon detecting gravitational waves with the ongoing experimental endeavors. The waveform catalogues resulting from time-dependent hydrodynamic and MHD simulations provide essential help to data analysis groups, since the chances for detection are noticeably enhanced for most sources through matched-filtering techniques.
In the following, we review the status of the numerical investigations in three distinctive astrophysical scenarios, which involve strong gravitational fields and, hence, relativistic physics: gravitational collapse, accretion onto black holes, and magnetohydrodynamic evolutions of neutron stars. Relativistic cosmology, another area where fundamental advances have been accomplished through numerical simulations, is not considered; the interested reader is directed to the Living Reviews article by Anninos [18] and references therein.

Gravitational collapse
The study of the gravitational collapse of massive stars has largely been pursued numerically over the years. This problem was already the main motivation when May and White built the first one-dimensional code [247,248] (see Section 2.1.1). Such codes are designed to integrate the coupled system of Einstein field equations and relativistic hydrodynamics, sidestepping Newtonian approaches.
Since early investigations the numerical study of gravitational collapse has been neatly divided Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 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 [146] 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 [78]. 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 [161] and references therein for details.

Core collapse supernovae
At the end of their thermonuclear evolution, massive stars in the (Main Sequence) mass range from 9 M to 30 M 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 ∼ 3 × 10 53 erg (M/M ) 2 (R/10 km) −1 , 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]  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.

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. The interested reader is directed to previous versions of this review for a chronological presentation of the most significant advances achieved with those early codes. The continuous numerical effort of the last 40 years has clearly shown the strong dependence of the explosion mechanism on the details of the post-bounce evolution: general relativity, the nuclear EOS and the properties of the nascent neutron star, the treatment of neutrino transport, and the neutrino-matter interaction. Recently, state-of-the-art treatments of neutrino transport have been achieved, even in the general-relativistic case, by solving the Boltzmann equation in connection with hydrodynamics, instead of using multigroup flux-limited diffusion [254,326,327,222,223,68,233]. The most recent results from the few spherically-symmetric simulations available show that neutrino-driven explosions may be produced in spherical symmetry for low-mass progenitors (below ∼ 10 M ) due to the particular density distribution in some of the burning layers of the star (the C-layer and the transition from this layer to the He-layer), which facilitates the re-expansion of the stalled prompt shock [192]. That is, however, not the case Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 for more massive progenitors, including the typical 15 M model. In this case, computationally expensive multidimensional 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 (see below).
From the numerical point of view, many of the earlier investigations of supernova core collapse 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, started in the late 1980s with the Newtonian simulations of Fryxell et al. [148], who developed an Eulerian second-order PPM code called PROMETHEUS (see [271] for a review of the present status). A comparison of the energetics of a stellar collapse simulation using this HRSC scheme and a standard May and White code with artificial viscosity was performed by Martí et al. [236] for the same initial model, grid size and EOS. These authors used Glaister's approximate Riemann solver [154] in a Lagrangian Newtonian hydrodynamic code. The outcome of the comparison was 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. [186] repeated this computation with a different EOS, using a PPM second-order Godunov-type scheme, disagreeing with Martí et al. [236]. The state-of-the-art implementation of the tensorial artificial viscosity terms in [186], together with the very fine numerical grids employed (unrealistic for three-dimensional simulations), could be the reason for the discrepancies.
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 [223].
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 [339] (see [424] 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 Γ = Γ min +η(log ρ−log ρ b ), where ρ b is the bounce density and η = 0 if ρ < ρ b and η > 0 otherwise [406]. A set of animations of the simulations presented in [339] is included in the four movies in compression is ∼ 1.0 M . The minimum value that the quantity α 2 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 α 2 ∼ 0.75).  [190]. 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. Visualization by Konstantinos Kifonidis. (To watch the movie, please go to the online version of this review article at http://www.livingreviews.org/lrr-2008-7.)

Axisymmetric Newtonian simulations.
Beyond spherical symmetry most of the existing simulations of core collapse supernova are still Newtonian. Axisymmetry allows for investigations of rotational core collapse, which are highly motivated nowadays by the prospects of direct detection of the gravitational waves emitted. The role played by rotation in core collapse progenitors is still primarily an open issue, being currently incorporated only in stellar evolution codes [172]. Observations of surface velocities support the fact that a large fraction of progenitor cores is rapidly rotating. However, magnetic torques can spin down the core by dynamo action, which couples to the outer layers of the star.
Axisymmetric core collapse simulations have been performed by [270,52,56] using a realistic EOS and some treatment of weak interaction processes, but neglecting neutrino transport, and by [185,268,181,144,145,208,67,273,408,68] employing some description of neutrino transport. The main challenge encountered in the incorporation of neutrino physics in multidimensional collapse simulations is that the resulting schemes are either very computationally expensive [67,273,408,68] or need to rely on simplifications of neutrino transport and its microphysics [208].
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 [243]). 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 15 M pre-collapse model of a blue super-giant star in [423]. The computations begin 20 ms after core bounce from a one-dimensional model of [65]. This model is obtained with the one-dimensional general-relativistic code mentioned previously [64], 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 Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 57 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 [190] 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 [184] and references therein). The SASI shows a preferential growth of low ( = 1 and = 2) shock-deformation modes at the time of shock revival, leading to large-amplitude bipolar oscillations, which have important implications for largescale explosion asymmetries, pulsar kicks, and the development of neutrino driven explosions.
A new line of investigation led by Burrows et al. [70] 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 g-mode oscillations. The simulations of [70] for an 11 M progenitor show that an = 1 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 [184].
In addition to understanding the explosion mechanism, one of the prime motivations for performing multidimensional core-collapse simulations is to compute theoretical gravitationalradiation waveforms, which may be confronted by experimental data currently being collected by a number of detectors. As early as 1982, Müller [270] obtained the first numerical evidence of the low gravitational wave efficiency of the core collapse scenario, with an energy smaller than 10 −6 M c 2 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 postbounce wave signal shows large amplitude oscillations associated with pulsations in the collapsed core [439,328] and (possibly) low−T /|W | rotational dynamic instabilities [306,305] (T and W 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 [273]. 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 (∼ 10 12 −10 13 G) or even larger at birth (∼ 10 14 −10 15 G). In extreme cases such as magnetars, the magnetic field can be so strong as to affect the internal structure of the star [51]. The presence of such intense magnetic fields Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 render magneto-rotational core collapse simulations mandatory. We note that as early as 1979 magneto-rotational core collapse was already proposed by [272] as a plausible supernova explosion mechanism.
In recent years, an increasing number of authors have performed axisymmetric magnetorotational 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 [209] and [426], who found differences with purely hydrodynamic models only for very strong initial fields (≥ 10 12 G). 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 (≤ 10 11 G, 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 (≥ 10 11 G) 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, [296] 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 [69]. 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.

Axisymmetric relativistic simulations.
Early investigations in general-relativistic rotational core collapse were mainly concerned with the question of black-hole formation under idealized conditions (see Section 5.1.2), but none of these studies have 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. In addition, as general relativity counteracts the stabilizing effect of rotation, a bounce caused by rotation has to occur at larger densities than in the Newtonian case, which further calls for studies based on a relativistic framework. Only since the early 2000s [93,94,95,96,97,379,365,367,74,98,305] have relativistic simulations of rotational core collapse become possible, although they still mostly employ idealized collapse models based on rotating polytropes.
Wilson [413] first computed neutron star bounces of Γ = 2 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 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 [418] (CFC hereafter). Correspondingly, the hydrodynamic equations were cast as the Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 59 first-order flux-conservative hyperbolic system described in Section 2.1.3. Details of the methodology and numerical code are given in [95] (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 [186]. 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 [439] -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 [96] 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 [439]. On average, thus, the CFC simulations of [96] show that the gravitational wave signals of relativistic models have similar amplitude but somewhat higher frequency than their Newtonian counterparts [439].
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 singlebounce 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 [96] 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. [96] is available at the MPA's website [243].
The movie presented in Figure 6 shows the time evolution of a multiple bounce model (model A2B4G1 in the notation of [95]), 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 [96] can be found at the MPA's website [243].
Relativistic simulations with improved dynamics and gravitational waveforms were reported by [74], 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 [96] a simple (modified polytrope) EOS was used, and the initial models were unmagnetized. The CFC simulations of [74] cover the basic morphology and dynamics of core collapse types studied by [96], 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 [74], 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 Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 the collapse to a black hole.) Comparisons of the CFC approach with fully general-relativistic simulations (employing the BSSN formulation) have been reported by [365] 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 [92], 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 [97]) 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 20 M 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 [351] together with an approximate (parameterized) neutrino treatment with deleptonization and neutrino pressure effects included, as recently suggested by [221]. In the simulations of [305,98] the electron fraction profiles are Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 61 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 [221] (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 [362], 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 (∼ 10 12 -10 13 G). Although these values are probably astrophysically not relevant, because the stellar evolution models of [172] predict a poloidal field strength of ∼ 10 6 G 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 λ crit ∼ 2πc A /Ω, where c A is the Alfvén speed and Ω is the angular velocity. The fastest growing MRI mode develops at length scales near λ crit on a typical time scale of τ MRI = 4π[ ∂ Ω] −1 , where = r sin θ 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 λ crit , which can be prohibitively small to capture for initially small magnetic fields (say, on the order of meters for a 10 10 G progenitor), especially for three-dimensional computations, where the affordable grid resolution is still limited. The results of [362] 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 (≤ 10 10 G), 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 T µν satisfy b µ b ν ρhu µ u ν + P g µν . 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 [76] presents the first relativistic simulations of magneto-rotational core collapse, which take into account the effects of a Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 microphysical EOS [351] and a simplified neutrino treatment. These effects have been incorporated into the code following the procedure employed by [305] and [98] 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. [76] 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 [76] 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 [362]). 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 simulations.
Most of the existing three-dimensional core collapse computations available in the literature have employed Newtonian physics (see, e.g., [56,328,63,147]), and only very recently have the first relativistic studies been performed [367,305]. In either case, no three-dimensional simulation, which includes neutrino transport, has yet been reported.
Concerning Newtonian studies, Bonazzola and Marck [56] performed pioneering 3D 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 [270]). 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, [328,63] have extended the work of [439] studying, with two different PPM hydrodynamics codes, the dynamic growth of nonaxisymmetric (bar mode) instabilities appearing in rotational post-collapse cores. A relativistic extension has been performed by [358,367] (see Section 5.3.1). The work of [147], based on an SPH supernova code with simple flux-limited diffusion to transport the neutrinos, shows that for rapidly rotating stars (unlikely from stellar evolution calculations), rotation modifies the convection above the proto-neutron star, but it is not fast enough to cause core fragmentation. These simulations also show that the protoneutron star does not spin fast enough to generate strong magnetic fields quickly after collapse, which hints that magnetic fields may not dominate the supernova explosion mechanism.
Three-dimensional and fully general-relativistic simulations of core collapse have been performed by [367] and [305]. 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 (0.01 ≤ T /|W | ≤ 0.02) and has a high degree of differential rotation, which seems unlikely according to stellar evolution models.
On the other hand, [305] have performed 3D simulations of general-relativistic core collapse within the BSSN framework and a finite-temperature nuclear EOS [351] and an approximate  (parameterized) neutrino treatment [221]. 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 [98]. 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 [98]. One of the models, labelled E20A in [305], 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 T /|W | ratio in the newly-formed proto-neutron star high enough to eventually develop low T /|W | nonaxisymmetric (bar-mode-like) deformations. This had already been found in a similar model using Newtonian physics in [306].
The left panel of Figure 7 displays the time evolution of the normalized mode amplitudes of the four lowest azimuthal density modes (∝ e imφ ) in the equatorial plane, and shows that the one-armed m = 1 mode (red curve) becomes dominant from t ∼ 20 ms after core bounce. Correspondingly, the right panel of Figure 7 plots the gravitational wave strains h + and h × 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, m = 4 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 m = 1, m = 2, and m = 3 features. The mode structure is nonstationary and exhibits strong dynamic variation with radius and time.
Ott et al. [305] find that the characteristic wave strain of this model peaks at a frequency of Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 ∼ 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 T /|W | 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.

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. 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.

Axisymmetric simulations.
Beyond spherical symmetry the investigations of blackhole formation started with the work of Nakamura [277], 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) [279]. The equations were finite-differenced using the donor cell scheme plus Friedrichs-Lax type viscosity terms. Nakamura used a "hypergeometric" slicing con-Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 dition (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 10 M "core" of a massive star with different amounts of rotational energy and an initial central density of 3 × 10 13 g cm −3 , up to the formation of a rotating black hole. However, due to limitations in computer resources, the affordable resolution (42 × 42 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 [278] (see also [281]) considered a configuration consisting of a neutron star (M = 1.09 M , ρ c = 10 15 g cm −3 ) with an accreted envelope of 0.81 M , which was thought to mimic mass fallback in a supernova explosion. Rotation and infall velocity were added to such configurations, simulating evolution depending on the prescribed rotation rates and rotation laws.
Bardeen, Stark, and Piran, in a series of papers [38,386,319,387], studied the collapse of rotating relativistic (Γ = 2) 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 1.9 × 10 15 (M/M ) −2 , and radius 6 GM/c 2 = 8.8 × 10 5 M/M cm. 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 ∆E/M c 2 < 7 × 10 −4 , and that the general waveform shape was nearly independent of the details of the collapse (but see [28,30] for new recent estimates of this energy figure).
Evans [118], building on previous work by Dykema [111], 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 (r = 0) and at the polar axis (θ = 0, π), where some fields diverge due to the coordinate singularities. Evans made important contributions toward regularizing the gravitational field equations in such situations [118]. 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 [385] and by Dimmelmeier et al. [97]. 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. [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 [354] investigated the effects of rotation on the criterion for prompt adiabatic collapse of rigidly and differentially rotating (Γ = 2) 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 [354] it was shown that the criterion for black-hole formation depends strongly on the amount of angular momentum, but Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 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 [369] 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, [110] performed long-term simulations of the gravitational collapse of rapidly rotating neutron stars to black holes in axisymmetry, employing the code described in [109] (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 J/M 2 < 1, J being the angular momentum and M the mass of the system. This result was also independently obtained by [356] 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, [110] 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 [419] 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 [378], where an axisymmetric Einstein-perfect fluid code is presented, and in [379] 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 [110], 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 lowdensity 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.

Three-dimensional simulations.
In full 3D, hydrodynamic simulations of the collapse of supermassive (Γ = 2) uniformly-rotating neutron stars to rotating black holes were first presented in [359], using the corresponding code discussed in Section 4.4. The simulations show no evidence for massive disk formation or outflows, which can be related to the moderate initial compactness of the stellar models (R/M ∼ 6). A proof-of-principle of the capabilities of the gr astro Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 67 code, mentioned in Section 4.4, to study black-hole formation was presented in [131], where the gravitational collapse of a spherical, unstable, relativistic polytrope was discussed. Similar tests with differentially rotating polytropes are given in [109] for an artificial viscosity-based, threedimensional general-relativistic hydrodynamics code, and test 3+1 simulations of the gravitational collapse of rapidly rotating neutron stars to black holes are mentioned in [110].
Supermassive black-hole formation was also studied by [437,438] using the cactus/whisky code. Alternative paths for formation were considered by following the general-relativistic evolution of a particular differentially-rotating N = 3 (Γ = 4/3) polytrope with such strong differential rotation as to have an initially toroidal shape. The study demonstrated the onset of nonaxisymmetric dynamic instabilities of supermassive stars in full general relativity. More precisely, it was found that the polytrope is unstable to nonaxisymmetric modes, which leads to a fragmentation into self-gravitating, collapsing components. Zink et al. [437] used adaptive-mesh-refinement techniques to follow the evolution to the formation of an apparent horizon centered on one of the fragments, signalling the formation of a black hole. A comprehensive study of the dependence of the fragmentation instability on parameters of the equilibrium model was reported in [438].
In [27] 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. Longterm 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 [28]. Despite good qualitative agreement being found with waveforms obtained in previous axisymmetric simulations [386], the newly computed amplitudes are about an order of magnitude smaller due to the more realistic rotation rates of the collapse progenitors. Figure 9: Still from a movie showing Gravitational radiation from the three-dimensional collapse of a neutron star to a rotating black hole [30]. The figure shows a snapshot in the evolution of the system once the black hole has been formed and the gravitational waves are being emitted in large amounts. See [30] for further details. Visualization made in AEI/ZIB. (To watch the movie, please go to the online version of this review article at http://www.livingreviews.org/lrr-2008-7.) More recently, [30] 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 ap-Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 proaches, 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 [30] 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 Q 20 . Moreover, in the central parts of the figure the collapse of the lapse function and the formation of a rotating black hole are visible. [30] 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.

Accretion on to black holes
The study of relativistic accretion and black-hole astrophysics is a very active field of research, both theoretically and observationally (see, e.g., [47,259] and references therein). On the one hand, advances in satellite instrumentation, such as the Rossi X-Ray Timing Explorer (RXTE), the Advanced Satellite for Cosmology and Astrophysics (ASCA), the Chandra X-ray Observatory, or the X-ray Multi-Mirror Mission-Newton, are great at stimulating and guiding theoretical research in accretion physics. The discovery of kHz quasi-periodic oscillations in low-mass X-ray binaries extends the frequency range over which these oscillations occur into timescales associated with the innermost regions of the accretion process (for a review, see [402]). Moreover, in extragalactic sources, spectroscopic evidence (broad iron emission lines) increasingly points to (rotating) black holes being the accreting central objects [394,207,62]. Important evidence has recently been collected thanks to the analysis of observations of relativistic X-ray emission lines from inner accretion disks around black holes carried out with Chandra, XMM-Newton, and Suzaku. Thick accretion disks (or tori) are orbiting the central black holes of many astrophysical objects such as quasars and other active galactic nuclei (AGNs), some X-ray binaries, and microquasars. In addition, they are believed to exist at the central engine of cosmic GRBs.
Disk accretion theory is primarily based on the study of (viscous) stationary flows and their stability properties through linearized perturbations thereof. A well-known example is the solution consisting of isentropic constant-angular-momentum tori, unstable to a variety of nonaxisymmetric global modes, discovered by Papaloizou and Pringle [313] (see [32] for a review of instabilities in astrophysical accretion disks). Since the pioneering work by Shakura and Sunyaev [348], thin disk models, parameterized by the αviscosity, in which the gas rotates with Keplerian angular momentum, which is transported radially by viscous stress, have been applied successfully to many astronomical objects. The thin disk model, however, is not valid for high luminosity systems, as it is unstable at high mass-accretion rates. In this regime, Abramowicz et al. [2] found the slim disk solution, which is stable against viscous and thermal instabilities. More recently, much theoretical work has been devoted to the problem of slow accretion, motivated by the discovery that many galactic nuclei are under-luminous (e.g., NGC 4258). Proposed accretion models involve the existence of advection-dominated accretion flows (ADAF solution; see, e.g., [285,284]) and advection-dominated inflow-outflow solutions (ADIOS solution [48]). The importance of convection for low values of the viscosity parameter α has also been discussed in the convection-dominated accretion flow (CDAF solution; see [179] and references therein). The importance of magnetic fields and their consequences for the stability properties of this solution are critically discussed in [31]. Magnetic effects have also been long advocated in connection with outflows from blackhole-accretion-disk systems, particularly when the mass accretion rate is much larger than the Eddington rate. Recent Chandra observations of the stellar-mass black-hole binary GRO J1655-40 [260] provide strong evidence supporting the magnetic nature of disk accretion onto black holes and its associated magnetic-powered winds. Theoretically, MHD outflows (and relativistic jets) are believed to be driven by energy flow through the accretion disk associated with the magnetic torque. The most convincing argument for the existence of the torque is the presence of the MRI, which not only drives the turbulent transport of angular momentum within the disk, self-regulating the accretion process, but also amplifies the strength of a seed vertical field on a dynamic timescale given by the inverse of the angular frequency of the disk. Extensive numerical work has been carried out in the last few years for a better understanding of such complex dynamics (see below).
For a wide range of accretion problems, the Newtonian theory of gravity is adequate for the description of the background gravitational forces (see, e.g., [191]). Extensive experience with Newtonian astrophysics has shown that explorations of the relativistic regime benefit from the use of model potentials. In particular, the Paczyński-Wiita pseudo-Newtonian potential for a Schwarzschild black hole [307] gives approximations of general-relativistic effects with an accuracy of 10 -20% outside the marginally stable radius (at r = 6 GM/c 2 ). Nevertheless, for comprehensive numerical work, a three-dimensional formalism is required, that is also able to cover the maximallyrotating black hole. In rotating spacetimes the gravitational forces cannot be captured fully with scalar potential formalisms. Additionally, geometric regions such as the ergosphere would be very hard to model without a metric description. Whereas the bulk of emission occurs in regions with almost Newtonian fields, only the observable features attributed to the inner region (as the relativistic X-ray emission lines) may crucially depend on the nature of the spacetime.
In the following we present a summary of representative time-dependent accretion simulations in relativistic hydrodynamics and MHD. We concentrate on multidimensional simulations. In the limit of spherical accretion, exact stationary solutions exist for both Newtonian gravity [57] and relativistic gravity [255]. Such solutions are commonly used to calibrate time-dependent hydrodynamic codes, by analyzing whether stationarity is maintained during a numerical evolution [171,237,117,339,36]. Correspondingly, similar tests have been proposed for the GRMHD equations, such as the spherically-symmetric accretion solution of a perfect fluid onto a Schwarzschild black hole in the presence of an (unphysical) radial magnetic field (b α = (b t , b r , 0, 0)) or the equatorial magnetized inflow solution in the region between a Kerr black-hole horizon and the marginally stable orbit, derived by [393] (see also [149]). These analytic solutions have been used in the literature to assess numerical codes [149,86,108,24,91].

Disk accretion
Pioneering numerical efforts in the study of black-hole accretion [411,171,167,168] made use of the frozen star paradigm of a black hole. In this framework, the time "slicing" of the spacetime is synchronized with that of asymptotic observers far from the black hole. Within this approach, Wilson [411] first investigated numerically the time-dependent accretion of inviscid matter onto a rotating black hole. This was the first problem to which his formulation of the hydrodynamic equations, as presented in Section 2.1.2, was applied. Wilson used an axisymmetric hydrodynamic code in cylindrical coordinates to study the formation of shock waves and X-ray emission in the strong-field regions close to the black-hole horizon, and was able to follow the formation of thick accretion disks during the simulations.
Wilson's formulation has been extensively used in time-dependent numerical simulations of thick disk accretion. In a system formed by a black hole surrounded by a thick disk, the gas flows in an effective (gravitational plus centrifugal) potential, whose structure is comparable to that of a close binary. The Roche torus surrounding the black hole has a cusp-like inner edge located at Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 the Lagrange point L 1 , where mass transfer driven by the radial pressure gradient is possible [3]. In [171] (see also [167]) Hawley and collaborators studied, in the test-fluid approximation and axisymmetry, the evolution and development of nonlinear instabilities in pressure-supported accretion disks formed as a consequence of the spiraling infall of fluid with some amount of angular momentum. The code used explicit second-order finite-difference schemes with a variety of choices to integrate the transport terms of the equations (i.e., those involving changes in the state of the fluid at a specific volume in space). The code also used a staggered grid (with scalars located at cell centers and vectors at the cell boundaries) for its suitability to difference the transport equations. Discontinuous solutions were stabilized with artificial viscosity terms.
With a three-dimensional extension of the axisymmetric code of Hawley et al. [170,171], Hawley [168] studied the global hydrodynamic nonaxisymmetric instabilities in thick, constant-angularmomentum accretion-gas tori orbiting around a Schwarzschild black hole. Such simulations showed that the Papaloizu-Pringle instability saturates in a strong spiral pressure wave, not in turbulence. In addition, the simulations confirmed that accretion flows through the torus could reduce and even halt the growth of the global instability. Extensions to Kerr spacetimes are reported in [84].
Igumenshchev and Beloborodov [180] performed two-dimensional relativistic hydrodynamic simulations of inviscid transonic disk accretion onto a Kerr black hole. The hydrodynamic equations follow Wilson's formulation, but the code avoids the use of artificial viscosity. The advection terms are evaluated with an upwind algorithm that incorporates the PPM scheme [80] to compute the fluxes. Their numerical work confirms analytic expectations regarding the strong dependence of the structure of the innermost disk region on the black-hole spin, and the functional form of the mass accretion rate with the potential barrier ∆W between the inner edge of the disk and the cusp,Ṁ ∝ (∆W ) Γ/(Γ−1) , with Γ being the adiabatic index.
Thick accretion disks orbiting black holes, on the other hand, may be subject to the runaway instability, as first proposed by Abramowicz et al. [1]. Starting with an initial disk filling its Roche lobe, so that mass transfer is possible through the cusp located at the L 1 Lagrange point, two evolutions are feasible when the mass of the black hole increases: either (i) the cusp moves inwards towards the black hole, which slows down the mass transfer, resulting in a stable situation, or (ii) the cusp moves deeper inside the disk material. In the latter case, the mass transfer speeds up, leading to runaway instability. This instability, whose existence is still a matter of debate (see, e.g., [129] and references therein), is an important issue for most models of cosmic GRBs, where the central engine responsible for the initial energy release is such a system consisting of a thick disk surrounding a black hole.
In [129] Font and Daigne presented time-dependent simulations of the runaway instability of constant-angular-momentum thick disks around black holes. The study was performed using a fully-relativistic hydrodynamics code based on HRSC schemes and the conservative formulation discussed in Section 2.1.3. The self-gravity of the disk was neglected and the evolution of the central black hole was assumed to be that of a sequence of Schwarzschild black holes of varying mass. The black-hole mass increase is determined by the mass accretion rate across the event horizon. In agreement with previous studies based on stationary models, [129] found that by allowing the mass of the black hole to grow, the disk becomes unstable. For all disk-to-blackhole mass ratios considered (between 1 and 0.05), the runaway instability appears very fast on a dynamic timescale of a few orbital periods (1 → 100), typically a few 10 ms and never exceeding 1 s, for a 2.5 M black hole and a large range of mass fluxes (ṁ ≥ 10 −3 M /s). Extensions of this work to marginally stable (or even stable) constant angular momentum disks were reported in Zanotti et al. [334] (animations can be visualized at the website listed in [330]).
An example of the simulations performed by [129] appears in Figure 10. This figure shows eight snapshots of the time-evolution of the rest-mass density, from t = 0 to t = 11.8 t orb . The contour levels are linearly spaced with ∆ρ = 0.1ρ 0 c , where ρ 0 c is the maximum value of the density at the center of the initial disk. In Figure 10 one can clearly follow the transition from a quasi-stationary Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Figure 10: Runaway instability of an unstable thick disk: contour levels of the rest-mass density ρ plotted at various times from t = 0 to t = 11.80 t orb , once the disk has almost entirely been destroyed. See [129] for details.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 accretion regime (panels (1) to (5)) to the rapid development of the runaway instability in about two orbital periods (panels (6) to (8)). At t = 11.80 t orb , the disk has almost entirely disappeared inside the black hole whose size has, in turn, noticeably grown.
On the other hand, recent simulations with nonconstant angular momentum disks and rotating black holes [128,82], show that the runaway instability is strongly suppressed in such cases. In the models built and evolved by [128,82] the angular momentum of the disks is assumed to increase outwards with the radial distance according to a power law, l = Kr α , where K and α are constant, such that K > 0 corresponds to prograde disks (with respect to the black-hole rotation) and K < 0 to retrograde disks. The results of [82] show that nonconstant angular momentum disks are dramatically stabilized for very small values of the angular momentum slope α, much smaller than the Keplerian value α = 1/2. In light of these results further ingredients need to be incorporated into the codes to finally discern the likelihood of the instability, namely magnetic fields and self-gravity.
The impact of magnetic fields on the dynamics of accretion disks around black holes has long been recognized, but only recently have the first systematic GRMHD simulations become possible [86,87,85,149,200,173,24,202,251,368,256,269]. We note, however, that as early as 1977 Wilson [412] developed the first numerical GRMHD scheme to study the axisymmetric evolution of plasma near a Kerr black hole. A key idea under intense numerical scrutiny currently is that the combined presence of a magnetic field and of a differentially rotating Keplerian disk, can lead to the magneto-rotational instability [32], generating effective viscosity and transporting angular momentum outwards through MHD turbulence. The dynamo generated by the MRI amplifies the magnetic field on dynamic timescales and, in turn, the associated magnetic torque is believed to be responsible for the generation of strong MHD outflows (see below). The sustained level of activity in the numerical modelling of magnetized accretion disks has also been accompanied by advances in the analytic front. In this respect, Komissarov [202] has derived an analytic solution for an axisymmetric, stationary, barotropic torus with constant distribution of specific angular momentum and a toroidal magnetic field configuration, which can be used to test multidimensional GRMHD codes in the strong gravity regime.
Yokosawa [429,430], using Wilson's formulation, studied the structure and dynamics of relativistic accretion disks and the transport of energy and angular momentum in MHD accretion on to a rotating black hole. In his code, the hydrodynamic equations are solved using the FCT scheme [60] (a second-order flux-limiter method that avoids oscillations near discontinuities by reducing the magnitude of the numerical flux), and the magnetic induction equation is solved using the constrained transport method [119]. The code contains a parameterized viscosity based on the α-model [348]. The simulations revealed different flow patterns inside the marginally-stable orbit, depending on the thickness, h, of the accretion disk. For thick disks with h ≥ 4r h , r h being the radius of the event horizon, the flow pattern becomes turbulent.
The evolution of the MRI in both the Schwarzschild and Kerr metrics has been numerically investigated, either in axisymmetry or in three dimensions, in a number of recent papers [86,149,87,173]. The simulations have been performed with two of the codes listed in Table 2, namely HARM [149] and the code of De Villiers and Hawley [86]. We note that both codes not only implement different formulations of the GRMHD equations but they also use entirely different numerical schemes (conservative and incomplete Riemann solver the former, nonconservative and artificial viscosity the latter). Furthermore, the code of [149] uses Kerr-Schild coordinates (modified to increase the resolution where the disk lies), regular at the event horizon (see below), while the code of [86] employs the more traditional Boyer-Lindquist coordinate system (singular at the horizon). We note that the initial models in all these investigations are purely hydrodynamic, constant-angular-momentum, thick disks in equilibrium, with an ad hoc poloidal magnetic field distribution superposed. The simulations, carried out for a large sample of parameter values (such as black-hole spin or initial magnetization of the plasma) reveal the following basic evolution of Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 the accretion process: first, there is an initial rapid growth of the toroidal magnetic field driven by shear. This is followed by the development of the MRI until saturation. The accretion process then continues through the quasi-stationary evolution of the disk by sustained MHD turbulence, which redistributes the density and angular momentum of the disk. An example of such an evolution is depicted in Figure 11, taken from [149]. This figure shows the first (left) and last snapshot (after about 8 orbital periods) of the logarithm of the density field for a magnetized, constant-angularmomentum torus around a Kerr black hole with a/M = 0.5. The exponential growth of the MRI happens within the first orbital period, after which the magnetic field has been amplified to high enough levels to distort the torus and initiate the accretion process.
De Villiers et al. [87] are able to identify five generic features in the accretion flow: a turbulent and dense disk body with a roughly constant H/R (H and R being the pressure scale height and radius, respectively), an inner disk consisting of an inner torus and a plunging inflow, a magnetized coronal envelope, an evacuated axial funnel, and a jet-like outflow along the funnel wall. Overall, they found a decrease in the accretion rate into the black hole with increasing spin parameter. The accretion rate is highly variable in all models, which makes the characteristics of the inner torus highly variable as well, which plays a key role in the accretion flow into the black hole. The outflows are found to be launched from a region in the coronal envelope above the inner torus, which is closer to the horizon as the black-hole spin increases, resulting in more powerful jets.
Three-dimensional numerical work has also addressed recently the evolution of accretion disks that are misaligned (tilted) with respect to the rotation axis of a Kerr black hole. Studies have been carried out both for hydrodynamic flows [139] and for magnetized flows [140]. These studies are motivated by the increasing evidence that misaligned black holes may actually exist in nature, evidence collected both from observational data as from theoretical arguments (see [140] and references therein). The simulations have been done using the cosmos and cosmos++ codes listed in Tables 1 and 2, respectively. For hydrodynamic flows it was found that Lense-Thirring precession causes the disks to warp differentially, yet the precession timescale is not shorter than other dynamic timescales. The late-time evolution of the disks showed solid-body precession at rates consistent with some low frequency quasi-periodic oscillations (QPOs). Similarly, the MRI turbulent tilted-disk simulations have not shown indication of a Bardeen-Petterson effect at large.
Recent work on disk accretion has also focused on improving the microphysics of the tori, Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 as in the axisymmetric magnetohydrodynamic simulations of [368] for neutrino-cooled accretion tori around rotating black holes. (See also [256] for comparisons between a simple ideal fluid EOS and the relativistic Synge EOS in GRMHD accretion flows.) The EOS used takes into account neutronization, the nuclear statistical equilibrium of a gas of free nucleons and α-particles, black body radiation, and a relativistic Fermi gas (neutrinos, electrons, and positrons), along with several cooling processes. It is found that the neutrino luminosity is ∼ 10 53 erg/s for a torus of mass ≥ 0.2 M and a black hole of 4 M , irrespective of its spin, a figure compatible with current estimates in short-hard GRBs. However, the luminosity decreases with decreasing torus mass, and, in light of recent state-of-the-art simulations of neutron-star binary mergers (see below), the formed black-hole-torus systems have generally small torus masses ( 0.1 M ). On the positive side, however, the simulations predict a strong variability in the neutrino luminosity (on timescales of a few ms), due to the MRI-turbulent accretion flow, which could help explain the observed variability of GRB light curves.

Jet formation
The most promising processes for producing relativistic jets like those observed in AGNs, microquasars, and GRBs involve the hydromagnetic centrifugal acceleration of material from the accretion disk [49], or the extraction of rotational energy from the ergosphere of a Kerr black hole by magnetic processes [316,50]. The Blandford-Payne model [49] is based upon the existence of a large-scale poloidal magnetic field passing through the disk, whose matter dynamics is dominated by magnetic tension. Following angular-momentum-conservation considerations, accelerated outflows are produced whenever the magnetic field lines have a sufficient angle with respect to the disk. Correspondingly, in the Blandford-Znajek mechanism [50] the black hole is seen as a magnetized rotating conductor whose rotational energy can be efficiently extracted by means of a magnetic torque, giving rise to a Poynting flux jet. On the other hand, in the magnetic Penrose process [316] the magnetic field lines across the ergosphere are twisted by frame dragging. The twist of the lines propagates outwards as a torsional Alfvén wave train carrying electromagnetic energy and making the total energy of the plasma near the black hole decrease to negative values. In addition, the accretion of this plasma by the black hole reduces the black-hole rotational energy. As in the Blandford-Znajek process, the type of outflows generated by this process is also a Poynting flux jet. The viability of all these various mechanisms has been investigated numerically in the last few years.
Koide et al. performed the first MHD simulations of jet formation in general relativity [196,195,197,194,288]. Their code uses the 3+1 formalism of general relativistic conservation laws of particle number, momentum, and energy, and Maxwell equations with infinite electric conductivity. The MHD equations are numerically solved in the test-fluid approximation (in the background geometry of Kerr spacetime) using a finite difference symmetric scheme [83]. The Kerr metric is described in Boyer-Lindquist coordinates, with a radial tortoise coordinate to enhance the resolution near the horizon.
In [197,194], in particular, the GRMHD behavior of a plasma flowing into a rapidly-rotating black hole (a = 0.99995) in a large-scale magnetic field was investigated numerically. The initial magnetic field is uniform and strong, B 0 = 10 ρ 0 c 2 , ρ 0 being the mass density. The initial Alfvén speed is v A = 0.953 c. The simulation shows how the rotating black hole drags the inertial frames around in the ergosphere. The azimuthal component of the magnetic field increases because of azimuthal twisting of the magnetic field lines, as depicted in Figure 12. This frame-dragging dynamo amplifies the magnetic field to a value which, by the end of the simulation, is three times larger than the initial one. The twist of the magnetic field lines propagates outwards as a torsional Alfvén wave. The magnetic tension torques the plasma inside the ergosphere in a direction opposite to that of the black-hole rotation. Therefore, the angular momentum of the plasma outside receives Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 75 Figure 12: Jet formation: the twisting of magnetic field lines around a Kerr black hole (black sphere). The yellow surface is the ergosphere. The red tubes show the magnetic field lines that cross into the ergosphere. Figure taken from [197] (used with permission).
Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 a net increase. Even though the plasma falls into the black hole, electromagnetic energy is ejected along the magnetic field lines from the ergosphere, due to the propagation of the Alfvén wave. By total-energy-conservation arguments, Koide et al. [197] conclude that the ultimate result of the generation of an outward Alfvén wave is the magnetic extraction of rotational energy from the Kerr black hole, a process the authors call the MHD Penrose process (but see [201]).
The first time-dependent GRMHD simulations of the magnetically-dominated monopole magnetospheres of black holes were performed by [200], who found that the numerical solution evolves towards a stable steady-state solution, which is very close to the corresponding force-free solution found by Blandford & Znajek [50]. A similar result was reported by [252] for weakly magnetized accretion disks around Kerr black holes, who found numerical solutions consistent with the Blandford-Znajek solution in the low-density funnel region around the black hole. The numerical simulations of [200] showed for the first time the development of an ultrarelativistic particle wind, Poynting dominated, from a rotating black hole.
Further numerical work on the Blandford-Znajek mechanism was pursued by Komissarov in [201] using initial models similar to those employed by Koide et al. [197,194]. The long-term solution shows properties that are significantly different from those of the initial transient phase studied by [197,194]. In particular, no regions of negative hydrodynamic 'energy at infinity' are seen inside the ergosphere and the MHD Penrose process does not operate. Similar conclusions were drawn by [252], which highlights the fact that the models of [197,194] were evolved for too short a time to observe unbound mass outflows along the funnel region (but see [288] for longer evolutions). The results of Komissarov [201] indicate that the rotational energy of the black hole continues to be extracted via the purely electromagnetic Blandford-Znajek mechanism. This effect also operates when the black-hole spin is the maximum allowed, a/M = 1 (extreme Kerr), as has been recently shown by [206]. Perhaps the most important result of these axisymmetric models has been not finding strong relativistic outflows from the black hole, contrary to observational data indicating the existence of strongly relativistic jets in numerous scenarios (but see below).
Ultrarelativistic outflows have not been found either in the comprehensive study of Hawley et al. [87,173,88,210,169]. Their third paper on the series [88], as well as [169], are devoted to studying the formation of unbound outflows in the simulations. Such unbound outflows, which are produced self-consistently from accretion disks that do not have large-scale magnetic fields as an initial condition, do not require rotation in the black hole, yet their strength is greatly enhanced by the black-hole spin. If the spin of the black hole is high, as in their almost extreme Kerr (Kepler disk) model KDJ with a/M = 0.99, the energy ejected in the Poynting-flux-dominated outflows can be tens of percent of the accreted rest mass. At low spin, kinetic energy and enthalpy of the matter dominate the outflow energetics. These basic results are in qualitative agreement with axisymmetric simulations performed by [252] using the conservative, shock-capturing scheme HARM. The jets formed in the three-dimensional simulations have two major components: a matter-dominated outflow, which moves at moderate speed (v/c ∼ 0.3) along the centrifugal barrier surrounding an evacuated axial funnel, and a highly-relativistic, Poynting-flux-dominated jet within the funnel. This jet is accelerated and collimated by magnetic and gas pressure forces in the inner torus and the surrounding corona. A critical discussion of the physical origin of the relativistic Poynting jet formed in these numerical simulations has been conducted by Punsly [325], who points out the importance of incorporating resistive MHD reconnection in the modelling to gain further insight.
Axisymmetric models of jet formation allow one to extend the computational time of the simulations at an affordable cost. This is the case for the simulations reported by [251], which extended the earlier work of [252] by studying jets from a collapsar GRB model (see below) until t ∼ 10 4 GM/c 3 and out to radial distances of r ∼ 10 4 GM/c 2 . For this study a new version of the GRMHD code HARM was used. It was found that, at such large spatial and temporal scales, the Poynting-flux-dominated jet is accelerated by continuous mass-loading from the disk to reach bulk Lorentz factors of ∼ 10 and maximum terminal Lorentz factors of ∼ 10 3 and it collimates to narrow Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity 77 half-opening angles of 5 • . This remarkable result is in concordance with the values observed in jets associated with GRBs, AGN, or X-ray binary systems.
At stellar scales, relativistic thick disks (or tori) orbiting around stellar-mass black holes may form after the gravitational collapse of the core of a rotating massive star (M > 25 M ) or after a neutron-star-binary merger. Such an accreting torus-plus-black-hole system is the most likely candidate to be the mechanism powering gamma-ray bursts. Part of the large amounts of energy released by accretion is deposited in the low-density region along the rotation axis, where the heated gas expands in a jet-like fireball. In connection with GRBs, van Putten and Levinson [405] have considered the theoretical evolution of an accretion torus in suspended accretion. These authors claim that the formation of baryon-poor outflows may be associated with a dissipative gap in a differentially-rotating magnetic-flux tube supported by an equilibrium magnetic moment of the black hole. Numerical simulations of nonideal MHD, incorporating radiative processes, are necessary to validate this picture. Numerical simulations of relativistic jets propagating through progenitor stellar models of collapsars have been presented in [9,435] (see also [266] for GRMHD simulations of jet formation in such context). The collapsar scenario, proposed by [421], is currently the most favored model for explaining long duration GRBs. The hydrodynamic simulations performed by [9] employ the three-dimensional code genesis [7] with a 2D spherical grid and equatorial-plane symmetry. The gravitational field of the black hole is described by the Schwarzschild metric, and the relativistic hydrodynamic equations are solved in the test-fluid approximation using a Godunov-type scheme. Aloy et al. [9] showed that the jet, initially formed by an ad hoc energy deposition of a few 10 50 erg/s within a 30 • cone around the rotation axis, reaches the surface of the collapsar progenitor intact, with a maximum Lorentz factor of ∼ 33.
Similarly, the genesis code has also been used to perform hydrodynamic simulations of rel-Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 ativistic jets in the neutron-star-binary merger scenario. Figure 13 shows a snapshot at time t = 0.05 s of an axisymmetric simulation of the launch and evolution of a relativistic outflow driven by energy deposition through neutrino-antineutrino annihilation in the vicinity of a blackhole-plus-accretion-torus system (model B1 of [8]). The left panel displays the logarithm of the Lorentz factor in the entire computational domain, while the right panel focuses on the innermost regions closer to the central black hole, depicting the logarithm of the rest-mass density. Since in this simulation the system is the remnant of a neutron-star-binary merger, it naturally provides a baryon-free channel along the rotation axis through which the outflow can reach ultrarelativistic velocities. Indeed, terminal Lorentz factors of about 1000 are attained for initial energy deposition rates of some 10 48 erg/s per stereo-radian. As shown in Figure 13 the outflow is strongly collimated due to the presence of the accretion torus, which, in turn, makes the gamma-ray emission very weak outside of the polar cones. The calculations reported by [8] predict that about one out of a hundred neutron-star mergers produces a detectable GRB, as its ultrarelativistic jet happens to be sent along the line of sight 2 . It is worth noting that the simulations performed by [8] have also been used by [11] to illustrate a new effect in relativistic hydrodynamics that can act as an efficient booster in jets and help explain the observed ultrarelativistic speeds. This effect, concealed in the subtleties of the relativistic Riemann problem, is purely hydrodynamic and occurs when large velocities tangential to a discontinuity are present in the flow, yielding Lorentz factors ∼ 10 2 − 10 3 or larger in flows with moderate initial Lorentz factors.

Wind accretion
The term "wind" or hydrodynamic accretion refers to the capture of matter by a moving object under the effects of the underlying gravitational field. The canonical astrophysical scenario in which matter is accreted in such a nonspherical way was suggested originally by Bondi, Hoyle, and Lyttleton [175,58], who studied, using Newtonian gravity, the accretion on to a gravitating point mass moving with constant velocity through a nonrelativistic gas of uniform density. The matter flow inside the accretion radius, after being decelerated by a conical shock, is ultimately captured by the central object. Such a process applies to describe mass transfer and accretion in compact X-ray binaries, in particular in the case in which the donor (giant) star lies inside its Roche lobe and loses mass via a stellar wind. This wind impacts on the orbiting compact star forming a bow-shaped shock front around it. This process is also believed to occur during the common envelope phase in the evolution of a binary system.
Numerical simulations have extended the simplified analytic models and have helped to develop a thorough understanding of the hydrodynamic accretion scenario, in its fully three-dimensional character (see, e.g., [341,42,125] and references therein). The numerical investigations have revealed the formation of accretion disks and the appearance of nontrivial phenomena such as shock waves and tangential (flip-flop) instabilities.
Most of the existing numerical work has used Newtonian hydrodynamics to study the accretion onto nonrelativistic stars [341]. More recently, Newtonian MHD studies have also investigated the spin-down of moving magnetars [399]. For compact accretors such as neutron stars or black holes, the correct numerical modeling requires a general-relativistic (magneto-)hydrodynamic description. Within the relativistic, frozen star framework, wind accretion onto "moving" black holes was first studied in axisymmetry by Petrich et al. [317]. In this work, Wilson's formulation of the hydrodynamic equations was adopted. The integration algorithm was borrowed from [387] with the transport terms finite-differenced following the prescription given in [171]. An artificial viscosity term of the form Q = aρ(∆v) 2 , with a a constant, was added to the pressure terms of the equations in order to stabilize the numerical scheme in regions of sharp pressure gradients. An extensive survey of the morphology and dynamics of relativistic wind accretion past a Schwarzschild black hole was performed by [133,132]. This investigation differs from [317] in both the use of a conservative formulation for the hydrodynamic equations and the use of Godunov-type HRSC schemes. Axisymmetric computations were compared to [317], finding major differences in the shock location, opening angle, and accretion rates of mass and momentum. The reasons for the discrepancies are related to the use of different formulations, numerical schemes, and grid resolution, much higher in [133,132]. Nonaxisymmetric two-dimensional studies, restricted to the equatorial plane of the black hole, were discussed in [132], motivated by the nonstationary patterns found in Newtonian simulations (see, e.g., [42,125]). The relativistic computations revealed that initially-asymptotic uniform flows always accrete onto the black hole in a stationary way that closely resembles the previous axisymmetric patterns.
In [308], Papadopoulos and Font presented a procedure that simplifies the numerical integration of the general-relativistic hydrodynamics equations near black holes. This procedure is based on identifying classes of coordinate systems in which the black-hole metric is free of coordinate singularities at the horizon (unlike the commonly adopted Boyer-Lindquist coordinates), independent of time, and admits a spacelike decomposition. With those coordinates the innermost radial boundary can be placed inside the horizon, allowing for an unambiguous treatment of the entire (exterior) physical domain. In [135,136] this approach was applied to the study of relativistic wind accretion onto rapidly-rotating black holes. The effects of the black-hole spin on the flow morphology were found to be confined to the inner regions of the black-hole potential well. Within this region, the black-hole angular momentum drags the flow, wrapping the shock structure around. An illustrative example is depicted in Figure 14. The left panel of this figure corresponds to a simulation employing the Kerr-Schild form of the Kerr metric, regular at the horizon. The right panel shows how the accretion pattern would look were the computation performed using the more common Boyer-Lindquist coordinates. The transformation induces a noticeable wrapping of the shock around the central black hole. The shock would wrap infinitely many times before reaching the horizon. As a result, the computation in these coordinates would be much more challenging than in Kerr-Schild coordinates. We note that such "horizon-penetrating" coordinates (or variations thereof) are currently employed in most of the state-of-the-art codes developed to study accretion onto black holes [149,201,140,396].

Gravitational radiation
Semi-analytical studies of finite-sized collections of dust, shaped in the form of stars or shells and falling isotropically onto a black hole are available in the literature [282,165,350,302,318]. These early investigations approximate gravitational collapse by a dust shell of mass m falling into a Schwarzschild black hole of mass M m. These studies have shown that for a fixed amount of infalling mass, the gravitational radiation efficiency is reduced compared to the point particle limit (∆E ∼ 0.0104 m 2 /M ), due to cancellations of the emission from distinct parts of the extended object. In [310], such conclusions were corroborated with numerical simulations 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 nonlinear evolution was integrated using a (nonlinear) hydrodynamics code (adopting the conservative formulation of Section 2.1.3 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., [77]). In the framework of the Newman-Penrose formalism, the equations for the perturbed Weyl tensor components decouple, and when written in the frequency domain, they even separate [397]. Papadopoulos and Font [310] used the limiting case for Schwarzschild black holes, i.e., the inhomogeneous Bardeen-Press equation [37]. The simulations showed the gradual excitation of the black-hole quasi-normal-mode frequency by sufficiently compact shells. Similar studies based on a metric formalism have recently been discussed in [275,276].
An example of the simulations of [310] appears in the movie of Figure 15. 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 parameterized according to ρ = ρ 0 + e −k(r * −r0) 2 sin 2 θ with k = 2, ρ 0 = 10 −2 and r 0 = 35 M . Additionally, r * 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 Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 the movie, one can clearly see how the burst of the emission coincides with the most dynamic accretion phase, when the shell crosses the peak of the potential and is subsequently captured by the black hole. The gravitational wave signal coincides with the quasinormal ringing frequency of the Schwarzschild black hole, 17 M . The existence of an initial burst, separated in time from the physical burst, is also noticeable in the movie. It just reflects the gravitational radiation content of the initial data (see [310] for a detailed explanation; this "junk" radiation is also a common inevitable feature in the recent black-hole-binary simulations [323]).
One-dimensional numerical simulations of a self-gravitating perfect fluid accreting onto a black hole are presented in [312], where the effects of mass accretion during the gravitational wave emission from a black hole of growing mass are explored. Using the conservative formulation outlined in Section 2.2.2 and HRSC schemes, Papadopoulos and Font [312] performed the simulations adopting an ingoing null foliation of a spherically-symmetric black-hole spacetime [311]. Such a foliation penetrates the black-hole horizon, allowing for an unambiguous numerical treatment of the inner boundary. The essence of nonspherical gravitational perturbations was captured by adding to the (characteristic) Einstein-perfect-fluid system, the evolution equation for a minimally-coupled massless scalar field. The simulations showed the familiar damped-oscillatory radiative decay, with both decay rate and frequencies being modulated by the mass-accretion rate. Any appreciable increase in the horizon mass during the emission reflects on the instantaneous signal frequency, f , which shows a prominent negative branch in theḟ (f ) evolution diagram. The features of the frequency evolution pattern reveal key properties of the accretion event, such as the total accreted mass and the accretion rate. To highlight the recent advances accomplished in numerical relativity it is worth it linking this result with the 3D simulations of black-hole formation performed by [30], which have provided information on the ring-down gravitational signal through which the dynamic black hole relaxes to its stationary state (see Section 5.1.2).
In a series of recent papers [434,433,269] it has also been shown that upon the introduction of perturbations, stable, relativistic, accretion tori manifest a long-term oscillatory behavior lasting for tens of orbital periods. The associated changes in the mass-quadrupole moment of such disks render these objects promising sources of high-frequency, detectable gravitational radiation for ground-based interferometers and advanced resonant bar detectors, particularly for Galactic systems and when the average disk density is close to nuclear matter density. If the disks are instead composed of low-density material stripped from the secondary star in Low-Mass X-Ray Binaries (LMXBs), their oscillations could help explain the high-frequency QPOs observed in the spectra of X-ray binaries. Indeed, such HFQPOs can be explained in terms of p-mode oscillations of a small-size torus orbiting around a stellar-mass black hole [332].
The studies reported in the existing papers have considered both Schwarzschild and Kerr black holes, as well as constant and nonconstant (power-law) distributions of the specific angular momentum of the disks. The most recent investigation [269] has accounted for magnetic fields in the tori, which, in addition to satisfying a polytropic EOS and having a constant distribution of the specific angular momentum, were built with a nonzero toroidal magnetic field component, for which equilibrium configurations exist [202]. A representative sample of initial models was built and the dependence of their dynamics on the strength of the magnetic field was investigated in long-term GRMHD evolutions in timescales of 100 orbital periods. Overall, the dynamics of the magnetized tori considered was found to be strikingly similar to that found in the purely hydrodynamic case by [434,433]. It should be noted that, since the axisymmetric initial models do not have a poloidal magnetic field component, the MRI, which could potentially alter the conclusions, can not develop.

Hydrodynamic and magnetohydrodynamic evolution of neutron stars
The numerical investigation of many interesting astrophysical processes involving neutron stars, such as the rotational evolution of proto-neutron stars (which can be affected by a dynamic barmode instability and by the Chandrasekhar-Friedman-Schutz instability), the appearance of MHD winds and outflows from pulsar magnetospheres, or the gravitational radiation from unstable pulsation modes or, more importantly, from the catastrophic coalescence and merger of neutron-star compact binaries, requires the ability of accurate, long-term hydrodynamic and MHD evolutions employing relativistic gravity. These scenarios are receiving increasing attention in recent years, as we discuss next.
In particular, the understanding of the magnetospheres of neutron stars has recently improved through time-dependent numerical simulations, which allow one to investigate the stability of steady-state analytic solutions of the pulsar equation. Existing attempts have used both ideal relativistic magnetohydrodynamics and resistive force-free electrodynamics (see [203] and references therein). The investigation by [203] permits one to rule out some stationary solutions, which are not naturally recovered by the time-dependent simulation. Further work is needed to account for resistive relativistic MHD computations.
Another remarkable investigation based on fully-relativistic MHD simulations by [205] has permitted us to understand the properties of the flow produced by a magnetized pulsar wind within a plerionic nebula. The results of the axisymmetric simulations have disclosed the complex dynamics of the post-shock flow, very different from the steady quasi-radial outflow assumed in earlier (spherically-symmetric) analytical models for plerions. In particular, the jet-torus pattern discovered by Chandra in the Crab nebula can be explained within the MHD approximation when the condition of spherical symmetry is no longer enforced, as shown by the simulated synchrotron X-ray images of the MHD numerical data.
General-relativistic MHD winds from rotating neutron stars have also been studied by [66]. Magnetically-dominated winds from stars are central to the angular momentum evolution of these objects, and GRMHD studies, such as those carried out by [66], may help us understand the observed spin-down of rotation-powered pulsars with strong magnetic fields.

Pulsations and instabilities of relativistic stars
The use of relativistic hydrodynamic codes to study the stability properties of neutron stars and to compute mode frequencies of oscillations of such objects has increased in recent years (see, e.g., the Living Reviews article by Stergioulas [388] and references therein). In the case of GRMHD codes, the first investigations are only just starting. An obvious advantage of the hydrodynamic approach is that it includes, by construction, nonlinear effects, which are important in situations where the linearized equations (commonly employed in calculations of the mode frequencies of pulsating stars) break down. In addition, while for nonrotating stars the frequencies of normal modes can be computed with perturbative methods and a theory of gravitational wave asteroseismology has already been formulated, there exist no accurate frequency determinations for rapidly rotating stars to date, nor has the theory of gravitational-wave asteroseismology been extended to include the effects of rotation on the oscillation frequencies. Due to the advanced status of current hydrodynamics codes, however, the computation of mode frequencies in rapidly rotating relativistic stars might be easier to achieve through nonlinear numerical evolutions than by using perturbative computations (see, e.g., the results in [131,355,364,389,99]).
Hydrodynamic evolutions of polytropic models of spherical neutron stars can be used as test-bed computations for multidimensional codes. Representative examples are the simulations by [157], with pseudo-spectral methods, and by [339] with HRSC schemes. These investigations adopted radial-gauge polar-slicing coordinates in which the general relativistic equations are expressed in a simple way that resembles Newtonian hydrodynamics. Gourgoulhon [157] used a numerical Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 code to detect, dynamically, the zero value of the fundamental mode of a neutron star against radial oscillations. Romero et al. [339] highlighted the accuracy of HRSC schemes by finding, numerically, a change in the stability behavior of two slightly different initial neutron star models: for a given EOS, a model with mass 1.94532 M is stable and a model of 1.94518 M is unstable. More recently, in [383] a method based on the nonlinear evolution of deviations from a background stationary-equilibrium star was applied to study nonlinear radial oscillations of a neutron star. The accuracy of the approach permitted a detailed investigation of nonlinear features such as quadratic and higher-order mode coupling and nonlinear transfer of energy.
Axisymmetric pulsations of rotating neutron stars can be excited in several scenarios, such as core collapse, crustquakes and corequakes, and binary mergers, and they could become detectable either via gravitational waves or high-energy radiation. An observational detection of such pulsations would yield valuable information about the EOS of relativistic stars. In recent years, the time evolution of the nonlinear equations governing the dynamics of matter and spacetime has been introduced as a promising new approach for computing mode frequencies [138,130,131,390,389,99]. For small amplitudes, the obtained frequencies are in excellent agreement with those expected by linear perturbation theory, while two-dimensional eigenfunctions can be obtained through a Fourier transform technique.
As a first step towards the study of pulsations of rapidly-rotating relativistic stars, Font, Stergioulas, and Kokkotas [138] developed an axisymmetric numerical code called ToniK (see Table 1) that integrates the GRHD equations in a fixed-background spacetime. The finite-difference code is based on a state-of-the-art approximate Riemann solver [102] and incorporates different second and third-order TVD and ENO numerical schemes. This code can accurately evolve rapidly rotating stars for many rotational periods, even for stars at the mass-shedding limit. The test simulations reported in [138] show that, for nonrotating stars, small amplitude oscillations have frequencies that agree to better than 1% with linear normal-mode frequencies (radial and nonradial) in the Cowling approximation (i.e., when the evolution of the spacetime variables is neglected). Axisymmetric modes of pulsating nonrotating stars are computed in [378], both in Cowling and in fully-coupled evolutions. Contrary to the 2+1 approach followed by [138], the code used in [378] evolves the relativistic stars on null spacetime foliations (see Section 2.2.2).
Until very recently (see below), the quasi-radial modes of rotating relativistic stars had been studied only under simplifying assumptions, such as in the slow-rotation approximation or in the relativistic Cowling approximation. An example of the latter is presented in [130], where a comprehensive study of all low-order = 0, 1, 2 and 3 axisymmetric modes of uniformly and rapidlyrotating relativistic stars was presented, using the code of [138]. This was done for a sequence of appropriately-perturbed stationary rotating stars, from the nonrotating limit to the mass-shedding limit. The frequencies of the axisymmetric modes are affected significantly by rotation only when the rotation rate exceeds about 50% of the maximum allowed. As expected, at large rotation rates, apparent mode crossings between different modes appear. A comparison of the mode-frequencies computed with ToniK was carried out by [188] with his pizaa code (see Table 1), which offers an improved treatment of quasistationary scenarios. The agreement found was satisfactory.
In [131], the first mode frequencies of uniformly-rotating stars in full general relativity and rapid rotation were obtained, using the three-dimensional code gr astro. Such frequencies were computed both in fixed spacetime evolutions and in coupled hydrodynamic and spacetime evolutions. The Cowling runs allowed a comparison with earlier results reported by [130], obtaining agreement at the 0.5% level. The fundamental mode frequencies and their first overtones obtained in fully coupled evolutions show a dependence on the increased rotation, which is similar to the one observed for the corresponding frequencies in the Cowling approximation [130].
Small-amplitude, nonlinear pulsations of both uniformly and differentially rotating neutron stars were studied by [389] employing the ToniK code. The centrifugal forces and the degree of differential rotation have significant effects on the mode eigenfunction and it was found that near Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 the mass-shedding limit, the pulsations are damped due to shocks forming at the surface of the star. This damping mechanism may set a small saturation amplitude for modes that are unstable to the emission of gravitational waves. It was also found that the fundamental quasi-radial mode is split, at least in the Cowling approximation and mainly in differentially rotating stars, into two different sequences. This study was revisited by [99] using the CoCoNuT code, which solves the GRHD equations for a conformally-flat three-metric (see Section 4.4). Differential rotation significantly shifts mode frequencies to smaller values, increasing the likelihood of detection by current gravitational wave interferometric detectors. An extended avoided crossing between the l = 0 and l = 4 first overtones was observed (previously known to exist from perturbative studies), which is important for correctly identifying mode frequencies in case of detection. For uniformly rotating stars near the mass-shedding limit, the existence of the mass-shedding-induced damping of pulsations was confirmed, even though the effect is not as strong as was previously found in the Cowling approximation [389].
Neutron stars following a core collapse supernova are rotating at birth and can be subject to various nonaxisymmetric instabilities (see, e.g., [388] for a review). Among those, if the rotation rate is high enough that the ratio of rotational kinetic energy T to gravitational potential energy W , β ≡ T /|W |, exceeds the critical value β d ∼ 0.27 inferred from studies with incompressible Maclaurin spheroids, the star is subject to a dynamic bar-mode (l = m = 2 f -mode) instability driven by hydrodynamics and gravity. Its study is highly motivated these days as such an instability bears important implications in the prospects of the detection of gravitational radiation from newlyborn rapidly-rotating neutron stars.
Newtonian simulations of the bar-mode instability from perturbed equilibrium models of rotating stars have shown that β d ∼ 0.27 and is quite independent of the stiffness of the EOS provided the star is not strongly differentially rotating. The relativistic simulations of [358] yield a value of β ∼ 0.24 − 0.25 for the onset of the instability, while the dynamics of the process closely resembles that found in Newtonian theory, i.e., unstable models with large enough β develop spiral arms following the formation of bars, ejecting mass and redistributing the angular momentum. As the degree of differential rotation becomes higher and more extreme, Newtonian simulations have also shown that β d can be as low as O(0.01).
Further relativistic hydrodynamic simulations of the dynamic bar-mode instability have been performed by [29] using the whisky code and extending the set of models considered by [358]. This study focused on investigating the persistence of the bar deformation once the instability has reached its saturation and on the precise determination of the parameter β d signalling the threshold for the onset of the instability. Generic nonlinear mode-coupling effects were found to appear during the development of the instability, which can severely limit the persistence of the bar deformation and even suppress the instability altogether. An animation of one such simulation can be seen in the movie of Figure 16.
A word of caution is needed here, as it remains unclear whether the requirements inferred from numerical simulations are at all met by core-collapse progenitors. As shown by [384] magnetic torques can spin down the core of the progenitor, which leads to slowly rotating neutron stars at birth (∼ 10 -15 ms). The most recent, state-of-the-art computations of the evolution of massive stars, which include angular momentum redistribution by magnetic torques and spin estimates of neutron stars at birth [172], lead to core-collapse progenitors, which do not seem to rotate fast enough to guarantee the unambiguous growth of the canonical bar-mode instability.
Regarding the CFS instability, relativistic hydrodynamic simulations of nonlinear r-modes are presented in [390], using a version of the gr astro code (see also [224] for Newtonian simulations). The gravitational radiation reaction-driven instability of the r-modes might have important astrophysical implications, provided the instability does not saturate at low amplitudes by nonlinear effects or by dissipative mechanisms. Regarding the latter, existing studies have shown that the mode amplitude could be limited by shear and bulk viscosity, by energy loss to a magnetic field Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 Figure 16: Still from a movie showing The time evolution of the bar-mode instability of one of the differentially rotating neutron-star models of [29]. The animation depicts the distribution of the rest-mass density. The exponential growth of the instability, which is followed by its saturation, the development of spiral arms, and the attenuation of the bar deformation, is all clearly visible in the movie. The last part of the dynamics is dominated by an almost axisymmetric configuration. Visualization developed at SISSA. Used with permission. (To watch the movie, please go to the online version of this review article at http://www.livingreviews.org/lrr-2008-7.) driven by differential rotation, by shock waves, or by the leak of the r-mode energy into some short wavelength oscillation modes (see [26] and references therein). The latter mechanism would dramatically limit the r-mode amplitude to a small value of O(10 −3 ), much smaller than those found in the early simulations of [390,224] for isentropic stars (see [388] for a complete list of references on the subject). Energy leak of the r-mode into other fluid modes was also considered by [160] through Newtonian hydrodynamic simulations, finding a catastrophic decay of the amplitude only once it has grown to a value larger than that reported by [26]. Detailed analysis showed that this breakdown is due to nonlinear 3-mode coupling of the r-mode to two other inertial modes. The saturation amplitude is of O(10 −2 ). It is believed that saturation amplitudes of this order are still appropriate to detect gravitational waves from the r-mode instability in LMXBs.
As mentioned in Section 5.1.2 the axisymmetric collapse of differentially-rotating magnetized neutron stars has been studied by [105,360] in full general relativity. The effects of magnetic fields on the magnetic braking of a magnetized differentially-rotating relativistic star have been studied by [114]. (Earlier works in the literature dealt with idealized, incompressible, uniformdensity stars.) In order to overcome the long evolution (Alfvén) timescale involved, of the order of 10 4 M for the models simulated, [114] adopted a perturbative metric approach and considered the limit of slow rotation and weak magnetic fields (still astrophysically realistic). This approach, however, permitted [114] to use a sufficiently high spatial resolution to track the development of the MRI and see its effects on the dynamics. The most important result found by [114] is that, independent of resolution, magnetic braking indeed drives the star toward uniform rotation within the timescales considered.

Neutron-star binary and black-hole-neutron-star coalescence
Many efforts in code development in numerical relativity and relativistic astrophysics aim to simulate the coalescence of compact binaries, composed of either two black holes, two neutron stars, or a black hole and a neutron star. The solution of the black-hole-binary problem has seen major breakthroughs in recent years, to the point that long-term stable evolutions for many orbits are now possible, and for various combinations of the parameter space of the problem, namely the spins and linear momentum of the black holes and their mass ratio (see [323] for a review). In this section we focus on simulations of neutron-star-binaries and neutron-star-black-hole configurations, for which progress has been equally dramatic.
Neutron-star binaries are among the most promising sources of gravitational radiation to be detected by the various ground-based interferometers worldwide. The computation of the gravitational waveform during the most dynamic phase of the coalescence and plunge depends crucially on hydrodynamic, finite-size effects. This phase begins once the stars, initially in quasi-equilibrium orbits of gradually smaller orbital radius (due to the emission of gravitational waves) reach the innermost stable circular orbit (ISCO). Approximately 10 8 years after formation of the binary system, the gravitational wave frequency enters the LIGO/VIRGO high frequency band. The final plunge of the two objects takes place on a dynamic timescale of a few ms. During the last 15 minutes before the stars finally merge, the frequency is inside the LIGO/VIRGO sensitivity range. About 16,000 cycles of waveform oscillation can be monitored, while the frequency gradually shifts from ∼ 10 Hz to ∼ 1 kHz. A perturbative treatment of gravitational radiation in the quadrupole approximation is valid as long as M/R 1 and M/r 1 simultaneously, M being the total mass of the binary, R the neutron star radius, and r the separation of the two stars. As the stars approach each other and merge, both inequalities are less valid and therefore, fully-relativistic hydrodynamic calculations become necessary.
Duez et al. [104] carried out the analytic modelling of the inspiral of corotational and irrotational relativistic neutron-star binaries as a sequence of quasi-equilibrium configurations, and computed the gravitational wavetrain from the last phase (a few hundred cycles) of the inspiral. These authors further showed a practical procedure to construct the entire wavetrain through coalescence by matching the late inspiral waveform to the one obtained by fully-relativistic hydrodynamic simulations. Detailed theoretical waveforms of the inspiral and plunge similar to those reported by [104] are crucial to enhance the chances of experimental detection in conjunction with matchedfiltering techniques.
The accurate simulation of a neutron-star-binary coalescence is, however, one of the most challenging tasks in numerical relativity. These scenarios involve strong gravitational fields, matter motion with (ultra) relativistic speeds, relativistic shock waves, and strong magnetic fields. The numerical difficulties are exacerbated by the intrinsic multidimensional character of the problem and by the inherent complexities in Einstein's theory of gravity, such as coordinate degrees of freedom and the possible formation of curvature singularities (e.g., collapse of matter configurations to black holes). It is thus not surprising that a large number of the available simulations have been attempted in the Newtonian (and post-Newtonian) framework (see [329] for a review). Many of these studies employ Lagrangian particle methods such as SPH, and only a few have considered (less viscous) high-order finite-difference methods such as PPM [342]. Notwithstanding the difficulties, major progress has been achieved recently in simulating neutron-star-binary mergers within a relativistic framework.
Concerning relativistic simulations, Wilson's formulation of the hydrodynamic equations (see Section 2.1.2) was used in [416,418,242]. Such investigations assumed a conformally-flat 3metric. These early simulations revealed the unexpected appearance of a "binary-induced collapse instability" of the neutron stars, prior to the eventual collapse of the final merged object. This effect was reduced, but not eliminated fully, in revised simulations [242], after Flanagan [124] pointed Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 out an error in the momentum-constraint equation as implemented by Wilson et al. [416,418]. A summary of this controversy can be found in [329]. Subsequent numerical simulations with the full set of Einstein equations (see below) did not find this effect.
The effect was neither found in the new set of CFC simulations performed by [300,299]. In the first investigation a new 3D SPH code was presented (see Section 4.4), tested, and applied to the coalescence of neutron-star binaries modelled by relativistic polytropes with Γ = 2.0, 2.6, and 3.0. Not only the system dynamics was studied but also gravitational-wave signals and luminosities were computed. The second investigation improved the microphysical treatment of the neutronstar models by using several nonzero temperature EOS, such as those of Shen et al. [351] and Lattimer-Swesty [214]. The results were compared with earlier investigations employing the BSSN equations performed by [375] (see below). It was found that the dynamics and outcome of the models depends strongly on the EOS and on the binary parameters (neutron-star masses and spins). Larger torus masses are found for asymmetric systems (up to 0.3 M for a neutron-star mass ratio of 0.55). Within tens of dynamic timescales only the Lattimer-Swesty EOS leads to a collapse of the postmerger remnant to a black hole. The gravitational wave emission from these simulations was analyzed in [298], where the authors concluded that the peak frequency of the post-merger signal is a sensitive indicator of the EOS, provided the total mass can be determined from the inspiral chirp signal. Nakamura et al. started a program in the late 1980's to simulate neutron-star-binary coalescence in general relativity (see, e.g., [280]). The group developed a three-dimensional code that solves the full set of Einstein and hydrodynamics equations through finite differences in a uniform Cartesian grid, using van Leer's scheme [403] with TVD flux limiters. Shockwaves are spread out using a tensor artificial-viscosity algorithm. The hydrodynamic equations follow Wilson's Eulerian formulation and the ADM formalism is adopted for the Einstein equations. This code has been tested with the study of the gravitational collapse of a rotating polytrope to a black hole (comparing to the axisymmetric computation of Stark and Piran [386]). A few results of neutron-star-binary simulations including gravitational radiation are reported in [304].
The headon collision of two neutron stars (a limiting case of a coalescence event) was considered by Miller et al. [261], who performed time-dependent relativistic simulations using the gr astro code. These simulations analyzed whether the collapse of the final object occurs in prompt timescales (a few milliseconds) or delayed (after neutrino cooling) timescales (a few seconds). In [349] it was argued that in a headon collision event, sufficient thermal pressure is generated to support the remnant in quasi-static equilibrium against (prompt) collapse, prior to slow cooling via neutrino emission (delayed collapse). In [261], prompt collapse to a black hole was found in the headon collision of two 1.4 M neutron stars modeled by a polytropic EOS with Γ = 2. The stars, initially separated by a proper distance of d = 44 km, were boosted toward one another at a speed of GM/d (the Newtonian infall velocity). The simulation employed a Cartesian grid of 192 3 points. The time evolution of this simulation can be followed in the Quicktime movie in Figure 17. This animation simultaneously shows the rest-mass density and the internal energy evolution during the on-axis collision. The formation of the black hole in prompt timescales is signalled by the sudden appearance of the apparent horizon at t = 0.16 ms (t = 63.194 in code units). The violet dotted circles indicate the trapped photons. The animation also shows a moderately-relativistic shockwave (Lorentz factor ∼ 1.2) appearing at t ∼ 36 (code units; ∼ 0.09 ms; yellow-white), which eventually is followed by two opposite moving shocks (along the infalling z direction) that propagate along the atmosphere surrounding the black hole. We note that, using the same code, critical phenomena has recently been found by [187] in the headon collison of two neutron stars.
Using the code of [109], Marronetti et al. [234] performed the first dynamic determination in relativistic gravity of the ISCO of a neutron-star binary (modelled as Γ = 2 polytropes). These authors were able to bracket the location of the ISCO by distinguishing stable circular orbits from Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7  [137,261]. The movie shows the evolution of the density and internal energy. The formation of the black hole in prompt timescales is demonstrated by the sudden appearance of the apparent horizon at t = 0.16 ms (t = 63.194 in code units), which is indicated by violet dotted circles representing the trapped photons. See [245] for download options of higher-quality versions of this movie. (To watch the movie, please go to the online version of this review article at http://www.livingreviews.org/lrr-2008-7.) unstable plunges. It was found that, for the simplified models considered, the ISCO frequency varies with compaction, but does not depend strongly on the stellar spin.
The most comprehensive study of neutron star coalescence in full general relativity has been performed by Shibata et al. [353,373,374,357,372,370]. The numerical code, briefly described in Section 4.4, allows the long-term simulation of the coalescences of both irrotational and corotational binaries, from the ISCO up to the formation and ringdown of the final collapsed object (either a black hole or a stable neutron star). The code also includes an apparent horizon finder along with microphysical EOS, and can extract the gravitational waveforms emitted in the collisions. The early simulations of Shibata and Uryu [374] accounted for a large sample of parameters of the binary system, such as the compactness of the (equal mass) neutron stars (0.12 ≤ M/R ≤ 0.16), the adiabatic index of the Γ-law EOS (1.8 ≤ Γ ≤ 2.5), and the maximum density, rest mass, gravitational mass, and total angular momentum. The initial data correspond to quasiequilibrium states, either corotational or irrotational, the latter being more realistic from considerations of viscous versus gravitational-radiation timescales. These initial data are obtained by solving the Einstein constraint equations and the equations for the gauge variables under the assumption of a conformallyflat 3-metric and the existence of a helicoidal Killing vector (see [374] for a detailed explanation). The binaries are chosen at the innermost orbits for which the Lagrange points appear at the inner edge of the neutron stars, and the plunge is induced by reducing the initial angular momentum by ∼ 2 -3%.
The comprehensive parameter-space survey carried out by [353,373,374] shows that the final outcome of the merger depends sensitively on the initial compactness of the neutron stars before plunge. Hence, depending on the stiffness of the EOS, controlled through the value of Γ, if the total rest mass of the system is ∼ 1.3 -1.7 times larger than the maximum rest mass of a spherical star in isolation, the end product is a black hole. Otherwise, a marginally-stable massive neutron star forms. In the latter case, the star is supported against self-gravity by rapid differential rotation.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 The image corresponds to the final snapshot of the evolution for the two models. The outcome of the merger for the more massive case is the delayed formation of a black hole, as signalled by the collapse of the lapse function shown in the lower panel of the image. In the less massive case, a long-lived hypermassive neutron star supported by rotation is formed. Simulations performed by [372]. (Used with permission.) The star could eventually collapse to a black hole once sufficient angular momentum has dissipated via neutrino emission or gravitational radiation. In turn, the different outcome of the merger is imprinted in the gravitational waveforms [374]. Therefore, future detection of high-frequency gravitational waves could help to constrain the maximum allowed mass of neutron stars. In addition, for prompt black-hole formation, a disk orbiting around the black hole forms, with a mass less than 1% the total rest mass of the system.
The input physics of those early simulations was later improved in [372] where hybrid EOS were adopted to mimic realistic nuclear equations of state. In this approach the EOS is divided into two parts, P = P cold + P th , where P cold is the cold part for which a fitting formula for realistic EOS of cold nuclear matter was assigned. The simulations used, in particular, the SLy and FPS EOS for which the maximum allowed ADM mass of cold and spherical neutron stars is ∼ 2.04 M and 1.80 M , respectively. Apart from improving the neutron-star models, the simulations of [372] also considered mergers of unequal-mass neutron stars. The underlying motivation was to find out whether massive accretion disks form depending on the initial mass ratio. Disk formation during neutron-star-binary coalescence, a fundamental issue for cosmological models of short duration GRBs, was indeed found to be enhanced for unequal-mass neutron stars, in which the less massive star is tidally disrupted during the evolution. The simulations also showed, in broad agreement with the simplified models, that depending on the threshold of the ADM mass of the system, which is also EOS dependent, the outcome is a black hole or a hypermassive neutron star with large ellipticity. For an ADM mass larger than the threshold, a black hole forms promptly in the merger irrespective of the mass ratio.
An example of such findings is shown in Figure 18. Both panels in this figure correspond to equal-mass neutron-star-binary mergers. The left panel is a 1.3 M − 1.3 M merger, while the right panel depicts the 1.35 M − 1.35 M case. Likewise, both panels display the velocity field and rest-mass isodensity contours in the equatorial plane at the final time of the evolution for the corresponding simulation. In addition, the lowest area of each panel shows the topology of the lapse function at the final time of the corresponding evolution. Despite the small difference of the initial neutron-star masses, the end product of each merger is remarkably different. In the former case, a hypermassive neutron star is formed, supported against collapse by centrifugal forces. In the latter case, the end product of this particular merger is the delayed formation of a rotating black hole. Animations of such simulations are available at [352].
Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 In the formation of a hypermassive neutron star with an ellipsoidal figure, quasiperiodic gravitational waves of large amplitude are emitted, with dimensionless strains of about h ∼ 6 − 7 × 10 −21 at a distance of 50 Mpc [357]. While these figures are within detectability for advanced laserinterferometric gravitational-wave detectors, the frequency of those waves, between 3 and 4 kHz, may, however, be too high. The eventual detection of such signals might help to further constrain the neutron star EOS.
A large number of neutron-star-binary merger simulations with an emphasis on the the blackhole formation case and on the resulting mass of the surrounding disk was presented by [370]. As in the simulations of [372], hybrid EOS were adopted to mimic realistic stiff nuclear EOS. In order to determine the mass of the disks that form, the simulations must be continued after black-hole formation, which [370] achieve with a simple excision technique. This technique, however, allows for long enough computations to even extract the ring-down gravitational waves associated with black-hole quasinormal modes. The resulting frequencies (too large) and amplitudes (too small at 50 Mpc) make unlikely the detection of this part of the signal by advanced laser interferometric detectors. One of the main results of the simulations of [370] was to find that the disk mass steeply increases with decreasing mass ratio for given ADM mass and EOS, suggesting that such small mass-ratio mergers are good candidates for producing the central engine of short-duration GRBs. Regarding the end product of the simulations, we note that the results of [370] do not agree with those of recent similar simulations of [297], at least for high-mass binaries. The reasons may be due to the different choice of nuclear EOS by both groups and also, perhaps more importantly, to the suppression of radiation-reaction effects on the gravitational waves (which rapidly dissipate angular momentum of the merged object) in the (approximate) CFC simulations of [297].
It is also worth noting the new neutron-star-binary merger simulations performed by [13,12], particularly the latter, for the innovative aspect of dealing with magnetized stars for the first time (which are, however, described by simple polytropes). The single simulation of a magnetized neutron star merger reported shows the feasibility of the numerical approach, leading to a strongly differentially-rotating object for the particular initial conditions used. Differences in the gravitational waveforms are found when comparing with a purely hydrodynamic merger. Those are motivated by the extremely large value of the magnetic field considered, which has a maximum magnitude of ∼ 10 16 G. In this respect we point out that Newtonian SPH simulations of magnetized neutron-star mergers carried out by [324] have revealed the amplification of the magnetic field beyond magnetar field strength. Such ultrastrong magnetic fields are caused by Kelvin-Helmholtz instabilities in the shear layer that forms between the neutron stars and lasts for a time scale of only 1 ms.
In recent years the first numerical simulations of the mixed binary system, i.e., that formed by a neutron star and a black hole, have also been presented [228,121,122,375,395,376,113]. Future simulations of such systems will benefit from the advances accomplished for the black-hole-binary case, and many are already adopting the same numerical treatment to achieve long-term stability in black-hole spacetimes, namely, the moving puncture approach (see [323] and references therein).
The simulations of [228] considered only the headon collision case, but thanks to the use of six levels of fixed mesh refinement, the evolutions were carried over up to times far beyond the collision, which provided time for the extraction of the gravitational wave signal and for estimating the radiated energy.
Fully three-dimensional simulations of the mixed binary system were first performed by [122]. However, two important simplifying assumptions were made in this study: first, the black hole was forced to remain fixed in space by making its mass sufficiently large compared to that of the neutron star, and, secondly, the treatment of the gravitational field equations was approximate due to the use of the conformal-flatness approximation. The first assumption (extreme massratio case) places the onset of tidal disruption of the neutron star outside the ISCO. The initial data are quasiequilibrium models for synchronized neutron-star polytropes generated as solutions Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 of the conformal thin-sandwich decomposition of the Einstein field equations (using the LORENE code [227]), from which relaxed SPH configurations were built and evolved using an SPH code. The main result of this work was to show that analytic models of mass transfer are not valid for neutron stars with low compactness, as they are highly unstable. The evolutions led to the formation of accretion disks around the black hole, which were discussed in [121] in the context of the mechanism that triggers short-hard GRBs.
The first fully self-consistent simulations of mixed compact binaries have been performed by Shibata and Uryu [375,376] and Shibata and Taniguchi [371]. These works use quasicircular states as initial conditions for the simulations, as the timescale of gravitational radiation reaction is a few times longer than the orbital period, and the black hole is modeled by a nonspinning moving puncture (see [323] and references therein for details). Furthermore, the first two works assume a corotating velocity field for the neutron star while the latter also considers the irrotational case. A Γ-law equation of state with Γ = 2 is used in all studies. High-order central schemes are used for the hydrodynamics and the Einstein equations are formulated and solved in the BSSN approach. The single simulation presented in [375] showed that a black-hole-plus-massive-disk (∼ 1 M ) system is not formed from nonspinning black-hole-neutron-star binaries, later confirmed in [376] for a slightly larger sample of models. A systematic study of black-hole-neutron-star-binary mergers was later reported in [371], focusing on the case that the neutron star is tidally disrupted. For all initial conditions, the neutron star was found to be tidally disrupted near the ISCO. The larger the size of the neutron star the more massive the resulting torus. The largest value found was ∼ 0.16 M for the case of a black hole of mass ∼ 4 M and a neutron star of mass ∼ 14.7 km.
The most recent simulations of mixed compact binaries are those of [113], and are also based on moving punctures. In this work, the assumptions adopted by the authors in their early work [122] are relaxed and the simulations are fully relativistic and self-consistent and do not only cover the extreme mass ratio case. They differ from those of [375,376,371] mainly in the way the initial data are built. Here, instead of quasicircular states, the conformal thin-sandwich formalism is used (see [395] for details). Other than this, the same BSSN formulation is used for the evolution of the gravitational field together with a conservative system for the hydrodynamics equations (that of [108]), which are solved using the HLL approximate-Riemann solver. The coupling between both sets of equations is based on the Cactus parallelization framework [244]. Despite different codes and initial data being used, the results of the simulations of [113] are reassuringly similar to those of [375,376,371]; the disks that form have masses, which may be too small to produce the required total energy output of a soft-hard GRB.

Additional Information
This last section of the review contains technical information concerning recent developments in the implementation of Riemann-solver-based numerical schemes in general-relativistic hydrodynamics.

Riemann problems in locally-Minkowskian coordinates
In [320], a procedure to integrate the general-relativistic hydrodynamics equations (as formulated in Section 2.1.3), taking advantage of the multitude of Riemann solvers developed in special relativity, was presented. The approach relies on a local change of coordinates in terms of which the spacetime metric is locally Minkowskian. This procedure allows, for 1D problems, the use of the exact solution of the special relativistic Riemann problem [239].
Such a coordinate transformation to locally Minkowskian coordinates at each numerical interface assumes that the solution of the Riemann problem is the one in special relativity and planar symmetry. This last assumption is equivalent to the approach followed in classical fluid dynamics, when using the solution of Riemann problems in slab symmetry for problems in cylindrical or spherical coordinates, as the solution breaks down near the singular points (e.g., the polar axis in cylindrical coordinates). In analogy to classical fluid dynamics, the numerical error depends on the magnitude of the Christoffel symbols, which might be large whenever huge gradients or large temporal variations of the gravitational field are present. Finer grids and improved time-advancing methods will be required in those circumstances.
Following [320], we illustrate the procedure for computing the second flux integral in Equation (76), which we call I. We begin by expressing the integral on a basis eα with e0 ≡ n µ and eî forming an orthonormal basis in the plane orthogonal to n µ with e1 normal to the surface Σ x 1 and e2 and e3 tangent to that surface. The vectors of this basis verify eα · eβ = ηαβ with ηαβ the Minkowski metric (in the following, caret superscripts will refer to vector components in this basis).
Denoting by x α 0 the coordinates at the center of the interface at time t, we introduce the following locally-Minkowskian coordinate system xα = Mα α (x α − x α 0 ), where the matrix Mα α is given by ∂ α = Mα α eα, calculated at x α 0 . In this system of coordinates the equations of generalrelativistic hydrodynamics transform into the equations of special relativistic hydrodynamics, in Cartesian coordinates, but with nonzero sources, and the flux integral reads (the caret symbol representing the numerical flux in Equation (76) is now removed to avoid confusion) with √ −ĝ = 1 + O(xα), where we have taken into account that, in the coordinates xα, Σ x 1 is described by the equation x1 − β1 α x0 = 0 (with βî = Mî i β i ), where the metric elements β 1 and α are calculated at x α 0 . Therefore, this surface is not at rest but moves with speed β1/α. At this point, all the theoretical work developed in recent years on special relativistic Riemann solvers can be exploited. The quantity in parentheses in Equation (103) represents the numerical flux across Σ x 1 , which can now be calculated by solving the special relativistic Riemann problem defined with the values at the two sides of Σ x 1 of two independent thermodynamic variables (namely, the rest mass density ρ and the specific internal energy ε) and the components of the velocity in the orthonormal spatial basis vî (vî = Mî i v i ).
Once the Riemann problem has been solved, we can take advantage of the self-similar character of the solution of the Riemann problem, which makes it constant on the surface Σ x 1 , simplifying the calculation of the above integral enormously: Living Reviews in Relativity http://www.livingreviews.org/lrr-2008-7 where the superscript (*) stands for the value on Σ x 1 obtained from the solution of the Riemann problem. Notice that the numerical fluxes correspond to the vector fields F 1 = {J, T·n, T·e1, T·e2, T · e3} and linearized Riemann solvers provide the numerical fluxes as defined in Equation (103).
Thus, the additional relation T · ∂ i = Mĵ i (T · eĵ) has to be used for the momentum equations. The integral in the right-hand side of Equation (104) is the area of the surface Σ x 1 and can be expressed in terms of the original coordinates as which can be evaluated for a given metric. The interested reader is directed to [320] for details on the testing and calibration of this procedure.

Acknowledgments
It is a pleasure to acknowledge José M. Ibáñez and José M. Martí for valuable suggestions and a careful reading of the manuscript. I am particularly grateful to Miguel Angel Aloy, Harald Dimmelmeier, Charles Gammie, Konstantinos Kifonidis, Shin Koide, Christian Ott, Luciano Rezzolla, José V. Romero, Masaru Shibata, and Wai-Mo Suen for providing some of the figures and animations contained in the article. This research has been supported in part by the Spanish Ministerio de Ciencia y Tecnología (grant AYA2004-08067-C03-01). The author has made use of SAO/NASA's Astrophysics Data System abstract service, which is gratefully acknowledged.