The first article in this series [1] contained a more-or-less complete and self-contained introduction to the underlying philosophy of the Gravitas computational framework, a summary of its core design principles, and an incomplete but high-level survey of the field of computational and numerical relativity as a whole, along with a detailed explanation of how Gravitas compared to the many other existing open source tools that are presently available. In the interests of brevity, we do not seek to repeat those introductory remarks here, and we instead encourage the reader to consult the previous article before proceeding with this one. The previous article also provided a comprehensive overview of Gravitas’s symbolic and analytical capabilities, including its tensor calculus and differential geometry functionality, its methodology for handling general Riemannian/pseudo-Riemannian/Lorentzian metrics and coordinate systems, and its algorithms for solving the Einstein field equations both analytically and numerically, including particular support for electromagnetic fields and the Einstein-Maxwell equations, as well as support for arbitrary non-gravitational fields via the stress-energy tensor. The present article will place particular emphasis upon the seamless integration between the aforementioned symbolic and analytical functionality and Gravitas’s powerful numerical subsystem. One of the primary design features of Gravitas is its synthesis of the capabilities of an abstract tensor calculus package that is comparable in many respects to xAct (a Wolfram Language framework developed by Martín-García et al., facilitating component calculus through xCoba, manipulation of scalar invariants of the Riemann tensor through Invar [2] [3] and canonicalization of tensor indices with respect to slot/permutation symmetries through xPerm [4]) or EinsteinPy [5] (a Python package facilitating elementary tensor manipulations, as well as the visualization and computation of geodesics within a small number of in-built spacetime geometries), with the capabilities of a full-scale numerical relativity framework such as Cactus [6]/McLachlan [7,8] (a component of the popular Einstein toolkit [9], facilitating the numerical solution of the Baumgarte-Shapiro-Shibata-Nakamura/BSSN formulation of the Einstein field equations [10] [11,12] using finite-difference approximations with up to eighth-order accuracy and block-structured adaptive mesh refinement/AMR), GRChombo [13,14] (a similar block-structured AMR code for solving both the BSSN and conformally-covariant Z4/CCZ4 formulations [15,16] of the Einstein field equations using fourth-order finite-difference methods) or SpEC [17] [18] (a code for solving the generalized harmonic formulation of the Einstein field equations using pseudospectral methods).

The principal advantage of this integration is the resulting ability of the user to configure initial conditions (i.e. initial spacelike hypersurface/Cauchy surface configurations), to define the desired evolution equations, to set up corresponding matter field configurations with appropriate coupling to the spacetime, to enforce appropriate gauge and coordinate constraints on the hyperbolic evolution, to run the resulting (automatically-optimized and automatically-parallelized) numerical simulations, and to analyze and visualize the raw outputs of those simulations (including the extraction of gravitational wave signals), all from within a single notebook environment. In particular, this eliminates the need to use separate computer algebra-based frameworks such as Kranc [19] in the Wolfram Language or SENR/NRPy+ [20] in Python for representing evolution equations and gauge conditions in an abstract tensorial form and then automatically synthesizing optimized low-level C or Fortran code from them, and indeed the need to use a separate tool such as SimulationTools in the Wolfram Language or kuibit [21] in Python for performing standard post-processing tasks and the extraction or analysis of key simulation data. Gravitas’s use of next-generation adaptive refinement algorithms based upon totally-unstructured hypergraph rewriting/Wolfram model evolution [22,23] [24] facilitates its handling of arbitrary (potentially overlapping) curvilinear coordinate systems and highly non-trivial spacetime topologies with ease, since one is not constrained to use a rigid Cartesian or spheroidal grid structure for the computational domain, and the hypergraph’s topology can instead adapt organically to both the topology and the dynamics of the problem being solved. This hypergraph rewriting formalism [25,26,27] (including its efficient hypergraph canonicalization system [28]), and the resulting numerical algorithms that it enables, have previously been used to investigate a variety of physical problems in discrete spacetimes, including binary black hole mergers [29], entanglement entropies in curved spacetime quantum field theories [30], singularity theorems in idealized gravitational collapse scenarios [31], amongst others.

In Section 2, we begin by introducing the abstract formalism of the ADM (or $3+1$) metric decomposition of a Riemannian/pseudo-Riemannian/Lorentzian manifold into a sequence of immersed hypersurfaces, and demonstrate how Gravitas is able to use the ADMDecomposition object (the most fundamental object in Gravitas’s numerical subsystem, analogous to MetricTensor in the symbolic/analytial subsystem) to represent such decompositions and their corresponding geometrical properties using only the specification of the initial Cauchy surface of the manifold and a suitable gauge choice. In particular, we show an excerpt of Gravitas’s in-built library of Cauchy initial data (including initial conditions for both black hole and cosmological spacetimes), as well as functionality for calculating normal vectors, “time vectors”, Gauss equations, Codazzi-Mainardi equations, ADM Hamiltonian constraints, ADM momentum constraints and (purely geometrical) ADM evolution equations from the underlying ADMDecomposition object. We also highlight functions such as ExtrinsicCurvatureTensor for computing the second fundamental form/extrinsic curvature of immersed hypersurfaces within the decomposition, and CottonTensor for determining conformal-flatness of immersed three-dimensional submanifolds. In Section 3, we proceed to demonstrate Gravitas’s library of supported ADM gauge conditions, including the geodesic, maximal, harmonic and $1+log$ slicing conditions on the lapse function, and the normal, harmonic, minimal distortion and pseudo-minimal distortion coordinate conditions on the shift vector. We also show how the (vacuum) ADM evolution equations and (vacuum) ADM Hamiltonian and momentum constraint equations can be solved, both analytically and numerically, using the SolveVacuumADMEquations function, with such solutions being represented symbolically as VacuumADMSolution objects. In Section 4, we illustrate how one can introduce a generic matter field into such a solution (and evolve it either numerically or analytically) using the ADMStressEnergyDecomposition function, which produces an ADM/$3+1$ decomposition of an arbitrary StressEnergyTensor object, thus allowing one to compute energy, momentum and stress projections of the tensor, as well as timelike (energy conservation) and spacelike (momentum conservation) projections of the stress-energy continuity equations, within any ADM decomposition and with any choice of gauge. We show how this can be combined with the aforementioned functionality in order to solve the full ADM evolution equations and full ADM Hamiltonian and momentum constraint equations, again both analytically and numerically, using the SolveADMEquations function, with such solutions being represented symbolically as ADMSolution objects. Finally, in Section 5, we demonstrate how the DiscreteHypersurfaceDecomposition function can be used to expose and visualize some of the internal details of how Gravitas obtains its numerical solutions, namely by means of a fourth-order Runge-Kutta time integration algorithm defined over hypergraphs of arbitrary topology, with a generalized local adaptive refinement algorithm that automatically refines and coarsens the hypergraph (using hypergraph rewriting methods) based upon the dynamics of the problem being solved. We show how all of these methods can be combined in order to simulate a non-trivial numerical relativistic problem in vacuum spacetime: the head-on collision of a binary black hole system, including the extraction of the resulting outgoing gravitational radiation. We conclude in Section 6 with some directions for future research and development.

