8.2 Present and future developments

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, and allow the ideal gas contributions to be arbitrarily degenerate and/or relativistic.

Depending on the problem to be simulated, effects due to heat conduction, diffusion, radiation transport, cooling, nuclear reactions, and viscosity may have to be considered, too. Including any of these effects is often a non-trivial task even in Newtonian hydrodynamics, as the differential operators describing advection and convection are of hyperbolic nature, while diffusion and conduction processes give rise to parabolic differential operators, and the treatment of constraints or self-gravity involves differential operators of elliptic type (see, e.g., the contributions in the book edited by Steiner and Gautschy [162]). There has been considerable development in the coupling of Newtonian HRSC methods to the nonhyperbolic terms arising in the equations from these physical processes using semi-implicit approaches, e.g., the predictor-corrector method [18]. Another example in this context provides the recent work of Howell and Grenough [127], who have coupled an explicit Newtonian Godunov-type integrator for the hyperbolic hydrodynamic equations to an implicit multigrid solver to describe effects of radiative diffusion on the flow and vice versa. We particularly mention this work here, as it also uses a block-structured adaptive mesh refinement algorithm (see Section 8.2.2). Although such sophisticated methods have not been applied in SRHD yet, they represent an important set of ideas that could provide a starting point for more elaborate SRHD simulations.

In the context of relativistic jets, Komissarov and Falle [148], and Scheck et al. [256] have considered a mixture of ideal, relativistic Boltzmann gases [272] hence allowing for jets with general (i.e., e, e⁺, p) composition. The usage of such a more general ideal EOS causes no special problem for the Riemann solvers although a higher nonlinearity is introduced in the process of the recovery of the primitive variables. In order to avoid this extra complexity, approximate expressions for the relativistic ideal gas EOS have been proposed [79, 266]. In case of the approximation proposed by Sokolov et al. [266], the recovery of the primitive variables is explicit. Moreover, the authors have developed an exact Riemann solver consistent with the approximate EOS.

An EOS describing matter consisting of a set of ideal, non-relativistic Boltzmann gases (nuclei in nuclear statistical equilibrium), a Fermi gas of electrons and photons was used in the simulations of relativistic jets from collapsars by Aloy et al. [7].

HRSC flow simulations involving elaborate microphysics may require the extension of the presently available relativistic Riemann solvers to handle general equations of state (see Section 9.1). This is the case for the Roe–Eulderink method, which can be extended following the procedure developed in the classical case by Glaister [103]. Methods based on exact solutions of the Riemann problem, like rPPM and rGlimm, can take advantage of the solution presented in Section 2.3 to cope with a general EOS. FCT based difference schemes used in simulations of relativistic heavy ion collisions (see Section 7.3) pose no specific numerical problem in handling a general EOS.

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., [233]). 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 [159], and in particular source terms due to reactions, have been proposed for classical flows [19, 132]. However, most HRSC codes still rely on operator splitting.

Peitz and Appl [222] 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. However, Peitz, and Appl have not yet implemented their model numerically.

A description of non-ideal hydrodynamics in general relativity is also the aim of Richardson and Chung’s work [243], although from a less formal basis. The authors introduce an approach (the so-called flow-field-dependent variation theory [54, 55] resting on the conservative Navier–Stokes system of equations for classical fluid dynamics) where local properties of the flow (advection, turbulence, or gravity dominated) are captured in terms of relevant parameters (measuring changes of the Lorentz factor, relativistic Reynolds and Froude numbers between adjacent numerical zones, respectively). These parameters are then used to produce a suitable numerical model (hyperbolic, parabolic, elliptic) which is subsequently discretized using finite difference or finite element methods. The latter approach has been applied by Richardson and Chung [243] for several test cases (mildly relativistic Riemann problem and general relativistic spherical dust infall).

8.2.2 Coupling of SRHD schemes with AMR

