Singularities have been discussed as an essential feature of general relativity more-or-less since its inception: the first non-trivial solution to the Einstein field equations (discovered independently by Schwarzschild [1] and Droste [2] within a year of Einstein’s original publication of the theory in 1915), namely the Schwarzschild metric for a spherically-symmetric matter distribution of mass M, was identified by Hilbert [3] as containing singularities (i.e., points at which certain components of the metric tensor would become divergent) at the coordinate values $r=0$ and $r=2M$. The singularity at $r=2M$ can easily be shown to be a byproduct of the choice of coordinate system, and hence non-physical, as first pointed out by Lemaître [4]. Specifically, although it is present in the case of Schwarzschild’s spherical coordinates, there exist many alternative coordinatizations of the Schwarzschild geometry (such as Eddington-Finkelstein coordinates or Gullstrand-Painlevé coordinates) for which the solution is perfectly regular at $r=2M$ or its equivalent. However, the singularity at $r=0$ appears somehow to be more fundamental: at this point, the Kretschmann scalar $K=R_{\mu \nu \rho \sigma }{R}^{\mu \nu \rho \sigma }$ also diverges, and since K is a quadratic scalar invariant of the Riemann curvature tensor (and thus invariant under arbitrary diffeomorphism transformations of spacetime), the singularity at $r=0$ cannot be merely a coordinate artifact, and must instead be considered an intrinsic feature of the geometry that is independent of coordinates. In 1939, Oppenheimer and Snyder [5] considered the case of an idealized, spherically-symmetric star undergoing continual gravitational collapse, showing that, to an outside observer in Schwarzschild coordinates, the stellar radius would appear to approach the coordinate singularity at $r=2M$ asymptotically, although for an infalling observer in, say, Eddington-Finkelstein coordinates, the stellar radius would appear to cross the coordinate horizon at $r=2M$ within finite proper time and thus continue collapsing indefinitely. Hence, at least within the Oppenheimer-Snyder model of gravitational collapse, the singularity at $r=0$ seems to be physically realized. Nor are such singularities purely a byproduct of spherically-symmetric spacetimes: the axially-symmetric Kerr geometry (as well as the more general Kerr-Newman class of electrovacuum solutions to the Einstein-Maxwell equations) exhibits a spatially-extended, “ring-like” singularity, as well as a pair of Schwarzschild-like coordinate singularities, at least in the oblate spheroidal coordinates of Boyer and Lindquist [6].

However, one might still have argued, perhaps perfectly reasonably, that the apparent existence of such “physical” singularities was nevertheless a consequence of the unphysically high degree of symmetry that these solutions possessed (either spherical in the case of Schwarzschild or axial in the case of Kerr), and that the singularity structure would not remain stable under perturbations away from axial symmetry. In 1964, Penrose [7] showed that this was not, in fact, the case. To begin, Penrose proposed the first mathematically rigorous definition of a generic spacetime singularity: a spacetime $\mathcal{M}$ contains a singularity if it is not geodesically complete, meaning that there exist timelike or null geodesics $\gamma$ in $\mathcal{M}$ that cannot be extended indefinitely either into the future, or into the past, or both (i.e., there exists a limit to how far one can smoothly continue the domain of definition of the time parameter $\tau$ of $\gamma$, for the case of timelike geodesics, or the affine parameter $\lambda$ of $\gamma$, for the case of null geodesics). Next, Penrose assumed a matter distribution exhibiting a sufficiently strong gravitational field that a trapped null/timelike surface (i.e., a compact, spacelike surface at which all outward-pointing null/timelike geodesics become convergent) is induced, but with no a priori assumptions regarding its symmetry. Penrose proceeded to prove that, so long as the spacetime is globally hyperbolic (and thus admits a foliation into non-overlapping spacelike hypersurfaces) and the matter distribution obeys either the null energy condition$T_{\mu \nu }{X}^{\mu }{X}^{\nu }\ge 0$ (for the case of null vector fields $\mathbf{X}$) or the strong energy condition$\left(T_{\mu \nu }-\frac{1}{2}T{g}_{\mu \nu }\right)X^{\mu }{X}^{\nu }\ge 0$ (for the case of timelike vector fields $\mathbf{X}$), the resulting spacetime is necessarily geodesically incomplete. Such energy conditions would naturally be satisfied by any perfect fluid distribution with non-negative density and pressure, as in the case of an idealized fluid model for a collapsing star. Penrose’s original singularity theorem principally concerned the case of future geodesic incompleteness relevant for gravitational collapse models; Hawking [8] subsequently considered the time-reversed case of past geodesic incompleteness relevant for the cosmology of expanding universes (in which the dominant energy condition must assumed - for a perfect fluid, this corresponds to the density being at least as large as the absolute value of the pressure). Many other singularity theorems, concerning different energy conditions or involving different assumptions on the global structure of the spacetime, are now known.

