Jonathan Gorard
email: gorardj@cardiff.ac.uk

Non-Vacuum Solutions, Gravitational Collapse and Discrete Singularity Theorems in Wolfram Model Systems

The celebrated geodesic congruence equation of Raychaudhuri, together with the resulting singularity theorems of Penrose and Hawking that it enabled, yield a highly general set of conditions under which a spacetime (or, more generically, a pseudo-Riemannian manifold) is expected to become geodesically incomplete. However, the proofs of these theorems traditionally depend upon a collection of assumptions about the continuum spacetime (and, in the physical case, the stress-energy distribution defined over it), including its global structure, its energy conditions, the existence of trapped null surfaces and the

various volume/intersection properties of geodesic congruences, that are inherently difficult to translate to the case of discrete spacetime formalisms, such as causal set theory or the Wolfram model. Some, such as the discrete analog of Raychaudhuri's equation for the volumes of geodesic congruences, are subtle to formulate due to intrinsic differences in the behavior of discrete vs. continuous geodesics; for others, such as the definition of a trapped null surface, no appropriate translation is known for the general discrete case due to a lack of a priori coordinate information. It is therefore a non-trivial question to ask whether (and to what extent) there exist equivalently general conditions under which one expects discrete spacetimes to become geodesically incomplete, and how these conditions might differ from those in the continuum. This article builds upon previous work, in which the conformal and covariant Z4 (CCZ4) formulation of the Cauchy problem for the Einstein field equations, with constraint-violation damping, was defined in

Introduction

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 [START_REF] Schwarzschild | Über das Gravitationsfeld eines Massenpunktes nach der Einsteinschen Theorie[END_REF] and Droste [START_REF] Droste | The field of a single centre in Einstein's theory of gravitation, and the motion of a particle in that field[END_REF] 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 [START_REF] Hilbert | Die Grundlagen der Physik[END_REF] 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 [START_REF] Lemaître | L'Univers en expansion[END_REF]. 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 µνρσ R µνρσ 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 [START_REF] Oppenheimer | On Continued Gravitational Contraction[END_REF] 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 [START_REF] Boyer | Maximal Analytic Extension of the Kerr Metric[END_REF].

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 [START_REF] Penrose | Gravitational Collapse and Space-Time Singularities[END_REF] 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 M contains a singularity if it is not geodesically complete, meaning that there exist timelike or null geodesics γ in 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 τ of γ, for the case of timelike geodesics, or the affine parameter λ of γ, 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 µν X µ X ν ≥ 0 (for the case of null vector fields X) or the strong energy condition T µν - 1 2 T g µν X µ X ν ≥ 0 (for the case of timelike vector fields 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 [START_REF] Hawking | Properties of expanding universes[END_REF] 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 [START_REF] Raychaudhuri | Relativistic Cosmology. I[END_REF] (much of it independently discovered by Landau [START_REF] Hawking | The Large Scale Structure of Space-Time[END_REF]) concerning the behavior of geodesic congruences (i.e. families of geodesics sharing the same affine parameter λ or time parameter τ ) in spacetime.

Raychaudhuri's equation relates the divergence of a geodesic congruence X to (amongst other things) the so-called Raychaudhuri scalar R µν X µ X ν (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 [START_REF] Matte | Sur de nouvelles solutions oscillatoires des equations de la gravitation[END_REF] [START_REF] Bel | Définition d'une densité d'énergie et d'un état de radiation totale généralisée[END_REF]).

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 µν X µ X ν ≥ 0 is equivalent to the satisfaction of the null energy condition T µν X µ X ν ≥ 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 λ: 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 [START_REF] Myers | Riemannian manifolds with positive mean curvature[END_REF], 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 [START_REF] Bombelli | Space-time as a causal set[END_REF][15] [START_REF] Sorkin | Spacetime and Causal Sets[END_REF] [START_REF] Sorkin | Forks in the Road on the Way to Quantum Gravity[END_REF] or the Wolfram model [START_REF] Wolfram | A New Kind of Science[END_REF][19] [START_REF] Gorard | Some Relativistic and Gravitational Properties of the Wolfram Model[END_REF] [START_REF] Gorard | Some Quantum Mechanical Properties of the Wolfram Model[END_REF], 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 [START_REF] Benincasa | Scalar Curvature of a Causal Set[END_REF],

or the discrete Einstein-Hilbert action from Ollivier-Ricci curvature in Wolfram model systems [START_REF] Gorard | Some Relativistic and Gravitational Properties of the Wolfram Model[END_REF][23] [START_REF] Gorard | Hypergraph Discretization of the Cauchy Problem in General Relativity via Wolfram Model Evolution[END_REF]), 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 SO + (1, 3), in the limit of infinite sprinkling density [START_REF] Bombelli | Discreteness without symmetry breaking: A theorem[END_REF]), 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 partiallyordered 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 [START_REF] Gorard | Some Relativistic and Gravitational Properties of the Wolfram Model[END_REF] [START_REF] Gorard | Hypergraph Discretization of the Cauchy Problem in General Relativity via Wolfram Model Evolution[END_REF]. 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 [START_REF] Bolognesi | Causal sets from simple models of computation[END_REF][27] [START_REF] Gorard | Algorithmic Causal Sets and the Wolfram Model[END_REF] (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 [START_REF] Gorard | Axiomatic Quantum Field Theory in Discrete Spacetime via Multiway Causal Structure: The Case of Entanglement Entropies[END_REF], 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 [START_REF] Dowker | Causal set d'Alembertians for various dimensions[END_REF], Sorkin [START_REF] Sorkin | Does Locality Fail at Intermediate Length-Scales?[END_REF] and Johnston [START_REF] Johnston | Particle propagators on discrete spacetime[END_REF] in the context of causal set theory. Following an ansatz proposed by Dowker et al. [START_REF] Dowker | Scalar Field Green Functions on Causal Sets[END_REF], as well as a more direct approach outlined by Johnston [START_REF] Johnston | Feynman Propagator for a Free Scalar Field on a Causal Set[END_REF][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 Wentzel [START_REF] Wentzel | Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanic[END_REF], Kramers [START_REF] Kramers | Wellenmechanik und halbzahlige Quantisierung[END_REF] and Brillouin [START_REF] Brillouin | La mécanique ondulatoire Schrödinger: une méthode générale de resolution par approximations successives[END_REF] 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 [START_REF] Lemaître | L'Univers en expansion[END_REF][38] [START_REF] Bondi | Spherically Symmetrical Models in General Relativity[END_REF]) in the spherically-symmetric case, or a spinning, collapsing disk of dust (described by the Weyl-Lewis-Papapetrou metric [START_REF] Weyl | Zur Gravitationstheorie[END_REF][41] [START_REF] Papapetrou | A Static Solution of the Equations of the Gravitational Field for an Arbitrary Charge-Distribution[END_REF]) 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. [START_REF] Alic | Conformal and covariant formulation of the Z4 system with constraint-violation damping[END_REF] and Bona et al. [START_REF] Bona | General-covariant evolution formalism for Numerical Relativity[END_REF] 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 [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF], 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 [START_REF] Gonc | Black hole formation from massive scalar fields[END_REF] 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 [START_REF] Neugebauer | The Einsteinian Gravitational Field of the Rigidly Rotating Disk of Dust[END_REF][48] [START_REF] Neugebauer | General Relativistic Gravitational Field of a Rigidly Rotating Disk of Dust: Solution in Terms of Ultraelliptic Functions[END_REF] (based on an earlier conjecture of Bardeen and Wagoner[50] [START_REF] Bardeen | Relativistic Disks. I. Uniform Rotation[END_REF]) 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. µ, ν, ρ, σ, 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.

Governing Equations, Discretization Scheme and Numerics

The starting point for our numerical implementation is the canonical 3 + 1 decomposition of the spacetime metric tensor g µν into the standard ADM line element of Arnowitt, Deser and Misner [START_REF] Arnowitt | Dynamical Structure and Definition of Energy in General Relativity[END_REF]:

ds 2 = -α 2 dt 2 + γ ij dx i + β i dt dx j + β j dt , (1)
where the lapse function α and the shift vector β are the ADM gauge variables (defining a foliation of a globally hyperbolic spacetime into non-overlapping spacelike hypersurfaces), and γ ij denotes the induced spatial metric tensor on the spacelike hypersurfaces of constant time. Within this foliation, the futurepointing (null) unit normal vector n to the spacelike hypersurfaces is given in terms of the contravariant derivative ∇ µ of the time coordinate t:

n µ = -α∇ µ t, (2)
with the corresponding time vector t given by:

t µ = αn µ + β µ . (3)
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 γ ij :

∂ t γ ij = -2αK ij + ∇ i β j + ∇ j β i , (4)
and 6 for the extrinsic curvature tensor on spacelike hypersurfaces K ij :

∂ t K ij = α R (3) ij -2K ik K k j + K ij K -∇ i ∇ j α-8π S ij - 1 2 γ ij (S -τ ) +β k ∂ k K ij +K ik ∂ j β k +K kj ∂ i β k . ( 5
)
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 γ ij along the normal vector n:

K ij = - 1 2 L n γ ij = -γ k i γ l j ∇ (4)
k n l = -∇

j n i , (3)
K denotes its trace (i.e. the mean curvature):

K = γ ij K ij , (7)
and S ij , S and τ are simple functions of the stress-energy tensor T µν :

τ = T µν n µ n ν , S ij = T µν γ µ i γ ν j , S = γ ij S ij . (8)
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:

H = R (3) -K ij K ij + K 2 = 0, (9)
and:

M i = γ jl ∂ l K ij -∂ i K jl -Γ m jl K mi + Γ m ji K ml = 0, (10)
respectively. Here, as elsewhere, the bracketed integers are used to distinguish between spatial vs. spacetime variants of expressions and operations: ∇ (3) designates a covariant derivative over a spacelike hypersurface, ∇ (4) designates a covariant derivative over the entire spacetime, R (3) 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ácek [START_REF] Bona | Symmetry-breaking mechanism for the Z4 general-covariant evolution system[END_REF], one can construct an explicitly covariant formulation of the Einstein field equations with Lagrangian density [START_REF] Bona | Action principle for Numerical Relativity evolution systems[END_REF]:

L = √ gg µν [R µν + 2∇ µ Z ν ] , (11)
for some 4-vector Z, such that ADM Hamiltonian and momentum constraints reduce to the trivial algebraic condition Z µ = 0, and the Einstein field equations themselves take the form:

R µν + ∇ µ Z ν + ∇ ν Z µ = 8π T µν - 1 2 T g µν . (12)
This is commonly known as the (fully covariant) Z4 formulation. The components of the 4-vector Z µ 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 [START_REF] Gundlach | Constraint damping in the Z4 formulation and harmonic gauge[END_REF], one can introduce damping terms into the covariant Z4 equations so as to formulate the evolution equations as a so-called "λ-system" based on the techniques of Brodbeck, Frittelli, Hübner and Reula [START_REF] Brodbeck | Einstein's Equations with Asymptotically Stable Constraint Propagation[END_REF]. More specifically, if one replaces the Z4 evolution equations with the damped form:

R µν + ∇ µ Z ν + ∇ ν Z µ -κ 1 [t µ Z ν + t ν Z µ -(1 + κ 2 ) g µν t σ Z σ ] = 8π T µν - 1 2 T g µν , (13)
where κ 1 ≥ 0 and κ 2 are real constants, and t µ 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 µ = 0 are satisfied will provably form a set of attractor solutions for the corresponding λ-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 β. The 6 evolution equations for the spatial metric tensor γ ij remain unchanged:

(∂ t -L β ) γ ij = -2αK ij , (14)
although the 6 evolution equations for the extrinsic curvature tensor K ij now inherit a non-trivial dependency on the Z 4-vector (as well as the damping terms κ 1 and κ 2 ):

(∂ t L β ) K ij = -∇ i α j + α R ij + ∇ i Z j + ∇ j Z i -2K l i K lj + (K -2Θ) K ij -κ 1 (1 + κ 2 ) Θγ ij -8πα S ij - 1 2 (S -τ ) γ ij , ( 15
)
where Θ is simply the projection of the Z 4-vector in the normal direction n:

Θ = n µ Z µ = αZ 0 . ( 16
)
This projection itself is subject to an evolution equation:

(∂ t -L β ) Θ = α 2 R + 2∇ j Z j + (K -2Θ) K -K ij K ij -2 Z j α j α -2κ 1 (2 + κ 2 ) Θ -16πτ , ( 17
)
as indeed is the overall Z 4-vector:

(∂ t -L β ) Z i = α ∇ j K j i -δ j i K + ∂ i Θ -2K j i Z j -Θ α i α -κ 1 Z i -8πS i , (18)
which together close the system. Here, the α i terms correspond to the constants α i = 0 that appear in the general definition of the λ-system for constraint quantities in the Einstein field equations:

∂λ ∂t = α 0 C -β 0 λ, ∂λ i ∂t = α 1 C i -β 1 λ i , ( 19
)
∂λ ij k ∂t = α 3 C ij k -β 3 λ ij k , ∂λ ij kl ∂t = α 4 C ij kl -β 4 λ ij kl , (20)
where the β i > 0 are also constants, C, C i , C ij k and C ij kl are scalar, vector, rank-3 tensor and rank-4 tensorvalued constraint quantities, respectively, and λ, λ i , λ ij k and λ ij kl are corresponding tensorial quantities with the same symmetries. Note also that the functions of the stress-energy tensor S ij , S i and τ 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):

S ij = T ij , S i = n ν T ν i , τ = n µ n ν T µν . (21)
The fully conformal and covariant formulation of Z4 (otherwise known as CCZ4) is thus obtained by performing a rescaling of the spatial metric tensor γ ij by a conformal factor φ[57] [START_REF] Hilditch | Compact binary evolutions with the Z4c formulation[END_REF]:

γij = φ 2 γ ij , (22)
chosen such as to guarantee unit determinant for the conformal spatial metric γij :

φ = (det (γ ij )) -1 6 (23)
The trace-free components Ãij of the total extrinsic curvature tensor K ij on spacelike hypersurfaces then become (after conformal rescaling):

Ãij = φ 2 (K ij ) T F = φ 2 K ij - 1 3 Kγ ij = φ 2 K ij - 1 3 K ij γ ij γ ij . (24)
By introducing the Christoffel symbols Γi jk associated to the conformal spatial metric tensor γ ij , as well as the pseudovector Γi derived from the contraction of this conformal connection, namely:

Γi jk = 1 2 γil (∂ j γkl + ∂ k γjl -∂ l γjk ) , (25)
and:

Γi = γjk Γi jk = γij γkl ∂ l γjk , (26)
respectively, we can further decompose the Ricci curvature tensor R

ij on spacelike hypersurfaces into a sum:

R (3) ij = R(3) ij + Rφ(3) ij ( 27
)
of a rank-2 tensor R

(3) ij depending only upon spatial derivatives of the conformal spatial metric tensor γij :

R(3) ij = - 1 2 γlm ∂ l ∂ m γij + γk(i ∂ j) Γk + Γk Γ(ij)k + γlm 2 Γk l(i Γ j)km + Γk im Γkjl , (28)
and a second rank-2 tensor Rφ(3) ij depending only upon derivatives of the conformal factor φ:

Rφ(3) ij = 1 φ 2 φ ∇i ∇j φ + γij ∇l ∇l φ -2γ ij ∇l φ ∇l φ . (29)
The fully conformal and covariant Z4 formulation now decomposes into a system of 6 evolution equations for the conformal spatial metric tensor γij (written here explicitly, without Lie derivatives):

∂ t γij = -2α Ãij + 2γ k(i ∂ j) β k - 2 3 γij ∂ k β k + β k ∂ k γij , (30)
6 evolution equations for the trace-free components of the (conformal) extrinsic curvature tensor Ãij :

∂ t Ãij = φ 2 [-∇ i ∇ j α + α (R ij + ∇ i Z j + ∇ j Z i -8πS ij )] T F + α Ãij (K -2Θ) -2α Ãil Ãl j + 2 Ãk(i ∂ j) β k - 2 3 Ãij ∂ k β k + β k ∂ k Ãij , (31)
and single evolution equations for the trace of the extrinsic curvature tensor K:

∂ t K = -∇ i ∇ i α + α R + 2∇ i Z i + K 2 -2ΘK + β j ∂ j K -3ακ 1 (1 + κ 2 ) Θ + 4πα (S -3τ ) , (32)
the conformal factor φ:

∂ t φ = 1 3 αφK - 1 3 φ∂ k β k + β k ∂ k φ, (33)
and the projection Θ of the Z 4-vector in the normal direction n:

∂ t Θ = 1 2 α R + 2∇ i Z i -Ãij Ãij + 2 3 K 2 -2ΘK -Z i ∂ i α + β k ∂ k Θ -ακ 1 (2 + κ 2 ) Θ -8πατ. ( 34
)
Finally, in order to close the system, we must introduce a set of 4 evolution equations for the components of the pseudovector Γi derived from the conformal connection coefficients Γi jk ; for simplicity, we represent these as equations for the components Γi , defined as:

Γi = Γi + 2γ ij Z j , (35)
namely:

