General Relativistic Hydrodynamics in Discrete Spacetime: Perfect Fluid Accretion onto Static and Spinning Black Holes

1. Introduction

The accretion of an idealized fluid onto a compact object (e.g. a neutron star or a black hole) remains one of the most widely-studied problems in astrophysics and cosmology, as it can be used as a minimal mathematical or numerical model for such a wide variety of phenomena, including the formation and growth of supermassive black holes at the centers of galaxies [1,2], the formation and growth of primordial black holes during the early universe [3,4], and the dynamics of pulsars, active galactic nuclei and other high-energy astrophysical phenomena [5,6]. Indeed, via modern observational techniques such as reverberation mapping [7] and measurements of quasi-periodic oscillations [8,9], the dynamics of the accretion region very close to a compact object (which are often reflected in the fast variability in the spectrum of the region’s X-ray emissions) can be used to provide a powerful and high-precision testbed for general relativity itself, for instance allowing one to test mathematical proposals such as the no-hair theorem experimentally, by determining the degree to which the exterior geometry surrounding a spinning black hole (or other compact object) appears to be well-described by the Kerr metric. This, in turn, presents the exciting possibility that deviations from the predictions of classical general relativity, for instance due to modifications in the microscopic structure of spacetime at or below the Planck scale arising from certain quantum gravity models, may become experimentally verifiable (or falsifiable) via astrophysical observations of such high-energy accretion phenomena in the near future. It would therefore be both useful and instructive to determine a robust and generic set of predictions regarding the effects of spacetime discreteness upon certain relevant accretion parameters, including the accretion rates of mass/energy and momentum onto the compact object; the lift and drag forces exerted upon the compact object (assuming that it possess a non-zero angular momentum value, and/or that the accretion is non-radial) due to the accretion flow; and the morphology, characteristics and dynamics of the accretion flow itself. The purpose of the present article is to commence the lengthy process of deriving such a set of predictions.

One of the prototypical idealized accretion cases conventionally studied is that of a compact object moving at a constant velocity through an ideal gas of uniform density (or, equivalently, an ideal gas with a uniform density and flow velocity accreting onto the compact object), commonly known as Bondi-Hoyle-Lyttleton accretion [10,11,12]. The standard interpolation formula for the rate of mass accretion in the Bondi-Hoyle-Lyttleton model is, in turn, derived from two important limiting cases: the case where the flow velocity is zero (known as Bondi accretion [11]), and the case where the flow velocity is supersonic (known as Hoyle-Lyttleton accretion [12]). In this article, we shall focus solely upon the case of (radial) Bondi accretion, and we leave the extension of these techniques to the supersonic Hoyle-Lyttleton case, and to the generalized Bondi-Hoyle-Lyttleton case, as an open research problem, ripe for future investigation. Bondi’s original analysis [11] considered the case of a spherically-symmetric distribution of ideal gas of initially uniform density (assumed to be of infinite extent), accreting onto a central point mass in pure Newtonian gravity. Michel [13] later extended Bondi’s analytic solution for spherical accretion to the general relativistic case of a static, uncharged, non-rotating black hole (as described by the Schwarzschild metric) as the central compact object, although the polytropic form of the ideal gas equation of state used within Michel’s analysis was previously shown by Taub [14] to be physically reasonable only in the strictly non-relativistic and strictly ultra-relativistic limits (with dimensionless gas temperature much less than unity, and much greater than unity, respectively), and not in the intermediate relativistic case (with dimensionless gas temperature approximately equal to unity). Rather surprisingly, Petrich, Shapiro and Teukolsky [15] were even able to extend this analysis beyond the spherically-symmetric spacetimes considered thus far to the axially-symmetric spacetime case, and hence to derive an analytic solution for the accretion of a stiff, ultra-relativistic fluid onto an uncharged but spinning black hole (as described by the Kerr metric) as the central compact object. Although the stiff, ultra-relativistic equation of state used therein is not expected to be physical (since it requires a perfect relativistic fluid whose local sound speed is equal to the speed of light), the exact solution of Petrich, Shapiro and Teukolsky nevertheless provides a useful benchmark for the testing of general relativistic hydrodynamics codes. Font and Ibáñez [16,17], and later Font, Ibáñez and Papadopoulos [18] later performed a detailed and systematic analysis of the (non-radial) case of Bondi-Hoyle-Lyttleton accretion of perfect relativistic fluids, obeying more general forms of the ideal gas equation of state, onto both static and spinning black holes by means of numerical simulations using the so-called $3+1$ “Valencia” formulation of the equations of general relativistic hydrodynamics in hyperbolic conservation law form, due originally to Banyuls, Font, Ibáñez, Martí and Miralles [19].

Some of the key insights yielded by this combination of analytical and numerical work included: the discovery of a systematic reduction in the mass accretion rate as a function of the black hole spin (an effect which becomes more significant for higher gas temperatures and higher values of the adiabatic exponent) [18]; the discovery of a drag force exerted on the black hole, either due to the presence of a downstream region of high fluid density caused by non-radial fluid motion in the case of Bondi-Hoyle-Lyttleton accretion, or redistribution of high fluid pressure regions caused by non-zero black hole spin, or both (which, in turn, result in increased gravitational forces exerted on the black hole by the fluid), with an absence of the “flip-flop” fluid instabilities typically seen in purely Newtonian accretion simulations [16,17]; and the discovery of a lift force, analogous to the Magnus effect in Newtonian fluid dynamics, exerted on spinning black holes by non-radially accreting fluids due to the asymmetry of the fluid pressure redistribution (with more pressure being redistributed onto the side of the black hole that is counter-rotating with the fluid) [18]. The primary objective of the present work is to begin the process of determining what kind of effect an underlying discretization of the background spacetime is expected to have on these types of black hole accretion phenomena, as well as on other related astrophysical processes. Since discreteness of the fundamental structure of spacetime is a generic feature of many proposed models of quantum gravity, including casual set theory [20,21,22,23], causal dynamical triangulations [24,25], loop quantum gravity [26,27,28], and the Wolfram model [29,30,31,32,33], it is hoped that such an investigation will eventually enable the observational investigation of certain classes of quantum gravity theories by means of astrophysical probes of near-black hole accretion regions. To this end, we make use of the Gravitas computational general relativity framework [34,35], which allows for the configuration, execution, visualization and analysis of complex numerical relativity simulations in both discrete and continuous spacetime settings, by combining a powerful tensor calculus and differential geometry framework on the analytical side, with a sophisticated hypergraph-based adaptive refinement system [36,37] on the numerical side. Most general relativistic simulations of black hole accretion consider a perfect relativistic fluid evolving on top of a fixed, time-independent spacetime metric (typically representing either a Schwarzschild geometry or a Kerr geometry), and thereby neglect any gravitational effects of the fluid density on the black hole itself. Since the mass densities of the fluids in question are usually much smaller than the mass of the central black hole, this is often not an unreasonable simplification to make (this is conventionally referred to as the “test-fluid” assumption within the relativistic hydrodynamics literature [38]). However, since many of the effects in which we are interested for the purposes of this article (such as drag forces exerted by a fluid upon a spinning black hole) depend crucially upon the two-way gravitational interaction between the black hole and the fluid, we do not make this assumption here. Instead, we use Gravitas to configure and run fully general relativistic two-way coupled simulations, evolving the fluid variables and the metric tensor together in parallel.

We begin in Section 2 with a brief overview of the purely hyperbolic $3+1$ “Valencia” formalism for general relativistic hydrodynamics of Banyuls, Font, Ibáñez, Martí and Miralles [19], together with a description of how a modified version of the formalism can be derived that is specifically adapted for numerical relativistic hydrodynamics in discrete spacetimes, by means of the discrete spacetime ADM formalism already implemented within the Gravitas framework [35]. The final result of this analysis will be the derivation of a complete and fully-coupled system of purely hyperbolic equations for the evolution of the components of the discrete spatial metric tensor and the discrete spacetime fluid variables jointly, together with a set of purely elliptic constraint equations for the discrete spacetime gauge, which can then be implemented directly into Gravitas. We proceed in Section 3 to present a weak, integral form of these equations that is amenable to direct numerical solution via finite-volume methods, and we validate the resulting numerical implementation against a standard special relativistic hydrodynamics shock tube problem (namely the mildly-relativistic blast wave problem of Donat, Font, Ibaéñez and Marquina [39]). Particular attention is paid to the validation of the implementation of the conservative-to-primitive variable reconstruction algorithm, which is generally a non-trivial operation in relativistic hydrodynamics and for which we follow the approach of Eulderink and Mellema [40] in deriving a one-dimensional iterative Newton-Raphson solver, which works generically for any ideal gas equation of state. Finally, in Section 4, we show the numerical results of our general relativistic hydrodynamics simulations, beginning with a simulation of radial (Bondi-type) accretion onto a static/Schwarzschild black hole, before proceeding to radial (Bondi-type) accretion onto spinning/Kerr black holes, with a variety of spin values, ranging from modest to near-extremal. The broad qualitative features (e.g. the shape of the density profile for the accretion region in the Schwarzschild case, or the splitting of the accretion region into several distinct “arms” in the rapidly-spinning Kerr case, etc.) of these simulations appear similar to those obtained from analogous general relativistic hydrodynamics simulations performed in continuous spacetime geometries, although a more rigorous quantitative analysis reveals certain notable discrepancies. In particular, we find that the rates of mass/energy and momentum accretion onto the black hole both appear to be monotonically-decreasing functions of the discretization scale of the underlying spacetime, with increased black hole spin values, higher fluid temperatures and larger values of the adiabatic exponent (i.e. stiffer equations of state) accentuating and amplifying this discretization effect. Moreover, we discover that the drag force exerted on the black hole exhibits a sensitive dependence upon the underlying discretization scale, with certain critical values of the discretization scale resulting in apparent instabilities in the spacetime structure, observable within the feedback effect of the fluid density onto the black hole geometry. All simulation results presented within this (and the previous) section are presented in both “horizon-adapted” and “non-horizon-adapted” coordinate systems, as proposed by Font, Ibáñez and Papadopoulos [41], so as to eliminate the possibility that any of these effects might simply be a byproduct of unphysical fluid behavior resulting from certain numerical divergences near the black hole horizon. We conclude in Section 5 with a brief discussion of potential astrophysical implications of these results, as well as directions for future research and investigation.

Note that all of the Gravitas functionality necessary to reproduce the results presented within this article can be found in the Gravitas GitHub repository, with extensive documentation available within both the Wolfram Function Repository (e.g. ADMDecomposition and StressEnergyTensor) and within the two previous articles [34,35]. This article follows all of the same notational and terminological conventions as these two previous articles; in particular, we assume geometric units with $c=G=\hslash =1$, we employ a metric signature of $\left(-,+,+,+\right)$ in all relevant cases, and the Einstein summation convention is assumed throughout (such that all repeated tensor indices are implicitly summed over).

2. General Relativistic Hydrodynamics in Discrete Spacetime

In order to derive a form of the equations of general relativistic hydrodynamics that is suitable for analysis within a discrete spacetime setting, we begin by considering the $3+1$ “Valencia” formulation of the curved spacetime hydrodynamics equations in conservation law form due to Banyuls, Font, Ibáñez, Martí and Miralles [19], which exploits the fundamentally hyperbolic character of the spacetime continuity equations. The equations of general relativistic hydrodynamics represent a mathematical encoding of two distinct physical laws, namely the law of conservation of energy-momentum, and the law of conservation of baryon number. Assuming a spacetime given by a smooth n-dimensional Lorentzian manifold $\left(\mathcal{M},g\right)$, the law of conservation of energy-momentum can be represented as a statement that the covariant divergence of the rank-2 stress-energy tensor $T^{\mu \nu }$ vanishes identically:

$${\nabla }_{\nu }{T}^{\mu \nu }=\frac{\partial }{\partial x^{\nu }}\left(T^{\mu \nu }\right)+{\Gamma }_{\nu \sigma }^{\mu }{T}^{\sigma \nu }+{\Gamma }_{\nu \sigma }^{\nu }{T}^{\mu \sigma }=0,$$

(1)

while the law of conservation of baryon number can be represented as a statement that the covariant divergence of the rank-1 (rest) mass current vector $J^{\mu }$ also vanishes identically:

$${\nabla }_{\mu }{J}^{\mu }=\frac{\partial }{\partial x^{\mu }}\left(J^{\mu }\right)+{\Gamma }_{\mu \sigma }^{\mu }{J}^{\sigma }=0.$$

(2)

In the above, the spacetime covariant derivative ${\nabla }_{\mu }$ is represented in terms of the coefficients of the Levi-Civita connection ∇ on the manifold $\left(\mathcal{M},g\right)$, namely the spacetime Christoffel symbols ${\Gamma }_{\mu \nu }^{\rho }$, themselves represented in terms of partial derivatives of the spacetime metric tensor $g_{\mu \nu }$:

$${\Gamma }_{\mu \nu }^{\rho }=\frac{1}{2}{g}^{\rho \sigma }\left(\frac{\partial }{\partial x^{\mu }}\left(g_{\sigma \nu }\right)+\frac{\partial }{\partial x^{\nu }}\left(g_{\mu \sigma }\right)-\frac{\partial }{\partial x^{\sigma }}\left(g_{\mu \nu }\right)\right).$$

(3)

For the specific case of a perfect relativistic fluid in equilibrium, obtained by neglecting all considerations of heat conduction, fluid viscosity and shear stress, the stress-energy tensor $T^{\mu \nu }$ and (rest) mass current vector $J^{\mu }$ take the forms:

$$T^{\mu \nu }=\rho h{u}^{\mu }{u}^{\nu }+P{g}^{\mu \nu },\mathrm{and}{J}^{\mu }=\rho u^{\mu },$$

(4)

respectively, where $\rho$ denotes the (rest) mass density of the fluid, P denotes its hydrostatic pressure, $u^{\mu }$ denotes its spacetime velocity, and h denotes its specific relativistic enthalpy:

$$h=1+\epsilon \left(\rho ,P\right)+\frac{P}{\rho },$$

(5)

where $\epsilon \left(\rho ,P\right)$ represents the specific internal energy of the fluid. The product $\rho h$ of the (rest) mass density $\rho$ and the specific relativistic enthalpy h constitutes the total mass-energy density of the fluid. In all of the above, the components $g^{\mu \nu }$ are components of the inverse metric tensor $g^{\mu \nu }={\left(g_{\mu \nu }\right)}^{-1}$, and the tensor indices $\mu ,\nu ,\rho ,\sigma$ range across all spacetime coordinate directions $\left\{0,\dots ,n-1\right\}$ (with $\left\{x^{\mu }\right\}$ being a local spacetime coordinate basis), in contrast to the $3+1$ decomposition formalism discussed below. The resulting system of equations may then be closed by defining an appropriate equation of state, allowing one either to calculate the specific internal energy as a function of the fluid density and hydrostatic pressure $\epsilon \left(\rho ,P\right)$, or, equivalently, to calculate the hydrostatic pressure as a function of the fluid density and specific internal energy $P\left(\rho ,\epsilon \right)$. The equation of state thus allows one to compute the local sound speed $c_s$ of the fluid as:

$$c_s=\frac{1}{\sqrt{h}}\sqrt{{\left.\left(\frac{\partial P}{\partial \rho }\right)\right|}_{\epsilon }+\left(\frac{P}{{\rho }^2}\right){\left.\left(\frac{\partial P}{\partial \epsilon }\right)\right|}_P},$$

(6)

where ${\left.\left(\frac{\partial P}{\partial \rho }\right)\right|}_{\epsilon }$ and ${\left.\left(\frac{\partial P}{\partial \epsilon }\right)\right|}_{\rho }$ denote partial derivatives assuming fixed internal energy $\epsilon$ and fixed fluid density $\rho$, respectively.

We now proceed to perform a “$3+1$ decomposition” (or “foliation”) of our n-dimensional spacetime $\left(\mathcal{M},g\right)$ into a time-ordered sequence of $\left(n-1\right)$-dimensional spacelike hypersurfaces of Riemannian signature, each with an induced/spatial metric tensor ${\gamma }_{\mu \nu }$, by means of the ADM formalism due originally to Arnowitt, Deser and Misner [42,43], and later adapted by York [44] into the form used for the purposes of this article. Within such a decomposition, the overall spacetime line element (or first fundamental form) $d{s}^2$, which normally takes the general form:

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu },$$

(7)

with $\mu ,\nu$ ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$ (and with $\left\{x^{\mu }\right\}$ taken to represent a local spacetime coordinate basis), can now be written instead as:

$$\begin{array}{c}d{s}^2=-{\alpha }^{2}d{t}^2+{\gamma }_{\mu \nu }\left(d{x}^{\mu }+{\beta }^{\mu }dt\right)\left(d{x}^{\nu }+{\beta }^{\nu }dt\right)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\left(-{\alpha }^2+{\gamma }_{\mu \sigma }{\beta }^{\sigma }{\beta }^{\mu }\right)d{t}^2+2{\gamma }_{\mu \sigma }{\beta }^{\sigma }dtd{x}^{\mu }+{\gamma }_{\mu \nu }d{x}^{\mu }d{x}^{\nu }\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\left(-{\alpha }^2+{\beta }_{\mu }{\beta }^{\mu }\right)d{t}^2+2{\beta }_{\mu }dtd{x}^{\mu }+{\gamma }_{\mu \nu }d{x}^{\mu }d{x}^{\nu },\end{array}$$

(8)

with $\mu ,\nu ,\sigma$ ranging across the spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only (and with $\left\{x^{\mu }\right\}$ now taken to represent a local spatial coordinate basis on each hypersurface), and where t designates a distinguished “time” coordinate. In the above, the scalar field $\alpha$ (known as the lapse function) and the $\left(n-1\right)$-dimensional vector field ${\beta }^{\mu }$ (known as the shift vector) correspond to the Lagrange multipliers of the ADM formalism, representing the proper time distance $d\tau$ between corresponding points on the neighboring spacelike hypersurfaces labeled by coordinate time values $t=t_0$ and $t=t_0+dt$:

$$d\tau \left(t_0,t_0+dt\right)=\alpha dt,$$

(9)

as measured in the direction $\mathbf{n}$ normal to the $t=t_0$ hypersurface, and the relabeling of the spatial coordinate basis $x^{\mu }\left(t_0\right)$ as one moves from the $t=t_0$ hypersurface to the neighboring $t=t_0+dt$ hypersurface:

$$x^{\mu }\left(t_0+dt\right)=x^{\mu }\left(t_0\right)-{\beta }^{\mu }dt,$$

(10)

respectively. The unit vector $\mathbf{n}$ that is normal to each spacelike hypersurface is given by the spacetime contravariant derivative ${ }^{\left(4\right)}{\nabla }^{\mu }$ of the distinguished time coordinate t:

$$n^{\mu }=-\alpha { }^{\left(4\right)}{\nabla }^{\mu }t=-\alpha g^{\mu \sigma }{ }^{\left(4\right)}{\nabla }_{\sigma }t=-\alpha g^{\mu \sigma }\frac{\partial }{\partial x^{\sigma }}\left(t\right),$$

(11)

while the “time vector” $\mathbf{t}$ that determines how points on the $t=t_0$ hypersurface map to corresponding points on the $t=t_0+dt$ hypersurface is given by:

$$t^{\mu }=\alpha n^{\left(\mu +1\right)}+{\beta }^{\mu }=-{\alpha }^2{g}^{\mu \sigma }\frac{\partial }{\partial x^{\sigma }}\left(t\right)+{\beta }^{\mu },$$

(12)

with $\mu ,\sigma$ ranging across the spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and where we have introduced the notational convention of using a bracketed “4” to designate spacetime quantities (such that ${ }^{\left(4\right)}{\nabla }_{\mu }$ denotes the spacetime covariant derivative, as defined above in terms of the spacetime Christoffel symbols ${\Gamma }_{\mu \nu }^{\rho }$, which are henceforth denoted ${ }^{\left(4\right)}{\Gamma }_{\mu \nu }^{\rho }$), in order to distinguish them from the corresponding spatial quantities, for which we use a bracketed “3” instead. Interpreting the ADM formalism as a Hamiltonian formulation of the Einstein field equations, we see that the components ${\gamma }_{\mu \nu }$ of the spatial metric tensor represent the dynamical variables of the theory, with the components $K_{\mu \nu }$ of the extrinsic curvature tensor (or second fundamental form) representing the corresponding conjugate momenta. These components can be obtained by computing the Lie derivative $\mathcal{L}$ of the spatial metric tensor ${\gamma }_{\mu \nu }$ in the direction of the normal vector $\mathbf{n}$ [45]:

$$K_{\mu \nu }=-\frac{1}{2}{\mathcal{L}}_{\mathbf{n}}{\gamma }_{\mu \nu },$$

(13)

which expands out to give, explicitly:

$$\begin{array}{c}{K}_{\mu \nu }=\frac{1}{2\alpha }\left({ }^{\left(3\right)}{\nabla }_{\nu }{\beta }_{\mu }+{ }^{\left(3\right)}{\nabla }_{\mu }{\beta }_{\nu }-\frac{\partial }{\partial t}\left({\gamma }_{\mu \nu }\right)\right)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{1}{2\alpha }\left(\frac{\partial }{\partial x^{\nu }}\left({\beta }_{\mu }\right)-{ }^{\left(3\right)}{\Gamma }_{\nu \mu }^{\sigma }{\beta }_{\sigma }+\frac{\partial }{\partial x^{\mu }}\left({\beta }_{\nu }\right)-{ }^{\left(3\right)}{\Gamma }_{\mu \nu }^{\sigma }{\beta }_{\sigma }-\frac{\partial }{\partial t}\left({\gamma }_{\mu \nu }\right)\right),\end{array}$$

(14)

where the spatial covariant derivative ${ }^{\left(3\right)}{\nabla }_{\mu }$ is represented in terms of the coefficients of the induced Levi-Civita connection ${ }^{\left(3\right)}\nabla$ on each spacelike hypersurface, namely the spatial Christoffel symbols ${ }^{\left(3\right)}{\Gamma }_{\mu \nu }^{\rho }$, themselves represented in terms of partial derivatives of the spatial metric tensor ${\gamma }_{\mu \nu }$:

$${ }^{\left(3\right)}{\Gamma }_{\mu \nu }^{\rho }=\frac{1}{2}{\gamma }^{\rho \sigma }\left(\frac{\partial }{\partial x^{\mu }}\left({\gamma }_{\sigma \nu }\right)+\frac{\partial }{\partial x^{\nu }}\left({\gamma }_{\mu \sigma }\right)-\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)\right).$$

(15)

Note also that the indices of the shift vector $\beta$ are raised and lowered using the spatial metric tensor ${\gamma }_{\mu \nu }$, and so, in particular, the covector form ${\beta }_{\mu }$ used above is given by:

$${\beta }_{\mu }={\gamma }_{\mu \sigma }{\beta }^{\sigma }.$$

(16)

In all of the above, $\mu ,\nu ,\sigma$ range across the spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only.

Just as one can perform a $3+1$/ADM decomposition of the overall spacetime metric $g_{\mu \nu }$ that appears on the left-hand-side of the Einstein field equations, one can equivalently perform a $3+1$/ADM decomposition of the overall spacetime stress-energy tensor $T^{\mu \nu }$ that appears on the right-hand-side of the Einstein field equations [46]. By projecting the continuity equations for the stress-energy tensor $T^{\mu \nu }$:

$${ }^{\left(4\right)}{\nabla }_{\nu }{T}^{\mu \nu }=\frac{\partial }{\partial x^{\nu }}\left(T^{\mu \nu }\right)+{ }^{\left(4\right)}{\Gamma }_{\nu \sigma }^{\mu }{T}^{\sigma \nu }+{ }^{\left(4\right)}{\Gamma }_{\nu \sigma }^{\nu }{T}^{\mu \sigma }=0,$$

(17)

with $\mu ,\nu ,\sigma$ ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$, in the purely timelike direction, we obtain the energy conservation equation:

$$\frac{\partial }{\partial t}\left(E\right)-{\mathcal{L}}_{\beta }E+\alpha \left({ }^{\left(3\right)}{\nabla }_{\mu }{p}^{\mu }-KE-K_{\mu \nu }{S}^{\mu \nu }\right)+2p^{\mu }{ }^{\left(3\right)}{\nabla }_{\mu }\alpha =0,$$

(18)

where the Lie derivative term ${\mathcal{L}}_{\beta }E$ expands to give:

$$\begin{array}{c}\frac{\partial }{\partial t}\left(E\right)-{\beta }^{\mu }\frac{\partial }{\partial x^{\mu }}\left(E\right)+\alpha \left({ }^{\left(3\right)}{\nabla }_{\mu }{p}^{\mu }-KE-K_{\mu \nu }{S}^{\mu \nu }\right)+2p^{\mu }{ }^{\left(3\right)}{\nabla }_{\mu }\alpha \hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{\partial }{\partial t}\left(E\right)-{\beta }^{\mu }\frac{\partial }{\partial x^{\mu }}\left(E\right)+\alpha \left(\frac{\partial }{\partial x^{\mu }}\left(p^{\mu }\right)+{ }^{\left(3\right)}{\Gamma }_{\mu \sigma }^{\mu }{p}^{\sigma }-KE-K_{\mu \nu }{S}^{\mu \nu }\right)+2p^{\mu }\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)=0,\end{array}$$

(19)

with $\mu ,\nu ,\sigma$ ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. On the other hand, projecting in the $\left(n-1\right)$ purely spacelike directions yields the momentum conservation equations:

$$\frac{\partial }{\partial t}\left(p_{\mu }\right)-{\mathcal{L}}_{\beta }{p}_{\mu }+\alpha { }^{\left(3\right)}{\nabla }_{\nu }{S}_{\mu }^{\nu }+S_{\mu \nu }{ }^{\left(3\right)}{\nabla }^{\nu }\alpha -\alpha K{p}_{\mu }+E{ }^{\left(3\right)}{\nabla }_{\mu }\alpha =0,$$

(20)

where the Lie derivative term ${\mathcal{L}}_{\beta }{p}_{\mu }$ expands, and the contravariant derivative operator ${ }^{\left(3\right)}{\nabla }^{\nu }$ may be replaced with a corresponding covariant derivative operator ${ }^{\left(3\right)}{\nabla }_{\sigma }$, to give:

$$\begin{array}{c}\frac{\partial }{\partial t}\left(p_{\mu }\right)-{\beta }^{\sigma }\frac{\partial }{\partial x^{\sigma }}\left(p_{\mu }\right)-p_{\sigma }\frac{\partial }{\partial x^{\mu }}\left({\beta }^{\sigma }\right)+\alpha { }^{\left(3\right)}{\nabla }_{\nu }{S}_{\mu }^{\nu }+S_{\mu \nu }{\gamma }^{\nu \sigma }{ }^{\left(3\right)}{\nabla }_{\sigma }\alpha -\alpha K{p}_{\mu }+E{ }^{\left(3\right)}{\nabla }_{\mu }\alpha \hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{\partial }{\partial t}\left(p_{\mu }\right)-{\beta }^{\sigma }\frac{\partial }{\partial x^{\sigma }}\left(p_{\mu }\right)-p_{\sigma }\frac{\partial }{\partial x^{\mu }}\left({\beta }^{\sigma }\right)+\alpha \left(\frac{\partial }{\partial x^{\nu }}\left(S_{\mu }^{\nu }\right)+{ }^{\left(3\right)}{\Gamma }_{\nu \sigma }^{\nu }{S}_{\mu }^{\sigma }-{ }^{\left(3\right)}{\Gamma }_{\nu \mu }^{\sigma }{S}_{\sigma }^{\nu }\right)\\ \hfill +S_{\mu \nu }{\gamma }^{\nu \sigma }\frac{\partial }{\partial x^{\sigma }}\left(\alpha \right)-\alpha K{p}_{\mu }+E\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)=0,\end{array}$$

(21)

with $\mu ,\nu ,\sigma$ again ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. Note that, in the above, E, $p_{\mu }$ and $S_{\mu \nu }$ denote the energy density, the momentum density (in covector form) and the (Cauchy) stress tensor, respectively, of the stress-energy distribution described by $T^{\mu \nu }$, as perceived by an observer moving in the direction $\mathbf{n}$ normal to the spacelike hypersurfaces, which can be calculated via the componentwise projections:

$$E=T_{\mu \nu }{n}^{\mu }{n}^{\nu },p_{\alpha }=-T_{\mu \nu }{n}^{\mu }{\perp }_{\alpha }^{\nu },\mathrm{and}{S}_{\alpha \beta }=T_{\mu \nu }{\perp }_{\alpha }^{\mu }{\perp }_{\beta }^{\nu },$$

(22)

with $\alpha ,\beta$ ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and $\mu ,\nu$ ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$, respectively. Here, ${\perp }_{\mu }^{\nu }$ are the components of the orthogonal projector (i.e. the projection operator in the normal direction $\mathbf{n}$):

$${\perp }_{\mu }^{\nu }={\delta }_{\left(\mu +1\right)}^{\nu }+n_{\left(\mu +1\right)}{n}^{\nu },$$

(23)

where ${\delta }_{\mu }^{\nu }$ is the identity tensor/Kronecker delta function, and the momentum vector $\mathbf{p}$ and (Cauchy) stress tensor $S_{\mu \nu }$ are raised and lowered using the spatial metric tensor ${\gamma }_{\mu \nu }$, and so, in particular, for the quantities $p^{\mu }$ (in vector form), $S_{\mu }^{\nu }$ (in mixed-index form) and $S^{\mu \nu }$ (in contravariant form) appearing in the equations above, one has:

$$p^{\mu }={\gamma }^{\mu \sigma }{p}_{\sigma },S_{\mu }^{\nu }={\gamma }_{\mu \sigma }{S}^{\sigma \nu }={\gamma }_{\mu \sigma }{\gamma }^{\lambda \nu }{S}_{\lambda }^{\sigma }={\gamma }^{\sigma \nu }{S}_{\mu \sigma },S^{\mu \nu }={\gamma }^{\mu \sigma }{S}_{\sigma }^{\nu }={\gamma }^{\sigma \nu }{S}_{\sigma }^{\mu }={\gamma }^{\mu \sigma }{\gamma }^{\lambda \nu }{S}_{\sigma \lambda },$$