Penrose’s proof is remarkably simple, in large part because it builds so heavily upon the earlier work of Raychaudhuri [9] (much of it independently discovered by Landau [10]) concerning the behavior of geodesic congruences (i.e., families of geodesics sharing the same affine parameter $\lambda$ or time parameter $\tau$) in spacetime. Raychaudhuri’s equation relates the divergence of a geodesic congruence $\mathbf{X}$ to (amongst other things) the so-called Raychaudhuri scalar$R_{\mu \nu }{X}^{\mu }{X}^{\nu }$ (otherwise known as the trace of the tidal tensor, or the trace of the electrogravitic tensor in the context of the Bel decomposition of the Riemann curvature tensor [11,12]). The essence of the proof is then as follows (stated here for the case of null geodesic congruences, but with a straightforward extension to timelike ones). The Raychaudhuri equation implies that if one starts with a congruence of (initially parallel) null geodesics emanating from some spacetime region, and that congruence starts to converge (i.e., the volume of the congruence begins to decrease), then it will continue to converge for as long as the Raychaudhuri scalar is non-negative. By the Einstein field equations, the non-negativity of the Raychaudhuri scalar $R_{\mu \nu }{X}^{\mu }{X}^{\nu }\ge 0$ is equivalent to the satisfaction of the null energy condition $T_{\mu \nu }{X}^{\mu }{X}^{\nu }\ge 0$ on the stress-energy tensor, which holds by assumption. Moreover, the initial convergence of the null geodesics is guaranteed by the assumed presence of the trapped null surface. Taken together, these imply that the null geodesic congruence will “collapse” to have zero volume within some finite value of the affine parameter $\lambda$: a crucial lemma known as the focusing theorem. This “collapse” of the null geodesic congruence implies that all neighboring null geodesics must intersect with one another in some way. However, if two null geodesics intersect, then they cannot lie on the boundary of the proper future of the initial spacetime region (indeed, one way to define the boundary of the proper future is that it is the collection of all null geodesics in the congruence that do not intersect). Thus, if all null geodesics in the congruence intersect with at least one other null geodesic, then the initial spacetime region must have no proper future boundary. Hence, the spacetime is (null) geodesically incomplete, which completes the proof. We see immediately that the Penrose singularity theorem is thus a kind of “physical analog” of the Bonnet-Myers theorem from Riemannian geometry [13], by which any Riemannian manifold whose Ricci curvature is bounded below by a positive constant must either be compact or geodesically incomplete.

In the case of discrete quantum gravity theories such as causal set theory [14,15,16,17] or the Wolfram model [18,19,20,21], in which the underlying structure of spacetime is given in terms of some fundamentally combinatorial data structure (such as a partially ordered set, a directed acyclic graph or a time-ordered sequence of hypergraphs), the situation is somewhat more complicated. Although in many instances agreement with standard (continuum) general relativity is established in cases where an appropriate continuum limit exists (for instance, via the Benincasa-Dowker action from discrete d’Alembertians in causal set theory [22], or the discrete Einstein-Hilbert action from Ollivier-Ricci curvature in Wolfram model systems [20,23,24]), far less is known about the relativistic properties of such theories at the sub-continuum scale, including the validity of the singularity theorems. Indeed, it is far from obvious that singularities are even a well-defined concept in a (finite) discrete spacetime model, since if the number of gravitational degrees of freedom is necessarily finite, then all components and projections of discrete metric and curvature tensors will also be finite. Even the supposedly idealized cases that do not require full singularity theorems in the continuum, such as the spherically-symmetric Oppenheimer-Snyder stellar collapse model, do not admit straightforward translations to the discrete spacetime case, due to the difficulties of enforcing the requisite symmetries. It is possible to construct discrete spacetimes that are compatible with continuous symmetry groups in the continuum limit (for instance, causal sets produced via Poisson sprinkling into a $3+1$-dimensional Minkowski spacetime are provably compatible with the action of the restricted Lorentz group $S{O}^{+}\left(1,3\right)$, in the limit of infinite sprinkling density [25]), but these symmetries do not hold generically when the process is truncated after a finite number of steps. Thus, even if one constructs the Cauchy initial data for the Einstein field equations by sprinkling into a Riemannian manifold exhibiting the requisite spherical symmetry, the inherent randomness of the sprinkling process will have the effect of perturbing the initial data away from perfect spherical symmetry, hence making the methods of Oppenheimer and Snyder fundamentally inapplicable. It is therefore an extremely important, yet highly non-trivial, question to ask whether a result akin to the Penrose singularity theorem can be proven in the discrete case, and what kinds of assumptions (analogous to energy conditions and global causality conditions in the continuum case) might be necessary for such a proof to be possible. Though we do not claim to be able to answer this question with any generality as of yet, this is nevertheless the question with which the present article is concerned.

The Wolfram model is a discrete spacetime formalism in which Cauchy initial data is specified in the form of a spatial hypergraph (i.e., the discrete analog of a spacelike hypersurface), with dynamics determined by hypergraph rewriting rules. The causal interactions between these hypergraph rewrites generate a partially-ordered set, which is typically represented as a directed acyclic graph known as a causal graph; subject to certain assumptions (such as causal invariance and asymptotic dimension preservation), the combinatorial structure of this causal graph is known to converge to the conformal structure of a Lorentzian manifold obeying the Einstein field equations in the continuum limit [20,24]. Causal set theory is a deeply related discrete spacetime model in which spacetime structure is also represented as a partially-ordered set; indeed, Wolfram model evolution may be interpreted as endowing causal set theory with an explicit algorithmic dynamics [23,26,27] (since the transitive reduction of a causal graph generated via Wolfram model evolution corresponds to the Hasse diagram of some corresponding causal set). The problem of defining singularities within such discrete spacetime models turns out to be surprisingly straightforward, since Penrose’s notion of (null) geodesic incompleteness may be imported more-or-less “wholesale”: spacetime geodesics become directed paths in the causal graph, with the question of geodesic completeness being the question of whether such paths can always be extended, or, more concretely, the question of future/past geodesic incompleteness is equivalent to the question of whether there exist vertices within some subgraph of the causal graph (corresponding to the geodesic congruence) whose out/in-components are empty. However, appropriate discrete translations for much of the remaining mathematical apparatus of the Penrose singularity theorem remain elusive, or at the very least obscure. For instance, geodesic “collisions” are far more common in discrete spacetimes than in continuous ones, even in cases where there is no net convergence of the geodesic congruence, for the simple reason that there are fewer degrees of freedom (i.e., there are fewer “spaces” for the geodesic to occupy within a discrete spacetime, so pairs of geodesics are inherently much more likely to occupy the same “space” by “accident”). This makes the intersection and volume properties of geodesic congruences in discrete spacetimes noticeably and qualitatively different than in continuous ones, rendering the appropriate discrete analog of the Raychaudhuri equation somewhat non-obvious. Although considerable attention has been paid to the structure of discrete vacuum spacetimes (especially in the Wolfram model case), far less is known about non-vacuum solutions, and there has as of yet been no systematic investigation of discrete analogs of standard general relativistic energy conditions. Finally, the lack of any a priori inner product structure (although such structures can be defined, albeit not uniquely) on causal graphs makes Penrose’s original definition of trapped null surfaces hard to utilize directly, though equivalent descriptions in terms of cross-sectional areas of null congruences may be more amenable to immediate discretization.