Note that, as with the previous article [1] in this series, much of the Gravitas functionality that is demonstrated and discussed here is fully-documented and available via the Wolfram Function Repository, including ADMDecomposition, SolveVacuumADMEquations and DiscreteHypersurfaceDecomposition. However, there are many more functions that have not yet been documented in this way, or which have received further development beyond their documented versions; as always, an up-to-date version of the Gravitas framework may be obtained from its GitHub repository (currently approximately 30,000 lines of Wolfram Language code, including experimental and research functionality). Note moreover that this article will employ the same notational and terminological conventions as the previous one [1]; for instance, we use geometric units with $c=G=\hslash =1$, we assume a metric signature of $\left(-,+,+,+\right)$ in all relevant cases, we adopt the Einstein summation convention of summing over all repeated tensor indices throughout, and Gravitas functions use the keywords “Reduced” and “Symbolic” within certain property names in order to signify that Gravitas should apply all known algebraic equivalences and return a given tensor expression in canonical form, or that Gravitas should leave all partial derivatives as purely symbolic/unevaluated operators, respectively.

At the core of the Gravitas framework’s numerical subsystem is the ADMDecomposition object, which plays a role analogous to the role played by the MetricTensor object within Gravitas’s symbolic subsystem [1], and which allows one to represent a decomposition or “foliation” of a given differentiable Riemannian or pseudo-Riemannian manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ into an ordered sequence of differentiable submanifolds $\left(\mathcal{N},g_{\mathcal{N}}\right)$, where $\mathcal{N}\subset \mathcal{M}$, each of codimension-1 [32,33]. In the context of general relativity, the “ambient” manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ is typically a pseudo-Riemannian/Lorentzian spacetime, and the submanifolds $\left(\mathcal{N},g_{\mathcal{N}}\right)$ are typically Riemannian spacelike hypersurfaces, ordered by a gauge-dependent “time” parameter [34]. Mathematically, if $i:\mathcal{N}\to \mathcal{M}$ denotes the embedding of the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$ into the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$, then the tangent bundle $\underset{\mathbf{x}\in \mathcal{N}}{⨆}{T}_{\mathbf{x}}^{\top }\mathcal{N}$ on the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$ is included in the tangent bundle $\underset{\mathbf{x}\in \mathcal{M}}{⨆}{T}_{\mathbf{x}}^{\top }\mathcal{M}$ on the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ by means of the pushforward construction, whose cokernel/quotient is then the normal bundle $\underset{\mathbf{x}\in \mathcal{N}}{⨆}{T}_{\mathbf{x}}^{\perp }\mathcal{N}$ on the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$, i.e. one has the following short exact sequence [35,36]:

In the above, $T_{\mathbf{x}}^{\top }\mathcal{N}$ and $T_{\mathbf{x}}^{\top }\mathcal{M}$ denote the tangent spaces of the manifolds $\left(\mathcal{N},g_{\mathcal{N}}\right)$ and $\left(\mathcal{M},g_{\mathcal{M}}\right)$ at points $\mathbf{x}\in \mathcal{N}$ and $\mathbf{x}\in \mathcal{M}$, respectively; ${\left.\left(\underset{\mathbf{x}\in \mathcal{M}}{⨆}{T}_{\mathbf{x}}^{\top }\mathcal{M}\right)\right|}_{i\left(\mathcal{N}\right)}$ designates the restriction of the tangent bundle $\underset{\mathbf{x}\in \mathcal{M}}{⨆}{T}_{\mathbf{x}}^{\top }\mathcal{M}$ on the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ to the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$, i.e., more formally, it designates the pullback $i^{*}\left(\underset{\mathbf{x}\in \mathcal{M}}{⨆}{T}_{\mathbf{x}}^{\top }\mathcal{M}\right)$ of the tangent bundle $\underset{\mathbf{x}\in \mathcal{M}}{⨆}{T}_{\mathbf{x}}^{\top }\mathcal{M}$ on the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ to a corresponding vector bundle on the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$; $T_{\mathbf{x}}^{\perp }\mathcal{N}$ denotes the normal space of the manifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$ at the point $\mathbf{x}\in \mathcal{N}$, i.e. the set of all vectors in the tangent space $T_{\mathbf{x}}^{\top }\mathcal{M}$ of the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ that are normal to every vector in the tangent space $T_{\mathbf{x}}^{\top }\mathcal{N}$ of the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$ [37]: such that the normal bundle $\underset{\mathbf{x}\in \mathcal{N}}{⨆}{T}_{\mathbf{x}}^{\perp }\mathcal{N}$ is simply given by the quotient bundle: and 0 designates the trivial (rank-0) vector bundle.

