Numerical Hydrodynamics in Special Relativity

Abstract

This review is concerned with a discussion of numerical methods for the solution of the equations of special relativistic hydrodynamics (SRHD). Particular emphasis is put on a comprehensive review of the application of high-resolution shock-capturing methods in SRHD. Results obtained with different numerical SRHD methods are compared, and two astrophysical applications of SRHD flows are discussed. An evaluation of the various numerical methods is given and future developments are analyzed.

1 Introduction

1.1 Current fields of research

Relativity is a necessary ingredient for describing astrophysical phenomena involving compact objects. Among these phenomena are core collapse supernovae, X-ray binaries, pulsars, coalescing neutron stars, black hole formations, micro-quasars, active galactic nuclei, superluminal jets and gamma-ray bursts. When strong gravitational fields are encountered as, for example, in the case of coalescing neutron stars or near black holes, general relativistic effects must be considered. Also the significant gravitational wave signal produced by some of these phenomena can only be understood in the framework of the general theory of relativity. There are, however, astrophysical phenomena which involve flows at relativistic speeds but no strong gravitational fields, and thus at least certain aspects of these phenomena can be described within the framework of special relativity alone, disregarding general relativistic effects.

Another field of research, where special relativistic “flows” are encountered, are present-day heavy-ion collision experiments taking place in large particle accelerators. The heavy ions are accelerated to ultra-relativistic velocities very close to the speed of light (∼ 99.998% [166]) to study the equation of state for hot dense nuclear matter.

1.2 Overview of the numerical methods

The first attempt to solve the equations of relativistic hydrodynamics (RHD) was made by Wilson [188, 189] and collaborators [28, 75] using an Eulerian explicit finite difference code with monotonic transport. The code relies on artificial viscosity techniques [185, 154] to handle shock waves. It has been widely used to simulate flows encountered in cosmology, axisymmetric relativistic stellar collapse, accretion onto compact objects and, more recently, collisions of heavy ions. Almost all the codes for numerical both special (SRHD) and general (GRHD) relativistic hydrodynamics developed in the eighties [142, 167, 126, 125, 127, 51] were based on Wilson’s procedure. However, despite its popularity it turned out to be uΠ_able to describe extremely relativistic flows (Lorentz factors larger than 2; see, e.g., [28]) accurately.

In the mid eighties, Norman & Winkler [131] proposed a reformulation of t he difference equations of SRHD with an artificial viscosity consistent with the relativistic dynamics of non-perfect fluids. The strong coupling introduced in the equations by the presence of the viscous terms in the definition of relativistic momentum and total energy densities required an implicit treatment of the difference equations. Accurate results across strong relativistic shocks with large Lorentz factors were obtained in combination with adaptive mesh techniques. However, no multidimensional version of this code was developed.

Attempts to integrate the RHD equations avoiding the use of artificial viscosity were performed in the early nineties. Dubal [45] developed a 2D code for relativistic magneto-hydrodynamics based on an explicit second-order Lax-Wendroff scheme incorporating a flux corrected transport (FCT) algorithm [20]. Following a completely different approach Mann [102] proposed a multidimensional code for general relativistic hydrodynamics based on smoothed particle hydrodynamics (SPH) techniques [121], which he applied to relativistic spherical collapse [104]. When tested against 1D relativistic shock tubes all these codes performed similar to the code of Wilson. More recently, Dean et al. [39] have applied flux correcting algorithms for the SRHD equations in the context of heavy ion collisions. Recent developments in relativistic SPH methods [30, 164] are discussed in Section 4.2.

A major break-through in the simulation of ultra-relativistic flows was accomplished when high-resolution shock-capturing (HRSC) methods, specially designed to solve hyperbolic systems of conservations laws, were applied to solve the SRHD equations [107, 106, 49, 50]. This review is intended to provide a comprehensive discussion of different HRSC methods and of related methods used in SRHD. Numerical methods for special relativistic MHD flows (MHD stands for magneto hydrodynamics) are not included, because they are beyond the scope of this review. However, we may include such a discussion in a future update of this article.

1.3 Plan of the review

The review is organized as follows: Section 2 contains a derivation of the equations of special relativistic (perfect) fluid dynamics, as well as a discussion of their main properties. In Section 3 the most recent developments in numerical methods for SRHD are reviewed paying particular attention to high-resolution shock-capturing methods. Other developments in special relativistic numerical hydrodynamics are discussed in Section 4. Numerical results obtained with different methods as well as analytical solutions for several test problems are presented in Section 6. Two astrophysical applications of SRHD are discussed in Section 7. An evaluation of the various numerical methods is given in Section 8 together with an outlook for future developments. Finally, some additional technical information is presented in Section 9.

The reader is assumed to have basic knowledge in classical [92, 35] and relativistic fluid dynamics [171, 6], as well as in finite difference / volume methods for partial differential equations [152, 132]. A discussion of modern finite volume methods for hyperbolic systems of conservation laws can be found, e.g., in [96, 98, 93]. The theory of spectral methods for fluid dynamics is developed in [24], and smoothed particle hydrodynamics (SPH) is reviewed in [121].

2 Special Relativistic Hydrodynamics

The equations of special relativistic (perfect) fluid dynamics are derived, and their main properties discussed. The derivation of the SRHD equations in 2.1 is supplemented by 9.1, which discusses algorithms to compute primitive variables, a procedure crucial in state-of-the-art SRHD codes. 2.2 reflects on the SRHD equations as a hyperbolic system of conservation laws, and 2.3 discusses the solution of the special relativistic Riemann problem, which is the basis for most modern numerical methods. This last subsection is completed by 9.3, where a FORTAN programme called RIEMANN for computing the solution of a special relativistic Riemann problem is provided for download.

2.1 Equations

Using the Einstein summation convention, the equations describing the motion of a relativistic fluid are given by the five conservation laws

$${(\rho {u^\mu })_{;\mu }}) = 0,$$

((1))

$${T^{\mu \nu }}_{;\nu } = 0,$$

((2))

where (μ, v = 0,…, 3), and where _;μ denotes the covariant derivative with respect to coordinate x^μ. Furthermore, ρ is the proper rest-mass density of the fluid, u^μ its four-velocity, and T^μv is the stress-energy tensor, which for a perfect fluid can be written as

$${T^{\mu \nu }} = \rho h{u^\mu }{u^\nu } + p{g^{\mu \nu }}.$$

((3))

Here g^μv is the metric tensor, p the fluid pressure, and and h the specific enthalpy of the fluid defined by

$$h = 1 + \varepsilon + \frac{p}{\rho },$$

((4))

where ɛ is the specific internal energy. Note that we use natural units (i.e., the speed of light c =1) throughout this review.

In Minkowski spacetime and Cartesian coordinates (t, x¹, x², x³), the conservation equations (1, 2) can be written in vector form as

$$\frac{{\partial {\mathrm{u}}}}{{\partial t}} + \frac{{\partial {{\mathrm{F}}^i}({\mathrm{u}})}}{{\partial {x^i}}} = 0,$$

((5))

where i = 1, 2, 3. The state vector u is defined by

$${\mathrm{u}} = {(D,\:{S^1},\:{S^2},\:{S^3},\:\tau )^{\mathrm{T}}},$$

((6))

and the flux vectors Fⁱ are given by

$${{\mathrm{F}}^i} = {(D{v^i},\:{S^1}{v^i} + p{\delta ^{1i}},\:{S^2}{v^i} + p{\delta ^{2i}},\:{S^3}{v^i} + p{\delta ^{3i}},\:{S^i} - D{v^i})^{\mathrm{T}}}.$$

((7))

The five conserved quantities D, S¹, S², S³ and τ are the rest-mass density, the three components of the momentum density, and the energy density (measured relative to the rest mass energy density), respectively. They are all measured in the laboratory frame, and are related to quantities in the local rest frame of the fluid (primitive variables) through

$$D = \rho W,$$

((8))

$${S^i} = \rho h{W^2}{v^i}\;\;\;\;(i = 1,2,3),$$

((9))

$$\tau = \rho h{W^2} - p - D,$$

((10))

where vⁱ are the components of the three-velocity of the fluid

$${v^i} = \frac{{{u^i}}}{{{u^0}}},$$

((11))

and W is the Lorentz factor

$$W = {u^0} = \frac{1}{{\sqrt {1 - {v^i}{v_i}} }}.$$

((12))

The system of equations (5) with definitions (6, 8, 9, 10, 11, 12) is closed by means of an equation of state (EOS), which we shall assume to be given in the form

$$p = p(\rho ,\:\varepsilon ).$$

((13))

In the non-relativistic limit (i.e., v ≪ 1, h → 1) D, Sⁱ and τ approach their Newtonian counterparts ρ, ρvⁱ and ρE = ρɛ + ρv²/2, and equations of system (5) reduce to the classical ones. In the relativistic case the equations of (5) are strongly coupled via the Lorentz factor and the specific enthalpy, which gives rise to numerical complications (see Section 2.3).

In classical numerical hydrodynamics it is very easy to obtain vⁱ from the conserved quantities (i.e., ρ and ρvⁱ). In the relativistic case, however, the task to recover (ρ, vⁱ, p) from (D, Sⁱ, τ) is much more difficult. Moreover, as state-of-the-art SRHD codes are based on conservative schemes where the conserved quantities are advanced in time, it is necessary to compute the primitive variables from the conserved ones one (or even several) times per numerical cell and time step making this procedure a crucial ingredient of any algorithm (see Section 9.1).

2.2 SRHD as a hyperbolic system of conservation laws

An important property of system (5) is that it is hyperbolic for causal EOS [6]. For hyperbolic systems of conservation laws, the Jacobians ∂Fⁱ(u)/∂u have real eigenvalues and a complete set of eigenvectors (see Section 9.2). Information about the solution propagates at finite velocities given by the eigenvalues of the Jacobians. Hence, if the solution is known (in some spatial domain) at some given time, this fact can be used to advance the solution to some later time (initial value problem). However, in general, it is not possible to derive the exact solution for this problem. Instead one has to rely on numerical methods which provide an approximation to the solution. Moreover, these numerical methods must be able to handle discontinuous solutions, which are inherent to non-linear hyperbolic systems.

The simplest initial value problem with discontinuous data is called a Riemann problem, where the one dimensional initial state consists of two constant states separated by a discontinuity. The majority of modern numerical methods, the so-called Godunov-type methods, are based on exact or approximate solutions of Riemann problems. Because of its theoretical and numerical importance, we discuss the solution of the special relativistic Riemann problem in the next subsection.

2.3 Exact solution of the Riemann problem in SRHD

Let us first consider the one dimensional special relativistic flow of an ideal gas with an adiabatic exponent γ in the absence of a gravitational field. The Riemann problem then consists of computing the breakup of a discontinuity, which initially separates two arbitrary constant states L (left) and R (right) in the gas (see Fig. 1 with L ≡ 1 and R ≡ 5). For classical hydrodynamics the solution can be found, e.g., in [35]. In the case of SRHD, the Riemann problem has been considered by Martí & Müller [108], who derived an exact solution generalizing previous results for particular initial data [173].

The solution to this problem is self-similar, because it only depends on the two constant states defining the discontinuity v_L and v_R, where v = (p, ρ, v), and on the ratio (x - x₀)/(t - t₀), where x₀ and t₀ are the initial location of the discontinuity and the time of breakup, respectively. Both in relativistic and classical hydrodynamics the discontinuity decays into two elementary nonlinear waves (shocks or rarefactions) which move in opposite directions towards the initial left and right states. Between these waves two new constant states v_L* and v_R* (note that v_L* ≡ 3 and v_R* ≡ 4 in Fig. 1) appear, which are separated from each other through a contact discontinuity moving with the fluid. Across the contact discontinuity the density exhibits a jump, whereas pressure and velocity are continuous (see Fig. 1). As in the classical case, the self-similar character of the flow through rarefaction waves and the Rankine-Hugoniot conditions across shocks provide the relations to link the intermediate states v_S* (S =L, R) with the corresponding initial states v_S. They also allow one to express the fluid flow velocity in the intermediate states v_S* as a function of the pressure p_S* in these states. Finally, the steadiness of pressure and velocity across the contact discontinuity implies

$${v_{{{\mathrm{L}}_*}}}({p_*}) = {v_{{{\mathrm{R}}_*}}}({p_*}),$$

((14))

where p_* = p_L* = P_R*, which closes the system. The functions v_S*(p) are defined by