Modeling astrophysical phenomena often involves an enormous range of length and time scales to be covered in the simulations (see, e.g., [205]). 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 rather 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 is the coupling of SRHD solvers with the adaptive mesh refinement (AMR) technique [21]. 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 constrain the computational costs less than without AMR.

For an overview of online information about AMR visit, e.g., the AMRA home page of Plewa [230], and for public domain AMR software, e.g., the AMRCLAW home page of LeVeque and Berger [161], the web page of the Lawrence Berkeley Lab dedicated to AMR [1], and the NASA Goddard Space Flight Center web page on PARAMESH [171]. Astrophysical applications based on PARAMESH can be found at the web site of the ASCI / Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago [2]. Although, as demonstrated by these web sites, there has been a considerable effort over the last few years in developing frameworks for block-structured adaptive mesh refinement, we will see that the application to SRHD is still in its infancy.

An SRHD simulation of a relativistic jet based on a combined HLL-AMR scheme was performed by Duncan and Hughes [78]. Plewa et al. [232, 231] have modeled the deflection of highly supersonic slab jets propagating through non-homogeneous environments using the HRSC scheme of Martí et al. [183] combined with the AMR implementation AMRA of Plewa [230]. A similar study, but in 3D, was performed by Hughes et al. [128] who studied the deflection and precession of cylindrical relativistic jets when impinging on an oblique density gradient using the SRHD code of Duncan and Hughes [78] extended to 3D and their own implementation of the AMR technique of Berger and Colella [21]. Komissarov and Falle [147] have combined their numerical scheme with the adaptive grid code Cobra, which has been developed by Mantis Numerics Ltd. for industrial applications [88], 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 (for a comprehensive review see Font [91]). All these attempts are based on the usage of linearized Riemann solvers [179, 84, 250, 15, 93]. In the most recent of these approaches, Font et al. [93] 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 [13] and has been implemented by Pons et al. [234]. 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 [24]. In order to include matter in such an scenario, Papadopoulos and Font [221] 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) reformulation 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 [180] 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 [173], general relativistic pseudo-spectral codes based on the (3+1) ADM formalism [11] for computing radial perturbations [112] and 3D gravitational collapse of neutron stars [32], general relativistic [172, 204] and post-Newtonian [12] SPH. 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.

Let us remark that, in the case of dynamic spacetimes, the equations of relativistic hydrodynamics are solved on the local (in space and time) background solution provided by the Einstein equations at every time step [91]. The solution of the Einstein gravitational field equations and its coupling with the hydrodynamic equations is the realm of Numerical Relativity, which is beyond the scope of this article (see, e.g., Lehner [207] for a recent review).

8.2.4 Relativistic magneto-hydrodynamics (RMHD)

The inclusion of magnetic effects is of great importance for many astrophysical phenomena. The formation and collimation process of (relativistic) jets (powering powerful extragalactic radio sources, galactic microquasars, and GRBs) 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 [41]. 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.

The aim of any (Newtonian or relativistic) MHD code is to evolve the induction equation to obtain the magnetic fields from which to calculate the Lorentz force. Magnetic fields are divergence free, i.e., . Hence, numerical schemes are required to maintain this constraint (if fulfilled for the initial data) during the evolution. A first step towards the development of a relativistic (in fact, general relativistic) MHD code was made by Evans and Hawley [86] who incorporated a numerical scheme for the induction equation (constrained transport), which maintained zero divergence of the magnetic field up to machine round-off error, into the axisymmetric, two-dimensional finite difference code of Hawley et al. [123]. In Evans and Hawley’s method the magnetic flux through cell interfaces is the fundamental “magnetic” variable. Their method is also based on the use of a staggered mesh (some quantities including the magnetic field components are defined at grid interfaces). Thus, even in classical MHD, the extension of the constrained transport method to Riemann-solver-based schemes (that rely on fluxes at cell interfaces derived from cell averaged quantities) is non-trivial [65, 253]. Tóth [280] reviews and compares strategies (namely the eight-wave formulation, several versions of the constrained transport, and the projection scheme) used in HRSC schemes in classical MHD to maintain the constraint numerically. His conclusions also apply to RMHD.