∂ t Γi = 2α Γi jk Ãjk -3 Ãjk ∂ j φ φ - 2 3 γij ∂ j K + 2γ ki α∂ k Θ -Θ∂ k α - 2 3 αKZ k -2 Ãij ∂ j α + γkl ∂ k ∂ l β i + 1 3 γik ∂ k ∂ l β l + 2 3 Γi ∂ k β k -Γk ∂ k β i + 2κ 3 2 3 γij Z j ∂ k β k -γjk Z j ∂ k β i + β k ∂ k Γi -2ακ 1 γij Z j -16παγ ij S j . ( 36
)
Note that the new κ 3 parameter in the evolution equations for the conformal connection coefficients Γi above is currently unconstrained; although κ 3 = 1 yields a fully covariant formulation (at the cost of numerical stability in the case of black hole spacetimes) and κ 3 = 1 2 yields a non-covariant formulation which is nevertheless numerically stable, we choose instead to follow the prescription of Alic, Kastaun and Rezzolla [START_REF] Alic | Constraint damping of the conformal and covariant formulation of the Z4 system in simulations of binary neutron stars[END_REF] and eliminate the numerical instabilities by choosing κ 1 → κ1 α as the choice of damping parameter, which preserves spatial covariance in all terms depending on κ 1 , and full spacetime covariance for all other terms. Moreover, for simulations involving spherical symmetry, we employ the maximal slicing gauge condition [START_REF] Alcubierre | Appearance of coordinate shocks in hyperbolic formalisms of general relativity[END_REF][61] on the lapse function α:

∇ 2 α = αK ij K ij , (37)
and the more general 1 + log slicing gauge condition [START_REF] Alcubierre | Simple excision of a black hole in 3 + 1 numerical relativity[END_REF] for the axially-symmetric case:

∂ t α = -2α (K -2Θ) + β k ∂ k α, (38)
due to their known stability (and singularity-avoidance) properties for black hole spacetimes. By the same token, we use the gamma-driver gauge conditions [START_REF] Anninos | Collision of two black holes[END_REF][64] on the components of the shift vector β i :

∂ t β i = η 1 B i + β k ∂ k β i , where ∂ t B i = µ β1 α µ β 2 ∂ t Γi -β k ∂ k Γi + β k ∂ k B i -η 2 B i , (39)
or, in the simplified case (without any advection terms):

∂ t β i = η 1 B i , where ∂ t B i = µ β1 α µ β 2 ∂ t Γi -η 2 B i (40)
where B i is an auxiliary vector field, and η 1 , η 2 , µ β1 and µ β2 are scalar parameters, chosen here to be η 1 = 3 4 , µ β1 = 1, µ β2 = 0 and η 2 = 1 (i.e. the standard hyperbolic driver conditions). The standard conformal BSSN formalism of Baumgarte, Shapiro, Shibata and Nakamura [START_REF] Nakamura | General Relativistic Collapse to Black Holes and Gravitational Waves from Black Holes[END_REF] is recovered in the special case where Θ = 0 and Z i = 0, the trace of the extrinsic curvature tensor K is replaced with K BSSN , defined by:

K BSSN = K -2Θ, (41)
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 (γ ij ) and tr Ãij ):