$${v_{{S_*}}}(p) = \left\{ {\begin{array}{*{20}{l}} {\;\;{{\mathcal R}^S}(p)}&{{\mathrm{if}}\:p \le {p_S},}\\ {\;\;{{\mathcal S}^S}(p)}&{{\mathrm{if}}\:p > {p_S},} \end{array}} \right.$$

((15))

where ${{\mathcal R}^S}(p)/{{\mathcal S}^S}(p)$ denotes the family of all states which can be connected through a rarefaction / shock with a given state v_S ahead of the wave.

Figure 1

Schematic solution of a Riemann problem in special relativistic hydrodynamics. The initial state at t = 0 (top figure) consists of two constant states (1) and (5) with p₁ > p₅, ρ₁ > ρ₅, and v₁ = v₂ = 0 separated by a diaphragm at x_D. The evolution of the flow pattern once the diaphragm is removed (middle figure) is illustrated in a spacetime diagram (bottom figure) with a shock wave (solid line) and a contact discontinuity (dashed line) moving to the right. The bundle of solid lines represents a rarefaction wave propagating to the left.

The fact that one Riemann invariant is constant through any rarefaction wave provides the relation needed to derive the function ${{\mathcal R}^S}$

$${{\mathcal R}^S}(p) = \frac{{(1 + {v_S}){A_ \pm }(p) - (1 - {v_S})}}{{(1 + {v_S}){A_ \pm }(p) + (1 - {v_S})}},$$

((16))

with

$${A_ \pm }(p) = {\left( {\frac{{\sqrt {\gamma - 1} - c(p)}}{{\sqrt {\gamma - 1} + c(p)}}\frac{{\sqrt {\gamma - 1} + {c_S}}}{{\sqrt {\gamma - 1} - {c_S}}}} \right)^{ \pm \frac{2}{{\sqrt {\gamma - 1} }}}},$$

((17))

the + / - sign of A± corresponding to S =L / S =R. In the above equation, c_S is the sound speed of the state v_S, and c(p) is given by

$$c(p) = {\left( {\frac{{\gamma (\gamma - 1)p}}{{(\gamma - 1){\rho _S}{{(p/{p_S})}^{1/\gamma }} + \gamma p}}} \right)^{1/2}}.$$

((18))

The family of all states ${{\mathcal S}^S}(p)$, which can be connected through a shock with a given state v_S ahead of the wave, is determined by the shock jump conditions. One obtains

$$\begin{array}{*{20}{l}} {{{\mathcal S}^S}(p)\;\;\,\, = \;\;\,\,({h_S}{W_S}{v_S} \pm \frac{{p - {p_S}}}{{j(p)\sqrt {1 - {V_ \pm }{{(p)}^2}} }})}\\ {\;\;\,\;\;\,\;\;\,\:\: = \;\;\,{{\left[ {{h_S}{W_S} + (p - {p_S})\left( {\frac{1}{{{\rho _S}{W_S}}} \pm \frac{{{v_S}}}{{j(p)\sqrt {1 - {V_ \pm }{{(p)}^2}} }}} \right)} \right]}^{ - 1}},} \end{array}$$

((19))

where the + / - sign corresponds to S =R / S =L. V_±(p) and j(p) denote the shock velocity and the modulus of the mass flux across the shock front, respectively. They are given by

$${V_ \pm }(p) = \frac{{\rho _S^2W_S^2{v_S} \pm j{{(p)}^2}\sqrt {1 + {{({\rho _S}/j(p))}^2}} }}{{\rho _S^2W_S^2 + j{{(p)}^2}}},$$

((20))

and

$$j(p) = \sqrt {\frac{{{p_S} - p}}{{\frac{{h_S^2 - h{{(p)}^2}}}{{{p_S} - p}} - \frac{{2{h_S}}}{{{\rho _S}}}}}} ,$$

((21))

where the enthalpy h(p) of the state behind the shock is the (unique) positive root of the quadratic equation

$$\left( {1 + \frac{{(\gamma - 1)({p_S} - p)}}{{\gamma p}}} \right){h^2} - \frac{{(\gamma - 1)({p_S} - p)}}{{\gamma p}}h + \frac{{{h_S}({p_S} - p)}}{{{\rho _S}}} - h_S^2 = 0,$$

((22))

which is obtained from the Taub adiabat (the relativistic version of the Hugoniot adiabat) for an ideal gas equation of state.

The functions v_L*(p) and v_R*(p) are displayed in Fig. 2 in a p-v diagram for a particular set of Riemann problems. Once p_* has been obtained, the remaining state quantities and the complete Riemann solution,

$${\mathrm{u}} = {\mathrm{u}}((x - {x_0})/(t - {t_0});{{\mathrm{u}}_{\mathrm{L}}},\:{{\mathrm{u}}_{\mathrm{R}}})),$$

((23))

can easily be derived.

Figure 2

Graphical solution in the p-v plane of the Riemann problems defined by the initial states (p_L = 10³, ρ_L = 1, v_L = 0.5), and (p ⁱ_R , ρ_R = 1, v_R = 0) with i = 1, 2, 3, 4, where p ¹_R = 10², p ²_R = 10, p ³_R = 1, and p ⁴_R = 10^-1, respectively. The adiabatic index of the fluid is 5/3 in all cases. Note the asymptotic behavior of the functions as they approach v = 1 (i.e., the speed of light).

In Section 9.3 we provide a FORTRAN program called RIEMANN, which allows one to compute the exact solution of an arbitrary special relativistic Riemann problem using the algorithm just described.

The treatment of multidimensional special relativistic flows is significantly more difficult than that of multidimensional Newtonian flows. In SRHD all components (normal and tangential) of the flow velocity are strongly coupled through the Lorentz factor, which complicates the solution of the Riemann problem severely. For shock waves, this coupling ‘only’ increases the number of algebraic jump conditions, which must be solved. However, for rarefactions it implies the solution of a system of ordinary differential equations [108].

3 High-Resolution Shock-Capturing Methods

The application of high-resolution shock-capturing (HRSC) methods caused a revolution in numerical SRHD. These methods satisfy in a quite natural way the basic properties required for any acceptable numerical method:

• (i) high order of accuracy,

• (ii) stable and sharp description of discontinuities, and

• (iii) convergence to the physically correct solution.

Moreover, HRSC methods are conservative, and because of their shock capturing property discontinuous solutions are treated both consistently and automatically whenever and wherever they appear in the flow.

As HRSC methods are written in conservation form, the time evolution of zone averaged state vectors is governed by some functions (the numerical fluxes) evaluated at zone interfaces. Numerical fluxes are mostly obtained by means of an exact or approximate Riemann solver. High resolution is usually achieved by using monotonic polynomials in order to interpolate the approximate solutions within numerical cells.

Solving Riemann problems exactly involves time-consuming computations, which are particularly costly in the case of multidimensional SRHD due to the coupling of the equations through the Lorentz factor (see Section 2.3). Therefore, as an alternative, the usage of approximate Riemann solvers has been proposed.

In this section we summarize how the numerical fluxes are computed in a number of methods for numerical SRHD. Methods based on exact Riemann solvers are discussed in Sections 3.1 and 3.2, while those based on approximate solvers are discussed in Sections 3.3, 3.4, 3.5, 3.6, and 3.7. Readers not familiar with HRSC methods are referred to Section 9.4, where the basic properties of these methods are described and an outline of the recent developments is given.

3.1 Relativistic PPM

Martí & Müller [109] have used the procedure discussed in Section 2.3 to construct an exact Riemann solver, which they then incorporated in an extension of the piecewise parabolic method (PPM) [33] for 1D SRHD. In their relativistic PPM method numerical fluxes are calculated according to

$${{\mathrm{\hat F}}^{{\mathrm{RPPM}}}} = {\mathrm{F}}({\mathrm{u}}(0;\;{{\mathrm{u}}_{\mathrm{L}}},\:\;{{\mathrm{u}}_{\mathrm{R}}})),$$

((24))

where u_L and u_R are approximations of the state vector at the left and right side of a zone interface obtained by a second-order accurate interpolation in space and time, and u(0; u_L, u_R) is the solution of the Riemann problem defined by the two interpolated states at the position of the initial discontinuity.

The PPM interpolation algorithm described in [33] gives monotonic conservative parabolic profiles of variables within a numerical zone. In the relativistic version of PPM, the original interpolation algorithm is applied to zone averaged values of the primitive variables v = (p, ρ, v), which are obtained from zone averaged values of the conserved quantities u. For each zone j, the quartic polynomial with zone-averaged values a_j-2, a_j-1, a_j, a_j+1, and a_j+2 (where a = ρ, p, v) is used to interpolate the structure inside the zone. In particular, the values of a at the left and right interface of the zone, a_L,j and a_R,j, are obtained this way. These reconstructed values are then modified such that the parabolic profile, which is uniquely determined by a_L,j, a_R,j, and a_j, is monotonic inside the zone.

Both, the non relativitic PPM scheme described in [33] and the relativistic approach of [109] follow the same procedure to compute the time-averaged fluxes at an interface j + 1/2 separating zones j and j + 1. They are computed from two spatially averaged states, ${{\mathrm{v}}_{j + \frac{1}{2},\,{\mathrm{L}}}}$ and ${{\mathrm{v}}_{j + \frac{1}{2},\,{\mathrm{R}}}}$ at the left and right side of the interface, respectively. These left and right states are constructed taking into account the characteristic information reaching the interface from both sides during the time step. The relativistic version of PPM uses the characteristic speeds and Riemann invariants of the equations of relativistic hydrodynamics in this procedure.

3.2 The relativistic Glimm method

Wen et al. [187] have extended Glimm’s random choice method [65] to 1D SRHD. They developed a first-order accurate hydrodynamic code combining Glimm’s method (using an exact Riemann solver) with standard finite difference schemes.

In the random choice method, given two adjacent states, u ⁿ_j and u ⁿ_j+1 , at time tⁿ, the value of the numerical solution at time t^n+1/2 and position x_j+1/2 is given by the exact solution u(x, t) of the Riemann problem evaluated at a randomly chosen point inside zone (j, j + 1), i.e.,

$${\mathrm{u}}_{j + \frac{1}{2}}^{n + \frac{1}{2}} = {\mathrm{u}}\left( {\frac{{(j + {\xi _n})\Delta x}}{{\left( {n + \frac{1}{2}} \right)\Delta t}};{\mathrm{u}}_j^n,\:{\mathrm{u}}_{j + 1}^n} \right),$$

((25))

where ξ_n is a random number in the interval [0,1].

Besides being conservative on average, the main advantages of Glimm’s method are that it produces both completely sharp shocks and contact discontinuities, and that it is free of diffusion and dispersion errors.

Chorin [29] applied Glimm’s method to the numerical solution of homogeneous hyperbolic conservation laws. Colella [31] proposed an accurate procedure of randomly sampling the solution of local Riemann problems and investigated the extension of Glimm’s method to two dimensions using operator splitting methods.

3.3 Two-shock approximation for relativistic hydrodynamics

This approximate Riemann solver is obtained from a relativistic extension of Colella’s method [31] for classical fluid dynamics, where it has been shown to handle shocks of arbitrary strength [31, 191]. In order to construct Riemann solutions in the two-shock approximation one analytically continues shock waves towards the rarefaction side (if present) of the zone interface instead of using an actual rarefaction wave solution. Thereby one gets rid of the coupling of the normal and tangential components of the flow velocity (see Section 2.3), and the remaining minor algebraic complications are the Rankine-Hugoniot conditions across oblique shocks. Balsara [8] has developed an approximate relativistic Riemann solver of this kind by solving the jump conditions in the shocks’ rest frames in the absence of transverse velocities, after appropriate Lorentz transformations. Dai & Woodward [36] have developed a similar Riemann solver based on the jump conditions across oblique shocks making the solver more efficient.

Table 1 gives the converged solution for the intermediate states obtained with both Balsara’s and Dai & Woodward’s procedure for the case of the Riemann problems defined in Section 6.2 (involving strong rarefaction waves) together with the exact solution. Despite the fact that both approximate methods involve very different algebraic expressions, their results differ by less than 2%. However, the discrepancies are much larger when compared with the exact solution (up to a 100% error in the density of the left intermediate state in Problem 2). The accuracy of the two-shock approximation should be tested in the ultra-relativistic limit, where the approximation can produce large errors in the Lorentz factor (in the case of Riemann problems involving strong rarefaction waves) with important implications for the fluid dynamics. Finally, the suitability of the two-shock approximation for Riemann problems involving transversal velocities still needs to be tested.

Table 1 Pressure p_*, velocity v_*, and densities ρ_L* (left), ρ_R* (right) for the intermediate state obtained for the two-shock approximation of Balsara [8] (B) and of Dai & Woodward [36] (DW) compared to the exact solution (Exact) for the Riemann problems defined in Section 6.2.

3.4 Roe-type relativistic solvers

Linearized Riemann solvers are based on the exact solution of Riemann problems of a modified system of conservation equations obtained by a suitable linearization of the original system. This idea was put forward by Roe [155], who developed a linearized Riemann solver for the equations of ideal (classical) gas dynamics. Eulderink at al. [49, 50] have extended Roe’s Riemann solver to the general relativistic system of equations in arbitrary spacetimes. Eulderink uses a local linearization of the Jacobian matrices of the system fulfilling the properties demanded by Roe in his original paper.

