2.2 Covariant approaches2 Equations of General Relativistic 2 Equations of General Relativistic

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 tex2html_wrap_inline4738, which describes the rate of advance of time along a timelike unit vector tex2html_wrap_inline4740 normal to a surface, and the spacelike shift vector tex2html_wrap_inline4742 that describes the motion of coordinates within a surface.

The line element is written as

  equation174

where tex2html_wrap_inline4734 is the 3-metric induced on each spacelike slice.

2.1.1 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 [172Jump To The Next Citation Point In The Article, 173Jump To The Next Citation Point In The Article]. Building on theoretical work by Misner and Sharp [185Jump To The Next Citation Point In The Article], 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

equation191

m being a radial (Lagrangian) coordinate, indicating the total rest-mass enclosed inside the circumference tex2html_wrap_inline4750 .

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

eqnarray193

In these coordinates the local conservation equation for the baryonic mass, Equation (2Popup Equation), can be easily integrated to yield the metric potential b :

equation200

The gravitational field equations, Equation (10Popup Equation), and the equations of motion, Equation (1Popup Equation), reduce to the following quasi-linear system of partial differential equations (see also [185Jump To The Next Citation Point In The Article]):

  equation207

with the definitions tex2html_wrap_inline4754 and tex2html_wrap_inline4756, satisfying tex2html_wrap_inline4758 . Additionally,

  equation240

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

Hydrodynamics codes based on the original formulation of May and White and on later versions (e.g., [293Jump To The Next Citation Point In The Article]) have been used in many nonlinear simulations of supernova and neutron star collapse (see, e.g., [184Jump To The Next Citation Point In The Article, 280Jump To The Next Citation Point In The Article] and references therein), as well as in perturbative computations of spherically symmetric gravitational collapse within the framework of the linearized Einstein equations [251, 252]. 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 [182]. 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 3+1 Eulerian formulation (Wilson) 

The use of Eulerian coordinates in multi-dimensional numerical relativistic hydrodynamics started with the pioneering work by Wilson [300Jump To The Next Citation Point In The Article]. Introducing the basic dynamical variables D, tex2html_wrap_inline4768, and E, representing the relativistic density, momenta, and energy, respectively, defined as

  equation258

the equations of motion in Wilson's formulation [300Jump To The Next Citation Point In The Article, 301Jump To The Next Citation Point In The Article] are

  equation267

  equation281

  equation306

with the ``transport velocity'' given by tex2html_wrap_inline4772 . We note that in the original formulation [301Jump To The Next Citation Point In The Article] the momentum density equation, Equation (21Popup Equation), is only solved for the three spatial components tex2html_wrap_inline4774, and tex2html_wrap_inline4776 is obtained through the 4-velocity normalization condition tex2html_wrap_inline4778 .

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 [300Jump To The Next Citation Point In The Article] 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 [295Jump To The Next Citation Point In The Article] 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 [195Jump To The Next Citation Point In The Article, 193Jump To The Next Citation Point In The Article, 199, 22Jump To The Next Citation Point In The Article, 276Jump To The Next Citation Point In The Article, 228Jump To The Next Citation Point In The Article, 79Jump To The Next Citation Point In The Article], accretion onto compact objects [122Jump To The Next Citation Point In The Article, 226Jump To The Next Citation Point In The Article], numerical cosmology [53, 54Jump To The Next Citation Point In The Article, 12] and, more recently, the coalescence of neutron star binaries [303Jump To The Next Citation Point In The Article, 304Jump To The Next Citation Point In The Article, 169Jump To The Next Citation Point In The Article]. This formalism has also been employed, in the special relativistic limit, in numerical studies of heavy-ion collisions [302, 175]. 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 (19Popup Equation), with the addition of a multiplicative tex2html_wrap_inline4738 factor (the lapse function) and the introduction of the Lorentz factor, tex2html_wrap_inline4782 :

  equation338

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 [301Jump To The Next Citation Point In The Article], building on a vacuum numerical relativity code specifically developed to investigate the head-on collision of two black holes [273]. The resulting code was axially symmetric and aimed to integrate the coupled set of equations in the context of stellar core collapse [82].

More recently, Wilson's formulation has been applied to the numerical study of the coalescence of binary neutron stars in general relativity [303Jump To The Next Citation Point In The Article, 304Jump To The Next Citation Point In The Article, 169Jump To The Next Citation Point In The Article] (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 tex2html_wrap_inline4734) is conformally flat . The line element, Equation (11Popup Equation), then reads

equation351

The curvature of the 3-metric is then described by a position dependent conformal factor tex2html_wrap_inline4786 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.

  

Click on thumbnail to view image

Figure 1: Results for the shock heating test of a cold, relativistically inflowing gas against a wall using the explicit Eulerian techniques of Centrella and Wilson [54Jump To The Next Citation Point In The Article]. The plot shows the dependence of the relative errors of the density compression ratio versus the Lorentz factor W for two different values of the adiabatic index of the flow, tex2html_wrap_inline4616 (triangles) and tex2html_wrap_inline4618 (circles) gases. The relative error is measured with respect to the average value of the density over a region in the shocked material. The data are from Centrella and Wilson [54Jump To The Next Citation Point In The Article] and the plot reproduces a similar one from Norman and Winkler [208Jump To The Next Citation Point In The Article].

