General relativistic smoothed particle hydrodynamics

TL;DR

This paper introduces a novel general relativistic smoothed particle hydrodynamics (GRSPH) method that accurately models relativistic fluids around black holes, incorporating shock capturing, Hamiltonian integration, and benchmarking against known solutions.

Contribution

The paper develops a new GRSPH framework with improved shock handling and a Hamiltonian integrator, enabling more accurate simulations of relativistic astrophysical phenomena.

Findings

Successfully captures Einstein and spin precession effects.

Accurately reproduces shock tube and accretion flow solutions.

Demonstrates robustness in 1D and 3D relativistic tests.

Abstract

We present a method for general relativistic smoothed particle hydrodynamics (GRSPH), based on an entropy-conservative form of the general relativistic hydrodynamic equations for a perfect fluid. We aim to replace approximate treatments of general relativity in current SPH simulations of tidal disruption events and accretion discs. We develop an improved shock capturing formulation that distinguishes between shock viscosity and conductivity in relativity. We also describe a new Hamiltonian time integration algorithm for relativistic orbital dynamics and GRSPH. Our method correctly captures both Einstein and spin-induced precession around black holes. We benchmark our scheme in 1D and 3D against mildly and ultra relativistic shock tubes, exact solutions for epicyclic and vertical oscillation frequencies, and Bondi accretion. We assume fixed background metrics (Minkowski, Schwarzschild and Kerr in Cartesian Boyer-Lindquist coordinates) but the method lays the foundation for future direct coupling with numerical relativity.

Figures

Equations

$$U^{\mu }=\frac{\mathrm{d}{x}^{\mu }}{\mathrm{d}\tau },$$

$$v^{\mu }\equiv \frac{\mathrm{d}{x}^{\mu }}{\mathrm{d}t}=\frac{U^{\mu }}{U^0},$$

$$U^0\equiv \frac{\mathrm{d}t}{\mathrm{d}\tau }=\frac{1}{\sqrt{-g_{\mu \nu }{v}^{\mu }{v}^{\nu }}}.$$

$${\rho }^{\ast }$$

$$p_i$$

$$e$$

$$w=1+u+\frac{P}{\rho }.$$

$${\rho }^{\ast }$$

$$p_i$$

$$e$$

$$d{s}^2=-{\alpha }^2\mathrm{d}{t}^2+{\gamma }_{ij}(\mathrm{d}{x}^i+{\beta }^i\mathrm{d}t)(\mathrm{d}{x}^j+{\beta }^j\mathrm{d}t).$$

$$e=w\Gamma \left(\alpha -V^i{\beta }_i\right)-\frac{\alpha P}{\Gamma \rho }.$$

$$\frac{\mathrm{d}{\rho }^{\ast }}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{p}_i}{\mathrm{d}t}$$

$$\frac{\mathrm{d}e}{\mathrm{d}t}$$

$$f_i\equiv$$

$$\Lambda \equiv$$

$$T^{\mu \nu }$$

$$=\rho w(U^0{)}^2{v}^{\mu }{v}^{\nu }+P{g}^{\mu \nu }.$$

$$\frac{\sqrt{-g}{T}^{\mu \nu }}{2{\rho }^{\ast }}=\frac{1}{2}\left[w{U}^0{v}^{\mu }{v}^{\nu }+\frac{P{g}^{\mu \nu }}{\rho U^0}\right].$$

$$P=\left({\gamma }_{\mathrm{ad}}-1\right)\rho u,$$

$$P=\frac{\rho k_{\mathrm{B}}T}{\mu m_{\mathrm{H}}},$$

$$\frac{\mathrm{d}{\rho }_a^{\ast }}{\mathrm{d}t}=$$

$$\frac{\mathrm{d}{p}_i^a}{\mathrm{d}t}=$$

$$+\frac{\sqrt{-g_b}{P}_b}{{\Omega }_b{\rho }_b^{\ast 2}}\frac{\partial W_{ab}(h_b)}{\partial x^i}]+f_i^a,$$

$$\frac{\mathrm{d}{e}_a}{\mathrm{d}t}=$$

$$+\frac{\sqrt{-g_b}{P}_b{v}_a^i}{{\Omega }_b{\rho }_b^{\ast 2}}\frac{\partial W_{ab}(h_b)}{\partial x^i}]+{\Lambda }_a,$$

$${\Omega }_a=1-\frac{\partial h_a}{\partial {\rho }_a^{\ast }}\sum _b{m}_b\frac{\partial W_{ab}(h_a)}{\partial h_a}.$$

$${\rho }_a^{\ast }=\sum _b{m}_b{W}_{ab}(h_a).$$

$$h_a=h_{\mathrm{fac}}{\left(\frac{m_a}{{\rho }_a^{\ast }}\right)}^{1/d},$$

$$\frac{\partial h_a}{\partial {\rho }_a^{\ast }}=-\frac{1}{d}\frac{h_a}{{\rho }_a^{\ast }}.$$

$$W_{ab}(h_a)=W(h_a,L_{ab})\equiv \frac{C_{\mathrm{norm}}}{{h_a}^d}f(q),$$

$$f(q)=\left\{\begin{array}[]{ll}(3-q)^{5}-6(2-q)^{5}+15(1-q)^{5},&0\leq q<1,\\ (3-q)^{5}-6(2-q)^{5},&1\leq q<2,\\ (3-q)^{5},&2\leq q<3,\\ 0,&q\geq 3,\end{array}\right.$$

$$L_{ab}\equiv \sqrt{{\eta }_{ij}{r}_{ab}^i{r}_{ab}^j},$$

$$\frac{\partial W_{ab}(h_a)}{\partial x^i}=N_i{F}_{ab}(h_a),$$

$$N_i\equiv {\eta }_{ij}{\hat{r}}_{ab}^j,$$

Full text

††pagerange: General relativistic smoothed particle hydrodynamics–D††pubyear: 2019

1 Introduction

Urgent motivation for relativistic hydrodynamics simulations arises from the coincident detection of an electromagnetic afterglow (Abbott et al., 2017b) alongside the first detection of gravitational waves from a binary neutron star merger (Abbott et al., 2017a). The coming decade should bring hundreds more such events as the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) reaches design sensitivity (Abbott et al., 2017a).

Smoothed particle hydrodynamics (SPH) (Gingold & Monaghan, 1977; Lucy, 1977; Rosswog, 2009; Price, 2012) is perfectly suited for neutron star merger simulations because there is no preferred geometry, resolution follows mass, and it avoids the limitations imposed by a background density floor (e.g. Oechslin et al., 2002). The main limitation of SPH in relativistic hydrodynamics to date is the approximate treatment of general relativity (Metzger, 2017).

Further motivation for a public GRSPH code comes from i) the race to find electromagnetic counterparts to binary black hole mergers detected by Advanced LIGO and the forthcoming Laser Interferometer Space Antenna (e.g. Armitage & Natarajan, 2002; Milosavljević & Phinney, 2005; Cerioli et al., 2016); ii) disc formation in tidal disruption events (e.g. Hayasaki et al., 2013, 2016; Bonnerot et al., 2016; Tejeda et al., 2017); iii) the recently discovered phenomena of ‘tearing’ in warped discs around spinning black holes (e.g. Nixon et al., 2012; Nealon et al., 2015) and the possible relation to state transitions and quasi periodic oscillations (Nixon & Salvesen, 2014); iv) forthcoming observations from the Event Horizon Telescope (Akiyama et al., 2017); and v) the need for relativistic SPH simulations of pulsar winds (Okazaki et al., 2011).

Kheyfets et al. (1990) were the first to derive relativistic SPH equations, followed by Mann (1991) and Laguna et al. (1993). However, their formulations were not in conservative form and had difficulty handling shocks in even mildly relativistic flows. Conservative forms of the relativistic SPH equations were derived by Chow & Monaghan (1997) (hereafter CM97); Siegler & Riffert (2000); Monaghan & Price (2001); Rosswog (2010a) and Rosswog (2010b) for both special and general relativity, but their applications to date have been limited to special relativity. Other recent studies have instead used post-Newtonian approximations within standard SPH codes to account for relativistic effects in fluid flows around black holes (e.g. Tejeda & Rosswog, 2013; Nealon et al., 2015; Bonnerot et al., 2016; Hayasaki et al., 2016). SPH applications to fluid flows in strong self gravity assume conformal flatness (e.g. Oechslin et al., 2002; Faber et al., 2004; Bauswein et al., 2010) but these papers do not include standardised tests in fixed metrics. Most recently, Tejeda et al. (2017) applied GRSPH to tidal disruption events, but also showed only limited standardised tests of their method.