Let ${\mathcal B} = \partial {\mathrm{F}}/\partial {\mathrm{u}}$ be the Jacobian matrix associated with one of the fluxes F of the original system, and u the vector of unknowns. Then, the locally constant matrix ${\tilde {\mathcal B}}$, depending on u_L and u_R (the left and right state defining the local Riemann problem) must have the following four properties:

• 1. It constitutes a linear mapping from the vector space u to the vector space F.

• 2. As ${{\mathrm{u}}_{\mathrm{L}}} \to {{\mathrm{u}}_{\mathrm{R}}} \to {\mathrm{u}},\;\tilde {\mathcal B}({{\mathrm{u}}_{\mathrm{L}}},\:{{\mathrm{u}}_{\mathrm{R}}}) \to {\mathcal B}({\mathrm{u}}).$.

• 3. For any ${{\mathrm{u}}_{\mathrm{L}}}{\mathrm{,}}\;{{\mathrm{u}}_{\mathrm{R}}},\;\tilde {\mathcal B}({{\mathrm{u}}_{\mathrm{L}}},\:{{\mathrm{u}}_{\mathrm{R}}})({{\mathrm{u}}_{\mathrm{R}}} - \:{{\mathrm{u}}_{\mathrm{L}}}) = {\mathrm{F(}}{{\mathrm{u}}_{\mathrm{R}}}) - {\mathrm{F(}}{{\mathrm{u}}_{\mathrm{R}}}).$.

• 4. The eigenvectors of $\tilde {\mathcal B}$ are linearly independent.

Conditions 1 and 2 are necessary if one is to recover smoothly the linearized algorithm from the nonlinear version. Condition 3 (supposing 4 is fulfilled) ensures that if a single discontinuity is located at the interface, then the solution of the linearized problem is the exact solution of the nonlinear Riemann problem.

Once a matrix $\tilde {\mathcal B}$ satisfying Roe’s conditions has been obtained for every numerical interface, the numerical fluxes are computed by solving the locally linear system. Roe’s numerical flux is then given by

$${\widehat {\mathrm{F}}^{{\mathrm{ROE}}}} = \frac{1}{2}\left[ {{\mathrm{F}}({{\mathrm{u}}_{\mathrm{L}}}) + {\mathrm{F}}({{\mathrm{u}}_{\mathrm{R}}}) - \sum\limits_p | \;{{\tilde \lambda }^{(p)}}|{{\tilde \alpha }^{(p)}}{{{\mathrm{\tilde r}}}^{(p)}}} \right],$$

((26))

with

$${\tilde \alpha ^{{{(p)}_ = }}}{{\mathrm{\tilde l}}^{(p)}}\;\cdot\;({{\mathrm{u}}_{\mathrm{R}}} - {{\mathrm{u}}_{\mathrm{L}}}),$$

((27))

where ̃λ^(p), ̃r^(p), and ̃l^(p) are the eigenvalues and the right and left eigenvectors of $\tilde {\mathcal B}$, respectively (p runs from 1 to the total number of equations of the system).

Roe’s linearization for the relativistic system of equations in a general space-time can be expressed in terms of the average state [49, 50]

$${\mathrm{\tilde w}} = \frac{{{{\mathrm{w}}_{\mathrm{L}}} + {{\mathrm{w}}_{\mathrm{R}}}}}{{{k_{\mathrm{L}}} + {k_{\mathrm{R}}}}},$$

((28))

with

$${\mathrm{\tilde w}} = (k{u^0},\:k{u^1},\:k{u^2},\:k{u^3},\:kp/(\rho h)),$$

((29))

and

$${k^2} = \sqrt { - g} \rho h,$$

((30))

where g is the determinant of the metric tensor g_μv. The role played by the density ρ in case of the Cartesian non-relativistic Roe solver as a weight for averaging, is taken over in the relativistic variant by k, which apart from geometrical factors tends to ρ in the non-relativistic limit. A Riemann solver for special relativistic flows and the generalization of Roe’s solver to the Euler equations in arbitrary coordinate systems are easily deduced from Eulderink’s work. The results obtained in 1D test problems for ultra-relativistic flows (up to Lorentz factors 625) in the presence of strong discontinuities and large gravitational background fields demonstrate the excellent performance of the Eulderink-Roe solver [50].

Relaxing condition 3 above, Roe’s solver is no longer exact for shocks but still produces accurate solutions, and moreover, the remaining conditions are fulfilled by a large number of averages. The 1D general relativistic hydrodynamic code developed by Romero et al. [157] uses flux formula (26) with an arithmetic average of the primitive variables at both sides of the interface. It has successfully passed a long series of tests including the spherical version of the relativistic shock reflection (see Section 6.1).

Roe’s original idea has been exploited in the so-called local characteristic approach (see, e.g., [198]). This approach relies on a local linearization of the system of equations by defining at each point a set of characteristic variables, which obey a system of uncoupled scalar equations. This approach has proven to be very successful, because it allows for the extension to systems of scalar nonlinear methods. Based on the local characteristic approach are the methods developed by Marquina et al. [106] and Dolezal & Wong [42], which both use high-order reconstructions of the numerical characteristic fluxes, namely PHM [106] and ENO [42] (see Section 9.4).

3.5 Falle and Komissarov upwind scheme

Instead of starting from the conservative form of the hydrodynamic equations, one can use a primitive-variable formulation in quasi-linear form

$$\frac{{\partial {\mathrm{v}}}}{{\partial t}} + {\mathcal A}\frac{{\partial {\mathrm{v}}}}{{\partial x}} = 0,$$

((31))

where v is any set of primitive variables. A local linearization of the above system allows one to obtain the solution of the Riemann problem, and from this the numerical fluxes needed to advance a conserved version of the equations in time.

Falle & Komissarov [55] have considered two different algorithms to solve the local Riemann problems in SRHD by extending the methods devised in [53]. In a first algorithm, the intermediate states of the Riemann problem at both sides of the contact discontinuity, v_L* and v_R*, are obtained by solving the system

$${{\mathrm{v}}_{{{\mathrm{L}}_*}}} = {{\mathrm{v}}_{\mathrm{L}}} + {b_{\mathrm{L}}}{\mathrm{r}}_{\mathrm{L}}^ - ,\:\;\;\;\;{{\mathrm{v}}_{{{\mathrm{R}}_*}}} = {{\mathrm{v}}_{\mathrm{R}}} + {b_{\mathrm{R}}}{\mathrm{r}}_{\mathrm{R}}^ + ,$$

((32))

where r ^-_L is the right eigenvector of ${\mathcal A}({{\mathrm{v}}_{\mathrm{L}}})$ associated with sound waves moving upstream and r ⁺_R is the right eigenvector of ${\mathcal A}({{\mathrm{v}}_{\mathrm{R}}})$ of sound waves moving downstream. The continuity of pressure and of the normal component of the velocity across the contact discontinuity allows one to obtain the wave strengths b_L and b_R from the above expressions, and hence the linear approximation to the intermediate state v_*(v_L, v_R).

In the second algorithm proposed by Falle & Komissarov [55], a linearization of system (31) is obtained by constructing a constant matrix $\tilde {\mathcal A}({{\mathrm{v}}_{\mathrm{L}}},{{\mathrm{v}}_{\mathrm{R}}}) = {\mathcal A}\left( {\frac{1}{2}{\mathrm{(}}{{\mathrm{v}}_{\mathrm{L}}}{\mathrm{ + }}{{\mathrm{v}}_{\mathrm{R}}})} \right)$. The solution of the corresponding Riemann problem is that of a linear system with matrix $\tilde {\mathcal A}$, i.e.,

$${{\mathrm{v}}_*} = {{\mathrm{v}}_{\mathrm{L}}} + \sum\limits_{\tilde \lambda (p) < 0} {{{\tilde \alpha }^{(p)}}} {{\mathrm{\tilde r}}^{(p)}},$$

((33))

or, equivalently,

$${{\mathrm{v}}_*} = {{\mathrm{v}}_{\mathrm{R}}} - \sum\limits_{\tilde \lambda (p) > 0} {{{\tilde \alpha }^{(p)}}} {{\mathrm{\tilde r}}^{(p)}},$$

((34))

with

$${\tilde \alpha ^{(p)}} = {{\mathrm{\tilde l}}^{(p)}}\;\cdot\;({{\mathrm{v}}_{\mathrm{R}}} - {{\mathrm{v}}_{\mathrm{L}}}),$$

((35))

where ̃λ^(p), ̃r^(p), and ̃l^(p) are the eigenvalues and the right and left eigenvectors of $\tilde {\mathcal A}$, respectively (p runs from 1 to the total number of equations of the system).

In both algorithms, the final step involves the computation of the numerical fluxes for the conservation equations

$${{\mathrm{\hat F}}^{{\mathrm{FK}}}}{\mathrm{ = F(u(}}{{\mathrm{v}}_{\mathrm{*}}}{\mathrm{(}}{{\mathrm{v}}_{\mathrm{L}}}{\mathrm{, }}{{\mathrm{v}}_{\mathrm{R}}}{\mathrm{)))}}{\mathrm{.}}$$

((36))

3.6 Relativistic HLL method

Schneider et al. [161] have proposed to use the method of Harten, Lax & van Leer [74], HLL hereafter, to integrate the equations of SRHD. This method avoids the explicit calculation of the eigenvalues and eigenvectors of the Jaco-bian matrices and is based on an approximate solution of the original Riemann problems with a single intermediate state

