Abstract

The current status of numerical solutions for the equations of ideal general relativistic hydrodynamics is reviewed. With respect to an earlier version of the article, the present update providesadditional information on numerical schemes, and extends the discussion of astrophysical simulations ingeneral relativistic hydrodynamics. Different formulations of the equations are presented, with specialmention of conservative and hyperbolic formulations well-adapted to advanced numerical methods. Alarge sample of available numerical schemes is discussed, paying particular attention to solution proceduresbased on schemes exploiting the characteristic structure of the equations through linearized Riemannsolvers. A comprehensive summary of astrophysical simulations in strong gravitational fields is presented.These include gravitational collapse, accretion onto black holes, and hydrodynamical evolutions of neutronstars. The material contained in these sections highlights the numerical challenges of various representativesimulations. It also follows, to some extent, the chronological development of the field, concerning advanceson the formulation of the gravitational field and hydrodynamic equations and the numerical methodologydesigned to solve them. 1 Introduction

The description of important areas of modern astronomy, such as high-energy astrophysics or gravitational wave astronomy, requires general relativity. High-energy radiation is often emitted by highly relativistic events in regions of strong gravitational fields near compact objects such as neutron stars or black holes. The production of relativistic radio jets in active galactic nuclei, explained by pure hydrodynamical effects as in the twin-exhaust model [35] , by hydromagnetic centrifugal acceleration as in the Blandford–Payne mechanism [34] , or by electromagnetic extraction of energy as in the Blandford–Znajek mechanism [36] , involves an accretion disk around a rotating supermassive black hole. 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 (see, e.g., [287] ). A relativistic description is also necessary in scenarios involving explosive collapse of very massive stars (
$\sim 30{M}_{\odot}$
) to a black hole (in the so-called collapsar and hypernova models), or during the last phases of the coalescence of neutron star binaries. These catastrophic events are believed to exist at the central engine of highly energetic
$\gamma $
-ray bursts (GRBs) [214, 200, 306, 215] . In addition, non-spherical gravitational collapse leading to black hole formation or to a supernova explosion, and neutron star binary coalescence are among the most promising sources of detectable gravitational radiation. Such astrophysical scenarios constitute one of the main targets for the new generation of ground-based laser interferometers, just starting their gravitational wave search (LIGO, VIRGO, GEO600, TAMA) [285, 205] .

A powerful way to improve our understanding of the above scenarios is through accurate, large scale, three-dimensional numerical simulations. Nowadays, computational general relativistic astrophysics is an increasingly important field of research. In addition to the large amount of observational data by high-energy Xand
$\gamma $
-ray satellites such as Chandra, XMM-Newton, or INTEGRAL, and the new generation of gravitational wave detectors, 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, the computational astrophysicist and the numerical relativist face a daunting task. In the most 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 (magneto-)hydrodynamic (MHD) equations and the Einstein gravitational field equations.

In many cases, the number of equations must be augmented to account for non-adiabatic processes, e.g., radiative transfer or sophisticated microphysics (realistic equations of state for nuclear matter, nuclear physics, magnetic fields, and so on).

Nevertheless, in some astrophysical situations of interest, e.g., accretion of matter onto compact objects or oscillations of relativistic stars, the “test fluid” approximation is enough to get an accurate description of the underlying dynamics. In this approximation the fluid self-gravity 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 (i.e., with the stress-energy tensor being that of a perfect fluid), is also a fairly standard choice in numerical astrophysics.

The main purpose of this review is to summarize the existing efforts to solve numerically the equations of (ideal) general relativistic hydrodynamics. To this aim, the most important numerical schemes will be presented first in some detail. Prominence will be given to the so-called Godunov-type schemes written in conservative form. Since [162] , it has been demonstrated gradually [93, 78, 243, 83, 21, 296, 228] that conservative methods exploiting the hyperbolic character of the relativistic hydrodynamic equations are optimally suited for accurate numerical integrations, even well inside the ultrarelativistic regime. The explicit knowledge of the characteristic speeds (eigenvalues) of the equations, together with the corresponding eigenvectors, provides the mathematical (and physical) framework for such integrations, by means of either exact or approximate Riemann solvers. The article includes, furthermore, a comprehensive description of “relevant” numerical applications in relativistic astrophysics, including gravitational collapse, accretion onto compact objects, and hydrodynamical evolution of neutron stars. Numerical simulations of strong-field scenarios employing Newtonian gravity and hydrodynamics, as well as possible post-Newtonian extensions, have received considerable attention in the literature and will not be covered in this review, which focuses on relativistic simulations. Nevertheless, we must emphasize that most of what is known about hydrodynamics near compact objects, in particular in black hole astrophysics, has been accurately described using Newtonian models. Probably the best known example is the use of a pseudo-Newtonian potential for non-rotating black holes that mimics the existence of an event horizon at the Schwarzschild gravitational radius [216] . This has allowed accurate interpretations of observational phenomena.

The organization of this article is as follows: Section 2 presents the equations of general relativistic hydrodynamics, summarizing the most relevant theoretical formulations that, to some extent, have helped to drive the development of numerical algorithms for their solution. Section 3 is mainly devoted to describing numerical schemes specifically designed to solve nonlinear hyperbolic systems of conservation laws. Hence, particular emphasis will be paid on conservative high-resolution shock-capturing (HRSC) upwind methods based on linearized Riemann solvers. Alternative schemes such as smoothed particle hydrodynamics (SPH), (pseudo-)spectral methods, and others will be briefly discussed as well. Section 4 summarizes a comprehensive sample of hydrodynamical simulations in strong-field general relativistic astrophysics. Finally, in Section 5 we provide additional technical information needed to build up upwind HRSC schemes for the general relativistic hydrodynamics equations. Geometrized units (
$G=c=1$
) are used throughout the paper except where explicitly indicated, as well as the metric conventions of [185] . Greek (Latin) indices run from 0 to 3 (1 to 3).

2 Equations of General Relativistic Hydrodynamics

The general relativistic hydrodynamic equations consist of the local conservation laws of the stress-energy tensor
${T}^{\mu \nu}$
(the Bianchi identities) and of the matter current density
${J}^{\mu}$
(the continuity equation):

As usual,
${\nabla}_{\mu}$
stands for the covariant derivative associated with the four-dimensional spacetime metric
${g}_{\mu \nu}$
. The density current is given by
${J}^{\mu}=\rho {u}^{\mu}$
,
${u}^{\mu}$
representing the fluid 4-velocity and
$\rho $
the rest-mass density in a locally inertial reference frame.

$$\begin{array}{c}{\nabla}_{\mu}{T}^{\mu \nu}=0,\end{array}$$ | (1) |

$$\begin{array}{c}{\nabla}_{\mu}{J}^{\mu}=0.\end{array}$$ | (2) |

The stress-energy tensor for a non-perfect fluid is defined as

where
$\varepsilon $
is the specific energy density of the fluid in its rest frame,
$p$
is the pressure, and
${h}^{\mu \nu}$
is the spatial projection tensor
${h}^{\mu \nu}={u}^{\mu}{u}^{\nu}+{g}^{\mu \nu}$
. In addition,
$\eta $
and
$\zeta $
are the shear and bulk viscosities.

$$\begin{array}{c}{T}^{\mu \nu}=\rho (1+\varepsilon ){u}^{\mu}{u}^{\nu}+(p-\zeta \theta ){h}^{\mu \nu}-2\eta {\sigma}^{\mu \nu}+{q}^{\mu}{u}^{\nu}+{q}^{\nu}{u}^{\mu},\end{array}$$ | (3) |

The expansion
$\theta $
, describing the divergence or convergence of the fluid world lines, is defined as
$\theta ={\nabla}_{\mu}{u}^{\mu}$
. The symmetric, trace-free, spatial shear tensor
${\sigma}^{\mu \nu}$
is defined by

and, finally,
${q}^{\mu}$
is the energy flux vector.

$$\begin{array}{c}{\sigma}^{\mu \nu}=\frac{1}{2}({\nabla}_{\alpha}{u}^{\mu}{h}^{\alpha \nu}+{\nabla}_{\alpha}{u}^{\nu}{h}^{\alpha \mu})-\frac{1}{3}\theta {h}^{\mu \nu},\end{array}$$ | (4) |

In the following we will neglect non-adiabatic effects, such as viscosity or heat transfer, assuming the stress-energy tensor to be that of a perfect fluid

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,
$\sqrt{-g}$
is the volume element associated with the 4-metric, with
$g=det\left({g}_{\mu \nu}\right)$
, and the
${\Gamma}_{\mu \lambda}^{\nu}$
are the 4-dimensional Christoffel symbols.

$$\begin{array}{c}{T}^{\mu \nu}=\rho h{u}^{\mu}{u}^{\nu}+p{g}^{\mu \nu},\end{array}$$ | (5) |

$$\begin{array}{c}h=1+\varepsilon +\frac{p}{\rho}.\end{array}$$ | (6) |

$$\begin{array}{c}\frac{\partial}{\partial {x}^{\mu}}\sqrt{-g}{J}^{\mu}=0,\end{array}$$ | (7) |

$$\begin{array}{c}\frac{\partial}{\partial {x}^{\mu}}\sqrt{-g}{T}^{\mu \nu}=-\sqrt{-g}{\Gamma}_{\mu \lambda}^{\nu}{T}^{\mu \lambda},\end{array}$$ | (8) |

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 thermodynamical quantities. In general, the EOS takes the form

Due to their simplicity, the most widely employed EOSs in numerical simulations are the ideal fluid EOS,
$p=(\Gamma -1)\rho \varepsilon $
, where
$\Gamma $
is the adiabatic index, and the polytropic EOS (e.g., to build equilibrium stellar models),
$p=K{\rho}^{\Gamma}\equiv K{\rho}^{1+1/N}$
,
$K$
being the polytropic constant and
$N$
the polytropic index.

$$\begin{array}{c}p=p(\rho ,\varepsilon ).\end{array}$$ | (9) |