det (γ ij ) = 1, and tr Ãij = 0, (42)
are enforced explicitly as part of the numerical integration algorithm at each time step (by subtracting out the trace from Ãij and rescaling the conformal spatial metric tensor γ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 Φ obeying the (massive) Klein-Gordon equation:

+ m 2 Φ = 0, (43)
or for a more general potential V (Φ) (with the massive case hence corresponding to V (Φ) = m 2 ΦΦ):

Φ + dV dΦ = 0, (44)
which, in curved spacetime, becomes:

g µν ∇ µ ∇ ν Φ = dV dΦ . (45)
The solution to the Klein-Gordon equation yields a stress-energy tensor T µν with components:

T µν = 2∇ µ Φ∇ ν Φ -g µν ∇ ρ Φ∇ ρ Φ -m 2 ΦΦ , (46)
in the massive case, or, in lowered-index form T µν :

T µν = ∇ µ Φ∇ ν Φ - 1 2 g µν (∇ ρ Φ∇ ρ Φ + 2V ) , (47)
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 Φ is then given by:

∂ t Φ = αΠ m + β i ∂ i Φ, ( 48
)
where Π m is the (negative of the) conjugate momentum of the scalar field given by:

Π m = 1 α ∂ t Φ -β i ∂ i Φ . (49)
However, we can immediately see that this "definition" of the (negative) conjugate momentum Π m is really just a trivial restatement of the "evolution equation" for the scalar field Φ; therefore, the "true" evolution equation for the scalar field must be given in terms of an evolution equation for its (negative) conjugate momentum Π m directly, instead, i.e:

∂ t Π m = β i ∂ i Π m + γ ij (α∂ j ∂ i Φ + ∂ j Φ∂ i α) + α KΠ m -γ ij Γ k ij ∂ k Φ + dV dΦ . (50)
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 [START_REF] Gorard | Hypergraph Discretization of the Cauchy Problem in General Relativity via Wolfram Model Evolution[END_REF], this can be achieved by applying a generalized bicubic:

f (x, y) = 3 i=0 3 j=0 a ij x i y j , (51)
or tricubic [START_REF] Lekien | Tricubic interpolation in three dimensions[END_REF]:

f (x, y, z) = 3 i=0 3 j=0 3 k=0 a ijk x i y j z k , (52)
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:

dΦ (x, t) dt = F (t, Φ (x, t)) , where Φ (x, t 0 ) = Φ 0 , (53)
where Φ (x, t) is a generic state vector, i.e:

Φ (x, t) = (φ 1 (x, t) , φ 2 (x, t) , φ 3 (x, t) , . . . ) , (54)
we can evolve the system forwards in units of a discrete time step ∆t (chosen so as to be compatible with the Courant condition, and therefore numerically stable) as [START_REF] Kreiss | On the stability definition of difference approximations for the initial boundary value problem[END_REF][68]:

Φ (x, t n+1 ) = Φ (x, t n ) + 1 6 ∆t (k 1 + 2k 2 + 2k 3 + k 4 ) , where t n+1 = t n + ∆t. ( 55
)
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:

∂Φ (x, t) ∂t = F (Φ (x, t)) , ( 56
)
where F is a non-linear operator acting on Φ, 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 x and a global time coordinate t for the hypergraphs, which can be obtained using the geometrical construction described above.

Figure 1: Interpolating a value for the endpoint of the red geodesic in a spatial hypergraph with a twodimensional Riemannian manifold-like limiting structure, as generated by the hypergraph rewriting rule {{x, x, y} , {x, z, w}} → {{w, w, v} , {v, w, y} , {z, y, v}}, using a generalized bicubic interpolation algorithm applied to the 7 vertices contained within the blue subhypergraph.

Figure 2: Computing an inner product of two discrete geodesics (shown in red and blue) in a spatial hypergraph with a two-dimensional Riemannian manifold-like limiting structure, as generated by the hypergraph rewriting rule {{x, x, y} , {x, z, w}} → {{w, w, v} , {v, w, y} , {z, y, v}}, by discrete projection/"dropping a perpendicular" (shown in green), yielding a normalized inner product value of 3 4 .

In the explicit time evolution formula above, we choose:

k 1 = F (t n , φ (x, t n )) , k 4 = F (t n + ∆t, Φ (x, t n ) + ∆tk 3 ) , (57)
and:

k 2 = F t n + ∆t 2 , Φ (x, t n ) + ∆t k 1 2 , k 3 = F t n + ∆t 2 , Φ (x, t n ) + ∆t k 2 2 , (58)
which completes the temporal discretization of the scheme. For spatial discretization, we choose a coordinate structure such that ∆x = x i+1 -x i , ∆y = y j+1 -y j and ∆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 Φ n i,j,k 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 [START_REF] Zlochower | Accurate black hole evolutions by fourth-order numerical relativity[END_REF], with first derivatives (in x) defined by:

∂ x F t n , Φ n i,j,k = 1 12∆x F t n , Φ n i-2,j,k -8F t n , Φ n i-1,j,k + 8F t n , Φ n i+1,j,k -F t n , Φ n i+2,j,k , (59)
and likewise for all other first derivatives. Note that, whenever advection terms (i.e. terms of the form

β i ∂ i F t n , Φ n i,j,k
) are present, one must instead use a fourth-order upwind finite-difference scheme:

∂ x F t n , Φ n i,j,k = 1 12∆x -F t n , Φ n i-3,j,k + 6F t n , Φ n i-2,j,k -18F t n , Φ n i-1,j,k +10F t n , Φ n i,j,k + 3F t n , Φ n i+1,j,k , (60)
whenever β x < 0, and:

∂ x F t n , Φ n i,j,k = 1 12∆x F t n , Φ n i+3,j,k -6F t n , Φ n i+2,j,k + 18F t n , Φ n i+1,j,k -10F t n , Φ n i,j,k -3F t n , Φ n i-1,j,k , (61)
whenever β 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:

∂ xx F t n , Φ n i,j,k = 1 12 (∆x) 2 -F t n , Φ n i+2,j,k + 16F t n , Φ n i+1,j,k -30F t n , Φ n i,j,k +16F t n , Φ n i-1,j,k -F t n , Φ n i-2,j,k , (62)
and:

∂ xy F t n , Φ n i,j,k = 1 144∆x∆y F t n , Φ n i-2,j-2,k -8F t n , Φ n i-1,j-2,k + 8F t n , Φ n i+1,j-2,k -F t n , Φ n i+2,j-2,k -8 F t n , Φ n i-2,j-1,k -8F t n , Φ n i-1,j-1,k + 8F t n , Φ n i+1,j-1,k -F t n , Φ n i+2,j-1,k + 8 F t n , Φ n i-2,j+1,k -8F t n , Φ n i-1,j+1,k + 8F t n , Φ n i+1,j+1,k -F t n , Φ n i+2,j+1,k -F t n , Φ n i-2,j+2,k -8F t n , Φ n i-1,j+2,k + 8F t n , Φ n i+1,j+2,k -F t n , Φ n i+2,j+2,k , (63)
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 [START_REF] Kreiss | Comparison of accurate methods for the integration of hyperbolic equations[END_REF] 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:

∂ t Φ i,j,k → ∂ t Φ i,j,k + σ 64∆x (Φ i+3,j,k -6Φ i+2,j,k +15Φ i+1,j,k -20Φ i,j,k + 15Φ i-1,j,k -6Φ i-2,j,k + Φ i-3,j,k ) , (64)
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 [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF], based on previous methods developed by Berger and Oliger[72] and Gropp [START_REF] Gropp | A Test of Moving Mesh Refinement for 2-D Scalar Hyperbolic Problems[END_REF]. Broadly speaking, we introduce a tagging function f (x, y, z) based on whether the L 2 norm of the change in some scalar field φ (which, for the simulations presented within this article, will always be chosen to be a Riemannian scalar curvature invariant) exceeds a pre-defined threshold σ (φ) across a given vertex or subhypergraph:

f (x, y, z) =          1, if 3 i=1 (φ (x + ∆xx i ) -φ (x -∆xx i )) 2 > σ (φ) , 0, otherwise, (65)
such that f (x, y, z) = 1 signifies that a given vertex or subhypergraph is tagged for refinement, and f (x, y, z) = 0 otherwise. The entire hypergraph is then partitioned into subhypergraphs, the signatures X (x), Y (y) and Z (z) are computed for each subhypergraph:

X (x) = f (x, y, z) dydz, Y (y) = f (x, y, z) dxdx, Z (z) = f (x, y, z) dxdy, (66)
and the Laplacians ∂ 2 x X (x), ∂ 2 y Y (y) and ∂ 2 z Z (z) 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 ∆ ∂ 2 i X i (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 (x i ) = 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 subhy-pergraph, 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:

Φ n+1 i,j,k = Φ n i,j,k - ∆t ∆x F Φ n i+ 1 2 ,j,k -F Φ n i-1 2 ,j,k - ∆t ∆y G Φ n i,j+ 1 2 ,k -G Φ n i,j-1 2 ,k - ∆t ∆x H Φ n i,j,k+ 1 2 -H Φ n i,j,k-1 2 , (67)
with F , G and H denoting the inter-vertex flux functions projected in the x, y and z directions, respectively, and with Φ i± 1 2 ,j,k , Φ i,j± 1 2 ,k and Φ i,j,k± 1 2 denoting the extrapolated values of the state vector Φ 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 Φ coarse (i.e. the state vector at refinement level l -1) are defined in terms of (conservative) averages of the values Φ f ine in the fine subhypergraph (i.e. the state vectors at refinement level l):

1 r 3 r-1 p=0 r-1 q=0 r-1 s=0 Φ f ine m+p,n+q,o+s → Φ coarse i,j,k , (68)
where i, j and k are the discrete coordinates of the coarse vertex, where the same vertex spans the discrete

intervals [m, m + r -1], [n, n + r -1] and [o, o + r -1]
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):

Φ i,j,k (t + ∆t coarse ) = Φ i,j,k - ∆t coarse ∆x F Φ i+ 1 2 ,j,k (t) - 1 r 3 r-1 q=0 r-1 p=0 r-1 s=0 F Φ m+ 1 2 ,n+p,o+s (t + q∆t f ine ) - ∆t coarse ∆y G Φ i,j+ 1 2 ,k (t) -G Φ i,j-1 2 ,k (t) - ∆t coarse ∆z H Φ i,j,k+ 1 2 (t) -H Φ i,j,k-1 2 (t) , (69)
with ∆t coarse and ∆t f ine 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:

∆t l ∆x l = ∆t l-1 ∆x l-1 = • • • = ∆t 1 ∆x 1 , (70)
and with ∆x, ∆y and ∆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 Φ coarse i,j,k (t + ∆t coarse ) and replaced with the corresponding fine ones. One starts by constructing a tensor δF of fluxes through all coarse subhypergraph boundaries that are also outer boundaries of fine subhypergraphs, i.e. in the x direction:

δF old Φ i+ 1 2 ,j,k = -F coarse Φ i+ 1 2 ,j,k , (71)
and likewise for the other coordinate directions. After each stable time step ∆t f ine in the fine subhypergraph, we then add a sum of all the fine subhypergraph fluxes through the boundary at discrete coordinates i + 1 2 , j, k :

δF new Φ i+ 1 2 ,j,k = δF old Φ i+ 1 2 ,j,k + 1 r 2 r-1 p=0 r-1 q=0 F Φ m+ 1 2 ,n+p,o+s , (72)
and likewise for the other discrete coordinate projections i, j + 1 2 , k and i, j, k + 1 2 . After r such time steps ∆t f ine have elapsed, the flux tensor δF new Φ i+ 1 coarse subhypergraph so as to match those computed using the modified time evolution formula presented above; for instance, for a vertex with discrete coordinates (i + 1, j, k) only one correction is needed:

Φ coarse(new) i+1,j,k = Φ coarse(old) i+1,j,k + ∆t coarse ∆x coarse δF new Φ i+ 1 2 ,j,k , (73)
whereas two corrections would be needed for the vertex at discrete coordinates (i + 2, j, k):

Φ coarse(new) i+2,j,k = Φ coarse(old) i+2,j,k + ∆t coarse ∆x coarse δF new Φ i+ 1 2 ,j,k - ∆t coarse ∆x coarse δF new Φ i+ 3 2 ,j,k , (74)
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 Φ, i.e. Φ i± 1 2 ,j,k , Φ i,j± 1 2 ,k and Φ i,j,k± 1 2 (and hence also to compute the inter-vertex flux functions

F Φ i± 1 2 ,j,k , G Φ i,j± 1 2 ,k and H Φ i,j,k± 1 2 ),
we choose to use a weighted, essentially non-oscillatory (WENO) spatial reconstruction scheme [START_REF] Jiang | Efficient Implementation of Weighted ENO Schemes[END_REF][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 {Ψ l } l=1 M + 1; the Laguerre polynomials j (x) constitute a natural such choice:

j (x) = 1≤m≤M +1,m =j x -x m x j -x m = (x -x 1 ) (x j -x 1 ) • • • (x -x j-1 ) (x j -x j-1 ) (x -x j+1 ) (x j -x j+1 ) • • • (x -x M +1 ) (x j -x M +1 ) , (75)
where 1 ≤ j ≤ M + 1, interpolating between the set of M + 1 nodal points (x j , y j , z j ) from (x 1 , y 1 , z 1 ) to

(x M +1 , y M +1 , z M +1
) 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 [0, 1] using the coordinate transformation (x, y, z) → (ξ, η, ζ):

ξ (x, i) = 1 ∆x i x -x i-1 2 , η (y, j) = 1 ∆y j y -y j-1 2 , ζ (z, k) = 1 ∆z k z -z k-1 2 . ( 76
)
We denote the M + 1 nodal points (x j , y j , z j ) using the shorthand

{x k } M +1
k=1 , such that:

ψ l (x k ) = δ lk , where l, k = 1, 2, . . . , M + 1. ( 77
)
Using the numerical stencils S i,j,k , defined for each of the Cartesian coordinate directions:

S s,x i,j,k = i+R e=i-L I e,j,k , S s,y i,j,k = j+R e=j-L I i,e,k , S s,z i,j,k = k+R e=k-L I i,j,e , (78)
where I i,j,k denotes the vertex with discrete coordinates (i, j, k), 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 fourthorder 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 = M 2 + 1, R = M 2 and s = 1, L = M 2 , R = M 2 + 1
, respectively), one left-sided (with s = 2, L = M , R = 0) and one rightsided (with s = 3, L = 0, R = M ). The spatial reconstruction procedure itself works by taking the secondorder boundary extrapolation step from the SLIC/MUSCL-Hancock finite-volume approaches [START_REF] Sweby | Godunov Methods[END_REF][77], namely (assuming extrapolation in the x coordinate direction):

Φ L i+ 1 1 2 φ + i-1 2 ,j,k Φ n i,j,k -Φ n i-1,j,k , (79)
for left boundary extrapolation and:

Φ R i+ 1 1 2 φ - i+ 3
the nodal basis polynomials Ψ l (ξ) as:

w s,x h (x, t n ) = M p=0 Ψ p (ξ) ŵn,s i,j,k,p = Ψ p (ξ) ŵn,s i,j,k,p , (81)
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 S s,x i,j,k , we obtain the following (linear) system of integral equations:

∀I e,j,k ∈ S s,x i,j,k , 1 ∆x e x e+ 1 2 x e-1 2 Ψ p (ξ (x)) ŵn,s i,j,k,p dx = Φn e,j,k , (82)
where Φn i,j,k denotes the (spatial) average of the solution vector Φ, when integrated over the vertex I i,j,k at time t n . The coefficients ŵn,s i,j,k,p 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:

w x h = Ψ p (ξ) ŵn i,j,k,p , where ŵn i,j,k,p = Ns s=1 ω s ŵn,s i,j,k,p , (83)
where we set N s = 4 since the polynomial degree M is odd, and the data-dependent non-linear weights ω s are defined by:

ω s = ωs q ωq , where ωs = λ s (σ s + ) r . (84)
In the above, we choose σ s to be an oscillation indicator function:

σ s = Σ lm wn,s l wn,s m , (85)
where Σ lm is an oscillation indicator matrix [START_REF] Dumbser | Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws[END_REF]:

Σ lm = M α=1 1 0 ∂ α Φ l (ξ) ∂ξ α ∂ α Ψ m (ξ) ∂ξ α dξ, (86)
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 λ s = 1 and λ 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 [START_REF] Gorard | Hypergraph Discretization of the Cauchy Problem in General Relativity via Wolfram Model Evolution[END_REF].

3 Massive Scalar Field Collapse to a Non-Rotating Schwarzschild

Black Hole

Our approach here will be to follow the analysis of Gonçalves and Moss [START_REF] Gonc | Black hole formation from massive scalar fields[END_REF], in which the WKB approximation of Wentzel [START_REF] Wentzel | Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanic[END_REF], Kramers [START_REF] Kramers | Wellenmechanik und halbzahlige Quantisierung[END_REF] and Brillouin [START_REF] Brillouin | La mécanique ondulatoire Schrödinger: une méthode générale de resolution par approximations successives[END_REF] is applied to the case of a minimally-coupled massive scalar field in spherical symmetry, in order to show that, in the limit of a scalar field of infinite mass, the resulting spacetime geometry is described by the Lemaître-Tolman-Bondi metric [START_REF] Lemaître | L'Univers en expansion[END_REF][38] [START_REF] Bondi | Spherically Symmetrical Models in General Relativity[END_REF] for a non-rotating inhomogeneous (collapsing or expanding) dust. From this analysis, it becomes possible to deduce sufficient conditions on the field parameters for a spherically-symmetric massive scalar field "bubble collapse" problem to yield the same idealized stellar collapse solution as that analyzed by Oppenheimer and Snyder [START_REF] Oppenheimer | On Continued Gravitational Contraction[END_REF]. These conditions may be derived analytically for the toy case of an idealized "top hat" initial density distribution of the scalar field, but for the more physically plausible case of an exponential initial density distribution, a numerical approach must be adopted instead. To a great extent, the analysis is simplified by the fact that Birkhoff's theorem guarantees that the exterior spacetime geometry for a spherically-symmetric ball of collapsing dust and the exterior spacetime geometry for the resulting uncharged, non-rotating black hole are identical: they are both described by the Schwarzschild metric. This allows us to neglect any considerations regarding the existence of a smooth transition between the relevant geometries (in contrast to the axially-symmetric case).

We begin by considering the line element for the most general possible spherically-symmetric metric, i.e.

the line element for a spacetime whose isometry group contains a subgroup of the special orthogonal group SO (3), such that the orbits of this group are all 2-spheres, namely:

ds 2 = -e 2Ψ(τ,r) dτ 2 + e -2Λ(τ,r) dr 2 + R 2 (τ, r) dΩ 2 , ( 87
)
where r is the ordinary radial coordinate, R (τ, r) is the circumferential radial coordinate (such that the proper circumference at radius R (τ, r) is always 2πR (τ, r)), τ is the proper time, Ψ (τ, r) and Λ (τ, r) are arbitrary functions, and dΩ 2 designates the induced metric on the 2-sphere, i.e:

dΩ 2 = dθ 2 + sin (θ) dφ 2 , (88)
for colatitude coordinate θ and longitude coordinate φ. If we rescale the proper time coordinate τ to correspond instead to the proper time t for an observer that is comoving with the radial coordinate r:

t = τ 0 e Ψ(τ ,r) dτ , (89)
then we can rewrite the metric line element in Gaussian polar coordinates as:

ds 2 = -dt 2 + e -2Λ(t,r) dr 2 + R 2 (t, r) dΩ 2 . ( 90
)
The non-zero components of the Einstein curvature tensor G µν for such a spacetime are then given, up to redundancies due to symmetry, by the time-time component G tt :

G tt = 1 R 2 (t, r) -R (t, r) e 2Λ(t,r) 2∂ r R (t, r) ∂ r Λ (t, r) + 2∂ rr R (t, r) + 1 R (t, r) (∂ r R (t, r)) 2 -2∂ t R (t, r) ∂ t Λ (t, r) R (t, r) + 1 + (∂ t R (t, r)) 2 , (91)
the radial-time component G rt :

G rt = - 2 R (t, r) [∂ rt R (t, r) + ∂ r R (t, r) ∂ t Λ (t, r)] , (92)
the radial-radial component G rr :

G rr = - 1 R 2 (t, r) e -2Λ(t,r) 2∂ tt R (t, r) R (t, r) + (∂ t R (t, r)) 2 + 1 -(∂ r R (t, r)) 2 , (93)
and the colatitude-colatitude component G θθ (or, equivalently, the longitude-longitude component G φφ ):

G θθ = G φφ sin 2 (θ) = R (t, r) ∂ t R (t, r) ∂ t Λ (t, r) + ∂ r R (t, r) ∂ r Λ (t, r) e 2Λ(t,r) + ∂ rr R (t, r) e 2Λ(t,r) -∂ tt R (t, r) + ∂ tt Λ (t, r) R (t, r) -(∂ t Λ (t, r)) 2 R (t, r) . ( 94
)
As stated previously, we take as our matter content a minimally-coupled, real-valued massive scalar field Φ (t, r) with mass µ, defined by the equation of motion:

-∇ µ ∇ µ + µ 2 Φ (t, r) = -+ µ 2 Φ (t, r) = 0, (95)
where designates the covariant wave operator:

= ∇ µ ∇ µ . ( 96
)
When written out explicitly in terms of our spherically-symmetric metric in Gaussian polar coordinates, this equation of motion becomes:

∂ tt Φ (t, r) -e 2Λ(t,r) ∂ rr Φ (t, r) + µ 2 Φ (t, r) + ∂ t Λ (t, r) + 2 R (t, r) ∂ t R (t, r) ∂ t Φ (t, r) -e 2Λ(t,r) ∂ r Λ (t, r) - 2 R (t, r) ∂ r R (t, r) ∂ r Φ (t, r) = 0, (97)
such that the non-zero components of the stress energy tensor T µν for the scalar field Φ (t, r) become, up to redundancies due to symmetry, the time-time component T tt :

T tt = 1 2 (∂ t Φ (t, r)) 2 + 1 2 e 2Λ(t,r) (∂ r Φ (t, r)) 2 + 1 2 µ 2 Φ (t, r) , (98)
the radial-time component T rt :

T rt = ∂ t Φ (t, r) ∂ r Φ (t, r) , (99)
the mixed-index radial-radial component T r r :

T r r = 1 2 (∂ t Φ (t, r)) 2 + 1 2 e 2Λ(t,r) (∂ r Φ (t, r)) 2 - 1 2 µ 2 Φ (t, r) , (100)
and the colatitude-colatitude component T θθ (or, equivalently, the longitude-longitude component T φφ ):

T θθ = T φφ sin 2 (θ) = 1 2 R 2 (t, r) (∂ t Φ (t, r)) 2 -e 2Λ(t,r) (∂ r Φ (t, r)) 2 -µ 2 Φ (t, r) . (101)
Following Gonçalves and Moss [START_REF] Gonc | Black hole formation from massive scalar fields[END_REF], we now recast two of the Einstein field equations (corresponding to the radial-time component G rt and the radial-radial component G rr of the Einstein curvature tensor, respectively) in terms of purely first derivatives of two scalar functions, namely k (t, r) and m (t, r), of the form:

k (t, r) = 1 -e 2Λ(t,r) (∂ r R (t, r)) 2 , and m (t, r) = 1 2 R (t, r) (∂ t R (t, r)) 2 + k (t, r) , (102)
thus yielding:

∂ t k (t, r) = 8πR (t, r) ∂ r R (t, r) T r t , (103)
and:

∂ t m (t, r) = 4πR 2 (t, r) ∂ r R (t, r) T r t -4πR 2 (t, r) ∂ t R (t, r) T r r , (104)
respectively. The overall system of equations can now be closed by also incorporating the aforementioned second-order equation of motion for the scalar field Φ (t, r) in Gaussian polar coordinates:

∂ tt Φ (t, r) -e 2Λ(t,r) ∂ rr Φ (t, r) + µ 2 Φ (t, r) + ∂ t Λ (t, r) + 2 R (t, r) ∂ t R (t, r) ∂ t Φ (t, r) -e 2Λ(t,r) ∂ r Λ (t, r) - 2 R (t, r) ∂ r R (t, r) ∂ r Φ (t, r) = 0. ( 105
)
A third Einstein field equation (corresponding to the time-time component G tt of the Einstein curvature tensor) can further be used to derive the following constraint on the spatial derivative of m (t, r):

∂ r m (t, r) = 4πR 2 (t, r) ∂ r R (t, r) T tt -4πR 2 (t, r) ∂ t R (t, r) T rt , (106)
from which we are therefore able to reconstruct the initial value of m (0, r), given the Cauchy initial data for the spacetime defined on a spacelike hypersurface with proper time t = 0.

From here, we seek to apply the standard WKB approximation of Wentzel, Kramers and Brillouin [START_REF] Wentzel | Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanic[END_REF][36] [START_REF] Brillouin | La mécanique ondulatoire Schrödinger: une méthode générale de resolution par approximations successives[END_REF] (as conventionally used in semiclassical approximations to the Schrödinger equation in non-relativistic quantum mechanics) in order to derive a wavelike solution to the equation of motion for the scalar field Φ (t, r).

Specifically, the WKB approximation applies whenever one has an n-th order differential equation with coefficients c i (x) depending on a spatial parameter x, and in which the n-th order derivative is multiplied by a scalar parameter

d n f (x) dx n + c n-1 (x) d n-1 f (x) dx n-1 + • • • + c 1 (x) df (x) dx + c 0 (x) f (x) = 0, ( 107
)
in which case, in the limit as → 0, one has the following ansatz for f (x) in the form of an asymptotic series expansion:

f (x) = lim δ→0 exp 1 δ ∞ i=0 δ i S i (x) , (108)
with expansion terms S i (x) yet to be determined, and in which the relative asymptotic scaling of the parameters δ and is determined by the equation in question. Therefore, if we let the parameter δ correspond now to the Compton wavelength of the scalar field Φ (t, r), i.e. 1 µ (for scalar field mass µ, assuming natural units with = c = 1), and if we assume that the function Φ (t, r) is of compact support, with a finite radius λ in which the value of the field Φ (t, r) is non-vanishing, then the WKB approximation applies whenever the scalar field is sufficiently massive, such that the Compton wavelength is much less than the radius of support:

1 µ λ. (109)
This allows us to perform an asymptotic series expansion for a scalar amplitude function Ψ (t, r), such that the following wavelike solution ansatz for Φ (t, r):

Φ (t, r) = 1 µ Ψ (t, r) cos (µt) , (110)
holds, with the non-zero components of the stress-energy tensor T µν for the scalar field Φ (t, r) therefore given (via simple differentiation) in terms of the amplitude function Ψ (t, r); up to redundancies due to symmetry, these are the time-time component T tt :

T tt = 1 2 Ψ 2 (t, r) - 1 2µ Ψ (t, r) ∂ t Ψ (t, r) sin (2µt) + 1 4µ 2 (∂ t Ψ (t, r)) 2 e 4Λ(t,r) (∂ r Ψ (t, r)) 2 [1 + cos (2µt)] , (111)
the radial-time component T rt :

T rt = - 1 2µ Ψ (t, r) ∂ r Ψ (t, r) sin (2µt) + 1 2µ 2 ∂ t Ψ (t, r) ∂ r Ψ (t, r) [1 + cos (2µt)] , (112)
the mixed-index radial-radial component T r r :

T r r = - 1 2 Ψ 2 (t, r) cos (2µt) - 1 2µ Ψ (t, r) ∂ t Ψ (t, r) sin (2µt) + 1 4µ 2 (∂ t Ψ (t, r)) 2 + e 4Λ(t,r) (∂ r Ψ (t, r)) 2 [1 + cos (2µt)] , (113)
and the colatitude-colatitude component T θθ (or, equivalently, the longitude-longitude component T φφ ):

T θθ = T φφ sin 2 (θ) = 1 2 R 2 (t, r) -Ψ 2 (t, r) cos (2µt) - 1 µ Ψ (t, r) ∂ t Ψ (t, r) sin (2µt) + 1 2µ 2 (∂ t Ψ (t, r)) 2 -e 4Λ(t,r) (∂ r Ψ (t, r)) 2 [1 + cos (2µt)] . (114)
If we now write out the asymptotic expansion for the scalar amplitude function Ψ (t, r) as an explicit trigonometric series of the general form:

Ψ (t, r) = Ψ 0 (t, r) + Ψ (cos) 01
(t, r) cos (µt) + Ψ (sin) 01

(t, r) sin (µt)

+ ∞ m=1 2m n=2 1 µ m Ψ (cos) mn (t, r) cos (nµt) + Ψ (sin) mn (t, r) sin (nµt) , (115)
for undetermined coefficient functions Ψ (cos) mn (t, r) and Ψ

(sin) mn (t, r), then, up to second-order in the Compton wavelength expansion parameter 1 µ , we obtain:

Ψ (t, r) = Ψ 01 (t, r) cos (µt) + 1 µ 2 Ψ (cos) 22 (t, r) cos (2µt) + O 1 µ 3 . ( 116
)
We can also write out the explicit trigonometric series expansions of the scalar functions k (t, r) and m (t, r)

(as derived from the radial-time component G rt and the radial-radial component G rr component of the Einstein curvature tensor, respectively) in much the same way, yielding:

k (t, r) = k 0 (t,
k (t, r) = k 0 (t, r) + 1 µ 2 k (cos) 22 (t, r) cos (2µt) + O 1 µ 3 , (119)
and:

m (t, r) = m 0 (t, r) + 1 µ m (sin) 12 (t, r) sin (2µt) + 1 µ 2 m (cos) 22 (t, r) cos (2µt) + O 1 µ 3 , (120)
respectively. Due to the known definitional relationship between the functions k (t, r) and the m (t, r), and the circumferential radial coordinate R (t, r), we can therefore deduce an analogous expansion for the R (t, r) coordinate up to second order in the expansion parameter 1 µ also:

R (t, r) = R 0 (t, r) + 1 µ 2 R (cos) 22 (t, r) cos (2µt) + O 1 µ 3 . ( 121
)
At this point, we are able to substitute these explicit trigonometric expansions for the scalar functions k (t, r) and m (t, r) into the first-order differential equations derived above, namely:

∂ t k (t, r) = 8πR (t, r) ∂ r R (t, r) T r t , (122)
and:

∂ t m (t, r) = 4πR 2 (t, r) ∂ r R (t, r) T r t -4πR 2 (t, r) ∂ t R (t, r) T r r , (123)
respectively, yielding simply, up to first-order in the expansion parameter 1 µ :

∂ t k 0 (t, r) = 0, and
∂ t m 0 (t, r) = 0, (124)
respectively. Moreover, we can rearrange the definition of the scalar function m (t, r) in order to obtain a first-order differential equation for the circumferential radial coordinate R (t, r) as well:

m (t, r) = 1 2 R (t, r) (∂ t R (t, r)) 2 + k (t, r) , =⇒ ∂ t R (t, r) = ± 2m (t, r) -k (t, r) R (t, r) R (t, r) , (125)
into which we can substitute its own trigonometric expansion in much the same way, obtaining (up to first-order in the expansion parameter 1 µ ):

∂ t R 0 (t, r) = ± 2m 0 (t, r) R (t, r) -k 0 (t, r). (126)
On the other hand, if we substitute the trigonometric expansion for the coordinate R (t, r) into the constraint equation for the spatial derivative of the function m (t, r) that is used in the reconstruction of the initial value of m (0, r) from the Cauchy data, namely:

∂ r m (t, r) = 4πR 2 (t, r) ∂ r R (t, r) T tt -4πR 2 (t, r) ∂ t R (t, r) T rt , (127)
then we also find (up to first-order in the expansion parameter 1 µ ):

∂ r m 0 (t, r) = 2πR 2 0 (t, r) ∂ r R 0 (t, r) Ψ 2 0 (t, r) . (128)
Therefore, we infer from the definitions of the scalar functions k (t, r) and m (t, r) that, up to first-order in 1 µ , the overall metric line element is of the following form (in Gaussian polar coordinates):

ds 2 = -dt 2 + (∂ r R 0 (t, r)) 2 1 -k 0 (t, r) dr 2 + R 2 0 (t, r) dΩ 2 , with dΩ 2 = dθ 2 + sin (θ) dφ 2 , (129)
with the scalar field amplitude parameter Ψ 0 (t, r) given by:

Ψ 0 (t, r) = ± ∂ r m 0 (t, r) 2πR 2 0 (t, r) ∂ r R 0 (t, r) . ( 130
)
As pointed out by Gonçalves and Moss [START_REF] Gonc | Black hole formation from massive scalar fields[END_REF], this metric has exactly the form of the Lemaître-Tolman-Bondi metric [START_REF] Lemaître | L'Univers en expansion[END_REF][38] [START_REF] Bondi | Spherically Symmetrical Models in General Relativity[END_REF] for a spherically-symmetric distribution of inhomogeneous dust that is either expanding or contracting uniformly due to gravity, namely:

ds 2 = -dt 2 + (∂ r R 0 (t, r)) 2 1 + 2E (t, r) dr 2 + R 2 0 (t, r) dΩ 2 , ( 131
)
where we have introduced the parameter E (t, r) to represent the specific energy (i.e. energy per unit mass) of the dust particles located at the comoving coordinate radius r at time t:

E (t, r) = - k 0 (t, r) 2 > - 1 2 , ( 132
)
obeying the same equation of motion derived above for the scalar field amplitude Ψ 0 (t, r), namely:

∂ t R 0 (t, r) = ± 2m 0 (t, r) R 0 (t, r) + 2E (t, r). (133)
We can verify explicitly that this indeed corresponds to a dust solution, since the pressure contribution to the stress-energy tensor vanishes identically (i.e. p (t, r) = 0), leaving only an energy density contribution ρ (t, r) of the form:

ρ (t, r) = ∂ r m 0 (t, r) 4πR 2 0 (t, r) ∂ r R 0 (t, r) , (134)
with the scalar function m 0 (t, r) now playing the role of the total gravitational mass contained within the sphere of comoving coordinate radius r at time t. The equation of motion for R 0 (t, r) can now be integrated analytically in terms of the parameter η, revealing that it has exactly three solutions, corresponding to the cases of hyperbolic evolution (i.e. E (t, r) > 0):

R 0 (t, r) = m 0 (t, r) 2E (t, r) (cosh (η) -1) , where (sinh (η) -η) = (2E (t, r)) 3 2 (t -t 0 (r)) m 0 (t, r) , (135)
parabolic evolution (i.e. E (t, r) = 0):

R 0 (t, r) = 9m 0 (t, r) (t -t 0 (r)) 2 2 1 3 , (136)
and elliptic evolution (i.e. E (t, r) < 0):

R 0 (t, r) = m 0 (t, r) 2 |E (t, r)| (1 -cos (η)) , where (η -sin (η)) (2 |E (t, r)|) 3 2 (t -t 0 (r)) m 0 (t, r) , (137)
respectively, where t 0 (r) is an arbitrary scalar function designating the initial time parameter for the world lines of dust particles located at the comoving coordinate radius r.

The Schwarzschild metric within a geodesic coordinate system can then be obtained as a limiting case of the Lemaître-Tolman-Bondi metric by setting the mass parameter m 0 (t, r) to be constant. For instance, if we take the specific energy of the dust to vanish (i.e. E (t, r) = 0), then we effectively transform the Schwarschild metric in the Schwarzschild coordinate system {t, r} (with Schwarzschild radius r s = 2m 0 ):

ds 2 = -1 - r s r dt 2 + dr 2 1 -rs r + r 2 dθ 2 + sin 2 (θ) dφ 2 , (138)
to the geodesic coordinate system {τ, ρ}:

dτ = dt + r s r 1 - r s r -1
dr, and

dρ = dt + r r s 1 - r s r -1 dr, (139)
in which the geodesics with constant ρ coordinate are all timelike, parametrized by proper time τ , corresponding to the trajectories of particles in free fall that begin at infinity with zero velocity, yielding the following form of the Schwarzschild metric in the Lemaître coordinate system [START_REF] Lemaître | L'Univers en expansion[END_REF]:

ds 2 = -dτ 2 + r s r dρ 2 + r 2 dθ 2 + sin 2 (θ) dφ 2 , where r = 3 2 (ρ -τ ) 2 3 r 1 3 s . ( 140
)
On the other hand, if we set the specific energy of the dust to be of the form:

E (t, r) = - 1 2 (1 + r 2 ) , (141)
then we obtain instead a coordinate system in which the Schwarzschild radial coordinate r is transformed to the Novikov radial coordinate r nov , dependent upon the initial radial height r max of the particle on its trajectory:

r nov = r max r s -1, (142)
with the time coordinate τ still corresponding to the proper time of the falling particle, yielding the Novikov form of the Schwarzschild metric, representing the trajectories of particles in free fall that begin at an arbitrary radial distance:

ds 2 = -dτ 2 + r nov + 1 r 2 ∂r ∂r nov 2 dr 2 nov + r 2 dθ 2 + sin 2 (θ) dφ 2 , (143)
from which the coordinates r and τ can be recovered by solving the following (implicitly time-dependent) pair of equations for η:

τ = m 0 r 2 nov + 1 3 2 (η + sin (η)) , and r = m 0 r 2 nov + 1 (η + cos (η)) . (144)
Similarly, the Friedmann-Lemaître-Robertson-Walker metric [START_REF] Friedmann | Über die Möglichkeit einer Welt mit konstanter negativer Krümmung des Raumes[END_REF][4] [START_REF] Robertson | Kinematics and World-Structure[END_REF][81] for a homogeneous and isotropic universe:

dτ 2 = -dt 2 + a 2 (t) dr 2 1 -kr 2 + r 2 dθ 2 + sin 2 (θ) dφ 2 , (145)
with dimensionless scale parameter a (t) and discrete curvature parameter k, is recovered in the limiting case of the Lemaître-Tolman-Bondi metric in which the initial time parameter t 0 (r) is constant, with the hyperbolic, parabolic and elliptic evolution cases (corresponding to E (t, r) > 0, E (t, r) = 0 and E (t, r) < 0, respectively) yielding discrete curvature values of k = -1, k = 0 and k = 1, respectively.

Above, we showed that the trigonometric expansions for the scalar functions k (t, r), m (t, r) and R (t, r), namely:

k (t, r) = k 0 (t, r) + 1 µ 2 k (cos) 22 (t, r) cos (2µt) + O 1 µ 3 , (146)
m (t, r) = m 0 (t, r) + 1 µ m (sin) 12 (t, r) sin (2µt) + 1 µ 2 m (cos) 22 (t, r) cos (2µt) + O 1 µ 3 , (147)
and:

R (t, r) = R 0 (t, r) + 1 µ 2 R (cos) 22 (t, r) cos (2µt) + O 1 µ 3 , (148)
respectively, could be applied to yield, up to first-order in the expansion parameter 1 µ , the following conditions on the partial derivatives ∂ t k 0 (t, r), ∂ t m 0 (t, r), ∂ t R 0 (t, r) and ∂ r m 0 (t, r):

∂ t k 0 (t, r) = 0, ∂ t m 0 (t, r) = 0, ∂ t R 0 (t, r) = ± 2m 0 (t, r) R (t, r) -k 0 (t, r), (149)
and:

∂ r m 0 (t, r) = 2πR 2 0 (t, r) ∂ r R 0 (t, r) Ψ 2 0 (t, r) , (150)
respectively. Moving now up to second-order in the expansion parameter 1 µ , we obtain the following subleading corrections to the functions k (t, r), m (t, r) and R (t, r):

k 22 (t, r) = 2π R 0 (t, r) ∂ r R 0 (t, r) (1 -k 0 (t, r)) Ψ 0 (t, r) ∂ r Ψ (t, r) , (151)
m 12 (t, r) = πR 2 0 (t, r) ∂ t R 0 (t, r) Ψ 2 0 (t, r) , (152)
and:

R 22 (t, r) = - 1 2 R 0 (t, r) Ψ 2 0 (t, r) , (153)
respectively. Without loss of generality, we shall assume the validity of the WKB approximation whenever the magnitudes of the sub-leading corrections k 22 (t, r), m 12 (t, r) and R 22 (t, r) are no greater than half of the magnitudes of the leading-order terms k 0 (t, r), m 0 (t, r) and R 0 (t, r) (the multiplicative factor of 1 2 can be modified by simply changing the value of the scalar field mass parameter µ). Since the single largest correction comes from the m 12 (t, r) contribution to m (t, r), we obtain the following region in the (t, r)-plane:

m 12 (t, r) µm 0 (t, r) ≤ 1 2 , ( 154
)
bounded by the curve:

∂ t R 0 (t, r) ∂ r m 0 (t, r) = µm 0 (t, r) ∂ r R 0 (t, r) , (155)
within which the WKB approximation holds. Thus, we henceforth restrict ourselves to points in our spacetime for which the entirety of the (interior) past light cone lies within this region.

Returning now to the equation of motion for the circumferential radial coordinate R 0 (t, r) (up to first order in the expansion parameter 1 µ ), expressed in terms of the scalar function k 0 (t, r):

∂ t R 0 (t, r) = ± 2m 0 (t, r) R 0 (t, r) -k 0 (t, r), (156)
we use the parametric analytic integration of R 0 (t, r) presented previously to deduce that:

t (η, r) = t 0 (r) + m 0 (t, r) k 3 2 0 (t, r) (η + sin (η)) , and R 0 (η, r) = 2m 0 (t, r) k 0 (t, r) cos 2 η 2 , (157)
subject to the hypothesis that ∂ r R 0 (t, r) = 0, thus ensuring that the constant-time "shells" of the dust solution do not intersect, and therefore that the evolution remains globally hyperbolic. If the Cauchy initial data are time-symmetric and of the general form: R 0 (0, r) = r, and

∂ t R 0 (t, r) = 0, ( 158
)
then this fixes the choice of radial coordinate r, as well as the choice of initial time parameter for world lines t 0 (r) and the scalar function k 0 (t, r): t 0 (r) = 0, and

k 0 (t, r) = 2m 0 (t, r) r , (159)
such that one has:

t (η, r) = r 3 8m 0 (t, r) (η + sin (η)) , and R 0 (η, r) = r cos 2 η 2 , (160)
where η ∈ [0, π]. Therefore, assuming an initial density profile ρ (0, r):

ρ (0, r) = T tt (0, r) , (161)
then the solution m 0 (t, r) can be reconstructed (for any time t) by means of the following integral:

m 0 (t, r) = 4π r 0 ρ (0, r ) (r ) 2 dr = 4π r 0 T tt (0, r ) (r ) 2 dr . (162)
Recalling that m 0 (t, r) designates the gravitational mass contained inside the comoving sphere of coordinate radius r at time t, if this gravitational mass parameter converges to a finite value in the limit as r → ∞, then this value consequently corresponds to the ADM mass M of the spacetime:

lim r→∞ [m 0 (t, r)] = M. (163)
For the purposes of the numerical tests presented within this article, we will generally employ an exponential initial density profile of the form:

ρ (0, r) = T tt (0, r) = ρ 0 exp - r λ 3 , ( 164
)
whose solution is therefore given by (after evaluating the integral for m 0 (t, r)):

m 0 (t, r) = M 1 -exp - r λ 3 , (165)
where ρ 0 denotes an initial density constant, and λ denotes (as previously) the radius of support, i.e. the radius within which the scalar field Φ (t, r) does not vanish.

We can determine the radial position of the resulting event horizon within this Lemaître-Tolman-Bondi solution by treating the parameter η as an explicit function of the radial coordinate, i.e. η (r), yielding the following equation for the radial lightlike geodesics:

dη (r) ∂r = 1 R 0 (t, r) ± r 2m -1 -1 2 ∂ r R 0 (t, r) - 1 2 γ (r) (η (r) + sin (η (r))) , (166)
where we have introduced a function γ (r) whose form depends upon the initial mass distribution m 0 (t, r):

γ (r) = 3 2 - r∂ r m 0 (t, r) 2m 0 (t, r) . (167)
Since the resulting black hole should have a Schwarzschild radius equal to twice the ADM mass r s = 2M , we impose the following boundary condition at spatial infinity:

lim r→∞ [R 0 (t, r)] = 2M ; (168)
moreover, we see that the mass distribution collapses to form a spacelike singularity at proper time t = π √ 8M . This solution allows us to compute the time and space derivatives of the circumferential radial coordinate R 0 (t, r) analytically as:

∂ t R 0 (t, r) = - 2m 0 (t, r) r tan η (r) 2 , (169)
and:

∂ r R 0 (t, r) = cos 2 η (r) 2 + 1 2 γ (r) (η (r) + sin (η (r))) tan η (r) 2 , (170)
respectively. We require that the initial mass distribution should satisfy γ (r) ≥ 0, since otherwise there exists the possibility that ∂ r R 0 (t, r) = 0, causing the constant-time "shells" of the dust solution to intersect, and thereby also causing the equation for the scalar field amplitude parameter Ψ 0 (t, r):

Ψ 0 (t, r) = ± ∂ r m 0 (t, r) 2πR 2 0 (t, r) ∂ r R 0 (t, r) , (171)
to diverge, due to a failure of global hyperbolicity.

It is instructive at this point to compare this solution against the case of a massive scalar field Φ (t, r)

with a constant initial density profile ρ (0, r) = T tt (0, r), as considered by Oppenheimer and Snyder [START_REF] Oppenheimer | On Continued Gravitational Contraction[END_REF], whose solution is therefore given by (after evaluating the integral for m 0 (t, r)):

m 0 (t, r) =        M r λ 3 , if 0 < r < λ, M, if r ≥ λ. (172)
This "top-hat" form of the potential causes the WKB approximation for the amplitude of the scalar field Ψ (t, r) to break down along the boundary r = λ due to the presence of the discontinuity in the solution, although this can be rectified by making the amplitude piecewise-linear near the boundary:

Ψ (t, r) =                Ψ 0 (t, r) , if r ≤ λ -1 µ , Ψ 1 (t, r) , if λ -1 µ < r < λ, 0, if r ≥ λ, (173)
such that Ψ 0 (t, r) is constant and Ψ 1 (t, r) is no steeper than linear (in r):

Ψ 0 (t, r) = O (1) , and ∂ r Ψ 1 (t, r) = O (1) . (174)
Recall that we previously ascertained the form of the curve bounding the region within which the WKB approximation remained valid, namely:

∂ t R 0 (t, r) ∂ r m 0 (t, r) = µm 0 (t, r) ∂ r R 0 (t, r) . (175)
Hence, if we now parametrize this boundary curve using the scalar function η * (r), then we can use the explicitly-computed time and space derivatives of the circumferential radial coordinate derived above, namely:

∂ t R 0 (t, r) = - 2m 0 (t, r) r tan η * (r) 2 , (176)
and:

∂ r R 0 (t, r) = cos 2 η * (r) 2 + 1 2 γ (r) (η * (r) + sin (η * (r))) tan η * (r) 2 , (177)
to obtain the following simple constraint on the parameter η * (r):

tan η * (r) 2 = µ 3 λ 3 2M cos 2 η * (r) 2 . (178)
From this, we can conclude that the parameter η * (r) must be constant, and hence the proper time coordinate t (η * , r) must also be correspondingly constant, inside the region of support for the scalar field Φ (t, r) (i.e. for all r < λ). We also know that the Lemaître-Tolman-Bondi solution holds outside the region of support (i.e. for all r > λ), and consequently the WKB approximation must break down somewhere along a curve with a fixed proper time coordinate t * and r ∈ [0, λ]; thus, if we wish (as indeed we do) to guarantee that the WKB approximation holds everywhere outside the event horizon of the resulting black hole, it is necessary that this curve should lie entirely within the Lemaître-Tolman-Bondi event horizon.

Outside the region of support for the scalar field Φ (t, r) (i.e. for r ≥ λ), we know that the event horizon of the Lemaître-Tolman-Bondi solution lies along the circumferential radial line R = 2M , and therefore, if the event horizon is described by the curve η EH (r), then we have:

∀r ≥ λ, cos 2 η EH (r) 2 = 2M r . ( 179
)
The endpoint of the curve along which the WKB approximation fails to hold will lie on the interior of the region bounded by η EH (r) whenever:

η * (r) = η EH (λ) ⇐⇒ tan η * (r) 2 > λ 2M -1. ( 180
)
From the aforementioned constraint on the parameter η * (r), namely:

tan η * (r) 2 = µ 3 λ 3 2M cos 2 η * (r) 2 , (181)
we can moreover deduce that:

tan η * (r) 2 1 + tan 2 η * (r) 2 = µ 3 λ 3 2M , (182)
yielding the desired breakdown of the WKB approximation on the interior of the Lemaître-Tolman-Bondi event horizon η EH (r) whenever the following inequality is satisfied:

µM > 3 2 1 - 2M λ , such that lim λ→∞ [µM ] > 3 2 . ( 183
)
Thus, any constant configuration of a massive scalar field Φ (t, r) defined within a region r < λ (with λ denoting the radius of support of the field), with total ADM mass M , must collapse to form a black hole whenever this inequality is satisfied.

In order to test this hypothesis numerically for the case of the exponential initial density profile described above, we place the outermost boundary of the computational domain at radius 60M , and enforce the Sommerfeld (radiative) boundary condition:

∂f (x i , t) ∂t = - vx i r ∂f (x i , t) ∂x i -v f (x i , t) -f 0 (x i , t) r , (184)
for scalar fields f , with r = x We confirm that the ADM mass of the overall spacetime (computed by integrating over a surface surrounding the boundary of asymptotic flatness) remains approximately constant, and that the linear and angular momenta of the resulting black hole converge to be approximately zero, as expected. The L 1 , L 2 and L ∞ -norms of the (non-vacuum) Hamiltonian constraint H:

H = R -K ij K ij + K 2 -16πρ, (185)
are computed using:

H 1 = i m i |H i | , H 2 = i m i |H 2 i |, H ∞ = max i [m i |H|] , (186)
respectively, with m i designating the fraction of the overall volume of the computational domain being occupied by the i-th partitioned subhypergraph:

Figure 9: Spatial hypergraphs corresponding to the final hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted using the local curvature in the Schwarzschild conformal factor ψ, and colored according to the value of the scalar field Φ (t, r).

Figure 10: Spatial hypergraphs corresponding to the initial hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 0M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted and colored using the local curvature in the Schwarzschild conformal factor ψ.

m i = V i V tot , (187)
and likewise for the (non-vacuum) momentum constraints M i

M i = γ jl ∂ l K ij ∂ i K jl -Γ m jl K mi + Γ m ji K ml -8πS i . ( 188
)
We also approximate the location of the event horizon of the resulting black hole in the moving puncture gauge with lapse α = 0.3, confirming that this enables us to excise the spacetime region with lapse α < 0.3 by setting H to be identically zero (although we stress that, due to the aforementioned singularity-avoidance properties of the maximal-slicing gauge condition [START_REF] Alcubierre | Appearance of coordinate shocks in hyperbolic formalisms of general relativity[END_REF][61], such an excision is strictly unnecessary). The convergence rates for the Hamiltonian constraint after time t = 4.5M , with respect to the L 1 , L 2 and L ∞norms, illustrating approximately fourth-order convergence of the finite-difference scheme, are shown in Table 1. Figure 15: Spatial hypergraphs corresponding to projections along the z-axis of the initial hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 0M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted using the local curvature in ψ, and colored according to the value of the scalar field Φ (t, r).

Vertices (L 1 ) (L 2 ) (L ∞ ) O (L 1 ) O (L 2 ) O (L ∞ ) 100
Figure 16: Spatial hypergraphs corresponding to projections along the z-axis of the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 1.5M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted using the local curvature in ψ, and colored according to the value of the scalar field Φ (t, r).

Figure 17: Spatial hypergraphs corresponding to projections along the z-axis of the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 3M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted using the local curvature in ψ, and colored according to the value of the scalar field Φ (t, r).

Figure 18: Spatial hypergraphs corresponding to projections along the z-axis of the final hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted using the local curvature in ψ, and colored according to the value of the scalar field Φ (t, r).

Figure 19: Spatial hypergraphs corresponding to projections along the z-axis of the initial hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 0M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted and colored using the local curvature in ψ.

Figure 20: Spatial hypergraphs corresponding to projections along the z-axis of the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 1.5M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted and colored using the local curvature in ψ.

Figure 21: Spatial hypergraphs corresponding to projections along the z-axis of the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 3M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted and colored using the local curvature in ψ.

Figure 22: Spatial hypergraphs corresponding to projections along the z-axis of the final hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted and colored using the local curvature in ψ.

Massive Scalar Field Collapse to a Maximally-Rotating (Extremal) Kerr Black Hole

With the case of a spherically-symmetric (non-rotating) scalar field distribution thus considered, we proceed now to consider the case of an axially-symmetric (rotating) scalar field distribution of the same general form. As mentioned previously, the primary complicating factor with the more general axially-symmetric case is that, unlike in the spherically-symmetric case (in which, by Birkhoff's theorem, all of the relevant geometries are asymptotically Schwarzschild), the exterior spacetime geometry for an axially-symmetric disk of uniformly-rotating, collapsing dust and the exterior spacetime geometry for the resulting maximallyrotating (extremal) Kerr black hole are distinct, being given by the Weyl-Lewis-Papepetrou metric [START_REF] Weyl | Zur Gravitationstheorie[END_REF][41] [START_REF] Papapetrou | A Static Solution of the Equations of the Gravitational Field for an Arbitrary Charge-Distribution[END_REF] and the Kerr metric, respectively. One is therefore forced to construct some kind of smooth transition between the two limiting geometries, at least if one wishes to obtain an analytical description of the collapse. All other aspects of the analysis of Gonçalves and Moss [START_REF] Gonc | Black hole formation from massive scalar fields[END_REF], including, in particular, the application of the WKB approximation:

Φ (t, r) = 1 µ Ψ (t, r) cos (µt) , (189)
in the limit of vanishing Compton wavelength 1 µ λ, in order to derive appropriate coordinate conditions for which the stress-energy tensor for an initially axially-symmetric massive scalar field Φ (t, r) may be approximated by the stress-energy tensor for a uniformly-rotating (and collapsing) disk of dust, may be employed without significant modification from their spherically-symmetric counterparts. Our approach to guaranteeing the necessary smooth transition between geometries will be based on the general methods developed by Neugebauer and Meinel [START_REF] Neugebauer | The Einsteinian Gravitational Field of the Rigidly Rotating Disk of Dust[END_REF][48] [START_REF] Neugebauer | General Relativistic Gravitational Field of a Rigidly Rotating Disk of Dust: Solution in Terms of Ultraelliptic Functions[END_REF], building upon the previous work of Bardeen and Wagoner [START_REF] Bardeen | Uniformly Rotating Disks in General Relativity[END_REF][51], which we briefly recapitulate here.

The line element for the spacetime region surrounding such an axially-symmetric distribution of dust can be expressed in the cylindrical Weyl-Lewis-Papapetrou coordinates (t, ρ, ϕ, ζ) as:

ds 2 = -e 2ν(ρ,ζ) dt 2 + W 2 (ρ, ζ) e -2ν(ρ,ζ) (dϕ -ω (ρ, ζ)) 2 + e 2α(ρ,ζ) dρ 2 + dζ 2 , ( 190
)
where the metric potential functions α (ρ, ζ), W (ρ, ζ), ν (ρ, ζ) and ω (ρ, ζ) depend upon the non-angular spatial coordinates ρ and ζ only (i.e. significantly they are independent of the angular coordinate ϕ); these are sometimes known as Weyl's canonical coordinates [START_REF] Weyl | Zur Gravitationstheorie[END_REF]. In the absence of any rotation, i.e. in the limit as the angular velocity parameter ω (ρ, ζ) → 0, this reduces to the standard line element for the static and axially-symmetric Weyl metric:

lim ω(ρ,ζ)→0 ds 2 = -e 2ν(ρ,ζ) dt 2 + e -2ν(ρ,ζ) ρ 2 dϕ 2 + e 2α(ρ,ζ)-2ν(ρ,ζ) dρ 2 + dζ 2 . ( 191
)
Along the axis of symmetry, i.e. in the limit as the radial coordinate ρ → 0, the required existence of a local isometry to Minkowski space (the elementary flatness condition) furthermore implies that:

lim ρ→0 W (ρ, ζ) ρ e -(α(ρ,ζ)+ν(ρ,ζ)) = 1, ( 192
)
which furnishes one with a natural set of boundary conditions on the metric potential functions α (ρ, ζ), W (ρ, ζ) and ν (ρ, ζ). The other natural set of boundary conditions arises from enforcing the asymptotic flatness condition, namely the requirement that, at spatial infinity, i.e. in the limit as ρ 2 + ζ 2 → ∞, the resulting line element must converge asymptotically to the line element of the Minkowski metric in the cylindrical coordinate system: lim

ρ 2 +ζ 2 →∞ ds 2 = -dt 2 + dρ 2 + dζ 2 + ρ 2 dϕ 2 , (193)
implying the explicit boundary conditions:

lim ρ 2 +ζ 2 →∞ [α (ρ, ζ)] = 0, lim ρ 2 +ζ 2 →∞ [W (ρ, ζ)] = ρ, lim ρ 2 +ζ 2 →∞ [ν (ρ, ζ)] = 0, (194) and
α (ρ, ζ) = k (ρ, ζ) -U (ρ, ζ) , and e 2ν(ρ,ζ) W (ρ, ζ) ± ω (ρ, ζ) = 1 W (ρ, ζ) e -2U (ρ,ζ) ∓ a (ρ, ζ) , ( 196
)
it becomes possible to rewrite the axially-symmetric spacetime line element in the equivalent form:

ds 2 = -e 2U (ρ,ζ) (dt + a (ρ, ζ) dϕ) 2 + e -2U (ρ,ζ) e 2k(ρ,ζ) dρ 2 + dζ 2 + W 2 (ρ, ζ) dϕ 2 . ( 197
)
In order to obtain the metric for a uniformly-rotating disk of dust, we set:

W (ρ, ζ) = ρ, (198)
such that the metric can be recovered from the (complex) Ernst potential f (ρ, ζ):

f (ρ, ζ) = e 2U (ρ,ζ) + iψ (ρ, ζ) . ( 199
)
In the above, we are assuming an oriented spacetime with a timelike Killing vector field X, such that the twist 1-form, denoted τ , is given by:

* τ = ξ ∧ dξ, where ξ = X i dx i , (200)
in Cartesian coordinates. The function ψ (ρ, ζ) in the imaginary part of the Ernst potential f (ρ, ζ) arises from the fact that, within such a spacetime, one can (locally) write the twist 1-form as:

τ = dψ (ρ, ζ) , (201)
for some constant ψ along X. The real part e 2U (ρ,ζ) of the Ernst potential f (ρ, ζ) is then simply given by:

e 2U (ρ,ζ) = X i X i . ( 202
)
The collapse of an axially-symmetric, uniformly-rotating disk of dust can therefore be solved as a boundary value problem for the non-linear Ernst equation:

[Re (f (ρ, ζ))] ∂ ρρ f (ρ, ζ) + ∂ ρ f (ρ, ζ) ρ + ∂ ζζ f (ρ, ζ) = (∂ ρ f (ρ, ζ)) 2 + (∂ ζ f (ρ, ζ)) 2 , ( 203
)
with the new metric potential functions k (ρ, ζ) and a (ρ, ζ) being calculable by means of the following integrals:

k (ρ, ζ) = ρ 0 ρ (∂ ρU (ρ, ζ)) 2 -(∂ ζ U (ρ, ζ)) 2 + e -4U ( ρ,ζ) 4 (∂ ρψ (ρ, ζ)) 2 -(∂ ζ ψ (ρ, ζ)) 2 dρ, (204)
and:

a (ρ, ζ) = ρ 0 ρe -4U ( ρ,ζ) ∂ ζ ψ (ρ, ζ) dρ, (205) respectively.
Although the Kerr solution is conventionally parametrized in terms of a total mass M and a formal "angular momentum" parameter J, we choose instead to parametrize the uniformly-rotating disk of dust in terms of a coordinate radius ρ 0 and a function µ of the angular velocity parameter Ω (not to be confused with the mass parameter for the scalar field, which will no longer be required):

µ = 2Ω 2 ρ 2 0 e -2V0 , where V 0 = U (0, 0) , (206)
with V 0 playing the role of a "surface potential". These parameters are chosen in part to avoid the otherwise counterintuitive relationship J ≥ M 2 between the angular momentum and mass of the disk (which, in the case of the Kerr solution, would indicate supercriticality). If we now introduce the following, purely imaginary, function H (K):

H (K) = µ log 1 + µ 2 1 + K 2 ρ 2 0 2 + µ 1 + K 2 ρ 2 0 πiρ 2 0 1 + µ 2 1 + K 2 ρ 2 0 2 , with Re (H (K)) = 0, ( 207
)
as well as the complex functions Z (K, ρ, ζ) and Z 1 (K, ρ, ζ):

Z (K, ρ, ζ) = (K + i (ρ + iζ)) (K -i (ρ -iζ)) (K 2 -K 2 1 ) (K 2 -K 2 2 ), (208)
and:

Z 1 (K, ρ, ζ) = (K + i (ρ + iζ)) (K -i (ρ -iζ)), (209)
such that Re (Z 1 (K, ρ, ζ)) < 0 whenever ρ and ζ lie strictly outside the disk, then we can write the Ernst potential f (ρ, ζ) purely in terms of hyperelliptic integrals as:

f (ρ, ζ) = exp Ka(ρ,ζ) K1 K 2 Z (K, ρ, ζ) dK + K b (ρ,ζ) K2 K 2 Z (K, ρ, ζ) dK - iρ0 -iρ0 H (K) Z 1 (K, ρ, ζ) K 2 dK , (210)
where the lower limits of integration K 1 and K 2 are related by complex conjugation:

K 1 = -K2 = ρ 0 i -µ µ , with Re (K 1 ) = -Re (K 2 ) < 0, (211)
and where the upper limits of integration K a (ρ, ζ) and K b (ρ, ζ) are specified in terms of ultraelliptic functions (i.e. hyperelliptic functions of two variables, ρ and ζ), defined via:

Ka(ρ,ζ) K1 dK Z (K, ρ, ζ) + K b (ρ,ζ) K2 dK Z (K, ρ, ζ) = iρ0 -iρ0 H (K) Z 1 (K, ρ, ζ) dK, (212)
and:

Ka(ρ,ζ) K1 K Z (K, ρ, ζ) dK + K b (ρ,ζ) K2 K Z (K, ρ, ζ) dK = iρ0 -iρ0 H (K) Z 1 (K, ρ, ζ) KdK, (213) respectively.
The problem of recovering the upper integration limits K a (ρ, ζ) and K b (ρ, ζ) from these integral relations is now an instance of the Jacobi inversion problem for Abelian integrals, which can be solved explicitly in terms of (hyperelliptic) theta functions [82][83]. Specifically, if we define a theta function θ (x, y; p, q, α) using:

θ (x, y; p, q, α) = ∞ m=-∞ ∞ n=-∞ (-1) m+n p m 2 q n 2 e 2mx+2ny+4mnα , (214)
then the expression for the Ernst potential f (ρ, ζ) in terms of hyperelliptic integrals may be reformulated as:

f (ρ, ζ) = θ (α 0 u (ρ, ζ) + α 1 v (ρ, ζ) -C 1 (ρ, ζ) , β 0 u (ρ, ζ) + β 1 v (ρ, ζ) -C 2 (ρ, ζ) ; p, q, α) θ (α 0 u (ρ, ζ) + α 1 v (ρ, ζ) + C 1 (ρ, ζ) , β 0 u (ρ, ζ) + β 1 v (ρ, ζ) + C 2 (ρ, ζ) ; p, a, α) × e -(γ0u(ρ,ζ)+γ1v(ρ,ζ)+µw(ρ,ζ)) , (215)
with functions u (ρ, ζ), v (ρ, ζ) and w (ρ, ζ) given by the integrals:

u (ρ, ζ) = iρ0 -iρ0 H (K) Z 1 (K, ρ, ζ) dK, v (ρ, ζ) = iρ0 -iρ0 H (K) Z 1 (K, ρ, ζ) KdK, (216)
and:

w (ρ, ζ) = iρ0 -iρ0 H (K) Z 1 (K, ρ, ζ) K 2 dK, (217)
respectively, and where the normalization parameters α 0 , α 1 , β 0 , β 1 , γ 0 and γ 1 , as well as the moduli of the theta function p, q and α, and the pair of functions C 1 (ρ, ζ) and C 2 (ρ, ζ), are all defined over the two-sheeted hyperelliptic Riemann surface associated to the complex function Z (K, ρ, ζ). More precisely, this is the hyperelliptic Riemann surface with cuts between the branch points K 1 and K1 , K 2 and K2 , and -i (ρ -iζ) and i (ρ + iζ), and with upper/lower sheets defined by Z (K, ρ, ζ) → µK 3 in the limit as K → ∞.

Thus, for consistency, the integrals from K 1 to K a (ρ, ζ) and from K 2 to K b (ρ, ζ) must be evaluated along the same curve on this Riemann surface. The first four normalization parameters α 0 , α 1 , β 0 and β 1 can be obtained from dω 1 and dω 2 , namely two normalized Abelian differentials of the first kind, i.e:

dω 1 = α 0 dK Z (K, ρ, ζ) + α 1 K Z (K, ρ, ζ)
dK, and

dω 2 = β 0 dK Z (K, ρ, ζ) + β 1 K Z (K, ρ, ζ) dK, (218)
given by:

am dω n = πiδ mn , with m, n ∈ {1, 2} , (219)
which can then be solved as a system of linear algebraic equations in order to recover α 0 , α 1 , β 0 and β 1 , represented in terms of closed curves (periods) a 1 and a 2 . Likewise, the remaining two normalization parameters γ 0 and γ 1 can be obtained from dω 3 , a normalized Abelian differential of the third kind, i.e:

dω 3 = γ 0 dK Z (K, ρ, ζ) + γ 1 K Z (K, ρ, ζ) dK + µ K 2 Z (K, ρ, ζ) dK, (220)
whose periods a 1 and a 2 vanish:

aj dω 3 = 0, with j ∈ {1, 2} , (221)
to yield another system of linear algebraic equations for γ 0 and γ 1 . Moreover, the moduli of the theta function p, q and α may be derived from the Riemann matrix:

(B ij ) =    log (p) 2α 2α log (q)    , (222)
given by:

B ij = bi dω j , with i, j ∈ {1, 2} , (223)
such that the real part of (B ij ) is negative definite, with another pair of closed curves (periods) b 1 and b 2 .

Finally, the pair of functions C 1 (ρ, ζ) and C 2 (ρ, ζ) can be computed by means of the integral:

C i (ρ, ζ) = - ∞ + -i(ρ-iζ) dω i , with i ∈ {1, 2} , (224)
with + denoting the upper sheet of the hyperelliptic Riemann surface.

We are now able to determine that the surface mass density of the uniformly rotating disk of dust, denoted σ (µ, ρ), and given by:

σ (µ, ρ) = - Ω 2πe V0(µ) ψ 0 µ 1 -ρ 2 ρ 2 0 e V0 µ 1-ρ 2 ρ 2 0 , (225)
is strictly positive, where, in the above, the surface potential parameter V 0 (µ) now depends only upon the angular velocity-dependent parameter µ:

V 0 (µ) = - 1 2 sinh -1    µ + 1 + µ 2 ℘ I (µ) ; 4 3 µ 2 -4, 8 3 µ 1 + µ 2 9 -2 3 µ    , (226)
with ℘ (x; g 2 , g 3 ) being the Weierstrass elliptic function:

∞ ℘(x;g2,g3) dt 4t 3 -g 2 t -g 3 = x, (227)
and with the single-parameter function I (µ) given by:

I (µ) = 1 π µ 0 log x + √ 1 + x 2 (1 + x 2 ) (µ -x) dx. (228)
Moreover, the imaginary part of the Ernst potential ψ (ρ, ζ) now also depends only upon µ, hence the use of the notation ψ 0 (µ):

ψ (ρ, ζ) = ψ 0 µ 1 - ρ 2 ρ 2 0 , with ζ = 0 + . (229)
If µ 0 denotes the first zero of the Weierstrass elliptic function expression in the denominator of the equation for the surface potential V 0 (µ), i.e:

℘ I (µ) ; 4 3 µ 2 -4, 8 3 µ 1 + µ 2 9 - 2 3 µ = 0, (230)
then we have that the solution is non-singular (i.e. regular) whenever:

0 < µ < µ 0 ≈ 4.63, ⇐⇒ 0 > V 0 (µ) > -∞. (231)
Note that both along the axis of symmetry (ρ = 0) and within the plane of rotation of the disk (ζ = 0), the hyperelliptic integrals for the Ernst potential:

f (ρ, ζ) = exp Ka(ρ,ζ) K1 K 2 Z (K, ρ, ζ) dK + K b (ρ,ζ) K2 K 2 Z (K, ρ, ζ) dK - iρ0 -iρ0 H (K) Z 1 (K, ρ, ζ) K 2 dK , (232)
as well as the corresponding hyperelliptic integrals for the integration limits, all reduce to ordinary elliptic ones. In the former case, with ρ = 0, the potential may be written out as the following ratio of elliptic integrals:

f (0, ζ) = 2π + 1 -1 β(x) ix-ζ ρ 0 dx 2π + 1 -1 α(x) ix-ζ ρ 0 dx , with ζ > 0, (233)
with the function β (x) defined as satisfying the following integral equation:

β (x) = (2µ) 3 2 e -V0 x 1 -x 2 -µ 2 1 -x 2 2 β (x) + 1 -x 2 π 2 1 -1 1 -(x ) 2 x -x dx 1 -1 β (x ) x -x dx , (234)
with the two improper integrals designating the Cauchy principal values, and with the function α (x) depending algebraically upon both x and β (x) as follows:

α (x) = (1 -w (x)) β (x) + i 4w 2 (x) e -4V0 (ψ 2 0 + 4Ω 2 ρ 2 0 x 2 ) -(e 4V0 + w 2 (x)) β 2 (x) iψ 0 -2Ωρ 0 x , (235)
where, for the sake of notational convenience, we have introduced the scalar function w (x) and the scalar constant ψ 0 as:

w (x) = 2Ω 2 ρ 2 0 1 -x 2 ,
and

ψ 0 = -1 -e 4V0 -4Ω 2 ρ 2 0 , (236)
respectively. Likewise in the latter case, with ζ = 0 (within the disk), the metric potential functions U (ρ, ζ), a (ρ, ζ) and k (ρ, ζ) may be computed purely in terms of ordinary elliptic functions as:

e 2U (ρ,ζ) = exp 2V 0 µ 1 - ρ 2 ρ 2 0 - µρ 2 2ρ 2 0 , (237)
(1 + Ωa (ρ, ζ)) e 2U (ρ,ζ) = e V0(µ) exp V 0 µ 1 - ρ 2 ρ 2 0 , (238)
and:

e 2k(ρ,ζ)-2U (ρ,ζ) = e -2V0(µ) exp    - µ µ 1-ρ 2 ρ 2 0 f 0 (μ) f 0 (μ) f 0 (μ) f0 (μ) dμ    , (239)
respectively, at least in the all cases for which ρ ≤ ρ 0 , where the function f 0 (µ) represents the value of the Ernst potential f (ρ, ζ) at ρ = 0 and ζ = 0 + :

f 0 (µ) = e 2V0(µ) + iψ 0 (µ) . (240)
The angular velocity parameter (considered as a function of both µ and the coordinate radius ρ 0 ) can hence be computed by combining the original definition of the µ parameter:

µ = 2Ω 2 (µ, ρ 0 ) ρ 2 0 e -2V0(µ) , (241)
with the expression for the surface potential parameter V 0 (µ) in terms of the Weierstrass elliptic function ℘ (x; g 2 , g 3 ):

V 0 (µ) = - 1 2 sinh -1    µ + 1 + µ 2 ℘ I (µ) ; 4 3 µ 2 -4, 8 3 µ 1 + µ 2 9 -2 3 µ    . (242)
The Newtonian limit of the uniformly-rotating disk of dust configuration now corresponds to the case in which the parameter µ 1, in which case the solution reduces to the so-called "Maclaurin disk" solution (i.e. a gaseous self-gravitating astrophysical disk under Newtonian gravity). We can see this directly by expanding the solution for the Ernst potential (written with an explicit dependence on µ) f (ρ, ζ) as a formal power series, with expansion parameter √ µ, around the point µ = 0:

f (µ, ρ, ζ) = 1 + ∞ n=1 f n (ρ, ζ) µ (n+1) 2 . ( 243
)
From the expressions for the integration limits K a (ρ, ζ) and K b (ρ, ζ), as well as for the Ernst potential f (ρ, ζ) itself, in terms of hyperelliptic integrals, we therefore obtain the following leading-order terms in the respective expansions: we introduce the elliptic coordinate system (ξ, η) using:

K a (ρ, ζ) -K 1 = O µ 3 2 , K b (ρ, ζ) -K 2 = O µ 3 2 , (244) and:
f (µ, ρ, ζ) = exp -µ i -i H (K -K 1 ) (K -K 2 ) Z (K, ρ, ζ) dK + O µ 3 . (245)
ρ = ρ 0 1 + ξ 2 1 -η 2 , with ζ = ρ 0 ξη, (246)
such that 0 ≤ ξ < ∞ and -1 ≤ η ≤ 1, then the integral for the first coefficient f 1 (ξ, η) may be evaluated as:

f 1 (ξ, η) = 1 πi i -i 1 + K 2 Z (K, ρ, ζ) dK = - 1 π 4 3 cot -1 (ξ) + ξ -ξ 2 + 1 3 cot -1 (ξ) 1 -3η 2 . ( 247
)
The Newtonian (Maclaurin disk) limit U M ac (µ, ρ, ζ) as µ → 0, given by (with an explicit dependence of the speed of light c reintroduced temporarily for the sake of dimensional clarity):

U M ac (µ, ρ, ζ) = c 2 µf 1 (ρ, ζ) 2 , (248)
can thus be recovered by considering the time-time component g tt of the metric tensor g µν :

g tt (µ, ρ, ζ) = -Re (f (µ, ρ, ζ)) = -(1 + µf 1 (ρ, ζ)) = -1 + 2U M ac (µ, ρ, ζ) c 2 , (249)
such that:

µ = - 2U M ac (µ, 0, 0) c 2 = 2Ω 2 (µ, ρ) ρ 2 0 c 2 , (250)
as required. On the other hand, the extremal Kerr black hole limit is recovered as µ → µ 0 ≈ 4.63, as we shall now see.

Recall first that the spacetime line element for a Kerr black hole in the Boyer-Lindquist coordinate system (r, θ, ϕ) (i.e. the three-dimensional orthogonal coordinate system on an oblate spheroid) is traditionally given as:

ds 2 = -1 - 2M r Σ dt 2 + Σ ∆ dr 2 + Σdθ 2 + r 2 + a 2 + 2M ra 2 Σ sin 2 (θ) sin 2 (θ) dϕ 2 - 4M ra sin 2 (θ) Σ dtdϕ, (251)
with a being the ratio of the black hole's angular momentum J to its ADM mass M :

a = J M , (252)
and with Σ and ∆ being two characteristic length scales, namely:

Σ = r 2 + a 2 cos 2 (θ) , and ∆ = r 2 -2M r + a 2 , (253)
respectively. However, by introducing the additional parameters W , ν and ω, defined by:

W 2 = ∆ sin 2 (θ) , e 2ν = ∆Σ (r 2 + a 2 ) 2 -∆a 2 sin 2 (θ) , ω = 2M ra ∆Σ e 2ν , (254)
this line element can be rewritten in the following compressed (but equivalent) form that is closer to the Weyl-Lewis-Papapetrou cylindrical metric:

ds 2 = -e 2ν dt 2 + Σ dr 2 ∆ + dθ 2 + W 2 e -2ν (dϕ -ωdt) 2 . (255)
The stationary and axially-symmetric nature of the metric can be inferred from the lack of explicit dependence of the line element on time t or the angular coordinate ϕ, respectively. The Schwarzschild event horizon (otherwise known as the interior horizon) appears at the point where the radial-radial component of the metric tensor g rr becomes singular, namely:

1 g rr = 0, ⇐⇒ r I = M + M 2 -a 2 , (256)
i.e. the largest root of the quadratic equation ∆ = 0. Likewise, the boundary of the ergosphere (otherwise known as the exterior horizon) appears at the point where the time-time component of the metric tensor g tt changes sign, namely:

g tt = 0, ⇐⇒ r E = M + M 2 -a 2 cos 2 (θ). (257)
Note that, in contrast to the metric for the uniformly-rotating disk of dust, this solution is only physical if:

a ≤ M, ⇐⇒ J ≤ M 2 . ( 258
)
Within the ergosphere, i.e. for r I < r < r E , the discrepancy between the non-zero value of the Boyer-Lindquist shift vector β Kerr i and the (approximately zero) value of the numerical shift vector β i gives rise to an off-diagonal distortion in the metric tensor, causing coordinates to be dragged in the angular direction θ, implying that all observers must rotate in the same direction as the black hole itself:

∂ϕ ∂t > 0, (259)
which we interpret physically as frame-dragging.

The trajectories of test particles in the Kerr geometry, with instantaneous 4-momentum (p t , p r , p θ , p ϕ ), can be described in terms of three conserved quantities, namely the total energy E = -p t , the angular momentum component parallel to the axis of symmetry L = p ϕ , and the more abstract quantity Q, defined by:

Q = p 2 θ + cos 2 (θ) a 2 µ 2 -p 2 t + p 2 ϕ sin 2 (θ) , (260)
with µ being the trivial fourth conserved quantity, i.e. the particle's rest mass (without loss of generality, we henceforth assume µ = 0 for photons and µ = 1 for all other particles). Q = 0 is both a necessary and sufficient condition for a particle moving within the "equatorial plane" θ = π 2 to remain within the equatorial plane indefinitely, while Q > 0 is a necessary condition for a particle to cross the equatorial plane. By solving the equation for Q first for the lower-index momentum quantities p µ , followed by the upper-index momentum quantities p µ , we obtain the following equations of motion for the orbital trajectories of the test particles:

Σ dt dλ = -a aE sin 2 (θ) -L + r 2 + a 2 T ∆ , Σ dr dλ = ± V r , , (261)
Σ dθ dλ = ± V θ , Σ dϕ dλ = -aE - L sin 2 (θ) + a T ∆ , (262)
with λ being a function of the proper time τ on the particle's world line τ :

λ = τ µ , (263)
which thus becomes an affine parameter in the limit µ → 0, with the quantity T being given by:

T = E r 2 + a 2 -La, (264)
and with quantities V r and V θ being "effective potentials" that govern the motions of particles in the r and θ coordinate directions, respectively:

V r = T 2 -∆ µ 2 r 2 + (L -aE) 2 + Q , V θ = Q -cos 2 (θ) a 2 µ 2 -E 2 + L 2 sin 2 (θ) . ( 265
)
For a test particle in a circular orbit within the equatorial plane θ = π 2 at a fixed radius r, the derivative dr dλ must vanish identically, yielding:

Σ dr dλ = ± V r , =⇒ V r (r) = V r (r) = 0, (266)
where these equations can be solved simultaneously to give the total energy E and angular momentum component L as:

E µ = r 3 2 -2M √ r ± a √ M r 3 4 r 3 2 -rM √ r ± 2a √ M , and L µ = ± √ M r 2 ∓ 2a √ M r + a 2 r 3 4 r 3 2 -3M √ r ± 2a √ M , (267)
respectively, implying that the coordinate angular velocity of the orbit Ω is simply:

Ω = dϕ dt = ± √ M r 3 2 ± a √ M , (268)
with the upper signs designating the direct orbits (i.e. corotating orbits with the black hole, such that L > 0), and with lower signs designating the retrograde orbits (i.e. counterrotating with the black hole, such that L < 0).

Evidently, such circular orbits can only exist if the denominator in the equations for E µ and L µ is real, i.e. if and only if the following inequality is satisfied.

r 3 2 -3M √ r ± 2a √ M ≥ 0. (269)
The limiting case of equality occurs only in the case of photons (since this is associated with infinite energy E per unit rest mass µ), whose orbits therefore form the innermost boundary, denoted r ph , of all the circular orbits:

r ph = 2M 1 + cos 2 3 cos -1 ∓ a M . (270)
For all other orbits (with r > r ph ), the trajectories of certain particles may diverge to infinity in an asymptotically hyperbolic fashion if they are given an infinitesimal outward perturbation; such unstable "unbound" orbits occur for E µ > 1 (with the instability arising from the fact that, although their trajectories are given by circular functions, their energies are now given by hyperbolic ones), with the limiting case of equality E µ = 1 occurring in the case of "marginally bound" orbits, denoted r mb , forming the innermost boundary of all orbits with parabolic energy functions:

r mb = 2M ∓ a + 2 √ M √ M ∓ a. (271)
Parabolic trajectories (which are a good approximation for all non-relativistic matter v c) that intersect the region r < r mb will therefore inevitably fall into the black hole. On the other hand, bound orbits with curvature in ψ), respectively. The convergence rates for the Hamiltonian constraint after time t = 4.5M , with respect to the L 1 , L 2 and L ∞ -norms, illustrating approximately fourth-order convergence of the finitedifference scheme, are shown in Table 2. We also confirm that the ADM mass of the spacetime remains approximately constant, that the linear momentum of the resulting extremal Kerr black hole converges to be approximately zero, and that the (ADM) angular momentum is approximately conserved (and therefore the angular momentum parameter J of the limiting extremal Kerr black hole is approximately equal to the initial angular momentum of the rotating disk), as required. Figure 23: Spatial hypergraphs corresponding to the initial hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 0M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted using the local curvature in the Boyer-Lindquist conformal factor ψ, and colored according to the value of the scalar field Φ (t, R).

Vertices (L 1 ) (L 2 ) (L ∞ ) O (L 1 ) O (L 2 ) O (L ∞ ) 100
Figure 24: Spatial hypergraphs corresponding to the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 1.5M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted using the local curvature in the Boyer-Lindquist conformal factor ψ, and colored according to the value of the scalar field Φ (t, R).

Figure 25: Spatial hypergraphs corresponding to the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 3M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted using the local curvature in the Boyer-Lindquist conformal factor ψ, and colored according to the value of the scalar field Φ (t, R).

Figure 26: Spatial hypergraphs corresponding to the final hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted using the local curvature in the Boyer-Lindquist conformal factor ψ, and colored according to the value of the scalar field Φ (t, R).

Figure 27: Spatial hypergraphs corresponding to the initial hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 0M , with resolutions of 200, 400 and 800 vertices, respectively. The hypegraphs have been adapted and colored using the local curvature in the Boyer-Lindquist conformal factor ψ. Figure 32: Spatial hypergraphs corresponding to projections along the z-axis of the initial hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 0M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Boyer-Lindquist conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted using the local curvature in ψ, and colored according to the value of the scalar field Φ (t, R).

Figure 33: Spatial hypergraphs corresponding to projections along the z-axis of the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 1.5M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Boyer-Lindquist conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted using the local curvature in ψ, and colored according to the value of the scalar field Φ (t, R).

Figure 34: Spatial hypergraphs corresponding to projections along the z-axis of the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 3M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Boyer-Lindquist conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted using the local curvature in ψ, and colored according to the value of the scalar field Φ (t, R).

Figure 35: Spatial hypergraphs corresponding to projections along the z-axis of the final hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Boyer-Lindquist conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted using the local curvature in ψ, and colored according to the value of the scalar field Φ (t, R).

Figure 36: Spatial hypergraphs corresponding to projections along the z-axis of the initial hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 0M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Boyer-Lindquist conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted and colored using the local curvature in ψ.

Figure 37: Spatial hypergraphs corresponding to projections along the z-axis of the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 1.5M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Boyer-Lindquist conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted and colored using the local curvature in ψ.

Figure 38: Spatial hypergraphs corresponding to projections along the z-axis of the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 3M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Boyer-Lindquist conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted and colored using the local curvature in ψ.

Figure 39: Spatial hypergraphs corresponding to projections along the z-axis of the final hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 vertices, respectively. The vertices have been assigned spatial coordinates according to the profile of the Boyer-Lindquist conformal factor ψ through a spatial slice perpendicular to the z-axis, and the hypergraphs have been adapted and colored using the local curvature in ψ.

Comparison with the Pure Wolfram Model Case

Our objective for this section is to demonstrate that the numerical results obtained via the Gravitas framework within the preceding two sections may be reproduced via a pure hypergraph rewriting system that provably satisfies the Einstein field equations in the continuum limit (with an event selection function determined by the ADM gauge conditions presented within Section 2), without recourse to discretizing an a priori continuous spacetime metric. Recall that the Wolfram model hypergraph rewriting formalism [START_REF] Wolfram | A New Kind of Science[END_REF][19] [20][21] may, at least in a large class of cases, be considered analogous to a full discretization of the Cauchy problem in general relativity [START_REF] Gorard | Hypergraph Discretization of the Cauchy Problem in General Relativity via Wolfram Model Evolution[END_REF], with each spacelike hypersurface being given by a (finite, usually directed) hypergraph H = (V, E), defined by:

E ⊆ P (V ) \ {∅} , (295)
for power set function P, i.e. each hyperedge e ∈ E connects an arbitrary (non-empty) subset of vertices in V . Such spatial hypergraphs may therefore be represented as finite collections of (ordered) relations between vertices, as illustrated in Figure 40. Once the Cauchy initial data have been specified through the construction of an appropriate initial spatial hypergraph (usually by Poisson sprinkling [START_REF] Gorard | Algorithmic Causal Sets and the Wolfram Model[END_REF]), the Cauchy surface may then be evolved forwards in time by applying an abstract rewriting rule R defined over hypergraphs, of the general form:

R :

H 1 = (V 1 , E 1 ) → H 2 = (V 2 , E 2 ) . (296)
Loosely speaking, such a rule specifies that a subhypergraph which matches the pattern defined by

H 1 = (V 1 , E 1 )
is to be replaced by a distinct subhypergraph matching the pattern defined by H 2 = (V 2 , E 2 ), as illustrated in the example shown in Figure 41; the hypergraph rewriting semantics of Wolfram model evolution can be made fully rigorous using the formalism of double-pushout rewriting over selective adhesive categories [START_REF] Gorard | ZX-Calculus and Extended Hypergraph Rewriting Systems I: A Multiway Approach to Categorical Quantum Information Theory[END_REF][85] [START_REF] Gorard | A Functorial Perspective on (Multi)computational Irreducibility[END_REF] If we proceed to apply this rule to every possible/matching (non-overlapping) subhypergraph at each step, it defines an effective dynamics for the Wolfram model system, as demonstrated in Figure 42 for the case of the rewriting rule {{x, y} , {y, z}} → {{w, y} , {y, z} , {z, w} , {x, w}}, applied to the "double self-loop" initial condition {{0, 0} , {0, 0}}. Since there will, in general, exist multiple such choices of maximal non-overlapping sets of subhypergraphs to which to apply the rule at each step, there exists a certain degree of freedom over the evolution history of the system; this corresponds to the gauge freedom inherent to the solution of the general relativistic Cauchy problem [START_REF] Gorard | Some Relativistic and Gravitational Properties of the Wolfram Model[END_REF][23] [START_REF] Gorard | Hypergraph Discretization of the Cauchy Problem in General Relativity via Wolfram Model Evolution[END_REF]. If the sequence of spatial hypergraphs represents the foliation of the discrete spacetime into spacelike hypersurfaces, then the causal (i.e. conformally invariant) structure of that discrete spacetime may be represented by means of a directed acyclic graph, known as a causal graph. Within such a causal graph, each vertex corresponds to an individual rewrite event, and each edge corresponds to a causal relationship between a pair of rewrite events, such that the directed edge A → B exists if and only if:

In (B) ∩ Out (A) = ∅, (297)
i.e. if and only if the input for rewrite event B makes use of hyperedges that were produced as part of the output of rewrite event A; again, the causal structure of Wolfram model systems may be made fully rigorous using the formalism of weak 2-categories and multiway evolution causal graphs [START_REF] Gorard | Fast Automated Reasoning over String Diagrams using Multiway Causal Structure[END_REF]. An example illustrating how such a causal graph may be constructed algorithmically is shown in Figure 43, and it is shown in Figure 44 how taking the transitive reduction of such a causal graph yields the Hasse diagram for a causal set, i.e. it represents the conformal structure of a kind of "skeletonized" version of a Lorentzian manifold [START_REF] Gorard | Algorithmic Causal Sets and the Wolfram Model[END_REF]. The condition of causal invariance guarantees that the causal graphs produced by different choices of updating sequence in the hypergraph (and therefore different choices of relativistic gauge) will all ultimately be isomorphic, which in turn guarantees compatibility with general covariance in the continuum limit. The procedure by which Wolfram model/hypergraph rewriting rules may be constructed that provably satisfy the Einstein field equations in the continuum limit is somewhat complex and has been outlined in great detail elsewhere [START_REF] Gorard | Some Relativistic and Gravitational Properties of the Wolfram Model[END_REF][23] [START_REF] Gorard | Hypergraph Discretization of the Cauchy Problem in General Relativity via Wolfram Model Evolution[END_REF], and so in the interest of concision we do not seek to reproduce those details here. The salient features of that construction for our present purposes are as follows: the condition of causal invariance allows one to define a compatible notion of discrete intrinsic curvature (based on the Ollivier-Ricci curvature construction for arbitrary metric-measure spaces [START_REF] Ollivier | Ricci curvature of metric spaces[END_REF][90]) on both causal graphs and hypergraphs, and hence to construct a discrete analog of the Einstein-Hilbert action by taking the sum of the discrete Ricci scalar over the entire causal graph; if one further assumes that the hypergraph rewriting dynamics are weakly ergodic, such that, in particular, the net flux of causal edges through any given hypersurface in the causal graph eventually converges to zero, then one is analytically justified in exchanging this discrete sum for an Figure 42: An elementary example of a 10-step evolution history for the Wolfram model/hypergraph rewriting system corresponding to the set substitution rule {{x, y} , {y, z}} → {{w, y} , {y, z} , {z, w} , {x, w}}, starting from the "double self-loop" initial condition {{0, 0} , {0, 0}}.

integral in the continuum limit; finally, if one assumes that the hypergraph rewriting rule is asymptotically dimension preserving, such that the dimension of the limiting hypergraph converges to some fixed, finite value, then it follows that the discrete Einstein-Hilbert action must be extremized in this limit. Thus, any Wolfram model rule which is causal invariant, weakly ergodic and asymptotically dimension preserving will be consistent with the Einstein field equations in the continuum limit. Note that here, as elsewhere, isomorphisms between causal graphs and spatial hypergraphs are determined by means of a generalized form of the "uniqueness trees" algorithm [START_REF] Gorard | Uniqueness Trees: A Possible Polynomial Approach to the Graph Isomorphism Problem[END_REF].

Figure 43:

An elementary example of a causal graph obtained after a 3-step evolution for the Wolfram model/hypergraph rewriting system corresponding to the set substitution rule {{x, y} , {x, z}} → {{x, y} , {x, w} , {y, w} , {z, w}}, starting from the "double self-loop" initial condition {{0, 0} , {0, 0}}. The complete causal graph obtained after a 5-step evolution for the Wolfram model/hypergraph rewriting system corresponding to the set substitution rule {{x, y} , {x, z}} → {{x, y} , {x, w} , {y, w} , {z, w}} (left), shown together with its transitive reduction (right), starting from the "double self-loop" initial condition {{0, 0} , {0, 0}}.

Although the majority of attention regarding the relativistic properties of the Wolfram model has thus far been focused on the case of vacuum spacetimes, in recent work [START_REF] Gorard | Axiomatic Quantum Field Theory in Discrete Spacetime via Multiway Causal Structure: The Case of Entanglement Entropies[END_REF] we showed how it is possible to equip an arbitrary casual graph (in particular, one generated via Wolfram model evolution) with the structure of a free, massless scalar field. In so doing, we broadly followed the approach of Johnston [START_REF] Johnston | Quantum Fields on Causal Sets[END_REF] in which one begins by applying the standard Fourier-analytic techniques described by Gel'fand and Shilov [START_REF] Gel | Generalized Functions[END_REF] and Egorov and Shubin [START_REF] Egorov | Partial Differential Equations II[END_REF] in order to derive retarded Green's functions for arbitrary d-dimensional (continuous) spacetimes; for instance, in d = 1, d = 2, d = 3 and d = 4, these are given by:

(G R ) (1) m (x) = θ (x) sin (mx) m , (G R ) (2) m (x) = θ x 0 θ τ 2 1 2 J 0 (mτ ) , (298)
(G R ) (3)
m (x) = θ x 0 θ τ 2 1 2π cos (mτ ) τ , (G R ) (4) m (x) = θ x 0 θ τ 2 1 2π δ τ 2 - m 4π J 1 (mτ ) τ , (299)
density profile for the scalar field Φ (t, R), embedded within a Poisson-sprinkled initial hypergraph with fixed sprinkling density ρ c ) for an axially-symmetric background metric in (cylindrical) Weyl-Lewis-Papapetrou coordinates: as well as the extraction of all associated curvature and stress-energy properties of the resulting discrete spacetimes, can be found in [START_REF] Gorard | Hypergraph Discretization of the Cauchy Problem in General Relativity via Wolfram Model Evolution[END_REF].

ds 2 = -e 2U (ρ,zeta) (dt + a (ρ, ζ) dϕ) 2 + e -2U (ρ,ζ) e 2k(ρ,ζ) dρ 2 + dζ 2 + W 2 (ρ, ζ) dϕ 2 , (322)
Vertices (L 1 ) (L 2 ) (L ∞ ) O (L 1 ) O (L 2 ) O (L ∞ ) 100

Concluding Remarks

This article has succeeded in demonstrating, through a combination of rigorous mathematical analysis and explicit numerical simulation using the Gravitas framework, that if one starts with an analytical solution to the massive scalar field "bubble collapse" problem in either spherical or axial symmetry (yielding convergence a non-rotating Schwarzschild black hole or a maximally-rotating/extremal Kerr black hole, respectively), and proceeds to apply microscopic perturbations away from these exact symmetries, at least of the kind that are consistent with discretization of the underlying spacetime, then the solution structure in both cases remains stable under appropriate assumptions. We have, in addition, shown that both solutions may be recovered as limiting cases of a pure Wolfram model/hypergraph rewriting evolution, without any a priori continuity assumptions on the underlying spacetime. To the best of our knowledge, this represents the first systematic analysis of the evolution of a Wolfram model system representing a true non-vacuum solution to the Einstein field equations in the continuum limit, as well as the first robust analysis of gravitational collapse and singularity-formation behavior within a reasonably generic discrete spacetime setting. However, the techniques employed here are clearly exceedingly limited in scope, depending as they do upon an unphysically high degree of a priori symmetry, as well as an extremely idealized matter model (i.e. a massive scalar field obeying the discrete Klein-Gordon equation) for the stress-energy tensor. Recent work investigating global homotopic aspects of Wolfram model systems [START_REF] Arsiwalla | Homotopies in Multiway (Non-Deterministic) Rewriting Systems as n-Fold Categories[END_REF][96] provides one plausible direction by which it might conceivably be possible to formulate a more global topological theorem regarding geodesic incompleteness in generic Wolfram model systems, in a similar style to Penrose's original 1964 singularity theorem [START_REF] Penrose | Gravitational Collapse and Space-Time Singularities[END_REF].

Moreover, there is an ongoing effort to integrate more sophisticated stress-energy models (e.g. relativistic dust, relativistic radiation, perfect fluids, electromagnetic fields, etc.) into the Gravitas framework, and thereby to make such stress-energy distributions compatible with Wolfram model evolution in broadly the same manner as we have done here with massive scalar fields; it would be highly interesting to investigate the analogs of relativistic energy conditions and equations of state within arbitrary Wolfram model causal graphs, and to probe the effects that such properties and constraints might have on the dynamics and singularity structure of the resulting discrete spacetimes. Such an extension might also afford one the opportunity to extend certain, presently highly idealized, quantum mechanical computations regarding (e.g.) black hole entropies [97][98] to the case of more astrophysically realistic black holes, such as those produced through the gravitational collapse of stellar matter. One might potentially even be able to use such techniques to study various open cosmological questions, such as the formation of supermassive primordial black holes within certain inflationary models [START_REF] Zel | The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model[END_REF], in the discrete spacetime setting. Finally, we note that it would be particularly exciting to study certain relativistic critical phenomena, such as the famous Choptuik scalar field collapse solution [START_REF] Choptuik | Universality and scaling in gravitational collapse of a massless scalar field[END_REF], by means of such methods, in order to assess whether and to what extent to which the scaling relations and universality properties observed within standard simulations of continuous spacetimes [START_REF] Gundlach | Critical Phenomena in Gravitational Collapse[END_REF] also extend to discrete ones.

Figure 3 :

 3 Figure 3: An illustration of the refinement/coarsening procedure for hypergraphs, wherein a hyperedge e = {u, v} is subdivided into a pair of hyperedges e 1 = {u, w}, e 2 = {w, v} for some newly-generated vertex w (i.e. refinement, on the left), or a pair of hyperedges e 1 = {u, w}, e 2 = {w, v} is smoothed into a single hyperedge e = {u, v} (i.e. coarsening, on the right).

Figure 4 :

 4 Figure 4: An illustration of the refinement procedure for (structured) hypergraphs, wherein a collection of twelve colored vertices in a two-dimensional quadrilateral mesh are marked for refinement (left), and are consequently replaced with finer (2-by-2) two-dimensional grids (right).

Figure 5 :

 5 Figure5: An illustration of which coarse vertices (highlighted in dark red) must have their finite-difference evolution scheme modified by virtue of being adjacent to interfaces with a fine grid/subhypergraph, following refinement of the two-dimensional quadrilateral mesh shown above.

2 1 + x 2 2 + x 2 3

 23 denoting the radial distance parameter, f 0 denoting the values of the specified scalar field at the boundary of the spacetime (which, for the purposes of the tests presented in this article, are taken to be the Minkowski space reference values), and v denoting the radiative velocity (henceforth, we take v = 1). We evolve the solution until a final time of t = 4.5M , with intermediate checks at times t = 1.5M and t = 3M ; the initial, first intermediate, second intermediate and final hypersurface configurations, with the hypergraphs adapted using the Schwarzschild conformal factor ψ and colored using the scalar field Φ (t, r), are shown in Figures6, 7, 8 and 9, respectively, with resolutions of 200, 400 and 800 vertices; similarly, Figures 10, 11, 12 and 13 show the initial, first intermediate, second intermediate and final hypersurface configurations, but with the hypergraphs both adapted and colored using the Schwarzschild conformal factor ψ, respectively. Figure 14 shows the discrete characteristic structure of the solutions after time t = 4.5M (using directed acyclic causal graphs to show discrete characteristic lines). Projections along the z-axis of the initial, first intermediate, second intermediate and final hypersurface configurations, with vertices assigned spatial coordinates according to the profile of the Schwarzschild conformal factor ψ, are shown in Figures 15, 16

  , 17 and 18 (with hypergraphs colored using the scalar field Φ (t, r)) and Figures 19, 20, 21 and 22 (with hypergraphs colored using the local curvature in ψ), respectively.

Figure 6 :

 6 Figure6: Spatial hypergraphs corresponding to the initial hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 0M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted using the local curvature in the Schwarzschild conformal factor ψ, and colored according to the value of the scalar field Φ (t, r).

Figure 7 :

 7 Figure7: Spatial hypergraphs corresponding to the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 1.5M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted using the local curvature in the Schwarzschild conformal factor ψ, and colored according to the value of the scalar field Φ (t, r).

Figure 8 :

 8 Figure8: Spatial hypergraphs corresponding to the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 3M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted using the local curvature in the Schwarzschild conformal factor ψ, and colored according to the value of the scalar field Φ (t, r).

Figure 12 :

 12 Figure12: Spatial hypergraphs corresponding to the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 3M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted and colored using the local curvature in the Schwarzschild conformal factor ψ.

Figure 13 :

 13 Figure13: Spatial hypergraphs corresponding to the final hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted and colored using the local curvature in the Schwarzschild conformal factor ψ.

Figure 14 :

 14 Figure14: Causal graphs corresponding to the discrete characteristic structure of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 hypergraph vertices, respectively.

: lim ρ 2

 2 +ζ 2 →∞ [ω (ρ, ζ)] = 0. (195) By introducing the alternative metric potential functions U (ρ, ζ), k (ρ, ζ) and a (ρ, ζ), relating them to the old functions α (ρ, ζ), W (ρ, ζ), ν (ρ, ζ) and ω (ρ, ζ) as follows:

By treating K = O ( 1 )

 1 within this integral and expanding the complex function Z (K, ρ, ζ), we can straightforwardly determine the expansion coefficients f n ρ ρ0 , ζ ρ0 in terms of elementary functions. For instance, if

Figure 28 :

 28 Figure28: Spatial hypergraphs corresponding to the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 1.5M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted and colored using the local curvature in the Boyer-Lindquist conformal factor ψ.

Figure 29 :

 29 Figure29: Spatial hypergraphs corresponding to the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 3M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted and colored using the local curvature in the Boyer-Lindquist conformal factor ψ.

Figure 30 :

 30 Figure30: Spatial hypergraphs corresponding to the final hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted and colored using the local curvature in the Boyer-Lindquist conformal factor ψ.

Figure 31 :

 31 Figure31: Causal graphs corresponding to the discrete characteristic structure of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 4.5M , with resolutions of 200, 400 and 800 hypergraph vertices, respectively.

Figure 40 :

 40 Figure 40: Two elementary examples of directed spatial hypergraphs, corresponding to the collections of ordered relations {{1, 2} , {1, 3} , {2, 3} , {4, 1}} and {{1, 2, 3} , {3, 4, 5}}, respectively.

Figure 41 :

 41 Figure 41: An elementary example of a hypergraph rewriting rule, corresponding to the Wolfram model/set substitution system {{x, y} , {y, z}} → {{w, y} , {y, z} , {z, w} , {x, w}}.

Figure 44 :

 44 Figure 44:The complete causal graph obtained after a 5-step evolution for the Wolfram model/hypergraph rewriting system corresponding to the set substitution rule {{x, y} , {x, z}} → {{x, y} , {x, w} , {y, w} , {z, w}} (left), shown together with its transitive reduction (right), starting from the "double self-loop" initial condition {{0, 0} , {0, 0}}.

Figure 47 :

 47 Figure 47: Spatial hypergraphs corresponding to the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 3M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the value of the scalar field Φ (t, r).

Figure 48 :

 48 Figure 48: Spatial hypergraphs corresponding to the final hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 4.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the value of the scalar field Φ (t, r).

Figure 49 :

 49 Figure49: Spatial hypergraphs corresponding to the initial hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 0M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the local curvature in the Schwarzschild conformal factor ψ.

Figure 50 :

 50 Figure 50: Spatial hypergraphs corresponding to the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 1.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the local curvature in the Schwarzschild conformal factor ψ.

Figure 51 :

 51 Figure51: Spatial hypergraphs corresponding to the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 3M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the local curvature in the Schwarzschild conformal factor ψ.

Figure 52 :

 52 Figure 52: Spatial hypergraphs corresponding to the final hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 4.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the local curvature in the Schwarzschild conformal factor ψ.

Figure 53 :

 53 Figure 53: Causal graphs corresponding to the discrete characteristic structure of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 4.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 hypergraph vertices, respectively.

Figure 54 :

 54 Figure54: Spatial hypergraphs corresponding to the initial hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 0M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the value of the scalar field Φ (t, R).

Figure 55 :

 55 Figure55: Spatial hypergraphs corresponding to the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 1.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the value of the scalar field Φ (t, R).

Figure 56 :

 56 Figure56: Spatial hypergraphs corresponding to the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 3M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the value of the scalar field Φ (t, R).

Figure 57 :

 57 Figure57: Spatial hypergraphs corresponding to the final hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 4.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the value of the scalar field Φ (t, R).

Figure 58 :

 58 Figure58: Spatial hypergraphs corresponding to the initial hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 0M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the local curvature in the Boyer-Lindquist conformal factor ψ.

Figure 59 :

 59 Figure 59: Spatial hypergraphs corresponding to the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 1.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the local curvature in the Boyer-Lindquist conformal factor ψ.

Figure 60 :

 60 Figure 60: Spatial hypergraphs corresponding to the second intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 3M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the local curvature in the Boyer-Lindquist conformal factor ψ.

Figure 61 :

 61 Figure61: Spatial hypergraphs corresponding to the final hypersurface configuration of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 4.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the local curvature in the Boyer-Lindquist conformal factor ψ.

Figure 62 :

 62 Figure 62: Causal graphs corresponding to the discrete characteristic structure of the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with a spinning (exponential) initial density distribution, at time t = 4.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 hypergraph vertices, respectively.

  r) + k

	(cos) 01	(t, r) cos (µt) + k	(sin) 01	(t, r) sin (µt)
				+	∞ m=1	2m n=2	1 µ m k (cos) mn (t, r) cos (nµt) + k (sin) mn (t, r) sin (nµt) , (117)
	for coefficient functions k	(cos) mn (t, r) and k mn (t, r), and: (sin)
	m (t, r) = m 0 (t, r) + m (cos) 01	(t, r) cos (µt) + m	(sin) 01	(t, r) sin (µt)
				+	∞ m=1	2m n=2	1 µ m m (cos)