The main result presented within this article is that, for a massive scalar field “bubble collapse” problem obeying one of two approximate spatial symmetries, when formulated and discretized as a Wolfram model evolution problem, the resulting discrete spacetime converges to one of two standard black hole spacetimes (namely either Schwarzschild or extremal Kerr). Specifically, if the initial scalar field distribution is approximately spherically-symmetric, then the resulting discrete spacetime is asymptotically Schwarzschild, whereas if the distribution is instead approximately axially-symmetric, then the resulting discrete spacetime is asymptotically extremal Kerr (reducing to Schwarzschild in degenerate cases). The significance of the word “approximately” in the above is in reference to the symmetry discretization problem referenced previously; although the Cauchy initial data is exactly spherically/axially-symmetric analytically, the discretization procedure inevitably has the effect of perturbing the initial hypersurface away from exact spatial symmetry, and so one must show that the resulting convergence remains stable under perturbations of this general form. We show this stable convergence property using both rigorous mathematical analysis and explicit numerical simulation (by means of a newly-developed, high-resolution, hypergraph-based numerical relativity code called Gravitas, featuring totally unstructured adaptive mesh refinement), and demonstrate the expected agreement between the analytical and numerical results. Our rationale for choosing a massive scalar field as the underlying stress-energy model (rather than the dust or perfect fluid models commonly used in relativistic astrophysics) is to enable more direct comparison with pure Wolfram model evolution: in recent work [28], we showed how a massless scalar field theory (obeying the discrete Klein-Gordon equation) could be defined over an arbitrary Wolfram model system, building upon the previous work of Dowker and Glaser [29], Sorkin [30] and Johnston [31] in the context of causal set theory. Following an ansatz proposed by Dowker et al. [32], as well as a more direct approach outlined by Johnston[33,34], the discrete massless Green’s functions can naturally be extended to the massive case. To the best of our knowledge, no comparable proposal has yet been made for equipping arbitrary causal sets/Wolfram model evolutions with matter fields consistent with either the dust or perfect fluid forms of the continuum stress-energy tensor. However, as we shall show through the course of this article, the WKB approximation of Wentze l [35], Kramers [36] and Brillouin [37] may nevertheless be applied to find a range of parameter values in which a collapsing massive scalar field bubble may be interpreted as a non-rotating, collapsing ball of dust (described by the Lemaître-Tolman-Bondi metric [4,38,39]) in the spherically-symmetric case, or a spinning, collapsing disk of dust (described by the Weyl-Lewis-Papapetrou metric [40,41,42]) in the axially-symmetric case.

In Section 2, we begin by introducing the fully covariant and conformally-invariant Z4 formulation (CCZ4) of the Einstein field equations due to Alic et al. [43] and Bona et al. [44] used by Gravitas in the formulation of the Cauchy problem for general relativity (along with the relevant gauge conditions, adapted for the spherically-symmetric and axially-symmetric spacetimes simulated in this article). We also outline the various numerical algorithms used in the discrete evolution of these equations, including the fourth-order Runge-Kutta scheme for the time evolution (with appropriate modification of the finite-difference stencils to make them suitable for general hypergraphs with totally unstructured topology), the generalized local adaptive mesh refinement (AMR) algorithm, based on the approach of Berger and Colella [45], for coarsening and refining the hypergraph topology, and the higher-order weighted essentially non-oscillatory (WENO) scheme for extrapolating boundary values. In Section 3, we follow the approach of Gonçalves and Moss [46] to show that the WKB approximation may be used to treat a minimally-coupled massive scalar field in spherical symmetry as a non-rotating inhomogeneous dust described by the Lemaître-Tolman-Bondi metric, at least in the limit as the mass of the field goes to infinity, and thus we rigorously derive a sufficient condition for the spherically-symmetric massive scalar field “bubble collapse” problem to yield the idealized stellar collapse solution of Oppenheimer and Snyder. This condition can be derived analytically for a “top hat” initial density distribution of the scalar field, and numerically for the (more physical) exponential initial density distribution. We also present the numerical solutions to the massive scalar field bubble collapse problem in spherical symmetry, showing agreement with the analytic predictions, within this section. However, the collapse of an axially-symmetric massive scalar field to an extremal Kerr black hole is inherently more complicated, since the exterior solution for a spinning disk of dust is not extremal Kerr, but rather the Weyl-Lewis-Papapetrou metric (in contrast to the spherically-symmetric case, in which the exterior solutions for both the collapsing dust and the resulting non-rotating black hole are described by the Schwarzschild geometry, by virtue of Birkhoff’s theorem), so one must attempt to construct some kind of smooth transition between the two distinct geometries. In Section 4, we follow the methods of Neugebauer and Meinel [47,48,49] (based on an earlier conjecture of Bardeen and Wagoner [50,51]) to prove that the collapse of a minimally-coupled massive scalar field bubble in axial symmetry (and thus, by the WKB approximation, the collapse of a uniformly-spinning disk of dust) may be formulated as a boundary-value problem for the Ernst equation, whose solution can be obtained as a special case of the Jacobi inversion problem for Abelian/hyperelliptic integrals. The full solution may be given in terms of so-called ultraelliptic functions (i.e., hyperelliptic functions of two variables) and ultraelliptic integrals, and ultimately depends upon both the the angular velocity of the disk and the relative redshift from the disk’s center, with the Maclaurin disk solution (i.e., a gaseous astrophysical disk under Newtonian gravity) and the extremal Kerr solution being obtained in the appropriate “Newtonian” and “exterior” limits, respectively. The numerical solutions to the massive scalar field bubble collapse problem in axial symmetry, illustrating agreement with the analytic solutions, are also presented within this section. In Section 5, we outline how the approach of Johnston may be used to equip an arbitrary causal graph produced by Wolfram model evolution with the Green’s function for a massive scalar field, thus allowing us to make a direct comparison between the numerical results obtained within the preceding sections and those produced via “pure” Wolfram model evolution (without any underlying PDE system), using a hypergraph rewriting rule that provably satisfies the Einstein field equations in the continuum limit. Finally, in Section 6, some potential directions for future research are proposed, including the extension to more complex matter/stress-energy models involving more sophisticated equations of state, the extension to more general classes of perturbations away from spherical or axial symmetry, and the more speculative possibility of extending the kinds of global topological techniques of Penrose in order to prove a fully generalized singularity theorem for arbitrary discrete spacetimes.