In the “test fluid” approximation, where the fluid self-gravity is neglected, the dynamics of the system are completely governed by Equations ( 1 ) and ( 2 ), together with the EOS ( 9 ). In those situations where 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 dynamical spacetime. A detailed description of the various numerical approaches to solve the Einstein equations is beyond the scope of the present article (see, e.g., Lehner [?] for a recent review). We briefly mention that the Einstein equations, in the presence of matter fields, can be formulated as an initial value (Cauchy) problem, using the so-called 3+1 decomposition of spacetime [15] . More details can be found in, e.g., [314] . Given a choice of gauge, the Einstein equations in the 3+1 formalism [15] split into evolution equations for the 3-metric
${\gamma}_{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. Long-term stable evolutions of the Einstein equations have recently been accomplished using various reformulations of the original 3+1 system (see, e.g., [25, 257, 6, 89] for simulations involving matter sources, and [5] and references therein for vacuum black-hole evolutions). Alternatively, a characteristic initial value problem formulation of the Einstein equations was developed in the 1960s by Bondi, van der Burg, and Metzner [45] , and Sachs [246] . This approach has gradually advanced to a state where long-term stable evolutions of caustic-free spacetimes in multi-dimensions are possible, even including matter fields (see [?] and references therein). A recent review of the characteristic formulation is presented in a Living Reviews article by Winicour [304] . Examples of this formulation in general relativistic hydrodynamics are discussed in various sections of the present article.

$$\begin{array}{c}{G}^{\mu \nu}=8\pi {T}^{\mu \nu},\end{array}$$ | (10) |

Traditionally, most of the approaches for numerical integrations of the general relativistic hydrodynamic equations have adopted spacelike foliations of the spacetime, within the 3+1 formulation.

Recently, however, covariant forms of these equations, well suited for advanced numerical methods, have also been developed. This is reviewed next in a chronological way.

2.1 Spacelike 3+1 approaches

In the 3+1 (ADM) formulation [15] , spacetime is foliated into a set of non-intersecting spacelike hypersurfaces. There are two kinematic variables describing the evolution between these surfaces:

the lapse function
$\alpha $
, which describes the rate of advance of time along a timelike unit vector
${n}^{\mu}$
normal to a surface, and the spacelike shift vector
${\beta}^{i}$
that describes the motion of coordinates within a surface.

The line element is written as

where
${\gamma}_{ij}$
is the 3-metric induced on each spacelike slice.

$$\begin{array}{c}d{s}^{2}=-({\alpha}^{2}-{\beta}_{i}{\beta}^{i})d{x}^{0}d{x}^{0}+2{\beta}_{i}d{x}^{i}d{x}^{0}+{\gamma}_{ij}d{x}^{i}d{x}^{j},\end{array}$$ | (11) |

2.1.1 $\mathbf{}$ 1+1 Lagrangian formulation (May and White)

The pioneering numerical work in general relativistic hydrodynamics dates back to the one-dimensional gravitational collapse code of May and White [171, 172] . Building on theoretical work by Misner and Sharp [184] , 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 3.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 worth describing its main features in some detail.

For a spherically symmetric spacetime, the line element can be written as

$m$
being a radial (Lagrangian) coordinate, indicating the total rest-mass enclosed inside the circumference
$2\pi R(m,t)$
.

$$\begin{array}{c}d{s}^{2}=-{a}^{2}(m,t)d{t}^{2}+{b}^{2}(m,t)d{m}^{2}+{R}^{2}(m,t)(d{\theta}^{2}+{sin}^{2}\theta d{\phi}^{2}),\end{array}$$ | (12) |

The co-moving character of the coordinates leads, for a perfect fluid, to a stress-energy tensor of the form

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 [184] ):

with the definitions
$u=1/a\cdot \partial R/\partial t$
and
$\Gamma =1/b\cdot \partial R/\partial m$
, satisfying
${\Gamma}^{2}=1-{u}^{2}-2M/R$
.

$$\begin{array}{ccc}{T}_{1}^{1}={T}_{2}^{2}={T}_{3}^{3}& =& -p,\end{array}$$ | (13) |

$$\begin{array}{ccc}{T}_{0}^{0}& =& (1+\varepsilon )\rho ,\end{array}$$ | (14) |

$$\begin{array}{ccc}{T}_{\nu}^{\mu}& =& 0\text{if}\mu \ne \nu .\end{array}$$ | (15) |

$$\begin{array}{c}b=\frac{1}{4\pi \rho {R}^{2}}.\end{array}$$ | (16) |

$$\begin{array}{c}\begin{array}{ccc}\frac{\partial p}{\partial m}+\frac{1}{a}\frac{\partial a}{\partial m}\rho h& =& 0,\\ \frac{\partial \varepsilon}{\partial t}+p\frac{\partial}{\partial t}\left(\frac{1}{\rho}\right)& =& 0,\\ \frac{\partial u}{\partial t}& =& -a\left(4\pi \frac{\partial p}{\partial m}{R}^{2}\frac{\Gamma}{h}+\frac{M}{{R}^{2}}+4\pi pR\right),\\ \frac{1}{\rho {R}^{2}}\frac{\partial \rho {R}^{2}}{\partial t}& =& -a\frac{\partial u/\partial m}{\partial R/\partial m},\end{array}\end{array}$$ | (17) |

Additionally,

represents the total mass interior to radius
$m$
at time
$t$
. The final system, Equations ( 17 ), is closed with an EOS of the form given by Equation ( 9 ).

$$\begin{array}{c}M={\int}_{0}^{m}4\pi {R}^{2}\rho (1+\varepsilon )\frac{\partial R}{\partial m}dm,\end{array}$$ | (18) |

Hydrodynamics codes based on the original formulation of May and White and on later versions (e.g., [292] ) have been used in many nonlinear simulations of supernova and neutron star collapse (see, e.g., [183, 279] and references therein), as well as in perturbative computations of spherically symmetric gravitational collapse within the framework of the linearized Einstein equations [250, 251] . In Section 4.1.1 below, some of these simulations are discussed in detail. An interesting analysis of the above formulation in the context of gravitational collapse is provided by Miller and Sciama [181] . By comparing the Newtonian and relativistic equations, these authors showed that the net acceleration of the infalling mass shells is larger in general relativity than in Newtonian gravity. The Lagrangian character of May and White's formulation, together with other theoretical considerations concerning the particular coordinate gauge, has prevented its extension to multi-dimensional calculations. However, 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 lack of numerical diffusion.

2.1.2 $\mathbf{}$ 3+1 Eulerian formulation (Wilson)

The use of Eulerian coordinates in multi-dimensional numerical relativistic hydrodynamics started with the pioneering work by Wilson [298] . Introducing the basic dynamical variables
$D$
,
${S}_{\mu}$
, and
$E$
, representing the relativistic density, momenta, and energy, respectively, defined as

the equations of motion in Wilson's formulation [298, 299] are

with the “transport velocity” given by
${V}^{\mu}={u}^{\mu}/{u}^{0}$
. We note that in the original formulation [299] the momentum density equation, Equation ( 21 ), is only solved for the three spatial components
${S}_{i}$
, and
${S}_{0}$
is obtained through the 4-velocity normalization condition
${u}_{\mu}{u}^{\mu}=-1$
.

$$\begin{array}{c}D=\rho {u}^{0},{S}_{\mu}=\rho h{u}_{\mu}{u}^{0},E=\rho \varepsilon {u}^{0},\end{array}$$ | (19) |

$$\begin{array}{c}\frac{1}{\sqrt{-g}}\frac{\partial}{\partial {x}^{0}}\left(D\sqrt{-g}\right)+\frac{1}{\sqrt{-g}}\frac{\partial}{\partial {x}^{i}}\left(D{V}^{i}\sqrt{-g}\right)=0,\end{array}$$ | (20) |

$$\begin{array}{c}\frac{1}{\sqrt{-g}}\frac{\partial}{\partial {x}^{0}}\left({S}_{\mu}\sqrt{-g}\right)+\frac{1}{\sqrt{-g}}\frac{\partial}{\partial {x}^{i}}\left({S}_{\mu}{V}^{i}\sqrt{-g}\right)+\frac{\partial p}{\partial {x}^{\mu}}+\frac{1}{2}\frac{\partial {g}^{\alpha \beta}}{\partial {x}^{\mu}}\frac{{S}_{\alpha}{S}_{\beta}}{{S}^{0}}=0,\end{array}$$ | (21) |

$$\begin{array}{c}\frac{\partial}{\partial {x}^{0}}\left(E\sqrt{-g}\right)+\frac{\partial}{\partial {x}^{i}}\left(E{V}^{i}\sqrt{-g}\right)+p\frac{\partial}{\partial {x}^{\mu}}\left({u}^{0}{V}^{\mu}\sqrt{-g}\right)=0,\end{array}$$ | (22) |

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 [298] 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 [294] to the relativistic regime (see Section 3.1.1 ).

Wilson's formulation has been widely used in hydrodynamical codes developed by a variety of research groups. Many different astrophysical scenarios were first investigated with these codes, including axisymmetric stellar core collapse [194, 192, 198, 22, 275, 227, 79] , accretion onto compact objects [122, 225] , numerical cosmology [53, 54, 11] and, more recently, the coalescence of neutron star binaries [302, 303, 168] . This formalism has also been employed, in the special relativistic limit, in numerical studies of heavy-ion collisions [301, 174] . We note that in most of these investigations, the original formulation of the hydrodynamic equations was slightly modified by re-defining the dynamical variables, Equation ( 19 ), with the addition of a multiplicative
$\alpha $
factor (the lapse function) and the introduction of the Lorentz factor,
$W\equiv \alpha {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 [299] , building on a vacuum numerical relativity code specifically developed to investigate the head-on collision of two black holes [272] . The resulting code was axially symmetric and aimed to integrate the coupled set of equations in the context of stellar core collapse [82] .

$$\begin{array}{c}D=\rho W,{S}_{\mu}=\rho hW{u}_{\mu},E=\rho \varepsilon W.\end{array}$$ | (23) |

More recently, Wilson's formulation has been applied to the numerical study of the coalescence of binary neutron stars in general relativity [302, 303, 168] (see Section 4.3.2 ). These studies adopted an approximation scheme for the gravitational field, by imposing the simplifying condition that the 3-geometry (the 3-metric
${\gamma}_{ij}$
) is conformally flat. The line element, Equation ( 11 ), then reads

The curvature of the 3-metric is then described by a position dependent conformal factor
${\phi}^{4}$
times a flat-space Kronecker delta. Therefore, in this approximation scheme all radiation degrees of freedom are removed, while 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 [139] 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.
Wilson's formulation showed some limitations in handling situations involving ultrarelativistic flows (
$W\gg 2$
), as first pointed out by Centrella and Wilson [54] . Norman and Winkler [207] performed a comprehensive numerical assessment of such formulation by means of special relativistic hydrodynamical simulations. Figure 1 reproduces a plot from [207] in which the relative error of the density compression ratio in the so-called 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 [54] . This figure shows that for Lorentz factors of about 2 (
$v\approx 0.86c$
), which is 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$
.

$$\begin{array}{c}d{s}^{2}=-({\alpha}^{2}-{\beta}_{i}{\beta}^{i})d{x}^{0}d{x}^{0}+2{\beta}_{i}d{x}^{i}d{x}^{0}+{\phi}^{4}{\delta}_{ij}d{x}^{i}d{x}^{j}.\end{array}$$ | (24) |

Norman and Winkler [207] concluded 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 Section 3.1.1 ), are only added to the pressure terms 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, [207] proposed to add the
$Q$
terms in a relativistically 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:

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
$\rho h$
. In general,
$Q$
is a nonlinear function of the velocity and, hence, the quantity
$Q{W}^{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 much more 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.

$$\begin{array}{c}{T}^{\mu \nu}=\rho (1+\varepsilon +(p+Q)/\rho ){u}^{\mu}{u}^{\nu}+(p+Q){g}^{\mu \nu}.\end{array}$$ | (25) |

$$\begin{array}{c}\frac{\partial}{\partial {x}^{0}}\left[\right(\rho h+Q\left){W}^{2}{V}_{j}\right]+\frac{\partial}{\partial {x}^{i}}\left[\right(\rho h+Q\left){W}^{2}{V}_{j}{V}^{i}\right]+\frac{\partial (p+Q)}{\partial {x}^{j}}=0.\end{array}$$ | (26) |

Their code, which in addition incorporates an adaptive grid, reproduces very accurate results even for ultrarelativistic flows with Lorentz factors of about 10 in one-dimensional, flat spacetime simulations.

Very recently, Anninos and Fragile [13] have compared state-of-the-art artificial viscosity schemes and high-order non-oscillatory central schemes (see Section 3.1.3 ) using Wilson's formulation for the former class of schemes and a conservative formulation (similar to the one considered in [220, 218] ; Section 2.2.2 ) for the latter. Using a three-dimensional Cartesian code, 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 about
$\sim 0.95c$
. On the other hand, results for the shock reflection problem with a second-order finite difference central scheme show the suitability of such a scheme to handle ultrarelativistic flows, the underlying reason being, most likely, the use of a conservative formulation of the hydrodynamic equations rather than the particular scheme employed (see Section 3.1.3 ).

Tests concerning spherical accretion onto a Schwarzschild black hole using both schemes yield the maximum relative errors near the event horizon, as large as
$\sim 24$
% for the central scheme.

2.1.3 $\mathbf{}$ 3+1 conservative Eulerian formulation (Ibán͂ez and coworkers)

In 1991, Martí, Ibán͂ez, and Miralles [162] presented a new formulation of the (Eulerian) general relativistic hydrodynamic equations. This formulation was aimed to take fundamental advantage of the hyperbolic and conservative character of the equations, contrary to the one discussed in the previous section. 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 the existing tools of classical fluid dynamics. This procedure departs from earlier approaches, most notably in avoiding the need for either artificial dissipation terms to handle discontinuous solutions [298, 299] , or implicit schemes as proposed in [207] .

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., [151] ). 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 [162] .