mn (t, r) cos (nµt) + m (sin) mn (t, r) sin (nµt) , (118) for coefficient functions m (cos) mn (t, r) and m (sin) mn (t, r), respectively, from which we obtain (up to second order in the expansion parameter 1 µ ):

Table 1 :

 1 Convergence rates for the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with respect to the L 1 , L 2 and L ∞ -norms for the Hamiltonian constraint H after time t = 4.5M , showing approximately fourth-order convergence.Figure11: Spatial hypergraphs corresponding to the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 1.5M , with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been adapted and colored using the local curvature in the Schwarzschild conformal factor ψ.

	-3	-	-	-

[START_REF] Penrose | Gravitational Collapse and Space-Time Singularities[END_REF]

.36 × 10 -3 9.05 × 10 -3 2.94 × 10

Table 2 :

 2 [START_REF] Raychaudhuri | Relativistic Cosmology. I[END_REF].43 × 10 -2 8.99 × 10 -2 8.14 × 10 -2 Convergence rates for the massive scalar field "bubble collapse" problem to a maximally-rotating (extremal) Kerr black hole test, with respect to the L 1 , L 2 and L ∞ -norms for the Hamiltonian constraint H after time t = 4.5M , showing approximately fourth-order convergence.

	-	-	-

  with ρ = R sin (θ) and ζ = R cos (θ). Once again, the initial, first intermediate, second intermediate and final hypersurface configurations (corresponding to times t = 0M , t = 1.5M , t = 3M and t = 4.5M , respectively), with hypergraphs colored using the scalar field Φ (t, R), are shown in Figures 54, 55, 56 and 57, respectively, with resolutions of 200, 400 and 800 vertices; similarly, Figures 58, 59, 60 and 61 show the initial, first intermediate, second intermediate and final hypersurface configurations, but with the hypergraphs colored using the Boyer-Lindquist conformal factor ψ, respectively. Figure 62 shows the corresponding causal graphs for the resulting Wolfram model evolutions, with the convergence rates for the Hamiltonian constraint after time t = 4.5M , with respect to the L 1 , L 2 and L ∞ -norms, illustrating approximately second-order convergence, shown in Table4. In both cases we confirm approximate conservation of the ADM mass, linear momentum and angular momentum of the overall spacetime for the duration of the evolution. Full details regarding the construction and evolution of general relativistic Cauchy data via Wolfram model systems,