The metric tensor $g_{\mathcal{M}}$ on the ambient manifold $\mathcal{M}$ has the effect of splitting this exact sequence, thus guaranteeing that the restriction of the tangent bundle on the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ to the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$ can be written as the following direct sum of tangent and normal bundles on the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$:

This splitting allows one, in turn, to write the Levi-Civita connection ${\nabla }^{\left(\mathcal{M}\right)}$ on the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ as a sum of tangential and normal components, i.e. one has: where $\mathbf{X}\in T_{\mathbf{x}}\mathcal{M}$ denotes an arbitrary tangent vector at the point $\mathbf{x}\in \mathcal{M}$, $\mathbf{Y}\in \mathsf{\Gamma }\left(\underset{\mathbf{x}\in \mathcal{M}}{⨆}{T}_{\mathbf{x}}\mathcal{M}\right)$ denotes an arbitrary vector field on the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$, $\mathsf{\Gamma }\left(\underset{\mathbf{x}\in \mathcal{M}}{⨆}{T}_{\mathbf{x}}\mathcal{M}\right)$ denotes the space of smooth sections of the tangent bundle $\underset{\mathbf{x}\in \mathcal{M}}{⨆}{T}_{\mathbf{x}}\mathcal{M}$ (i.e. the space of vector fields definable on the ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$), the tangential component $\top \left({\nabla }_{\mathbf{X}}^{\left(\mathcal{M}\right)}\mathbf{Y}\right)$ denotes the orthogonal projection of the covariant derivative ${\nabla }_{\mathbf{X}}^{\left(\mathcal{M}\right)}\mathbf{Y}$ onto the tangent bundle $\underset{\mathbf{x}\in \mathcal{N}}{⨆}{T}_{\mathbf{x}}^{\top }\mathcal{N}$, and the normal component $\perp \left({\nabla }_{\mathbf{X}}^{\left(\mathcal{M}\right)}\mathbf{Y}\right)$ denotes the orthogonal projection of the covariant derivative ${\nabla }_{\mathbf{X}}^{\left(\mathcal{M}\right)}\mathbf{Y}$ onto the normal bundle $\underset{\mathbf{x}\in \mathcal{N}}{⨆}{T}_{\mathbf{x}}^{\perp }\mathcal{N}$. The induced Levi-Civita connection ${\nabla }^{\left(\mathcal{N}\right)}$ on the submanifold $\left(\mathcal{N},g_{\mathcal{N}}\right)$ is then given by the orthogonal projection of the ambient Levi-Civita connection ${\nabla }^{\left(\mathcal{M}\right)}$ onto the tangent bundle: and the second fundamental form/extrinsic curvature tensor K (a symmetric, vector-valued differential form taking values in the normal bundle $\underset{\mathbf{x}\in \mathcal{N}}{⨆}{T}_{\mathbf{x}}^{\perp }\mathcal{N}$) is given by the orthogonal projection of the ambient Levi-Civita connection ${\nabla }^{\left(\mathcal{M}\right)}$ onto the normal bundle [36]: i.e. one has, overall:

Returning now to the specific case of general relativity, in which our ambient manifold $\left(\mathcal{M},g_{\mathcal{M}}\right)$ is a Lorentzian spacetime, we proceed to define a universal “time function” $t:\mathcal{M}\to \mathbb{R}$ mapping each point/event to its corresponding coordinate time value; the submanifolds $\left(N,g_{\mathcal{N}}\right)$ are then given by the level surfaces of this universal time function, i.e:

We thus obtain a time-ordered sequence of purely spacelike hypersurfaces (each of codimension-1, with purely Riemannian signature) $\left({\mathcal{N}}_{t_0},g_{{\mathcal{N}}_{t_0}}\right)$ which covers the entire ambient spacetime $\left(\mathcal{M},g_{\mathcal{M}}\right)$:

If our ambient spacetime $\left(\mathcal{M},g_{\mathcal{M}}\right)$ is, moreover, globally hyperbolic, then these spacelike hypersurfaces are also strictly non-intersecting, i.e:

Henceforth, we shall denote the overall metric tensor on the ambient spacetime manifold $\mathcal{M}$ (assumed to be of dimension n) by g, and the “induced” metric tensor on the spacelike hypersurfaces ${\mathcal{N}}_{t_0}$ (assumed to be of dimension $n-1$) by $\gamma$. Concretely, ADMDecomposition represents such spacetime “foliations” by means of the ADM formalism first proposed by Arnowitt, Deser and Misner [38,39] (although the particular mathematical formulation of the $3+1$ decomposition that Gravitas uses is due to York [40]), wherein one defines a metric tensor $\gamma$ on some “initial” spacelike hypersurface $\left({\mathcal{N}}_{t_0},\gamma \right)$, as well as the Lagrange multipliers $\alpha$ (a scalar field known as the lapse function) and ${\beta }^{\mu }$ (an $\left(n-1\right)$-dimensional vector field known as the shift vector, with $\mu \in \left\{0,\cdots ,n-2\right\}$, yet in which every such ${\beta }^{\mu }$ depends upon all n spacetime coordinates), which play the role of the gauge variables and consequently define how neighboring hypersurfaces are related to one another geometrically. More precisely, the lapse function $\alpha$ determines the proper time distance $d\tau$ between corresponding points on the spacelike hypersurfaces at coordinate time values $t=t_0$ and $t=t_0+dt$, as measured along the unit normal direction $\mathbf{n}$ to the hypersurface at $t=t_0$ (i.e. the direction followed by normal observers): whereas the shift vector (field) ${\beta }^{\mu }$ determines how the spatial coordinates $x^{\mu }\left(t_0\right)$ get relabeled as one moves from the spacelike hypersurface at coordinate time value $t=t_0$ to the spacelike hypersurface at coordinate time value $t=t_0+dt$, i.e:

The overall line element/first fundamental form $d{s}^2$ for the ambient (spacetime) MetricTensor object $g_{\mu \nu }$ is therefore given in terms of the components of the induced (spatial) MetricTensor object ${\gamma }_{\mu \nu }$, the lapse function $\alpha$ and the components of the shift vector (field) ${\beta }^{\mu }$ as: with $\mu ,\nu ,\sigma$ ranging across all $\left\{0,\cdots ,n-2\right\}$ (i.e across spatial coordinate indices only). Note that the spatial line element/first fundamental form $d{l}^2$ for the induced (spatial) MetricTensor object ${\gamma }_{\mu \nu }$ constituting the initial spacelike hypersurface of the Schwarzschild metric [41], as developed independently by Droste [42] (representing, for instance, an uncharged, non-rotating black hole with mass M in Schwarzschild or spherical polar coordinates $\left(t,r,\theta ,\varphi \right)$) is given by: whereas, for the initial spacelike hypersurface of the Kerr metric [43] (representing, for instance, an uncharged, spinning black hole with mass M and angular momentum J in Boyer-Lindquist [44] or oblate spheroidal coordinates $\left(t,r,\theta ,\varphi \right)$), the spatial line element is given instead by: with $\mu ,\nu$ ranging across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only) in both cases. The corresponding ADMDecomposition objects, including representations of both the induced/spatial and ambient/spacetime MetricTensor objects, for the Schwarzschild and Kerr metrics are shown in Figure 1, assuming the most general choice of gauge consisting of the lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and the shift vector (field) $\left({\beta }^1\left(t,r,\theta ,\varphi \right),{\beta }^2\left(t,r,\theta ,\varphi \right),{\beta }^3\left(t,r,\theta ,\varphi \right)\right)$. By default, the most general forms of the ADM gauge variables are always assumed, and appropriate formal symbols are chosen for the various parameters of the decomposition (such as mass M and angular momentum J), as well as for the coordinates (such as time coordinate t and radial coordinate r), although these defaults may always be overridden by passing additional arguments to ADMDecomposition, as illustrated in Figure 2 for the case of modified shift vectors $\left(\beta \left(t,r,\theta ,\varphi \right),0,0\right)$ (for Schwarzschild) and $\left(0,0,\beta \left(t,r,\theta ,\varphi \right)\right)$ (for Kerr).

Note that several other Cauchy surface geometries are built into ADMDecomposition’s small (though continually growing) library of initial spacelike hypersurface configurations, including Cauchy initial data for the Reissner-Nordström metric [45,46], as developed independently by Weyl [47] and Jeffery [48] (representing, for instance, a charged, non-rotating black hole with mass M and electric charge Q in Schwarzschild or spherical polar coordinates $\left(t,r,\theta ,\varphi \right))$ and the Kerr-Newman metric [49,50] (representing, for instance, a charged, spinning black hole with mass M, angular momentum J and electric charge Q in Boyer-Lindquist or oblate spheroidal coordinates $\left(t,r,\theta ,\varphi \right)$), with spatial line elements given by: and: respectively, as shown in Figure 3. Cosmological Cauchy initial data is also supported, such as for the Friedmann-Lemaître-Robertson-Walker or FLRW metric [51,52,53,54] (representing, for instance, a perfectly homogeneous and isotropic universe with scale factor $a\left(t\right)$ and curvature k in spherical polar coordinates $\left(t,r,\theta ,\varphi \right)$), with spatial line element given by: as shown in Figure 4, both for the default choice of gauge with general lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and general shift vector (field) $\left({\beta }^1\left(t,r,\theta ,\varphi \right),{\beta }^2\left(t,r,\theta ,\varphi \right),{\beta }^3\left(t,r,\theta ,\varphi \right)\right)$, and for the trivial choice of gauge with lapse function identically equal to 1 (i.e. $\alpha \left(t,r,\theta ,\varphi \right)=1$) and shift vector (field) identically equal to zero (i.e. $\beta \left(t,r,\theta ,\varphi \right)=\left(0,0,0\right)$). In the above, as before, $\mu$ and $\nu$ represent spatial indices, and therefore range across $\left\{0,\cdots ,n-2\right\}$ only.

In all that follows, we will use a bracketed “3” to indicate that a given object is restricted to spacelike hypersurfaces, with a bracketed “4” indicating instead that a given object is defined over the entire ambient spacetime (for instance, ${}^{\left(3\right)}{\nabla }_{\mu }$ denotes the covariant derivative operator on spacelike hypersurfaces, whose corresponding Christoffel symbols ${}^{\left(3\right)}{\mathsf{\Gamma }}_{\mu \nu }^{\rho }$ are defined in terms of derivatives of the induced/spatial metric tensor ${\gamma }_{\mu \nu }$, while ${}^{\left(4\right)}{\nabla }_{\mu }$ denotes the covariant derivative operator on the ambient spacetime, whose corresponding Christoffel symbols ${}^{\left(4\right)}{\mathsf{\Gamma }}_{\mu \nu }^{\rho }$ are defined in terms of derivatives of the overall spacetime metric tensor $g_{\mu \nu }$); note that we use these particular numbers purely for illustrative purposes (since the formalism is often referred to as “$3+1$ formalism”) - the objects themselves are defined within Gravitas for arbitrary n-dimensional spacetimes, which are then foliated into $\left(n-1\right)$-dimensional spacelike hypersurfaces. The future-pointing, timelike unit vector $\mathbf{n}$ that is normal to each spacelike hypersurface can now be computed as the spacetime contravariant derivative ${}^{\left(4\right)}{\nabla }^{\mu }$ of the distinguished time coordinate t [34]: whereas the “time vector” $\mathbf{t}$ determines how corresponding points on neighboring hypersurfaces are related (and so, in particular, will differ from the normal vector $\mathbf{n}$ if the corresponding observer does not move in a normal direction to the hypersurfaces, which in turn corresponds to the case of a non-vanishing shift vector $\beta \ne \mathbf{0}$):