(24)

respectively, where $\mu ,\nu ,\sigma ,\lambda$ range across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and where ${\gamma }^{\mu \nu }$ are components of the inverse spatial metric tensor ${\gamma }^{\mu \nu }={\left({\gamma }_{\mu \nu }\right)}^{-1}$. Moreover, the indices of the stress-energy tensor $T^{\mu \nu }$ and the normal vector $\mathbf{n}$ are raised and lowered using the spacetime metric tensor $g_{\mu \nu }$, and so, in particular, for the covariant forms $T_{\mu \nu }$ and $n_{\mu }$ appearing above, one has:

$$T_{\mu \nu }=g_{\mu \sigma }{g}_{\lambda \nu }{T}^{\sigma \lambda }=g_{\sigma \nu }{T}_{\mu }^{\sigma }=g_{\mu \sigma }{T}_{\nu }^{\sigma },n_{\mu }=g_{\mu \sigma }{n}^{\sigma },$$

(25)

with $\mu ,\nu ,\sigma ,\lambda$ ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$. We have also introduced the notation K to indicate the trace of the extrinsic curvature tensor $K_{\mu \nu }$, i.e:

$$K=K_{\mu }^{\mu }={\gamma }^{\mu \nu }{K}_{\mu \nu }.$$

(26)

Upon comparing the decomposition of the stress-energy tensor $T^{\mu \nu }$ to the decomposition of the spacetime metric tensor $g_{\mu \nu }$, we see that the energy density E plays the same as the lapse function $\alpha$, the momentum density covector $p_{\mu }$ plays the same role as the shift vector ${\beta }^{\mu }$, and the (Cauchy) stress tensor $S_{\mu \nu }$ plays the same role as the induced/spatial metric tensor ${\gamma }_{\mu \nu }$.

Returning now from considerations of the general ADM formalism to the specific case of general relativistic hydrodynamics, we proceed to consider the (spatial) fluid velocity $\mathbf{v}$, as perceived by an observer moving in the direction $\mathbf{n}$ normal to the spacelike hypersurfaces, namely:

$$v^{\mu }=\frac{u^{\left(\mu +1\right)}}{\alpha u^0}+\frac{{\beta }^{\mu }}{\alpha },$$

(27)

where $\alpha u^0$ represents the Lorentz factor of the fluid:

$$\alpha u^0=-u_{\left(\mu +1\right)}{n}^{\left(\mu +1\right)}=\frac{1}{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}},$$

(28)

with $\mu ,\nu$ in all of the above ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. We shall henceforth treat the fluid (rest) mass density $\rho$, the (spatial) fluid velocity components for a normal observer $v^{\mu }$, and the fluid pressure P, as the primitive variables of our forthcoming system of hyperbolic partial differential equations in conservation law form. For a perfect relativistic fluid, the energy conservation equation obtained from taking a timelike projection of the stress-energy continuity equations becomes:

$$\begin{array}{c}\frac{1}{\sqrt{-det\left(g_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right)\right.\hfill \\ \left.+\frac{\partial }{\partial x^{\rho }}\left(\sqrt{-det\left(g_{\mu \nu }\right)}\left(\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)+P{v}^{\rho }\right)\right)\right]\\ \hfill =\alpha \left(T^{\mu 0}\frac{\partial }{\partial x^{\mu }}\left(log\left(\alpha \right)\right)-T^{\mu \nu \left(4\right)}{\Gamma }_{\nu \mu }^0\right),\end{array}$$

(29)

with $\mu ,\nu ,\rho$ on the left-hand-side of the equation ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and $\mu ,\nu$ on the right-hand-side of the equation ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$. We have introduced the notation $det\left(g_{\mu \nu }\right)$ and $det\left({\gamma }_{\mu \nu }\right)$ in the above to represent determinants of the spacetime and spatial metric tensors $g_{\mu \nu }$ and ${\gamma }_{\mu \nu }$, respectively (regarded here as explicit matrices in covariant form); within this notation, the indices $\mu$ and $\nu$ should therefore be thought of as being purely “structural”. Likewise, the momentum conservation equations obtained from taking $\left(n-1\right)$ spacelike projections of the stress-energy continuity equations become:

$$\begin{array}{c}\frac{1}{\sqrt{-det\left(g_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\right)\right.\hfill \\ \left.+\frac{\partial }{\partial x^{\rho }}\left(\sqrt{-det\left(g_{\mu \nu }\right)}\left(\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)+P{\delta }_{\sigma }^{\rho }\right)\right)\right]\\ \hfill =T^{\mu \nu }\left(\frac{\partial }{\partial x^{\mu }}\left(g_{\nu \left(\sigma +1\right)}\right)-{ }^{\left(4\right)}{\Gamma }_{\nu \mu }^{\lambda }{g}_{\lambda \left(\sigma +1\right)}\right),\end{array}$$

(30)

with $\mu ,\nu ,\rho$ on the left-hand-side of the equation again ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, $\mu ,\nu ,\lambda$ on the right-hand-side of the equation ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$, and with $\sigma$ on both sides ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. Finally, the baryon number continuity equation:

$${ }^{\left(4\right)}{\nabla }_{\mu }{J}^{\mu }=\frac{\partial }{\partial x^{\mu }}\left(J^{\mu }\right)+{ }^{\left(4\right)}{\Gamma }_{\mu \sigma }^{\mu }{J}^{\sigma }=0,$$

(31)

yields:

$$\begin{array}{c}\frac{1}{\sqrt{-det\left(g_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right)\right.\hfill \\ \hfill \left.+\frac{\partial }{\partial x^{\rho }}\left(\sqrt{-det\left(g_{\mu \nu }\right)}\left(\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)\right)\right)\right]=0,\end{array}$$

(32)

with $\mu ,\nu ,\rho$ ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. In the above, the indices of the (spatial) fluid velocity vector $\mathbf{v}$ are raised and lowered using the spatial metric tensor ${\gamma }_{\mu \nu }$, and so, in particular, one has the covector form:

$$v_{\mu }={\gamma }_{\mu \sigma }{v}^{\sigma },$$

(33)

with $\mu ,\sigma$ ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. The conserved quantity appearing within the baryon number continuity equation represents the (rest) mass density D of the fluid as measured by an observer moving in the normal direction $\mathbf{n}$:

$$D=\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}=-J_{\mu }{n}^{\mu },$$

(34)

with $\mu ,\nu$ on the left-hand-side ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and $\mu$ on the right-hand-side ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$; the conserved quantity appearing within the energy conservation equation is the difference between the energy density E measured by a normal observer and the (rest) mass density D measured by that same observer:

$$E-D=\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}=T_{\mu \nu }{n}^{\mu }{n}^{\nu }-J_{\mu }{n}^{\mu },$$

(35)

with $\mu ,\nu$ on the left-hand-side ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and $\mu ,\nu$ on the right-hand-side ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$; while, finally, the conserved quantities appearing within the momentum conservation equations are simply the components of the momentum density $p_{\mu }$ (represented in covector form) measured by a normal observer:

$$p_{\sigma }=\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}=-T_{\mu \nu }{n}^{\mu }{\perp }_{\sigma }^{\nu },$$

(36)

with $\mu ,\nu$ on the left-hand-side ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, $\mu ,\nu$ on the right-hand-side ranging across across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$, and with $\sigma$ on both sides ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. Note that the source terms appearing on the right-hand-sides of the energy and momentum conservation equations do not contain any derivatives of the primitive variables $\rho$, $v^{\mu }$ and P, and therefore the hyperbolic character of the overall system of equations is preserved.

However, observe also that the source terms for the energy and momentum conservation equations currently depend upon the overall spacetime metric tensor $g_{\mu \nu }$, its partial derivatives and its corresponding Christoffel symbols ${ }^{\left(4\right)}{\Gamma }_{\mu \nu }^{\rho }$. Moreover, the hyperbolic equations themselves involve a dependence on the spacetime metric determinant $det\left(g_{\mu \nu }\right)$. For many simulations in general relativistic hydrodynamics, this does not present a problem, since a time-independent (and often analytic) spacetime metric, such as the Schwarzschild metric for a static black hole or the Kerr metric for a spinning one, is assumed to be fixed in advance, and then a relativistic fluid is simply evolved on top of it [47]. In such cases, all of the necessary spacetime metric components, spacetime metric derivatives and spacetime Christoffel symbols (along with the spacetime metric determinant) may be precalculated, and their analytical forms (or some appropriate numerical approximations to them) can then be incorporated into the overall simulation code. Although this is an entirely reasonable idealization to use in cases where the gravitational influence of the fluid on the underlying metric may be safely neglected (i.e. the “test-fluid” assumption [38]), which is often true in the case of black hole accretion simulations, this is clearly unsatisfactory for our present purposes, since we intend to evolve the fluid parameters and the spatial metric tensor together in a fully-coupled fashion, in order to determine the effects of spacetime discretization on both the fluid morphology and the resulting spacetime geometry. Eliminating the dependence of the equations on the spacetime metric determinant $det\left(g_{\mu \nu }\right)$ is straightforward since, due to the geometry of the ADM decomposition, this determinant can be directly related to the spatial metric determinant $det\left({\gamma }_{\mu \nu }\right)$ by means of the lapse function $\alpha$:

$$\sqrt{-det\left(g_{\mu \nu }\right)}=\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}.$$

(37)

On the other hand, by means of a somewhat more involved calculation, we can rewrite the source terms for the energy and momentum conservation equations purely in terms of components of the stress-energy tensor $T^{\mu \nu }$, the primitive variables of the fluid $\rho$, $v^{\mu }$ (or equivalently $v_{\mu }$) and P, the ADM gauge variables $\alpha$ and ${\beta }^{\mu }$, the spatial metric tensor components ${\gamma }_{\mu \nu }$, and the extrinsic curvature tensor components $K_{\mu \nu }$, as follows:

$$\begin{array}{c}\alpha \left(T^{\mu 0}\frac{\partial }{\partial x^{\mu }}\left(log\left(\alpha \right)\right)-T^{\mu \nu }{ }^{\left(4\right)}{\Gamma }_{\nu \mu }^0\right)=T^{00}\left({\beta }^{\mu }{\beta }^{\nu }{K}_{\mu \nu }-{\beta }^{\mu }\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)\right)\hfill \\ \hfill +T^{0\left(\mu +1\right)}\left(-\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)+2{\beta }^{\nu }{K}_{\mu \nu }\right)+T^{\left(\mu +1\right)\left(\nu +1\right)}{K}_{\mu \nu },\end{array}$$

(38)

with $\mu ,\nu$ on the left-hand-side ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$ and $\mu ,\nu$ on the right-hand side ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and:

$$\begin{array}{c}{T}^{\mu \nu }\left(\frac{\partial }{\partial x^{\mu }}\left(g_{\nu \left(\sigma +1\right)}\right)-{ }^{\left(4\right)}{\Gamma }_{\nu \mu }^{\lambda }{g}_{\lambda \left(\sigma +1\right)}\right)=T^{00}\left(\frac{1}{2}{\beta }^{\mu }{\beta }^{\nu }\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)-\alpha \frac{\partial }{\partial x^{\sigma }}\left(\alpha \right)\right)+T^{0\left(\mu +1\right)}{\beta }^{\nu }\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)\hfill \\ \hfill +\frac{1}{2}{T}^{\left(\mu +1\right)\left(\nu +1\right)}\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)+\frac{\rho h{v}_{\rho }}{\alpha \left(1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }\right)}\frac{\partial }{\partial x^{\sigma }}\left({\beta }^{\rho }\right),\end{array}$$

(39)

with $\mu ,\nu ,\lambda$ on the left-hand-side ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$, $\mu ,\nu ,\rho$ on the right-hand side ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and with $\sigma$ on both sides ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only.

With these new modifications put in place, our hyperbolic system of equations governing the evolution of a perfect relativistic fluid on an arbitrary (and potentially dynamically-evolving) spacetime now consists of the following form of the energy conservation law:

$$\begin{array}{c}\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right)\right.\hfill \\ \left.+\frac{\partial }{\partial x^{\rho }}\left(\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)+P{v}^{\rho }\right)\right)\right]\\ \hfill =T^{00}\left({\beta }^{\mu }{\beta }^{\nu }{K}_{\mu \nu }-{\beta }^{\mu }\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)\right)+T^{0\left(\mu +1\right)}\left(-\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)+2{\beta }^{\nu }{K}_{\mu \nu }\right)+T^{\left(\mu +1\right)\left(\nu +1\right)}{K}_{\mu \nu },\end{array}$$

(40)

with $\mu ,\nu ,\rho$ now ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only for the whole equation, the following form of the momentum conservation law:

$$\begin{array}{c}\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\right)\right.\hfill \\ \left.+\frac{\partial }{\partial x^{\rho }}\left(\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)+P{\delta }_{\sigma }^{\rho }\right)\right)\right]\\ =T^{00}\left(\frac{1}{2}{\beta }^{\mu }{\beta }^{\nu }\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)-\alpha \frac{\partial }{\partial x^{\sigma }}\left(\alpha \right)\right)+T^{0\left(\mu +1\right)}{\beta }^{\nu }\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)\\ \hfill +\frac{1}{2}{T}^{\left(\mu +1\right)\left(\nu +1\right)}\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)+\frac{\rho h{v}_{\rho }}{\alpha \left(1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }\right)}\frac{\partial }{\partial x^{\sigma }}\left({\beta }^{\rho }\right),\end{array}$$

(41)

with $\mu ,\nu ,\rho ,\sigma$ now ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only for the whole equation, and the following form of the baryon number conservation law:

$$\begin{array}{c}\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right)\right.\hfill \\ \hfill \left.+\frac{\partial }{\partial x^{\rho }}\left(\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)\right)\right)\right]=0,\end{array}$$

(42)

with $\mu ,\nu ,\rho$ now ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only for the whole equation. The characteristic wave speeds of the fluid system can now be calculated by performing an eigendecomposition of its corresponding $\left(n+1\right)$-dimensional Jacobian matrices ${\mathbf{B}}^{\rho }$, namely:

$${\mathbf{B}}^{\rho }=\alpha \frac{\partial \left[\begin{array}{c}\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)+P{v}^{\rho }\\ \left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)+P{\delta }_{\sigma }^{\rho }\\ \left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)\end{array}\right]}{\partial \left[\begin{array}{c}\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\\ \frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\\ \frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\end{array}\right]},$$

(43)

with one Jacobian matrix ${\mathbf{B}}^{\rho }$ associated with each spatial coordinate direction $x^{\rho }$. As first calculated by Anile [48], Eulderink and Mellema [40], and later Banyuls et al. [19], the $\left(n+1\right)$ eigenvalues of each Jacobian matrix ${\mathbf{B}}^{\rho }$ can be grouped into those corresponding to the material wave speeds, i.e:

$${\lambda }_0^{\rho }=\alpha v^{\rho }-{\beta }^{\rho },$$

(44)

which have algebraic multiplicity n, and those corresponding to the acoustic wave speeds, i.e:

$${\lambda }_{\pm }^{\rho }=\frac{\alpha }{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }{c}_s^2}\left[v^{\rho }\left(1-c_s^2\right)\pm c_s\sqrt{\left(1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }\right)\left[{\gamma }^{\rho \rho }\left(1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }{c}_s^2\right)-v^{\rho }{v}^{\rho }\left(1-c_s^2\right)\right]}\right]-{\beta }^{\rho },$$

(45)

which each have algebraic multiplicity 1.

All that remains for us now is to consider the equations governing the dynamics of the discrete spacetime geometry itself. We start from the full Einstein field equations (including arbitrary stress-energy source terms), which are of a mixed hyperbolic-elliptic character:

$${ }^{\left(4\right)}{G}_{\mu \nu }+\Lambda g_{\mu \nu }={ }^{\left(4\right)}{R}_{\mu \nu }-\frac{1}{2}{ }^{\left(4\right)}R{g}_{\mu \nu }+\Lambda g_{\mu \nu }=8\pi T_{\mu \nu },$$

(46)

where ${ }^{\left(4\right)}{G}_{\mu \nu }$ is the spacetime Einstein tensor:

$${ }^{\left(4\right)}{G}_{\mu \nu }={ }^{\left(4\right)}{R}_{\mu \nu }-\frac{1}{2}{ }^{\left(4\right)}R{g}_{\mu \nu },$$

(47)

${ }^{\left(4\right)}{R}_{\mu \nu }$ is the spacetime Ricci tensor, obtained by contraction (${ }^{\left(4\right)}{R}_{\mu \nu }={ }^{\left(4\right)}{R}_{\mu \sigma \nu }^{\sigma }$) of the spacetime Riemann tensor ${ }^{\left(4\right)}{R}_{\sigma \mu \nu }^{\rho }$:

$${ }^{\left(4\right)}{R}_{\sigma \mu \nu }^{\rho }=\frac{\partial }{\partial x^{\mu }}\left({ }^{\left(4\right)}{\Gamma }_{\sigma \nu }^{\rho }\right)-\frac{\partial }{\partial x^{\nu }}\left({ }^{\left(4\right)}{\Gamma }_{\mu \sigma }^{\rho }\right)+{ }^{\left(4\right)}{\Gamma }_{\mu \lambda }^{\rho }{ }^{\left(4\right)}{\Gamma }_{\sigma \nu }^{\lambda }-{ }^{\left(4\right)}{\Gamma }_{\lambda \nu }^{\rho }{ }^{\left(4\right)}{\Gamma }_{\mu \sigma }^{\lambda },$$

(48)

${ }^{\left(4\right)}R$ is the spacetime Ricci scalar, obtained as the trace (i.e. ${ }^{\left(4\right)}R={ }^{\left(4\right)}{R}_{\mu }^{\mu }=g^{\mu \nu }{ }^{\left(4\right)}{R}_{\mu \nu }$) of the spacetime Ricci tensor ${ }^{\left(4\right)}{R}_{\mu \nu }$, and $\Lambda$ is the cosmological constant (essentially taken to be an arbitrary integration constant for our present purposes). We can decompose the ten independent components of the full Einstein field equations (assuming a four-dimensional spacetime manifold $\left(\mathcal{M},g\right)$, otherwise the number of independent components is equal to $\frac{1}{2}n\left(n+1\right)$ for n-dimensional spacetimes) into a system of six purely hyperbolic evolution equations (or $\frac{1}{2}n\left(n-1\right)$ evolution equations, for n-dimensional spacetimes) and a collection of four purely elliptic constraint equations (or n constraint equations, for n-dimensional spacetimes), with the latter (constraint) equations arising from the contracted Bianchi identities, which assert that the covariant divergence of the spacetime Einstein tensor ${ }^{\left(4\right)}{G}_{\mu \nu }$ must vanish identically:

$${ }^{\left(4\right)}{\nabla }_{\nu }{ }^{\left(4\right)}{G}^{\mu \nu }=\frac{\partial }{\partial x^{\nu }}\left({ }^{\left(4\right)}{G}^{\mu \nu }\right)+{ }^{\left(4\right)}{\Gamma }_{\nu \sigma }^{\mu }{ }^{\left(4\right)}{G}^{\sigma \nu }+{ }^{\left(4\right)}{\Gamma }_{\nu \sigma }^{\nu }{ }^{\left(4\right)}{G}^{\mu \sigma }=0.$$

(49)

In all of the above, $\mu ,\nu ,\rho ,\sigma ,\lambda$ range across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$. Upon performing an ADM decomposition of the spacetime metric, the hyperbolic evolution equations take the form:

$$\begin{array}{c}\frac{\partial }{\partial t}\left(K_{\nu }^{\mu }\right)=\alpha { }^{\left(3\right)}{R}_{\nu }^{\mu }-{ }^{\left(3\right)}{\nabla }_{\rho }\left({ }^{\left(3\right)}{\nabla }_{\nu }\alpha \right){\gamma }^{\rho \mu }+\alpha K{K}_{\nu }^{\mu }+{\beta }^{\rho }{ }^{\left(3\right)}{\nabla }_{\rho }{K}_{\nu }^{\mu }\hfill \\ \hfill +K_{\rho }^{\mu }{ }^{\left(3\right)}{\nabla }_{\nu }{\beta }^{\rho }-K_{\nu }^{\rho }{ }^{\left(3\right)}{\nabla }_{\rho }{\beta }^{\mu }-\alpha \left(8\pi T_{\left(\rho +1\right)\left(\nu +1\right)}{\gamma }^{\rho \mu }-4\pi T{\delta }_{\nu }^{\mu }\right)-\alpha \left(\frac{2\Lambda }{n-2}{\gamma }_{\rho \nu }\right){\gamma }^{\rho \mu },\end{array}$$

(50)

which then expand out to give:

$$\begin{array}{c}\frac{\partial }{\partial t}\left(K_{\nu }^{\mu }\right)=\alpha { }^{\left(3\right)}{R}_{\nu }^{\mu }-\left(\frac{\partial }{\partial x^{\rho }}\left(\frac{\partial }{\partial x^{\nu }}\left(\alpha \right)\right)-{ }^{\left(3\right)}{\Gamma }_{\rho \nu }^{\sigma }\left(\frac{\partial }{\partial x^{\sigma }}\left(\alpha \right)\right)\right){\gamma }^{\rho \mu }+\alpha K_{\nu }^{\mu }\hfill \\ +{\beta }^{\rho }\left(\frac{\partial }{\partial x^{\rho }}\left(K_{\nu }^{\mu }\right)+{ }^{\left(3\right)}{\Gamma }_{\rho \sigma }^{\mu }{K}_{\nu }^{\sigma }-{ }^{\left(3\right)}{\Gamma }_{\rho \nu }^{\sigma }{K}_{\sigma }^{\mu }\right)+K_{\rho }^{\mu }\left(\frac{\partial }{\partial x^{\nu }}\left({\beta }^{\rho }\right)+{ }^{\left(3\right)}{\Gamma }_{\nu \sigma }^{\rho }{\beta }^{\sigma }\right)\\ \hfill -K_{\nu }^{\rho }\left(\frac{\partial }{\partial x^{\rho }}\left({\beta }^{\mu }\right)+{ }^{\left(3\right)}{\Gamma }_{\rho \sigma }^{\mu }{\beta }^{\sigma }\right)-\alpha \left(8\pi T_{\left(\rho +1\right)\left(\nu +1\right)}{\gamma }^{\rho \mu }-4\pi T{\delta }_{\nu }^{\mu }\right)-\alpha \left(\frac{2\Lambda }{n-2}{\gamma }_{\rho \nu }\right){\gamma }^{\rho \mu },\end{array}$$

(51)

where ${ }^{\left(3\right)}{R}_{\mu \nu }$ is the spatial Ricci tensor, obtained by contraction (i.e. ${ }^{\left(3\right)}{R}_{\mu \nu }={ }^{\left(3\right)}{R}_{\mu \sigma \nu }^{\sigma }$) of the spatial Riemann tensor ${ }^{\left(3\right)}{R}_{\sigma \mu \nu }^{\rho }$:

$${ }^{\left(3\right)}{R}_{\sigma \mu \nu }^{\rho }=\frac{\partial }{\partial x^{\mu }}\left({ }^{\left(3\right)}{\Gamma }_{\sigma \nu }^{\rho }\right)-\frac{\partial }{\partial x^{\nu }}\left({ }^{\left(3\right)}{\Gamma }_{\mu \sigma }^{\rho }\right)+{ }^{\left(3\right)}{\Gamma }_{\mu \lambda }^{\rho }{ }^{\left(3\right)}{\Gamma }_{\sigma \nu }^{\lambda }-{ }^{\left(3\right)}{\Gamma }_{\lambda \nu }^{\rho }{ }^{\left(3\right)}{\Gamma }_{\mu \sigma }^{\lambda },$$

(52)

T is the trace of the stress-energy tensor (i.e. $T=g^{\mu \nu }{T}_{\mu \nu }$), and the indices of the extrinsic curvature tensor $K_{\mu \nu }$ and spatial Ricci tensor ${ }^{\left(3\right)}{R}_{\mu \nu }$ are raised and lowered using the spatial metric tensor ${\gamma }_{\mu \nu }$, and so, in particular, for the quantities $K_{\mu }^{\nu }$ and ${ }^{\left(3\right)}{R}_{\mu }^{\nu }$ (both in mixed-index form) appearing above, one has:

$$K_{\mu }^{\nu }={\gamma }_{\mu \sigma }{K}^{\sigma \nu }={\gamma }_{\mu \sigma }{\gamma }^{\lambda \nu }{K}_{\lambda }^{\sigma }={\gamma }^{\sigma \nu }{K}_{\mu \sigma },\mathrm{and}{R}_{\mu }^{\nu }={\gamma }_{\mu \sigma }{R}_{\sigma \nu }={\gamma }_{\mu \sigma }{\gamma }^{\lambda \nu }{R}_{\lambda }^{\sigma }={\gamma }^{\sigma \nu }{R}_{\mu \sigma }.$$

(53)

Finally, the elliptic constraint equations resulting from the contracted Bianchi identities may be decomposed into a timelike projection, yielding the Hamiltonian constraint equation:

$$\mathcal{H}={ }^{\left(3\right)}R+K^2-K_{\nu }^{\mu }{K}_{\mu }^{\nu }-16\pi {\alpha }^2{T}^{00}-2\Lambda =0$$

(54)

where $\left(3\right)R$ is the spatial Ricci scalar, obtained as the trace (i.e. ${ }^{\left(3\right)}R={ }^{\left(3\right)}{R}_{\mu }^{\mu }={\gamma }^{\mu \nu }{ }^{\left(3\right)}{R}_{\mu \nu }$) of the spatial Ricci tensor ${ }^{\left(3\right)}{R}_{\mu \nu }$, and into a collection of $\left(n-1\right)$ spacelike projections, yielding the momentum constraint equations:

$${\mathcal{M}}_{\mu }={ }^{\left(3\right)}{\nabla }_{\nu }{K}_{\mu }^{\mu }-{ }^{\left(3\right)}{\nabla }_{\mu }K-8\pi T_{\left(\mu +1\right)}^0=0,$$

(55)

which then expand out to give:

$${\mathcal{M}}_{\mu }=\frac{\partial }{\partial x^{\mu }}\left(K_{\mu }^{\nu }\right)+{ }^{\left(3\right)}{\Gamma }_{\nu \sigma }^{\nu }{K}_{\mu }^{\sigma }-{ }^{\left(3\right)}{\Gamma }_{\nu \mu }^{\sigma }{K}_{\sigma }^{\nu }-\frac{\partial }{\partial x^{\mu }}\left(K\right)-8\pi T_{\left(\mu +1\right)}^0=0.$$

(56)

In all of the above, $\mu ,\nu ,\rho ,\sigma ,\lambda$ range across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. Although Gravitas solves the elliptic constraint equations automatically when running numerical simulations (typically by means of an iterative solver), we have also validated the algorithms employed within this article by using violations of the constraint equations (and, in particular, the propagation of certain constraint-violating modes) as a means of measuring and quantifying the robustness of the relevant numerical schemes. Note that, analytically, due to the Einstein field equations, the Hamiltonian and momentum constraint equations on the spacetime are satisfied identically whenever the energy and momentum conservation equations on the stress-energy distribution are satisfied, and vice versa.

3. Numerical Validation: Special Relativistic Hydrodynamics

In order to render the discrete spacetime general relativistic hydrodynamics equations derived within the previous section in a form that is more directly amenable to explicit numerical solution, we begin by subdividing our overall n-dimensional spacetime $\left(\mathcal{M},g\right)$ into a collection of simply-connected n-dimensional submanifolds $\Omega \subseteq \mathcal{M}$ (known as “control volumes”, “computational cells”, or “nodes” within Gravitas’s terminology), each of which has a closed $\left(n-1\right)$-dimensional boundary $\partial \Omega$. We can then integrate over each of these submanifolds in turn, yielding:

$$\begin{array}{c}{\int }_{\Omega }\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right)\right]d\Omega \hfill \\ +{\int }_{\Omega }\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial x^{\rho }}\left(\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)+P{v}^{\rho }\right)\right)\right]d\Omega \\ \hfill ={\int }_{\Omega }\left[T^{00}\left({\beta }^{\mu }{\beta }^{\nu }{K}_{\mu \nu }-{\beta }^{\mu }\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)\right)+T^{0\left(\mu +1\right)}\left(-\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)+2{\beta }^{\nu }{K}_{\mu \nu }\right)+T^{\left(\mu +1\right)\left(\nu +1\right)}{K}_{\mu \nu }\right]d\Omega ,\end{array}$$

(57)

for the energy conservation equation;

$$\begin{array}{c}{\int }_{\Omega }\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\right)\right]d\Omega \hfill \\ +{\int }_{\Omega }\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial x^{\rho }}\left(\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)+P{\delta }_{\sigma }^{\rho }\right)\right)\right]d\Omega \\ ={\int }_{\Omega }\left[T^{00}\left(\frac{1}{2}{\beta }^{\mu }{\beta }^{\nu }\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)-\alpha \frac{\partial }{\partial x^{\sigma }}\left(\alpha \right)\right)+T^{0\left(\mu +1\right)}{\beta }^{\nu }\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)\right.\\ \hfill \left.+\frac{1}{2}{T}^{\left(\mu +1\right)\left(\nu +1\right)}\frac{\partial }{\partial x^{\sigma }}\left({\gamma }_{\mu \nu }\right)+\frac{\rho h{v}_{\rho }}{\alpha \left(1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }\right)}\frac{\partial }{\partial x^{\sigma }}\left({\beta }^{\rho }\right)\right]d\Omega ,\end{array}$$

(58)

for the momentum conservation equations; and:

$$\begin{array}{c}{\int }_{\Omega }\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right)\right]d\Omega \hfill \\ \hfill {\int }_{\Omega }\frac{1}{\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left[\frac{\partial }{\partial x^{\rho }}\left(\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\left(v^{\rho }-\frac{{\beta }^{\rho }}{\alpha }\right)\right)\right)\right]=0,\end{array}$$

(59)