Table 3 :

 3 Convergence rates for the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with respect to the L 1 , L 2 and L ∞ -norms for the Hamiltonian constraint H after time t = 4.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), showing approximately second-order convergence.Figure45: Spatial hypergraphs corresponding to the initial hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 0M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the value of the scalar field Φ (t, r).

	-2	-	-	-

1.50 × 10 -1 1.79 × 10 -1 1.76 × 10

Table 4 :

 4 Convergence rates for the massive scalar field "bubble collapse" to a maximally-rotating (extremal) Kerr black hole test, with respect to the L 1 , L 2 and L ∞ -norms for the Hamiltonian constraint H after time t = 4.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), showing approximately second-order convergence.Figure46: Spatial hypergraphs corresponding to the first intermediate hypersurface configuration of the massive scalar field "bubble collapse" to a non-rotating Schwarzschild black hole test, with an exponential initial density distribution, at time t = 1.5M , produced via pure Wolfram model evolution (with set substitution rule {{x, y} , {y, z} , {z, w} , {w, v}} → {{y, u} , {u, v} , {w, x} , {x, u}}), with resolutions of 200, 400 and 800 vertices, respectively. The hypergraphs have been colored according to the value of the scalar field Φ (t, r).

,j,k can be used to correct the solutions over the

,j,k = Φ n i,j,k +