Wilson's formulation showed some limitations in handling situations involving ultrarelativistic flows (tex2html_wrap_inline4794), as first pointed out by Centrella and Wilson [54Jump To The Next Citation Point In The Article]. Norman and Winkler [208Jump To The Next Citation Point In The Article] performed a comprehensive numerical assessment of such formulation by means of special relativistic hydrodynamical simulations. Figure  1 reproduces a plot from [208Jump To The Next Citation Point In The Article] 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 [54Jump To The Next Citation Point In The Article]. This figure shows that for Lorentz factors of about 2 (tex2html_wrap_inline4798), 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 .

Norman and Winkler [208Jump To The Next Citation Point In The Article] 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 (21Popup Equation), and at the divergence of the velocity in the source of the energy equation, Equation (22Popup Equation). However, [208Jump To The Next Citation Point In The Article] 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:

  equation373

In this way, for instance, the momentum equation takes the following form (in flat spacetime):

  equation380

In Wilson's original formulation, Q is omitted in the two terms containing the quantity tex2html_wrap_inline4808 . In general, Q is a nonlinear function of the velocity and, hence, the quantity tex2html_wrap_inline4812 in the momentum density of Equation (26Popup Equation) 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. 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 [13Jump To The Next Citation Point In The Article] 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 [221Jump To The Next Citation Point In The Article, 218Jump To The Next Citation Point In The Article]; 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 tex2html_wrap_inline4814 . 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 tex2html_wrap_inline4816 % for the central scheme.

2.1.3 3+1 conservative Eulerian formulation (Ibáñez and coworkers) 

In 1991, Martí, Ibáñez, and Miralles [163Jump To The Next Citation Point In The Article] 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 [300Jump To The Next Citation Point In The Article, 301Jump To The Next Citation Point In The Article], or implicit schemes as proposed in [208Jump To The Next Citation Point In The Article].

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., [152Jump To The Next Citation Point In The Article]). 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 [163Jump To The Next Citation Point In The Article].

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 [10Jump To The Next Citation Point In The Article] and references therein) by applying Friedrichs' definition of hyperbolicity [100] to a quasi-linear form of the system of hydrodynamic equations,

  equation406

where tex2html_wrap_inline4820 are the Jacobian matrices of the system and tex2html_wrap_inline4822 is a suitable set of primitive (physical) variables (see below). The system (27Popup Equation) is hyperbolic in the time direction defined by the vector field tex2html_wrap_inline4824 with tex2html_wrap_inline4826, if the following two conditions hold: (i) tex2html_wrap_inline4828 and (ii) for any tex2html_wrap_inline4700 such that tex2html_wrap_inline4832, tex2html_wrap_inline4834, the eigenvalue problem tex2html_wrap_inline4836 has only real eigenvalues tex2html_wrap_inline4838 and a complete set of right-eigenvectors tex2html_wrap_inline4840 . 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 tex2html_wrap_inline4842 . In Font et al. [93Jump To The Next Citation Point In The Article] those calculations were extended to an arbitrary reference frame in which the motion of the fluid is described by the 4-velocity tex2html_wrap_inline4844 .

The approach followed in [93] for the equations of special relativistic hydrodynamics was extended to general relativity in [21Jump To The Next Citation Point In The Article]. The choice of evolved variables (conserved quantities) in the 3+1 Eulerian formulation developed by Banyuls et al. [21Jump To The Next Citation Point In The Article] differs slightly from that of Wilson's formulation [300Jump To The Next Citation Point In The Article]. It comprises the rest-mass density (D), the momentum density in the j -direction (tex2html_wrap_inline4850), 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 [21Jump To The Next Citation Point In The Article] for more complete definitions and geometrical foundations.

In terms of the so-called primitive variables tex2html_wrap_inline4854, the conserved quantities are written as

  equation451

where the contravariant components tex2html_wrap_inline4856 of the 3-velocity are defined as

equation458

and W is the relativistic Lorentz factor tex2html_wrap_inline4860 with tex2html_wrap_inline4862 .

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

  equation471

with tex2html_wrap_inline4864 satisfying tex2html_wrap_inline4866 with tex2html_wrap_inline4868 . The state vector is given by

  equation491

with tex2html_wrap_inline4870 . The vector of fluxes is

  equation496

and the corresponding sources tex2html_wrap_inline4872 are

  equation514

The local characteristic structure of the previous system of equations was presented in [21Jump To The Next Citation Point In The Article]. 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 (30Popup Equation) 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 [21Jump To The Next Citation Point In The Article] (see also [96Jump To The Next Citation Point In The Article]).

The range of applications considered so far in general relativity employing the above formulation of the hydrodynamic equations, Equation (30Popup Equation, 31Popup Equation, 32Popup Equation, 33Popup Equation), 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., [167, 8Jump To The Next Citation Point In The Article]). The first applications in general relativity were performed, in one spatial dimension, in [163Jump To The Next Citation Point In The Article], 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 [39Jump To The Next Citation Point In The Article]. These results are discussed in [161Jump To The Next Citation Point In The Article, 38Jump To The Next Citation Point In The Article]. 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 [129Jump To The Next Citation Point In The Article, 244Jump To The Next Citation Point In The Article, 67Jump To The Next Citation Point In The Article]. 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. [96Jump To The Next Citation Point In The Article] (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 [96Jump To The Next Citation Point In The Article] the spectral decomposition (eigenvalues and right-eigenvectors) of the general relativistic hydrodynamic equations, valid for general spatial metrics, was derived, extending earlier results of [21Jump To The Next Citation Point In The Article] for non-diagonal metrics. A complete set of left-eigenvectors was presented by Ibáñ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 approaches2 Equations of General Relativistic 2 Equations of General Relativistic

image Numerical Hydrodynamics in General Relativity
José A. Font
http://www.livingreviews.org/lrr-2003-4
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to livrev@aei-potsdam.mpg.de