Our goal in this paper is to present a detailed description of a fully general relativistic SPH code, able to perform hydrodynamic simulations in any given fixed background metric, with the potential to be applied to dynamically evolving metrics. Our requirements for a modern GRSPH code include being able to capture relativistic shocks, a precise and accurate treatment of orbital dynamics, and the ability to work in Kerr geometry. To this end we derive an improved shock capturing method and a new Hamiltonian time integration algorithm suitable for GRSPH and relativistic orbital dynamics. We describe our complete method and its implementation in depth, before showing the results of various benchmark tests. We implement our method in Phantom (Price et al., 2018), a publicly available SPH code.

2 Relativistic hydrodynamics

We follow Monaghan & Price (2001) who obtained the equations of general relativistic smoothed particle hydrodynamics from a Lagrangian. In the following, Greek indices are summed over ([math],$1$,$2$,$3$) while Latin indices are summed over ($1$,$2$,$3$). The subscripts $a$ and $b$ are reserved for particle labels.

2.1 Definitions

The four velocity of a fluid at a given spacetime coordinate $x^{\mu}$ is given by

[TABLE]

where $\tau$ is the proper time, and $x^{\mu}\equiv(t,x^{i})$. The coordinate velocity $v^{\mu}$ is defined with respect to the coordinate time $t$ such that

[TABLE]

where from the normalisation condition $U^{\mu}U_{\mu}=-1$ we have

[TABLE]

We assume standard relativistic units where $c=G=1$ and a metric signature ($-$,$+$,$+$,$+$).

2.2 Conserved variables

The conserved density, momentum and energy are respectively

[TABLE]

where $g$ is the determinant of the covariant metric $g_{\mu\nu}$ and we use $w$ to denote the specific enthalpy

[TABLE]

The specific internal energy, pressure, and density in the rest frame of the fluid are denoted by $u,P$ and $\rho$, respectively. Equivalently, given a $3+1$ decomposition of the metric (Arnowitt et al., 2008) we can express the conserved quantities as

[TABLE]

where $\Gamma=\left(1-V^{i}V_{i}\right)^{-1/2}$ is the generalised Lorentz factor, and $V^{i}=\left(v^{i}+\beta^{i}\right)/\alpha$ is the fluid velocity in the frame of a local Eulerian observer. The lapse function, shift vector, and spatial three-metric ($\alpha$, $\beta_{i}$, $\gamma_{ij}$) are related to the four-metric $ds^{2}=g_{\mu\nu}\mathrm{d}x^{\mu}\mathrm{d}x^{\nu}$ via

[TABLE]

Spatial vectors can then have their indices raised and lowered by the three-metric, e.g. $v^{i}v_{i}=\gamma_{ij}v^{i}v^{j}=\gamma^{ij}v_{i}v_{j}$. We note that $\gamma$ is the determinant of $\gamma_{ij}$, and that Equation 2.10 can be written alternatively in terms of the enthalpy and pressure as

[TABLE]

2.3 Equations of relativistic hydrodynamics

From Monaghan & Price (2001), the equations of relativistic hydrodynamics in Lagrangian conservative form, assuming an ideal fluid, are given by

[TABLE]

where the source terms containing derivatives of the metric are defined according to

[TABLE]

For an ideal fluid the stress-energy tensor is given by

[TABLE]

Thus, the term required to compute the source terms is

[TABLE]

We close the equation set with the ideal gas equation of state

[TABLE]

where ${\gamma_{\mathrm{ad}}}$ is the adiabatic index. For our calculations we assume an ideal atomic gas, ${\gamma_{\mathrm{ad}}}=5/3$. The internal energy $u$ can be related to the gas temperature $T$ through the ideal gas law

[TABLE]

where $k_{\mathrm{B}}$ is Boltzmann’s constant, $\mu$ is the mean molecular weight, and $m_{\mathrm{H}}$ is the mass of a Hydrogen atom.

3 SPH Equations of relativistic hydrodynamics

The discrete form of Equations 2.13, 2.14, and 2.15 in the absence of dissipation, are given respectively by (e.g. Siegler & Riffert, 2000; Monaghan & Price, 2001; Rosswog, 2010b)

[TABLE]

where $W_{ab}$ is the interpolating kernel, $h_{a}$ is the smoothing length, and $\Omega_{a}$ is a term related to the gradient of the smoothing length given by (Monaghan, 2002; Springel & Hernquist, 2002)

[TABLE]

The conserved density can also be computed directly via

[TABLE]

We choose to compute $\rho^{*}$ by summation, rather than by integration of Equation 3.1. Given equal mass particles, the smoothing length is adapted by simultaneously solving 3.5 and

[TABLE]

where $d$ is the number of dimensions, and $h_{\mathrm{fac}}$ is a numerical parameter specifying the smoothing length in terms of the mean particle spacing. The term required in 3.4 is then

[TABLE]

For details on solving 3.5 and 3.6 simultaneously we refer the reader to Price & Monaghan (2007).

3.1 Interpolating kernel

We use a kernel in the form

[TABLE]

where $L_{ab}$ is the inter-particle spacing, $C_{\mathrm{norm}}$ is a normalisation constant, and $q=L_{ab}/h_{a}$. The standard choice for the dimensionless function $f(q)$ in SPH is the cubic spline. With $h_{\mathrm{fac}}=1.2$ this results in 57.9 neighbours per particle on average in 3D. In general, more accurate results can be obtained by using a smoother spline such as the quintic, given by

[TABLE]

for which $C_{\mathrm{norm}}=[1/120,1/(120\pi)]$ in 1D and 3D. This comes at the cost of a greater number of neighbours per particle; 113 on average in 3D when $h_{\mathrm{fac}}=1$. For consistency, we choose to use the quintic spline everywhere, however the cubic spline gives satisfactory results also. The neighbour finding procedure in more than one dimension is non trivial. We employ a $k$d-tree for this process, the details of which can be found in the Phantom paper, Price et al. (2018) (see their section 2.1.7).

For the inter-particle spacing $L_{ab}$ we assume local flatness and thus use the Euclidean distance. We define this distance in Cartesian coordinates as

[TABLE]

where $\eta_{ij}$ is the spatial part of the Minkowski metric, and $r^{i}_{ab}\equiv x^{i}_{a}-x^{i}_{b}$ is the line of sight vector pointing from particle $a$ to particle $b$ also in Cartesian coordinates. Computing the kernel and its derivative is then identical to what is done in standard non-relativistic SPH. In Cartesian coordinates, since $\eta_{ij}$ is a constant, we can write the kernel derivative as

[TABLE]

where $F_{ab}\leq 0$ is the scalar part of the kernel gradient. We define the line of sight projection operator,

[TABLE]

where,

[TABLE]

is the line of sight unit vector.

We emphasise that the choice between using the Euclidean distance norm and proper lengths for neighbour distances is arbitrary. Using a different distance metric would simply require a different particle placement to obtain the same conserved density. So this choice reflects how one wishes to place particles to resolve features in the spacetime. We argue that the Euclidean norm is the most natural choice because the particle placement is independent of the metric. For example, a uniformly spaced lattice of particles produces a constant $\rho^{*}$. Using the proper length for neighbour distances would require stretching or squeezing the initial particle arrangement to compensate for the determinant of the metric in the volume element. It would also require computing a path integral in the given metric between all neighbouring particles.

One possible advantage to using the proper distance would be that smoothing spheres would never cross the event horizon. That is, neighbours close to the event horizon would become infinitely distant. However, it is no longer obvious how to set up the particles to achieve a particular $\rho^{*}$ in this case, and using different metrics would require different particle setups. We demonstrate in Section 8.4 that using the Euclidean distance norm indeed produces the correct conserved densities close to the black hole so long as the horizon is sufficiently well resolved.

4 Shock capturing

It is necessary to introduce dissipative terms in order to capture shocks. As in CM97 and Monaghan (1997) the construction of the dissipative terms is motivated by approximate Riemann solvers used by finite volume codes (Martí et al., 1991). To demonstrate this, consider the equation set written in conservative form and in one spatial dimension (for simplicity)

[TABLE]