Note that all of the code necessary to reproduce all of the results presented within this article is fully open source and freely available as part of the Gravitas package on GitHub. In particular, Gravitas features in-built functionality for performing the canonical $3+1$ metric decomposition and computing the corresponding evolution equations, constraint equations, gauge satisfaction equations, etc. (e.g., through ADMDecomposition), discretizing the resulting evolution using adaptive hypergraph methods (e.g., through DiscreteHypersurfaceDecomposition), coupling arbitrary spacetimes to arbitrary matter/stress-energy distributions (e.g., through StressEnergyTensor) and solving the Einstein field equations, both numerically and analytically (e.g., through SolveEinsteinEquations and SolveVacuumEinsteinEquations). Gravitas also contains an extensive library of in-built metrics, geometries, gauge choices and coordinate systems. Several functions are also fully-documented and exposed through the Wolfram Function Repository, including MetricTensor for representing arbitrary spacetime metrics, RicciTensor for computing Ricci curvature tensors (and their projections) for arbitrary spacetimes, MultiwaySystem for evolving arbitrary Wolfram model (multiway) systems and computing their corresponding causal graphs, etc. Note also that, throughout this article, we adopt the general convention that Greek indices (i.e., $\mu$, $\nu$, $\rho$, $\sigma$, etc.) correspond to spacetime coordinates, while Latin indices (i.e., i, j, k, l, etc.) correspond to spatial coordinates. Einstein summation convention is assumed throughout, unless otherwise specified.

The starting point for our numerical implementation is the canonical $3+1$ decomposition of the spacetime metric tensor $g_{\mu \nu }$ into the standard ADM line element of Arnowitt, Deser and Misner [52]: where the lapse function $\alpha$ and the shift vector $\beta$ are the ADM gauge variables (defining a foliation of a globally hyperbolic spacetime into non-overlapping spacelike hypersurfaces), and ${\gamma }_{ij}$ denotes the induced spatial metric tensor on the spacelike hypersurfaces of constant time. Within this foliation, the future-pointing (null) unit normal vector $\mathbf{n}$ to the spacelike hypersurfaces is given in terms of the contravariant derivative ${\nabla }^{\mu }$ of the time coordinate t: with the corresponding time vector $\mathbf{t}$ given by:This decomposition allows us to formulate the Einstein field equations as a Cauchy initial-value problem on spacelike hypersurfaces, defined by a system of 12 independent hyperbolic partial differential equations: 6 for the spatial metric tensor ${\gamma }_{ij}$: and 6 for the extrinsic curvature tensor on spacelike hypersurfaces $K_{ij}$:In the above, $K_{ij}$, i.e., the extrinsic curvature tensor, may be defined abstractly in terms of the Lie derivative of the spatial metric tensor ${\gamma }_{ij}$ along the normal vector $\mathbf{n}$: K denotes its trace (i.e., the mean curvature): and $S_{ij}$, S and $\tau$ are simple functions of the stress-energy tensor $T_{\mu \nu }$:The resulting system is closed by imposing appropriate conditions on the ADM Hamiltonian/energy constraint H and the three ADM momentum constraints $M_i$, given by: and: respectively. Here, as elsewhere, the bracketed integers are used to distinguish between spatial vs. spacetime variants of expressions and operations: ${\nabla }^{\left(3\right)}$ designates a covariant derivative over a spacelike hypersurface, ${\nabla }^{\left(4\right)}$ designates a covariant derivative over the entire spacetime, $R^{\left(3\right)}$ designates the Ricci scalar curvature restricted to a spacelike hypersurface, etc. Unbracketed expressions and operations are assumed to refer to their spacetime variants unless otherwise specified.