Special relativistic 2D MHD test problems with Lorentz factors up to  3 have been investigated by Dubal [77] with a code based on FCT techniques (see Section 4).

Van Putten [286, 287, 290] 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 the method to 2D blast waves [289] and astrophysical jets [288, 291].

In a series of papers, Koide and coworkers [138, 136, 210, 211, 139] have investigated relativistic magnetized jets using a symmetric TVD scheme (see Section 3). Koide, Nishikawa, and Mutel [138] simulated a 2D RMHD slab jet, whereas Koide [136] investigated the effect of an oblique magnetic field on the propagation of a relativistic slab jet. Nishikawa et al. [210, 211] 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). Concerning higher order symmetric non-oscillatory schemes, the very recent work by Del Zanna et al. [72] has to be mentioned. Their third order scheme produces results which are competitive with those obtained by Riemann-solver based methods (see next paragraph) but avoiding all the complexity associated with the spectral decomposition into characteristic fields (particularly the degeneracies). Its high order and its simplicity make this approach very appealing.

Steps towards the extension of linearized Riemann solvers to ideal RMHD have already been taken. All theoretical aspects (RMHD as a quasi-linear hyperbolic system, spectral decomposition of the Jacobian of the flux vector in covariant form, study of simple waves and shock waves) are compiled in the book by Anile [8], augmenting previous work of Lichnerowicz [163]. Romero [251] derived an analytic expression for the spectral decomposition of the Jacobian matrix of the flux vector in the case of a planar relativistic flow field permeated by a transversal magnetic field (nonzero field component only orthogonal to flow direction). Anile and Pennisi [9] and Van Putten [292] studied the characteristic structure of the RMHD equations in (constraint free) covariant form. Finally, Balsara [14] and Komissarov [143] have developed robust, second-order accurate (in space and time), Godunov-type schemes for 1D and 2D RMHD, respectively. Both start from the spectral decomposition of the RMHD system of (ten) equations in covariant form, extract the relevant information (wave speeds, jumps in the characteristic variables) for the (seven) physical waves, and analyze the cases of degeneracy (i.e., cases where several wave speeds corresponding to different waves become degenerate). Komissarov’s RMHD scheme is an extension of the scheme developed for RHD [89] described in Section 3.5, which avoids the use of the left eigenvectors (in [14] they are computed numerically). In its multi-dimensional version, Komissarov’s code enforces by employing the integral form of the induction equation. This code has been used to study the propagation of light, highly relativistic jets carrying toroidal magnetic fields [144].

Koide, Shibata, and Kudoh [139] performed simulations of magnetically driven axisymmetric jets from black hole accretion disks. Their GRMHD code [140] is an extension of the special relativistic MHD code developed by Koide et al. [138, 136, 210]. 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. The authors have recently extended their code to Kerr spacetimes [141] and performed simulations of axisymmetric jets formed by extracting rotational energy from a black hole [137, 142]. Finally, using a 3D GRMHD code, Nishikawa et al. [212] have investigated the dynamics of a freely falling corona and of a Keplerian accretion disk around a Schwarzschild black hole to form bipolar relativistic jets assuming axisymmetry as in previous simulations.

With the pioneering work of Koide and collaborators, numerical simulations have entered into the realm of GRMHD. However, despite their success, present simulations only cover a tiny fraction of dynamical time scales (about 2 rotational periods of the accretion disk) and jets are formed with very small terminal speeds (Lorentz factors less than 2). Hence, the quest for robust codes able to follow the formation of steady relativistic jets is still open. Given their success in SRHD, the extension of Riemann-solver based HRSC methods is an obvious option to bear in mind. Again, the third-order symmetric HRSC algorithms developed recently [72] represent a very interesting alternative.