where $\mathbf{u}=(\rho^{*},p_{i},e)$ is the vector of conservative variables and $\mathbf{F}(\mathbf{u})$ is the flux. Finite volume schemes require computing the numerical flux between adjacent states $\mathbf{u}_{\mathrm{L}},\mathbf{u}_{\mathrm{R}}$. The simplest approach utilises the ‘local Lax-Friedrichs’ or ‘Rusanov’ flux given by

[TABLE]

where $v_{\mathrm{sig}}$ is the maximum characteristic speed. Importantly, the dissipative term (proportional to $v_{\mathrm{sig}}$) acts on jumps in the conservative variables.

Although the exact relationship between the dissipation introduced in approximate Riemann solvers and physical dissipation terms is unclear in relativity, since the correct way to formulate viscosity and heat conduction in relativistic hydrodynamics is not well understood (e.g. Andersson & Comer, 2007), the approach used above is standard in GRMHD codes (e.g. Martí et al., 1991; Gammie et al., 2003; Zhang & MacFadyen, 2006).

Our SPH shock capturing formulation is based on approximate Riemann solvers, following CM97. If one considers adjacent particles in SPH as left and right states, an appropriate form of the dissipative terms requires a jump in the conservative variables and a maximum signal speed $v_{\mathrm{sig}}$ between the particles. CM97 proposed corresponding terms in the momentum and energy equations in SPH given by

[TABLE]

where $p^{*}$ and $e^{*}$ are computed using only velocities along the line of sight ($v^{*}$), $\bar{\rho}^{*}_{ab}=(\rho_{a}^{*}+\rho_{b}^{*})/2$, and $\alpha_{\mathrm{AV}}\simeq 1$ is a numerical parameter to control the amount of dissipation.

The problem with the formulation given by CM97 is that, while they proved that the entropy increase was positive definite for cold gas ($w=1$), we found this to no longer be the case when the enthalpy differs from unity. To remedy this, we adapted their formulation to not use the full jump in momentum or total energy

[TABLE]

where starred quantities again denote the use of the line of sight velocity $v^{*}$. Instead, in order to have a positive definite change in entropy we use a jump only in the ‘kinetic terms’, $\Gamma_{a}^{*}V_{a}^{*}-\Gamma_{b}^{*}V_{b}^{*}$. For the momentum equation we multiply the corresponding term by an average enthalpy, and for the energy equation we multiply by an average of the term $wv$, in order to recover the correct dimensions. That is, we define

[TABLE]

where $\overline{w}_{ab}$ denotes a suitable average of $w$ over particles $a$ and $b$.

A second problem with the CM97 formulation was that it does not distinguish between viscosity and conductivity. This leads to over-smoothing of contact discontinuities (see Section 8.1.1). Showing how to split the dissipation into viscosity and conductivity in relativity is non-trivial. We demonstrate in Section 8.1.1 and Appendix A.1 that the relevant splitting can be achieved by defining the energy jump in two parts as follows

[TABLE]

where $\bar{v}^{*}_{ab}$ denotes an average of $v^{*}$ over particles $a$ and $b$ (see below). This is equivalent to defining a relativistic thermal energy

[TABLE]

Our third change is more cosmetic. Rather than using separate averages for $\overline{w}_{ab}$ and $\bar{v}_{ab}^{*}$, we average them along with the density, signal speed and kernel derivative in the manner of Price & Federrath (2010) following the non-relativistic terms in Phantom (Price et al., 2018). That is, we write the viscosity as a $P+q$ term in the momentum and energy equations.

4.1 Artificial viscosity

Equations 3.2 and 3.3 with our dissipation become

[TABLE]

where $q_{a}$ and $q_{b}$ encapsulate the artificial viscosity, and $\Pi^{a}_{\mathrm{cond}}$ the artificial conductivity. For added simplicity in notation, we have defined

[TABLE]

Our artificial viscosity terms are

[TABLE]

and involve a fluid velocity that is along the line of sight between particles, similar to CM97, i.e.

[TABLE]

If the particles are separating, $N_{i}(v_{a}^{i}-v_{b}^{i})\geq 0$, we set $q_{a}=q_{b}=0$.

4.2 Artificial conductivity

Our artificial conductivity term is given by

[TABLE]

where again for simplicity we define

[TABLE]

which in Cartesian coordinates is

[TABLE]

Note that the terms $\alpha_{\mathrm{AV}}$ and $\alpha_{u}$ allow us to control the amount of artificial viscosity and conductivity, respectively. These should not be confused with $\alpha_{a}$ and $\alpha_{b}$ which are the lapse function at the positions of particles $a$ and $b$.

4.3 Signal speed — viscosity

When constructing an approximation for the maximum signal speed in the artificial viscosity, $v_{\mathrm{sig}}$, we impose two requirements. First, the signal speed must not exceed the speed of light. Second, the signal speed should reduce to the form used in standard SPH codes in the non-relativistic limit, namely

[TABLE]

where $c_{\mathrm{s}}$ is the sound speed, $\mathbf{v}_{ab}=\mathbf{v}_{a}-\mathbf{v}_{b}$ is the relative velocity, and $\mathbf{\hat{r}}_{ab}$ is the unit vector along the line of sight. The numerical parameters $[\alpha_{\mathrm{AV}},\beta_{\mathrm{AV}}]$ control the amount of dissipation. Most importantly, the signal speed should involve a relative velocity between the particles to avoid excessive dissipation when there is a bulk velocity and individual particle speeds are high.

Because of these requirements, we do not use the signal speeds proposed by CM97, since in some of their formulations the maximum signal speed may exceed the speed of light. We also do not use the signal speed proposed by Rosswog (2010b), based on the eigenvalues of the Euler equations. Their formulation, although limited by the speed of light, involves the absolute velocity of the SPH particles. This means their viscosity term is no longer Galilean invariant. We therefore require a signal speed that uses only the relative velocity between particles. To this end we propose a signal speed given by

[TABLE]

where the relative velocity between particles, considered along the line of sight, is defined as

[TABLE]

and the relativistic sound speed is given by

[TABLE]

We use the velocities in the frame of the local Eulerian observer $V^{i}$ rather than the coordinate velocities $v^{i}$, since only $V^{i}$ is a physical velocity that is limited by the speed of light, that is $|V_{i}V^{i}|\leq 1$. The coordinate velocity on the other hand, is not strictly under the same restriction.

4.4 Signal speed — conductivity

For simplicity one may assume the same signal speed for both viscosity and conductivity. That is, $v^{u}_{\mathrm{sig}}=v_{\mathrm{sig}}$. In general however, one should use different speeds since the conductivity is mainly required at contact discontinuities which travel at a different speed to the shock. Price (2008) proposed using a signal speed based on the relative pressure jump, designed to eliminate the pressure blip at contact discontinuities. We generalise this to the relativistic case using

[TABLE]

For cases where equilibrium pressure gradients exist due to other physics (e.g. gravity) it is better to simply use the relative speeds (Price, 2008; Wadsley et al., 2008), giving

[TABLE]

when generalised to the relativistic case. In the non-relativistic limit, both eq. 4.25 and 4.26 reduce to the non-relativistic implementation in Phantom (Price et al., 2018), where the choice is based on whether or not gravity and/or external forces are used.

The main advantage of the above signal speeds is that conductivity becomes $\mathcal{O}(h^{2})$ in the numerical code, rather than first order. For most of the tests in this paper, we simply use $v^{u}_{\mathrm{sig}}=v_{\mathrm{sig}}$ with a reduced value of $\alpha_{u}$. In section 8.1.6 we demonstrate the benefit of using Eq. 4.26 in the case of a self-gravitating, oscillating polytrope.

5 Evolving entropy

Instead of the total specific energy $e$, we evolve an entropy variable (cf. Springel & Hernquist 2002) in order to guarantee positivity of the thermal energy. We define

[TABLE]

Taking the time derivative and using Equation 2.21 we find

[TABLE]

Comparing this to the second law of thermodynamics

[TABLE]

we see that in the absence of dissipation

[TABLE]

Furthermore, using Equation 2.22, $K$ can be related to the specific entropy $s$

[TABLE]

More generally, we can rewrite Equation 5.3 in terms of time derivatives of the conserved variables ($\rho^{*}$, $p_{i}$, $e$), giving

[TABLE]

which, after substituting Equations 3.1, 4.10, and 4.11, gives

[TABLE]