As first proposed by Bona, Ledvinka, Palenzuela and Z̆ác̆ek [53], one can construct an explicitly covariant formulation of the Einstein field equations with Lagrangian density [54]: for some 4-vector $\mathbf{Z}$, such that ADM Hamiltonian and momentum constraints reduce to the trivial algebraic condition $Z_{\mu }=0$, and the Einstein field equations themselves take the form:This is commonly known as the (fully covariant) Z4 formulation. The components of the 4-vector $Z_{\mu }$ may therefore be thought of as quantifying the degree to which the ADM constraints are violated; as shown by Gundlach, Martín-García, Calabrese and Hinder[55], one can introduce damping terms into the covariant Z4 equations so as to formulate the evolution equations as a so-called “$\lambda$-system” based on the techniques of Brodbeck, Frittelli, Hübner and Reula[56]. More specifically, if one replaces the Z4 evolution equations with the damped form: where ${\kappa }_1\ge 0$ and ${\kappa }_2$ are real constants, and $t^{\mu }$ is any non-vanishing null vector field (this effectively translates to a weak constraint on the permitted choices of ADM gauge conditions), then the surface in phase space on which the constraint equations $Z_{\mu }=0$ are satisfied will provably form a set of attractor solutions for the corresponding $\lambda$-system, thus guaranteeing symmetric hyperbolicity of the (constrained) evolution equations. The explicit time evolution equations for the (damped) Z4 system, in terms of the canonical $3+1$/ADM decomposition, can now be expressed reasonably succinctly via Lie derivatives along the shift vector $\beta$. The 6 evolution equations for the spatial metric tensor ${\gamma }_{ij}$ remain unchanged: although the 6 evolution equations for the extrinsic curvature tensor $K_{ij}$ now inherit a non-trivial dependency on the $\mathbf{Z}$ 4-vector (as well as the damping terms ${\kappa }_1$ and ${\kappa }_2$): where $\Theta$ is simply the projection of the $\mathbf{Z}$ 4-vector in the normal direction $\mathbf{n}$:This projection itself is subject to an evolution equation: as indeed is the overall $\mathbf{Z}$ 4-vector: which together close the system. Here, the ${\alpha }_i$ terms correspond to the constants ${\alpha }_i\ne 0$ that appear in the general definition of the $\lambda$-system for constraint quantities in the Einstein field equations: where the ${\beta }_i>0$ are also constants, $\mathcal{C}$, ${\mathcal{C}}^i$, ${\mathcal{C}}_k^{ij}$ and ${\mathcal{C}}_{kl}^{ij}$ are scalar, vector, rank-3 tensor and rank-4 tensor-valued constraint quantities, respectively, and $\lambda$, ${\lambda }^i$, ${\lambda }_k^{ij}$ and ${\lambda }_{kl}^{ij}$ are corresponding tensorial quantities with the same symmetries. Note also that the functions of the stress-energy tensor $S_{ij}$, $S_i$ and $\tau$ are subtly different for the covariant Z4 system from their counterparts in conventional ADM (indeed, one of them has been now been promoted from a scalar to a covector):

The fully conformal and covariant formulation of Z4 (otherwise known as CCZ4) is thus obtained by performing a rescaling of the spatial metric tensor ${\gamma }_{ij}$ by a conformal factor $\varphi$ [57,58]: chosen such as to guarantee unit determinant for the conformal spatial metric ${\tilde{\gamma }}_{ij}$:The trace-free components ${\tilde{A}}_{ij}$ of the total extrinsic curvature tensor $K_{ij}$ on spacelike hypersurfaces then become (after conformal rescaling):By introducing the Christoffel symbols ${\tilde{\Gamma }}_{jk}^i$ associated to the conformal spatial metric tensor ${\gamma }_{ij}$, as well as the pseudovector ${\tilde{\Gamma }}^i$ derived from the contraction of this conformal connection, namely: and: respectively, we can further decompose the Ricci curvature tensor $R_{ij}^{\left(3\right)}$ on spacelike hypersurfaces into a sum: of a rank-2 tensor $R_{ij}^{\left(3\right)}$ depending only upon spatial derivatives of the conformal spatial metric tensor ${\tilde{\gamma }}_{ij}$: and a second rank-2 tensor ${\tilde{R}}_{ij}^{\varphi \left(3\right)}$ depending only upon derivatives of the conformal factor $\varphi$:The fully conformal and covariant Z4 formulation now decomposes into a system of 6 evolution equations for the conformal spatial metric tensor ${\tilde{\gamma }}_{ij}$ (written here explicitly, without Lie derivatives): 6 evolution equations for the trace-free components of the (conformal) extrinsic curvature tensor ${\tilde{A}}_{ij}$: and single evolution equations for the trace of the extrinsic curvature tensor K: the conformal factor $\varphi$: and the projection $\Theta$ of the $\mathbf{Z}$ 4-vector in the normal direction $\mathbf{n}$:Finally, in order to close the system, we must introduce a set of 4 evolution equations for the components of the pseudovector ${\tilde{\Gamma }}^i$ derived from the conformal connection coefficients ${\tilde{\Gamma }}_{jk}^i$; for simplicity, we represent these as equations for the components ${\widehat{\Gamma }}^i$, defined as: namely:

Note that the new ${\kappa }_3$ parameter in the evolution equations for the conformal connection coefficients ${\widehat{\Gamma }}^i$ above is currently unconstrained; although ${\kappa }_3=1$ yields a fully covariant formulation (at the cost of numerical stability in the case of black hole spacetimes) and ${\kappa }_3=\frac{1}{2}$ yields a non-covariant formulation which is nevertheless numerically stable, we choose instead to follow the prescription of Alic, Kastaun and Rezzolla [59] and eliminate the numerical instabilities by choosing ${\kappa }_1\to \frac{{\kappa }_1}{\alpha }$ as the choice of damping parameter, which preserves spatial covariance in all terms depending on ${\kappa }_1$, and full spacetime covariance for all other terms. Moreover, for simulations involving spherical symmetry, we employ the maximal slicing gauge condition [60,61] on the lapse function $\alpha$: and the more general $1+log$ slicing gauge condition[62] for the axially-symmetric case: due to their known stability (and singularity-avoidance) properties for black hole spacetimes. By the same token, we use the gamma-driver gauge conditions [63,64] on the components of the shift vector ${\beta }^i$: or, in the simplified case (without any advection terms): where $B^i$ is an auxiliary vector field, and ${\eta }_1$, ${\eta }_2$, ${\mu }_{{\beta }_1}$ and ${\mu }_{{\beta }_2}$ are scalar parameters, chosen here to be ${\eta }_1=\frac{3}{4}$, ${\mu }_{{\beta }_1}=1$, ${\mu }_{{\beta }_2}=0$ and ${\eta }_2=1$ (i.e., the standard hyperbolic driver conditions). The standard conformal BSSN formalism of Baumgarte, Shapiro, Shibata and Nakamura [65] is recovered in the special case where $\Theta =0$ and $Z^i=0$, the trace of the extrinsic curvature tensor K is replaced with $K^{BSSN}$, defined by: and the Hamiltonian constraint is used to eliminate any explicit dependence on the (spatial) Ricci scalar R from the evolution equation for the trace of the extrinsic curvature tensor K. The consistency conditions on the two additional scalar fields introduced by the conformal rescaling of the covariant Z4 system (namely $\det \left({\tilde{\gamma }}_{ij}\right)$ and $\mathrm{tr}\left({\tilde{A}}_{ij}\right)$): are enforced explicitly as part of the numerical integration algorithm at each time step (by subtracting out the trace from ${\tilde{A}}_{ij}$ and rescaling the conformal spatial metric tensor ${\tilde{\gamma }}_{ij}$ accordingly). The singularity-avoidance properties of our chosen gauge conditions remove any need for us to employ techniques such as excision or punctures for the case of the symmetrical black hole spacetimes simulated within this article.

With regards to the governing equations for our numerical simulations, all that remains is to specify how our spacetime may be equipped with a minimally-coupled (massive) scalar field $\Phi$ obeying the (massive) Klein-Gordon equation: or for a more general potential $V\left(\Phi \right)$ (with the massive case hence corresponding to $V\left(\Phi \right)=m^2\overline{\Phi }\Phi$): which, in curved spacetime, becomes:The solution to the Klein-Gordon equation yields a stress-energy tensor $T^{\mu \nu }$ with components: in the massive case, or, in lowered-index form $T_{\mu \nu }$: in the case of a general potential. Within the initial-value formulation of the Einstein field equations afforded by the ADM decomposition, the (hyperbolic) evolution equation for this scalar field $\Phi$ is then given by: where ${\Pi }_m$ is the (negative of the) conjugate momentum of the scalar field given by:

However, we can immediately see that this “definition” of the (negative) conjugate momentum ${\Pi }_m$ is really just a trivial restatement of the “evolution equation” for the scalar field $\Phi$; therefore, the “true” evolution equation for the scalar field must be given in terms of an evolution equation for its (negative) conjugate momentum ${\Pi }_m$ directly, instead, i.e:

In order to evolve the resulting system of equations using an explicit finite-difference time-stepping method, it is first necessary to introduce a means of defining finite-difference numerical stencils over hypergraphs of arbitrary topology (i.e., hypergraphs with no a priori spatial coordinate structure). As outlined in previous work [24], this can be achieved by applying a generalized bicubic: or tricubic [66]: interpolation scheme to the neighborhoods at the endpoints of hypergraph geodesics, depending upon whether the hypergraph has a limiting 2-dimensional or 3-dimensional structure, respectively (and where the number of “a” coefficients is therefore either $2^4$ or $2^6$), as shown in Figure 1. By equipping the hypergraph with a local inner product structure using discrete projections, as shown in Figure 2 (with the axioms of linearity and conjugate symmetry being satisfied whenever the hypergraph converges to a Riemannian manifold in the continuum limit), one can therefore define a local set of orthonormal coordinate axes on the hypergraph, over which a typical finite-difference scheme can be formulated. For the simulation results presented within this article, we choose to use an explicit fourth-order Runge-Kutta numerical scheme, which, for an initial value problem of the form: where $\Phi \left(\mathbf{x},t\right)$ is a generic state vector, i.e: we can evolve the system forwards in units of a discrete time step $\Delta t$ (chosen so as to be compatible with the Courant condition, and therefore numerically stable) as [67,68]:

Such an initial-value ODE problem can be obtained from a system of non-linear (hyperbolic) PDEs, such as the Einstein field equations, of the general form: where $\mathcal{F}$ is a non-linear operator acting on $\Phi$, by computing a new set of fluxes between vertices in the hypergraph (thus linearizing the problem) at each time step. Note that this assumes a set of local spatial coordinates $\mathbf{x}$ and a global time coordinate t for the hypergraphs, which can be obtained using the geometrical construction described above.

In the explicit time evolution formula above, we choose: and: which completes the temporal discretization of the scheme. For spatial discretization, we choose a coordinate structure such that $\Delta x=x_{i+1}-x_i$, $\Delta y=y_{j+1}-y_j$ and $\Delta z=z_{k+1}-z_k$ (i.e., the three discrete coordinates are indexed by i, j and k, respectively), leading to a discrete form of the state vector ${\Phi }_{i,j,k}^n$ defined at each vertex and for each discrete time step. The spatial derivatives may then be computed by means of the fourth-order centered finite-difference stencils of Zlochower, Baker, Campanelli and Lousto [70], with first derivatives (in x) defined by: and likewise for all other first derivatives. Note that, whenever advection terms (i.e., terms of the form ${\beta }^i{\partial }_{i}F\left(t_n,{\Phi }_{i,j,k}^n\right)$) are present, one must instead use a fourth-order upwind finite-difference scheme: whenever ${\beta }^x<0$, and: whenever ${\beta }^x>0$, and likewise for all other first derivatives. The finite-difference stencils for second derivatives (both for derivatives purely in x, and for mixed derivatives in x and y) can be deduced by applying the first derivative stencils sequentially, in any order, for instance yielding: and: respectively, for the fourth-order centered case, and likewise for the fourth-order upwind case (and for all other second derivatives). We also incorporate a dissipation term of Kreiss-Oliger type [71] into the definition of the (discrete) time derivative for the fourth-order finite-difference scheme, in order to reduce the probability of spurious high-frequency modes destabilizing the numerical solution: for derivatives in the x-direction, and likewise for all other directional derivatives.