for the baryon number continuity equation, with $\mu ,\nu ,\rho ,\sigma$ in all of the above ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. This so-called “weak” integral form of the conservation equations can then be solved directly using Gravitas’s hypergraph-based finite-volume numerical algorithms [35,36]. As a means of validating this numerical implementation, we first consider the simplified case of four-dimensional special relativistic hydrodynamics, in which the spacetime metric tensor $g_{\mu \nu }$ is simply the four-dimensional Minkowski metric ${\eta }_{\mu \nu }$, i.e. $g_{\mu \nu }={\eta }_{\mu \nu }=\mathrm{diag}\left(-1,1,1,1\right)$, the spatial metric tensor ${\gamma }_{\mu \nu }$ is the three-dimensional Euclidean metric ${\delta }_{\mu \nu }$, i.e. ${\gamma }_{\mu \nu }={\delta }_{\mu \nu }=\mathrm{diag}\left(1,1,1\right)$, and we select trivial gauge conditions in which the lapse function obeys the geodesic slicing condition (i.e. $\alpha =1$) and the shift vector obeys the normal coordinate conditions (i.e. $\beta =\mathbf{0}$). The energy conservation equation now reduces to:

$$\frac{\partial }{\partial t}\left(\frac{\rho h}{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)+\frac{\partial }{\partial x^{\rho }}\left(\left(\frac{\rho h}{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)v^{\rho }+P{v}^{\rho }\right)=0,$$

(60)

with the corresponding weak form:

$$\begin{array}{c}{\int }_{\Omega }\left[\frac{\partial }{\partial t}\left(\frac{\rho h}{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right]d\Omega \hfill \\ \hfill +{\int }_{\Omega }\left[\frac{\partial }{\partial x^{\rho }}\left(\left(\frac{\rho h}{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)v^{\rho }+P{v}^{\rho }\right)\right]d\Omega =0\end{array}$$

(61)

the momentum conservation equations reduce to:

$$\frac{\partial }{\partial t}\left(\frac{\rho h{v}_{\sigma }}{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)+\frac{\partial }{\partial x^{\rho }}\left(\left(\frac{\rho h{v}_{\sigma }}{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)v^{\rho }+P{\delta }_{\sigma }^{\rho }\right)=0,$$

(62)

with their corresponding weak forms being:

$${\int }_{\Omega }\left[\frac{\partial }{\partial t}\left(\frac{\rho h{v}_{\sigma }}{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\right]d\Omega +{\int }_{\Omega }\left[\frac{\partial }{\partial x^{\rho }}\left(\left(\frac{\rho h{v}_{\sigma }}{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)v^{\rho }+P{\delta }_{\sigma }^{\rho }\right)\right]d\Omega =0$$

(63)

and the baryon number continuity equation reduces to:

$$\frac{\partial }{\partial t}\left(\frac{\rho }{\sqrt{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)+\frac{\partial }{\partial x^{\rho }}\left(\left(\frac{\rho }{\sqrt{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)v^{\rho }\right)=0,$$

(64)

with its corresponding weak form being:

$${\int }_{\Omega }\left[\frac{\partial }{\partial t}\left(\frac{\rho }{\sqrt{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right]d\Omega +{\int }_{\Omega }\left[\left(\left(\frac{\rho }{\sqrt{1-{\delta }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)v^{\rho }\right)\right]d\Omega ,$$

(65)

and where $\mu ,\nu ,\rho ,\sigma$ in all of the above range, again, across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only.

Our primary objective here is to validate our conservative-to-primitive variable reconstruction algorithm. Since the conservative variables are the quantities that are actually evolved by our numerical algorithm, yet the primitive variables are the quantities that are required for the computation of the flux function on the next time step, it is necessary to perform a conversion from the conservative variables, consisting of the (rest) mass density of the fluid D measured by an observer moving in the normal direction $\mathbf{n}$

$$D=\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}=-J_{\mu }{n}^{\mu },$$

(66)

with $\mu ,\nu$ on the left-hand-side ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and $\mu$ on the right-hand-side ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$; the components of the momentum density covector $p_{\sigma }$ of the fluid measured by that normal observer:

$$p_{\sigma }=\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}=-T_{\mu \nu }{n}^{\mu }{\perp }_{\sigma }^{\nu },$$

(67)

with $\mu ,\nu$ on the left-hand-side ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, $\mu ,\nu$ on the right-hand-side ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$, and with $\sigma$ on both sides ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only; and the difference between the energy density E and (rest) mass density D of the fluid measured by that same normal observer:

$$E-D=\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P-\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}=T_{\mu \nu }{n}^{\mu }{n}^{\nu }-J_{\mu }{n}^{\mu },$$

(68)

with $\mu ,\nu$ on the left-hand-side ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only, and $\mu ,\nu$ on the right-hand-side ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$, to the primitive variables, consisting of the (rest) mass density of the fluid $\rho$, the components of the spatial velocity vector of the fluid measured by a normal observer $v^{\mu }$, and the hydrostatic pressure of the fluid P. In non-relativistic hydrodynamics, such a conversion can typically be performed purely algebraically, but in both special and general relativistic hydrodynamics (at least assuming a reasonably generic equation of state), the components of the momentum density covector $p_{\sigma }$ are not algebraically independent of one other due to the presence of the Lorentz factor $\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}$, and therefore no closed form expression for the primitive variables in terms of the conservative ones is known in general [49]. One notable exception to this is the case of stiff, ultra-relativistic fluids [50], in which such a closed form expression does exist (and, indeed, the existence of such a simplification is closely related to why Petrich, Shapiro and Teukolsky [15] were able to derive an analytic solution for accretion for the accretion of such fluids onto spinning black holes in general axial symmetry), although this particular property of the stiff, ultra-relativistic equation of state certainly does not generalize, as we shall discuss later on in this article.

For this reason, we choose to follow the approach proposed by Eulderink and Mellema [40], and apply an iterative, non-linear root-finding algorithm (namely the one-dimensional Newton-Raphson method) in order to approximate the roots to the following quartic polynomial in $\xi$ numerically:

$${\alpha }_4{\xi }^3\left(\xi -\eta \right)+{\alpha }_2{\xi }^2+{\alpha }_1\xi +{\alpha }_0=0,$$

(69)

where we have defined the variable $\xi$ and the constant $\eta$ to be given by:

$$\xi =\frac{\sqrt{-g_{\mu \nu }{T}^{0\mu }{T}^{0\nu }}}{\rho h{u}^0}=\frac{\left(\sqrt{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\mu }}\right)}^2}\right)\left(\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}{\rho h},$$

(70)

and:

$$\eta =\frac{2\rho u^0\left(\Gamma -1\right)}{\left(\sqrt{-g_{\mu \nu }{T}^{0\mu }{T}^{0\nu }}\right)\Gamma }=\frac{2\rho \left(\Gamma -1\right)}{\left(\sqrt{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}\right)\left(\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)\Gamma },$$

(71)

respectively, and where the coefficients ${\alpha }_4$, ${\alpha }_2$, ${\alpha }_1$ and ${\alpha }_0$ in front of the terms in the quartic are given by:

$${\alpha }_4=\frac{{\left(T^{00}\right)}^2}{g^{00}{g}_{\mu \nu }{T}^{0\mu }{T}^{0\nu }}-1=\frac{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2}{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}-1,$$

(72)

$$\begin{array}{c}{\alpha }_2=\left(\frac{\Gamma -2}{\Gamma }\right)\left(\frac{{\left(T^{00}\right)}^2}{g^{00}{g}_{\mu \nu }{T}^{0\mu }{T}^{0\nu }}-1\right)+1+\left(\frac{{\left(\rho u^0\right)}^2}{g_{\mu \nu }{T}^{0\mu }{T}^{0\nu }}\right){\left(\frac{\Gamma -1}{\Gamma }\right)}^2\hfill \\ =\left(\frac{\Gamma -2}{\Gamma }\right)\left(\frac{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2}{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}-1\right)+1\\ \hfill -\left(\frac{{\rho }^2}{\left({\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2\right)\left(1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }\right)}\right){\left(\frac{\Gamma -1}{\Gamma }\right)}^2,\end{array}$$

(73)

$${\alpha }_1=-\frac{2\rho u^0\left(\Gamma -1\right)}{\left(\sqrt{-g_{\mu \nu }{T}^{0\mu }{T}^{0\nu }}\right){\Gamma }^2}=-\frac{2\rho \left(\Gamma -1\right)}{\left(\sqrt{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}\right)\left(\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right){\Gamma }^2},$$

(74)

and:

$${\alpha }_0=-\frac{1}{{\Gamma }^2},$$

(75)

respectively. In all of the equations above, $\mu ,\nu$ on the left-hand-side range across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$, while $\mu ,\nu ,\sigma$ on the right-hand-side range across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. Since the variable $\xi$ itself represents the fluid quantity:

$$\xi =\frac{\left(\sqrt{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}\right)\left(\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}{\rho h},$$

(76)

we may proceed to use its approximate numerical value as a starting point for computing all of the other fluid quantities of interest, such as the new value of the Lorentz factor $W_{new}$:

$$\begin{array}{c}{W}_{new}=\frac{1}{2}\left(\frac{\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P}{\sqrt{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}}\right)\xi \hfill \\ \hfill \times \left(1+\sqrt{1+4\left(\frac{\Gamma -1}{\Gamma }\right)\left(\frac{1-\frac{\rho \xi }{\left(\sqrt{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}\right)\left(\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}}{\frac{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2{\xi }^2}{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}}\right)}\right),\end{array}$$

(77)

and, from it, the new value of the fluid (rest) mass density ${\rho }_{new}$:

$${\rho }_{new}=\frac{\rho }{\left(\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)W_{new}},$$

(78)

the new value of the specific relativistic enthalpy of the fluid $h_{new}$ (from which the new value of its hydrostatic pressure $P_{new}$ can then be recovered by means of the equation of state):

$$h_{new}=\frac{\left(\sqrt{{\left(\frac{\rho h}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}-P\right)}^2-{\left(\frac{\rho h{v}_{\sigma }}{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}^2}\right)\left(\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\right)}{\rho {\xi }^2},$$

(79)

and finally the components of the new spatial velocity covector of the fluid, as measured by a normal observer, $v_{\sigma }^{new}$:

$$v_{\sigma }^{new}=\frac{\rho h{v}_{\sigma }}{\left(1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }\right){\rho }_{new}{h}_{new}{W}_{new}^2}.$$

(80)

In all of the equations above, $\mu ,\nu ,\sigma$ range across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. As proved by Eulderink and Mellema [40], this quartic possesses exactly two real roots, with only one obeying the physicality condition $\xi >0$, and the iterative Newton-Raphson method is guaranteed to converge at least quadratically to this solution. In all of the above, we have assumed that the fluid obeys the ideal gas equation of state [51], with specific relativistic enthalpy given by

$$h=1+\frac{P}{\rho }\left(\frac{\Gamma }{\Gamma -1}\right),$$

(81)

where $\Gamma$ is the adiabatic exponent of the fluid, such that the local speed of sound $c_s$ is simply:

$$c_s=\sqrt{\frac{\Gamma P}{\rho \left(1+\left(\frac{P}{\rho }\right)\left(\frac{\Gamma }{\Gamma -1}\right)\right)}}.$$

(82)

We validate against a standard one-dimensional special relativistic shock tube problem, namely the mildly-relativistic blast wave problem proposed by Donat, Font, Ibáñez and Marquina [39], and we compare against the numerical solution of Del Zanna and Bucciantini [52], and, since Riemann problems in one-dimensional special relativistic hydrodynamics admit exact solutions, we compare also against the exact solution derived by Pons, Martí and Müller [53]. The fluid on the left-hand-side of the shock tube is of a high temperature and pressure, with $\rho =10$ and $P=13.3$; the fluid on the right-hand-side of the shock tube is of a low temperature and negligible pressure, with $\rho =1$ and $P=0$ (in some older papers using more unstable numerical methods, $P={10}^{-6}$ is chosen instead). Although this is a one-dimensional problem, since Gravitas’s hypergraph-based numerical algorithms work optimally in higher-dimensional geometries, we choose instead to evolve it as a three-dimensional problem in spherical symmetry, with the high temperature/high pressure fluid contained initially within a small spherical region in the center of the domain. It is then trivial to interpolate from a solution to this higher-dimensional spherically-symemtric problem back to a solution to the original one-dimensional Riemann problem. Since it will facilitate certain comparisons that we intend to perform later on in the case of general relativistic hydrodynamics in black hole spacetimes (which will be simulated in either spherical or axial symmetry), we choose to evolve this problem within both a spherically-symmetric Minkowksi spacetime parameterized by spherical polar coordinates $\left(t,r,\theta ,\varphi \right)$:

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-d{t}^2+d{r}^2+r^{2}d{\theta }^2+r^2{sin}^2\left(\theta \right)d{\varphi }^2,$$

(83)

with the initial hypersurface geometry:

$$d{l}^2={\gamma }_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=d{r}^2+r^{2}d{\theta }^2+r^2{sin}^2\left(\theta \right)d{\varphi }^2,$$

(84)

as well as within a standard rectangular Minkowski spacetime parameterized by Cartesian coordinates $\left(t,x,y,z\right)$:

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-d{t}^2+d{x}^2+d{y}^2+d{z}^2,\mathrm{with}d{l}^2={\gamma }_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=d{x}^2+d{y}^2+d{z}^2.$$

(85)

We run all simulations with a hypergraph resolution of 10,000 vertices. Internally, Gravitas uses an adaptive fourth-order Runge-Kutta algorithm [35,36] with a hypergraph rewriting/canonicalization algorithm based on [54]. An ideal gas equation of state with adiabatic exponent $\Gamma =\frac{5}{3}$ is assumed throughout. The initial ($t=0$) configurations of the domains in both cases are shown in Figure 1, with vertices colored based on fluid density and with vertex coordinates assigned using a two-dimensional projection of the spatial coordinates, yielding two-dimensional visualizations of the respective simulation domains. We can construct three-dimensional visualizations of the domains by also assigning a third vertex coordinate based on the fluid density, as shown in Figure 2. Finally, to give an indication of the hypergraph topology produced by Gravitas’s adaptive hypergraph refinement algorithm (applied here as a preconditioning step for the initial data), we show the initial hypergraphs without any vertex coordinate information assigned in Figure 3. We shall make use of all three modes of visualization throughout the remainder of this article. The exact solution to the mildly-relativistic blast wave Riemann problem is known to consist of three waves: a slow-moving rarefaction wave, a contact discontinuity, and a fast-moving shock wave. In Figure 4, Figure 5 and Figure 6 (showing the solution at coordinate time $t=0.4$ with two-dimensional spatial coordinates, three-dimensional spatial and fluid density coordinates, and no vertex coordinates, respectively), we see that all three waves are resolved correctly in both the spherically-symmetric and rectangular cases.

4. Black Hole Accretion in Discrete Spacetime

Following the recent historical exposition of Aguayo-Ortiz, Tejeda, Sarbach and López-Cámara [55], we note that Bondi [11] originally analyzed the case of an infinite, spherically-symmetric distribution of ideal gas with adiabatic exponent $\Gamma$, with initially uniform density $\rho$ and pressure P, accreting radially onto a compact object of mass M in Newtonian gravity. Under these assumptions, the steady-state accretion flow in spherical polar coordinates $\left(t,r,\theta ,\varphi \right)$ is governed by the continuity and radial Euler equations, namely:

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\rho v_r\right)=0,\mathrm{and}{v}_r\frac{d{v}_r}{dr}+\frac{1}{\rho }\frac{dP}{dr}+\frac{M}{r^2}=0,$$

(86)

respectively, where $v_r$ designates the radial component of the fluid velocity, i.e:

$$v_r=\left|\frac{dr}{dt}\right|.$$

(87)

The initial conditions of the fluid distribution are given in terms of its density and pressure at infinite radial distance, i.e. ${\rho }_{\infty }$ and $P_{\infty }$, respectively, with the fluid assumed to be at rest at infinite radial distance (i.e. ${\left(v_r\right)}_{\infty }=0$) and with the (local) adiabatic speed of sound $c_s$ given by the following partial derivative assuming fixed internal energy $\epsilon$:

$$c_s=\sqrt{{\left.\frac{\partial P}{\partial \rho }\right|}_{\epsilon }}.$$

(88)

Subject to the additional assumption that the flow is transonic, i.e. that there exists a radius $r=r_{trans}$ such that the radial velocity $v_r=v_{trans}$ is equal to the (local) adiabatic speed of sound, which can be calculated to be:

$$r_{trans}=\frac{M}{2v_{trans}^2 },\mathrm{where}{v}_{trans}=c_s=c_{\infty }\sqrt{\frac{2}{5-3\Gamma }},$$

(89)

and where $c_{\infty }$ is the adiabatic speed of sound at infinite radial distance, Bondi showed that there exists a unique analytic solution to the continuity and radial Euler equations which maximizes the rate at which mass is accreted onto the compact object, namely:

$$\frac{dM}{dt}=\pi {\left(\frac{2}{5-3\Gamma }\right)}^{\frac{5-3\Gamma }{2\left(\Gamma -1\right)}}{M}^2\frac{{\rho }_{\infty }}{c_{\infty }^3}.$$

(90)

Bondi’s purely Newtonian analysis was subsequently extended by Michel [13] to the case of an infinite, spherically-symmetric distribution of ideal gas, again with an initially uniform density $\rho$ and pressure P, accreting radially onto a static, uncharged and non-rotating black hole (as described by a Schwarzschild geometry) of mass M in general relativity. The Schwarzschild metric is taken to be given in Schwarzschild/spherical polar coordinates $\left(t,r,\theta ,\varphi \right)$ by [56,57]:

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-\left(1-\frac{2M}{r}\right)d{t}^2+{\left(1-\frac{2M}{r}\right)}^{-1}d{r}^2+r^2\left(d{\theta }^2+{sin}^2\left(\theta \right)d{\varphi }^2\right),$$

(91)

and, within such a spacetime, the spherical symmetry of the problem allows us to reformulate the energy-momentum conservation equations as a single (radial) ordinary differential equation:

$${\nabla }_{\nu }{T}^{\mu \nu }=\frac{\partial }{\partial x^{\nu }}\left(T^{\mu \nu }\right)+{\Gamma }_{\nu \sigma }^{\mu }{T}^{\sigma \nu }+{\Gamma }_{\nu \sigma }^{\nu }{T}^{\mu \sigma }=0,\Rightarrow \frac{d}{dr}\left(r^2\rho h{u}_t{u}^r\right)=0$$

(92)

and the baryon number continuity equation may be reformulated similarly:

$${\nabla }_{\mu }\left(\rho u^{\mu }\right)=\frac{\partial }{\partial x^{\mu }}\left(\rho u^{\mu }\right)+{\Gamma }_{\mu \sigma }^{\mu }\left(\rho u^{\sigma }\right)=0,\Rightarrow \frac{d}{dr}\left(r^2\rho u^r\right)=0,$$

(93)

where h is, as usual, the specific relativistic enthalpy of the fluid:

$$h=1+\epsilon \left(\rho ,P\right)+\frac{P}{\rho }.$$

(94)

In all other respects, the specification of the initial conditions of the fluid in terms of its density ${\rho }_{\infty }$, pressure $P_{\infty }$ and specific relativistic enthalpy $h_{\infty }$ at infinite radial distance (with the fluid again assumed to be at rest at this point, i.e. ${\left(u_r\right)}_{\infty }=0$) is directly analogous to the Newtonian case, with the (local) adiabatic speed of sound $c_s$ now given by the following partial derivative, assuming a fixed hydrostatic pressure P:

$$c_s=\sqrt{\frac{\rho }{h}{\left(\left.\frac{\partial h}{\partial \rho }\right)\right|}_P}.$$

(95)

Once again, we make the assumption that the flow is transonic, with the radial velocity $u_r={\left(u_r\right)}_{trans}$ at the transonic radius $r=r_{trans}$ now being such that the norm of the spatial velocity of the fluid is measured by any local static observer as being equal to the (local) adiabatic speed of sound:

$$r_{trans}=\frac{M}{2{\left(u_r\right)}_{trans}},\mathrm{where}{\left(u_r\right)}_{trans}=\sqrt{\frac{\frac{1}{3}\left(\frac{h_{trans}^2}{h_{\infty }^2}-1\right)}{\left(\frac{h_{trans}^2}{h_{\infty }^2}\right)}},$$

(96)

where the transonic value of the specific relativistic enthalpy $h_{trans}$ is calculated as:

$$h_{trans}=2h_{\infty }\sqrt{\Gamma -\frac{2}{3}}sin\left(\frac{1}{3}arccos\left(\frac{3\left(\Gamma -1\right)}{2h_{\infty }{\left(\sqrt{\Gamma -\frac{2}{3}}\right)}^3}\right)\right).$$

(97)

Once again, there exists a unique analytic solution to the energy-momentum and baryonic number conservation equations within this setting, such that the flow satisfies the steady-state condition and regularity across the event horizon of the black hole is preserved (as proved by Chaverra and Sarbach [58], and later Chaverra, Mach and Sarbach [59]), with the accretion rate of mass onto the black hole given by:

$$\frac{dM}{dt}=\pi {\left(\frac{h_{trans}}{h_{\infty }}\right)}^{\frac{3\Gamma -2}{\Gamma -1}}{\left(\frac{\sqrt{\frac{1}{3}{\left(\frac{h_{trans}^2}{h_{\infty }^2}-1\right)}^2}}{c_{\infty }}\right)}^{\frac{5-3\Gamma }{\Gamma -1}},$$

(98)

where $c_{\infty }$ is, again, the adiabatic speed of sound at infinite radial distance.

Font, Ibáñez and Papadopoulos [41] have stressed the importance of using “horizon-adapted” coordinates, i.e. coordinate systems (such as Kerr-Schild coordinates) which remain regular at the black hole event horizon, when performing black hole accretion studies, so as to avoid unphysical fluid behavior near the horizon resulting from coordinate divergences. For this reason, we shall perform our radial accretion simulations onto static, uncharged, non-rotating black holes (i.e. Schwarzschild black holes) expressed in both the Schwarzschild/spherical polar coordinate system $\left(t,r,\theta ,\varphi \right)$, which is not horizon-adapted, with the metric of the initial spacelike hypersurface given by:

$$d{l}^2={\gamma }_{\mu \nu }d{x}^{\mu }d{x}^{\nu }={\left(1-\frac{2M}{r}\right)}^{-1}d{r}^2+r^2\left(d{\theta }^2+si{n}^2\left(\theta \right)d{\varphi }^2\right),$$

(99)

and the Kerr-Schild/Cartesian coordinate system $\left(t,x,y,z\right)$ [60], which is horizon-adapted, with the metric of the initial spacelike hypersurface given by:

$$d{l}^2={\gamma }_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=d{x}^2+d{y}^2+d{z}^2+\mathcal{F}{\left(l_{\mu }d{x}^{\mu }\right)}^2,$$

(100)

where:

$$\mathcal{F}=\frac{2M}{r},\mathrm{and}{l}_{\mu }d{x}^{\mu }=\frac{x}{r}dx+\frac{y}{r}dy+\frac{z}{r}dz,$$

(101)

and with r being the usual radial coordinate:

$$r=\sqrt{x^2+y^2+z^2},$$

(102)

in order that we be able to compare the hydrodynamic results obtained across the two coordinate schemes. As in the special relativistic hydrodynamics cases from before, we run all simulations presented in this section with a hypergraph resolution of 10,000 vertices, using an ideal gas equation of state:

$$h=1+\frac{P}{\rho }\left(\frac{\Gamma }{\Gamma -1}\right),\mathrm{and}{c}_s=\sqrt{\frac{\Gamma P}{\rho \left(1+\left(\frac{P}{\rho }\right)\left(\frac{\Gamma }{\Gamma -1}\right)\right)}},$$

(103)

with adiabatic exponent $\Gamma =\frac{5}{3}$. The geometries of the initial ($t=0$) hypersurface configurations in both the Schwarzschild/spherical polar coordinate system and the Kerr-Schild/Cartesian coordinate system are shown in Figure 7, with vertices colored based on the value of the extrinsic curvature tensor, and with vertex coordinates assigned using a two-dimensional spatial projection. Three-dimensional visualizations are shown in Figure 8, with the third vertex coordinate assigned based on extrinsic curvature. Finally, coordinate-free representations of the pure hypergraph topologies are shown in Figure 9.

In order to evolve this spatial metric (together with the fluid variables defined on top of it) forwards in time using Gravitas, we must first select an appropriate set of gauge conditions. For the lapse function $\alpha$, we choose to use the maximal slicing condition initially developed by Lichnerowicz [61] and later developed into a directly usable form by York [44]:

$${ }^{\left(3\right)}\Delta \alpha =\alpha K^{\mu \nu }{K}_{\mu \nu }-\frac{\partial }{\partial t}\left(K\right),$$

(104)

where ${ }^{\left(3\right)}\Delta$ denotes the connection Laplacian on spacelike hypersurfaces, defined for arbitrary scalar fields $\varphi$ as:

$${ }^{\left(3\right)}\Delta \varphi ={ }^{\left(3\right)}{\nabla }^{\mu }\left({ }^{\left(3\right)}{\nabla }_{\mu }\varphi \right)={\gamma }^{\mu \sigma }{ }^{\left(3\right)}{\nabla }_{\sigma }\left({ }^{\left(3\right)}{\nabla }_{\mu }\varphi \right)={\gamma }^{\mu \sigma }\left(\frac{\partial }{\partial x^{\sigma }}\left(\frac{\partial }{\partial x^{\mu }}\left(\varphi \right)\right)-{ }^{\left(3\right)}{\Gamma }_{\sigma \mu }^{\lambda }\left(\frac{\partial }{\partial x^{\lambda }}\left(\varphi \right)\right)\right),$$

(105)

which, using the contraction properties of the (spatial) Christoffel symbols ${ }^{\left(3\right)}{\Gamma }_{\sigma \mu }^{\lambda }$, becomes:

$${ }^{\left(3\right)}\Delta \varphi =\frac{1}{\sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left(\frac{\partial }{\partial x^{\mu }}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left({\gamma }^{\mu \nu }\frac{\partial }{\partial x^{\nu }}\left(\varphi \right)\right)\right)\right),$$

(106)

leading to the following explicit form of the maximal slicing condition, with the lapse function $\alpha$ being treated as a scalar field defined over spacelike hypersurfaces:

$$\frac{1}{\sqrt{det\left({\gamma }_{\mu \nu }\right)}}\left(\frac{\partial }{\partial x^{\mu }}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left({\gamma }^{\mu \nu }\frac{\partial }{\partial x^{\mu }}\left(\alpha \right)\right)\right)\right)=\alpha K^{\mu \nu }{K}_{\mu \nu }-\frac{\partial }{\partial t}\left(K\right).$$

(107)

Note, as before, that the indices of the extrinsic curvature tensor $K_{\mu \nu }$ are raised and lowered using the spatial metric tensor ${\gamma }_{\mu \nu }$, and so, in particular, for the contravariant form $K^{\mu \nu }$, one has:

$$K^{\mu \nu }={\gamma }^{\mu \sigma }{K}_{\sigma }^{\nu }={\gamma }^{\sigma \nu }{K}_{\sigma }^{\mu }={\gamma }^{\mu \sigma }{\gamma }^{\lambda \nu }{K}_{\sigma \lambda }.$$

(108)

The maximal slicing condition seeks to maximize the spatial volume of each spacelike hypersurface by reducing the evolution rate in high-curvature regions and increasing it in low-curvature regions, thus equipping it with highly favorable singularity-avoidance properties that make it ideal for simulating black hole spacetimes. For the shift vector $\beta$, we choose to use the minimal distortion coordinate conditions of Smarr and York [62], subsequently adapted into the form employed here by Brady, Creighton and Thorne [63]:

$$\begin{array}{c}{ }^{\left(3\right)}{\nabla }^{\mu }\left({ }^{\left(3\right)}{\nabla }_{\mu }{\beta }^{\nu }\right)+{ }^{\left(3\right)}{\nabla }^{\nu }\left({ }^{\left(3\right)}{\nabla }_{\mu }{\beta }^{\mu }\right)-2{ }^{\left(3\right)}{\nabla }_{\mu }\left(\alpha K^{\mu \nu }\right)\hfill \\ \hfill ={\gamma }^{\mu \sigma }{ }^{\left(3\right)}{\nabla }_{\sigma }\left({ }^{\left(3\right)}{\nabla }_{\mu }{\beta }^{\nu }\right)+{\gamma }^{\nu \sigma }{ }^{\left(3\right)}{\nabla }_{\sigma }\left({ }^{\left(3\right)}{\nabla }_{\mu }{\beta }^{\mu }\right)-2{ }^{\left(3\right)}{\nabla }_{\mu }\left(\alpha K^{\mu \nu }\right)=0,\end{array}$$

(109)

which expands out to give:

$$\begin{array}{c}{\gamma }^{\mu \sigma }\left(\frac{\partial }{\partial x^{\sigma }}\left(D_{\mu }^{\nu }\right)+{ }^{\left(3\right)}{\Gamma }_{\sigma \lambda }^{\nu }{D}_{\mu }^{\lambda }-{ }^{\left(3\right)}{\Gamma }_{\sigma \mu }^{\lambda }{D}_{\lambda }^{\mu }\right)+{\gamma }^{\nu \sigma }\left(\frac{\partial }{\partial x^{\sigma }}\left(D_{\mu }^{\mu }\right)+{ }^{\left(3\right)}{\Gamma }_{\sigma \lambda }^{\mu }{D}_{\mu }^{\lambda }-{ }^{\left(3\right)}{\Gamma }_{\sigma \mu }^{\lambda }{D}_{\lambda }^{\mu }\right)\hfill \\ \hfill -2\left(\frac{\partial }{\partial x^{\mu }}\left(\alpha K^{\mu \nu }\right)+{ }^{\left(3\right)}{\Gamma }_{\mu \sigma }^{\mu }\left(\alpha K^{\sigma \nu }\right)+{ }^{\left(3\right)}{\Gamma }_{\mu \sigma }^{\nu }\left(\alpha K^{\mu \sigma }\right)\right)=0,\end{array}$$

(110)

with the rank-2 tensor $D_{\mu }^{\nu }$ consisting of (spatial) covariant derivatives of the shift vector components ${\beta }^{\nu }$:

$$D_{\mu }^{\nu }={ }^{\left(3\right)}{\nabla }_{\nu }{\beta }^{\mu }=\frac{\partial }{\partial x^{\mu }}\left({\beta }^{\nu }\right)+{ }^{\left(3\right)}{\Gamma }_{\mu \sigma }^{\nu }{\beta }^{\sigma }.$$

(111)

The minimal distortion coordinate conditions seek to minimize the distortion (or “strain”) in the spatial coordinates as one evolves from one hypersurface to the next, which makes it preferable for the case of hydrodynamics simulations, in which one generally wishes for the spatial coordinate system within which the fluid is evolved to remain as consistent as possible between time steps. In all of the above, $\mu ,\nu ,\sigma ,\lambda$ range across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. We initialize our simulation with the (non-dimensional) gas temperature at infinite radial distance set to be ${\Theta }_{\infty }=0.1$, from which the density and pressure at infinite radial distance (i.e. ${\rho }_{\infty }$ and $P_{\infty }$, respectively) can be computed using the relation ${\Theta }_{\infty }=\frac{P_{\infty }}{{\rho }_{\infty }}$. In Figure 10, Figure 11 and Figure 12 (showing the solution at coordinate time $t=100M$ with two-dimensional spatial coordinates, three-dimensional spatial and fluid density coordinates, and no vertex coordinates, respectively), we see that the fluid eventually evolves to a steady-state configuration with a high-density spherical accretion region surrounding the black hole event horizon, with no substantive difference between the solutions seen in the Schwarzschild/spherical polar and Kerr-Schild/Cartesian coordinate systems. The rates of mass/energy accretion onto the black hole are found to be slightly lower than in the analytic solution of Michel [13], with this discrepancy vanishing in the limit as the discretization scale goes to zero.

The analyses of Bondi [11] and Michel [13] described above both made extensive use of the spherical symmetry of the accretion problem. However, the treatment of an infinite, initially spherically-symmetric distribution of ideal gas, with initially uniform density $\rho$ and pressure P, accreting radially onto an uncharged but spinning black hole (as described by a Kerr geometry) of mass M and spin J in full general relativity requires breaking this spherical symmetry, and replacing it with a more general axially-symmetric spacetime geometry. The Kerr metric [64] is taken to be given in Boyer-Lindquist/oblate spheroidal coordinates [65]$\left(t,r,\theta ,\varphi \right)$ by:

$$\begin{array}{c}d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-\left(1-\frac{2M}{\left(r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)\right)}\right)d{t}^2+\left(\frac{r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)}{r^2-2M+{\left(\frac{J}{M}\right)}^2}\right)d{r}^2\hfill \\ +\left(r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)\right)d{\theta }^2+\left(r^2+{\left(\frac{J}{M}\right)}^2+\frac{2J^2{sin}^2\left(\theta \right)}{M\left(r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)\right)}\right){sin}^2\left(\theta \right)d{\varphi }^2\\ \hfill -\left(\frac{4J{sin}^2\left(\theta \right)}{r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)}\right)dtd\varphi ,\end{array}$$

(112)

and, following Petrich, Shapiro and Teukolsky [15], we assume an ultra-relativistic equation of state in which the rest-mass energy of the fluid is negligible when compared to its internal energy:

$$P=\left(\Gamma -1\right)\rho ,$$

(113)

and we assume, moreover, that the fluid is stiff in the sense that $\Gamma =2$ and therefore $P=\rho$ identically. Subject to the additional assumption that the flow obtains a steady-state configuration that is non-rotational, the resulting fluid equations can be expressed purely in terms of the gradient of a certain stream function (or scalar potential) $\Phi$:

$$h{u}_{\mu }=\frac{\partial }{\partial x^{\mu }}\left(\Phi \right),$$

(114)

which, due to the normalization convention $u_{\mu }{u}^{\mu }=1$ for the spacetime velocity vector $\mathbf{u}$, implies that the specific relativistic enthalpy h can be written purely in terms of the stream function gradient:

$$h=\sqrt{-\frac{\partial }{\partial x^{\mu }}\left(\Phi \right)g^{\mu \sigma }\frac{\partial }{\partial x^{\sigma }}\left(\Phi \right)}.$$

(115)

Substituting the stream function gradient equation into the baryon number conservation equation:

$${ }^{\left(4\right)}{\nabla }_{\mu }\left(\rho u^{\mu }\right)=\frac{\partial }{\partial x^{\mu }}\left(\rho u^{\mu }\right)+{ }^{\left(4\right)}{\Gamma }_{\mu \sigma }^{\mu }\left(\rho u^{\sigma }\right)=0,$$

(116)

yields the following harmonic equation for $\Phi$:

$${ }^{\left(4\right)}\Delta \Phi ={ }^{\left(4\right)}{\nabla }^{\mu }\left({ }^{\left(4\right)}{\nabla }_{\mu }\Phi \right)=g^{\mu \sigma }{ }^{\left(4\right)}{\nabla }_{\sigma }\left({ }^{\left(4\right)}{\nabla }_{\mu }\Phi \right)=0,$$

(117)

where ${ }^{\left(4\right)}\Delta$ denotes the connection Laplacian on spacetime, which expands out to give:

$$g^{\mu \sigma }\left(\frac{\partial }{\partial x^{\sigma }}\left(\frac{\partial }{\partial x^{\mu }}\left(\Phi \right)\right)-{ }^{\left(4\right)}{\Gamma }_{\sigma \mu }^{\lambda }\left(\frac{\partial }{\partial x^{\lambda }}\left(\Phi \right)\right)\right)=0,$$

(118)

which, using the contraction properties of the spacetime Christoffel symbols ${ }^{\left(4\right)}{\Gamma }_{\sigma \mu }^{\lambda }$, becomes:

$$\frac{1}{\sqrt{-det\left(g_{\mu \nu }\right)}}\left(\frac{\partial }{\partial x^{\sigma }}\left(\sqrt{-det\left(g_{\mu \nu }\right)}\left(g^{\mu \sigma }\frac{\partial }{\partial x^{\mu }}\left(\Phi \right)\right)\right)\right)=0.$$

(119)

In all of the above, $\mu ,\nu ,\lambda ,\sigma$ range across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$.

By imposing the boundary condition that the fluid should be at rest at infinite radial distance (i.e. $v_{\infty }^{\mu }=0$), and also that the fluid distribution should be uniform at this distance (i.e. that ${\rho }_{\infty }$, $P_{\infty }$ and $h_{\infty }$ should all be constant), there exists an analytic solution to this equation for the stream function, as derived by Aguayo-Ortiz, Sarbach and Tejeda [66], namely:

$$\Phi =h_{\infty }\left[-t+2Mlog\left(\frac{r-M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}}{2\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}}\right)\right].$$

(120)

This, in turn, yields the following values for the timelike, radial, polar and azimuthal projections of the spacetime velocity vector $\mathbf{u}$:

$$\frac{h}{h_{\infty }}{u}^t=1+\frac{2Mr}{r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)}\left(\frac{r+M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}}{r-M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}}\right),$$

(121)

$$\frac{h}{h_{\infty }}{u}^r=-\frac{2M\left(M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}\right)}{r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)},$$

