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).

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: 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: 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 }$: 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: 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: 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: 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: 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: 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$: 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: 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: 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: 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]: which expands out to give, explicitly: 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 }$: 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: 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 }$: 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: where the Lie derivative term ${\mathcal{L}}_{\beta }E$ expands to give: 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: 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: 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: 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}$): 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: 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: 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: 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: where $\alpha u^0$ represents the Lorentz factor of the fluid: 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: 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: 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: yields: 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: 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}$: 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: 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: 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$: 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: 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: 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: 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: 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: 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: 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: which have algebraic multiplicity n, and those corresponding to the acoustic wave speeds, i.e: 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: where ${ }^{\left(4\right)}{G}_{\mu \nu }$ is the spacetime Einstein tensor: ${ }^{\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$ 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: 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: which then expand out to give: 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 }$: 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: Finally, the elliptic constraint equations resulting from the contracted Bianchi identities may be decomposed into a timelike projection, yielding the Hamiltonian constraint equation: 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: which then expand out to give: 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.

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: for the energy conservation equation; for the momentum conservation equations; and: 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: with the corresponding weak form: the momentum conservation equations reduce to: with their corresponding weak forms being: and the baryon number continuity equation reduces to: with its corresponding weak form being: and where $\mu ,\nu ,\rho ,\sigma$ in all of the above range, again, across spatial coordinate indices $\left\{0,\dots ,n-2\right\}$ only.