Thus, the evolution of $K$ is governed purely by the dissipative terms. An important distinction between SPH and finite-volume schemes is that when evolving $K$, energy remains exactly conserved (to the accuracy of the time integration algorithm), since Eq 5.7 holds to round-off error. In Appendix A we demonstrate that Eq. 5.8 is always positive (for special relativity) thus ensuring that entropy always increases.

6 Time Integration

Relativistic hydrodynamics, in the absence of dissipation, forms a Hamiltonian system where the canonical variables are the positions $x^{i}$ and momenta $p_{i}$. Ideally therefore one should use a geometric integrator that preserves the Hamiltonian properties. The main complication is that the Hamiltonian is non-separable — it cannot be separated into kinetic, thermal, and potential energy terms.

6.1 Generalised leapfrog method

Leimkuhler & Reich (2005) suggest the ‘generalised leapfrog method’ for the symplectic integration of Hamilton’s equations with a non-separable Hamiltonian. With $p_{i}$ and $x^{i}$ as the conjugate variables, their scheme is given by

[TABLE]

where the superscripts $n+m$ refer to quantities evaluated at times $t^{n+m}=t^{n}+m\Delta t$. It may be seen that this reduces to the standard Kick-Drift-Kick leapfrog scheme when $v^{i}\equiv p_{i}$ and when the acceleration ${\mathrm{d}}p_{i}/{\mathrm{d}}t$ depends only on position. It is also time-reversible and symplectic. However, a general property of geometric integrators for non-separable Hamiltonians is that they are implicit in one or more of the steps (e.g. Yoshida, 1993). In the above integrator this occurs in the first ‘Kick’ step, where the momentum at the half step appears on both sides of the equation, and in the ‘Drift’ step where the position at the full step appears on both sides.

6.2 Alternative generalised leapfrog method

For GRSPH, evaluating ${\mathrm{d}}p_{i}/{\mathrm{d}}t(p_{i},x^{i})$ means performing a sum over particles, while evaluation of ${{\mathrm{d}}x^{i}}/{{\mathrm{d}}t}(p_{i},x^{i})$ means performing a primitive variable solve (see Section 7) which can be performed independently for each particle. Hence the previous integrator would require two calls to the summation routines, since ${\mathrm{d}}p_{i}/{\mathrm{d}}t$ used in the second Kick step is different from the one required for the subsequent first Kick step at the start of the next timestep. For this reason we propose a similar integrator where the symmetrisation occurs in the drift step. Our alternative integrator is given by

[TABLE]

This integrator retains the time symmetry of the previous integrator, however it does not exactly conserve angular momentum. The difference is that the implicit steps are now in the Drift step, as previously, and in the second Kick step. Importantly, the derivative of the momentum used in the second Kick step is the same as that used at the beginning of the subsequent step, meaning that it can be re-used. The integrator is also easier to start, because it merely requires an evaluation of the momentum derivative involving position and momentum from the initial conditions.

The implicit solve in the Drift step (in the position derivative) involves only the primitive variable solver, so is independent for each particle. The implicit nature of the second Kick step (in the momentum derivative) in principle involves multiple loops over neighbouring particles, but this is similar to our existing approach to dealing with viscosity and other velocity-dependent terms in non-relativistic SPH (see below). With a good prediction of $p^{n+1}_{i}$, only one call to the summation routines is required.

6.3 Hybrid leapfrog, and implementation in GRSPH

Since our modified integrator does not exactly conserve angular momentum, we implement a hybrid of the two algorithms in order to retain the benefits of both. Taking an operator splitting approach similar to that of the reversible reference system propagator algorithm (RESPA) derived by Tuckerman et al. (1992) we split the acceleration into ‘long range’ and ‘short range’ forces (see Price et al. 2018 for an implementation in SPH). In our case the ‘long range’ forces are the SPH forces, while the ‘short range’ forces are the external forces arising from the curvature of the spacetime. These correspond to the pressure gradient term and $f_{i}$ in the momentum equation, respectively. i.e.

[TABLE]

We use the ‘long range’ forces $f_{i}^{\,\mathrm{sph}}$ to update the momentum of each particle and use the ‘short range’ forces $f_{i}^{\,\mathrm{ext}}$ in updating their positions. Our complete hybrid algorithm is therefore

[TABLE]

[TABLE]

[TABLE]

where the position update is done over $m$ substeps (indicated by the equations in braces), such that $m=\mathrm{int}(\Delta t_{\mathrm{ext}}/\Delta t_{\mathrm{sph}}+1)$ and is the minimum over all particles. We denote the updating of a variable by an arrow, some of which first require solving and implicit equation, namely Equations 6.10, 6.12, and 6.15.

The combination of the two algorithms ensures that one of the SPH force calculations can be reused between steps, reducing computational cost, and that angular momentum is conserved in the case where the SPH forces can be neglected. More importantly it ensures that our method is not restricted by the external force timestep, since computing $f_{i}^{\,\mathrm{ext}}$ is cheap whereas $f_{i}^{\,\mathrm{sph}}$ is expensive.

It is useful to point out that in the absence of metric gradients the hybrid method reduces to our proposed scheme in section 6.2, while for pressureless fluid it reduces to the method of Leimkuhler & Reich (2005) in section 6.1.

6.3.1 Evaluating implicit steps

To evaluate Equations 6.10 and 6.12, we employ two first-order predictions which we refer to as $p^{*}_{i}$ and $x^{i,*}$. That is, for Eq 6.10 we do

[TABLE]

and for evaluating Eq 6.12 we do

[TABLE]

To complete the implicit step we check that

[TABLE]

after Equations 6.18 and 6.20, respectively, where $\epsilon_{x}$ and $\epsilon_{p}$ are small tolerances. If these criteria are violated we iterate Equation 6.18 (or Eq. 6.20) by a fresh evaluation of the corresponding derivative with the updated value, that is, simply replacing $p_{i}^{*}$ with $\tilde{p}_{i}$, and $p_{i}$ with $p_{i}^{*}$ on the right hand side of Eq. 6.18 (and similarly with $x^{i,*}$ and $x^{i}$ in Eq. 6.20).

Equation 6.15 is evaluated in the same manner as 6.10 but with $f_{i}^{\,\mathrm{sph}}$ and $\Delta t_{\mathrm{sph}}$ instead of $f_{i}^{\,\mathrm{ext}}$ and $\Delta t_{\mathrm{ext}}$. We use the error criterion

[TABLE]

where we define $|p|\equiv\sqrt{p_{i}p_{j}\eta^{ij}}$ and $\Delta p_{i}\equiv\tilde{p}_{i}-p_{i}^{*}$, consistent with the implementation of the equivalent implicit step in Phantom (cf. Price et al., 2018). The maximum is taken over all the particles in the numerator, while the denominator evaluates the root mean square, also over all the particles.

In general, iterations of Equations 6.10 and 6.12 are cheap, while iterations of Equation 6.15 are expensive, so $\epsilon_{p}$ and $\epsilon_{x}$ can be set arbitrarily small while $\epsilon_{p}^{\mathrm{sph}}$ must allow for an efficient code. By default the tolerance is set to $\epsilon_{p}^{\mathrm{sph}}=10^{-2}$ in Phantom. To try to prevent iterations from occurring, we constrain the timestep by

[TABLE]

since an iteration is as expensive as halving the timestep.

In terms of storage, we require only one position and momenta for each particle, but storage of two derivatives (current and previous) of both $v^{i}$ and $f_{i}$.

6.4 Timestep constraint

We adopt the usual stability conditions, namely

[TABLE]

where $C_{\mathrm{C}}$ and $C_{\mathrm{f}}$ are the usual safety factors, set to 0.3 and 0.25, respectively by default in Phantom, the maximum $v_{\mathrm{sig}}$ is over neighbours to particle $a$, and $|f|\equiv\sqrt{f_{i}f_{j}\eta^{ij}}$. In our tests we use global timestepping where the global timestep $\Delta t_{\mathrm{sph}}$ is set to the minimum value of $\Delta t_{a}^{\mathrm{sph}}$ on the particles, although our scheme is simple to generalise to the case where particles have individual timesteps. The timestep constraint for external force substeps is similarly

[TABLE]

and the global $\Delta t_{\mathrm{ext}}$ is the minimum of $\Delta t_{a}^{\mathrm{ext}}$ on the particles.

7 Recovery of primitive variables

Computing the right hand side of Equations 3.1, 3.2, and 3.3 explicitly requires the primitive variables $(\rho,u,v^{i},P)$. We must therefore invert the set of equations describing the conserved variables (Eq. 2.8, 2.9, and 2.10) at every timestep. To do this we use the recovery scheme given by Tejeda (2012).