In both of the above, $\mu ,\sigma$ range across all $\left\{0,\cdots ,n-1\right\}$ (i.e. across all spacetime coordinate indices). Figure 5 and Figure 6 show the components of the future-pointing unit normal vector $n^{\mu }$ and of the future-pointing unit “time vector” $t^{\mu }$, computed directly from the respective ADMDecomposition objects for the Schwarzschild and Kerr metrics, assuming a fully generic choice of gauge with lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and shift vector $\left({\beta }^1\left(t,r,\theta ,\varphi \right),{\beta }^2\left(t,r,\theta ,\varphi \right),{\beta }^3\left(t,r,\theta ,\varphi \right)\right)$ in both cases. Recall from the discussion above that the second fundamental form/extrinsic curvature tensor K is simply the symmetric, vector-valued differential form given by the orthogonal projection ⊥ of the ambient/spacetime Levi-Civita connection ${}^{\left(4\right)}\nabla$ onto the normal bundle; we can express this in component form by taking the Lie derivative $\mathcal{L}$ of the spatial metric tensor ${\gamma }_{\mu \nu }$ in the normal direction $\mathbf{n}$ [55]: i.e., in expanded form: where $\mu ,\nu ,\sigma$ range across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only), as shown in Figure 7 for the case of ExtrinsicCurvatureTensor objects computed directly from the ADMDecomposition objects for the Schwarzschild and Kerr metrics, assuming again a fully generic choice of gauge with lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and shift vector $\left({\beta }^1\left(t,r,\theta ,\varphi \right),{\beta }^2\left(t,r,\theta ,\varphi \right),{\beta }^3\left(t,r,\theta ,\varphi \right)\right)$. Note that both the shift vector ${\beta }^{\mu }$ and the extrinsic curvature tensor $K_{\mu \nu }$ are purely spatial objects, and so their indices are raised and lowered using the spatial metric tensor ${\gamma }_{\mu \nu }$ rather than full spacetime metric tensor $g_{\mu \nu }$, i.e: for the case of the rank-1 shift vector/covector (field) $\beta$, and: for the case of the rank-2 extrinsic curvature tensor (field) K. Their covariant derivatives ${}^{\left(3\right)}{\nabla }_{\mu }$ (i.e. their derivatives along tangent vectors) are therefore represented in terms of the coefficients of the induced/spatial Levi-Civita connection ${}^{\left(3\right)}\nabla$, i.e. the induced/spatial Christoffel symbols ${}^{\left(3\right)}{\mathsf{\Gamma }}_{\mu \nu }^{\rho }$, which are themselves represented in terms of partial derivatives of the spatial metric tensor ${\gamma }_{\mu \nu }$: namely: for the case of the rank-1 shift vector/covector (field) $\beta$, and: for the case of the rank-2 extrinsic curvature tensor (field) K. This is demonstrated in Figure 8, in which all covariant derivatives of the ExtrinsicCurvatureTensor object in lowered-index/covariant form (i.e. $K_{\mu \nu }$, the default case), as well as all covariant derivatives of the ExtrinsicCurvatureTensor object in raised-index/contravariant form (i.e. $K^{\mu \nu }$) are computed for the ADMDecomposition objects for the Schwarzschild and Kerr metrics, respectively, assuming generic gauge choice.