Theoretical advances on 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 magneto-hydrodynamics (MHD) was exhaustively studied by Anile and collaborators (see [10] and references therein) by applying Friedrichs' definition of hyperbolicity [100] to a quasi-linear form of the system of hydrodynamic equations,

where
${\mathcal{A}}^{\mu}$
are the Jacobian matrices of the system and
$\mathbf{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
$\xi $
with
${\xi}_{\mu}{\xi}^{\mu}=-1$
, if the following two conditions hold: (i)
$det\left({\mathcal{A}}^{\mu}{\xi}_{\mu}\right)\ne 0$
and (ii) for any
$\zeta $
such that
${\zeta}_{\mu}{\xi}^{\mu}=0$
,
${\zeta}_{\mu}{\zeta}^{\mu}=1$
, the eigenvalue problem
${\mathcal{A}}^{\mu}({\zeta}_{\mu}-\lambda {\xi}_{\mu})\mathbf{r}=0$
has only real eigenvalues
$\{{\lambda}_{n}{\}}_{n=1,\cdots ,5}$
and a complete set of right-eigenvectors
$\{{\mathbf{r}}_{n}{\}}_{n=1,\cdots ,5}$
. Besides verifying the hyperbolic character of the relativistic hydrodynamic equations, Anile and collaborators [10] obtained the explicit expressions for the eigenvalues and eigenvectors in the local rest frame, characterized by
${u}^{\mu}={\delta}_{0}^{\mu}$
. In Font et al. [93] those calculations were extended to an arbitrary reference frame in which the motion of the fluid is described by the 4-velocity
${u}^{\mu}=W(1,{v}^{i})$
.

$$\begin{array}{c}{\mathcal{A}}^{\mu}(\mathbf{w})\frac{\partial \mathbf{w}}{\partial {x}^{\mu}}=0,\end{array}$$ | (27) |

The approach followed in [93] for the equations of special relativistic hydrodynamics was extended to general relativity in [21] . The choice of evolved variables (conserved quantities) in the 3+1 Eulerian formulation developed by Banyuls et al. [21] differs slightly from that of Wilson's formulation [298] .

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 are the natural extension (for a generic spacetime) of the Eulerian observers in classical fluid dynamics. Interested readers are directed to [21] for more complete definitions and geometrical foundations.

In terms of the so-called primitive variables
$\mathbf{w}=(\rho ,{v}_{i},\varepsilon )$
, the conserved quantities are written as

where the contravariant components
${v}^{i}={\gamma}^{ij}{v}_{j}$
of the 3-velocity are defined as

and
$W$
is the relativistic Lorentz factor
$W\equiv \alpha {u}^{0}=(1-{v}^{2}{)}^{-1/2}$
with
${v}^{2}={\gamma}_{ij}{v}^{i}{v}^{j}$
.

$$\begin{array}{c}\begin{array}{ccc}D& =& \rho W,\\ {S}_{j}& =& \rho h{W}^{2}{v}_{j},\\ E& =& \rho h{W}^{2}-p,\end{array}\end{array}$$ | (28) |

$$\begin{array}{c}{v}^{i}=\frac{{u}^{i}}{\alpha {u}^{0}}+\frac{{\beta}^{i}}{\alpha},\end{array}$$ | (29) |

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-energy tensor components. More precisely, the first-order flux-conservative hyperbolic system, well suited for numerical applications, reads

with
$g\equiv det\left({g}_{\mu \nu}\right)$
satisfying
$\sqrt{-g}=\alpha \sqrt{\gamma}$
with
$\gamma \equiv det\left({\gamma}_{ij}\right)$
. The state vector is given by

with
$\tau \equiv E-D$
. The vector of fluxes is

and the corresponding sources
$\mathbf{S}(\mathbf{w})$
are

The local characteristic structure of the previous system of equations was presented in [21] .

$$\begin{array}{c}\frac{1}{\sqrt{-g}}\left(\frac{\partial \sqrt{\gamma}\mathbf{U}(\mathbf{w})}{\partial {x}^{0}}+\frac{\partial \sqrt{-g}{\mathbf{F}}^{i}(\mathbf{w})}{\partial {x}^{i}}\right)=\mathbf{S}(\mathbf{w}),\end{array}$$ | (30) |

$$\begin{array}{c}\mathbf{U}(\mathbf{w})=(D,{S}_{j},\tau ),\end{array}$$ | (31) |

$$\begin{array}{c}{\mathbf{F}}^{i}(\mathbf{w})=\left(D\left({v}^{i}-\frac{{\beta}^{i}}{\alpha}\right),{S}_{j}\left({v}^{i}-\frac{{\beta}^{i}}{\alpha}\right)+p{\delta}_{j}^{i},\tau \left({v}^{i}-\frac{{\beta}^{i}}{\alpha}\right)+p{v}^{i}\right),\end{array}$$ | (32) |

$$\begin{array}{c}\mathbf{S}(\mathbf{w})=\left(0,{T}^{\mu \nu}\left(\frac{\partial {g}_{\nu j}}{\partial {x}^{\mu}}-{\Gamma}_{\nu \mu}^{\delta}{g}_{\delta j}\right),\alpha \left({T}^{\mu 0}\frac{\partial ln\alpha}{\partial {x}^{\mu}}-{T}^{\mu \nu}{\Gamma}_{\nu \mu}^{0}\right)\right).\end{array}$$ | (33) |

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. System ( 30 ) satisfies, hence, the definition of hyperbolicity. As it will become apparent in Section 3.1.2 below, the knowledge of the spectral information is essential in order to construct HRSC schemes based on Riemann solvers. This information can be found in [21] (see also [96] ).

The range of applications considered so far in general relativity employing the above formulation of the hydrodynamic equations, Equation ( 30 , 31 , 32 , 33 ), is still small and mostly devoted to the study of stellar core collapse and accretion flows onto black holes (see Sections 4.1.1 and 4.2 below). In the special relativistic limit this formulation is being successfully applied to simulate the evolution of (ultra-)relativistic extragalactic jets, using numerical models of increasing complexity (see, e.g., [166, 8] ). The first applications in general relativity were performed, in one spatial dimension, in [162] , using a slightly different form of the equations. Preliminary investigations of gravitational stellar collapse were attempted by coupling the above formulation of the hydrodynamic equations to a hyperbolic formulation of the Einstein equations developed by [39] . These results are discussed in [160, 38] . More recently, successful evolutions of fully dynamical spacetimes in the context of adiabatic stellar core collapse, both in spherical symmetry and in axisymmetry, have been achieved [129, 243, 67] . These investigations are considered in Section 4.1.1 below.

An ambitious three-dimensional, Eulerian code which evolves the coupled system of Einstein and hydrodynamics equations was developed by Font et al. [96] (see Section 3.3.2 ). The formulation of the hydrodynamic equations in this code follows the conservative Eulerian approach discussed in this section. The code is constructed for a completely general spacetime metric based on a Cartesian coordinate system, with arbitrarily specifiable lapse and shift conditions. In [96] the spectral decomposition (eigenvalues and right-eigenvectors) of the general relativistic hydrodynamic equations, valid for general spatial metrics, was derived, extending earlier results of [21] for non-diagonal metrics. A complete set of left-eigenvectors was presented by Ibán͂ez et al. [127] . 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 5.2 of this article.

2.2 Covariant approaches

General (covariant) conservative formulations of the general relativistic hydrodynamic equations for ideal fluids, i.e., not restricted to spacelike foliations, have been reported in [78] and, more recently, in [220, 218] . 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 we describe the existing covariant formulations in some detail.

2.2.1 Eulderink and Mellema

Eulderink and Mellema [76, 78] first derived a covariant formulation of the general relativistic hydrodynamic 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 on a generalization of Roe's approximate Riemann solver.

Such a solver was first applied to the non-relativistic Euler equations in [241] 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=(\Gamma -1)\rho \varepsilon $
,
$\Gamma $
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 re-written in flux-conservative form. The flow variables are then expressed in terms of a parameter vector
$\omega $
as

where
${\omega}^{\alpha}\equiv K{u}^{\alpha}$
,
${\omega}^{4}\equiv Kp/\left(\rho h\right)$
and
${K}^{2}\equiv \sqrt{-g}\rho h=-{g}_{\alpha \beta}{\omega}^{\alpha}{\omega}^{\beta}$
. The vector
${\mathbf{F}}^{0}$
represents the state vector (the unknowns), and each vector
${\mathbf{F}}^{i}$
is the corresponding flux in the coordinate direction
${x}^{i}$
.

$$\begin{array}{c}{\mathbf{F}}^{\alpha}=\left(\left[K-\frac{\Gamma}{\Gamma -1}{\omega}^{4}\right]{\omega}^{\alpha},{\omega}^{\alpha}{\omega}^{\beta}+K{\omega}^{4}{g}^{\alpha \beta}\right),\end{array}$$ | (34) |

Eulderink and Mellema computed the exact “Roe matrix” [241] for the vector ( 34 ) 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 state
$\stackrel{~}{\omega}=({\omega}_{L}+{\omega}_{R})/({K}_{L}+{K}_{R})$
, where L and R denote the left and right states in a Riemann problem (see Section 3.1.2 ). Further technical details can be found in [76, 78] .

The performance of this general relativistic Roe solver was tested in a number of one-dimensional problems for which exact solutions are known, including non-relativistic shock tubes, special relativistic shock tubes, and 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 [77] . However, no astrophysical applications in strong-field general relativistic flows have yet been attempted with this formulation.

2.2.2 Papadopoulos and Font

In this formulation [220] , the spatial velocity components of the 4-velocity,
${u}^{i}$
, together with the rest-frame density and internal energy,
$\rho $
and
$\varepsilon $
, 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 a five dimensional space
$\mathbf{w}=(\rho ,{u}^{i},\varepsilon )$
. 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,
$\mathbf{U}$
, individually denoted
$(D,{S}^{i},E)$
:

Note that the state vector variables slightly differ from previous choices (see, e.g., Equations ( 19 ), ( 28 ), and ( 34 )). With those definitions the equations of general relativistic hydrodynamics take the standard conservation law form,

with
$A=(0,i,4)$
. The flux vectors
${\mathbf{F}}^{j}$
and the source terms
$\mathbf{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
$\mathbf{U}$
or
$\mathbf{w}$
, and each one can be obtained from the other via the definitions ( 35 ) and the use of the normalization condition for the 4-velocity,
${g}_{\mu \nu}{u}^{\mu}{u}^{\nu}=-1$
. The local characteristic structure of the above system of equations was presented in [220] , where the formulation proved well suited for the numerical implementation of HRSC schemes. The formulation presented in this section has been developed for a perfect fluid EOS. Extensions to account for generic EOS are given in [218] .

$$\begin{array}{c}\begin{array}{ccc}D& =& {\mathbf{U}}^{0}={J}^{0}=\rho {u}^{0},\\ {S}^{i}& =& {\mathbf{U}}^{i}={T}^{0i}=\rho h{u}^{0}{u}^{i}+p{g}^{0i},\\ E& =& {\mathbf{U}}^{4}={T}^{00}=\rho h{u}^{0}{u}^{0}+p{g}^{00}.\end{array}\end{array}$$ | (35) |

$$\begin{array}{c}\frac{\partial \left(\sqrt{-g}{\mathbf{U}}^{A}\right)}{\partial {x}^{0}}+\frac{\partial \left(\sqrt{-g}{\mathbf{F}}^{j}\right)}{\partial {x}^{j}}=\mathbf{S},\end{array}$$ | (36) |

$$\begin{array}{c}{\mathbf{F}}^{j}=({J}^{j},{T}^{ji},{T}^{j0})=(\rho {u}^{j},\rho h{u}^{i}{u}^{j}+p{g}^{ij},\rho h{u}^{0}{u}^{j}+p{g}^{0j}),\end{array}$$ | (37) |

$$\begin{array}{c}\mathbf{S}=(0,-\sqrt{-g}{\Gamma}_{\mu \lambda}^{i}{T}^{\mu \lambda},-\sqrt{-g}{\Gamma}_{\mu \lambda}^{0}{T}^{\mu \lambda}).\end{array}$$ | (38) |

This reference further contains a comprehensive analysis of general relativistic hydrodynamics in conservation form.

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
$\mathbf{U}$
. After this update the vector of primitive quantities
$\mathbf{w}$
must be re-evaluated, as those are needed in the Riemann solver (see Section 3.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. Specific details on this issue for each formulation of the equations can be found in Refs. [21, 78, 220] . 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 solving a quartic polynomial. Therefore, iterative methods are still preferred [78] . 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.

Lightcone hydrodynamics: A comprehensive numerical study of the formulation of the hydrodynamic equations discussed in this section was presented in [220] , 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, dynamical spacetime evolutions of relativistic polytropes (i.e., stellar models satisfying the so-called 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 have been presented before in [133] (see [31] for a recent implementation). This approach is geared towards smooth isentropic flows. A Lagrangian method, limited to spherical symmetry, was developed by [180] . Recent work in [71] includes an Eulerian, non-conservative, formulation for general fluids in null hypersurfaces and spherical symmetry, including their matching to a spacelike section.

The general formalism laid out in [220, 218] is currently being systematically applied to astrophysical problems of increasing complexity. Applications in spherical symmetry include the investigation of accreting dynamic black holes, which can be found in [220, 221] . Studies of the gravitational collapse of supermassive stars are discussed in [155] , and studies of the interaction of scalar fields with relativistic stars are presented in [269] . Axisymmetric neutron star spacetimes have been considered in [268] , as part of a broader program aimed at the study of relativistic stellar dynamics and gravitational collapse using characteristic numerical relativity. We note that there has been already a proof-of-principle demonstration of the inclusion of matter fields in three dimensions [31] .

2.3 Going further: Non-ideal hydrodynamics

Formulations of the equations of non-ideal hydrodynamics in general relativity are also available in the literature. The term “non-ideal” accounts for additional physics such as viscosity, magnetic fields, and radiation. These non-adiabatic effects can play a major role in some astrophysical systems as, such as accretion disks or relativistic stars.

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 [191] and references therein). These extended fluid theories, however, remain unexplored, numerically, in astrophysical systems. The reason may be the lack of an appropriate formulation well-suited for numerical studies. Work in this direction was done by Peitz and Appl [223] who provided a 3+1 coordinate-free representation of different types of dissipative relativistic fluid theories which possess, in principle, the potentiality of being well adapted to numerical applications.

The inclusion of magnetic fields and the development of formulations for the MHD equations, attractive to numerical studies, is still very limited in general relativity. Numerical approaches in special relativity are presented in [143, 290, 20] . In particular, Komissarov [143] , and Balsara [20] developed two different upwind HRSC (or Godunov-type) schemes, providing the characteristic information of the corresponding system of equations, and proposed a battery of tests to validate numerical MHD codes. 3+1 representations of relativistic MHD can be found in [271, 81] . In [312] the transport of energy and angular momentum in magneto-hydrodynamical accretion onto a rotating black hole was studied adopting Wilson's formulation for the hydrodynamic equations (conveniently modified to account for the magnetic terms), and the magnetic induction equation was solved using the constrained transport method of [81] . Recently, Koide et al. [141, 142] performed the first MHD simulation, in general relativity, of magnetically driven relativistic jets from an accretion disk around a Schwarzschild black hole (see Section 4.2.2 ). These authors used a second-order finite difference central scheme with nonlinear dissipation developed by Davis [61] . Even though astrophysical applications of Godunov-type schemes (see Section 3.1.2 ) in general relativistic MHD are still absent, it is realistic to believe this situation may change in the near future.

The interaction between matter and radiation fields, present in different levels of complexity in all astrophysical systems, is described by the equations of radiation hydrodynamics. The Newtonian framework is highly developed (see, e.g., [179] ; the special relativistic transfer equation is also considered in that reference). Pons et al. [229] discuss a hyperbolic formulation of the radiative transfer equations, paying particular attention to the closure relations and to extend HRSC schemes to those equations. General relativistic formulations of radiative transfer in curved spacetimes are considered in, e.g., [236] and [315] (see also references therein).

3 Numerical Schemes

We turn now to describe the numerical schemes, mainly those based on finite differences, specifically designed to solve nonlinear hyperbolic systems of conservation laws. As discussed in the previous section, the equations of general relativistic hydrodynamics fall in this category, irrespective of the formulation. Even though we also consider schemes based on artificial viscosity techniques, the emphasis is on the so-called 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., [151, 152, 286, 128] ). For this reason only the most relevant features will be covered here, addressing the reader to the appropriate literature for further details. In particular, an excellent introduction to the implementation of HRSC schemes in special relativistic hydrodynamics is presented in the Living Reviews article by Martí and Müller [165] . Alternative techniques to finite differences, such as smoothed particle hydrodynamics, (pseudo-)spectral methods and others, are briefly considered last.

3.1 Finite difference schemes

Any system of equations presented in the previous section 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
$\mathbf{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 time step constrained by the Courant–Friedrichs–Lewy (CFL) condition.

The hydrodynamic equations (either in Newtonian physics or in general relativity) constitute a nonlinear hyperbolic system and, hence, smooth initial data can transform into discontinuous data (the crossing of characteristics in the case of shocks) in a finite time during the evolution. As a consequence, classical finite difference schemes (see, e.g., [151, 286] ) 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 as the grid is refined. 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.

3.1.1 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 [294] 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
$\Delta x$
of the space variable. In their original work, von Neumann and Richtmyer proposed the following expression for the viscosity term:

$$Q=\{\begin{array}{cc}-\alpha \frac{\partial v}{\partial x}& \text{if}\frac{\partial v}{\partial x}<0\text{or}\frac{\partial \rho}{\partial t}>0,\\ 0& \text{otherwise},\end{array}$$

with
$\alpha =\rho (k\Delta x{)}^{2}\partial v/\partial x$
,
$v$
being the fluid velocity,
$\rho $
the density,
$\Delta x$
the spatial interval, and
$k$
a constant parameter whose value is conveniently 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 conjuction 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 [171] the artificial viscosity term, obtained in analogy with the one proposed by von Neumann and Richtmyer [294] , is introduced in the equations accompanying the pressure, in the form

$$Q=\{\begin{array}{cc}\frac{\rho}{\Gamma}{\left(\frac{a\Delta m}{{R}^{2}}\right)}^{2}\frac{\partial {R}^{2}u}{\partial m}& \text{if}\frac{\partial \rho}{\partial t}>0,\\ 0& \text{otherwise}.\end{array}$$

Further examples of similar expressions for the artificial viscosity terms, in the context of Wilson formulation, can be found in, e.g., [298, 123] . A state-of-the-art formulation of the artificial viscosity approach is reported in [13] .
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 [207, 13] . 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 in the 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 [206] .

3.1.2 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
$\alpha $
is used, the global error of the numerical solution has to tend to zero as
$\mathcal{O}(\Delta x{)}^{\alpha}$
as the cell width
$\Delta x$
tends to zero). For hyperbolic systems of conservation laws, schemes written in conservation form are preferred since, according to the Lax–Wendroff theorem [150] , they guarantee that the convergence, if it exists, is to one of the so-called 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
${\mathcal{C}}^{1}$
classical solutions (continuous and differentiable) in regions where they are continuous and have a finite number of discontinuities.

For the sake of simplicity let us consider in the following an initial value problem for a one-dimensional scalar hyperbolic conservation law,

and let us introduce a discrete numerical grid of space-time 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

where
$\hat{f}$
is a consistent numerical flux function (i.e.,
$\hat{f}(u,u,\cdots ,u)=f\left(u\right)$
) of
$p+q+1$
arguments and
$\Delta t$
and
$\Delta x$
are the time step and cell width respectively. Furthermore,
${u}_{j}^{n}$
is an approximation of the average of
$u(x,t)$
within the numerical cell
$[{x}_{j-1/2},{x}_{j+1/2}]$
$({x}_{j\pm 1/2}=({x}_{j}+{x}_{j\pm 1})/2)$
:

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
$\eta \to 0$
, of the “viscous version” of the scalar conservation law, Equation ( 39 ):

Mathematically, the solution one needs to search for is characterized by the so-called entropy condition (in the language of fluids, the condition that the entropy of any fluid element should increase when running into a discontinuity). The characterization of these entropy-satisfying solutions for scalar equations was given by Oleinik [211] . For hyperbolic systems of conservation laws it was developed by Lax [149] .

$$\begin{array}{c}\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}=0,u(x,t=0)={u}_{0}\left(x\right),\end{array}$$ | (39) |

$$\begin{array}{c}{u}_{j}^{n+1}={u}_{j}^{n}-\frac{\Delta t}{\Delta x}\left(\hat{f}\right({u}_{j-p}^{n},{u}_{j-p+1}^{n},\cdots ,{u}_{j+q}^{n})-\hat{f}({u}_{j-p-1}^{n},{u}_{j-p}^{n},\cdots ,{u}_{j+q-1}^{n}\left)\right),\end{array}$$ | (40) |

$$\begin{array}{c}{u}_{j}^{n}\approx \frac{1}{\Delta x}{\int}_{{x}_{j-1/2}}^{{x}_{j+1/2}}u(x,{t}^{n})dx.\end{array}$$ | (41) |

$$\begin{array}{c}\frac{\partial u}{\partial t}+\frac{\partial f\left(u\right)}{\partial x}=\eta \frac{{\partial}^{2}u}{\partial {x}^{2}}.\end{array}$$ | (42) |

The Lax–Wendroff theorem [150] 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., [240] ). Along this direction, the notion of total-variation 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
$\left({u}^{n}\right)$
, is defined as

A numerical scheme is said to be TV-stable if TV
$\left({u}^{n}\right)$
is bounded for all
$\Delta t$
at any time for each initial data. In the case of nonlinear, scalar conservation laws it can be proved that TV-stability is a sufficient condition for convergence [151] , 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 so-called total variation diminishing (TVD) schemes [115] (see below).

$$\begin{array}{c}\text{TV}\left({u}^{n}\right)={\sum}_{j=0}^{+\infty}|{u}_{j+1}^{n}-{u}_{j}^{n}|.\end{array}$$ | (43) |

Let us now consider the specific system of hydrodynamic equations as formulated in Equation ( 30 ), and let us consider a single computational cell of our discrete spacetime. Let
$\Omega $
be a region (simply connected) of a given four-dimensional manifold
$\mathcal{\mathcal{M}}$
, bounded by a closed three-dimensional surface
$\partial \Omega $
. We further take the 3-surface
$\partial \Omega $
as the standard-oriented hyper-parallelepiped made up of two spacelike surfaces
$\{{\Sigma}_{{x}^{0}},{\Sigma}_{{x}^{0}+\Delta {x}^{0}}\}$
plus timelike surfaces
$\{{\Sigma}_{{x}^{i}},{\Sigma}_{{x}^{i}+\Delta {x}^{i}}\}$
that join the two temporal slices together. By integrating system ( 30 ) over a domain
$\Omega $
of a given spacetime, the variation in time of the state vector
$\mathbf{U}$
within
$\Omega $
is given – keeping apart the source terms – by the fluxes
${\mathbf{F}}^{i}$
through the boundary
$\partial \Omega $
. The integral form of system ( 30 ) is

which can be written in the following conservation form, well-adapted to numerical applications:

where

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.

$$\begin{array}{c}{\int}_{\Omega}\frac{1}{\sqrt{-g}}\frac{\partial \sqrt{\gamma}\mathbf{U}}{\partial {x}^{0}}d\Omega +{\int}_{\Omega}\frac{1}{\sqrt{-g}}\frac{\partial \sqrt{-g}{\mathbf{F}}^{i}}{\partial {x}^{i}}d\Omega ={\int}_{\Omega}\mathbf{S}d\Omega ,\end{array}$$ | (44) |

$$\begin{array}{ccc}& & (\overline{\mathbf{U}}\Delta V{)}_{{x}^{0}+\Delta {x}^{0}}-(\overline{\mathbf{U}}\Delta V{)}_{{x}^{0}}=\end{array}$$ |

$$\begin{array}{ccc}& & -\left({\int}_{{\Sigma}_{\text{}x{\text{}}^{1}+\Delta \text{}x{\text{}}^{1}}}\sqrt{-g}{\mathbf{F}}^{1}d{x}^{0}d{x}^{2}d{x}^{3}-{\int}_{{\Sigma}_{\text{}x{\text{}}^{1}}}\sqrt{-g}{\mathbf{F}}^{1}d{x}^{0}d{x}^{2}d{x}^{3}\right)\end{array}$$ |

$$\begin{array}{ccc}& & -\left({\int}_{{\Sigma}_{\text{}x{\text{}}^{2}+\Delta \text{}x{\text{}}^{2}}}\sqrt{-g}{\mathbf{F}}^{2}d{x}^{0}d{x}^{1}d{x}^{3}-{\int}_{{\Sigma}_{\text{}x{\text{}}^{2}}}\sqrt{-g}{\mathbf{F}}^{2}d{x}^{0}d{x}^{1}d{x}^{3}\right)\end{array}$$ |

$$\begin{array}{ccc}& & -\left({\int}_{{\Sigma}_{\text{}x{\text{}}^{3}+\Delta \text{}x{\text{}}^{3}}}\sqrt{-g}{\mathbf{F}}^{3}d{x}^{0}d{x}^{1}d{x}^{2}-{\int}_{{\Sigma}_{\text{}x{\text{}}^{3}}}\sqrt{-g}{\mathbf{F}}^{3}d{x}^{0}d{x}^{1}d{x}^{2}\right)\end{array}$$ |

$$\begin{array}{ccc}& & +{\int}_{\Omega}\mathbf{S}d\Omega ,\end{array}$$ | (45) |

$$\begin{array}{c}\overline{\mathbf{U}}=\frac{1}{\Delta V}{\int}_{\Delta V}\sqrt{\gamma}\mathbf{U}d{x}^{1}d{x}^{2}d{x}^{3},\end{array}$$ | (46) |

$$\begin{array}{c}\Delta V={\int}_{{x}^{1}}^{{x}^{1}+\Delta {x}^{1}}{\int}_{{x}^{2}}^{{x}^{2}+\Delta {x}^{2}}{\int}_{{x}^{3}}^{{x}^{3}+\Delta {x}^{3}}\sqrt{\gamma}d{x}^{1}d{x}^{2}d{x}^{3}.\end{array}$$ | (47) |

The computation of the time integrals of the interface fluxes appearing in Equation ( 45 ) is one of the distinctive features of upwind HRSC schemes. One needs first to define the so-called numerical fluxes, which are recognized as approximations to the time-averaged fluxes across the cell interfaces, which depend on the solution at those interfaces,
$\mathbf{U}({x}^{j}+\Delta {x}^{j}/2,{x}^{0})$
, during a time step,

At the cell interfaces the flow can be discontinuous and, following the seminal idea of Godunov [108] , 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 so-called Godunov-type methods [117, 75] .

$$\begin{array}{c}{\hat{\mathbf{F}}}_{j+\frac{1}{2}}\approx \frac{1}{\Delta t}{\int}_{{t}^{n}}^{{t}^{n+1}}\mathbf{F}(\mathbf{U}({x}^{j+1/2},{x}^{0}\left)\right).\end{array}$$ | (48) |

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: a shock wave, a rarefaction wave, and a contact discontinuity.

The complete structure of the Riemann problem can be solved analytically (see [108] for the solution in Newtonian hydrodynamics and [163, 230] in special relativistic hydrodynamics) and, accordingly, used to update the solution forward in time.
For reasons of numerical efficiency and, particularly in multi-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. After extensive experimentation, the results achieved with approximate Riemann solvers are comparable to those obtained with the exact solver (see [286] for a comprehensive overview of Godunov-type methods, and [165] for an excellent summary of approximate Riemann solvers in special relativistic hydrodynamics).

In the frame of the local characteristic approach, the numerical fluxes appearing in Equation ( 45 ) 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 [241] it adopts the following functional form:

where
${\mathbf{w}}_{L}$
and
${\mathbf{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.

$$\begin{array}{c}{\hat{\mathbf{F}}}_{j+\frac{1}{2}}=\frac{1}{2}\left(\mathbf{F}\left({\mathbf{w}}_{R}\right)+\mathbf{F}\left({\mathbf{w}}_{L}\right)-{\sum}_{n=1}^{5}\left|{\stackrel{~}{\lambda}}_{n}\right|\Delta {\stackrel{~}{\omega}}_{n}{\stackrel{~}{\mathbf{r}}}_{n}\right),\end{array}$$ | (49) |

The way these 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., [286, 152] ) most commonly used for TVD schemes [115] are the second order, piecewise-linear reconstruction, introduced by van Leer [289] in the design of the MUSCL scheme (Monotonic Upstream Scheme for Conservation Laws), and the third order, piecewise parabolic reconstruction developed by Colella and Woodward [58] 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 developed. Among those, we mention the total variation bounded (TVB) schemes [267] and the essentially non-oscillatory (ENO) schemes [116] .

Alternatively, high-order methods for nonlinear hyperbolic systems have also been designed using flux limiters rather than slope limiters, as in the FCT scheme of Boris and Book [46] . 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,
$\hat{\mathbf{F}}={\hat{\mathbf{F}}}_{h}-(1-\Phi )({\hat{\mathbf{F}}}_{h}-{\hat{\mathbf{F}}}_{l})$
with the limiter
$\Phi \in [0,1]$
(see [286, 152] for further details).

The last term in the flux-formula, Equation ( 49 ), 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 ( 49 ),
$\{{\stackrel{~}{\lambda}}_{n},{\stackrel{~}{\mathbf{r}}}_{n}{\}}_{n=1,...,5}$
are the eigenvalues and right-eigenvectors of the Jacobian matrix of the flux vector, respectively. Correspondingly, the quantities
$\{\Delta {\stackrel{~}{\omega}}_{n}{\}}_{n=1,...,5}$
are the jumps of the so-called characteristic variables across each characteristic field. They are obtained by projecting the jumps of the state vector variables with the left-eigenvectors matrix:

The “tilde” in Equations ( 49 ) and ( 50 ) indicates that the corresponding fields are averaged at the cell interfaces from the left and right (reconstructed) values.

$$\begin{array}{c}\mathbf{U}\left({\mathbf{w}}_{R}\right)-\mathbf{U}\left({\mathbf{w}}_{L}\right)={\sum}_{n=1}^{5}\Delta {\stackrel{~}{\omega}}_{n}{\stackrel{~}{\mathbf{r}}}_{n}.\end{array}$$ | (50) |

During the last few years most of the standard Riemann solvers developed in Newtonian fluid dynamics have been extended to relativistic hydrodynamics: Eulderink [78] , as discussed in Section 2.2.1 , explicitly derived a relativistic Roe Riemann solver [241] ; Schneider et al. [249] carried out the extension of Harten, Lax, van Leer, and Einfeldt's (HLLE) method [117, 75] ; Martí and Müller [164] extended the PPM method of Woodward and Colella [305] ; Wen et al. [296] extended Glimm's exact Riemann solver; Dolezal and Wong [68] put into practice Shu–Osher ENO techniques; Balsara [19] extended Colella's two-shock approximation, and Donat et al. [69] extended Marquina's method [70] . Recently, much effort has been spent concerning the exact special relativistic Riemann solver and its extension to multi-dimensions [163, 230, 237, 238] . The interested reader is addressed to [165] and references therein for a comprehensive description of such solvers in special relativistic hydrodynamics.

3.1.3 High-order central schemes

The use of high-order non-oscillatory symmetric (central) TVD schemes for solving hyperbolic systems of conservation laws emerged at the mid 1980s [61, 242, 310, 204] (see also [311] and [286] and references therein) as an alternative approach to upwind HRSC schemes. Symmetric schemes are based on either high-order schemes (e.g., Lax–Wendroff ) with additional dissipative terms [61, 242, 310] , or on non-oscillatory high-order extensions of first-order central schemes (e.g., Lax–Friedrichs) [204] .

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 multi-dimensional problems. Its use is, thus, not limited to those systems of equations where the characteristic information is explicitly known or to systems where the solution of the Riemann problem is prohibitively expensive to compute. This approach has gradually developed during the last decade to reach a mature status where a number of straightforward 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 [286] . An up-to-date summary of the status and applications of this approach is discussed in [286, 145, 280] .

In the context of special and general relativistic MHD, Koide et al. [141, 142] applied a second-order central scheme with nonlinear dissipation developed by Davis [61] to the study of black hole accretion and formation of relativistic jets. One-dimensional test simulations in special relativistic hydrodynamics performed by the author and coworkers [92] using the SLIC (slope limiter centred) scheme (see [286] for details) showed its capabilities to yield satisfactory results, comparable to those obtained by HRSC schemes based on Riemann solvers, even well inside the ultrarelativistic regime. The slopes of the original central scheme were limited using high-order reconstruction procedures such as PPM [58] , which was essential to keep the inherent diffusion of central schemes at discontinuities at reasonable levels. Very recently, Del Zanna and Bucciantini [63] assessed a third-order convex essentially non-oscillatory central scheme in multi-dimensional special relativistic hydrodynamics. Again, 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 assessed by [13] in one-dimensional special and general relativistic hydrodynamics, where similar results to those of [63] are reported. These authors also validate their central scheme in simulations of spherical accretion onto a Schwarzschild black hole, and further provide a comparison with an artificial viscosity based scheme.

It is worth emphasizing that early pioneer approaches in the field of relativistic hydrodynamics [207, 54] 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\sim 2$
). A key feature lacking in those investigations was the use of a conservative approach for both the system of equations and the numerical schemes. A posteriori, and in the light of the results reported by [63, 13, 92] , it appears that this was the ultimate reason preventing the extension of the early computations to the ultrarelativistic regime.

The alternative of using high-order central schemes instead of upwind HRSC schemes becomes apparent when the spectral decomposition of the hyperbolic system under study is not known. The straightforwardness of a central scheme makes its use very appealing, especially in multi-dimensions where computational efficiency is an issue. Perhaps the most important example in relativistic astrophysics is the system of (general) relativistic MHD equations. Despite some progress in recent years (see, e.g., [20, 143] ), much more work is needed concerning their solution with HRSC schemes based upon Riemann solvers. Meanwhile, an obvious choice is the use of central schemes [141, 142] .

3.1.4 Source terms

Most “conservation laws” involve source terms on the right hand side of the equations. In hydrodynamics, for instance, those terms arise when considering external forces such as gravity, which make the right hand side 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) or ionization (resulting in a collection of non-homogeneous 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, addressing the interested reader to [286, 152] and references therein for further details.

There are, essentially, two basic 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 ( 45 )):

The temporal order of this basic algorithm can be improved by introducing successive sub-steps 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,
${\mathbf{U}}^{n+1}={\mathcal{\mathcal{L}}}_{s}^{\Delta t}{\mathcal{\mathcal{L}}}_{f}^{\Delta t}{\mathbf{U}}^{n}$
, the operator
${\mathcal{\mathcal{L}}}_{f}^{\Delta t}$
solves the homogeneous PDE in the first step to yield the intermediate value
${\mathbf{U}}^{*}$
. Then, in the second step, the operator
${\mathcal{\mathcal{L}}}_{s}^{\Delta t}$
solves the ordinary differential equation
${\partial}_{t}{\mathbf{U}}^{*}=\mathbf{S}\left({\mathbf{U}}^{*}\right)$
to yield the final value
${\mathbf{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, where
${\mathbf{U}}^{n+1}={\mathcal{\mathcal{L}}}_{s}^{\Delta t/2}{\mathcal{\mathcal{L}}}_{f}^{\Delta t}{\mathcal{\mathcal{L}}}_{s}^{\Delta t/2}{\mathbf{U}}^{n}$
. Therefore, this method advances by a half time step the solution for the ODE containing the source terms, and by a full time step the conservation law (the PDE) in between each source step.

$$\begin{array}{c}{\mathbf{U}}_{j}^{n+1}={\mathbf{U}}_{j}^{n}-\frac{\Delta t}{\Delta x}\left({\hat{\mathbf{F}}}_{j+\frac{1}{2}}-{\hat{\mathbf{F}}}_{j-\frac{1}{2}}\right)+\Delta t{\mathbf{S}}_{j}^{n}.\end{array}$$ | (51) |

We note that in some cases the source terms may become stiff, as in phenomena that either occur on much faster timescales than the hydrodynamic time step
$\Delta t$
, or act over much smaller spatial scales than the grid resolution
$\Delta x$
. Stiff source terms may easily lead to numerical difficulties.

The interested reader is directed to [152] (and references therein) for further information on various approaches to overcome the problems of stiff source terms.

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

3.2.1 Smoothed particle hydrodynamics

The Lagrangian particle method SPH, derived independently by Lucy [157] and Gingold and Monaghan [104] , has shown successful performance to model 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 hydrodynamical variables, the finite extent of the particles being determined by a smoothing function (the kernel) containing a characteristic length scale
$h$
. The main advantage of this 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 multi-dimensional applications, when, e.g., modelling regions containing large voids. Experience in Newtonian hydrodynamics shows that SPH produces very accurate results with a small number of particles (
$\approx {10}^{3}$
or even less). However, if more than
${10}^{4}$
particles have to be used for certain problems, and self-gravity has to be included, the computational power, which grows as the square of the number of particles, may exceed the capabilities of current supercomputers.

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.

Reviews of the classical SPH equations are abundant in the literature (see, e.g., [186, 190] and references therein). The reader is addressed to [190] for a summary of comparisons between SPH and HRSC methods.

Recently, implementations of SPH to handle (special) relativistic (and even ultrarelativistic) flows have been developed (see, e.g., [57] and references therein). However, SPH has been scarcely applied to simulate relativistic flows in curved spacetimes. Relevant references include [137, 146, 147, 270] .

Following [146] , let us describe the implementation of an SPH scheme in general relativity. Given a function
$f(\mathbf{x})$
, its mean smoothed value
$\langle f(\mathbf{x})\rangle ,(\mathbf{x}=(x,y,z\left)\right)$
can be obtained from

where
$W$
is the smoothing kernel,
$h$
the smoothing length, and
$\sqrt{{g}^{\prime}}{d}^{3}{x}^{\prime}$
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 [105] . The kernel is required to satisfy a normalization condition,

which is assured by choosing
$W(\mathbf{x},{\mathbf{x}}^{\prime};h)=\xi (\mathbf{x})\Omega \left(v\right)$
, with
$v=|\mathbf{x}-{\mathbf{x}}^{\prime}|/h$
,
$\xi (\mathbf{x})$
being a normalization function, and
$\Omega \left(v\right)$
a standard spherical kernel.

$$\begin{array}{c}\langle f(\mathbf{x})\rangle \equiv \int W(\mathbf{x},{\mathbf{x}}^{\prime};h)f\left({\mathbf{x}}^{\prime}\right)\sqrt{{g}^{\prime}}{d}^{3}{x}^{\prime},\end{array}$$ | (52) |

$$\begin{array}{c}\int W(\mathbf{x},{\mathbf{x}}^{\prime};h)\sqrt{{g}^{\prime}}{d}^{3}{x}^{\prime}=1,\end{array}$$ | (53) |

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
$\langle n(\mathbf{x})\rangle \equiv {\sum}_{a=1}^{N}\delta (\mathbf{x}-{\mathbf{x}}_{a})/\sqrt{g}$
, with
$\{{\mathbf{x}}_{a}{\}}_{a=1,...,N}$
being the collection of
$N$
particles where the functions are known. The previous representations then read:

with
${\Omega}_{ab}\equiv \Omega ({\mathbf{x}}_{a},{\mathbf{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 [146] . The time evolution of the final system of ODEs is performed in [146] using a second-order Runge–Kutta time integrator with adaptive time step. 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 [159] .

$$\begin{array}{c}\langle \nabla f(\mathbf{x})\rangle =\nabla \langle f(\mathbf{x})\rangle -\langle f(\mathbf{x})\rangle \nabla ln\xi (\mathbf{x}),\end{array}$$ | (54) |

$$\begin{array}{c}\langle \nabla \cdot \mathbf{A}(\mathbf{x})\rangle =\nabla \cdot \langle \mathbf{A}(\mathbf{x})\rangle -\langle \mathbf{A}(\mathbf{x})\rangle \cdot ln\xi (\mathbf{x}).\end{array}$$ | (55) |

$$\begin{array}{ccc}\langle f\left({\mathbf{x}}_{a}\right)\rangle & =& \xi \left({\mathbf{x}}_{a}\right){\sum}_{b=1}^{N}\frac{f\left({\mathbf{x}}_{b}\right)}{\langle n\left({\mathbf{x}}_{b}\right)\rangle}{\Omega}_{ab},\end{array}$$ | (56) |

$$\begin{array}{ccc}\langle \nabla f\left({\mathbf{x}}_{a}\right)\rangle & =& \xi \left({\mathbf{x}}_{a}\right){\sum}_{b=1}^{N}\frac{f\left({\mathbf{x}}_{b}\right)}{\langle n\left({\mathbf{x}}_{b}\right)\rangle}{\nabla}_{{\mathbf{x}}_{a}}{\Omega}_{ab},\end{array}$$ | (57) |

$$\begin{array}{ccc}\langle \nabla \cdot \mathbf{A}\left({\mathbf{x}}_{a}\right)\rangle & =& \xi \left({\mathbf{x}}_{a}\right){\sum}_{b=1}^{N}\frac{\mathbf{A}\left({\mathbf{x}}_{b}\right)}{\langle n\left({\mathbf{x}}_{b}\right)\rangle}\cdot {\nabla}_{{\mathbf{x}}_{a}}{\Omega}_{ab},\end{array}$$ | (58) |

Recently, Siegler and Riffert [270] have developed a Lagrangian conservation form of the general relativistic hydrodynamic equations for perfect fluids (with artificial viscosity) in arbitrary background spacetimes. Within that formulation these authors [270] 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., [159, 137, 146] ). The conservative character of their scheme has allowed the modelling of ultrarelativistic flows including shocks with Lorentz factors as large as 1000.

3.2.2 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 [110, 52] .

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 have been obtained in one-dimensional supernova collapse simulations, both in the framework of a tensor-scalar theory of gravitation [208, 210] and in general relativity [209] .

Following [41] we illustrate the main ideas of spectral methods considering the quasi-linear one-dimensional scalar equation:

with
$u=u(t,x)$
, and
$\lambda $
a constant parameter. In the linear case (
$\lambda =0$
), and assuming the function
$u$
to be periodic, spectral methods expand the function into a Fourier series:

From the numerical point of view, the series is truncated for a suitable value of
$k$
. Hence, Equation ( 59 ), with
$\lambda =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)$
.

$$\begin{array}{c}\frac{\partial u}{\partial t}=\frac{{\partial}^{2}u}{\partial {x}^{2}}+\lambda u\frac{\partial u}{\partial x},t\ge 0,x\in [0,1],\end{array}$$ | (59) |

$$\begin{array}{c}u(x,t)={\sum}_{k=0}^{\infty}\left[{a}_{k}\right(t)cos(2\pi kx)+{b}_{k}(t)sin(2\pi kx\left)\right].\end{array}$$ | (60) |

$$\begin{array}{c}\frac{d{a}_{k}}{dt}=-{k}^{2}{a}_{k}\left(t\right),\frac{d{b}_{k}}{dt}=-{k}^{2}{b}_{k}\left(t\right).\end{array}$$ | (61) |

In the nonlinear case,
$\lambda \ne 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
$\partial u/\partial 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 ( 60 ) 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 [110, 52] 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 the study of compact objects, as well as the associated violent phenomena of gravitational collapse and supernova explosion. They have developed a fully object-oriented library (based on the C++ computer language) called LORENE [156] to implement (multi-domain) spectral methods within spherical coordinates. A comprehensive summary of applications in general relativistic astrophysics is presented in [41] . The most recent ones deal with the computation of quasi-equilibrium configurations of either synchronized or irrotational binary neutron stars in general relativity [112, 283, 282] . Such initial data are currently being used by fully relativistic, finite difference time-dependent codes (see Section 3.3.2 ) to simulate the coalescence of binary neutron stars.

3.2.3 Flow field-dependent variation method

Richardson and Chung [239] have recently proposed 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 general relativistic hydrodynamic equations are expanded in a special form of the Taylor series:

with
${s}_{a}$
and
${s}_{b}$
denoting the first-order and second-order variation parameters. Using the above expressions, the time update then reads:

Combining the conservation form of the equations and the rearranged Taylor series expansion, allows us to rewrite the time update into standard matrix (residual) form, which can then be discretized using either standard finite difference or finite element methods [239] .

$$\begin{array}{c}\begin{array}{ccc}{\mathbf{U}}^{n+1}& =& {\mathbf{U}}^{n}+\Delta t\frac{\partial {\mathbf{U}}^{n+{s}_{a}}}{\partial t}+\frac{\Delta {t}^{2}}{2}\frac{{\partial}^{2}{\mathbf{U}}^{n+{s}_{b}}}{\partial {t}^{2}}+\mathcal{O}(\Delta {t}^{3}),\\ \frac{\partial {\mathbf{U}}^{n+{s}_{a}}}{\partial t}& =& \frac{\partial {\mathbf{U}}^{n}}{\partial t}+{s}_{a}\frac{\partial \Delta {\mathbf{U}}^{n+1}}{\partial t},& 0<{s}_{a}<1,\\ \frac{{\partial}^{2}{\mathbf{U}}^{n+{s}_{b}}}{\partial {t}^{2}}& =& \frac{{\partial}^{2}{\mathbf{U}}^{n}}{\partial {t}^{2}}+{s}_{b}\frac{{\partial}^{2}\Delta {\mathbf{U}}^{n+1}}{\partial {t}^{2}},& 0<{s}_{b}<1,\end{array}\end{array}$$ | (62) |

$$\begin{array}{ccc}{\mathbf{U}}^{n+1}={\mathbf{U}}^{n}& +& \Delta t\left(\frac{\partial {\mathbf{U}}^{n}}{\partial t}+{s}_{a}\frac{\partial \Delta {\mathbf{U}}^{n+1}}{\partial t}\right)\end{array}$$ |

$$\begin{array}{ccc}& +& \frac{\Delta {t}^{2}}{2}\left(\frac{{\partial}^{2}{\mathbf{U}}^{n}}{\partial {t}^{2}}+{s}_{b}\frac{{\partial}^{2}\Delta {\mathbf{U}}^{n+1}}{\partial {t}^{2}}\right)+\mathcal{O}(\Delta {t}^{3}).\end{array}$$ | (63) |

The physical interpretation of the coefficients
${s}_{a}$
and
${s}_{b}$
is the foundation of the FDV method.

The first-order parameter
${s}_{a}$
is proportional to
${\mathbf{a}}_{i}\partial {\mathbf{U}}^{n+1}/\partial {x}_{i}$
, where
${\mathbf{a}}_{i}$
is the convection Jacobian representing the change of convective motion. If the Lorentz factor remains constant in space and time, then
${s}_{a}=0$
. However, if the Lorentz factor between adjacent zones is large,
${s}_{a}=1$
.

Similar assessments in terms of the Reynolds number can be provided for the diffusion and diffusion gradients, while the Froude number is used for the source term Jacobian
$\partial \mathbf{S}/\partial \mathbf{U}$
. Correspondingly, the second-order FDV parameters
${s}_{b}$
are chosen to be exponentially proportional to the first-order ones.

Obviously, the main drawback of the FDV method is the dependence of the solution procedure on a large number of problem-dependent parameters, a drawback shared to some extent by artificial viscosity schemes. Richardson and Chung [239] have implemented the FDV method in a C++ code called GRAFSS (General Relativistic Astrophysical Flow and Shock Solver). The test problems they report are the special relativistic shock tube (problem 1 in the notation of [164] ) and the Bondi accretion onto 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. Applications to non-ideal hydrodynamics and relativistic MHD are in progress [239] .

3.2.4 Going further

The finite element method is the preferred choice to solve multi-dimensional structural engineering problems since the late 1960s [317] . First steps in the development of the finite element method for modeling astrophysical flows in general relativity are discussed by Meier [175] . 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 [175] 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 to accurately compute the equilibrium stellar structure of Newtonian polytropes, either spherical or rotating. The main limitation of the finite element method, as Meier points out [175] , is that it is generally fully implicit, which results in severe computer time and memory limitations for threeand four-dimensional problems.

3.3 State-of-the-art three-dimensional codes

The most advanced time-dependent, finite difference, three-dimensional Cartesian codes to solve the system of Einstein and hydrodynamics equations are those developed by Shibata [257] and by the Washington University/NCSA/AEI-Golm Numerical Relativity collaboration [96, 89] . The main difference between both codes lies in the numerical methods used to solve the relativistic hydrodynamic equations: artificial viscosity in Shibata's code [257] , and upwind HRSC schemes in the code of [96, 89] . We note, however, that very recently Shibata has upgraded his code to incorporate HRSC schemes (in particular, a Roe-type approximate Riemann solver and piecewise parabolic cell-reconstruction procedures) [259] . Further 3D codes are currently being developed by a group in the U.S. (Duez et al. [73] ) and by a E.U. Research Training Network collaboration [119] .

3.3.1 Shibata

In Shibata's code [257] , the hydrodynamics equations are formulated following Wilson's approach (Section 2.1.2 ) for a conformal-traceless reformulation of the spacetime variables (see below). An important difference with respect to the original system, Equations ( 20 , 21 , 22 ), is that an equation for the entropy is solved instead of the energy equation. The hydrodynamic equations are integrated using van Leer's [289] second order finite difference scheme with artificial viscosity, following the approach of a precursor code developed by [213] .

The ADM Einstein equations are reformulated into a conformal traceless system, an idea originally introduced by Shibata and Nakamura [263] (see also [196] ), and further developed by Baumgarte and Shapiro [25] . This “BSSN” formulation, which shows enhanced long-term stability compared to the original ADM system, makes use of a conformal decomposition of the 3-metric,
${\stackrel{~}{\gamma}}_{ij}={e}^{-4\phi}{\gamma}_{ij}$
and the trace-free part of the extrinsic curvature,
${A}_{ij}={K}_{ij}-{\gamma}_{ij}K/3$
, with the conformal factor
$\phi $
chosen to satisfy
${e}^{4\phi}={\gamma}^{1/3}\equiv det({\gamma}_{ij}{)}^{1/3}$
. In this formulation, as shown by [25] , in addition to the evolution equations for the conformal 3–metric
${\stackrel{~}{\gamma}}_{ij}$
and the conformal-traceless extrinsic curvature variables
${\stackrel{~}{A}}_{ij}$
, there are evolution equations for the conformal factor
$\phi $
, the trace of the extrinsic curvature
$K$
and the “conformal connection functions”
${\stackrel{~}{\Gamma}}^{i}$
. Further details can be found in [25, 257] .

The code uses an approximate maximal slicing condition, which amounts to solving a parabolic equation for
$ln\alpha $
, and a minimal distortion gauge condition for the shift vector. It also admits
$\pi $
-rotation symmetry around the
$z$
-axis, as well as plane symmetry with respect to the
$z=0$
plane, allowing computations in a quadrant region. In addition, Shibata has also implemented in the code the “cartoon” method proposed by the AEI Numerical Relativity group [7] , which makes possible axisymmetric computations with a Cartesian grid. “Approximate” outgoing boundary conditions are used at the outer boundaries; these do not completely eliminate numerical errors due to spurious back reflection of gravitational waves^{
$\text{1}$
}
. A staggered leapfrog method is used for the time evolution of the metric quantities. Correspondingly, the hydrodynamic equations are updated using a second-order two-step Runge–Kutta scheme. In each time step, the staggered metric quantities needed for the hydrodynamics update are properly extrapolated to intermediate time levels to reach the desired order of accuracy.

The code developed by Shibata [257, 259] has been tested in a variety of problems, including spherical collapse of dust to a black hole, signalled by the appearance of the apparent horizon (comparing 1D and 3D simulations), stability of spherical stars and computation of the radial oscillation period, quadrupole oscillations of perturbed spherical stars and computation of the associated gravitational radiation, preservation of the rotational profile of (approximate) rapidly rotating stars, and the preservation of a co-rotating binary neutron star in a quasi-equilibrium state (assuming a conformally flat 3-metric) for more than one orbital period. Improvements of the hydrodynamical schemes have been considered very recently [259] , in order to tackle problems in which shocks play an important role, e.g., stellar core collapse to a neutron star. 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, and, perhaps most importantly, the coalescence of binary neutron stars, a long-standing problem in numerical relativistic hydrodynamics. These applications are discussed in Section 4 . The most recent simulations of binary neutron star coalescence [266] have been performed on a FACOM VPP5000 computer with typical grid sizes of (505, 505, 253) in
$(x,y,z)$
.

^{
$\text{1}$
}
We note however that all codes based on the 3+1 formalism share this problem since the outer boundaries are located at finite radii.

Further work on the development of more sophisticated boundary conditions is needed to solve this problem. Alternative solutions are to follow the light-cone approach developed by Winicour et al. [304] or the conformal formalism of Friedrich [99]

3.3.2 CACTUS/GR_ASTRO

A three-dimensional general relativistic hydrodynamics code developed by the Washington University/NCSA/AEI-Golm collaboration for the NASA Neutron Star Grand Challenge Project [295] is discussed in Refs. [96, 89] . The code is built upon the Cactus Computational Toolkit [170] . A version of this code that passed the milestone requirement of the NASA Grand Challenge project was released to the community. This code has been benchmarked at over 140 GFlop/sec on a
$1024$
node Cray T3E, with a scaling efficiency of over 95%, showing the potential for large scale 3D simulations of realistic astrophysical systems.

The hydrodynamics part of the code uses the conservative formulation discussed in Section 2.1.3 .

A variety of state-of-the-art Riemann solvers are implemented, including a Roe-like solver [241] and Marquina's flux formula [70] . The code uses slope-limiter methods to construct second-order TVD schemes by means of monotonic piecewise linear reconstructions of the cell-centered quantities for the computation of the numerical fluxes. The standard “minmod” limiter and the “monotonized central-difference” limiter [288] are implemented. Both schemes provide the desired second-order accuracy for smooth solutions, while still satisfying the TVD property. In addition, third-order piecewise parabolic (PPM) reconstruction has been recently implemented and used in [278] .

The Einstein equations are solved using the following different approaches: (i) the standard ADM formalism, (ii) a hyperbolic formulation developed by [40] , and (iii) the BSSN formulation of [196, 263, 25] . The time evolution of both the ADM and the BSSN systems can be performed using several numerical schemes [96, 6, 89] . Currently, a leapfrog (non-staggered in time), and an iterative Crank–Nicholson scheme have been coupled to the hydrodynamics solver. The code is designed to handle arbitrary shift and lapse conditions, which can be chosen as appropriate for a given spacetime simulation. The AEI Numerical Relativity group [169] is currently developing robust gauge conditions for (vacuum) black hole spacetimes (see, e.g., [5] and references therein). Hence, all advances accomplished here can be incorporated in future versions of the code for non-vacuum spacetime evolutions. Similarly, since it is a general purpose code, a number of different boundary conditions can be imposed for either the spacetime variables or for the hydrodynamical variables.

We refer the reader to [96, 6, 89] for additional details.

The code has been subjected to a series of convergence tests [96] , with many different combinations of the spacetime and hydrodynamics finite differencing schemes, demonstrating the consistency of the discrete equations with the differential equations. The simulations performed in [96] include, among others, the evolution of equilibrium configurations of compact stars (solutions to the TOV equations) and the evolution of relativistically boosted TOV stars (
$v=0.87c$
) transversing diagonally across the computational domain – a test for which an exact solution exists. In [6, 4] the improved stability properties of the BSSN formulation of the Einstein equations were compared to the ADM system. In particular, in [6] a number of strongly gravitating systems were simulated, including weak and strong gravitational waves, black holes, boson stars and relativistic stars. While the error growth-rate can be decreased by going to higher grid resolutions, the BSSN formulation requires grid resolutions higher than the ones needed in the ADM formulation to achieve the same accuracy. Furthermore, it was shown in [89] that the code successfully passed stringent long-term evolution tests, such as the evolution of both spherical and rapidly rotating, stationary stellar configurations, and the formation of an apparent horizon during the collapse of a relativistic star to a black hole. The high accuracy of the hydrodynamical schemes employed has allowed the detailed investigation of the frequencies of radial, quasi-radial and quadrupolar oscillations of relativistic stellar models, and use them as a strong assessment of the accuracy of the code. The frequencies obtained have been compared to the frequencies computed with perturbative methods and in axisymmetric nonlinear evolutions [88] . In all of the cases considered, the code reproduces these results with excellent accuracy and is able to extract the gravitational waveforms that might be produced during non-radial stellar pulsations.

4 Hydrodynamical Simulations in Relativistic Astrophysics

With the exception of the vacuum two-body problem (i.e., the coalescence of two black holes), all realistic astrophysical systems and sources of gravitational radiation involve matter. Not surprisingly, 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 body of numerical investigations in the literature dealing with hydrodynamical integrations in static background spacetimes. Most of those 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 has begun in more recent years.

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

Three-dimensional simulations have only started more recently. Much effort is nowadays focused on the study of the coalescence of compact neutron star binaries (as well as the vacuum black hole binary counterpart). These theoretical investigations are driven by the emerging possibility of soon detecting gravitational waves with the different experimental efforts currently underway.

The waveform catalogues resulting from time-dependent hydrodynamical simulations may provide some help to data analysis groups, since the chances for detection may be enhanced through matched-filtering techniques.

In the following, we review the status of the numerical investigations in three astrophysical scenarios all involving strong gravitational fields and, hence, relativistic physics: gravitational collapse, accretion onto black holes, and hydrodynamical evolution 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 [12] and references therein.

4.1 Gravitational collapse

The study of gravitational collapse of massive stars has been largely pursued numerically over the years. This problem was already the main motivation when May and White built the first one-dimensional code [171, 172] . Such codes are designed to integrate the coupled system of Einstein field equations and relativistic hydrodynamics, sidestepping Newtonian approaches.

By browsing through the literature one realizes that the numerical study of gravitational collapse has been neatly divided between two main communities since the early investigations.

First, the computational astrophysics community has traditionally focused on the physics of the (type II/Ib/Ic) supernova mechanism, i.e., collapse of unstable stellar iron cores followed by a bounce and explosion. Nowadays, numerical simulations are highly developed, as they include rotation and detailed state-of-the-art microphysics (e.g., EOS and sophisticated treatments for neutrino transport). These studies are currently performed, routinely, in multi-dimensions with advanced HRSC schemes. Progress in this direction has been achieved, however, at the expense of accuracy in the treatment of the hydrodynamics and the gravitational field, by using Newtonian physics. In this approach the emission of gravitational radiation is computed through the quadrupole formula. For reviews of the current status in this direction see [190, 134] and references therein.

On the other hand, the numerical relativity community has been more interested, since the beginning, in highlighting those aspects of the collapse problem having to do with relativistic gravity, in particular, the emission of gravitational radiation from non-spherical collapse (see [205] for a recent Living Reviews article on gravitational radiation from gravitational collapse). Much of the effort has focused on developing reliable numerical schemes to solve the gravitational field equations and much less emphasis, if any, has been given to the details of the microphysics of the collapse. In addition, this approach has mostly considered gravitational collapse leading to black hole formation, employing relativistic hydrodynamics and gravity. Quite surprisingly, both approaches have barely interacted over the years, except for a handful of simulations in spherical symmetry. Nevertheless, numerical relativity is steadily reaching a state in which it is not unrealistic to expect a much closer interaction in the near future (see, e.g., [262, 258, 89, 259] and references therein).

4.1.1 Core collapse supernovae

At the end of their thermonuclear evolution, massive stars in the (Main Sequence) mass range from
$9{M}_{\odot}$
to
$30{M}_{\odot}$
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
$\sim 3\times {10}^{53}erg(M/{M}_{\odot}{)}^{2}(R/10km{)}^{-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. The accepted scenario, which has emerged from numerical simulations, contains the following two necessary ingredients for any plausible explosion mechanism (see [134] for a review): (i) the existence of the shock wave together with neutrino heating that re-energizes it (in the so-called delayed-mechanism by which the stalled prompt supernova shock is reactivated by an increase in the post-shock pressure due to neutrino energy deposition [300, 30] ); and (ii) the convective overturn which rapidly transports energy into the shocked region [59, 29] (and which can lead to large-scale deviations from spherically symmetric explosions).
However, the path to reach such conclusions has been mostly delineated from numerical simulations in one spatial dimension. Fully relativistic simulations of microphysically detailed core collapse off spherical symmetry are still absent and they might well introduce some modifications. The broad picture described above has been demonstrated in multi-dimensions using sophisticated Newtonian models, as is documented in [190] .

Spherically symmetric simulations. May and White's formulation and their corresponding one-dimensional code formed the basis of most spherically symmetric codes built to study the collapse problem. Wilson [297] investigated the effect of neutrino transport on stellar collapse, concluding that, in contrast to previous results, heat conduction by neutrinos does not produce the ejection of material [60, 14] . The code solved the coupled set of hydrodynamic and Einstein equations, supplemented with a Boltzmann transport equation to describe the neutrino flow. Van Riper [292] used a spherically symmetric general relativistic code to study adiabatic collapse of stellar cores, considering the purely hydrodynamical bounce as the preferred explosion mechanism. The important role of general relativistic effects to produce collapses otherwise absent in Newtonian simulations was emphasized. Bruenn [49] developed a code in which the neutrino transport was handled with an approximation, the multigroup flux-limited diffusion. Baron et al. [24] obtained a successful and very energetic explosion for a model in which the value of the incompressibility modulus of symmetric nuclear matter at zero temperature was particularly small,
${K}_{0}^{sym}=180MeV$
(its precise value, nowadays preferred around
$240MeV$
[107] , is still a matter of debate). Mayle et al. [173] computed neutrino spectra from stellar collapse to stable neutron stars for different collapse models using, as in [49] , multigroup flux-limited diffusion for the neutrino transport.

Van Riper [293] and Bruenn [50] verified that a softer supranuclear EOS, combined with general relativistic hydrodynamics, increases the chances for a prompt explosion. In [248, 247] and [177] the neutrino transport was first handled without approximation by using a general relativistic Boltzmann equation. Whereas in [248, 247] the neutrino transport is described by moments of the general relativistic Boltzmann equation in polar slicing coordinates [22] , [177] used maximal slicing coordinates. Miralles et al. [183] , employing a realistic EOS based upon a Hartree–Fock approximation for Skyrme-type nuclear forces at densities above nuclear density, found that the explosion was favoured by soft EOS, a result qualitatively similar to that of Baron et al. [24] , who used a phenomenological EOS. Swesty et al. [279] also focused on the role of the nuclear EOS in stellar collapse on prompt timescales. Contrary to previous results they found that the dynamics of the shock is almost independent of the nuclear incompressibility once the EOS is not unphysically softened as in earlier simulations (e.g., [292, 24, 293, 50, 183] ). Swesty and coworkers used a finite temperature compressible liquid drop model EOS for the nuclear matter component [148] .

For the shock to propagate promptly to a large radius they found that the EOS must be very soft at densities just above nuclear densities, which is inconsistent with the
$1.44{M}_{\odot}$
neutron star mass constraint imposed by observations of the binary pulsar PSR 1913+16.

From the above discussion it is clear that numerical simulations demonstrate a strong sensitivity of the explosion mechanism to the details of the post-bounce evolution: general relativity, the nuclear EOS and the properties of the nascent neutron star, the treatment of the neutrino transport, and the neutrino–matter interaction. Recently, state-of-the-art treatments of the neutrino transport have been achieved, even in the general relativistic case, by solving the Boltzmann equation in connection with the hydrodynamics, instead of using multigroup flux-limited diffusion [176, 231, 232, 153] .

The results of the few spherically symmetric simulations available show, however, that the models do not explode, neither in the Newtonian or in the general relativistic case. Therefore, computationally expensive multi-dimensional simulations with Boltzmann neutrino transport, able to account for convective effects, are needed to draw further conclusions on the viability of the neutrino-driven explosion mechanism.

From the numerical point of view, many of the above investigations used artificial viscosity techniques to handle shock waves. Together with a detailed description of the microphysics, the correct numerical modeling of the shock wave is the major issue in stellar collapse. In this context, the use of HRSC schemes, specifically designed to capture discontinuities, is much more recent, starting in the late 1980s with the Newtonian simulations of Fryxell et al. [103] , which used an Eulerian second-order PPM scheme (see [190] for a review of the present status). There are only a few relativistic simulations, so far restricted to spherical symmetry [129, 243, 210] .

Martí et al. [161] implemented Glaister's approximate Riemann solver [106] in a Lagrangian Newtonian hydrodynamical code. They performed a comparison of the energetics of a stellar collapse simulation using this HRSC scheme and a standard May and White code with artificial viscosity, for the same initial model, grid size and EOS. They found that the simulation employing a Godunov-type scheme produced larger kinetic energies than that corresponding to the artificial viscosity scheme, with a factor of two difference in the maximum of the infalling velocity. Motivated by this important difference, Janka et al. [136] repeated this computation with a different EOS, using a PPM second-order Godunov-type scheme, disagreeing with Martí et al. [161] . The state-of-the-art implementation of the tensorial artificial viscosity terms in [136] , together with the very fine numerical grids employed (unrealistic for three-dimensional simulations), could be the reason for the discrepancies.

As mentioned in Section 2.1.3 , Godunov-type methods were first used to solve the equations of general relativistic hydrodynamics in [162] , where the characteristic fields of the one-dimensional (spherically symmetric) system were derived. The first astrophysical application was stellar collapse.

In [38] the hydrodynamic equations were coupled to the Einstein equations in spherical symmetry.

The field equations were formulated as a first-order flux-conservative hyperbolic system for a harmonic gauge [39] , somehow “resembling” the hydrodynamic equations. HRSC schemes were applied to both systems simultaneously (only coupled through the source terms of the equations).