Go to previous page Go up Go to next page

2.1 Spacelike (3+1) approaches

In the 3+1 (ADM) formulation [25158], spacetime is foliated into a set of nonintersecting spacelike hypersurfaces. There are two kinematic variables describing the evolution between these surfaces: the lapse function α, which describes the rate of advance of time along a timelike unit vector nμ normal to a surface, and the spacelike shift vector i β that describes the motion of coordinates within a surface.

The line element is written as

ds2 = − (α2 − β βi)dx0dx0 + 2β dxidx0 + γ dxidxj, (11 ) i i ij
where γij is the 3-metric induced on each spacelike slice.

2.1.1 1+1 Lagrangian formulation (May and White)

For historical reasons it is worth beginning the overview of the formulations with the pioneering numerical work in general relativistic hydrodynamics, namely the one-dimensional gravitational collapse code of May and White [247Jump To The Next Citation Point248Jump To The Next Citation Point]. Building on theoretical work by Misner and Sharp [263Jump To The Next Citation Point], May and White developed a time-dependent numerical code to solve the evolution equations describing adiabatic spherical collapse in general relativity. This code was based on a Lagrangian finite-difference scheme (see Section 4.1), in which the coordinates are co-moving with the fluid. Artificial viscosity terms were included in the equations to damp the spurious numerical oscillations caused by the presence of shock waves in the flow solution. May and White’s formulation became the starting point of a large number of numerical investigations in subsequent years and, hence, it is useful to describe its main features in some detail.

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

ds2 = − a2(m, t)dt2 + b2(m, t)dm2 + R2 (m, t)(dθ2 + sin2 θdφ2), (12 )
m being a radial (Lagrangian) coordinate, indicating the total rest mass enclosed inside the sphere 43 πR (m, t)3.

The co-moving character of the coordinates leads, for a perfect fluid, to a stress-energy tensor whose nonvanishing components are 1 2 3 0 T1 = T 2 = T3 = − p,T 0 = (1 + ɛ)ρ. In these coordinates the local conservation equation for the baryonic mass, Equation (2View Equation), can be easily integrated to yield the metric potential b:

--1---- b = 4πρR2 . (13 )

The gravitational field equations, Equation (10View Equation), and the equations of motion, Equation (1View Equation), reduce to the following quasi-linear system of partial differential equations (see also [263]):

∂p 1 ∂a ----+ -----ρh = 0, (14 ) ∂m a ∂m
( ) ∂ɛ- -∂- 1- ∂t + p∂t ρ = 0, (15 )
∂u ( ∂p Γ M ) --- = − a 4π----R2 --+ --2 + 4πpR , (16 ) ∂t ∂m h R
1 ∂ρR2 ∂u∕ ∂m ---2------= − a-------, (17 ) ρR ∂t ∂R ∕∂m
with the definitions u = 1∂R- a∂t and Γ = 1 ∂R- b∂m, satisfying Γ 2 = 1 − u2 − 2M R. Additionally,
∫ m 2 ∂R M = 4πR ρ (1 + ɛ)----dm, (18 ) 0 ∂m
represents the total mass interior to radius m at time t. The final system, Equations (14View Equation) – (17View Equation), is closed with an EOS of the form given by Equation (9View Equation).

Hydrodynamics codes based on the original formulation of May and White and on later versions (e.g., [406Jump To The Next Citation Point]) have been used in many nonlinear simulations of supernova and neutron-star collapse (see, e.g., [262391] and references therein), as well as in perturbative computations of spherically-symmetric gravitational collapse within the framework of the linearized Einstein equations [346347]. However, the Lagrangian character of May and White’s formulation, together with other theoretical considerations concerning the particular coordinate gauge, has prevented its extension to multiple-dimensional calculations. In spite of this, for one-dimensional problems, the Lagrangian approach adopted by May and White has considerable advantages with respect to an Eulerian approach with spatially-fixed coordinates, most notably the reduction of numerical diffusion.

2.1.2 3+1 Eulerian formulation (Wilson)