(122)

$$\frac{h}{h_{\infty }}{u}^{\theta }=0,$$

(123)

and:

$$\frac{h}{h_{\infty }}{u}^{\varphi }=\frac{2Jr}{\left(r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)\right)\left(r-M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}\right)},$$

(124)

respectively, as well as the following relation for the fluid (rest) mass density $\rho$ and/or specific relativistic enthalpy h:

$$\frac{\rho }{{\rho }_{\infty }}=\frac{h}{h_{\infty }}=\sqrt{1+\left(\frac{2M}{r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)}\right)\left(\frac{r\left(r+M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}\right)+2M\left(M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}\right)}{r-M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}}\right)}.$$

(125)

One consequently recovers the analytic solution of Petrich, Shapiro and Teukolsky [15] for the accretion rate of mass onto the spinning black hole:

$$\frac{dM}{dt}=8\pi M\left(M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}\right){\rho }_{\infty }=4\pi \left({\left(M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}\right)}^2+{\left(\frac{J}{M}\right)}^2\right){\rho }_{\infty }.$$

(126)

As in the Schwarzschild case described previously, we perform our radial accretion simulations onto uncharged, spinning black holes (i.e. Kerr black holes) expressed in both the Boyer-Lindquist/oblate spheroidal coordinate system $\left(t,r,\theta ,\varphi \right)$, which is not horizon-adapted, with the metric of the initial spacelike hypersurface given by:

$$\begin{array}{c}d{l}^2={\gamma }_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=\left(\frac{r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)}{r^2-2Mr+{\left(\frac{J}{M}\right)}^2}\right)d{r}^2+\left(r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)\right)d{\theta }^2\hfill \\ \hfill +\left(r^2+{\left(\frac{J}{M}\right)}^2+\frac{2J^2{sin}^2\left(\theta \right)}{M\left(r^2+{\left(\frac{J}{M}\right)}^2{cos}^2\left(\theta \right)\right)}\right){sin}^2\left(\theta \right)d{\varphi }^2,\end{array}$$

(127)

and the Kerr-Schild/Cartesian coordinate system $\left(t,x,y,z\right)$, which is horizon-adapted, with the metric of the initial spacelike hypersurface given by:

$$d{l}^2={\gamma }_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=d{x}^2+d{y}^2+d{z}^2+\mathcal{F}{\left(l_{\mu }d{x}^{\mu }\right)}^2,$$

(128)

where:

$$\mathcal{F}=\frac{2M{r}^2}{r^4+{\left(\frac{J}{M}\right)}^2{z}^2},\mathrm{and}{l}_{\mu }d{x}^{\mu }=\frac{z}{r}dz+\frac{r}{r^2+{\left(\frac{J}{M}\right)}^2}\left(xdx+ydy\right)-\frac{{\left(\frac{J}{M}\right)}^2}{r^2+{\left(\frac{J}{M}\right)}^2}\left(xdy-ydx\right),$$

(129)

with r no longer being the usual radial coordinate, but rather being defined implicitly as a (positive, real) solution to the following algebraic equation:

$$\frac{x^2+y^2}{r^2+{\left(\frac{J}{M}\right)}^2}+\frac{z^2}{r^2}=1.$$

(130)

We initialize our simulations, as before, with an ideal gas equation of state with adiabatic exponent $\Gamma =\frac{5}{3}$, and the non-dimensional gas temperature at infinite radial distance (which determines both ${\rho }_{\infty }$ and $P_{\infty }$) set to be ${\Theta }_{\infty }=0.1$. We begin by considering a black hole with only a modest spin value of $J=0.6M$, and in Figure 13, Figure 14 and Figure 15 (showing the solution at coordinate time $t=100M$ with, as before, two-dimensional spatial coordinates, three-dimensional spatial and fluid density coordinates, and no vertex coordinates, respectively), we see that the fluid in this case evolves to a steady-state configuration with a single high-density “swirl” surrounding the black hole horizon. For a rapidly-spinning black hole with a spin value of $J=0.9M$, as shown in Figure 16, Figure 17 and Figure 18, we see that this high-density “swirl” effectively splits into two distinct “arms”, while for a black hole spinning close to the threshold of extremality with $J=0.99M$, as shown in Figure 19,–Figure 21, we see splitting of the “swirl” into three “arms” instead. We see evidence of some slight boundary effects around the edges and corners of the domain in Kerr-Schild/Cartesian coordinates, and no evidence of unphysical fluid behavior close to the horizon in Boyer-Lindquist/oblate spheroidal coordinates (which we attribute to our robust and singularity-avoiding choice of gauge). The rates of mass/energy accretion onto the spinning black holes are, again, found to be slightly lower than in the analytic solution of Petrich, Shapiro and Teukolsky [15], with this discrepancy vanishing in the limit as the discretization scale and the black hole spin both go to zero.

In order to conduct a more quantitative analysis of these simulation results, we begin by extracting the rates of mass/energy and momentum accretion onto the black hole in each case. Following the approach of Petrich, Shapiro, Stark and Teukolsky [50], the total fluid (rest) mass M contained within a given simply-connected spatial volume V of dimension $\left(n-1\right)$, with a closed $\left(n-2\right)$-dimensional boundary $\partial V$, can be computed by means of the following integral of the fluid (rest) mass density:

$$M={\int }_V\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)dV,$$

(131)

which, upon application of the Leibniz integral rule, allows us to write the mass accretion rate within that volume as:

$$\frac{dM}{dt}={\int }_V\frac{\partial }{\partial t}\left(\sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\right)\right)dV,$$

(132)

with $\mu ,\nu$ in the above ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. Since, by the definitions of the Lorentz factor $\alpha u^0$, the (rest) mass current vector $J^{\mu }$ and the spacetime metric determinant $det\left(g_{\mu \nu }\right)$, we have:

$$\alpha u^0=-u_{\left(\mu +1\right)}{n}^{\left(\mu +1\right)}=\frac{1}{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}},J^{\mu }=\rho u^{\mu },\mathrm{and}\sqrt{-det\left(g_{\mu \nu }\right)}=\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)},$$

(133)

respectively, with $\mu$ ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$ in $u_{\left(\mu +1\right)}{n}^{\left(\mu +1\right)}$, $J^{\mu }$ and $\rho u^{\mu }$, and with $\mu ,\nu$ ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only everywhere else, we can rewrite these two integrals as:

$$M={\int }_V{J}^0\sqrt{-det\left(g_{\mu \nu }\right)}dV,\mathrm{and}\mathrm{therefore}\frac{dM}{dt}={\int }_V\frac{\partial }{\partial t}\left(J^0\sqrt{-det\left(g_{\mu \nu }\right)}\right)dV,$$

(134)

respectively. The time derivative in the latter integral can now be rewritten as a difference between a spacetime covariant derivative and a spatial partial derivative, namely:

$${\int }_V\frac{\partial }{\partial t}\left(J^0\sqrt{-det\left(g_{\mu \nu }\right)}\right)dV={\int }_V\left[{ }^{\left(4\right)}{\nabla }_{\mu }\left(J^{\mu }\sqrt{-det\left(g_{\mu \nu }\right)}\right)-\frac{\partial }{\partial x^{\nu }}\left(J^{\left(\nu +1\right)}\sqrt{-det\left(g_{\mu \nu }\right)}\right)\right]dV,$$

(135)

with $\mu$ ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$ and $\nu$ ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. However, the covariant derivative term vanishes by virtue of the conservation of baryon number:

$$\begin{array}{c}{ }^{\left(4\right)}{\nabla }_{\mu }{J}^{\mu }=\frac{\partial }{\partial x^{\mu }}\left(J^{\mu }\right)+{ }^{\left(4\right)}{\Gamma }_{\mu \sigma }^{\mu }{J}^{\sigma }=0,\hfill \\ \hfill \Rightarrow { }^{\left(4\right)}{\nabla }_{\mu }\left(J^{\mu }\sqrt{-det\left(g_{\mu \nu }\right)}\right)=\sqrt{-det\left(g_{\mu \nu }\right)}\left({ }^{\left(4\right)}{\nabla }_{\mu }{J}^{\mu }\right)=\sqrt{-det\left(g_{\mu \nu }\right)}\left(\frac{\partial }{\partial x^{\mu }}\left(J^{\mu }\right)+{ }^{\left(4\right)}{\Gamma }_{\mu \sigma }^{\mu }{J}^{\sigma }\right)=0,\end{array}$$

