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 [14], 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_inline3615, which describes the rate of advance of time along a timelike unit vector tex2html_wrap_inline3617 normal to a surface, and the spacelike shift vector tex2html_wrap_inline3619 that describes the motion of coordinates within a surface.

The line element is written as

equation148

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

2.1.1 May and White 

The pioneering numerical work in general relativistic hydrodynamics dates back to the one-dimensional gravitational collapse code of May and White [132Jump To The Next Citation Point In The Article, 133Jump To The Next Citation Point In The Article]. Building on theoretical work by Misner and Sharp [142Jump To The Next Citation Point In The Article], May and White developed a numerical code to solve the evolution equations describing adiabatic spherical collapse in general relativity. This code was based on a Lagrangian finite difference scheme, 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

equation163

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

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

eqnarray165

In these coordinates the local conservation equation for the baryonic mass, Eq. (2Popup Equation), can be easily integrated:

equation171

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

equation178

equation186

equation194

equation205

with the definitions

equation213

satisfying

equation221

Additionally,

  equation227

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

Codes based on the original formulation of May and White and on later versions (e.g. [221Jump To The Next Citation Point In The Article]) have been used in many non-linear simulations of supernova and neutron star collapse (see, e.g., [141Jump To The Next Citation Point In The Article, 214Jump To The Next Citation Point In The Article] and references therein), as well as in perturbative computations of spherically symmetric gravitational collapse employing the linearized Einstein theory [191, 193, 192]. In Section  4.1.1 some of these simulations are reviewed in detail. The Lagrangian character of May and White's code, together with other theoretical considerations concerning the particular coordinate gauge, has prevented its extension to multidimensional 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, mainly the lack of numerical diffusion.

2.1.2 Wilson 

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

  equation243

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

  equation252

  equation266

  equation291