Use the equation of state to solve for the enthalpy as a function of density and pressure only,

[TABLE] 2. 2.

Express $\rho$, $P$ and $v^{i}$ as functions of only $w$ and the conserved variables,

[TABLE]

where $v_{i}=\gamma_{ij}~{}v^{j}$, and $\Gamma(w)$ can be calculated from Eq. 2.9 as,

[TABLE] 3. 3.

Given Equations 7.2 and 7.3, solve Equation 7.1 for the enthalpy using a standard root-finding algorithm. (See Appendix B for more detail). 4. 4.

Compute $v_{i}$ using Equation 7.4, and raise the index using the three-metric $v^{i}=\gamma^{ij}v_{j}$.

8 Numerical tests

We split the testing of our numerical method into three stages.

First, we developed a dedicated 1D special relativistic SPH code to test the shock capturing terms in section 4, and the conservative-to-primitive algorithm in section 7. For the time integration in this code we used the standard 2nd order Runge-Kutta method (RK2), with time steps constrained by Equation 6.25). We use RK2 instead of the algorithm described in section 6 to show that the particular timestepping algorithm used is not important for shock tube tests. 2. 2.

Next, we developed a 3D N-body code to test the force terms involving the metric derivatives $f_{i}$. This involved integrating Equation 2.14 with $P=0$, i.e.

[TABLE]

Although this is equivalent to solving the geodesic equation, it allowed us to test our time integration scheme (section 6.3), as well as the conservative-to-primitive algorithm, again, but with more complicated metrics than the Minkowski metric. We implemented the Schwarzschild metric in both spherical and Cartesian coordinates, along with the Kerr metric in Boyer-Lindquist coordinates (both spherical and Cartesian). For this code we used a fixed time step. 3. 3.

Finally, after testing individual aspects of our method in the previous two codes, we incorporated them all into Phantom, a fully 3-dimensional SPH code (Price et al., 2018). In this code, we implemented the full, hybrid-Leapfrog method for the time integration, as described in section 6.3. The timestep is constrained according to Equations 6.24, 6.25, & 6.26.

For SPH simulations we always use equal mass particles, and assume an adiabatic gas with ${\gamma_{\mathrm{ad}}}=\nicefrac{{5}}{{3}}$. For the artificial dissipation we use $\alpha_{\mathrm{AV}}=1$ and $\alpha_{u}=0.1$, unless otherwise stated.

We quantify the error in our simulations with the dimensionless $L_{2}$ error

[TABLE]

where $N$ is the number of points and $y_{\mathrm{exact}}$ is the exact solution at point $i$.

8.1 Special Relativity

8.1.1 Mildly relativistic shock

In Figures 2 & 2 we show the results for a mildly relativistic shock in both 1D and 3D, respectively. This is ‘Problem 1’ in Martí & Müller (2003) with a maximum fluid Lorentz factor of $\Gamma=1.38$. For the initial conditions, we set

[TABLE]

The particle spacing to the left of the interface is 0.0005, and 0.005 to the right. For the initial particle arrangement in 3D, we follow Price et al. (2018). We employ a uniform close-packed lattice of $1000\times 26\times 26$ particles for $x\in[-0.5,0]$ and $434\times 12\times 12$ particles for $x\in[0,0.5]$, with periodic boundaries in the $y$ and $z$ directions. Of the four hydrodynamic variables plotted ($\rho^{*}$, $v^{x}$, $u$, $P$), the error is largest in $v^{x}$ with $L_{2}=4.0\times 10^{-2}$ in 1D, and $L_{2}=3.2\times 10^{-2}$ in 3D.

The most noticeable difference between the 1D and 3D simulations is the slight overshoot at the left end of the rarefaction fan $x\approx-0.2$. This is caused by the lack of dissipation applied to the small discontinuity at the end of the rarefaction fan. It can be eliminated in 1D by applying the viscosity term for both approaching and receding particles. In 3D, even with viscosity applied only to approaching particles there is some dissipation applied at this location because the velocity differences are computed more isotropically. The small excess velocity at $x\approx 0.25$ seen in 3D is caused by particles rearranging behind the shock front.

Figure 3 shows the same shock but with $\alpha_{u}=0$. Without artificial conductivity we find a ‘blip’ in the pressure profile at the contact discontinuity. The same feature is discussed by Price (2008) and Price (2012) in the context of non-relativistic SPH simulations. In both cases, the ‘blip’ can be smoothed out by including artificial conductivity. Simulations by CM97 do not contain such a ‘blip’, but instead display excessive thermal diffusion at the contact discontinuity, since artificial conductivity is included implicitly within their dissipation terms (compare our results to figures 9 and 10 in their paper).

8.1.2 Strong blast

In Figures 5 & 5 we show the results for a strong blast wave in both 1D and 3D, respectively. This is ‘Problem 2’ in Martí & Müller (2003) with a maximum fluid Lorentz factor of $\Gamma=3.6$. We follow the same setup as for the mildly-relativistic shock but with differing initial conditions

[TABLE]

The particle spacing is 0.001 on both sides of the interface, with $1000\times 12\times 12$ particles for the whole domain in 3D. Importantly, we use unsmoothed initial conditions, unlike in CM97 and Rosswog (2010b). We found this was only possible when we evolved the entropy rather than total energy as the conserved variable, otherwise negative pressures could result in the absence of smoothing. The error is again largest in $v^{x}$ for both cases, with $L_{2}=2.5\times 10^{-1}$ in 1D and $L_{2}=3.1\times 10^{-1}$ in 3D. The peak of the density spike is also within 2.7% of the exact solution in both the 1D and 3D simulations. The same differences between 1D and 3D simulations occur as in the previous test, and for the same reasons.

Figure 6 shows the details of the density spike at different resolutions. At low resolutions the shock speed is overestimated, consistent with the results of CM97. With 500 particles the centre of the spike is ahead of the exact solution by approximately one spike width. As the resolution is increased, our simulations approach the exact solution.

8.1.3 Sinusoidally perturbed shock

In Figure 7 we show the results in 1D and 3D for a shock tube with a sinusoidally perturbed right density state (e.g. Del Zanna & Bucciantini, 2002; Rosswog, 2010b) where the fluid reaches a maximum Lorentz factor of $\Gamma=1.15$. The initial conditions for this setup are

[TABLE]

The particle spacing is 0.0004 on the left side of the interface, and on the right such that the particle masses are the same as on the left. We find that the sine wave perturbation is easily transported across the shock front in both cases, with no noticeable loss in accuracy. The peaks of the propagated sine wave reach the correct levels of the two limiting solutions; to within 1.2% in 1D, and within 1.3% in 3D. The left end of the rarefaction fan in 1D displays an overshoot in the density, for the same reasons as discussed in the previous tests.

8.1.4 Shocks with transverse velocity

We consider two problems where the shock contains a non-zero transverse velocity to the main flow. This is important in relativistic flows since the magnitude of the flow speed is limited by the speed of light, and the transverse velocity is coupled to the evolution of the other hydrodynamic variables through the Lorentz factor.

Figure 9 shows our results in 1D for a shocktube with initial conditions given by

[TABLE]

The initial particle spacing is 0.0001 on the left side of the interface. Of the four hydrodynamic variables plotted ($\rho$, $v^{x}$, $v^{y}$, $P$), the error is greatest in $\rho$ with $L_{2}=3.7\times 10^{-2}$, and a $5\%$ overshoot at the contact discontinuity. The fluid has a maximum Lorentz factor of $\Gamma=7.09$.

A much more challenging test is the shocktube with initial conditions given by

[TABLE]

As pointed out by Zhang & MacFadyen (2006) this tests requires a ‘very high resolution to resolve the complicated structure of the transverse velocity’. Figure 9 shows the results for this problem in 1D. We find that this test is particularly challenging since the transverse velocity changes rapidly, requiring very high spatial resolution. The most obvious artefact is a spike in $v^{y}$ at $x\approx 0.2$. This is also observed in Zhang & MacFadyen (2006) but is worse in SPH since this feature occurs in a low density part of the simulation, where the resolution is lowest. The largest error is in $\rho$ with $L_{2}=4.0\times 10^{-1}$. The fluid reaches a maximum Lorentz factor of $\Gamma=44.5$.

8.1.5 3D spherical blast wave