,j,k = Φ n i+1,j,k -

,j,k Φ n i+2,j,k -Φ i+1,j,k ,(80)for right boundary extrapolation, with (diagonal) limiter functions φ ± , and replacing the linearized approximation scheme on the right-hand side with a higher-degree spatial reconstruction polynomial w s,x h (x, t n ). This polynomial (again assuming reconstruction in the x coordinate direction) may be expanded in terms of

Acknowledgments

The author would like to thank the DAMTP GR journal club at the University of Cambridge (and particularly Justin Ripley) for providing the opportunity to present an early version of these results and receive useful feedback, Salvatore Vultaggio for clarifying discussions regarding the physical significance of scalar fields in general relativity and scalar-tensor gravity theories, José M. Martín-García for providing enlightening pointers to recent literature on relativistic critical phenomena in gravitational collapse, and Stephen Wolfram for general encouragement and helpful discussions throughout.

r > r mb will be stable if and only if V r (r) ≤ 0, yielding the equivalent conditions:

Thus, the "marginally stable" orbits, denoted r ms , forming the innermost boundary of all stable orbits r > r ms , will be given by:

where, for compactness, we have introduced the constants Z 1 and Z 2 :

The two limiting cases for the Kerr black hole consequently correspond to the non-rotating a = 0 (J = 0) case in which the solution simply reduces to the Schwarzschild metric, since the interior and exterior horizon radii now both coincide with the Schwarzschild radius (i.e. r I = r E = 2M ) and so the ergosphere disappears, and the maximally-rotating a = M (J = M 2 ) case, which we refer to here as the "extremal Kerr" solution.