As stated previously, the orthogonal projection ⊤ of the ambient/spacetime Levi-Civita connection ${}^{\left(4\right)}\nabla$ onto the tangent bundle yields the induced/spatial Levi-Civita connection ${}^{\left(3\right)}\nabla$, while the orthogonal projection ⊥ of the ambient/spacetime Levi-Civita connection ${}^{\left(4\right)}\nabla$ onto the normal bundle yields the symmetric (vector-valued) second fundamental form/extrinsic curvature tensor K. The statement that this is necessarily the case is often referred to as the Gauss relation [56], which asserts the following general relationship between (projections of) the ambient/spacetime Riemann tensor ${}^{\left(4\right)}R$ and the induced/spatial Riemann tensor ${}^{\left(3\right)}R$ (assuming that, as before, that $\left(\mathcal{M},g_{\mathcal{M}}\right)$ denotes the ambient/spacetime manifold and $\left(\mathcal{N},g_{\mathcal{N}}\right)$ denotes the induced/spatial submanifold): or, in explicit (component-based) form: where ${\perp }_{\mu }^{\nu }$ designates the orthogonal projector, i.e. the projection operator in the normal direction $\mathbf{n}$: with ${\delta }_{\mu }^{\nu }$ being the identity tensor (i.e. the Kronecker delta function). In the above, $\rho ,\sigma ,\mu ,\nu$ range across all $\left\{0,\cdots ,n-1\right\}$ (i.e. across all spacetime coordinate indices), while $\alpha ,\beta ,\gamma ,\delta$ range across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only), with ${}^{\left(4\right)}{R}_{\sigma \mu \nu }^{\rho }$ being the components of the ambient/spacetime Riemann tensor: where $\lambda$ here ranges across all spacetime coordinate indices $\left\{0,\cdots ,n-1\right\}$, ${}^{\left(3\right)}{R}_{\beta \gamma \delta }^{\alpha }$ being the components of the induced/spatial Riemann tensor: where $\lambda$ here ranges across spatial coordinate indices $\left\{0,\cdots ,n-2\right\}$ only, and with ${}^{\left(4\right)}{\mathsf{\Gamma }}_{\mu \nu }^{\rho }$ being the ambient/spacetime Christoffel symbols, represented in terms of partial derivatives of the overall spacetime metric tensor $g_{\mu \nu }$: The Gauss equations necessarily hold identically for any ADMDecomposition object computed directly from an (initial) spatial MetricTensor object in the manner described previously, as illustrated in Figure 9 for the cases of the Schwarzschild metric (representing, for instance, an uncharged, non-rotating black hole with mass M in Schwarzschild or spherical polar coordinates $\left(t,r,\theta ,\varphi \right)$) and the Kerr metric (representing, for instance, an uncharged, spinning black hole with mass M and angular momentum J in Boyer-Lindquist or oblate spheroidal coordinates $\left(t,r,\theta ,\varphi \right)$), assuming a restricted choice of gauge consisting of the lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and the modified shift vectors $\left(\beta \left(t,r,\theta ,\varphi \right),0,0\right)$ (for Schwarzschild) and $\left(0,0,\beta \left(t,r,\theta ,\varphi \right)\right)$ (for Kerr). On the other hand, the Codazzi-Mainardi relation [36,56] asserts the following general relationship between (orthogonal projections of) the ambient/spacetime Ricci tensor and the covariant derivatives of the second fundamental form/extrinsic curvature tensor: or, in explicit (component-based) form: i.e., in expanded form: with $\rho ,\sigma$ ranging across all $\left\{0,\cdots ,n-1\right\}$ (i.e. across all spacetime coordinate indices), while $\mu ,\nu ,\lambda$ range across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only). In the above, ${}^{\left(4\right)}{R}_{\mu \nu }$ designates the ambient/spacetime Ricci tensor, obtained through appropriate contraction of the ambient/spacetime Riemann tensor ${}^{\left(4\right)}{R}_{\sigma \mu \nu }^{\rho }$ (i.e. ${}^{\left(4\right)}{R}_{\mu \nu }={}^{\left(4\right)}{R}_{\mu \sigma \nu }^{\sigma }$), and K designates the trace of the extrinsic curvature tensor $K_{\mu \nu }$ (i.e. $K=K_{\mu }^{\mu }={\gamma }^{\mu \nu }{K}_{\mu \nu }$). Just as for the Gauss equations, the Codazzi-Mainardi equations also necessarily hold identically for any ADMDecomposition object computed directly from an (initial) spatial MetricTensor object in this way, as demonstrated in Figure 10 for the cases of the Schwarzschild and Kerr metrics with a restricted choice of gauge consisting, as before, of the lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and the modified shift vectors $\left(\beta \left(t,r,\theta ,\varphi \right),0,0\right)$ (for Schwarzschild) and $\left(0,0,\beta \left(t,r,\theta ,\varphi \right)\right)$ (for Kerr). In some sense, the Gauss and Codazzi-Mainardi equations therefore represent a set of geometrical consistency conditions that must necessarily be satisfied in order for the hypersurfaces in the foliation to “plumb together” in some appropriate way.

In much the same vein as the Gauss and Codazzi-Mainardi equations above, it is also possible to relate the components of the ambient/spacetime Ricci tensor ${}^{\left(4\right)}{R}_{\mu \nu }$ and the induced/spatial Ricci tensor ${}^{\left(3\right)}{R}_{\mu \nu }$; by making the time derivative of the extrinsic curvature tensor $\frac{\partial }{\partial t}\left(K_{\nu }^{\mu }\right)$ the subject of these equations, one can then proceed to formulate them as a system of time evolution equations for the components of the extrinsic curvature tensor $K_{\mu \nu }$ (regarded here as the conjugate momenta of the components of the spatial metric tensor ${\gamma }_{\mu \nu }$, which are in turn regarded here as the dynamical variables of the theory): i.e., in expanded form: with $\mu ,\nu ,\rho ,\sigma$ ranging across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only). Since these are purely geometrical identities relating components of ${}^{\left(4\right)}{R}_{\mu \nu }$ to components of ${}^{\left(3\right)}{R}_{\mu \nu }$, these evolution equations necessarily hold identically for any ADMDecomposition object computed from an (initial) spatial MetricTensor object in the usual way, as shown in Figure 11 for the cases of the Schwarzschild and Kerr metrics with the same restricted choice of gauge as before, consisting of the lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and the modified shift vectors $\left(\beta \left(t,r,\theta ,\varphi \right),0,0\right)$ (for Schwarzschild) and $\left(0,0,\beta \left(t,r,\theta ,\varphi \right)\right)$ (for Kerr). However, as we shall see shortly, once the Einstein field equations are used to replace the ambient/spacetime Ricci tensor terms ${}^{\left(4\right)}{R}_{\mu \nu }$ with stress-energy tensor terms $T_{\mu \nu }$, these purely geometrical relations instead become bona fide evolution equations, describing how the components of the extrinsic curvature tensor $K_{\mu \nu }$ on spacelike hypersurfaces change as one moves forwards (or backwards) in time. By taking appropriate contractions of the Gauss equations, we also obtain a timelike projection of the contracted Bianchi identities, known as the Hamiltonian constraint $\mathcal{H}$: or, more succinctly: with $\mu ,\nu$ ranging across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only), where we have used the definition of the ambient/spacetime Einstein tensor ${}^{\left(4\right)}{G}_{\mu \nu }$: with $\mu ,\nu$ here ranging across all $\left\{0,\cdots ,n-1\right\}$ (i.e. across all spacetime coordinate indices), and where ${}^{\left(3\right)}R$ and ${}^{\left(4\right)}R$ denote the induced/spatial and ambient/spacetime Ricci scalars: respectively. Likewise, one can obtain the corresponding spacelike projections of the contracted Bianchi identities by taking appropriate contractions of the Codazzi-Mainardi equations, thus yielding the so-called momentum constraints ${\mathcal{M}}_{\mu }$ (represented here in covector form): i.e., in expanded form: which, upon simplification again using the ambient/spacetime Einstein tensor ${}^{\left(4\right)}{G}_{\mu \nu }$, become: with $\mu ,\nu ,\sigma$ here ranging across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only). Just as with the evolution equations above, these forms of the Hamiltonian $\mathcal{H}$ and momentum ${\mathcal{M}}_{\mu }$ constraints are purely geometrical identities, and therefore vanish identically for any ADMDecomposition object obtained from an (initial) spatial MetricTensor object in the same way, as illustrated in Figure 12 and Figure 13, again for the cases of the Schwarzschild and Kerr metrics with modified shift vectors $\left(\beta \left(t,r,\theta ,\varphi \right),0,0\right)$ and $\left(0,0,\beta \left(t,r,\theta ,\varphi \right)\right)$, respectively. The equations asserting that the Hamiltonian and momentum constraints vanish identically, i.e. that $\mathcal{H}=0$ and ${\mathcal{M}}_{\mu }=0$, are consequently referred to as the Hamiltonian and momentum constraint equations, and they can be computed and verified directly within Gravitas, as shown in Figure 14 and Figure 15. Much like the evolution equations, once the Einstein field equations are applied and used to transform the ambient/spacetime Einstein tensor terms ${}^{\left(4\right)}{G}_{\mu \nu }$ into stress-energy tensor terms $T_{\mu \nu }$, the Hamiltonian and momentum constraint equations will cease to be purely geometrical identities and will instead become bona fide physical constraints on the evolution and choice of gauge. Indeed, in (for instance) the $3+1$-dimensional case, the ten independent Einstein field equations may be projected in the space-space directions (thus yielding the six evolution equations), the time-time direction (thus yielding the one Hamiltonian constraint equation) and the time-space/space-time directions (thus yielding the three momentum constraint equations), with the latter two corresponding to the four redundant degrees of freedom arising from the contracted Bianchi identities.