with the ``transport velocity'' given by tex2html_wrap_inline3637 . Notice that the momentum density equation, Eq. (26Popup Equation), is only solved for the three spatial components, tex2html_wrap_inline3639, and tex2html_wrap_inline3641 is obtained through the 4-velocity normalization condition tex2html_wrap_inline3643 .

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 non-linear hyperbolic systems of equations, namely the preservation of their conservation form . This is a necessary feature to guarantee correct evolution in regions of sharp entropy generation (i.e. shocks). As a consequence, some amount of numerical dissipation must be used to stabilize the solution across discontinuities. The first attempt to solve the equations of general relativistic hydrodynamics in the original Wilson scheme [226Jump To The Next Citation Point In The Article] used a combination of finite difference upwind techniques with artificial viscosity terms. Such terms extended the classic treatment of shocks introduced by von Neumann and Richtmyer [223Jump To The Next Citation Point In The Article] into 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 [149Jump To The Next Citation Point In The Article, 148Jump To The Next Citation Point In The Article, 152, 18Jump To The Next Citation Point In The Article, 211Jump To The Next Citation Point In The Article, 177Jump To The Next Citation Point In The Article, 65Jump To The Next Citation Point In The Article], accretion onto compact objects [94Jump To The Next Citation Point In The Article, 175Jump To The Next Citation Point In The Article], numerical cosmology [47Jump To The Next Citation Point In The Article, 48Jump To The Next Citation Point In The Article, 12Jump To The Next Citation Point In The Article] and, more recently, the coalescence and merger of neutron star binaries [230Jump To The Next Citation Point In The Article, 231Jump To The Next Citation Point In The Article]. This formalism has also been extensively employed, in the special relativistic limit, in numerical studies of heavy-ion collisions [229, 135]. We note that in these investigations, the original formulation of the hydrodynamic equations was slightly modified by re-defining the dynamical variables, Eq. (24Popup Equation), with the addition of a multiplicative tex2html_wrap_inline3615 factor and the introduction of the Lorentz factor, tex2html_wrap_inline3647 (the ``relativistic gamma''):

  equation322

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, this was first considered in [227Jump 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 [209]. The resulting code was axially symmetric and aimed to integrate the coupled set of equations in the context of stellar core collapse [67].

More recently, Wilson's formulation has also been applied to the numerical study of the coalescence of binary neutron stars in general relativity [230Jump To The Next Citation Point In The Article, 231Jump To The Next Citation Point In The Article] (see Section  4.3). An approximation scheme for the gravitational field has been adopted in these studies, by imposing the simplifying condition that the three-geometry (the three metric) is conformally flat. The line element then reads

equation332

The curvature of the three metric is then described by a position dependent conformal factor tex2html_wrap_inline3649 times a flat-space Kronecker delta. Therefore, in this approximation scheme all radiation degrees of freedom are thrown away, and 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 approximation is identical to Einstein's theory, in more general situations it has the same accuracy as the first post-Newtonian approximation [108].

  

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 [48Jump To The Next Citation Point In The Article]: 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_inline3491 (triangles) and tex2html_wrap_inline3493 (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 [48Jump To The Next Citation Point In The Article], and the plot reproduces a similar one from Norman and Winkler [159Jump To The Next Citation Point In The Article].

Wilson's formulation showed some limitations in handling situations involving ultrarelativistic flows, as first pointed out by Centrella and Wilson [48Jump To The Next Citation Point In The Article]. Norman and Winkler [159Jump To The Next Citation Point In The Article] performed a comprehensive numerical study of such formulation by means of special relativistic hydrodynamical simulations. Fig.  1 reproduces a plot from [159Jump To The Next Citation Point In The Article] in which the relative error of the density compression ratio in the relativistic shock reflection problem - the heating of a cold gas which impacts at relativistic speeds with a solid wall and bounces back - is displayed as a function of the Lorentz factor W of the incoming gas. The source of the data is [48Jump To The Next Citation Point In The Article]. This figure shows that for Lorentz factors of about 2 (tex2html_wrap_inline3659), 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 [159Jump 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, called collectively Q (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 and at the divergence of the velocity in the source of the energy equation. However, [159Jump To The Next Citation Point In The Article] proposed to add the Q terms globally 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:

  equation352

In this way, in flat spacetime, the momentum equation takes the form:

  equation361

In Wilson's formulation Q is omitted in the two terms containing the quantity tex2html_wrap_inline3669 . In general Q is a non-linear function of the velocity and, hence, the quantity tex2html_wrap_inline3673 in the momentum density of Eq. (31Popup Equation) is a highly non-linear 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 to describe more accurately such coupling. Their code, which incorporates an adaptive grid, reproduces very accurate results even for ultrarelativistic flows with Lorentz factors of about 10 in one-dimensional flat spacetime simulations.

2.1.3 Valencia 

In 1991 the Valencia group  [125Jump 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. Numerically, the hyperbolic and conservative nature of the relativistic Euler equations allows for the use of schemes based on the characteristic fields of the system, translating existing tools of classical fluid dynamics to relativistic hydrodynamics. This procedure departs from earlier approaches, most notably in avoiding the need for artificial dissipation terms to handle discontinuous solutions [227], as well as implicit schemes as proposed in [159Jump To The Next Citation Point In The Article].

A numerical scheme written in conservation form automatically guarantees the correct Rankine-Hugoniot (jump) conditions across discontinuities (the shock-capturing property). 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 (HRSC in the following) schemes from classical fluid dynamics into the realm of relativity [125Jump To The Next Citation Point In The Article].

Theoretical advances on the mathematical character of the relativistic hydrodynamic equations were achieved studying the special relativistic limit. In Minkowski spacetime, the hyperbolic character of relativistic (magneto-) hydrodynamics 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 [78] to a quasi-linear form of the system of hydrodynamic equations,

  equation381

where tex2html_wrap_inline3675 are the Jacobian matrices of the system and tex2html_wrap_inline3677 are a suitable set of primitive variables (see below). System (32Popup Equation) will be hyperbolic in the time-direction defined by the vector field tex2html_wrap_inline3679 with tex2html_wrap_inline3681, if the following two conditions hold: (i) tex2html_wrap_inline3683 and (ii) for any tex2html_wrap_inline3587 such that tex2html_wrap_inline3687, tex2html_wrap_inline3689, the eigenvalue problem tex2html_wrap_inline3691 has only real eigenvalues tex2html_wrap_inline3693, and a complete set of right-eigenvectors tex2html_wrap_inline3695 . 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_inline3697 . In Font et al. [72Jump To The Next Citation Point In The Article] those calculations were extended to an arbitrary reference frame in which the motion of the fluid was described by the 4-velocity tex2html_wrap_inline3699 .

The extension to general relativity of the approach, followed in [72Jump To The Next Citation Point In The Article] for special relativity, was accomplished in [17Jump To The Next Citation Point In The Article]. We will refer to the formulation of the general relativistic hydrodynamic equations presented in [17Jump To The Next Citation Point In The Article] as the Valencia formulation . The choice of evolved variables (conserved quantities) in this formulation differs slightly from Wilson's formulation. It comprises the rest-mass density (D), the momentum density in the j -direction (tex2html_wrap_inline3705), 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 addressed to [17Jump To The Next Citation Point In The Article] for their definition and geometrical foundations.

In terms of the primitive variables tex2html_wrap_inline3709, the conserved quantities are written as:

  equation426

  equation429

  equation432

where the contravariant components tex2html_wrap_inline3711 of the three-velocity are defined as

equation436

and W is the relativistic Lorentz factor tex2html_wrap_inline3715 with tex2html_wrap_inline3717 .

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, coming from the spacetime geometry, which do not contain derivatives of stress-energy tensor components. More precisely, the first-order flux-conservative hyperbolic system, well suited for numerical applications, reads:

  equation449

with tex2html_wrap_inline3719 satisfying tex2html_wrap_inline3721 with tex2html_wrap_inline3723 . The state vector is given by

equation469

with tex2html_wrap_inline3725 . The vector of fluxes is

equation473

and the corresponding sources tex2html_wrap_inline3727 are

equation490

The local characteristic structure of the previous system of equations was presented in [17Jump 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) and a complete set of right-eigenvectors exists. System (37Popup Equation) satisfies, hence, the definition of hyperbolicity. As discussed 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 [17Jump To The Next Citation Point In The Article] (see also [75Jump To The Next Citation Point In The Article]).

The range of applications considered so far in general relativity employing this formulation is still small and mostly devoted to the study of accretion flows onto black holes (see Section  4.2.2 below). In the special relativistic limit this formulation is being successfully applied to model the evolution of (ultra-) relativistic extragalactic jets (see, e.g., [129, 9]). The first numerical studies in general relativity were performed, in one spatial dimension, in [125Jump 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 Valencia formulation to a hyperbolic formulation of the Einstein equations developed by [34Jump To The Next Citation Point In The Article]. Some discussion of these results can be found in [123Jump To The Next Citation Point In The Article, 33Jump To The Next Citation Point In The Article]. More recently, successful evolutions of fully dynamical spacetimes in the context of adiabatic spherically symmetric stellar core-collapse have been achieved [101Jump To The Next Citation Point In The Article, 184Jump To The Next Citation Point In The Article]. We will come back to these issues in Section  4.1.1 below.

Recently, a three-dimensional, Eulerian, general relativistic hydrodynamical code, evolving the coupled system of the Einstein and hydrodynamic equations, has been developed [75Jump To The Next Citation Point In The Article]. The formulation of the hydrodynamic equations follows the Valencia approach. The code is constructed for a completely general spacetime metric based on a Cartesian coordinate system, with arbitrarily specifiable lapse and shift conditions. In [75Jump 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, correcting earlier results of [17Jump To The Next Citation Point In The Article] for non-diagonal metrics. A complete set of left-eigenvectors has been recently presented by Ibáñez et al. [99]. This information is summarized in Section  5.2 .

The formulation of the coupled set of equations and the numerical code reported in [75Jump To The Next Citation Point In The Article] were used for the construction of the milestone code ``GR3D'' for the NASA Neutron Star Grand Challenge project. For a description of the project see the website of the Washington University Gravity Group [1]. A public domain version of the code has recently been released to the community at the same website, the source and documentation of this code can be downloaded at [2].



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-2000-2
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to livrev@aei-potsdam.mpg.de