As our adaptive refinement algorithm, we use the hypergraph generalization of the local adaptive mesh refinement (AMR) algorithm proposed by Berger and Colella [45], based on previous methods developed by Berger and Oliger [72] and Gropp [73]. Broadly speaking, we introduce a tagging function $f\left(x,y,z\right)$ based on whether the $L_2$ norm of the change in some scalar field $\varphi$ (which, for the simulations presented within this article, will always be chosen to be a Riemannian scalar curvature invariant) exceeds a pre-defined threshold $\sigma \left(\varphi \right)$ across a given vertex or subhypergraph: such that $f\left(x,y,z\right)=1$ signifies that a given vertex or subhypergraph is tagged for refinement, and $f\left(x,y,z\right)=0$ otherwise. The entire hypergraph is then partitioned into subhypergraphs, the signatures$X\left(x\right)$, $Y\left(y\right)$ and $Z\left(z\right)$ are computed for each subhypergraph: and the Laplacians ${\partial }_x^{2}X\left(x\right)$, ${\partial }_y^{2}Y\left(y\right)$ and ${\partial }_z^{2}Z\left(z\right)$ of these signatures determine the axis that will separate the tagged and untagged vertices in the direction orthogonal to the signature function (i.e., the axis of partition), since this will be the local coordinate direction which maximizes $\Delta \left({\partial }_i^2{X}_i\right)$ (and hence which induces an appropriate inflection point in the Laplacian). If the signature is identically zero in any direction (i.e., if $X_i\left(x_i\right)=0$), then that direction is chosen as the axis of partition instead. So long as any given subhypergraph satisfies the two principal axioms of hierarchical grid/hypergraph nesting (i.e., that fine subhypergraphs must always be contained within the neighborhood of a vertex in the next coarsest subhypergraph, and that a vertex that is not at the boundary of the domain at refinement level l must be separated from a vertex at refinement level $l-2$ by at least one vertex at refinement level $l-1$, in any direction in the neighborhood), and the ratio of tagged vertices to total vertices is at least $ϵ$ (for some predefined $0<ϵ<1$), the partitioning procedure terminates; otherwise, it continues on recursively subdividing subhypergraphs into further subhypergraphs and the procedure starts again. The refinement/coarsening procedure for arbitrary vertices/hypergraphs, which can be enacted via the operations of hyperedge subdivision and hyperedge smoothing, respectively, is illustrated in Figure 3.

Following refinement, the time evolution scheme for the finite-difference method, which may be represented in a very explicit (and manifestly conservative) form as: with F, G and H denoting the inter-vertex flux functions projected in the x, y and z directions, respectively, and with ${\Phi }_{i\pm \frac{1}{2},j,k}$, ${\Phi }_{i,j\pm \frac{1}{2},k}$ and ${\Phi }_{i,j,k\pm \frac{1}{2}}$ denoting the extrapolated values of the state vector $\Phi$ at the x-boundaries, y-boundaries and z-boundaries, respectively, must now be modified so as to account for the presence of new boundaries between coarse and fine subhypergraphs. Whenever a coarse vertex (i.e., a vertex at refinement level $l-1$) is overlaid by a refined subhypergraph (i.e., a subhypergraph at refinement level l), the values of the coarse state vector ${\Phi }^{coarse}$ (i.e., the state vector at refinement level $l-1$) are defined in terms of (conservative) averages of the values ${\Phi }^{fine}$ in the fine subhypergraph (i.e., the state vectors at refinement level l): where i, j and k are the discrete coordinates of the coarse vertex, where the same vertex spans the discrete intervals $\left[m,m+r-1\right]$, $\left[n,n+r-1\right]$ and $\left[o,o+r-1\right]$ within the fine subhypergraph. On the other hand, whenever a coarse vertex (at refinement level $l-1$) is adjacent to an interface with a fine subhypergraph (at refinement level l), it is necessary to modify the time evolution scheme to be of the following general form (assuming that the coarse-to-fine interface is projected in the x direction, with the obvious modifications for the other coordinate directions): with $\Delta t_{coarse}$ and $\Delta t_{fine}$ denoting the stable time steps at refinement levels $l-1$ and l, respectively, scaled using the same refinement ratios as the vertex sizes themselves, i.e: and with $\Delta x$, $\Delta y$ and $\Delta z$ designating the spatial size of the vertices in the coarse subhypergraph.

This modification can be implemented as a corrector step to the naïve flux computation in the coarse subhypergraph, whereby the coarse fluxes are subtracted from ${\Phi }_{i,j,k}^{coarse}\left(t+\Delta t_{coarse}\right)$ and replaced with the corresponding fine ones. One starts by constructing a tensor $\delta F$ of fluxes through all coarse subhypergraph boundaries that are also outer boundaries of fine subhypergraphs, i.e., in the x direction: and likewise for the other coordinate directions. After each stable time step $\Delta t_{fine}$ in the fine subhypergraph, we then add a sum of all the fine subhypergraph fluxes through the boundary at discrete coordinates $\left(i+\frac{1}{2},j,k\right)$: and likewise for the other discrete coordinate projections $\left(i,j+\frac{1}{2},k\right)$ and $\left(i,j,k+\frac{1}{2}\right)$. After r such time steps $\Delta t_{fine}$ have elapsed, the flux tensor $\delta F^{new}\left({\Phi }_{i+\frac{1}{2},j,k}\right)$ can be used to correct the solutions over the coarse subhypergraph so as to match those computed using the modified time evolution formula presented above; for instance, for a vertex with discrete coordinates $\left(i+1,j,k\right)$ only one correction is needed: whereas two corrections would be needed for the vertex at discrete coordinates $\left(i+2,j,k\right)$: etc. An example of a refinement of a two-dimensional spatial hypergraph with a regular quadrilateral grid structure into a collection of 2-by-2 two-dimensional grids is shown in Figure 4, along with the corresponding configuration of coarse vertices for which the finite-difference evolution scheme must be modified to account for the adjacencies between coarse and fine subhypergraphs in Figure 5.