In Figure 10 we show the results of a spherically symmetric relativistic blast in 3D, analogous to the non-relativistic Sedov-Taylor blast wave. This is the same problem as done by Del Zanna & Bucciantini (2002) and Zhang & MacFadyen (2006) (see their Figures 8 and 17, respectively).

For the initial conditions we set $640\times 740\times 784$ particles on a closepacked lattice within a uniform periodic box $x,y,z\in[-1,1]$. The density is initially set to $\rho^{*}=\rho=1$ everywhere. We inject a sphere of uniform pressure $P_{\mathrm{in}}=1000$ at the origin with radius $r_{0}=0.4$, and set the pressure outside the sphere to $P_{\mathrm{out}}=0$. Our setup and resolution is equivalent to Zhang & MacFadyen (2006), however we simulate the whole box $x,y,z\in[-1,1]$, instead of just the quadrant $x,y,z\in[0,1]$ which is only possible in spherical symmetry.

Following Chow & Monaghan (1997) and Rosswog (2010b) we slightly smoothed the initial discontinuity in the pressure using

[TABLE]

where for the characteristic transition length $\Delta r$ we used half the initial particle spacing.

We compare our results to the 1D high resolution results of Zhang & MacFadyen (2006). Our solution shows qualitative agreement with theirs, with the most notable difference being an $\sim 8\%$ increase in the speed of our shock front. There is also a small amount of smoothing at the left end of the rarefaction fan ($x\approx 0.05$), which is present in the 3D simulations of Del Zanna & Bucciantini (2002) and Zhang & MacFadyen (2006). Figure 11 shows the (primitive) density slice through the $z=0$ plane, demonstrating our code’s ability to preserve radial symmetry and resolve the structure of the thin shell. The fluid reaches a maximum Lorentz factor of $\Gamma=19.56$.

8.1.6 Oscillating polytrope

A concern with artificial conductivity (expressed by the referee of this paper) is that it may introduce unwanted side effects. For example, spuriously applied conductivity may seriously disrupt a stellar equilibrium state. To address this concern, Figure 12 shows the results of a 3D oscillating polytropic star with Newtonian self-gravity in the Minkowski metric. The initial conditions are a relaxed star consisting of $65752$ $\left(\approx 4\pi 25^{3}/3\right)$ particles and polytropic equation of state $P(r)=K\rho(r)^{(n+1)/n}$ with $K=10^{-10}$ and $n=1.5$ (i.e. ${\gamma_{\mathrm{ad}}}=5/3$).

To relax the star we apply an additional force to the particles in the momentum equation which depends on the particle velocity (Gingold & Monaghan, 1977)

[TABLE]

and turn off all dissipative terms. That is, we set $\alpha_{\mathrm{AV}}=0$ and $\alpha_{u}=0$. We integrate until $t=20000$ when all particle motions have been sufficiently damped.

To induce radial oscillations in the star, we perturb the relaxed star profile with a small velocity perturbation

[TABLE]

Figure 12 compares the results of three different simulations, one with no dissipation, a second with conductivity turned on using $v^{u}_{\mathrm{sig}}=v_{\mathrm{sig}}$, and a third using conductivity with $v^{u}_{\mathrm{sig}}$ as in Eq. 4.26.

With no dissipation the amplitude of the stellar oscillation does not decay with time since energy is conserved exactly (black line in both panels). However with $v^{u}_{\mathrm{sig}}=v_{\mathrm{sig}}$, the stellar profile slowly becomes disturbed, and the central density increases. The results of the third simulation demonstrate a way to remedy this issue. The signal speed in Eq. 4.26 gives conductivity of $\mathcal{O}(h^{2})$, and which vanishes when the star is in hydrostatic equilibrium. This allows us to reproduce the same undamped oscillations as in the dissipation-free simulation.

8.2 Orbital dynamics — Schwarzschild Metric

In order to test our ability to perform orbital dynamics accurately, we use the Schwarzschild metric in spherical Schwarzschild coordinates. We confirmed that the same results can be reproduced using the equivalent Cartesian representation (see Appendix C.1). We use a black hole mass of $M=1$ in all simulations.

8.2.1 Circular orbits

The orbital frequency $\Omega$ of a circular geodesic in the Schwarzschild metric is identical to the frequency of a Newtonian orbit

[TABLE]

Figure 13 (left) shows a test particle orbit with initial position $(r,\theta,\phi)=(10M,\pi/2,0)$ and initial velocity $(v^{r},v^{\theta},v^{\phi})=(0,0,\Omega)$, simulated for $15$ orbital periods with constant time step $\Delta t=0.01M$ throughout. Energy and angular momentum $(e,p_{\phi})$ are conserved to machine precision ($\Delta E/E\approx 10^{-15}$), and the radius also remains constant to machine precision.

8.2.2 Radial geodesics

The velocity of a radially in-falling particle beginning at rest in the Schwarzschild metric is given by

[TABLE]

where $r_{0}$ is the starting radius. Figure 15 shows our calculations of radially in-falling test particles from several radii, compared to the exact solution. We integrate with constant time step $\Delta t=0.01M$ until the particles reach $r=2.001M$. The largest $L_{2}$ error out of the 10 simulations occurs for the $r_{0}=4M$ case, with $L_{2}=6.7\times 10^{-8}$ when using data points every $31$ time steps. This confirms that radial infall is accurately captured by our code.

8.2.3 Apsidal advance

Figure 16 (left) shows the trajectory of a precessing orbit. The orbit begins at $r_{0}=90M$ with tangential velocity $v_{\mathrm{t}}=0.0521157$ such that the pericentre distance is $\approx 10M$, i.e.

[TABLE]

The trajectory of our simulation follows the trajectory as computed from a direct integration of the geodesic equation with coordinate time derivatives (see e.g. Equation 4.2 of Tejeda et al. 2017). We find a change in azimuthal angle of approximately $82.3^{\circ}$ between the first apocentre and the starting position. This is consistent with the exact value of $82.4^{\circ}$ (see e.g. Equation 2 of Wegg 2012) to within measurement error. A more accurate value can be obtained with higher temporal resolution at the apocentre. This confirms that the apsidal advance is as expected from theory.

8.3 Orbital dynamics — Kerr Metric

We perform further tests of orbital dynamics using the Kerr metric in spherical Boyer-Lindquist coordinates, however the same results can be reproduced using the equivalent Cartesian representation (see Appendix C.2). We again use a black hole mass of $M=1$ in all simulations.

8.3.1 Circular orbits

Circular geodesics are allowed in the Kerr metric when the plane of orbit is perpendicular to the spin axis of the black hole (i.e. $\theta=\pi/2$ in Boyer-Lindquist coordinates). The orbital frequency for a circular orbit in the Kerr metric is given by (Abramowicz et al., 1978)

[TABLE]

for orbits going in the positive $\phi$ direction, where $|a|\leq M$ is the angular momentum parameter of the black hole. Prograde orbits correspond to $a>0$, whilst retrograde orbits are when $a<0$. Figure 13 (right) shows a circular orbit at a radius of $2M$ for a maximally rotating black hole ($a/M=1$), evolved for 15 orbital periods with constant $\Delta t=0.01M$ and compared to the exact orbit. Again, energy and angular momentum $(e,p_{\phi})$ are conserved to machine precision, and the radius also remains constant to machine precision.

8.3.2 Epicyclic motion

An otherwise circular orbit given a small, linear, radial perturbation, will undergo epicyclic motion in the Kerr metric. The epicyclic frequency $\kappa$ is given by (Kato, 1990; Lubow et al., 2002)

[TABLE]

We compute $\kappa$ for a range of values of $a$, by setting a series of test particles on epicyclic orbits at different radii. To induce the epicyclic motion we set the initial angular velocity to be $v_{\phi}=1.00001\Omega$ (i.e. a 0.001% increase from the circular velocity). For simplicity and because test particles do not interact, simulations with the same $a$ were computed simultaneously. Each simulation was run for $t=30000M$ to ensure that epicycles with the lowest frequency completed at least $20$ cycles. We then measure the frequency, as described in Appendix D. Figure 14 (left) shows the results of all simulations for $a\in[-1,1]$ and $M=1$. The radii at which $\kappa=0$ correspond to the innermost stable circular orbit (ISCO). The typical errors are of order $\Delta\kappa/\kappa\approx 10^{-3}$. However, the errors can be made arbitrarily small by taking a longer time series, since the smallest frequency measurement possible with a Fourier transform is $\delta\Omega=2\pi/t_{\mathrm{max}}$.