The use of Eulerian coordinates in multiple-dimensional numerical-relativistic hydrodynamics started with the pioneering work of Wilson [411Jump To The Next Citation Point417Jump To The Next Citation Point]. Introducing the basic dynamic variables D, S μ, and E, representing the relativistic density, momenta, and generalized internal energy, respectively, and defined as

0 0 0 D = ρu , Sμ = ρhu μu , E = ρɛu , (19 )
the equations of motion in Wilson’s formulation [411Jump To The Next Citation Point413Jump To The Next Citation Point] are:
1 ∂ √ --- 1 ∂ √ --- √---- --0(D − g) + √-------i(DV i − g) = 0, (20 ) − g ∂x − g∂x
1 ∂ √ --- 1 ∂ √ --- ∂p 1 ∂gαβ S S √-------0(Sμ − g) + √-------i(SμV i − g) +---μ + -----μ--α-0β = 0, (21 ) − g∂x − g∂x ∂x 2 ∂x S
∂ √ --- ∂ √ --- ∂ √ --- ---0(E − g) + ---i(EV i − g ) + p-μ(u0V μ − g) = 0, (22 ) ∂x ∂x ∂x
with the “transport velocity” given by V μ = u μ∕u0. We note that, in the original formulation [413Jump To The Next Citation Point], the momentum density equation, Equation (21View Equation), is solved only for the three spatial components, Si, and S0 is obtained through the 4-velocity normalization condition u u μ = − 1 μ.

A direct inspection of the system shows that the equations are written as a coupled set of advection equations. In doing so, the terms containing derivatives (in space or time) of the pressure are treated as source terms. This approach, hence, sidesteps an important guideline for the formulation of nonlinear hyperbolic systems of equations, namely the preservation of their conservation form. This is a necessary condition to guarantee correct evolution in regions of sharp entropy generation (i.e., shocks). Furthermore, some amount of numerical dissipation must be used to stabilize the solution across discontinuities. In this spirit, the first attempt to solve the equations of general-relativistic hydrodynamics in the original Wilson’s scheme [411Jump To The Next Citation Point] 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 [407Jump To The Next Citation Point] to the relativistic regime (see Section 4.1.1).

Wilson’s formulation has been (and still is) widely used in hydrodynamic codes developed by a variety of research groups. (The latest extensions made to incorporate magnetic fields are discussed in Section 3.1.1). Many different astrophysical scenarios were first investigated with these codes, including axisymmetric stellar core collapse [279Jump To The Next Citation Point277Jump To The Next Citation Point28338Jump To The Next Citation Point386Jump To The Next Citation Point319Jump To The Next Citation Point118Jump To The Next Citation Point], accretion onto compact objects [170Jump To The Next Citation Point317Jump To The Next Citation Point], numerical cosmology [7273Jump To The Next Citation Point17] and, more recently, the coalescence of neutron-star binaries [416Jump To The Next Citation Point418Jump To The Next Citation Point242Jump To The Next Citation Point] (see [417Jump To The Next Citation Point] for an up-to-date summary of such studies). This formalism has also been employed, in the special-relativistic limit, in numerical studies of heavy-ion collisions [415249]. We note that in most of these investigations, the original formulation of the hydrodynamic equations was slightly modified by redefining the dynamic variables, Equation (19View Equation), with the addition of a multiplicative α factor (the lapse function) and the introduction of the Lorentz factor, 0 W ≡ αu:

D = ρW, Sμ = ρhW u μ, E = ρɛW. (23 )

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

More recently, Wilson’s formulation has been applied to the numerical study of the coalescence of neutron-star binaries in general relativity [416Jump To The Next Citation Point418Jump To The Next Citation Point242Jump To The Next Citation Point] (see Section 5.3.2). These studies adopt an approximation scheme for the gravitational field, by imposing the simplifying condition that the three-geometry (the 3-metric γij) is conformally flat. The curvature of the three metric is then described by a position-dependent conformal factor φ4 times a flat-space Kronecker delta, and the line element, Equation (11View Equation), reads

2 2 i 0 0 i 0 4 i j ds = − (α − βiβ )dx dx + 2βidx dx + φ δijdx dx . (24 )