(136)

due to the metric compatibility of the spacetime covariant derivative operator (i.e. since ${ }^{\left(4\right)}{\nabla }_{\nu }{g}^{\mu \nu }=\mathbf{0}$ identically), with $\mu ,\nu ,\sigma$ ranging across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$. Thus, the mass accretion rate reduces to the following surface integral over the boundary $\partial V$, with surface element $d{S}_{\mu }$ and corresponding area element $dA$:

$$\frac{dM}{dt}=-{\int }_{\partial V}{J}^{\left(\mu +1\right)}\sqrt{-det\left(g_{\mu \nu }\right)}d{S}_{\mu },\mathrm{with}d{S}_{\mu }={\gamma }_{\mu \nu }{n}^{\left(\nu +1\right)}dA,$$

(137)

which, using the definitions of the (rest) mass current vector $J^{\mu }$ (i.e. $J^{\mu }=\rho u^{\mu }$) and the (spatial) fluid velocity vector $v^{\mu }$ perceived by an observer moving in the normal direction $\mathbf{n}$:

$$v^{\mu }=\frac{u^{\left(\mu +1\right)}}{\alpha u^0}+\frac{{\beta }^{\mu }}{\alpha }=\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}\left(u^{\left(\mu +1\right)}\right)+\frac{{\beta }^{\mu }}{\alpha },$$

(138)

becomes:

$$\begin{array}{c}-{\int }_{\partial V}{J}^{\left(\mu +1\right)}\sqrt{-det\left(g_{\mu \nu }\right)}d{S}_{\mu }=-{\int }_{\partial V}\rho u^{\left(\mu +1\right)}\sqrt{-det\left(g_{\mu \nu }\right)}d{S}_{\mu }\hfill \\ \hfill =-{\int }_{\partial V}\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\left(v^{\mu }-\frac{{\beta }^{\mu }}{\alpha }\right)\right)d{S}_{\mu },\end{array}$$

(139)

with $\mu ,\nu$ in all of the above ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only. In order to evaluate this integral numerically, we exploit the radial nature of the accretion flow by sampling the integrand uniformly in the angular range $\theta \in \left[0,\frac{\pi }{2}\right]$ and then evaluating:

$$\frac{dM}{dt}=4\pi {\int }_0^{\frac{\pi }{2}}\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(\frac{\rho }{\sqrt{1-{\gamma }_{\mu \nu }{v}^{\mu }{v}^{\nu }}}\left(v^r-\frac{{\beta }^r}{\alpha }\right)\right)d\theta ,$$

(140)

using standard methods of numerical quadrature. By means of a similar argument, we follow the approach of Petrich, Shapiro, Stark and Teukolsky [50] and approximate the linear momentum accretion rate (assuming only local gravitational effects, and assuming evaluation of the integral over an asymptotically-flat boundary of spacetime) as:

$$\frac{d{P}^{\mu }}{dt}=-{\int }_{\partial V}\alpha \sqrt{det\left({\gamma }_{\mu \nu }\right)}\left(T^{\left(\mu +1\right)\left(\nu +1\right)}\right)d{S}_{\nu },$$

(141)

with $\mu ,\nu$ ranging across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only.

Next, we determine the rate of decrease in the spin value J of the black hole due to dynamical/tidal friction effects (also known as the “drag force” of the fluid) exerted by the fluid onto the underlying spacetime geometry. Within a general Bondi-Hoyle-Lyttleton accretion setup, the asymmetry of the fluid pressure distribution surrounding a spinning black hole causes an increase in fluid pressure on the side of the black hole that is counter-rotating with the fluid and a decrease on the side that is co-rotating with it, which, in turn distorts the metric in such a way as to induce a gravitational field opposing the direction of spin of the black hole. In the purely radial (Bondi) accretion setup, this represents the general relativistic analog of the viscous dissipation of angular momentum due to tidal deformation that occurs in Newtonian gravity. In order to analyze this effect quantitatively, we follow the approach of Hawking and Hartle [67], and later Hartle [68], by exploiting the known relationship between the surface area A of the horizon of a Kerr black hole and its corresponding spin value J:

$$A=8\pi M\left(M+\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}\right),$$

(142)

thus implying that the rate of decrease in the spin of the black hole $\frac{dJ}{dt}$ is related directly to the rate of decrease of the surface area of the black hole horizon $\frac{dA}{dt}$

$$\frac{dJ}{dt}=-\left(\frac{\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}}{8\pi {\left(\frac{J}{M}\right)}^2}\right)\left(\frac{dA}{dt}\right).$$

(143)

If we now choose an orthonormal tetrad such that the vector $\mathbf{l}$ is normal to the horizon of the black hole, with:

$$\frac{\partial }{\partial x^{\mu }}\left(t\right)l^{\mu }=1,$$

(144)

and the vector $\mathbf{m}$ (along with its complex conjugate $\overline{\mathbf{m}}$) is such that:

$$m^{\mu }{l}_{\mu }=0,\mathrm{and}{m}^{\mu }{\overline{m}}_{\mu }=-1,$$

(145)

then the rate of decrease of the surface area of the black hole can be rewritten as:

$$\frac{dA}{dt}=-2{\int }_S\left({ }^{\left(4\right)}{\nabla }_{\nu }{l}_{\mu }\right)m^{\mu }{\overline{m}}^{\nu }dA=-2{\int }_S\left(\frac{\partial }{\partial x^{\nu }}\left(l_{\mu }\right)-{ }^{\left(4\right)}{\Gamma }_{\nu \mu }^{\sigma }{l}_{\sigma }\right)m^{\mu }{\overline{m}}^{\nu }dA,$$

(146)

where the $\left(n-2\right)$-dimensional surface S represents the intersection of the event horizon of the black hole with a constant-time hypersurface, and the integrand represents the convergence of null geodesic generators of the event horizon in the Newman-Penrose formalism [69]. Hence, the overall rate of decrease in the spin of the black hole is given by:

$$\frac{dJ}{dt}=\left(\frac{\sqrt{M^2-{\left(\frac{J}{M}\right)}^2}}{4\pi {\left(\frac{J}{M}\right)}^2}\right)\left({\int }_S\left(\frac{\partial }{\partial x^{\nu }}\left(l_{\mu }\right)-{ }^{\left(4\right)}{\Gamma }_{\nu \mu }^{\sigma }{l}_{\sigma }\right)m^{\mu }{\overline{m}}^{\nu }dA\right).$$

(147)

Equivalently, the rate of decrease can be calculated as:

$$\frac{d\left(\frac{J}{M}\right)}{dt}=-\frac{A}{4J}{\int }_S{\left|\left({ }^{\left(4\right)}{\nabla }_{\nu }{l}_{\mu }\right)m^{\mu }{m}^{\nu }\right|}^{2}dA=-\frac{A}{4J}{\int }_S{\left|\left(\frac{\partial }{\partial x^{\nu }}\left(l_{\mu }\right)-{ }^{\left(4\right)}{\Gamma }_{\nu \mu }^{\sigma }{l}_{\sigma }\right)m^{\mu }{m}^{\nu }\right|}^{2}dA,$$

(148)

where now the integrand represents the square of the shear of null geodesic generators of the perturbed black hole event horizon, evaluated to the first-order of perturbation theory. In all of the above, $\mu ,\nu ,\sigma$ range across all spacetime coordinate indices $\left\{0,\dots ,n-1\right\}$. Within our numerical validation tests, we did not find any substantive difference in the results obtained through these two mathematical approaches, so we generally opt to use the latter due to the reduced algorithmic complexity of its implementation within Gravitas. It is instructive to compare this angular drag force experienced by the black hole in the purely radial/Bondi accretion case with the linear drag force ${\mathbf{F}}_{\infty }^D$ experienced by the black hole in the general Bondi-Hoyle-Lyttleton accretion case, as calculated analytically by Ostriker [70], yielding:

$${\mathbf{F}}_{\infty }^D=-\left(\frac{1}{2}log\left(\frac{1+{\mathcal{M}}_{\infty }}{1-{\mathcal{M}}_{\infty }}\right)-{\mathcal{M}}_{\infty }\right)\left(\frac{4\pi M^2{\rho }_{\infty }}{v_{\infty }^2}\right){\widehat{\mathbf{v}}}_{\infty },$$

(149)

in the case of subsonic flow at radial infinity (i.e. ${\mathcal{M}}_{\infty }<1$), and:

$${\mathbf{F}}_{\infty }^D=-\left(\frac{1}{2}log\left({\mathcal{M}}_{\infty }^2-1\right)+log\left(\Lambda \right)\right)\left(\frac{4\pi M^2{\rho }_{\infty }}{v_{\infty }^2}\right){\widehat{\mathbf{v}}}_{\infty },$$

(150)

in the case of supersonic flow at radial infinity (i.e. ${\mathcal{M}}_{\infty }>1$). In the above, ${\widehat{\mathbf{v}}}_{\infty }$ is the unit vector representing the relative velocity between the fluid flow (at radial infinity) and the black hole, ${\mathcal{M}}_{\infty }$ designates the Mach number at radial infinity, ${\rho }_{\infty }$ and $v_{\infty }$ denote (as usual) the fluid density and velocity at radial infinity, respectively, and $log\left(\Lambda \right)$ represents the Coulomb logarithm for particle collisions, with a typical numerical estimate (for instance derived by Chapon, Mayer and Teyssier [71] using numerical simulations of supermassive black hole binary mergers) of $log\left(\Lambda \right)=3.2$.

We find, in each of the cases simulated, that the rates of mass/energy accretion $\frac{dM}{dt}$ and linear momentum accretion $\frac{d{P}^{\mu }}{dt}$ decrease monotonically as the discretization scale of the underlying spacetime increases. This effect becomes progressively more pronounced as one increases the black hole spin value J, the dimensionless gas temperature at radial infinity ${\Theta }_{\infty }=\frac{P_{\infty }}{{\rho }_{\infty }}$, and the adiabatic exponent $\Gamma$. For instance, the rate of decrease for a black hole spinning close to extremality, with spin parameter $J=0.99M$, experiences a mass/energy accretion rate decrease that is approximately six times as rapid (as a function of discretization scale) with a stiff equation of state (i.e. $\Gamma =2$), and approximately two times as rapid with an ultra-relativistic equation of state (i.e. $\Gamma =\frac{4}{3}$), as the rate of decrease for a non-rotating black hole with the equivalent equations of state. These effects disappear in the non-relativistic limit as ${\Theta }_{\infty }\to 0$, and become divergent in the ultra-relativistic limit as $\Theta \to \infty$. We do not find any correspondingly systematic relationship between the discretization scale and the angular drag force exerted on a spinning black hole, although we do find evidence of an advective-acoustic instability within the fluid, which then propagates to become an instability in the underlying disrcete spacetime structure, that appears for certain critical values of the discretization scale, and becomes more pronounced at higher Mach numbers. This may be a purely numerical artefact (since similar such instabilities were found by Beckmann, Slyz and Devriendt [72] in simulations of supermassive black holes using the RAMSES code, and were discovered to be dependent upon numerical resolution), or may be due to some more physical “inverse energy cascade” effect caused by a truncation of fluid interactions at short length-scales. As a consequence, we treat this result as necessarily more tentative than the mass/energy and momentum accretion rate results, and believe that it warrants further and more systematic investigation.

5. Concluding Remarks

In this article, we have derived and numerically validated a new formulation of the equations of general relativistic hydrodynamics that is amenable to analysis within arbitrary discrete spacetime settings, and have implemented the resulting formalism into the Gravitas computational general relativity framework. We then proceeded to simulate radial (Bondi-type) accretion of a perfect relativistic fluid obeying an ideal gas equation of state onto both static and spinning black holes, as described by the Schwarzschild and Kerr metrics respectively, with a variety of black hole spin parameters and in a variety of different coordinate systems. Our simulations suggest that there exists a fairly robust (namely an inverse, monotonic) relationship between the mass/energy and momentum accretion rates onto the black hole and the discretization scale of the underlying spacetime. We have also found preliminary evidence of a possible advective-acoustic instability in the angular drag force exerted on the black hole by the fluid, that becomes significantly more pronounced at certain key values of the discretization scale, although further numerical experiments will be required before the extent to which this corresponds to a physically realistic effect can be confidently determined. These results provide tentative evidence that there may exist astrophysically observable effects of the underlying discreteness of spacetime arising within certain quantum gravity models that are reflected in the dynamics of the fluid region close to radially-accreting black holes, especially those whose spin values are approaching extremality, and onto whom the accretion flow is ultra-relativistic. However, in order to render the analysis both mathematically and computationally tractable, we have had to make several strong assumptions which limit the physical reasonableness and generality of our results, including the assumption of a form of the ideal gas equation of state that was shown by Taub [14] only to be physical in either the strictly non-relativistic limit (i.e. with ${\Theta }_{\infty }\to 0$ and $\Gamma =\frac{5}{3}$) or in the strictly ultra-relativistic limit (i.e. with ${\Theta }_{\infty }\to \infty$ and $\Gamma =\frac{4}{3}$), and crucially not in the relativistic case with dimensionless gas temperatures on the order of unity (i.e. ${\Theta }_{\infty }\sim 1$). Therefore, extension of the simulation results presented within this article to the case of fluid accretion involving more physically reasonable equations of state, and from the highly idealized case of purely radial/Bondi-type accretion to the more astrophysically relevant case of non-radial/Bondi-Hoyle-Lyttleton-type accretion, remains a particular priority.

Many other directions exist for future research, including the inclusion of the effects arising from certain quantum gravitational, quantum field-theoretic and/or quantum information-theoretic properties of black hole event horizons in discrete spacetimes [73,74,75] (especially within discrete black holes spinning close to extremality) into simulations of the resulting accretion dynamics, as well as the effects of certain features of the global spacetime topology that are characteristic of discrete/emergent spacetime theories [76,77]. It is also highly likely that extending these simulations into more complex astrophysical and cosmological settings involving strong relativistic field dynamics, such as the accretion of matter onto a merging binary black hole system [78,79,80], would reveal yet more intricate physics that is peculiar to the discrete spacetime setting, although such an analysis would require significant advances in the numerical algorithms employed within the Gravitas framework, certainly well beyond the capabilities of the algorithms used within this article. On the more observational side, a more complete and systematic survey of the parameter space of black hole spin values, discretization scales and equation of state parameters would be necessary in order to determine which (if any) of the effects discussed within this article might realistically be detectable within the X-ray emission spectra of black hole accretion regions in the near-term, as would the incorporation of electromagnetic effects (which would, in turn, facilitate investigations of the impact of spacetime discreteness on phenomena such as the Blandford-Znajek mechanism [81] and the formation of astrophysical jets within active galactic nuclei [82], for example) into the spacetime description. Finally, it would be particularly exciting to extend the mathematical and numerical methods developed here to other general relativistic scenarios involving the two-way interaction between perfect fluid matter and a discrete underlying spacetime in strong gravity, such as the inspiral and collision of binary neutron star systems [83,84].