In order to extrapolate the boundary values of the state vector $\Phi$, i.e., ${\Phi }_{i\pm \frac{1}{2},j,k}$, ${\Phi }_{i,j\pm \frac{1}{2},k}$ and ${\Phi }_{i,j,k\pm \frac{1}{2}}$ (and hence also to compute the inter-vertex flux functions $F\left({\Phi }_{i\pm \frac{1}{2},j,k}\right)$, $G\left({\Phi }_{i,j\pm \frac{1}{2},k}\right)$ and $H\left({\Phi }_{i,j,k\pm \frac{1}{2}}\right)$), we choose to use a weighted, essentially non-oscillatory (WENO) spatial reconstruction scheme [74,75] with the required fourth-order accuracy. For this purpose, we first choose a set of $M+1$ linearly-independent polynomials in the nodal basis, denoted ${\left\{{\Psi }_l\right\}}_{l=1}M+1$; the Laguerre polynomials ${\ell }_j\left(x\right)$ constitute a natural such choice: where $1\le j\le M+1$, interpolating between the set of $M+1$nodal points$\left(x_j,y_j,z_j\right)$ from $\left(x_1,y_1,z_1\right)$ to $\left(x_{M+1},y_{M+1},z_{M+1}\right)$ in such a way that no two $x_j$ are ever the same. In this context, the term nodal basis indicates that the polynomials must all be of degree M (as the Laguerre polynomials are), and that the basis functions have been rescaled so as to fit within the unit reference interval $\left[0,1\right]$ using the coordinate transformation $\left(x,y,z\right)\to \left(\xi ,\eta ,\zeta \right)$:We denote the $M+1$ nodal points $\left(x_j,y_j,z_j\right)$ using the shorthand ${\left\{x_k\right\}}_{k=1}^{M+1}$, such that:Using the numerical stencils ${\mathcal{S}}_{i,j,k}$, defined for each of the Cartesian coordinate directions: where $I_{i,j,k}$ denotes the vertex with discrete coordinates $\left(i,j,k\right)$, and where L and R designate the spatial extent of the numerical stencil to the left and right, respectively, it is now possible to perform the required spatial reconstruction. Since we wish in this particular case for the spatial reconstruction to be fourth-order accurate (i.e., $M+1=4$), this implies that the basis polynomials M will all be of odd degree (i.e., $M=3$), and so four numerical stencils are used: two centered (with $s=0$, $L=⌊\frac{M}{2}⌋+1$, $R=⌊\frac{M}{2}⌋$ and $s=1$, $L=⌊\frac{M}{2}⌋$, $R=⌊\frac{M}{2}⌋+1$, respectively), one left-sided (with $s=2$, $L=M$, $R=0$) and one right-sided (with $s=3$, $L=0$, $R=M$). The spatial reconstruction procedure itself works by taking the second-order boundary extrapolation step from the SLIC/MUSCL-Hancock finite-volume approaches [76,77], namely (assuming extrapolation in the x coordinate direction): for left boundary extrapolation and: for right boundary extrapolation, with (diagonal) limiter functions ${\varphi }^{\pm }$, and replacing the linearized approximation scheme on the right-hand side with a higher-degree spatial reconstruction polynomial ${\mathbf{w}}_h^{s,x}\left(x,t^n\right)$. This polynomial (again assuming reconstruction in the x coordinate direction) may be expanded in terms of the nodal basis polynomials ${\Psi }_l\left(\xi \right)$ as: where we have made use of the Einstein summation convention within the second equality. By imposing the weak (i.e., integral) form of the conservation equations across every vertex within the numerical stencil ${\mathcal{S}}_{i,j,k}^{s,x}$, we obtain the following (linear) system of integral equations: where ${\overline{\Phi }}_{i,j,k}^n$ denotes the (spatial) average of the solution vector $\Phi$, when integrated over the vertex $I_{i,j,k}$ at time $t^n$. The coefficients ${\widehat{\mathbf{w}}}_{i,j,k,p}^{n,s}$ of the reconstruction polynomials for each stencil may thus be obtained by solving the resulting linear system, and so the reconstruction for the entire vertex can be performed by taking the following non-linear combination of the polynomials for each stencil: where we set $N_s=4$ since the polynomial degree M is odd, and the data-dependent non-linear weights ${\omega }_s$ are defined by:In the above, we choose ${\sigma }_s$ to be an oscillation indicator function: where ${\Sigma }_{lm}$ is an oscillation indicator matrix [78]: so as to damp the effects of spurious numerical oscillations in the solution (hence making the resulting scheme essentially non-oscillatory, as required). Based on the outcomes of numerical experimentation, for the remainder of this article we choose ${\lambda }_s=1$ and ${\lambda }_s={10}^5$ for one-sided and centered stencils, respectively, $ϵ={10}^{-4}$ (to prevent division by zero in the case of smooth solutions) and $r=8$. For a more complete description of the hypergraph-based numerical algorithms used within this article, we invite the reader to consult [24].