Therefore, in this approximation scheme all radiation degrees of freedom are removed. Moreover, under the maximal-slicing condition (K = 0), the field equations reduce to a set of five Poisson-like elliptic equations in flat spacetime for the lapse, the shift vector, and the conformal factor. While in spherical symmetry this approach is no longer an approximation, being identical to Einstein’s theory, beyond spherical symmetry its quality degrades. In [193] it was shown by means of numerical simulations of extremely relativistic disks of dust that it has the same accuracy as the first post-Newtonian approximation. In less extreme situations, however, as in rotational stellar-core collapse, the approximation yields results comparable to those obtained using full general relativity (see [74Jump To The Next Citation Point92Jump To The Next Citation Point365Jump To The Next Citation Point]).

Wilson’s formulation showed some limitations in handling situations involving ultrarelativistic flows (W ≥ 2), as first pointed out by Centrella and Wilson [73Jump To The Next Citation Point]. Norman and Winkler [291Jump To The Next Citation Point] performed a comprehensive numerical assessment of such formulation by means of special-relativistic hydrodynamics simulations. Figure 1View Image reproduces a plot from [291Jump To The Next Citation Point] 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 [73Jump To The Next Citation Point]. This figure shows that for Lorentz factors of about 2 (v ≈ 0.86 c), the threshold of the ultrarelativistic limit, the relative errors are between 5% and 7% (depending on the adiabatic exponent of the gas), showing a linear growth with W.

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 [73Jump To The Next Citation Point]. The plot shows the dependence of the relative errors of the density compression ratio versus the Lorentz factor W for two different values of the adiabatic index of the flow, Γ = 4∕3 (triangles) and Γ = 5 ∕3 (circles) gases. The relative error is measured with respect to the average value of the density over a region in the shocked material. The data are from Centrella and Wilson [73Jump To The Next Citation Point] and the plot reproduces a similar one from Norman and Winkler [291Jump To The Next Citation Point].