In all dimensions $n\ge 4$, the vanishing of all components of the unique, trace-free, rank-4 Weyl tensor $C_{\sigma \mu \nu }^{\rho }$ [57], obtained by subtracting out all trace components (i.e. all Ricci tensor components $R_{\mu \nu }$) from the full Riemann tensor $R_{\sigma \mu \nu }^{\rho }$ via the Ricci decomposition theorem [58]: is both a necessary and sufficient condition for the underlying manifold $\left(\mathcal{M},g\right)$ to be conformally-flat (by the Weyl-Schouten theorem [59]), since the Weyl tensor $C_{\sigma \mu \nu }^{\rho }$ is invariant under all conformal transformations of the form $\widetilde{g_{\mu \nu }}={\mathsf{\Omega }}^2{g}_{\mu \nu }$ with real conformal factor $\mathsf{\Omega }\in \mathbb{R}$; the Weyl tensor $C_{\sigma \mu \nu }^{\rho }$ is also significant from a physical standpoint, since it governs the propagation of gravitational radiation within vacuum spacetimes. However, the fact that the Weyl tensor $C_{\sigma \mu \nu }^{\rho }$ vanishes identically in dimension $n=3$, as demonstrated in Figure 16 for the cases of spatial MetricTensor objects extracted from ADMDecomposition objects for the Schwarzschild and Kerr metrics (with the most general choice of gauge, although this makes no difference to the spacelike hypersurface geometry) presents problems in the case of $3+1$ decomposition and the ADM formalism, since it makes it more difficult to determine whether a given 3-dimensional spacelike hypersurface is conformally-flat (or, indeed, to extract gravitational wave data from 3-dimensional spacelike hypersurfaces directly); indeed, as seen here, the WeylTensor objects will (incorrectly) report both manifolds to be conformally-flat in these cases. To this end, it is instructive instead to consider the rank-3 Cotton tensor $C_{\rho \mu \nu }$ [60], defined by: i.e., in expanded form: for which vanishing of all components $C_{\rho \mu \nu }=0$ is both a necessary and sufficient condition for conformal-flatness in dimension $n=3$ (in dimensions $n>3$, the vanishing of the Cotton tensor $C_{\rho \mu \nu }$ is a necessary but not sufficient condition). Figure 17 illustrates the use of the CottonTensor function in Gravitas to demonstrate that the spatial MetricTensor object extracted from the ADMDecomposition object is conformally-flat for the case of the Schwarzschild metric, but not for the Kerr metric (since it is provable that there do not exist any conformally-flat spatial slices of the Kerr metric [61]). In fact, in the particular case of dimension $n=3$, one can simplify the Cotton tensor $C_{\rho \mu \nu }$ even further by using the Hodge dual operation ★ in the associated Grassmann algebra to reduce the rank-3 Cotton tensor to the rank-2 Cotton-York tensor (density) $C_{\mu }^{\nu }$, defined by: i.e., in expanded form: where ${\epsilon }_{\rho \sigma \nu }$ designates the totally-antisymmetric Levi-Civita symbol of rank-3, and therefore the Cotton-York tensor (density) $C_{\mu }^{\nu }$ is definable only in three dimensions (in which, like the full Cotton tensor $C_{\rho \mu \nu }$, the vanishing of all components of the Cotton-York tensor $C_{\mu }^{\nu }=0$ is both a necessary and sufficient condition for conformal-flatness). Figure 18 illustrates the use of the CottonYorkTensor function in Gravitas to demonstrate, once again, that the spatial MetricTensor object extracted from the ADMDecomposition object is conformally-flat in the Schwarzschild case but not in the Kerr case. The geometrical significance of the Cotton tensor $C_{\rho \mu \nu }$ lies in its very simple transformation law under conformal rescalings of the form $\widetilde{g_{mu\nu }}={\mathsf{\Omega }}^2{g}_{\mu \nu }$ with $\mathsf{\Omega }\in \mathbb{R}$, namely:

In order to be able to formulate the Einstein field equations as a (well-posed) initial value problem, it is necessary not only to specify the (initial) Cauchy data in the form of the (initial) spacelike hypersurface/spatial metric tensor ${\gamma }_{\mu \nu }$, but also to give an appropriate prescription for the gauge, i.e. a sufficient set of mathematical conditions for one to be able to determine (uniquely) the values of the lapse function $\alpha$ and the shift vector field ${\beta }^{\mu }$ at every point in the ambient spacetime. The Gravitas framework makes the terminological distinction between slicing conditions (which are constraints on the lapse function $\alpha$, and which therefore determine geometrically how the spacetime is foliated into a time-ordered sequence of spacelike hypersurfaces) and coordinate conditions (which are constraints on the shift vector ${\beta }^{\mu }$, and which therefore determine how the spatial coordinates are relabeled as one moves between neighboring spacelike hypersurfaces). The most trivial choice of gauge would involve imposing the geodesic slicing condition, in which $\alpha =1$ everywhere, along with the normal coordinate conditions, in which ${\beta }^{\mu }=0$ everywhere; these conditions are extremely straightforward to impose, although they have a tendency to produce coordinate pathologies (and, in particular, make no effort to avoid singularities) and are therefore of limited utility for numerical simulations. In addition to geodesic slicing and normal coordinates, ADMDecomposition includes a small library of in-built gauge conditions (including most standard slicing and coordinate conditions), with many more planned for future inclusion. For example, a very commonly-used slicing condition in numerical relativity (due to its tendency to avoid singularities, and therefore its utility when simulating black hole spacetimes) is the maximal slicing condition, due originally to Lichnerowicz [62] and subsequently developed further by York [40]: where ${}^{\left(3\right)}\mathsf{\Delta }$ designates the induced (connection) Laplacian on spacelike hypersurfaces, where the connection Laplacian $\mathsf{\Delta }$ is defined abstractly (for arbitrary tensor fields T) as the trace of the second covariant derivative ${\nabla }^2$, i.e. $\mathsf{\Delta }T=\mathrm{tr}\left({\nabla }^{2}T\right)$, where the second covariant derivative is itself given abstractly by:

In explicit component form, the action of the induced (connection) Laplacian ${}^{\left(3\right)}\mathsf{\Delta }$ on an arbitrary scalar field $\varphi$ (i.e. a rank-0 tensor field) can therefore be written as: or, equivalently: hence allowing us to express the maximal slicing condition in the following (explicit, expanded) form, by treating the lapse function $\alpha$ as a pure scalar field defined over each spacelike hypersurface: with $\mu ,\nu ,\sigma ,\lambda$ in all of the above ranging across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only). Intuitively, the maximal slicing condition seeks to maximize the spatial volume of each hypersurface by essentially “slowing” the evolution in regions of high curvature and “accelerating” it in regions of low curvature. Figure 19 shows the maximal slicing gauge condition, as computed directly from the ADMDecomposition objects for the Schwarzschild metric (representing, for instance, an uncharged, non-rotating black hole with mass M in Schwarzschild or spherical polar coordinates $\left(t,r,\theta ,\varphi \right)$) and the Kerr metric (representing, for instance, an uncharged, spinning black hole with mass M and angular momentum J in Boyer-Lindquist or oblate spheroidal coordinates $\left(t,r,\theta ,\varphi \right)$), assuming a restricted choice of gauge consisting of the lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and the modified shift vectors $\left(\beta \left(t,r,\theta ,\varphi \right),0,0\right)$ (for Schwarzschild) and $\left(0,0,\beta \left(t,r,\theta ,\varphi \right)\right)$ (for Kerr). Another slicing condition that is frequently used whenever one wishes to cast the Einstein field equations in a strongly hyperbolic form is the harmonic slicing condition first investigated by Bona and Massó [63,64], though its singularity-avoidance properties were discovered by Alcubierre and Massó [65] to be somewhat less favorable than those of maximal slicing, in which one simply asserts that the (spacetime) Laplacian of the time coordinate t (considered now as a pure scalar field on the spacetime) vanishes identically: where ${}^{\left(4\right)}\mathsf{\Delta }$ designates the ambient (connection) Laplacian on spacetime, i.e., in expanded form, one has: with $\mu ,\sigma ,\lambda$ here ranging across all $\left\{0,\cdots ,n-1\right\}$ (i.e. across all spacetime coordinate indices). Intuitively, the harmonic slicing condition replaces the typical elliptic slicing condition (as is the case with, say, maximal slicing) with a corresponding hyperbolic one. Figure 20 demonstrates the harmonic slicing gauge condition, as computed directly from the ADMDecomposition objects for the Schwarzschild and Kerr metrics, with the same restricted choice of gauge (with lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and modified shift vectors $\left(\beta \left(t,r,\theta ,\varphi \right),0,0\right)$ and $\left(0,0,\beta \left(t,r,\theta ,\varphi \right)\right)$, respectively). Finally, a slicing condition that is frequently used in the simulation of compact binary inspirals and mergers [66] [67] (due to its empirical stability properties and comparable singularity-avoidance properties to maximal slicing) is the $1+log$ slicing condition: so-named because it can be rewritten in the slightly expanded form: which, in the particular case where normal coordinate conditions are enforced and in which the shift vector therefore vanishes everywhere, i.e. $\beta =\mathbf{0}$, reduces to: with $\sigma ,\rho$ in all of the above ranging across all $\left\{0,\cdots ,n-2\right\}$ (i.e. across spatial coordinate indices only). $1+log$ slicing is itself a special case of the general Bona-Massó slicing condition presented in  [63] [64]. Figure 21 illustrates the $1+log$ slicing gauge condition, as computed directly from the ADMDecomposition objects for the Schwarzschild and Kerr metrics, with the same restricted choice of gauge (with lapse function $\alpha \left(t,r,\theta ,\varphi \right)$ and modified shift vectors $\left(\beta \left(t,r,\theta ,\varphi \right),0,0\right)$ and $\left(0,0,\beta \left(t,r,\theta ,\varphi \right)\right)$, respectively).