For the extremal a = M case, the interior and exterior horizons are maximally separated, so the ergosphere becomes maximally extended, since r I = M and, in the equatorial plane with θ = π 2 , r E = 2M . The equations for the specific total energy E µ and specific angular momentum L µ given above reduce in this case to:

respectively. For the Schwarzschild (a = 0) case, the photon orbit radius is given by r ph = 3M , the marginally bound orbit radius is given by r mb = 4M , and the marginally stable orbit radius is given by r ms = 6M . On the other hand, for the extremal Kerr (a = M ) case, we can infer from the reduced equations above that the photon orbit radius is given by either r ph = M (for direct orbits) or r ph = 4M (for retrograde orbits), the marginally bound orbit radius is given by r mb = M (for direct orbits) or

(for retrograde orbits), and the marginally stable orbit radius is given by r ms = M (for direct orbits) or r ms = 9M (for retrograde orbits). Switching momentarily from the cylindrical coordinate system (ρ, ζ) to the spherical coordinate system (R, θ):

this allows us to write the complex Ernst potential f (ρ, ζ) for the extremal Kerr solution in the following, rather elegant, form:

assuming a strictly positive radial coordinate value R > 0.

In order to see how this extremal Kerr form of the complex Ernst potential is recovered in the limit as µ → µ 0 ≈ 4.63 from the more general Ernst potential f (ρ, ζ) for a uniformly-rotating disk of dust:

with upper integration limits K a (ρ, ζ) and K b (ρ, ζ) defined as ultraelliptic functions via:

and:

respectively, it is helpful to take K b (ρ, ζ) to be on the other sheet of the Riemann surface, as compared to the previous setup. This yields the following equivalent form for the Ernst potential:

with the corresponding equivalent forms for the integration limits:

and:

respectively. From the original definitions of the µ parameter:

and the surface potential parameter V 0 (µ):

we obtain, in the limit as µ → µ 0 (since µ 0 is, by definition, the first zero of the elliptic Weierstrass function

appearing in the denominator), the following limiting expressions for the integrals above in terms of the spherical coordinate system (R, θ):

iρ0 -iρ0

modulo considerations of periodicity, and under the hypothesis that the radial coordinate is strictly positive, i.e. R > 0. In the limit as µ → µ 0 , both K 1 and K 2 (i.e. the lower limits of integration) vanish:

under the hypothesis that Re (K 1 ) = -Re (K 2 ) < 0, and therefore the complex function Z (K, ρ, ζ) in the above integrals simplifies to:

The unique solution to the above (now elementary) integral equations is therefore given by the complex Ernst potential for the extremal Kerr solution:

as required. This completes the proof.

We now proceed to simulate the evolution of an axially-symmetric massive scalar field "bubble collapse" problem numerically, corresponding to the collapse of a uniformly-rotating disk of dust to a maximallyrotating (extremal) Kerr black hole, using the same techniques as for the spherically-symmetric case analyzed previously. As before, we assume an exponential initial density profile for the massive scalar field Φ (t, R):

for initial density constant ρ 0 and radius of support λ, albeit replacing the spherically-symmetric background metric in Gaussian polar coordinates:

with the axially-symmetric background metric in (cylindrical) Weyl-Lewis-Papapetrou coordinates:

τ denotes proper time distance:

J α denote Bessel functions of the first kind (of order α):

assuming validity of the series expansion around x = 0, and δ is the standard (1-dimensional) Dirac delta function, with the corresponding advanced Green's functions simply given by a slight modification of the boundary conditions:

By definition, such Green's functions correspond to solutions to the (massive) Klein-Gordon equation for a (real) scalar field:

and:

When used in conjunction with Johnston's highly evocative "hops and stops" formalism for causal set

Feynman propagators [START_REF] Johnston | Quantum Fields on Causal Sets[END_REF], these massless retarded Green's functions in the continuum allow one to derive massless discrete propagator matrices (K)

0 (x, y), such as:

in the 1 + 1-dimensional case, derived via summing over chains (hence the C subscript), or:

in the 3 + 1-dimensional case, derived via summing over paths (hence the P subscript). In the above, C 0 (x, y)

designates the causal matrix :

i.e. the matrix designating which pairs of events are connected by directed edges in the transitive closure of the causal graph (or, equivalently, which pairs of events are connected by directed paths in the full causal graph), whilst L 0 (x, y) designates the link matrix :

i.e. the matrix designating which pairs of events are connected by directed edges in the transitive reduction of the causal graph.

If one opts instead not to take the massless (i.e. m → 0) limit, then one obtains a collection of slightly more complicated expressions for the massive discrete propagator matrices (K)

M (x, y), for instance:

in the 1 + 1-dimensional/chain-summing case, or:

in the 3 + 1-dimensional/path-summing case. In the above, the matrix exponents (C 0 (x, y)) k and (L 0 (x, y)) k are defined in terms of a discrete convolution operation A * B (which reduces simply to matrix multiplication in the case of causal and link matrices):

These general forms were proposed by Dowker et al. [START_REF] Dowker | Scalar Field Green Functions on Causal Sets[END_REF] as a generic ansatz for obtaining massive causal set Green's functions in arbitrary (integer) numbers of dimensions, based on the formal expansion of massive scalar field Green's functions G

with the continuum convolution operation A * B here being given by:

Moreover, the mass parameter M here is a "normalized" version of the scalar field mass m, defined by

, where ρ c is the sprinkling density of the causal graph/causal set, i.e. the ratio ρ c = n v of the expected number of vertices in a causal subgraph n to the spacetime volume of that subgraph v in the appropriate continuum limit. As shown by Johnston [START_REF] Johnston | Particle propagators on discrete spacetime[END_REF][33], in the limit ρ c → ∞ of infinite sprinkling density, these massive discrete propagator matrices converge (modulo multiplicative constants) to the appropriate massive scalar field Green's functions:

as required. Note that, for the sake of simplicity, we have presented the above analysis for the case of flat/Minkowski space, although (following Birrell and Davies [START_REF] Birrell | Quantum Fields in Curved Space[END_REF] and Fulling [START_REF] Fulling | Aspects of Quantum Field Theory in Curved Space-Time[END_REF]) by modifying the form of the Klein-Gordon equation such that the d-dimensional Green's function G

where ζ ∈ R is a real constant, is the curved spacetime analog of the d'Alembertian operator on scalar fields φ (x):

and R (x) is the (spacetime) Ricci scalar evaluated at point x, we straightforwardly obtain the corresponding construction for massive scalar field Green's functions in arbitrary curved spacetimes (or at least those that permit Riemann normal neighborhoods). By replacing the continuum Ricci scalar R (x) with the discrete (Ollivier-Ricci) scalar curvature for Wolfram model hypergraphs/causal graphs, we obtain a corresponding construction for arbitrary Wolfram model systems. Full details may be found in [START_REF] Gorard | Axiomatic Quantum Field Theory in Discrete Spacetime via Multiway Causal Structure: The Case of Entanglement Entropies[END_REF].

Using this overall approach, we construct an exponential initial density profile for a massive scalar field Φ (t, r) defined over Wolfram model hypergraphs:

with initial density constant ρ 0 and radius of support λ. The Cauchy initial data is set up by taking a spherically-symmetric background metric in Gaussian polar coordinates: shown in Table 3. Furthermore, we perform the analogous construction (consisting of an exponential initial