$${{\mathrm{u}}^{{\mathrm{HLL}}}}(x/t;\;{{\mathrm{u}}_{\mathrm{L}}},\;{{\mathrm{u}}_{\mathrm{R}}}) = \left\{ {\begin{array}{*{20}{c}} {{{\mathrm{u}}_{\mathrm{L}}}}\\ {{{\mathrm{u}}_{\mathrm{*}}}}\\ {{{\mathrm{u}}_{\mathrm{R}}}} \end{array}} \right.\begin{array}{*{20}{c}} {{\mathrm{for }}\;x\; < \;{a_{\mathrm{L}}}t\;\;\;\;\;\;\;\;\;}\\ {\;\;{\mathrm{for }}\;{a_{\mathrm{L}}}t \le x \le \;{a_{\mathrm{R}}}t\;\;\;,}\\ {{\mathrm{for }}\;x\; > \;{a_{\mathrm{R}}}t\;\;\;\;\;\;\;\;\;} \end{array}$$

((37))

where a_L and a_R are lower and upper bounds for the smallest and largest signal velocities, respectively. The intermediate state u_* is determined by requiring consistency of the approximate Riemann solution with the integral form of the conservation laws in a grid zone. The resulting integral average of the Riemann solution between the slowest and fastest signals at some time is given by

$${{\mathrm{u}}_{\mathrm{*}}}{\mathrm{ = }}\frac{{{a_{\mathrm{R}}}{{\mathrm{u}}_{\mathrm{R}}}{\mathrm{ - }}{a_{\mathrm{L}}}{{\mathrm{u}}_{\mathrm{L}}}{\mathrm{ - F(}}{{\mathrm{u}}_{\mathrm{R}}}{\mathrm{) + F(}}{{\mathrm{u}}_{\mathrm{L}}}{\mathrm{)}}}}{{{a_{\mathrm{R}}}{\mathrm{ - }}{a_{\mathrm{L}}}}}{\mathrm{,}}$$

((38))

and the numerical flux by

$${{\mathrm{\hat F}}^{{\mathrm{HLL}}}} = \frac{{a_{\mathrm{R}}^ + {\mathrm{F(}}{{\mathrm{u}}_{\mathrm{L}}}) - a_{\mathrm{L}}^ - {\mathrm{F(}}{{\mathrm{u}}_{\mathrm{R}}}) + a_{\mathrm{R}}^ + a_{\mathrm{L}}^ - {\mathrm{(}}{{\mathrm{u}}_{\mathrm{R}}} - {{\mathrm{u}}_{\mathrm{L}}})}}{{a_{\mathrm{R}}^ + - a_{\mathrm{L}}^ - }},$$

((39))

where

$$a_{\mathrm{L}}^ - = {\mathrm{min\{ 0,}}\;{a_{\mathrm{L}}}\} ,\;\;\;\;a_{\mathrm{R}}^ + = {\mathrm{max\{ 0,}}\;{a_{\mathrm{R}}}\} .$$

((40))

An essential ingredient of the HLL scheme are good estimates for the smallest and largest signal velocities. In the non-relativistic case, Einfeldt [48] proposed to calculate them based on the smallest and largest eigenvalues of Roe’s matrix. This HLL scheme with Einfeldt’s recipe is a very robust upwind scheme for the Euler equations and possesses the property of being positively conservative. The method is exact for single shocks, but it is very dissipative, especially at contact discontinuities.

Schneider et al. [161] have presented results in 1D ultra-relativistic hydrodynamics using a version of the HLL method with signal velocities given by

$${a_{\mathrm{R}}} = (\bar v + {\bar c_s})/(1 + \bar v\,{\bar c_s}),$$

((41))

$${a_{\mathrm{L}}} = (\bar v - {\bar c_s})/(1 + \bar v\,{\bar c_s}),$$

((42))

where c_s is the relativistic sound speed, and where the bar denotes the arithmetic mean between the initial left and right states. Duncan & Hughes [46] have generalized this method to 2D SRHD and applied it to the simulation of relativistic extragalactic jets.

3.7 Marquina’s flux formula

Godunov-type schemes are indeed very robust in most situations although they fail spectacularly on occasions. Reports on approximate Riemann solver failures and their respective corrections (usually a judicious addition of artificial dissipation) are abundant in the literature [153]. Motivated by the search for a robust and accurate approximate Riemann solver that avoids these common failures, Donat & Marquina [44] have extended to systems a numerical flux formula which was first proposed by Shu & Osher [163] for scalar equations. In the scalar case and for characteristic wave speeds which do not change sign at the given numerical interface, Marquina’s flux formula is identical to Roe’s flux. Otherwise, the scheme switches to the more viscous, entropy satisfying local Lax-Friedrichs scheme [163]. In the case of systems, the combination of Roe and local-Lax-Friedrichs solvers is carried out in each characteristic field after the local linearization and decoupling of the system of equations [44]. However, contrary to Roe’s and other linearized methods, the extension of Marquina’s method to systems is not based on any averaged intermediate state.

Martí et al. have used this method in their simulations of relativistic jets [110, 111]. The resulting numerical code has been successfully used to describe ultra-relativistic flows in both one and two spatial dimensions with great accuracy (a large set of test calculations using Marquina’s Riemann solver can be found in Appendix II of [111]). Numerical experimentation in two dimensions confirms that the dissipation of the scheme is sufficient to eliminate the carbuncle phenomenon [153], which appears in high Mach number relativistic jet simulations when using other standard solvers [43].

Aloy et al. [3] have implemented Marquina’s flux formula in their three dimensional relativistic hydrodynamic code GENESIS.

Font et al. [59] have developed a 3D general relativistic hydro code where the matter equations are integrated in conservation form and fluxes are calculated with Marquina’s formula.

3.8 Symmetric TVD schemes with nonlinear numerical dissipation

The methods discussed in the previous subsections are all based on exact or approximate solutions of Riemann problems at cell interfaces in order to stabilize the discretization scheme across strong shocks. Another successful approach relies on the addition of nonlinear dissipation terms to standard finite difference methods. The algorithm of Davis [38] is based on such an approach. It can be interpreted as a Lax- Wendroff scheme with a conservative TVD (total variation diminishing) dissipation term. The numerical dissipation term is local, free of problem dependent parameters and does not require any characteristic information. This last fact makes the algorithm extremely simple when applied to any hyperbolic system of conservation laws.

A relativistic version of Davis’ method has been used by Koide et al. [82, 81, 129] in 2D and 3D simulations of relativistic magneto-hydrodynamic jets with moderate Lorentz factors. Although the results obtained are encouraging, the coarse grid zoning used in these simulations and the relative smallness of the beam flow Lorentz factor (4.56, beam speed ≈ 0.98c) does not allow for a comparison with Riemann-solver-based HRSC methods in the ultra-relativistic limit.

4 Other Developments

In this Section we summarize some recent developments in numerical RHD based on non-HRSC methods. The corresponding methods have been shown to be capable of simulating high Lorentz factor flows with shock waves. Van Put-ten’s approach, described in 4.1, was originally developed for numerical RMHD. 4.2 is devoted to outline recent relativistic extensions of SPH methods (originally developed for Newtonian hydrodynamics). Finally, 4.3 describes the main properties of the relativistic version of the beam scheme, a method based on the numerical solution of the equilibrium limit of the non-relativistic Boltzmann equation.

4.1 Van Putten’s approach

Relying on a formulation of Maxwell’s equations as a hyperbolic system in divergence form, van Putten [179] has devised a numerical method to solve the equations of relativistic ideal MHD in flat spacetime [181]. Here we only discuss the basic principles of the method in one spatial dimension. In van Putten’s approach, the state vector u and the fluxes F of the conservation laws are decomposed into a spatially constant mean (subscript 0) and a spatially dependent variational (subscript 1) part

$${\mathrm{u}}(t,\:x) = {{\mathrm{u}}_0}(t) + {{\mathrm{u}}_1}(t,\:x)\:,\:{\mathrm{F}}(t,\:x) = {{\mathrm{F}}_0}(t) + {{\mathrm{F}}_1}(t,\:x).$$

((43))

The RMHD (for relativistic MHD) equations then become a system of evolution equations for the integrated variational parts u_i^*, which reads

$$\frac{{\partial {{\mathrm{u}}_1}*}}{{\partial t}} + {{\mathrm{F}}_1} = 0,$$

((44))

together with the conservation condition

$$\frac{{d{{\mathrm{F}}_0}}}{{dt}} = 0.$$

((45))

The quantities u₁^* are defined as

$${{\mathrm{u}}_1}*(t,\:x) = \int_{}^x {{{\mathrm{u}}_{\mathrm{1}}}} (t,\:y)dy.$$

((46))

They are continuous, and standard methods can be used to integrate the system (44). Van Putten uses a leapfrog method.

The new state vector u(t, x) is then obtained from u₁^*(t, x) by numerical differentiation. This process can lead to oscillations in the case of strong shocks and a smoothing algorithm should be applied. Details of this smoothing algorithm and of the numerical method in one and two spatial dimensions can be found in [180] together with results on a large variety of tests.

Van Putten has applied his method to simulate relativistic hydrodynamic and magneto hydrodynamic jets with moderate flow Lorentz factors (< 4.25) [182, 184].

4.2 Relativistic SPH

Besides finite volume schemes, another completely different method is widely used in astrophysics for integrating the hydrodynamic equations. This method is Smoothed Particle Hydrodynamics, or SPH for short [100, 63, 121]. The fundamental idea of SPH is to represent a fluid by a Monte Carlo sampling of its mass elements. The motion and thermodynamics of these mass elements is then followed as they move under the influence of the hydrodynamics equations. Because of its Lagrangian nature there is no need within SPH for explicit integration of the continuity equation, but in some implementations of SPH this is done nevertheless for certain reasons. As both the equation of motion of the fluid and the energy equation involve continuous properties of the fluid and their derivatives, it is necessary to estimate these quantities from the positions, velocities and internal energies of the fluid elements, which can be thought of as particles moving with the flow. This is done by treating the particle positions as a finite set of interpolating points where the continuous fluid variables and their gradients are estimated by an appropriately weighted average over neighboring particles. Hence, SPH is a free-Lagrange method, i.e., spatial gradients are evaluated without the use of a computational grid.

A comprehensive discussion of SPH can be found in the reviews of Hernquist & Katz [76], Benz [12] and Monaghan [120, 121]. The non-relativistic SPH equations are briefly discussed in Section 9.5. The capabilities and limits of SPH are explored, e.g., in [169, 172], and the stability of the SPH algorithm is investigated in [170].

The SPH equations for special relativistic flows have been first formulated by Monaghan [120]. For such flows the SPH equations given in Section 9.5 can be taken over except that each SPH particle a carries v_a baryons instead of mass m_a [120, 30]. Hence, the rest mass of particle a is given by m_a = m₀v_a, where m₀ is the baryon rest mass (if the fluid is made of baryons). Transforming the notation used in [30] to ours, the continuity equation, the momentum and the total energy equations for particle a are given by (unit of velocity is c)

$$\frac{{d{N_a}}}{{dt}} = - \sum\limits_b {{\nu _b}} ({{\rm{v}}_a} - {{\rm{v}}_b})\:\cdot\:{\nabla _a}{W_{ab}},$$

((47))

$$\frac{{d{{\widehat {\rm{S}}}_a}}}{{dt}} = - \sum\limits_b {{\nu _b}} \left( {\frac{{{p_a}}}{{N_a^2}} + \frac{{{p_b}}}{{N_b^2}} + {\Pi _{ab}}} \right)\:\cdot\:{\nabla _a}{W_{ab}},$$

((48))

and

$$\frac{{d\hat \tau }}{{dt}} = - \sum\limits_b {{\nu _b}} \left( {\frac{{{p_a}{{\rm{v}}_a}}}{{N_a^2}} + \frac{{{p_b}{{\rm{v}}_b}}}{{N_b^2}} + {\Omega _{ab}}} \right)\:\cdot\:{\nabla _a}{W_{ab}},$$

((49))

respectively. Here, the summation is over all particles other than particle a, and d/dt denotes the Lagrangian time derivative.

$$N = \frac{D}{{{m_0}}}$$

((50))

is the baryon number density,

$${\mathrm{\hat S}} \equiv \frac{{\mathrm{S}}}{N} = {m_0}hW{\mathrm{v}}$$

((51))

the momentum per particle, and

$$\hat \tau \equiv \frac{\tau }{N} + {m_0} = {m_0}hW - \frac{p}{N}$$

((52))

the total energy per particle (all measured in the laboratory frame). The momentum density S = (S¹, S², S³)^T, the energy density τ (measured in units of the rest mass energy density), and the specific enthalpy h are defined in Section 2.1. Π_ab and Ω_ab are the SPH dissipation terms, and ∇_aW_ab denotes the gradient of the kernel W_ab (see Section 9.5 for more details).

Special relativistic flow problems have been simulated with SPH by [90, 80, 102, 104, 30, 164]. Extensions of SPH capable of treating general relativistic flows have been considered by [80, 89, 164]. Concerning relativistic SPH codes the artificial viscosity is the most critical issue. It is required to handle shock waves properly, and ideally it should be predicted by a relativistic kinetic theory for the fluid. However, unlike its Newtonian analogue, the relativistic theory has not yet been developed to the degree required to achieve this. For Newtonian SPH Lattanzio et al. [94] have shown that in high Mach number flows a viscosity quadratic in the velocity divergence is necessary. They proposed a form of the artificial viscosity such that the viscous pressure could be simply added to the fluid pressure in the equation of motion and the energy equation. Because this simple form of the artificial viscosity has known limitations, they also proposed a more sophisticated form of the artificial viscosity terms, which leads to a modified equation of motion. This artificial viscosity works much better, but it cannot be generalized to the relativistic case in a consistent way. Utilizing an equation for the specific internal energy both Mann [102] and Laguna et al. [89] use such an inconsistent formulation. Their artificial viscosity term is not included into the expression of the specific relativistic enthalpy. In a second approach, Mann [102] allows for a time-dependent smoothing length and SPH particle mass, and further proposed a SPH variant based on the total energy equation. Lahy [90] and Siegler & Riffert [164] use a consistent artificial viscosity pressure added to the fluid pressure. Siegler & Riffert [164] have also formulated the hydrodynamic equations in conservation form.

Monaghan [122] incorporates concepts from Riemann solvers into SPH. For this reason he also proposes to use a total energy equation in SPH simulations instead of the commonly used internal energy equation, which would involve time derivatives of the Lorentz factor in the relativistic case. Chow & Monaghan [30] have extended this concept and have proposed an SPH algorithm, which gives good results when simulating an ultra-relativistic gas. In both cases the intention was not to introduce Riemann solvers into the SPH algorithm, but to use them as a guide to improve the artificial viscosity required in SPH.

In Roe’s Riemann solver [155], as well as in its relativistic variant proposed by Eulerdink [49, 50] (see Section 3.4), the numerical flux is computed by solving a locally linear system and depends on both the eigenvalues and (left and right) eigenvectors of the Jacobian matrix associated to the fluxes and on the jumps in the conserved physical variables (see Eqs. (26) and (27)). Monaghan [122] realized that an appropriate form of the dissipative terms Π_ab and Ω_ab for the interaction between particles a and b can be obtained by treating the particles as the equivalent of left and right states taken with reference to the line joining the particles. The quantity corresponding to the eigenvalues (wave propagation speeds) is an appropriate signal velocity v_sig (see below), and that equivalent to the jump across characteristics is a jump in the relevant physical variable. For the artificial viscosity tensor, Π_ab, Monaghan [122] assumes that the jump in velocity across characteristics can be replaced by the velocity difference between a and b along the line joining them.

With these considerations in mind Chow & Monaghan [30] proposed for Π_ab in the relativistic case the form

$${\Pi _{ab}} = - \frac{{K{v_{sig}}({\mathrm{\hat S}}_a^* - {\mathrm{\hat S}}_b^*)\;\cdot\;{\mathrm{j}}}}{{{{\bar N}_{ab}}}},$$

((53))

when particles a and b are approaching, and Π_ab = 0 otherwise. Here K = 0.5 is a dimensionless parameter, which is chosen to have the same value as in the non-relativistic case [122]. ̅N_ab = (N_a + N_b)/2 is the average baryon number density, which has to be present in (53), because the pressure terms in the summation of (90) have an extra density in the denominator arising from the SPH interpolation. Furthermore,

$${\mathrm{j}} = \frac{{{{\mathrm{r}}_{ab}}}}{{|{{\mathrm{r}}_{ab}}|}}$$

((54))

is the unit vector from b to a, and

$${{\mathrm{\hat S}}^{\mathrm{*}}} = {m_0}h{W^*}{\mathrm{v}},$$

((55))

where

$${W^*} = \frac{1}{{\sqrt {1 - {{({\mathrm{v}} \cdot {\mathrm{j}})}^2}} }}.$$

((56))

Using instead of ̂S (see Eq. (51)) the modified momentum ̂S^*, which involves the line of sight velocity v · j, guarantees that the viscous dissipation is positive definite [30].

The dissipation term in the energy equation is derived in a similar way and is given by [30]

$${\Omega _{ab}} = - \frac{{K{v_{sig}}(\widehat {\tau _a^*} - \widehat {\tau _b^*}){\mathrm{j}}}}{{{{\bar N}_{ab}}}},$$

((57))

if a and b are approaching, and Ω_ab = 0 otherwise. Ω_ab involves the energy $\widehat {{\tau ^*}}$, which is identical to ̂τ (see Eq. (52)) except that W is replaced by W*.

To determine the signal velocity Chow & Monaghan [30] (and Monaghan [122] in the non-relativistic case) start from the (local) eigenvalues, and hence the wave velocities (v ± c_s)/(1 ± vc_s) and v of one-dimensional relativistic hydro-dynamic flows. Again considering particles a and b as the left and right states of a Riemann problem with respect to motions along the line joining the particles, the appropriate signal velocity is the speed of approach (as seen in the computing frame) of the signal sent from a towards b and that from b to a. This is the natural speed for the sharing of physical quantities, because when information about the two states meets it is time to construct a new state. This speed of approach should be used when determining the size of the time step by the Courant condition (for further details see [30]).

Chow & Monaghan [30] have demonstrated the performance of their Riemann problem guided relativistic SPH algorithm by calculating several shock tube problems involving ultra-relativistic speeds up to v = 0.9999. The algorithm gives good results, but finite volume schemes based on Riemann solvers give more accurate results and can handle even larger speeds (see Section 6).

4.3 Relativistic beam scheme

Sanders & Prendergast [159] proposed an explicit scheme to solve the equilibrium limit of the non-relativistic Boltzmann equation, i.e., the Euler equations of Newtonian fluid dynamics. In their so-called beam scheme the Maxwellian velocity distribution function is approximated by several Dirac delta functions or discrete beams of particles in each computational cell, which reproduce the appropriate moments of the distribution function. The beams transport mass, momentum, and energy into adjacent cells, and their motion is followed to first-order accuracy. The new (i.e., time advanced) macroscopic moments of the distribution function are used to determine the new local non-relativistic Maxwell distribution in each cell. The entire process is then repeated for the next time step. The Courant-Friedrichs-Levy (CFL) stability condition requires that no beam of gas travels farther than one cell in one time step. This beam scheme, although being a particle method derived from a microscopic kinetic description, has all the desirable properties of modern characteristic-based wave propagating methods based on a macroscopic continuum description.

The non-relativistic scheme of Sanders & Prendergast [159] has been extended to relativistic flows by Yang et al. [194]. They replaced the Maxwellian distribution function by its relativistic analogue, i.e., by the more complex Jüttner distribution function, which involves modified Bessel functions. For three-dimensional flows the Jüttner distribution function is approximated by seven delta functions or discrete beams of particles, which can be viewed as dividing the particles in each cell into seven distinct groups. In the local rest frame of the cell these seven groups represent particles at rest and particles moving in ±x, ±y and ±z directions, respectively.

Yang et al. [194] show that the integration scheme for the beams can be cast in the form of an upwind conservation scheme in terms of numerical fluxes. They further show that the beam scheme not only splits the state vector but also the flux vectors, and has some entropy-satisfying mechanism embedded as compared with approximate relativistic Riemann solver [42, 161] based on Roe’s method [155]. The simplest relativistic beam scheme is only first-order accurate in space, but can be extended to higher-order accuracy in a straightforward manner. Yang et al. consider three high-order accurate variants (TVD2, ENO2, ENO3) generalizing their approach developed in [195, 196] for Newtonian gas dynamics, which is based on the essentially non-oscillatory (ENO) piecewise polynomial reconstruction scheme of Harten at al. [73].

Yang et al. [194] present several numerical experiments including relativistic one-dimensional shock tube flows and the simulation of relativistic two-dimensional Kelvin-Helmholtz instabilities. The shock tube experiments consist of a mildly relativistic shock tube, relativistic shock heating of a cold flow, the relativistic blast wave interaction of Woodward & Colella [191] (see Section 6.2.3), and the perturbed relativistic shock tube flow of Shu & Osher [163].

5 Summary of Methods

This section contains a summary of all the methods reviewed in the two preceding sections as well as several FCT and artificial viscosity codes. The main characteristic of the codes (dissipation algorithm, spatial and temporal orders of accuracy, reconstruction techniques) are listed in two tables (Table 2 for HRSC codes; Table 3 for other approaches).

Table 2 High-resolution shock-capturing methods. All the codes rely on a conservation form of the RHD equations with the exception of ref. [187].

Table 3 Code characteristics.

6 Test Bench

In this section we compare the performance of the numerical methods described in the previous sections based on a couple of test problems which have an analytical solution. In 6.1 we compare how the different methods handle the relativistic shock heating of a cold gas in different geometries based on previously published data. In Table 4 we summarize the results and give for every numerical method both the highest Lorentz factor achieved for this problem and the mean error in the computation of the post-shock density. The results obtained with different numerical methods for two Riemann problems involving shock waves and relativistic velocities appear in Section 6.2.1 (mildly relativistic Riemann problem) and Section 6.2.2 (highly relativistic Riemann problem), respectively. The performance of the methods is summarized in Tables 6 and 7. Finally, a challenging test problem based on the collision of two relativistic blast waves is discussed in Section 6.2.3.

Table 4 Summary of relativistic shock heating test calculations by various authors in planar (α = 0), cylindrical (α = 1), and spherical (α = 2) geometry. W_max and (Terror are the maximum inflow Lorentz factor and compression ratio error extracted from tables and figures of the corresponding reference. W_max should only be considered as indicative of the maximum Lorentz factor achievable by every method. The methods are described in Sections 3 and 4 and their basic properties summarized in Section 5 (Tables 2, 3).

6.1 Relativistic shock heating in planar, cylindrical, and spherical geometry

Shock heating of a cold fluid in planar, cylindrical or spherical geometry has been used since the early developments of numerical relativistic hydrodynamics as a test case for hydrodynamic codes, because it has an analytical solution ([18] in planar symmetry; [111] in cylindrical and spherical symmetry), and because it involves the propagation of a strong relativistic shock wave.

In planar geometry, an initially homogeneous, cold (i.e., ε ≈ 0) gas with coordinate velocity v₁ and Lorentz factor W₁ is supposed to hit a wall, while in the case of cylindrical and spherical geometry the gas flow converges towards the axis or the center of symmetry. In all three cases the reflection causes compression and heating of the gas as kinetic energy is converted into internal energy. This occurs in a shock wave, which propagates upstream. Behind the shock the gas is at rest (v₂ = 0). Due to conservation of energy across the shock the gas has a specific internal energy given by

$${\varepsilon _2} = {W_1} - 1.$$

((58))

The compression ratio of shocked and unshocked gas, σ, follows from

$$\sigma = \frac{{\gamma + 1}}{{\gamma - 1}} + \frac{\gamma }{{\gamma - 1}}{\varepsilon _2},$$

((59))

where γ is the adiabatic index of the equation of state. The shock velocity is given by

$${V_s} = \frac{{(\gamma - 1){W_1}|{v_1}|}}{{{W_1} + 1}}.$$

((60))

In the unshocked region (r ∈ [V_st, ∞[) the pressure-less gas flow is self-similar and has a density distribution given by

$$\rho (t,\:r) = {\left( {1 + \frac{{|{v_1}|t}}{r}} \right)^\alpha }{\rho _0},$$

((61))

where α = 0, 1, 2 for planar, cylindrical or spherical geometry, and where ρ₀ is the density of the inflowing gas at infinity (see Fig. 3).

Figure 3

Schematic solution of the shock heating problem in spherical geometry. The initial state consists of a spherically symmetric flow of cold (p = 0) gas of unit rest mass density having a coordinate inflow velocity |v₁| = 1 everywhere. A shock is generated at the center of the sphere, which propagates upstream with constant speed. The post-shock state is constant and at rest. The pre-shock state, where the flow is self-similar, has a density which varies as ρ = (1 + t/r)² with time t and radius r.

In the Newtonian case the compression ratio σ of shocked and unshocked gas cannot exceed a value of σ_max = (γ + 1)/(γ - 1) independently of the inflow velocity. This is different for relativistic flows, where σ grows linearly with the flow Lorentz factor and becomes infinite as the inflowing gas velocity approaches to speed of light.

The maximum flow Lorentz factor achievable for a hydrodynamic code with acceptable errors in the compression ratio σ is a measure of the code’s quality. Table 4 contains a summary of the results obtained for the shock heating test by various authors.

Explicit finite-difference techniques based on a non-conservative formulation of the hydrodynamic equations and on non-consistent artificial viscosity [28, 75] are able to handle flow Lorentz factors up to ≈ 10 with moderately large errors (σ_error ≈ 1 - 3%) at best [190, 113]. Norman & Winkler [131] got very good results (σ_error ≈ 0.01%) for a flow Lorentz factor of 10 using consistent artificial viscosity terms and an implicit adaptive-mesh method.

The performance of explicit codes improved significantly when numerical methods based on Riemann solvers were introduced [107, 106, 49, 161, 50, 109, 55]. For some of these codes the maximum flow Lorentz factor is only limited by the precision by which numbers are represented on the computer used for the simulation [42, 187, 3].

Schneider et al. [161] have compared the accuracy of a code based on the relativistic HLL Riemann solver with different versions of relativistic FCT codes for inflow Lorentz factors in the range 1.6 to 50. They found that the error in σ was reduced by a factor of two when using HLL.

Within SPH methods, Chow & Monaghan [30] have obtained results comparable to those of HRSC methods (σ_error < 2 ¿ 10^-3) for flow Lorentz factors up to 70, using a relativistic SPH code with Riemann solver guided dissipation. Sieglert & Riffert [164] have succeeded in reproducing the post-shock state accurately for inflow Lorentz factors of 1000 with a code based on a consistent formulation of artificial viscosity. However, the dissipation introduced by SPH methods at the shock transition is very large (10 – 12 particles in the code of ref. [164]; 20 – 24 in the code of ref. [30]) compared with the typical dissipation of HRSC methods (see below).

The performance of a HRSC method based on a relativistic Riemann solver is illustrated by means of an MPEG movie (Fig. 4) for the planar shock heating problem for an inflow velocity v₁ = -0.99999 c (W₁ ≈ 223). These results are obtained with the relativistic PPM code of [109], which uses an exact Riemann solver based on the procedure described in Section 2.3.

Figure 4

mpg-Movie (711 KB) Still from an MPEG movie showing the evolution of the density distribution for the shock heating problem with an inflow velocity v₁ = -0.99999 c in Cartesian coordinates. The reflecting wall is located at x = 0. The adiabatic index of the gas is 4/3. For numerical reasons, the specific internal energy of the inflowing cold gas is set to a small finite value (ɛ₁ = 10^-7W_i). The figure also shows the analytical solution (blue lines). The simulation has been performed on an equidistant grid of 100 zones. (For video see appendix)

The shock wave is resolved by three zones and there are no post-shock numerical oscillations. The density increases by a factor ≈ 900 across the shock. Near x = 0 the density distribution slightly undershoots the analytical solution (by ≈ 8%) due to the numerical effect of wall heating. The profiles obtained for other inflow velocities are qualitatively similar. The mean relative error of the compression ratio σ_error < 10^-3, and, in agreement with other codes based on a Riemann solver, the accuracy of the results does not exhibit any significant dependence on the Lorentz factor of the inflowing gas.

Some authors have considered the problem of shock heating in cylindrical or spherical geometry using adapted coordinates to test the numerical treatment of geometrical factors [157, 111, 187]. Aloy et al. [3] have considered the spherically symmetric shock heating problem in 3D Cartesian coordinates as a test case for both the directional splitting and the symmetry properties of their code GENESIS. The code is able to handle this test up to inflow Lorentz factors of the order of 700.

In the shock reflection test conventional schemes often give numerical approximations which exhibit a consistent O(1) error for the density and internal energy in a few cells near the reflecting wall. This ‘overheating”, as it is known in classical hydrodynamics [130], is a numerical artifact which is considerably reduced when Marquina’s scheme is used [44]. In passing we note that the strong overheating found by Noh [130] for the spherical shock reflection test using PPM (Fig. 24 in [130]) is not a problem of PPM, but of his implementation of PPM. When properly implemented PPM gives a density undershoot near the origin of about 9% in case of a non-relativistic flow. PLM gives an undershoot of 14% in case of ultra-relativistic flows (e.g., Tab. 1 and Fig. 1 in [157]).

6.2 Propagation of relativistic blast waves

Riemann problems with large initial pressure jumps produce blast waves with dense shells of material propagating at relativistic speeds (see Fig. 5). For appropriate initial conditions, both the speed of the leading shock front and the velocity of the shell material approach the speed of light producing very narrow structures. The accurate description of these thin, relativistic shells involving large density contrasts is a challenge for any numerical code. Some particular blast wave problems have become standard numerical tests. Here we consider the two most common of these tests. The initial conditions are given in Table 5.

Figure 5

Generation and propagation of a relativistic blast wave (schematic). The large pressure jump at a discontinuity initially located at r = 0.5 gives rise to a blast wave and a dense shell of material propagating at relativistic speeds. For appropriate initial conditions both the speed of the leading shock front and the velocity of the shell approach the speed of light producing very narrow structures.

Table 5 Initial data (pressure p, density ρ, velocity v) for two common relativistic blast wave test problems. The decay of the initial discontinuity leads to a shock wave (velocity v_shock, compression ratio σ_shock) and the formation of a dense shell (velocity v_shell, time-dependent width w_shell) both propagating to the right. The gas is assumed to be ideal with an adiabatic index γ = 5/3.

Problem 1 was a demanding problem for relativistic hydrodynamic codes in the mid eighties [28, 75], while Problem 2 is a challenge even for today’s state-of-the-art codes. The analytical solution of both problems can be obtained with program the RIEMANN (see Section 9.3). Generation and propagation of relativistic blast waves

6.2.1 Problem 1

In Problem 1, the decay of the initial discontinuity gives rise to a dense shell of matter with velocity v_shell = 0.72 (W_shell = 1.38) propagating to the right. The shell trailing a shock wave of speed v_shock = 0.83 increases its width, w_shell, according to wshell = 0.11 t, i.e., at time t = 0.4 the shell covers about 4% of the grid (0 ≤ x ≤ 1). Tables 6 and 7 give a summary of the references where this test was considered for non-HRSC and HRSC methods, respectively.

Table 6 Non-HRSC methods — Summary of references where the blast wave Problem 1 (defined in Table 5) has been considered in 1D, 2D, and 3D, respectively. The methods are described in Sections 3 and 4 and their basic properties summarized in Section 5 (Tables 2, 3). Note: CD stands for contact discontinuity.

Table 7 HRSC methods — Summary of references where the blast wave Problem 1 (defined in Table 5) has been considered in 1D, 2D, and 3D, respectively. The methods are described in Sections 3 and 4 and their basic properties summarized in Section 5 (Tables 2, 3). Note: CD stands for contact discontinuity.

Using artificial viscosity techniques, Centrella & Wilson [28] were able to reproduce the analytical solution with a 7% overshoot in v_shell, whereas Hawley et al. [75] got a 16% error in the shell density.

The results obtained with early relativistic SPH codes [102] were affected by systematic errors in the rarefaction wave and the constant states, large amplitude spikes at the contact discontinuity and large smearing. Smaller systematic errors and spikes are obtained with Laguna et al.’s (1993) code [89]. This code also leads to a large overshoot in the shell’s density. Much cleaner states are obtained with the methods of Chow & Monaghan (1997) [30] and Siegler & Riffert (1999) [164], both based on conservative formulations of the SPH equations. In the case of Chow & Monaghan’s (1997) method [30], the spikes at the contact discontinuity disappear but at the cost of an excessive smearing. Shock profiles with relativistic SPH codes are more smeared out than with HRSC methods covering typically more than 10 zones.

Van Putten has considered a similar initial value problem with somewhat more extreme conditions (v_shell ≈ 0.82 c, σ_shock ≈ 5.1) and with a transversal magnetic field. For suitable choices of the smoothing parameters his results are accurate and stable, although discontinuities appear to be more smeared than with typical HRSC methods (6 – 7 zones for the strong shock wave; ≈ 50 zones for the contact discontinuity).

An MPEG movie (Figure 6) shows the Problem 1 blast wave evolution obtained with a modern HRSC method (the relativistic PPM method introduced in Section 3.1). The grid has 400 equidistant zones, and the relativistic shell is resolved by 16 zones. Because of both the high order accuracy of the method in smooth regions and its small numerical diffusion (the shock is resolved with 4–5 zones only) the density of the shell is accurately computed (errors less than 0.1%). Other codes based on relativistic Riemann solvers [50] give similar results (see Table 7). The relativistic HLL method [161] underestimates the density in the shell by about 10% in a 200 zone calculation.

Figure 6

mpg-Movie (436 KB) Still from an MPEG movie showing the evolution of the density distribution for the relativistic blast wave Problem 1 (defined in Table 5). This figure also shows the analytical solution (blue lines). The simulation has been performed with relativistic PPM on an equidistant grid of 400 zones. (For video see appendix)

6.2.2 Problem 2

Problem 2 was first considered by Norman & Winkler [131]. The flow pattern is similar to that of Problem 1, but more extreme. Relativistic effects reduce the post-shock state to a thin dense shell with a width of only about 1% of the grid length at t = 0.4. The fluid in the shell moves with v_shell = 0.960 (i.e., W_shell = 3.6), while the leading shock front propagates with a velocity v_shock = 0.986 (i.e., W_shock = 6.0). The jump in density in the shell reaches a value of 10.6. Norman & Winkler [131] obtained very good results with an adaptive grid of 400 zones using an implicit hydro-code with artificial viscosity. Their adaptive grid algorithm placed 140 zones of the available 400 zones within the blast wave thereby accurately capturing all features of the solution.

Several HRSC methods based on relativistic Riemann solvers have used Problem 2 as a standard test [107, 106, 109, 55, 187, 43]. Table 8 gives a summary of the references where this test was considered.

Table 8 Summary of references where the blast wave Problem 2 (defined in Table 5) has been considered. Methods are described in Sections 3 and 4 and their basic properties summarized in Section 5 (Tables 2, 3).

An MPEG movie (Fig. 7) shows the Problem 2 blast wave evolution obtained with the relativistic PPM method introduced in Section 3.1 on a grid of 2000 equidistant zones. At this resolution the relativistic PPM code yields a converged solution. The method of Falle & Komissarov [55] requires a seven-level adaptive grid calculation to achieve the same, the finest grid spacing corresponding to a grid of 3200 zones. As their code is free of numerical diffusion and dispersion, Wen et al. [187] are able to handle this problem with high accuracy (see Fig 8). At lower resolution (400 zones) the relativistic PPM method only reaches 69% of the theoretical shock compression value (54% in case of the second-order accurate upwind method of Falle & Komissarov [55]; 60% with the code of Donat et al. [43]).

Figure 7

mpg-Movie (3.64 MB) Still from an MPEG movie showing the evolution of the density distribution for the relativistic blast wave Problem 2 (defined in Table 5). This figure also shows the analytical solution (blue lines). The simulation has been performed with relativistic PPM on an equidistant grid of 2000 zones. (For video see appendix)

Figure 8

Results from [187] for the relativistic blast wave Problems 1 (left column) and 2 (right column), respectively. The relativistic Glimm method is only used in regions with steep gradients. Standard finite difference schemes are applied in the smooth remaining part of the computational domain. In the above plots, Lax and LW stand respectively for Lax and Lax-Wendroff methods; G refers to the pure Glimm method.

Chow & Monaghan [30] have considered Problem 2 to test their relativistic SPH code. Besides a 15% overshoot in the shell’s density, the code produces a non-causal blast wave propagation speed (i.e., v_shock > 1).

6.2.3 Collision of two relativistic blast waves

The collision of two strong blast waves was used by Woodward & Colella [191] to compare the performance of several numerical methods in classical hydrodynamics. In the relativistic case, Yang et al. [194] considered this problem to test the high-order extensions of the relativistic beam scheme, whereas Martí & Müller [109] used it to evaluate the performance of their relativistic PPM code. In this last case, the original boundary conditions were changed (from reflecting to outflow) to avoid the reflection and subsequent interaction of rarefaction waves, allowing for a comparison with an analytical solution. In the following we summarize the results on this test obtained by Martí & Müller in [109].

The initial data corresponding to this test, consisting of three constant states with large pressure jumps at the discontinuities separating the states (at x = 0.1 and x = 0. 9), as well as the properties of the blast waves created by the decay of the initial discontinuities, are listed in Table 9. The propagation velocity of the two blast waves is slower than in the Newtonian case, but very close to the speed of light (0.9776 and -0.9274 for the shock wave propagating to the right and left, respectively). Hence, the shock interaction occurs later (at t = 0. 420) than in the Newtonian problem (at t = 0.028). The top panel in Fig. 9 shows four snapshots of the density distribution, including the moment of the collision of the blast waves at t = 0.420 and x = 0.5106. At the time of collision the two shells have a width of Δx = 0.008 (left shell) and Δx = 0.019 (right shell), respectively, i.e., the whole interaction takes place in a very thin region (about 10 times smaller than in the Newtonian case, where Δx ≈ 0.2).

Table 9 Initial data (pressure p, density ρ, velocity v) for the test problem of two colliding relativistic blast waves. The de cay of the initial discontinuities (at x = 0.1 and x = 0.9) produces two shock waves (velocities v_shock, compression ratios σ_shock) moving in opposite directions followed by two trailing dense shells (velocities v_shell, time-dependent widths w_shell). The gas is assumed to be ideal with an adiabatic index γ = 1.4.

Figure 9

The top panel shows a sequence of snapshots of the density profile for the colliding relativistic blast wave problem up to the moment when the waves begin to interact. The density profile of the new states produced by the interaction of the two waves is shown in the bottom panel (note the change in scale on both axes with respect to the top panel).

The collision gives rise to a narrow region of very high density (see lower panel of Fig. 9), bounded by two shocks moving at speeds 0.088 (shock at the left) and 0.703 (shock at the right) and large compression ratios (7.26 and 12.06, respectively) well above the classical limit for strong shocks (6.0 for γ = 1.4). The solution just described applies until t = 0. 430 when the next interaction takes place.

The complete analytical solution before and after the collision up to time t = 0.430 can be obtained following Appendix II in [109].

An MPEG movie (Fig. 10) shows the evolution of the density up to the time of shock collision at t = 0. 4200. The movie was obtained with the relativistic PPM code of Martí & Müller [109]. The presence of very narrow structures with large density jumps requires very fine zoning to resolve the states properly. For the movie a grid of 4000 equidistant zones was used. The relative error in the density of the left (right) shell is always less than 2.0% (0.6%), and is about 1.0% (0.5%) at the moment of shock collision. Profiles obtained with the relativistic Godunov method (first-order accurate, not shown) show relative errors in the density of the left (right) shell of about 50% (16%) at t = 0. 20. The errors drop only slightly to about 40% (5%) at the time of collision (t = 0.420).

Figure 10

mpg-Movie (2.00 MB) Still from an MPEG movie showing the evolution of the density distribution for the colliding relativistic blast wave problem up to the interaction of the waves. This figure also shows the analytical solution (blue lines). The computation has been performed with relativistic PPM on an equidistant grid of 4000 zones. (For video see appendix)

An MPEG movie (Fig. 11) shows the numerical solution after the interaction has occurred. Compared to the other MPEG movie (Fig. 10) a very different scaling for the x-axis had to be used to display the narrow dense new states produced by the interaction. Obviously, the relativistic PPM code resolves the structure of the collision region satisfactorily well, the maximum relative error in the density distribution being less than 2.0%. When using the first-order accurate Godunov method instead, the new states are strongly smeared out and the positions of the leading shocks are wrong.

Figure 11

mpg-Movie (697 KB) Still from an MPEG movie showing the evolution of the density distribution for the colliding relativistic blast wave problem around the time of interaction of the waves at an enlarged spatial scale. This figure also shows the analytical solution (blue lines). The computation has been performed with relativistic PPM on an equidistant grid of 4000 zones. (For video see appendix)

7 Applications

7.1 Astrophysical jets

The most compelling case for a special relativistic phenomenon are the ubiquitous jets in extragalactic radio sources associated with active galactic nuclei. In the commonly accepted standard model [10], flow velocities as large as 99% of the speed of light (in some cases even beyond) are required to explain the apparent superluminal motion observed in many of these sources. Models which have been proposed to explain the formation of relativistic jets, involve accretion onto a compact central object, such as a neutron star or stellar mass black hole in the galactic micro-quasars GRS 1915+105 [118] and GRO J1655-40 [174], or a rotating super massive black hole in an active galactic nucleus, which is fed by interstellar gas and gas from tidally disrupted stars.

Inferred jet velocities close to the speed of light suggest that jets are formed within a few gravitational radii of the event horizon of the black hole. Moreover, very-long-baseline interferometric (VLBI) radio observations reveal that jets are already collimated at subparsec scales. Current theoretical models assume that accretion disks are the source of the bipolar outflows which are further collimated and accelerated via MHD processes (see, e.g., [16]). There is a large number of parameters which are potentially important for jet powering: the black hole mass and spin, the accretion rate and the type of accretion disk, the properties of the magnetic field and of the environment.

At parsec scales the jets, observed via their synchrotron and inverse Compton emission at radio frequencies with VLBI imaging, appear to be highly collimated with a bright spot (the core) at one end of the jet and a series of components which separate from the core, sometimes at superluminal speeds. In the standard model [17], these speeds are interpreted as a consequence of relativistic bulk motions in jets propagating at small angles to the line of sight with Lorentz factors up to 20 or more. Moving components in these jets, usually preceded by outbursts in emission at radio wavelengths, are interpreted in terms of traveling shock waves.

Finally, the morphology and dynamics of jets at kiloparsec scales are dominated by the interaction of the jet with the surrounding extragalactic medium, the jet power being responsible for dichotomic morphologies (the so called Fanaroff-Riley I and II classes [56], FR I and FR II, respectively). Current models [14, 91] interpret FR I morphologies as the result of a smooth deceleration from relativistic to non-relativistic, transonic speeds on kpc scales due to a slower shear layer. For the most powerful radio galaxies (FR II) and quasars on the other hand, the observation of flux asymmetries between jet and counter-jet indicates that in these sources relativistic motion extends up to kpc scales, although with smaller values of the overall bulk speeds [21].

Although MHD and general relativistic effects seem to be crucial for a successful launch of the jet (for a review see, e.g., [23]), purely hydrodynamic, special relativistic simulations are adequate to study the morphology and dynamics of relativistic jets at distances sufficiently far from the central compact object (i.e., at parsec scales and beyond). The development of relativistic hydro-dynamic codes based on HRSC techniques (see Sections 3 and 4) has triggered the numerical simulation of relativistic jets at parsec and kiloparsec scales.

At kiloparsec scales, the implications of relativistic flow speeds and / or relativistic internal energies for the morphology and dynamics of jets have been the subject of a number of papers in recent years [112, 46, 110, 111, 86]. Beams with large internal energies show little internal structure and relatively smooth cocoons allowing the terminal shock (the hot spot in the radio maps) to remain well defined during the evolution. Their morphologies resemble those observed in naked quasar jets like 3C273 [37]. Fig. 12 shows several snapshots of the time evolution of a light, relativistic jet with large internal energy. The dependence of the beam’s internal structure on the flow speed suggests that relativistic effects may be relevant for the understanding of the difference between slower, knotty BL Lac jets and faster, smoother quasar jets [60].

Figure 12

Time evolution of a light, relativistic (beam flow velocity equal to 0.99) jet with large internal energy. The logarithm of the proper rest-mass density is plotted in grey scale, the maximum value corresponding to white and the minimum to black.

Highly supersonic models, in which kinematic relativistic effects due to high beam Lorentz factors dominate, have extended over-pressured cocoons. These over-pressured cocoons can help to confine the jets during the early stages of their evolution [110] and even cause their deflection when propagating through non-homogeneous environments [148]. The cocoon overpressure causes the formation of a series of oblique shocks within the beam in which the synchrotron emission is enhanced. In long term simulations (see Fig. 13), the evolution is dominated by a strong deceleration phase during which large lobes of jet material (like the ones observed in many FR IIs, e.g., Cyg A [25]) start to inflate around the jet’s head. These simulations reproduce some properties observed in powerful extragalactic radio jets (lobe inflation, hot spot advance speeds and pressures, deceleration of the beam flow along the jet) and can help to constrain the values of basic parameters (such as the particle density and the flow speed) in the jets of real sources.

Figure 13

Logarithm of the proper rest-mass density and energy density (from top to bottom) of an evolved, powerful jet propagating through the intergalactic medium. The white contour encompasses the jet material responsible for the synchrotron emission.

The development of multidimensional relativistic hydrodynamic codes has allowed, for the first time, the simulation of parsec scale jets and superluminal radio components [68, 85, 117]. The presence of emitting flows at almost the speed of light enhances the importance of relativistic effects in the appearance of these sources (relativistic Doppler boosting, light aberration, time delays). Hence, one should use models which combine hydrodynamics and synchrotron radiation transfer when comparing with observations. In these models, moving radio components are obtained from perturbations in steady relativistic jets. Where pressure mismatches exist between the jet and the surrounding atmosphere reconfinement shocks are produced. The energy density enhancement produced downstream from these shocks can give rise to stationary radio knots as observed in many VLBI sources. Superluminal components are produced by triggering small perturbations in these steady jets which propagate at almost the jet flow speed. One example of this is shown in Fig. 14 (see also [68]), where a superluminal component (apparent speed ≈ 7 times the speed of light) is produced from a small variation of the beam flow Lorentz factor at the jet inlet. The dynamic interaction between the induced traveling shocks and the underlying steady jet can account for the complex behavior observed in many sources [67].

Figure 14

Computed radio maps of a compact relativistic jet showing the evolution of a superluminal component (from left to right). Two resolutions are shown: present VLBI resolution (white contours) and resolution provided by the simulation (black/white images).

The first magnetohydrodynamic simulations of relativistic jets have been already undertaken in 2D [82, 81] and 3D [128, 129] to study the implications of ambient magnetic fields in the morphology and bending properties of relativistic jets. However, despite the impact of these results in specific problems like, e.g., the understanding of the misalignment of jets between pc and kpc scales, these 3D simulations have not addressed the effects on the jet structure and dynamics of the third spatial degree of freedom. This has been the aim of the work undertaken by Aloy et al. [2].

Finally, Koide et al. [83] have developed a general relativistic MHD code and applied it to the problem of jet formation from black hole accretion disks. Jets are formed with a two-layered shell structure consisting of a fast gas pressure driven jet (Lorentz factor ≈ 2) in the inner part and a slow magnetically driven outflow in the outer part, both of which are being collimated by the global poloidal magnetic field penetrating the disk.

7.2 Gamma-Ray Bursts (GRBs)

A second phenomenon which involves flows with velocities very close to the speed of light are gamma-ray bursts (GRBs). Although known observationally for over 30 years, until recently their distance ("local” or “cosmological”) has been, and their nature still is, a matter of controversial debate [57, 115, 143, 144]. GRBs do not repeat except for a few soft gamma-ray repeaters. They are detected with a rate of about one event per day, and their duration varies from milliseconds to minutes. The duration of the shorter bursts and the temporal substructure of the longer bursts implies a geometrically small source (less than ∼ c ¿ 1 msec ∼ 100 km), which in turn points towards compact objects, like neutron stars or black holes. The emitted gamma-rays have energies in the range 30 keV to 2 MeV.

Concerning the distance of GRB sources major progress has occurred through the observations by the BATSE detector on board the Compton Gamma-Ray Observatory (GRO), which have proven that GRBs are distributed isotropically over the sky [114]. Even more important the detection and the rapid availability of accurate coordinates (∼ arc minutes) of the fading X-ray counterparts of GRBs by the BeppoSAX spacecraft beginning in 1997 [34, 146], has allowed for subsequent successful ground based observations of faint GRB afterglows at optical and radio wavelength. In the case of GRB 990123 the optical, X-ray and gamma-ray emission was detected for the first time almost simultaneously (optical observations began 22 seconds after the onset of the GRB) [22, 1]. From optical spectra thus obtained, redshifts of several gamma-ray bursts have been determined, e.g., GRB 970508 (z = 0.835 [116, 141]), GRB 971214 (z = 3.42 [87]), GRB 980703 (z = 0.966 [41]), and GRB 990123 (1.60 ≤ z ≤ 2.05 [5]), which confirm that (at least some) GRBs occur at cosmological distances. Assuming isotropic emission the inferred total energy of cosmological GRBs emitted in form of gamma-rays ranges from several 10⁵¹ erg to 3 ¿ 10⁵³ erg (for GRB 971214) [26], and exceeds 10⁵⁴ erg for GRB 990123 [5, 22]. Updated information on GRBs localized with BeppoSAX, BATSE / RXTE (PCA) or BATSE / RXTE (ASM) can be obtained from a web site maintained by Greiner [71].

The compact nature of the GRB source, the observed flux, and the cosmo-logical distance taken together imply a large photon density. Such a source has a large optical depth for pair production. This is, however, inconsistent with the optically thin source indicated by the non-thermal gamma-ray spectrum, which extends well beyond the pair production threshold at 500 keV. This problem can be resolved by assuming an ultra-relativistic expansion of the emitting region, which eliminates the compactness constraint. The bulk Lorentz factors required are then W > 100 (see, e.g., [144]).

In April 1998 the pure cosmological origin of GRBs was challenged by the detection of the Type Ib/c supernova SN 1998bw [61, 62] within the 8 arc minute error box of GRB 980425 [165, 140]. Its explosion time is consistent with that of the GRB, and relativistic expansion velocities are derived from radio observations of SN 1998bw [88]. BeppoSAX detected two fading X-ray sources within the error box, one being positionally consistent with the supernova and a fainter one not consistent with the position of SN 1998bw [140]. Taken together these facts suggest a relationship between GRBs and SNe Ib/c, i.e., core collapse supernovae of massive stellar progenitors which have lost their hydrogen and helium envelopes [62, 78, 193]. As the host galaxy ESO 184-82 of SN 1998bw is only at a redshift of z = 0.0085 [175] and as GRB 980425 was not extraordinarily bright, GRB-supernovae are more than four orders of magnitude fainter (E_tot γ = 7 ¿ 10⁴⁷ erg for GRB 980425 [26]) than a typical cosmological GRB. However, the observation of the second fading X-ray source within the error box of GRB 980425 and unrelated with SN 1998bw still causes some doubts on the GRB supernova connection, although the probability of chance coincidence of GRB 980425 and SN 1998bw is extremely low [140].

In order to explain the energies released in a GRB various catastrophic collapse events have been proposed including neutron-star/neutron-star mergers [134, 69, 47], neutron-star/black-hole mergers [119], collapsars [192, 101], and hypernovae [135]. These models all rely on a common engine, namely a stellar mass black hole which accretes several solar masses of matter from a disk (formed during a merger or by a non-spherical collapse) at a rate of ∼ 1M_⊙ s^-1 [151]. A fraction of the gravitational binding energy released by accretion is converted into neutrino and anti-neutrino pairs, which in turn annihilate into electron-positron pairs. This creates a pair fireball, which will also include baryons present in the environment surrounding the black hole. Provided the baryon load of the fireball is not too large, the baryons are accelerated together with the e+ e- pairs to ultra-relativistic speeds with Lorentz factors > 10² [27, 145, 144]. The existence of such relativistic flows is supported by radio observations of GRB 980425 [88]. It has been further argued that the rapid temporal decay of several GRB afterglows is inconsistent with spherical (isotropic) blast wave models, and instead is more consistent with the evolution of a relativistic jet after it slows down and spreads laterally [160]. Independent of the flow pattern the bulk kinetic energy of the fireball then is thought to be converted into gamma-rays via cyclotron radiation and/or inverse Compton processes (see, e.g., [115, 144]).

One-dimensional numerical simulations of spherically symmetric relativistic fireballs have been performed by several authors to model GRB sources [145, 137, 136]. Multi-dimensional modeling of ultra-relativistic jets in the context of GRBs has for the first time been attempted by Aloy et al. [4]. Using a collapsar progenitor model of MacFadyen & Woosley [101] they have simulated the propagation of an axisymmetric jet through the mantle and envelope of a collapsing massive star (10M_⊙) using the GENESIS special relativistic hydrodynamic code [3]. The jet forms as a consequence of an assumed energy deposition of 10⁵¹ erg/sec within a 30 degree cone around the rotation axis. At break-out, i.e., when the jet reaches the surface of the stellar progenitor, the maximum Lorentz factor of the jet flow is about 20. The latter fact implies that Newtonian simulations of this phenomenon [101] are clearly inadequate.

8 Conclusion

8.1 Evaluation of the methods

An assessment of the quality of the numerical methods should consider, at least, the following aspects:

• (i) accuracy and robustness in describing high Lorentz factor flows with strong shocks;

• (ii) effort required to extend to multi dimensions;

• (iii) effort required to extend to RMHD and GRHD.

In Table 10 we have summarized these aspects of numerical methods for SRHD.

Table 10 Evaluation of numerical methods for SRHD. Methods have been categorized for clarity.

Since their introduction in numerical RHD at the beginning of nineties, HRSC methods have demonstrated their ability to describe accurately (stable and without excessive smearing) relativistic flows of arbitrarily large Lorentz factors and strong discontinuities, reaching the same quality as in classical hydrodynamics. In addition (as it is the case for classical flows, too), HRSC methods show the best performance compared to any other method (e.g., artificial viscosity, FCT or SPH).

Despite of the latter fact, a lot of effort has been put into improving these non-HRSC methods. Using a consistent formulation of artificial viscosity has significantly enhanced the capability of finite difference schemes [131] as well as of relativistic SPH [164] to handle strong shocks without spurious post-shock oscillations. However, this comes at the price of a large numerical dissipation at shocks. Concerning relativistic SPH, recent investigations using a conservative formulation of the hydrodynamic equations [30, 164] have reached an unprecedented accuracy with respect to previous simulations, although some issues still remain. Besides the strong smearing of shocks, the description of contact discontinuities and of thin structures moving at ultra-relativistic speeds needs to be improved (see Section 6.2).

Concerning FCT techniques, those codes based on a conservative formulation of the RHD equations have been able to handle relativistic flows with discontinuities at all flow speeds, although the quality of the results is below that of HRSC methods in all cases [161].

The extension to multi-dimensions is simple for most relativistic codes. Finite difference techniques are easily extended using directional splitting. Note, however, that HRSC methods based on exact solutions of the Riemann problem [109, 187] first require the development of a multidimensional version of the relativistic Riemann solver. The adapting-grid, artificial viscosity, implicit code of Norman & Winkler [131] and the relativistic Glimm method of Wen et al. [187] are restricted to one dimensional flows. Note that Glimm’s method produces the best results in all the tests analyzed in Section 6.

The symmetric TVD scheme proposed by Davis [38] and extended to GRMHD (see below) by Koide et al. [82] combines several characteristics making it very attractive. It is written in conservation form and is TVD, i.e., it is converging to the physical solution. In addition, it is independent of spectral decompositions, which allows for a simple extension to RMHD. Quite similar statements can be made about the approach proposed by van Putten [181]. In contrast to FCT schemes (which are also easily extended to general systems of equations), both Koide et al.’s and van Putten’s methods are very stable when simulating mildly relativistic flows (maximum Lorentz factors ≈ 4) with discontinuities. Their only drawback is an excessive smearing of the latter. A comparison of Davis’ method with Riemann solver based methods would be desirable.

8.2 Further developments

The directions of future developments in this field of research are quite obvious. They can be divided into four main categories:

8.2.1 Incorporation of realistic microphysics

Up to now most astrophysical SRHD simulations have assumed matter whose thermodynamic properties can be described by an inviscid ideal equation of state with a constant adiabatic index. This simplification may have been appropriate in the first generation of SRHD simulations, but it clearly must be given up when aiming at a more realistic modeling of astrophysical jets, gamma-ray burst sources or accretion flows onto compact objects. For these phenomena a realistic equation of state should include contributions from radiation (γ = 4/3-“fluid”), allow for the formation of electron-positron pairs at high temperatures, allow the ideal gas contributions to be arbitrarily degenerate and/or relativistic. Depending on the problem to be simulated, effects due to heat conduction, radiation transport, cooling, nuclear reactions, and viscosity may have to be considered, too. To include any of these effects is often a non trivial task even in Newtonian hydrodynamics (see, e.g., the contributions in the book edited by Steiner & Gautschy [168]).

When simulating relativistic heavy ion collisions, the use of a realistic equation of state is essential for an adequate description of the phenomenon. However, as these simulations have been performed with FCT based difference schemes (see, e.g., [166]), this poses no specific numerical problem. The simulation of flows obeying elaborated microphysics with HRSC methods needs in some cases the extension of the present relativistic Riemann solvers to handle general equations of state. This is the case of the Roe-Eulderink method (extensible by the procedure developed in the classical case by Glaister [64]), and rPPM and rGlimm both relying on an exact solution of the Riemann problem for ideal gases with constant adiabatic exponent (which can also be extended following the procedure of Colella & Glaz [32] for classical hydrodynamics). We expect the second generation of SRHD codes to be capable of treating general equations of state and various source/sink terms routinely.

Concerning the usage of complex equations of state (EOS) a limitation must be pointed out which is associated with the Riemann solvers used in HRSC methods, even in the Newtonian limit. These problems are especially compounded in situations where there are phase transitions present. In this case the EOS may have a discontinuous adiabatic exponent and may even be non-convex. The Riemann solver of Colella & Glaz [32] often fails in these situations, because it is derived under the assumption of convexity in the EOS. When convexity is not present the character of the solution to the Riemann problem changes. Situations where phase transitions cause discontinuities in the adiabatic index or non-convexity of the EOS are encountered, e.g., in simulations of neutron star formation, of the early Universe, and of relativistic heavy ion collisions.

Another interesting area that deserves further research is the application of relativistic HRSC methods in simulations of reactive multi-species flows, especially as such flows still cause problems for the Newtonian CFD community (see, e.g., [149]). The structure of the solution to the Riemann problem becomes significantly more complex with the introduction of reactions between multiple species. Riemann solvers that incorporate source terms [97], and in particular source terms due to reactions, have been proposed for classical flows [11, 79], but most HRSC codes still rely on operator splitting.

8.2.2 Coupling of SRHD schemes with AMR

Modeling astrophysical phenomena often involves an enormous range of length scales and time scales to be covered in the simulations (see, e.g., [124]). In two and definitely in three spatial dimensions many such simulations cannot be performed with sufficient spatial resolution on a static equidistant or non-equidistant computational grid, but they will require dynamic, adaptive grids. In addition, when the flow problem involves stiff source terms (e.g., energy generation by nuclear reactions) very restrictive time step limitations may result. A promising approach to overcome these complications will be the coupling of SRHD solvers with the adaptive mesh refinement (AMR) technique [13]. AMR automatically increases the grid resolution near flow discontinuities or in regions of large gradients (of the flow variables) by introducing a dynamic hierarchy of grids until a prescribed accuracy of the difference approximation is achieved. Because each level of grids is evolved in AMR on its own time step, time step restrictions due to stiff source terms are constraining the computational costs less than without AMR. For an overview of online information about AMR visit, e.g., the AMRA home page of Plewa [147], and for public domain AMR software, e.g., the AMRCLAW home page of LeVeque & Berger [99], and the AMRCART home page of Walder [186].

A SRHD simulation of a relativistic jet based on a combined HLL-AMR scheme was performed by Duncan & Hughes [46]. Plewa et al. [148] have modeled the deflection of highly supersonic jets propagating through non-homogeneous environments using the HRSC scheme of Martí et al. [111] combined with the AMR implementation AMRA of Plewa [147]. Komissarov & Falle [85] have combined their numerical scheme with the adaptive grid code Cobra, which has been developed by Mantis Numerics Ltd. for industrial applications [54], and which uses a hierarchy of grids with a constant refinement factor of two between subsequent grid levels.

8.2.3 General relativistic hydrodynamics (GRHD)

Up to now only very few attempts have been made to extend HRSC methods to GRHD and all of these have used linearized Riemann solvers [107, 50, 157, 9, 59]. In the most recent of these approaches Font et al. [59] have developed a 3D general relativistic HRSC hydrodynamic code where the matter equations are integrated in conservation form and fluxes are calculated with Marquina’s formula.

A very interesting and powerful procedure was proposed by Balsara [8] and has been implemented by Pons et al. [150]. This procedure allows one to exploit all the developments in the field of special relativistic Riemann solvers in general relativistic hydrodynamics. The procedure relies on a local change of coordinates at each zone interface such that the spacetime metric is locally flat. In that locally flat spacetime any special relativistic Riemann solver can be used to calculate the numerical fluxes, which are then transformed back. The transformation to an orthonormal basis is valid only at a single point in spacetime. Since the use of Riemann solvers requires the knowledge of the behavior of the characteristics over a finite volume, the use of the local Lorentz basis is only an approximation. The effects of this approximation will only become known through the study of the performance of these methods in situations where the structure of the spacetime varies rapidly in space and perhaps time as well. In such a situation finer grids and improved time advancing methods will definitely be required. The implementation is simple and computationally inexpensive.

Characteristic formulations of the Einstein field equations are able to handle the long term numerical description of single black hole spacetimes in vacuum [15]. In order to include matter in such an scenario, Papadopoulos & Font [138] have generalized the HRSC approach to cope with the hydrodynamic equations in such a null foliation of spacetime. Actually, they have presented a complete (covariant) re-formulation of the equations in GR, which is also valid for spacelike foliations in SR. They have extensively tested their method calculating, among other tests, shock tube problem 1 (see Section 6.2.1), but posed on a light cone and using the appropriate transformations of the exact solution [108] to account for advanced and retarded times.

Other developments in GRHD in the past included finite element methods for simulating spherically symmetric collapse in general relativity [103], general relativistic pseudo-spectral codes based on the (3+1) ADM formalism [7] for computing radial perturbations [70] and 3D gravitational collapse of neutron stars [19], and general relativistic SPH [102]. The potential of these methods for the future is unclear, as none of them is specifically appropriate for ultra-relativistic speeds and strong shock waves which are characteristic of most astrophysical applications.

Peitz & Appl [139] have addressed the difficult issue of non-ideal GRHD, which is of particular importance, e.g., for the simulation of accretion discs around compact objects, rotating relativistic fluid configurations, and the evolution of density fluctuations in the early universe. They have accounted for dissipative effects by applying the theory of extended causal thermodynamics, which eliminates the causality violating infinite signal speeds arising from the conventional Navier-Stokes equation. Peitz & Appl have not implemented their model numerically yet.

8.2.4 Relativistic magneto-hydrodynamics (RMHD)

The inclusion of magnetic effects is of great importance in many astrophysical flows. The formation and collimation process of (relativistic) jets most likely involves dynamically important magnetic fields and occurs in strong gravitational fields. The same is likely to be true for accretion discs around black holes. Magneto-relativistic effects even play a non-negligible role in the formation of proto-stellar jets in regions close to the light cylinder [23]. Thus, relativistic MHD codes are a very desirable tool in astrophysics. The non-trivial task of developing such a kind of code is considerably simplified by the fact that because of the high conductivity of astrophysical plasmas one must only consider ideal RMHD in most applications.

Evans & Hawley [52] extended the second-order accurate, Newtonian, artificial-viscosity transport method of Hawley et al. [75] to the evolution of the MHD induction equation. Special relativistic 2D MHD test problems with Lorentz factors up to ∼ 3 have been investigated by Dubal [45] with a code based on FCT techniques (see Section 4).

In a series of papers Koide and coworkers [82, 81, 128, 129, 83] have investigated relativistic magnetized jets using a symmetric TVD scheme (see Section 3). Koide, Nishikawa & Mutel [82] simulated a 2D RMHD slab jet, whereas Koide [81] investigated the effect of an oblique magnetic field on the propagation of a relativistic slab jet. Nishikawa et al. [128, 129] extended these simulations to 3D and considered the propagation of a relativistic jet with a Lorentz factor W = 4.56 along an aligned and an oblique external magnetic field. The 2D and 3D simulations published up to now only cover the very early propagation of the jet (up to 20 jet radii) and are performed with moderate spatial resolution on an equidistant Cartesian grid (up to 101 zones per dimension, i.e., 5 zones per beam radius).

Van Putten [180, 181] has proposed a method for accurate and stable numerical simulations of RMHD in the presence of dynamically significant magnetic fields in two dimensions and up to moderate Lorentz factors. The method is based on MHD in divergence form using a 2D shock-capturing method in terms of a pseudo-spectral smoothing operator (see Section 4). He applied this method to 2D blast waves [183] and astrophysical jets [182, 184].

Steps towards the extension of linearized Riemann solvers to ideal RMHD have already been taken. Romero [158] has derived an analytical expression for the spectral decomposition of the Jacobian in the case of a planar relativistic flow field permeated by a transversal magnetic field (nonzero field component only orthogonal to flow direction). Van Putten [178] has studied the characteristic structure of the RMHD equations in (constraint free) divergence form. Finally, Komissarov [84] has presented a robust Godunov-type scheme for RMHD, which is based on a linear Riemann solver, has second-order accuracy in smooth regions, enforces magnetic flux conservation, and which can cope with ultra-relativistic flows.

We end with the simulations performed by Koide, Shibata & Kudoh [83] on magnetically driven axisymmetric jets from black hole accretion disks. Their GRMHD code is an extension of the special relativistic MHD code developed by Koide et al. [82, 81, 128]. The necessary modifications of the code were quite simple, because in the (nonrotating) black hole’s Schwarzschild spacetime the GRMHD equations are identical to the SRMHD equations in general coordinates, except for the gravitational force terms and the geometric factors of the lapse function. With the pioneering work of Koide, Shibata & Kudoh the epoch of exciting GRMHD simulations has just begun.