8.3.3 Vertical-oscillation frequency

A circular orbit given a small, linear, vertical perturbation, will undergo vertical oscillations in the Kerr metric. The vertical oscillation frequency $\Omega_{z}$ is given by (Kato, 1990; Lubow et al., 2002)

[TABLE]

Figure 14 (right) shows our calculations of the vertical oscillation frequency, following a similar setup as above (see Section 8.3.2). To induce vertical oscillations we perturb test particles up out of the $\theta=\pi/2$ plane by 0.001% (i.e. $\theta=\left(1-0.00001\right)\pi/2$). The exact solutions for each value of $a$ are truncated at the radius corresponding to the ISCO. Typical errors are similar to those for the epicyclic frequency ($\Delta\Omega_{z}/\Omega_{z}\approx 10^{-3}$), and can again be made smaller by taking a longer time series.

8.3.4 Apsidal advance

Figure 16 (middle & right) demonstrates the effect of black hole spin on apsidal advance. We use similar initial conditions as in Section 8.2.3 however with $a/M=\{-0.1,0,0.1\}$ in the Kerr metric and performed these calculations in Cartesian-like coordinates. We chose relatively small values of spin in order to distinguish the direction of orbital advance, however there is no limit on the values of spin that the code can handle. The initial conditions for this orbit are such that in the x–y plane the initial starting position and velocity is the same for all three orbits since $r^{2}\neq x^{2}+y^{2}+z^{2}$ in Boyer-Lindquist coordinates, i.e.

[TABLE]

where again we use $r_{0}=90M$ and $v_{\mathrm{tan}}=0.0521157$. The angle of precession decreases for prograde orbits, and increases for retrograde orbits, as expected. Combined with the high accuracy of our epicyclic frequency measurements, this confirms that we capture Kerr orbital dynamics accurately.

8.3.5 Conservation of angular momentum

Figure 17 shows the conservation of angular momentum in our timestepping algorithm described in Section 6.3. We compute the precessing orbit in the Schwarzschild metric in Cartesian coordinates (Section 8.2.3) using both double (64 bit) and quad (128 bit) precision, whilst varying the number of steps per orbit (step size), and the tolerance ($\epsilon_{x}$, $\epsilon_{p}$) of the implicit solve. The conserved angular momentum in the Schwarzschild metric is

[TABLE]

We define the orbital period for this trajectory to be the time between consecutive apocentres ($t_{\mathrm{orbit}}\approx 2390M$ with our setup). With $10^{4}$ steps per orbit, we find that we conserve angular momentum to round off error provided $\epsilon_{p}$ and $\epsilon_{x}$ are sufficiently small. Using tolerances of $10^{-15}$ (double precision) and $10^{-30}$ (quad precision) angular momentum is conserved to round-off error. Therefore in our 3D code, which is double precision, we set $\epsilon_{p}=\epsilon_{x}=10^{-15}$ by default.

8.4 General relativity - Bondi accretion

To test our ability to perform hydrodynamics in a curved spacetime and not just special relativity, we consider spherically symmetric flow in the Schwarzschild metric. Hawley et al. (1984b) give the analytic solution, which is a generalisation of non-relativistic Parker-Bondi flow. To simulate the accretion solution, we inject shells (geodesic spheres) of SPH particles from a particular radius at discrete, regular intervals. The particles are given density, velocity and thermal energy corresponding to the exact solution at the radius which they are injected. For details of the algorithm we refer the reader to the publicly available subroutine in Phantom used in applications of stellar winds (Toupin et al. 2015a, b; Siess et al. in prep).

8.4.1 Geodesic flow

The solution for marginally bound, spherically symmetric, radial flow of ‘pressure-less’ gas is given by (Hawley et al., 1984b)111Note that equation (67b) of Hawley et al. (1984b) is written incorrectly

[TABLE]

where $\rho_{0}$ and $u_{0}$ are normalisation constants. The solution assumes that $u\ll 1$ everywhere, hence we must use $u_{0}\ll 1$ to be consistent. The solution is also the same for inflow and outflow, aside from the sign of $v^{r}$. For our simulations we use normalisation constants $u_{0}=1\times 10^{-9}$ and $\rho_{0}=1$.

Following Hawley et al. (1984a), we inject fluid radially in toward the black hole from a radius of $r_{\mathrm{in}}=18.1M$. To have a smooth boundary condition at $r_{\mathrm{in}}$, we keep 10 shells of particles, uniformly spaced in time, outside of the injection radius at all times. The hydrodynamic quantities $v^{r}$, $\rho$ and $u$ are fixed to the exact solution while the shells remain outside the injection radius. As each shell moves inward and crosses $r_{\mathrm{in}}$, it is ‘released’ and left to evolve on its own. Once a shell is released, a new boundary shell is added to the outside of the remaining boundary shells in order to maintain a constant number of boundary shells. For the inner boundary, we simply remove particles from the simulation when they cross $r=2.1M$.

Figures 19 & 19 show our results for geodesic Bondi accretion with $8412$ particles per shell, a time spacing of $\delta t=0.602M$ between consecutive shells, and individual particle masses of $m_{a}=9\times 10^{-4}M$. We find that we achieve best results with artificial viscosity turned off i.e. $\alpha_{\mathrm{AV}}=0$ (Fig. 19). After evolving the simulation to $t=90M$ the $L_{2}$ errors are $5.4\times 10^{-2}$, $1.8\times 10^{-4}$ and $4.4\times 10^{-2}$ for $\rho^{*}$, $v^{r}$ and $u$, respectively.

Figure 19 shows the same simulation but with the standard dissipation parameter $\alpha_{\mathrm{AV}}=1$ for comparison. Although there are no shocks present, viscous heating leads to excess thermal energy (third panel of Figure 19). The excess heating is small but easily overwhelms the negligible thermal energy assumed in the test problem. To reduce viscous heating away from shocks, the standard approach in non-relativistic SPH is to implement shock detection switches. The most widely used is the switch proposed by (Cullen & Dehnen, 2010), however this is non-trivial to generalise to relativity. We experimented with a modified form of the older Morris & Monaghan (1997) switch but found it ineffective for the Bondi problem since it relies on $(-\partial v^{i}/\partial x^{i})$ as the shock detector, which does not switch off for adiabatic contraction.

There is a small turnover in density at the inner boundary in both cases. This is merely a consequence of our inner boundary condition and does not reflect the accuracy of our density calculation in a strong field regime. The drop off occurs because particles close to the accretion radius have fewer neighbours and thus a lower density. This effect is not unique to GRSPH and would occur wherever the inner boundary was chosen to be.

8.4.2 Sonic point flow

Dropping the assumption $u\ll 1$ means that the pressure between fluid elements is no longer negligible. The solution for this flow (Michel, 1972; Hawley et al., 1984b) can be written in terms of a temperature variable $T\equiv P/\rho=({\gamma_{\mathrm{ad}}}-1)u$, not to be confused with the gas temperature as defined in equation 2.22,

[TABLE]

where $n\equiv 1/({\gamma_{\mathrm{ad}}}-1)$ is the polytropic index and $K_{0}$ is an entropy normalisation constant. $T(r)$ is determined by the implicit equation

[TABLE]

The constants $C_{1}$ and $C_{2}$ can be both uniquely determined through the choice of a critical point $r_{c}$,

[TABLE]

where

[TABLE]

To write the solution in terms of $v^{r}(r)$ and $\rho^{*}(r)$, we use equations 2.2 & 2.4 along with

[TABLE]

as determined by equation 2.3 for the Schwarzschild metric. For ${\gamma_{\mathrm{ad}}}=5/3$, equation 8.27 yields two real solutions, one for inflow and the other for outflow (Michel, 1972). For our simulation we take care to choose the correct branch, corresponding to an inflow. We choose a critical value $r_{c}=8M$ and normalise the entropy to $K_{0}=1$.

Figures 21 & 21 show the results of our simulation. We follow the same procedure as for the geodesic flow simulation, injecting shells of particles from $r_{\mathrm{in}}=18.1M$, however this time we fill the domain all the way to the accretion radius at $r=2.1M$ before begining the simulation. This reduces the time taken to reach a steady state. Each shell consists of 8412 particles, with time intervals of $\delta t=1.404M$ between each consecutive shell, and individual particle masses of $m_{a}=4\times 10^{-4}M$.