Norman and Winkler [291Jump To The Next Citation Point] conclude that those large errors were mainly due to the way in which the artificial viscosity terms are included in the numerical scheme in Wilson’s formulation. These terms, commonly called Q in the literature (see Section 4.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 (21View Equation), and at the divergence of the velocity in the source of the energy equation, Equation (22View Equation). However, [291Jump To The Next Citation Point] proposes that one add the Q terms in a consistent way, in order to consider the artificial viscosity as a real viscosity. Hence, the hydrodynamic equations should be rewritten for a modified stress-energy tensor of the following form:

μν μ ν μν T = ρ(1 + ɛ + (p + Q )∕ρ)u u + (p + Q )g . (25 )
In this way, for instance, the momentum equation takes the following form (in flat spacetime):
-∂-- 2 -∂-- 2 i ∂(p +-Q)- ∂x0 [(ρh + Q )W Vj ] + ∂xi [(ρh + Q )W VjV ] + ∂xj = 0. (26 )
In Wilson’s original formulation, Q is omitted in the two terms containing the quantity ρh. In general, Q is a nonlinear function of the velocity and, hence, the quantity QW 2V in the momentum density of Equation (26View 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 more strongly coupled than their Newtonian counterparts. As a result, Norman and Winkler proposed the use of implicit schemes as a way to describe more accurately such coupling. Their code, which in addition incorporates an adaptive grid, reproduces very accurate results even for ultrarelativistic flows with Lorentz factors of about ten in one-dimensional, flat spacetime simulations.

Recently, Anninos and Fragile [19Jump To The Next Citation Point] have compared state-of-the-art artificial viscosity schemes and high-order nonoscillatory central schemes (see Section 4.1.3) using Wilson’s formulation for the former class of schemes and a conservative formulation for the latter (similar to the one considered in [311Jump To The Next Citation Point309Jump To The Next Citation Point]; see Section 2.2.2). Using the three-dimensional Cartesian code cosmos, these authors found that earlier results for artificial viscosity schemes in shock tube tests or shock reflection tests are not improved, i.e. the numerical solution becomes increasingly unstable for shock velocities greater than ∼ 0.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 (see Figure 8 of [19Jump To The Next Citation Point]), the underlying reason being the use of a conservative formulation of the hydrodynamic equations rather than the particular scheme employed (see Section 4.1.3).

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

In 1991, Martí, Ibáñez, and Miralles [237Jump To The Next Citation Point] presented a new formulation of the (Eulerian) general-relativistic hydrodynamics equations. This formulation was designed to take fundamental advantage of the hyperbolic and conservative character of the equations, contrary to the formulation discussed in the previous Section 2.1.2. From the numerical point of view, the hyperbolic and conservative nature of the relativistic Euler equations allows for the use of schemes based on the characteristic fields of the system, bringing to relativistic hydrodynamics existing tools of computational classical fluid dynamics. This procedure departs from earlier approaches, most notably in avoiding the need for artificial dissipation terms to handle discontinuous solutions [411Jump To The Next Citation Point413Jump To The Next Citation Point], as well as implicit schemes as proposed in [291Jump To The Next Citation Point].

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., [218Jump To The Next Citation Point]). 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 [237Jump To The Next Citation Point].

Theoretical advances in the mathematical character of the relativistic hydrodynamic equations were first achieved studying the special relativistic limit. In Minkowski spacetime, the hyperbolic character of relativistic hydrodynamics and MHD is exhaustively studied by Anile and collaborators (see [15Jump To The Next Citation Point] and references therein) by applying Friedrichs’ definition of hyperbolicity [143] to a quasi-linear form of the system of hydrodynamic and MHD equations,

μ ∂w-- 𝒜 (w )∂xμ = 0, (27 )
where μ 𝒜 are the Jacobian matrices of the system and w is a suitable set of primitive (physical) variables (see below). The system (27View Equation) is hyperbolic in the time direction defined by the vector field ξ with ξμξ μ = − 1, if the following two conditions hold: (i) det (𝒜 μξμ) ⁄= 0 and (ii) for any ζ such that ζμξμ = 0, ζμζμ = 1, the eigenvalue problem 𝒜 μ(ζμ − λ ξμ)r = 0 has only real eigenvalues {λ } n n=1,⋅⋅⋅,5 and a complete set of right-eigenvectors {r } n n=1,⋅⋅⋅,5. Besides verifying the hyperbolic character of the relativistic hydrodynamic equations, Anile and collaborators [15Jump To The Next Citation Point] obtained the explicit expressions for the eigenvalues and eigenvectors in the local rest frame, characterized by uμ = δμ0. In Font et al. [134Jump To The Next Citation Point] those calculations were extended to an arbitrary reference frame in which the motion of the fluid is described by the 4-velocity u μ = W (1,vi).

The approach followed in [134] for the equations of special-relativistic hydrodynamics was extended to general relativity in [36Jump To The Next Citation Point]. The choice of state vector (conserved quantities) in the 3+1 Eulerian formulation developed by Banyuls et al. [36Jump To The Next Citation Point] differs slightly from that of Wilson’s formulation [411Jump To The Next Citation Point]. It comprises the rest-mass density (D), the momentum density in the j-direction (S j), and the total energy density (E), measured by a family of observers, which is the natural extension (for a generic spacetime) of the Eulerian observers in classical fluid dynamics. Interested readers are directed to [36Jump To The Next Citation Point] for more complete definitions and geometrical foundations.

In terms of the primitive variables w = (ρ,vi,ɛ), the conserved quantities are written as:

D = ρW, (28 )
Sj = ρhW 2vj, (29 )
E = ρhW 2 − p, (30 )
where the contravariant components i ij v = γ vj of the three-velocity are defined as
i ui βi v = ---0 + ---, (31 ) αu α
and W is the relativistic Lorentz factor 0 2 −1∕2 W ≡ αu = (1 − v ) with 2 i j v = γijv v.

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 resulting first-order flux-conservative hyperbolic system, well suited for numerical applications, reads:

( √ -- √ --- i ) √-1-- ∂--γU-(w-)+ ∂--−-gF-(w-)- = S(w ), (32 ) − g ∂x0 ∂xi
with g ≡ det(gμν) satisfying √--- √ -- − g = α γ with γ ≡ det(γij). The state vector is given by
U (w) = (D, Sj,τ ), (33 )
with τ ≡ E − D. The vector of fluxes is
( ( ) ( ) ( ) ) i i βi- i βi- i i βi- i F (w ) = D v − α ,Sj v − α + pδj,τ v − α + pv , (34 )
and the corresponding sources S(w ) are
( ( ) ( )) μν ∂gνj δ μ0∂lnα μν 0 S (w ) = 0,T ---μ-− Γνμgδj ,α T ---μ-− T Γ νμ . (35 ) ∂x ∂x

The local characteristic structure of the previous system of equations was presented in [36Jump To The Next Citation Point]. The eigenvalues (characteristic speeds) of the corresponding Jacobian matrices are all real (but not distinct, one showing a threefold degeneracy as a result of the assumed directional splitting approach) and a complete set of right eigenvectors exists. More precisely, for a fluid moving along the x-direction, the eigenvalues read:

λ = αvx − βx (triple), (36 ) 0 { ∘ ------------------------------------} λ± = ----α--- vx (1 − c2) ± cs (1 − v2)[γxx(1 − v2c2) − vxvx (1 − c2)] 1 − v2c2s s s s − βx, (37 )
where cs is the speed of sound. We note that the Minkowskian limit of these expressions is recovered properly (see [101Jump To The Next Citation Point]) as well as the Newtonian one (λ0 = vx, λ± = vx ± cs). Hence, system (32View Equation) satisfies the definition of hyperbolicity. As will become apparent in Section 4.1.2, the knowledge of the spectral information of the flux-vector Jacobian matrices is essential in order to construct (upwind) HRSC schemes based on Riemann solvers. This information can be found in [36Jump To The Next Citation Point] (see also [137Jump To The Next Citation Point]).

Three-dimensional, Eulerian codes, which evolve the coupled system of Einstein and hydrodynamics equations following the conservative Eulerian approach discussed in this section have been developed by Font et al. [137Jump To The Next Citation Point] and by Baiotti et al. [27Jump To The Next Citation Point] (see Section 4.4 for further details). In particular, in [137Jump To The Next Citation Point] the spectral decomposition (eigenvalues and right-eigenvectors) of the general-relativistic hydrodynamics equations, valid for general spatial metrics, was derived, extending earlier results of [36Jump To The Next Citation Point] for nondiagonal metrics. A complete set of left-eigenvectors was presented by Ibáñez et al. [176Jump To The Next Citation Point]. Due to the paramount importance of the characteristic structure of the equations in the design of upwind HRSC schemes based upon Riemann solvers, we summarize all necessary information in Section 6.2 of the article.

The range of applications considered so far in relativistic astrophysics employing the above formulation of the hydrodynamics equations, Equations (32View Equation) – (35View Equation), includes the study of gravitational collapse (both stellar-core collapse to neutron stars and black-hole formation), accretion flows on to black holes, as well as neutron-star pulsations and neutron-star–binary mergers (see Section 5). The first applications in general relativity were performed, in one spatial dimension, in [237Jump To The Next Citation Point], using a slightly different form of the equations. Preliminary investigations of gravitational stellar collapse were attempted in [23553] by coupling the above formulation of the hydrodynamic equations to a hyperbolic formulation of the Einstein equations developed by [54]. Evolutions of fully-dynamic spacetimes in the context of stellar-core collapse, both in spherical symmetry and in axisymmetry, have also been done [178Jump To The Next Citation Point339Jump To The Next Citation Point96Jump To The Next Citation Point97Jump To The Next Citation Point74Jump To The Next Citation Point]. These investigations are discussed in Section 5.1.1. In the special-relativistic limit this formulation has successfully been applied to simulate the evolution of (ultra-) relativistic extragalactic jets, using numerical models of increasing complexity (see, e.g., [2417Jump To The Next Citation Point]).

  Go to previous page Go up Go to next page