Figure 21 shows our results with a reduced artificial viscosity parameter of $\alpha_{\mathrm{AV}}=0.1$. After evolving the simulation to $t=360M$ in order to reach a steady state, the $L_{2}$ errors are $4.2\times 10^{-2}$, $2.5\times 10^{-2}$ and $3.2\times 10^{-2}$ for $\rho^{*}$, $v^{r}$ and $u$, respectively.

In Figure 21 we show the same simulation but with $\alpha_{\mathrm{AV}}=1$ for comparison. As with the geodesic simulations above, we find excessive viscous heating when the standard artificial viscosity parameter is used. A small drop off in density can again be seen at the inner boundary, for the same reasons as discussed in the geodesic simulation. Figure 22 shows a rendered cross section of the conserved density $\rho^{*}$ and the specific thermal energy $u$ close to the inner boundary for the same simulation as Fig. 21. This demonstrates that the boundary condition is not a serious concern, so long as the event horizon is well resolved.

9 Conclusions

We have developed and benchmarked a 3D GRSPH code for hydrodynamics in strong field gravity. The novel aspects are:

Shock capturing: We demonstrated how to split viscous and thermal conduction components of the shock dissipation terms in relativistic SPH. This is important for accurate treatment of contact discontinuities (Price, 2008). 2. 2.

Entropy: We show how to formulate the equations with entropy as the conserved variable. This enables us to robustly simulate strong blast waves with no smoothing of the initial conditions without the code failing due to negative pressures. 3. 3.

Timestepping: We showed how to generalise the Leapfrog time integration scheme to relativistic SPH. We combined a symplectic and reversible method for non-separable Hamiltonians from Leimkuhler & Reich (2005) with a modified algorithm that is efficient for SPH. This allows for symplectic and reversible integration of the orbital dynamics and reversible integration of the SPH equations in the absence of dissipation. Our algorithm generalises the RESPA scheme (Tuckerman et al., 1992) employed by Price et al. (2018). 4. 4.

Benchmarking: We benchmark both our algorithm and codes against a series of standardised tests in 1D and 3D. We showed results for mildly and ultra-relativistic shocks in 1D and 3D, circular and eccentric orbits, radial geodesics, vertical and epicyclic oscillation about circular orbits in Schwarzschild and Kerr metrics, 3D relativistic spherical blast, and generalised Bondi flow. 5. 5.

Public code: We implemented our algorithms into the public SPH code Phantom (Price et al., 2018). Our intention is to make GR Phantom publicly available in due course.

We showed how to use SPH, for flows in curved metrics of space. We captured the shocks, and adjusted the clocks. Our tests, they all worked in each case.

Acknowledgments

We thank Paul Lasky, Yuri Levin, Lionel Siess, Kimitake Hayasaki and Clement Bonnerot for useful discussions and feedback. DP thanks Walter Dehnen and Stephan Rosswog for useful discussions. We also thank the anonymous referee for comments which have improved the paper. We acknowledge CPU time on OzSTAR, funded by Swinburne University and the Australian Government. DL is funded through the Australian government RTP (Research Training Program) stipend. DP is grateful for funding from the Australian Research Council via DP130102078 and FT130100034.

Appendix A Positivity of the entropy change

A.1 Conductivity

First, we will prove that the artificial conductivity contributes to a positive definite change to the total entropy. Taking Equation 5.8 and relating it back to the specific entropy

[TABLE]

The change in total entropy is

[TABLE]

and thus the contribution from artificial conductivity is

[TABLE]

where we have used the relation $U^{0}=\Gamma/\alpha$, and defined the negative, symmetric term

[TABLE]

Splitting Eq. A.4 into two equal sums, interchanging the labels on one, and then recombining into a single sum, we get

[TABLE]

We note that $U^{0}$ is a strictly positive quantity, and $u\propto T$. Hence, if $u_{a}>u_{b}$, we must have $T_{a}>T_{b}$, and

[TABLE]

Thus the contribution to the total entropy, due to artificial conductivity, is always positive in both special and general relativity.

A.2 Viscosity

Next we consider the change in entropy per particle due to artificial viscosity in Cartesian coordinates

[TABLE]

where $v^{*}=N_{i}\,v^{i}$ is a line of sight velocity, and we defined the positive quantity

[TABLE]

In the case that $V^{*}$ is a monotonically increasing function of $v^{*}$ in a given spacetime, if $v_{a}^{*}>v_{b}^{*}$ we must have $V_{a}^{*}>V_{b}^{*}$ as well as $\Gamma_{a}^{*}>\Gamma_{b}^{*}$, and hence

[TABLE]

The condition that $V^{*}$ be monotonic is easily satisfied in the Minkowski metric, where $V^{i}=v^{i}$, since there is no shift, $\beta^{i}=0$, and the lapse $\alpha=1$. Thus in the case of special relativity, the artificial viscosity will only ever increase the entropy of a particle, however this cannot be guaranteed in all spacetimes.

Appendix B Recovery of primitive variables

Step 3 of section 7 involves solving the equation

[TABLE]

for $w$, where $\rho$ and $P$ are functions of $w$ (Eq. 7.2 and 7.3). For this we use the standard Newton-Raphson algorithm

[TABLE]

which requires the derivative $f^{\prime}(w)$. Here, the subscript $n$ corresponds to the iteration number. For the initial prediction, $w_{0}$, we use the enthalpy computed at the previous timestep. Iterations of Equation B.2 are repeated until the following convergence criteria is satisfied

[TABLE]

for some sufficiently small number $\epsilon_{w}$. In this paper we use $\epsilon_{w}=10^{-12}$. Given the ideal gas equation of state (Eq. 2.21), we can write

[TABLE]

and hence the derivative can be written as

[TABLE]

where $\rho$, $P$, and $\Gamma$ are all functions of $w$ (Eq. 7.2, 7.3 and 7.5). Equivalently, when evolving the entropy variable $K$

[TABLE]

where instead of Equation 7.3 for the pressure, we use

[TABLE]

Appendix C Metrics in Cartesian Coordinates

C.1 Schwarzschild metric

The Schwarzschild metric line element, in standard spherical coordinates $(t,r,\theta,\phi)$ is

[TABLE]

If we make the coordinate transformation

[TABLE]

we can write the components of the covariant metric in Cartesian coordinates $(t,x,y,z)$

[TABLE]

as well as the components of the contravariant metric

[TABLE]

where $r=\sqrt{x^{2}+y^{2}+z^{2}}$. In these coordinates $\sqrt{-g}=1$. The derivatives $\partial g_{\mu\nu}/\partial x^{i}$ can be found analytically by differentiating $g_{\mu\nu}$, or through a finite difference approximation.

C.2 Kerr metric

The Kerr metric line element in standard Boyer-Lindquist coordinates $(t,r,\theta,\phi)$ is

[TABLE]

where

[TABLE]

and $|a|\leq M$ is the spin parameter. If we make the coordinate transformation

[TABLE]

we can write the components of the covariant metric in Cartesian-like coordinates $(t,x,y,z)$

[TABLE]

and the components of the contravariant metric

[TABLE]

where

[TABLE]

We again have $\sqrt{-g}=1$ in the new coordinates and the derivatives $\partial g_{\mu\nu}/\partial x^{i}$ can be found analytically by differentiating $g_{\mu\nu}$, or through a finite difference approximation.

C.3 Derivatives of the metric

We compute the derivatives of the Schwarzschild metric by analytically differentiating $g_{\mu\nu}$, however for the Kerr metric we choose to compute the derivatives through a 2nd order centred finite difference approximation

[TABLE]

for $\epsilon\ll 1$. In our simulations we use $\epsilon=10^{-8}$.

C.4 Inverse

The contravariant metric $g^{\mu\nu}$ can also be found by inverting the matrix $g_{\mu\nu}$ at every point, since

[TABLE]

Given the determinant of the covariant metric $g$, we simply compute the adjugate matrix of $g_{\mu\nu}$, and then

[TABLE]

In practice we find that this conserves Eq.C.19 better than the analytic contravariant terms near the event horizon. As such, we compute the contravariant metric in this manner by default.

Appendix D Computing oscillation frequencies

To measure the oscillation frequencies described in Sections 8.3.2 and 8.3.3, we take the maximum value of the 1-D power spectrum of $r(t)$, using numpy.fft.fft in Python. Figure 23 shows the time series (left) and corresponding power spectrum (right) for the $a=1$, $r=1.2$ case.