Spherically symmetric black hole spacetimes on hyperboloidal slices

Gravitational radiation and some global properties of spacetimes can only be unambiguously measured at future null infinity (ℐ⁺). This motivates the interest in reaching it within simulations of coalescing compact objects, whose waveforms are extracted for gravitational wave modeling purposes. One promising method to include future null infinity in the numerical domain is the evolution on hyperboloidal slices: smooth spacelike slices that reach future null infinity. The main challenge in this approach is the treatment of the compactified asymptotic region at ℐ⁺. Evolution on a hyperboloidal slice of a spacetime including a black hole entails an extra layer of difficulty in part due to the finite coordinate distance between the black hole and future null infinity. Spherical symmetry is considered here as the simplest setup still encompassing the full complication of the treatment along the radial coordinate. First, the construction of constant-mean-curvature hyperboloidal trumpet slices for Schwarzschild and Reissner-Nordström black hole spacetimes is reviewed from the point of view of the puncture approach. Then, the framework is set for solving hyperboloidal-adapted hyperbolic gauge conditions for stationary trumpet initial data, providing solutions for two specific sets of parameters. Finally, results of testing these initial data in evolution are presented.

1. Introduction

The accurate numerical treatment of black holes (BHs) and their emitted gravitational wave (GW) signals is primordial for the field of GW astronomy. BHs are the most common participants in the compact binary coalescences observed so far [1, 2] but are challenging to model numerically due to the presence of the physical singularity inside of their horizon. GWs, as radiation propagating at the speed of light, are only unambiguously defined at future null infinity ℐ⁺, the collection of the end points of future-directed null geodesics. Future null infinity also corresponds to the idealized location of observers of astrophysical events [3–5], such as GW interferometers, where GWs signals should ideally be extracted from simulations.

Two main descriptions of BHs are common in numerically simulated spacetimes. Excision [6] involves setting an artificial timelike inner boundary inside the BH horizon to avoid the slices from reaching the physical singularity. This exploits the fact that no physical information is allowed to exit the BH, but the need to know the location of the apparent horizon at all times makes this approach technically difficult for generic spacetimes. Nevertheless, it has been successfully used to produce the largest, longest, and most accurate binary BH waveform catalog currently available [7]. In the puncture method, a specific singularity-avoidant slice of the BH spacetime is considered. This slice can have the topology of a wormhole, where the asymptotically flat end at the other side of the BH is compactified and represents the BH's location [8–10]. In the evolution of wormhole puncture initial data with the “moving puncture” gauge [11, 12], the initial slice numerically detaches from the asymptotically flat end beyond the horizon and its topology becomes close to that of a compactified trumpet [13–15], where the proper distance becomes infinite while reaching toward the symmetric point to future timelike infinity i⁺. For embedding diagrams of the wormhole and trumpet geometries, see figures 1 and 2 in Hannam et al. [15]. The construction of maximal trumpet slices has been tackled for Schwarzschild [16–20], Reissner-Nordström (RN) [21], and Kerr [22, 23]. Asymptotically, the slices considered in those studies are spacelike Cauchy and thus reach spatial infinity i⁰. The trumpet puncture approach is also chosen in the present study, for its simpler technical implementation and for the possibility to reach a portion inside of the horizon. While the latter is not required for GW extraction, it can provide insights into the numerical behavior of slices inside of the horizon, e.g., useful for the construction of Penrose diagrams of dynamical scenarios [24].

Including future null infinity within the numerical integration domain is possible by evolving on a suitable choice of foliation. The most straightforward option are characteristic slices, which can provide considerable simplifications in the equation used [25] but are prone to the development of caustics. Cauchy-characteristic matching [25–27] joins an inner Cauchy spacelike slice to an outer characteristic one along a timelike boundary. However, compatible formulations of the Einstein equations for each domain are required. In Cauchy-characteristic evolution [28–31], the same setup is used, but the Cauchy evolution is performed independently and then used as inner boundary data for the characteristic evolution. A more flexible and elegant alternative is the evolution on hyperboloidal [32–35] slices, which are spacelike and reach null infinity. An advantage that hyperboloidal evolution is expected to have and that has been achieved with Cauchy-characteristic evolution is resolving GW memory [36]. A radial compactification on hyperboloidal slices allows to include future null infinity in a finite domain. Unlike a compactification of Cauchy slices where radiation traveling out is slowed down and becomes underresolved, the outward propagation speed of signals on compactified hyperboloidal slices is of order unity and they reach ℐ⁺ at a finite coordinate time without any loss of resolution. Figure 5 in Vañó-Viñuales et al. [37] illustrates this effect with a scalar field perturbation.

Conformal compactification [38] is one method to tackle compactified hyperboloidal slices, which allow us to reach ℐ⁺ with a finite value of the coordinates. The core idea is that instead of working with the physical metric ${\tilde{g}}_{ab}$ that diverges at infinity when the coordinates are compactified, the Einstein equations are instead expressed in terms of a finite conformally rescaled metric ḡ_ab, related to the physical one by a conformal factor Ω that vanishes at ℐ⁺ at the appropriate rate

$\begin{array}{ll}{\overline{g}}_{ab}={\Omega }^2{\tilde{g}}_{ab}.& (1)\end{array}$

One of the most difficult aspects of this approach to the hyperboloidal initial value problem [35, 39] is the regularization of the resulting formally singular equations [see (2) in 2] in a way that works numerically and avoids instabilities arising from the continuum equations, in particular for hyperbolic-free evolution schemes considered here. At the analytical level, the equations were shown to be manifestly regular at ℐ⁺ [32, 35], however that specific formulation does not treat BHs in a straightforward way and suffers from continuum instabilities [40]. In contrast to the conformal approach, the dual foliation method [41, 42], a generalization of the dual coordinate frame method used in Scheel et al. [43], aims to minimize the divergent terms in the equations, making them as regular as possible. The present implementation follows Zenginoğlu's approach [44–47] to conformal compactification, using free evolution and a time-independent conformal factor Ω. Stable evolutions in spherical symmetry of regular initial data that do not form BHs were presented in Vañó-Viñuales et al. [37], while Vañó-Viñuales and Husa [48] covers experiments with suitable hyperbolic gauge conditions.

Evolving a hyperboloidal slice of spacetime including BHs is particularly challenging¹, especially in the puncture approach where both regions inside the horizons and the asymptotic far field are compactified. Constant-mean-curvature (CMC) foliations, where the trace of the physical extrinsic curvature takes a constant value, are well known in the literature, e.g., for the Schwarzschild [50–52] and RN [53, 54] spacetimes. Of special interest are those specific CMC slices that correspond to trumpet slices in their corresponding BH geometry: in a certain way, these are generalizations of the maximal trumpet slices mentioned above. The difference is that CMC slices with non-vanishing trace of the extrinsic curvature asymptotically reach null infinity and thus can be used as hyperboloidal trumpet slices suitable for evolving a BH spacetime all the way to future null infinity.

Several studies have considered hyperboloidal initial data including BHs. Configurations in spherical and in axial symmetry were presented in Schneemann [55], while Schinkel et al. [56] considered axisymmetric CMC slices for Kerr and [57] studied perturbed Kerr initial data on asymptotic CMC slices. The generalization of Bowen-York initial data to hyperboloidal slices for binaries of boosted and spinning BHs was carried out in Buchman et al. [58], whereas properties such as the Bondi-Sachs energy and momentum of the above setups were presented in Bardeen and Buchman [59]. The binary BH scenario was also studied in Schinkel [60]. However, these studies were designed with the aim to treat the BHS via excision and thus not a lot of effort was put into regularizing the slices beyond the BH horizon.

The description of BHs via punctures requires a careful treatment of the hyperboloidal slices inside the BHs as well. In a previous study [49, 61, 62], the evolution of hyperboloidal CMC Schwarzschild trumpet initial data was considered as well as the collapse into a BH of a scalar field perturbation on a regular spacetime. The trumpet dynamics was found to be highly dependent on the choice of gauge conditions. CMC trumpet initial data stationary with respect to the given gauge conditions are very desired as the evolution of any perturbation on these initial data would be easier to identify and study. Imposing stationarity is the approach suggested in Ohme et al. [63] although the slicing condition considered there is most likely not appropriate for numerical evolutions. In previous numerical experiments, a stationary solution is reached by the evolution at late times (such as that on the right of figure 2 in Vañó-Viñuales and Husa [62]), be it with BH trumpet or collapsing scalar field initial data, for at least some choices of gauge conditions. It thus makes sense to consider stationary solutions of the gauge conditions as candidates for an initial hyperboloidal trumpet slice.

The aims of this study are to review hyperbolic CMC trumpet BH initial data suitable for numerical evolutions with the puncture approach in mind (4) and to set the basic infrastructure in terms of initial data and gauge conditions to calculate stationary trumpet slices (5). An example of such a stationary configuration is solved for a specific choice of gauge, and basic evolutions for both CMC and solved-for initial data are performed on hyperboloidal slices (6). For this purpose, the already non-trivial hyperboloidal initial value problem in spherical symmetry is considered as it still contains the critical part of the radial treatment.

This article is organized as follows: In 2, the used formulation of the conformally compactified Einstein equations is briefly described, and the gauge conditions considered here are covered in 3. Initial data including a BH is treated in the following two sections: as constant-mean-curvature (CMC) in 4, while an example of solving hyperboloidal-adapted gauge conditions is provided in 5. Section 6 presents basic evolution results, and final thoughts on this study are summarized in the conclusions. The appendix collects an equation used in 5. Sections 2, 3, and 4 cover previously treated material, while sections 5 and 6 present new research.

The chosen metric signature is (−, +, +, +) and, as is customary, the fundamental constants are set to G = c = 1. The convention for the sign of the extrinsic curvature is that of Misner, Thorne, and Wheeler [64], meaning that a negative² value means expansion of the normals. Notation for the metrics is the same one as used in Vañó-Viñuales et al. [37]: the four-dimensional physical spacetime metric is denoted as $\tilde{g}$, the 4D conformal metric is denoted as ḡ, the 3D conformal spatial metric (induced by ḡ) is denoted as $\overline{\gamma }$, the 3D twice conformal metric is denoted as γ, and the 3D twice conformal background metric is denoted as $\widehat{\gamma }$.

2. Formulation

The emphasis in this study is on hyperboloidal BH initial data, so only a brief review of the formulation of the evolved system with corresponding references is given. Expressed in terms of the rescaled metric ḡ_ab as defined in (1), the 4D Einstein equations take the form as follows:

$\begin{array}{ll}{G}_{ab}[\overline{g}]=8\pi T_{ab}-\frac{2}{\Omega }\left({\overline{\nabla }}_a{\overline{\nabla }}_b\Omega -{\overline{g}}_{ab}{\overline{\nabla }}^c{\overline{\nabla }}_c\Omega \right)-\frac{3}{{\Omega }^2}{\overline{g}}_{ab}({\overline{\nabla }}_c\Omega ){\overline{\nabla }}^c\Omega .& (2)\end{array}$

The well-posed formulations considered are either the generalized BSSN [65–68] or a similar conformal version of the Z4 [69–71], the Z4c equations [72, 73]. The full derivation of these equations in terms of the conformally rescaled metric is described in Vañó-Viñuales et al. [37] and in Chapter 2 of Vañó-Viñuales [49]. The equations used in the simulations are those included in appendix C of Vañó-Viñuales et al. [37] (or appendix A in Vañó-Viñuales and Husa [48]) and again in Chapter 2 of Vañó-Viñuales [49]. There is a modification related to the evolution of BH spacetimes: a constraint damping term of the form

$\begin{array}{ll}-{\kappa }_0\frac{Z_r}{r},& (3)\end{array}$

where Z_r is a Z4 variable and κ₀ a freely specifiable parameter, can be added to ${\dot{\Lambda }}^r$'s right-hand-side (RHS). This term helps suppress instabilities if extrapolating boundary conditions are used at r = 0 (this was not necessary for evolutions of regular spacetimes as parity conditions could be imposed at the origin).

The evolution variables are the 3D conformally rescaled spatial metric:

$\begin{array}{ll}{\gamma }_{ab}=\chi {\overline{\gamma }}_{ab},& (4)\end{array}$

where ${\overline{\gamma }}_{ab}$ is the spatial metric induced from ḡ_ab, and χ is the spatial conformal factor. The conformal extrinsic curvature tensor ${\overline{K}}_{ab}$ is decomposed into its conformal trace-free part as follows:

$\begin{array}{ll}{A}_{ab}=\chi {\overline{K}}_{ab}-\frac{1}{3}{\gamma }_{ab}\overline{K},\text{ with }\overline{K}={\overline{K}}_{ab}{\overline{\gamma }}^{ab}\equiv K_{ab}{\gamma }^{ab},& (5)\end{array}$

and (in this formulation) its physical trace, mixed with the physical Z4 variable $\tilde{\Theta }$,

$\begin{array}{ll}\tilde{K}=\Omega \overline{K}-\frac{3{\beta }^a{\partial }_a\Omega }{\alpha }-2\tilde{\Theta }.& (6)\end{array}$

Evolved are A_ab, and $\tilde{K}$'s variation with respect to its initial value $\Delta \tilde{K}=\tilde{K}-{\tilde{K}}_0=\tilde{K}-K_{CMC}$ (this last parameter will be explained in 4.3). The quantity $\tilde{\Theta }$ is evolved as well if using the Z4 formulation. The Z4 variable Z_a is absorbed into the vector as follows:

$\begin{array}{ll}{\Lambda }^a={\gamma }^{bc}\left({\Gamma }_{bc}^a-{\widehat{\Gamma }}_{bc}^a\right)+2{\gamma }^{ab}{Z}_b,& (7)\end{array}$

where ${\Gamma }_{bc}^a$ are the Christoffel symbols calculated from γ_ab and ${\widehat{\Gamma }}_{bc}^a$ the ones built from a time-independent background metric ${\widehat{\gamma }}_{ab}$. The latter is chosen to be the flat spatial metric in spherical coordinates, and its explicit components (following an equivalent notation to that in (8) are given in (11). The evolved gauge variables are the conformal lapse α and the shift βⁱ.

2.1. Spherically symmetric reduction variables

The following spherically symmetric ansatz is used for the spherically symmetric line element in the conformally compactified domain (with σ² ≡ dθ² + sin² θ dϕ²):

$\begin{array}{l}d{s}^2=-\left({\alpha }^2-{\chi }^{-1}{\gamma }_{rr}{\beta }^r{}^2\right)d{t}^2\\ +{\chi }^{-1}\left[2{\gamma }_{rr}{\beta }^{r}dtdr+{\gamma }_{rr}d{r}^2+{\gamma }_{\theta \theta }{r}^{2}d{\sigma }^2\right],& (8)\end{array}$

where the freedom introduced by the spatial conformal factor χ is fixed by eliminating ${\gamma }_{\theta \theta }={\gamma }_{rr}^{-1/2}$.

In spherical symmetry, the only independent component of the trace-free part of the conformal extrinsic curvature A_ab after explicitly imposing its trace-freeness is A_rr. Only the radial component of the quantities Λ^a, β^a, and Z_a (denoted by Λ^r, β^r, and Z_r respectively) remains non-zero. The evolution variables of the spherically symmetric reduced system are χ, γ_rr, A_rr, $\Delta \tilde{K}$, Λ^r, α, β^r, and $\tilde{\Theta }$.

The conformal factor Ω is set to be a time-independent function of the compactified radial coordinate r as

$\begin{array}{ll}\Omega (r)=\left(-K_{CMC}\right)\frac{r_{ℐ}^2-r^2}{6r_{ℐ}},& (9)\end{array}$

with r_ℐ being the coordinate location of future null infinity (set to r_ℐ = 1 in the implementation without restricting generality) and K_CMC being a negative parameter described in 4.3. This expression satisfies that Ω(r) is a regular function that becomes zero at ℐ, with non-vanishing derivative there (compare e.g. [55, 74]). The origin of this expression is explained in 4.5.

3. Gauge conditions

Hyperboloidal constrained evolutions [75–77] have used suitable gauges imposed via the resolution of elliptic constraint equations. In this study, the free evolution approach is employed for its faster performance in simulations, and it requires the use of hyperbolic gauge conditions. The gauge quantities, lapse α and shift βⁱ, control the behavior of the coordinates, and they are critically important for a successful and efficient evolution. Bad choices will easily lead simulations to crash at an earlier or later time. An example of the effects of gauge choices in this hyperboloidal work is that they can induce deformations in propagating signals, as is illustrated by the (deformed) scalar field signals at ℐ⁺ in figure 2 in Vañó-Viñuales and Husa [48], that are to be corrected in post-processing.

For vanishing cosmological constant and vacuum or compact support matter sources, future null infinity is an ingoing null hypersurface. This means that no information is allowed to enter the domain from outside, making it a natural boundary for the numerical integration domain, where no boundary conditions need to be imposed—radiation just needs to be allowed to leave the spacetime. It is possible to fix ℐ⁺ to a specific coordinate location in the numerical grid (r_ℐ (9) in the present setup) for compactified hyperboloidal slices. This procedure is called scri-fixing [39, 44].

The background behavior of hyperboloidal slices differs from Cauchy ones in that the trace of the physical extrinsic curvature $\tilde{K}$ is non-zero asymptotically. This requires a modification of the usual slicing conditions commonly used in numerical simulations. See for example, the generalizations of the Bona-Massó family of slicing conditions [78] and the modifications of the Gamma-driver shift [79] and harmonic shift conditions [80] included in Vañó-Viñuales and Husa [48]. The basic idea behind those modifications is the addition of specific non-principal-part source terms to the gauge evolution equations to ensure that a hyperboloidal slice of Minkowski spacetime will be a stationary solution of the gauge equations. This is described in the next subsection.

An optimal prescription for hyperbolic gauge conditions for the conformally compactified hyperboloidal approach is still to be found. Experimentation with possible gauge source functions has provided several successful working examples. They are being further studied and extended by including a BH in the spacetime, and elsewhere by being tested in the full 3D case [81]. Work toward finding suitable gauge conditions [82] is also being tackled from the dual foliation approach.

3.1. Hyperbolic gauge conditions tested with BH spacetimes

When applying the gauge conditions discussed in Vañó-Viñuales and Husa [48] to a BH spacetime, one important aspect is to recognize that harmonic slicing is only marginally singularity avoiding, which means that a singularity is reached in an infinite coordinate time. Harmonic slicing is thus not a good choice in the neighborhood of a BH if excision is not used. However, near ℐ⁺, the physical propagation speeds of harmonic lapse (and shift) ensure that no unknown gauge information enters the numerical domain through future null infinity. Thus, the optimal scenario is to use harmonic slicing near ℐ⁺ and something different close to the BH. A condition that has provided successful evolutions using trumpet initial data is, with ${}^{\cdot }\equiv {\partial }_t$ and ${}^{\prime }\equiv {\partial }_r$ ,

$\begin{array}{ll}\dot{\alpha }={\beta }^r{\alpha }^{\prime }-{\hat{\beta }}^r{\hat{\alpha }}^{\text{'}}-\frac{(n_{cK}+{\alpha }^2)△\tilde{K}}{\Omega }+\frac{{\Omega }^{\prime }}{\Omega }({\hat{\beta }}^r\hat{\alpha }-{\beta }^r\alpha )+\frac{{\xi }_{cK}(\hat{\alpha }-\alpha )}{\Omega },& (10)\end{array}$

where ξ_cK is a parameter used to damp the behavior of the lapse at ℐ⁺. This equation is equivalent to (20) in Vañó-Viñuales and Husa [48] with ${\xi }_1\widehat{\alpha }={\xi }_{cK}$, ξ₂ = 0, and ${\alpha }^{2}f=n_{cK}+{\alpha }^2$, later setting n_cK³ to be proportional to Ω. Note that the coefficient in front of $\Delta \tilde{K}$ is similar to the shock-avoiding slicing condition [20, 83, 84]. This form was chosen for the following considerations. The α² part provides physical propagation speeds for the gauge modes (the first three lines listed in Figure 1), as mentioned above. This is desired at ℐ⁺ because then all propagation speeds are either positive or zero, and there are no incoming modes there. However, near the location of the trumpet inside of the BH's horizon, the physical propagation speeds become zero (as α = 0 at the location of the trumpet). The effect is that any signals that have entered the BH region and travel along the infinitely long cylinder of the trumpet slice will propagate slower and slower, soon becoming underresolved, which can lead to numerical instabilities. Increasing the gauge propagation speeds allows perturbations to leave the domain in a finite time and provides more stable evolutions in general and also gives smoother stationary values for the evolution quantities at the trumpet. Examples of modified propagation speeds for the lapse and shift conditions are shown in Figure 1.

FIGURE 1

Figure 1. Zero speed $c_0=-{\beta }^r$ and incoming (−) and outgoing (+) lightspeeds $c_{\pm }=-{\beta }^r\pm \alpha \sqrt{\frac{\chi }{{\gamma }_{rr}}}$ plotted for a CMC slice of a Schwarzschild BH with M = 1, using K_CMC = −1 and r_ℐ = 1. The vertical line indicates the radial position of the horizon, where c₀, c₋ < 0 and c₊ = 0. At ℐ⁺, we have c₀, c₊ > 0 and c₋ = 0. Except for c₀, outgoing speeds are shown in black and ingoing ones in blue. The two sets of dashed lines correspond to the incoming and outgoing modified propagation speeds related to the slicing condition as in (10) with the choices of parameter m_cK used in (10). While the speeds go to zero at the location of the trumpet (r = 0) and coincide whith the lightspeeds at ℐ⁺, their values are different in the rest of the domain. The two sets of dash-dot curves show the characteristic speeds associated with the shift condition with λ = 1, for shift choices (12) and (13). In the second case, where the shift advection terms have been dropped, the incoming speed is made to be zero at ℐ⁺, but that forced the outgoing one to be equal to c₀ (instead of c₊) there compared to a Minkowski equivalent (with different parameter choices) in figure 1 in Vañó-Viñuales and Husa [48].

These gauge source functions are designed to make a hyperboloidal CMC slice of Minkowski (encoded in the hatted quantities), a stationary solution of the slicing equation: $\stackrel{.}{\alpha }=0\leftrightarrow \alpha =\widehat{\alpha }$, ${\beta }^r={\hat{\beta }}^{{}^r}$, $\Delta \tilde{K}=0$. The components of the background conformally compactified metric (following an ansatz like that of (8)) that appear in (10), and are used to calculate ${\widehat{\Gamma }}_{bc}^a$ in (7), are

$\begin{array}{l}\widehat{\chi }={\hat{\gamma }}_{rr}={\hat{\gamma }}_{\theta \theta }=1,\widehat{\alpha }=\sqrt{{\Omega }^2+{\left(\frac{K_{CMC}r}{3}\right)}^2}\text{ and}\\ {\widehat{\beta }}^r=\frac{K_{CMC}r}{3}.& (11)\end{array}$

They are obtained (4.3) from (21) or (24) setting $A(\frac{r}{\overline{\Omega }})=1$, M = 0, C_CMC = 0, and $\overline{\Omega }=\Omega$.

For the shift condition, two different options are considered. One is a variant of the integrated Gamma-driver [79] adapted to hyperboloidal slices:

$\begin{array}{l}\stackrel{.}{{\beta }^r}={\beta }^r{\beta }^r{}^{\prime }-{\hat{\beta }}^r{\hat{\beta }}^r{}^{\prime }+\left(\lambda (r_{ℐ}^2-r^2)+\frac{3}{4}{\alpha }^2\chi \right){\Lambda }^r+\eta ({\hat{\beta }}^r-{\beta }^r)\\ +{\xi }_{{\beta }^r}\left(\frac{{\hat{\beta }}^r}{\Omega }-\frac{{\beta }^r}{\Omega }\right),& (12)\end{array}$

mostly the same as (26) in Vañó-Viñuales and Husa [48]. The coefficient in front of Λ^r is chosen in such a way that the associated propagation speeds will be the physical ones at ℐ⁺. The positive parameter λ will only increase the speeds near the trumpet in a similar fashion as n_cK for the slicing condition above. This is shown in Figure 1. The other shift option is to have an expression purely proportional to Λ^r: the resulting system is still hyperbolic, and it will have conformally flat initial data as a stationary solution (more on this in 5). However, dropping the advection terms modifies the characteristic propagation associated to the shift condition. In order to ensure that the related ingoing speed at ℐ⁺ is still zero, the coefficient in front of Λ^r is modified as

$\begin{array}{ll}\stackrel{.}{{\beta }^r}=\left(\lambda (r_{ℐ}^2-r^2)+3{\alpha }^2\chi +\frac{3}{2}{\gamma }_{rr}{\beta }^r{}^2+\frac{9}{2}\sqrt{{\gamma }_{rr}\chi }\alpha {\beta }^r\right){\Lambda }^r.& (13)\end{array}$

The resulting outgoing propagation speed is also modified: it is smaller (although still positive) at ℐ⁺ (see Figure 1), and it would be positive even inside of the horizon if the λ term was not present. The choice of the coefficient in (13) giving zero ingoing speed is not unique, but it has been used here for its good behavior in numerical evolutions, both at the origin and at ℐ⁺.

Hyperboloidal CMC trumpet initial data [(24) and (25) as derived in 4, or any initial data satisfying the relations in 5.2] are a stationary solution of the Einstein equations as described in 2. However, if they are evolved together with gauge conditions whose source functions are filled with hyperboloidal CMC Minkowski data (11) as described above, the right-hand-sides of the gauge evolution equations will not be zero. Thus, some gauge dynamics will take place in which the trumpet slice readjusts and settles into a new stationary solution, see the plot on the right of figure 2 in Vañó-Viñuales and Husa [62]. The change in the slices is easier to understand when depicted as a Carter-Penrose diagram, as in Figure 11B. While this scenario is satisfactory in the sense that a long-term solution is found, the initial dynamics does not allow to decouple any potential perturbations of the system from the trumpet dynamics. Naively, a way to try to obtain the desired outcome—trumpet initial data that are a stationary solution of the gauge conditions—is to fill in the gauge source terms in the gauge conditions with (24), the same data as the one given initially. This has been tested (see section 8.2 in Vañó-Viñuales [49]), with the result that a slow exponential growth appeared in the evolutions, causing the simulations to a crash in finite time. The conclusion of these tests is that the chosen trumpet initial data are a stationary but not a stable solution for the gauge conditions with trumpet source terms (more on this in 5). However, the growth in these simulations is slow enough to study small scalar field perturbations, as presented in Vañó-Viñuales and Husa [61]. Whether a different choice of trumpet slice or form of the gauge conditions would not cause the growth is an open question.⁴ Meanwhile, an attempt to combine stability and stationarity together is described in 5, where a solution for the gauge conditions with hyperboloidal Minkowski source functions is determined for a specific setup.

4. Constant-mean-curvature initial data

4.1. Main ingredients of hyperboloidal conformal compactification

At the core of the hyperboloidal approach is the foliation of spacetime along hyperboloidal slices, which can be characterized as the level sets of a specific parameter. This parameter is taken to be the hyperboloidal time coordinate t, and it is related to the usual time coordinate $\tilde{t}$ via the height function $h(\tilde{r})$ [85, 86] as

$\begin{array}{ll}t=\tilde{t}-h(\tilde{r}).& (14)\end{array}$

The height function satisfies $dh/d\tilde{r}<1$ everywhere except asymptotically, where $dh/d\tilde{r}{|}_{ℐ}=1$ holds, thus characterizing the hyperboloidal slices as spacelike but extending to ℐ⁺.

In order to reach future null infinity with a finite value of the spatial coordinates, the radial coordinate $\tilde{r}$ on a hyperboloidal slice is compactified into a new r using a compactification factor $\overline{\Omega }(r)$

$\begin{array}{ll}\tilde{r}=\frac{r}{\overline{\Omega }(r)}.& (15)\end{array}$

Following (1), the line element is conformally rescaled by the conformal factor Ω, to provide regular metric components all the way to ℐ⁺

$\begin{array}{ll}d{s}^2={\Omega }^{2}d{\tilde{s}}^2.& (16)\end{array}$

The compactification factor $\overline{\Omega }$ is not to be confused with the conformal factor Ω, as they are a priori different quantities. While the conformal compactification method relies in both having the same (or at least proportional) behavior near ℐ⁺, their behavior in other parts of the domain (especially at the location of the BHs) can be chosen to be very different. For the spherically symmetric data considered here, an example of this is illustrated in Figure 5.

4.2. Spherically symmetric conformally compactified hyperboloidal slices

A suitable starting point to derive spherically symmetric vacuum initial data on a hyperboloidal slice is the following line element on an uncompactified Cauchy slice:

$\begin{array}{llll}d{\tilde{s}}^2& =& -A(\tilde{r})d{\tilde{t}}^2+\frac{1}{A(\tilde{r})}d{\tilde{r}}^2+{\tilde{r}}^{2}d{\sigma }^2& (\mathrm{17a})\end{array}$

$\begin{array}{l}=-A(\tilde{r})d{t}^2-2A(\tilde{r})h^{\prime }(\tilde{r})dtd\tilde{r}+\frac{1-A(\tilde{r}){\left(h^{\prime }(\tilde{r})\right)}^2}{A(\tilde{r})}d{\tilde{r}}^2\\ +{\tilde{r}}^{2}d{\sigma }^2,& (\mathrm{17b})\end{array}$

expressed first in terms of the usual time $\tilde{t}$ and then in terms of the hyperboloidal time coordinate t after using (14). This ansatz for the initial metric is general enough to consider flat spacetime, the Schwarzschild and Reissner-Nordström (RN) spacetimes, and the addition of a non-vanishing cosmological constant. After applying the radial compactification (15) and conformal rescaling (16), it becomes

$\begin{array}{l}d{s}^2=-A{\Omega }^{2}d{t}^2+\frac{{\Omega }^2}{{\overline{\Omega }}^2}[-2A{h}^{\prime }(\overline{\Omega }-r{\overline{\Omega }}^{\prime })dtdr\\ +\frac{\left[1-{\left(A{h}^{\prime }\right)}^2\right]}{A}\frac{{(\overline{\Omega }-r{\overline{\Omega }}^{\prime })}^2}{{\overline{\Omega }}^2}d{r}^2+r^{2}d{\sigma }^2],& (18)\end{array}$

where A and h′ are functions of $\frac{r}{\overline{\Omega }}$, while $\overline{\Omega }$ and Ω depend on r.

4.3. Hyperboloidal CMC slices

A convenient way of slicing spacetime is doing it along constant-mean-curvature (CMC) slices, on which the trace of the physical extrinsic curvature ($\tilde{K}$) is a constant. A special case of CMC slices are maximal slices [87], where $\tilde{K}=0$ and spatial infinity is reached asymptotically. Maximal Schwarzschild trumpet slices are analytically described in Baumgarte and Naculich [16]. Generalizations for slices with a non-vanishing $\tilde{K}$ are considered in Refs. [85, 86, 88, 89] and including the critical case for trumpets in Malec and O'Murchadha [50] and Buchman et al. [58], while a study of CMC slices in the RN geometry has been performed in Tuite and O'Murchadha [53].

This derivation of a height function providing CMC slices follows [85, 86], see subsection 3.2.2 in Vañó-Viñuales [49] for a more detailed derivation. The basic procedure is to express the unit normal ñ^a to the hypersurface in terms of the metric (17b) and use it to calculate the expression for the trace of the physical extrinsic curvature:

$\begin{array}{ll}\tilde{K}=-\frac{1}{\sqrt{-\tilde{g}}}{\partial }_a\left(\sqrt{-\tilde{g}}{ñ}^a\right)=-\frac{1}{r^2}{\partial }_r\left[r^2\frac{A^{3/2}(\tilde{r})h^{\prime }(\tilde{r})}{\sqrt{1-{(A(\tilde{r})h^{\prime }(\tilde{r}))}^2}}\right].& (19)\end{array}$

Setting it equal to a constant value of $\tilde{K}=K_{CMC}<0$ and introducing C_CMC as an integration constant, the first derivative of the height function is isolated to give

$\begin{array}{ll}{h}^{\prime }(\tilde{r})=-\frac{\frac{K_{CMC}\tilde{r}}{3}+\frac{C_{CMC}}{{\tilde{r}}^2}}{A(\tilde{r})\sqrt{A(\tilde{r})+{(\frac{K_{CMC}\tilde{r}}{3}+\frac{C_{CMC}}{{\tilde{r}}^2})}^2}}.& (20)\end{array}$

The expression for the flat spacetime case is obtained by setting $A(\tilde{r})=1$ and C_CMC = 0, and the height function can be integrated to $h(\tilde{r})=\sqrt{{(3/K_{CMC})}^2+{\tilde{r}}^2}$.

Comparing our line element of initial data (18) with our metric ansatz (8) and substituting (20), we assign the following initial values to our metric components, where the notation X₀ ≡ X(t = 0) is used:

$\begin{array}{ll}{\gamma }_{\theta \theta 0}=1,{\chi }_0=\frac{{\overline{\Omega }}^2}{{\Omega }^2},{\gamma }_{rr0}=\frac{{(\overline{\Omega }-r{\overline{\Omega }}^{\prime })}^2}{{\tilde{\alpha }}_0^2{\overline{\Omega }}^2},{\alpha }_0=\Omega {\tilde{\alpha }}_0,& (\mathrm{21a})\end{array}$

$\begin{array}{l}{\beta }_0^r=\frac{\left(\frac{K_{CMC}r}{3\overline{\Omega }}+\frac{C_{CMC}{\overline{\Omega }}^2}{r^2}\right){\tilde{\alpha }}_0{\overline{\Omega }}^2}{(\overline{\Omega }-r{\overline{\Omega }}^{\prime })},\text{ with}\\ {\tilde{\alpha }}_0=\sqrt{A(\frac{r}{\overline{\Omega }})+{\left(\frac{K_{CMC}r}{3\overline{\Omega }}+\frac{C_{CMC}{\overline{\Omega }}^2}{r^2}\right)}^2}.& (\mathrm{21b})\end{array}$

A height function determined by imposing CMC is not the only suitable choice. It could also be only asymptotically CMC [57]. Other possible options, e.g., [42], where $h^{\prime }(\tilde{r})$ is chosen to provide unit outgoing radial coordinate lightspeed, and [90], where $h(\tilde{r})$ is introduced as part of the minimal gauge [91].

4.4. Hyperboloidal CMC trumpet slices

The choice of the integration constant C_CMC is a relevant matter as a critical value exists that provides trumpet [14] CMC data [58] in an equivalent way as done in the non-hyperboloidal case [16]. This critical value of C_CMC depends on M, Q, and K_CMC and is calculated [16, 87] by setting to zero the discriminant of the denominator (6th order polynomial) of ⁵

$\begin{array}{ll}{\gamma }_{\tilde{r}\tilde{r}0}(\tilde{r})=\frac{1}{A(\tilde{r})+{\left(\frac{K_{CMC}\tilde{r}}{3}+\frac{C_{CMC}}{{\tilde{r}}^2}\right)}^2}.& (22)\end{array}$

With this critical choice of C_CMC (see (3.42) in Vañó-Viñuales [49] for the explicit expression for RN), the denominator now has a double real root at $\tilde{r}=R_0$. This finite value of the radial coordinate $\tilde{r}$ is where the slice (reaching ℐ⁺ in its outer end) finishes, corresponding to the location of the trumpet [58]. However, in terms of proper distance, the inner end of the slice is infinitely far away from the singularity. For instance, in the Schwarzschild case and for maximal K_CMC = 0, the double root is R₀ = 3M/2 [16, 87], whereas for K_CMC → −∞, it tends to R₀ → 2M. The dependence of the double root R₀ on the charge Q and the value of K_CMC is shown in Figure 2. Note that in the extreme Reissner-Nordström case (Q = M), R₀/M is always unity.

FIGURE 2

Figure 2. Innermost value of the Schwarzschild-like radial coordinate (the double root R₀) reached by the outer CMC trumpet slice, as a function of K_CMC and Q (for zero cosmological constant). Taken from Vañó-Viñuales [49].

The effect of C_CMC's value on the CMC slices is depicted in figure 1 in Buchman et al. [58] and in Tuite and O'Murchadha [53] and illustrated in Figure 3 in the form of Carter-Penrose diagrams of the Schwarzschild spacetime. For a value of C_CMC smaller than the critical one, the denominator of (22) has two different real roots (R₁, R₂) for $\tilde{r}$, the outer one R₂ corresponding to the location of the minimal surfaces mentioned in Buchman et al. [58]. If $C_{CMC}<-\frac{1}{3}{K}_{CMC}{\left(M+\sqrt{M^2-Q^2}\right)}^3$ (such as the example shown in Figure 3A), the slices reach inside of the white hole, while for a larger value of C_CMC (Figure 3B), they enter the BH. Quantities become complex for $\tilde{r}\in (R_1,R_2)$; the corresponding part of the diagrams is left in white. As mentioned above, for the critical value of C_CMC, a double root appears and complete CMC trumpet slices (Figure 3C) exist, joining either ℐ⁺ to the symmetric point to future timelike infinity (i⁺) (the outer slices) or the singularity to i⁺ (the inner ones). For C_CMC larger than the critical one (Figure 3D), there is no root to the polynomial in (22) and the CMC slices extend between null infinity and the singularity. Examples of the outer CMC Schwarzschild trumpet slices for critical C_CMC for different values of K_CMC is given in figure 3.5 in Vañó-Viñuales [49], showing the maximal case $\tilde{K}=0$ corresponding to the usual trumpet slices [16] in the first subfigure. For a positive value of K_CMC (in the current sign convention), the hyperboloidal slices reach past null infinity ℐ⁻ instead of ℐ⁺. Equivalent penrose diagrams depicting CMC slices for the RN spacetime (with $A(\tilde{r})=1-\frac{2M}{\tilde{r}}+\frac{Q^2}{{\tilde{r}}^2}$) are shown in figures 3.11, 3.12, 3.13, and 3.14 in Vañó-Viñuales [49]. The non-extremal case is illustrated by the choice Q = 0.9M. The trumpet case (with critical C_CMC) should compare to panel 2 in figure 2 in Tuite and O'Murchadha [53]; the slices have a different profile because theirs were probably not numerically determined. The extremal case Q = M is shown in 3.14 in Vañó-Viñuales [49], and it has the feature that all slices are trumpet ones with R₀ = M always, as can also be seen in Figure 2. For the over-extreme case (Q > M), no critical value of C_CMC is found; thus, no trumpet slices can be constructed (at least with this method).

FIGURE 3

Figure 3. Carter-Penrose diagrams representing hyperboloidal CMC Schwarzschild foliations with M = 1 and K_CMC = −1 and the choices of C_CMC of the height functions depicted in Figure 4. They correct the versions included in figures 3.3 and 3.4 in Vañó-Viñuales [49]. (A) Incomplete slices entering the white hole. (B) Incomplete slices entering BH. (C) Trumpet slices. (D) Slices reaching the singularity. In figures (A, B), the height function is imaginary between R₁ and R₂. For the critical C_CMC (C), the slices that reach the singularity and those that reach ℐ⁺ are separated by a thick line located at $\tilde{r}=R_0$ that corresponds to the double root. In the C_CMC = 4 case on (D), the hyperboloidal slices go all the way from ℐ⁺ into the singularity. This last diagram is qualitatively the same (with different values of the parameters) as the Penrose diagram in figure 10 in Zenginoğlu [44]. The present numerical experiments use the outer slices (those reaching ℐ⁺) of (C) as the critical value of C_CMC allows to map the trumpet slice to the whole of the radial coordinate range r ∈ (0, r_ℐ].

The effect of C_CMC can also be seen in the Schwarzschild examples of the height function $h(\tilde{r})$ shown in Figure 4. They were used to construct the respective Penrose diagrams in Figure 3. The height function, integrated numerically from (20), if expressed in Schwarzschild coordinates has a coordinate singularity at the location of the horizon ($\tilde{r}=2M$): it diverges downward for C_CMC = 2 because it enters the white hole and upwards for the cases crossing the BH horizon. The region between (R₁, R₂) is complex for the subcritical values of C_CMC (absence of curves), so it is not straightforward how to join the inner and outer real parts of the corresponding height functions. For the critical case, the outer part of the height function goes to −∞ at the root $\tilde{r}=R_0=1.905$. For the supercritical C_CMC = 4, the height function attains a finite value at the singularity $\tilde{r}=0$. The integration constants for each case have been set in such a way that as $\tilde{r}\to \infty$, the Schwarzschild height functions approach the flat spacetime one $h(\tilde{r})=\sqrt{{(3/K_{CMC})}^2+{\tilde{r}}^2}$.

FIGURE 4

Figure 4. Example of numerically integrated height functions for the Schwarzschild spacetime with M = 1 and K_CMC = −1. The values of C_CMC correspond to those used in Figure 3. The integration constant has been set in such a way that as $\tilde{r}\to \infty$, the Schwarzschild height function approaches the flat spacetime one ($h(\tilde{r})=\sqrt{{(3/K_{CMC})}^2+{\tilde{r}}^2}$). In blue the parts of the height function that are imaginary (for C_CMC, smaller than the critical one).

Initial data developed in refs. [56–58] aims for its evolution using excision. For that purpose, any CMC slices with $C_{CMC}>-\frac{1}{3}{K}_{CMC}{\left(M+\sqrt{M^2-Q^2}\right)}^3$ (intersecting the BH horizon), are probably suitable, as in any case they will be cut before reaching the singularity. The puncture approach pursued here, however, compactifies the full slice, so the most suitable choice is the CMC trumpet slice for the critical value of C_CMC. This is the case that will be considered from now on.

4.5. Compactification $\overline{\Omega }$ and conformal rescaling Ω for CMC slices

It is convenient to determine the compactification factor $\overline{\Omega }$ by imposing a conformally flat initial spatial metric, which is in a way equivalent to transforming to the isotropic radial coordinate:

$\begin{array}{ll}{\gamma }_{rr0}=\frac{{(\overline{\Omega }-r{\overline{\Omega }}^{\prime })}^2}{{\overline{\Omega }}^2\left[A(\frac{r}{\overline{\Omega }})+{\left(\frac{K_{CMC}r}{3\overline{\Omega }}+\frac{C_{CMC}{\overline{\Omega }}^2}{r^2}\right)}^2\right]}\equiv 1,& (23)\end{array}$

as this is also a simple choice compatible with an initial zero (32). The factor $\overline{\Omega }$ is expected to vanish at the same rate as the conformal factor Ω at ℐ⁺, but it is allowed to have a different behavior elsewhere. Expression (23) is solved numerically for $\overline{\Omega }$; a suitable procedure is described in subsection 6.6.1 of Vañó-Viñuales [49]. The result for the Schwarzschild case is shown as the solid line in Figure 5, with the dotted line representing the linear behavior of $\overline{\Omega }$ near r = 0, corresponding to the asymptotic behavior near the trumpet $\tilde{r}=R_0$ (see (5) in Baumgarte and Naculich [16]). The resulting $\overline{\Omega }$ will only extend to r = 0 if the critical value of C_CMC is used. If C_CMC is larger than the critical value, $\overline{\Omega }$ will diverge as r → 0, while using a value smaller than the critical C_CMC will provide a $\overline{\Omega }$ that does not reach the origin (and does not vanish at the smaller r it gets to).

FIGURE 5

Figure 5. Conformal (Ω) and compactification ($\overline{\Omega }$) factors, as well as the behavior of $\overline{\Omega }~r/R_0$ near the origin (with R₀ the location of the trumpet in uncompactified Schwarzschild radial coordinate), for K_CMC = −1 and critical C_CMC. This figure is also included in Vañó-Viñuales and Husa [61].

In flat spacetime ($A(\frac{r}{\overline{\Omega }})=1$), condition (23) can be solved analytically resulting in expression (9), where $\overline{\Omega }$ is substituted by Ω. This result is commonly seen in the literature [55, 74] and has been used in a preceding study [37, 48] as conformal and compactification factors for the flat spacetime case. In order to ensure that the slices reach ℐ⁺ (instead of ending at spacelike infinity i⁰), (9) satisfies that Ω|_ℐ = 0 and ∇_aΩ|_ℐ ≠ 0. Then, its behavior at ℐ⁺ is unaffected by the choices of the parameters M and C_CMC, and it is well-behaved at the origin of the coordinate system. It is thus a good candidate to be used as a time-independent conformal factor Ω in general, so this is indeed also the choice in this study dealing with BH spacetimes, as already mentioned toward the end of 2.1. The profile of this choice of Ω appears in Figure 5 as a dashed line.

A comparison between compactified spherically symmetric BHs for various values of the charge Q is shown in Figure 6. In Figure 6A, the CMC trumpet compactification factors $\overline{\Omega }$ corresponding to Schwarzschild, to RN with Q = 0.9M and to extreme RN (Q = M) are presented. An interesting effect of the extremality of the Q = M case is that the cylindrical infinity of the trumpet and the BH horizon are mapped to the same point r = 0 of the isotropic radius r. This can be easily recognized by looking at the profiles of the CMC trumpet initial values of the shift β^r, displayed on the right in Figure 6B. In the Schwarzschild and non-extreme RN cases, the shift is positive at the horizon (mapped to r_Schw ≈ 0.13 and r_RN+ ≈ 0.071 for Q = 0.9M respectively), but in the extreme case, the shift never becomes positive. The curves corresponding to the Schwarzschild case of the shift β^r in Figure 6B (see plot on the left figure 3.10 in Vañó-Viñuales [49] for a representation of the lapse) compared to figure 5 in Hannamet al. [13] and figure 2 in Baumgarte and Naculich [16], with the difference that here the data are compactified on a hyperboloidal slice instead of a spacelike one (with vanishing $\tilde{K}$).

FIGURE 6

Figure 6. (A) CMC trumpet compactification factor. (B) CMC trumpet values for shift β^r. A version of these figures was originally included in Vañó-Viñuales [49]. (A) Profiles of the compactification factor $\overline{\Omega }$ with K_CMC = −1 and critical C_CMC for Schwarzschild, RN with Q = 0.9M and extreme RN CMC trumpet geometries. (B) Trumpet values for the shift β^r for Schwarzschild, RN with Q = 0.9M and extreme RN geometries. The detail in β^r's plot shows how the shift in the extreme RN case is never positive, a consequence of the trumpet and the horizon being mapped to the same point of the compactified radial coordinate.

4.6. Vacuum initial data: Schwarzschild spacetime

From now on, the described CMC trumpet BH initial data, suitable for evolutions using the formalism in 2, will be restricted to the Schwarzschild ($A(\tilde{r})=1-\frac{2M}{\tilde{r}}$) spacetime. The Reissner-Nordström case works equivalently. The consideration of a non-vanishing cosmological constant is left for future studies. After imposing conformal flatness (23), ${\overline{\Omega }}^{\prime }$ can be isolated from there and introduced into (21), yielding for the metric components

$\begin{array}{ll}{\chi }_0=\frac{{\overline{\Omega }}^2}{{\Omega }^2},{\gamma }_{rr0}={\gamma }_{\theta \theta 0}=1,& (\mathrm{24a})\end{array}$

$\begin{array}{l}{\alpha }_0=\Omega \sqrt{\left(1-\frac{2M\overline{\Omega }}{r}\right)+{\left(\frac{K_{CMC}r}{3\overline{\Omega }}+\frac{C_{CMC}{\overline{\Omega }}^2}{r^2}\right)}^2},\\ {\beta }_0^r=\frac{K_{CMC}r}{3}+\frac{C_{CMC}{\overline{\Omega }}^3}{r^2},& (\mathrm{24b})\end{array}$

and for the extrinsic curvature ($\tilde{K}=K_{CMC}$) and rest of the variables [calculated from (30), (31), and (32) using (24)]

$\begin{array}{ll}{\Lambda }_0^r={\tilde{\Theta }}_0=Z_{r0}=0,A_{rr0}=-\frac{2C_{CMC}{\overline{\Omega }}^3}{r^3\Omega }\text{ and }△{\tilde{K}}_0=0.& (25)\end{array}$

The background metric ${\widehat{\gamma }}_{ab}$ is chosen to be the initial value of the evolved one, that is, conformally flat ${\widehat{\gamma }}_{rr}={\widehat{\gamma }}_{\theta \theta }=1$ as indicated in (11). The profiles of the initial values of the variables (24) and (25) are shown in Figure 7. As mentioned in 3.1, CMC initial data (24) and (25) are a stationary solution of the Einstein equations but not of the gauge conditions described in 3.

FIGURE 7

Figure 7. CMC Schwarzschild trumpet BH initial data for the evolution variables with K_CMC = −1, M = 1, and critical C_CMC. The vertical line denotes the location of the horizon. As is expected for trumpet data, the lapse α and spatial conformal factor χ are zero at the puncture. The shift component β^r is positive at the horizon and negative at ℐ⁺. Compare to solved-for stationary data shown in Figure 10.

5. Stable stationary initial data for given gauge conditions

In this study, “stationary” initial data means data that correspond to a stationary solution of the Einstein equations and also of the chosen gauge conditions. Only if all RHSs (right-hand sides) of the evolution equations are zero for the initial data, the solution will be stationary. This is valid for both physical and gauge dynamics. In 4, an example was given describing how CMC initial data were a stationary solution of the gauge conditions with CMC-constructed source terms. However, the evolution would start diverging from those data as soon as the simulation started, and variable values would grow exponentially, leading the simulation to crash. This stationary solution was unstable under the small discretization errors naturally arising in a numerical code.

In this same context, “stable” describes an “attractor-type” solution to the system. As explained in 4, gauge source functions calculated from CMC Minkowski data give a stable stationary end state (after some trumpet gauge dynamics) for some choices of gauge conditions and parameters. Initial and final trumpet states of an instance of that evolution are included in figure 2 in Vañó-Viñuales and Husa [62]), where the latter also coincides with the final state of a collapsed scalar field creating a BH with the same total mass⁶ (for the same gauge configuration used). Ideally, we want stable stationary hyperboloidal trumpet initial data, which will remain a solution of the system even when small initial perturbations are present.

The way to set up initial data for a given spherically symmetric Schwarzschild BH spacetime (for some chosen values of M and K_CMC) is to impose the conditions (listed in 5.2) that satisfy the Einstein equations in the non-dynamical regime and then solve the gauge conditions for stationarity. The latter means setting $\stackrel{.}{\alpha }=0$ and ${\dot{\beta }}^r=0$, which correspond, respectively, to the choice of hyperboloidal trumpet slicing and the compactification factor and solving for the remaining degrees of freedom. There are several options to tackle the last part:

• Solve first the slicing condition on the uncompactified domain (using the uncompactified radial coordinate $\tilde{r}$ or equivalent) and afterwards use the shift condition to determine a suitable compactification for the radial coordinate, in the form of (15). The advantages of this approach are that the steps are performed separately and only one equation is to be solved at a time. The disadvantages are that the first integration is to be performed up to infinite values of the uncompactified radial coordinate, which will introduce considerable errors near ℐ⁺ unless some type of compactification is performed, and it also requires determining the location of the trumpet, which is not trivial. This optional partial compactification is difficult to deal with in the second step (imposing stationarity on the shift condition to obtain the compactification) as the solving procedure has to be built on top of it consistently – all of this assuming that a solution to the equation indeed exists.

• Solve both slicing and shift conditions [e.g., in the forms (10) and (12)] for stationarity at the same time. This would in principle allow us great freedom in the gauge conditions that we choose to solve for, provided they give an existing final stable stationary solution for a trumpet after some gauge dynamics. However, there has been no success so far despite numerous attempts. The main difficulty seems to lie within the form of the shift condition. The stable stationary trumpet state reached at late times in experiments sometimes shows a non-smooth profile of the field Λ^r at ℐ⁺ (see fourth panels in figure 8.31 in Vañó-Viñuales [49]). Relation (32) needs to hold in the stationary regime and maybe that condition is incompatible with the presence of advection terms in the shift condition (12). In general, the final stable stationary state is very sensitive to the choice of gauge conditions.

• Solve the slicing condition for the trumpet geometry and impose an initial conformally flat metric to obtain the compactification, in an equivalent way as done for CMC data with (23). The advantage is that both equations can be solved at the same time (the compactification allowing to solve all the way to ℐ⁺). The slicing condition [here, (10) is considered] is to be chosen carefully as it plays an important role in the trumpet geometry [19]. The disadvantage of this approach is that the shift condition needs to be modified in order to keep the obtained initial data stationary, namely, either dropping the advection and “η” terms or having its gauge source functions filled in by the stationary solution found.

The three options have been attempted, with only some success for some specific cases in the third way of proceeding. This best choice will be described in 5.3 and solved for two example configurations.

5.1. Comparison of metric quantities in physical and conformal domains

For the purpose of clarifying the relations between physical and conformally compactified quantities in the metric, let us introduce the following ansatz for the spherically symmetric line element in the uncompactified physical domain in terms of the hyperboloidal time t (again with dσ² ≡ dθ² + sin² θ dϕ²):

$\begin{array}{ll}d{\tilde{s}}^2=-\left({\tilde{\alpha }}^2-{\tilde{X}}_{\tilde{r}\tilde{r}}{{\beta }^{\tilde{r}}}^2\right)d{t}^2+2{\tilde{X}}_{\tilde{r}\tilde{r}}{\beta }^{\tilde{r}}dtd\tilde{r}+{\tilde{X}}_{\tilde{r}\tilde{r}}d{\tilde{r}}^2+{\tilde{X}}_{\theta \theta }{\tilde{r}}^{2}d{\sigma }^2.& (26)\end{array}$

The values of the metric components for a hyperboloidal slice can be read off by comparing this metric ansatz to (17b) [or for a Cauchy slice if using $\tilde{t}$ instead and relating to (17a)].

For the conformally compactified version, as will be used in initial data calculations in 5.3 and relates to the physical one as in (16), set

$\begin{array}{l}d{s}^2={\Omega }^{2}d{\tilde{s}}^2\\ =-\left({\alpha }^2-X_{rr}{\beta }^r{}^2\right)d{t}^2+2X_{rr}{\beta }^{r}dtdr+X_{rr}d{r}^2+X_{\theta \theta }{r}^{2}d{\sigma }^2.& (27)\end{array}$

It relates to the line element (8) used in the evolution formalism by X_rr ≡ γ_rr/χ and X_θθ ≡ γ_θθ/χ, but for convenience and clarity, it is written in terms of the X quantities. Its components can be read off from (18), as done for (8) in (21) after substitution of (20).

The relations between the physical and conformally rescaled metric quantities are the following (most listed in (2.39) and (2.69) in Vañó-Viñuales [49])

$\begin{array}{l}\tilde{\alpha }=\frac{\alpha }{\Omega },\tilde{\chi }={\Omega }^2\overline{\chi },\text{ with }\tilde{\chi }=\frac{{\tilde{\gamma }}_{\theta \theta }}{{\tilde{X}}_{\theta \theta }}\text{ equivalently as in }(8)\text{ and}\\ \overline{\chi }=\frac{\chi }{{\overline{\Omega }}^2},\text{ so that }\tilde{\chi }=\frac{{\Omega }^2}{{\overline{\Omega }}^2}\chi .& (28)\end{array}$

The quantity $\overline{\chi }$ accounts purely for the conformal rescaling (1) on the spatial conformal factor. However, χ is conformally rescaled and also includes the compactification of the radial coordinate (15).

The shift does not change due to the 4D conformal rescaling, but its radial components change under a transformation in the radial coordinate. The changes for the X metric components include both effects

$\begin{array}{l}{\beta }^{\tilde{r}}=\left(\frac{\overline{\Omega }-r{\overline{\Omega }}^{\prime }}{{\overline{\Omega }}^2}\right){\beta }^r,{\tilde{X}}_{\tilde{r}\tilde{r}}=\frac{{\overline{\Omega }}^2}{{\Omega }^2}{\left(\frac{\overline{\Omega }}{\overline{\Omega }-r{\overline{\Omega }}^{\prime }}\right)}^2{X}_{rr},\text{ and }\\ {\tilde{X}}_{\theta \theta }=\frac{{\overline{\Omega }}^2}{{\Omega }^2}{X}_{\theta \theta }.& (29)\end{array}$

5.2. Relations holding in the stationary regime

Subsection 5.3 will tackle the derivation of stationary initial data in relation to the gauge conditions from the same starting point as [63]. However, a slicing condition that has been tested experimentally is considered, and the calculations will take place in the conformally compactified domain instead of the physical one. The procedure will require knowing the conditions that stationarity puts on the metric, which is taken to be the Schwarzschild one from now onward.

Imposing stationarity on the evolution equation of the metric components allows to find the desired time-independent expressions for the trace of the extrinsic curvature, given in (30) and (31). Setting the RHS of (2.82a) in Vañó-Viñuales [49] to zero together with the evolved Z4 constraint $\tilde{\Theta }=0$ gives the following expression for the quantity $\Delta \tilde{K}=\tilde{K}-K_{CMC}$ (the variation of the physical trace of the extrinsic curvature with respect to the background value K_CMC) in terms of the conformally compactified metric components:

$\begin{array}{l}\Delta \tilde{K}=-K_{CMC}+\frac{\Omega }{\alpha }\left({\beta }^r{}^{\prime }-\frac{3{\beta }^r{\chi }^{\prime }}{2\chi }+\frac{{\beta }^r{\gamma }_{\theta \theta }{}^{\prime }}{{\gamma }_{\theta \theta }}+\frac{{\beta }^r{\gamma }_{rr}{}^{\prime }}{2{\gamma }_{rr}}+\frac{2{\beta }^r}{r}\right)\\ -\frac{3{\beta }^r{\Omega }^{\prime }}{\alpha }.& (30)\end{array}$

In essence, the relation above is equivalent to (7) in Ohme et al. [63]; only here, different variables are used, and the relation holds in the conformally compactified domain. Setting now the RHS of (2.82b) in Vañó-Viñuales [49] to zero provides the following expression for A_rr to hold in the stationary regime:

$\begin{array}{ll}{A}_{rr}=\frac{1}{3\alpha }\left({\beta }^r{\gamma }_{rr}^{\prime }-\frac{{\beta }^r{\gamma }_{rr}{\gamma }_{\theta \theta }^{\prime }}{{\gamma }_{\theta \theta }}+2{\gamma }_{rr}{{\beta }^r}^{\prime }-\frac{2{\beta }^r{\gamma }_{rr}}{r}\right).& (31)\end{array}$

This expression will not be used in further derivations but is provided here for completeness. The stationary expression for Λ^r in terms of the spatial metric components obtained from the Z4 constraint (2.81c) in Vañó-Viñuales [49] is

$\begin{array}{l}{\Lambda }^r=\frac{2{\hat{\gamma }}_{\theta \theta }}{{\hat{\gamma }}_{rr}{\gamma }_{\theta \theta }r}-\frac{2}{{\gamma }_{rr}r}+\frac{{\gamma }_{rr}^{\text{'}}}{2{\gamma }_{rr}^2}-\frac{{\gamma }_{\theta \theta }^{\text{'}}}{{\gamma }_{rr}{\gamma }_{\theta \theta }}+\frac{{\widehat{\gamma }}_{\theta \theta }^{\text{'}}}{{\widehat{\gamma }}_{rr}{\gamma }_{\theta \theta }}-\frac{{\widehat{\gamma }}_{rr}^{\text{'}}}{2{\widehat{\gamma }}_{rr}{\gamma }_{rr}}\\ =\frac{2}{r}\left(\sqrt{{\gamma }_{rr}}-\frac{1}{{\gamma }_{rr}}\right)+\frac{{\gamma }_{rr}^{\text{'}}}{{\gamma }_{rr}^2}.& (32)\end{array}$

After the second equality above, ${\widehat{\gamma }}_{rr}$ = ${\widehat{\gamma }}_{\theta \theta }=1$ have been set and the substitution ${\gamma }_{\theta \theta }={\gamma }_{rr}^{-1/2}$ has been imposed. The latter is required by the introduction of the spatial conformal factor χ in the formulation, and it is to be applied to (30) and (31) as well. The second expression for Λ^r in the stationary regime (32), with its formal divergence as r → 0, is not straightforward to solve. To find finite values of Λ^r for a γ_rr that neither diverges nor goes to zero at the origin fine-tuning is required (with the exception of γ_rr = 1 that gives Λ^r = 0). This condition (32) possibly puts stringent limitations on the possible solutions of the shift equation (12). A way to simplify the problem is to, instead of solving for the shift condition, impose Λ^r = 0 initially and accordingly choose conformally flat initial data, as was mentioned at the beginning of the section and will be used in 5.3.

The Schwarzschild spacetime is a static solution of the Einstein equations, given by (17) with $A(\tilde{r})=1-\frac{2M}{\tilde{r}}$. The following relations between metric components in the uncompactified physical domain hold:

$\begin{array}{ll}{\tilde{\alpha }}^2-{\beta }^{{\tilde{r}}^2}{\tilde{X}}_{\tilde{r}\tilde{r}}={\tilde{\alpha }}^2-{\tilde{\beta }}^2=1-\frac{2M}{\tilde{r}},& (\mathrm{33a})\end{array}$

$\begin{array}{ll}{\tilde{X}}_{rr}=\frac{1}{{\tilde{\alpha }}^2},{\tilde{X}}_{\theta \theta }=1.& (\mathrm{33b})\end{array}$

The first condition is, e.g., (9), in Ohme et al. [63] and corresponds to using the Killing lapse and shift. The second one is used in (8) also in Ohme et al. [63], introducing

$\begin{array}{l}\tilde{\beta }=\sqrt{{\tilde{X}}_{\tilde{r}\tilde{r}}}{\tilde{\beta }}^{\tilde{r}}\equiv \sqrt{\frac{{\gamma }_{\tilde{r}\tilde{r}}}{\tilde{\chi }}}{\beta }^{\tilde{r}}\text{ that transforms via }\tilde{\beta }=\frac{\beta }{\Omega }\text{ to}\\ \beta =\sqrt{X_{rr}}{\beta }^r\equiv \sqrt{\frac{{\gamma }_{rr}}{\chi }}{\beta }^r.& (34)\end{array}$

The last expression in (33b) means that the physical areal radius remains constant.

The equivalent expressions in the conformally compactified domain, obtained from (33) using the transformation relations in (28) and (29), are

$\begin{array}{ll}{\alpha }^2-{\beta }^r{}^2\frac{{\gamma }_{rr}}{\chi }\equiv {\alpha }^2-{\beta }^r{}^2{X}_{rr}={\alpha }^2-{\beta }^2=c^2{\Omega }^2\left(1-\frac{2M\overline{\Omega }}{r}\right),& (\mathrm{35a})\end{array}$

$\begin{array}{ll}\frac{{\gamma }_{rr}}{\chi }\equiv X_{rr}=\frac{c^2}{{\alpha }^2}{\left(\frac{{\Omega }^2}{{\overline{\Omega }}^2}(\overline{\Omega }-r{\overline{\Omega }}^{\prime })\right)}^2,\frac{{\gamma }_{\theta \theta }}{\chi }\equiv X_{\theta \theta }=\frac{{\Omega }^2}{{\overline{\Omega }}^2}.& (\mathrm{35b})\end{array}$

The compactification factor $\overline{\Omega }$ cannot be chosen freely (the conformal factor Ω can), but it is to be substituted from the second relation in (35b), a consequence of the physical areal radius chosen to be constant in (33b). Note the presence of the c factor (completely unrelated to the speed of light) in the first two relations: it corresponds to a constant rescaling of the hyperboloidal time coordinate, t → ct, in the same way as in the transformation to hyperboloidal time (55) in Panosso Macedo [90]. The constant c will take a different value in the stationary regime depending on the gauge equations chosen and the values of their parameters. Relations (35) have been checked experimentally in spherically symmetric hyperboloidal trumpet evolutions. The quantity c can be absorbed into the value of K_CMC (in the conformal factor) in the above expressions, but this approach will not be followed here. It is indeed possible that the c rescaling of the time coordinate in the evolution is a consequence of a rescaling of the conformal factor Ω (chosen to be time-independent) that takes place as a result of the change in the hyperboloidal slices. Understanding this effect is left for future work.

5.3. Solving the slicing condition and imposing an explicitly conformally flat metric

A delicate evolution variable in the BSSN/Z4 formulations is Λⁱ, which is usually set initially to zero corresponding to a choice of conformally flat initial metric, i.e., γ_ij = η_ij. As pointed out after introducing (32), for this spherically symmetric setup where ${\gamma }_{\theta \theta }={\gamma }_{rr}^{-1/2}$ holds and ${\widehat{\gamma }}_{rr}={\widehat{\gamma }}_{\theta \theta }=1$, the problem is considerably simplified if γ_rr = 1 is set as initial condition, ensuring a well-behaved initial Λ^r = 0. This is indeed a conformally compactified version of isotropic coordinates for puncture data [18, 19]. Thus, from now on, the initial metric will be chosen to be conformally flat explicitly. This mimics the procedure done for CMC slices in 4, where the compactification factor is also determined imposing conformal flatness via (23). The difference is that now the slicing will not be CMC but determined by stationarity of the slicing condition (10).

The coupled system of equations to solve is the left equation in (35b) and (10)'s RHS set to zero with $n_{cK}=-m_{cK}\frac{6rℐ}{K_{CMC}}\Omega =m_{cK}(r_{ℐ}^2-r^2)$, where m_cK is a non-zero constant. The reason why n_cK is set to be proportional to the conformal factor is to ensure that it vanishes at ℐ⁺. In this way, the gauge propagation speed associated with the slicing condition will be the physical one at future null infinity (as is the case for the harmonic slicing for n_cK = 0), as is shown in Figure 1. The other relations in 5.2 are used to substitute all quantities in those two equations in terms of β [defined in (34)], the following rescaling of the compactification factor:

$\begin{array}{ll}\overline{\omega }=\frac{\overline{\Omega }}{r\Omega },& (36)\end{array}$

which is expected to be finite and non-zero everywhere in the integration domain, and the constant c, which will be used as parameter to shoot-and-match on during the solving procedure. The profiles of β and $\overline{\omega }$ for the CMC case, with c = 1, are depicted in Figure 8. The explicit substitutions to be performed are

$\begin{array}{l}\overline{\Omega }=r\Omega \overline{\omega },{\beta }^r=\beta \sqrt{\frac{\chi }{{\gamma }_{rr}}},{\gamma }_{\theta \theta }={\gamma }_{rr}^{-1/2},{\gamma }_{rr}=1,\chi =r^2{\overline{\omega }}^2,\\ \alpha =\sqrt{{\beta }^2+c^2{\Omega }^2(1-2M\overline{\omega }\Omega )}.& (37)\end{array}$

FIGURE 8

Figure 8. (A) Profiles for β, with δβ in inset. (B) Profiles for $\overline{\omega }$. CMC trumpet profiles together with solutions to the slicing condition (10) and (38) (with minus sign) for two different sets of parameters. The corresponding values of c are included in the main text.

The expression from χ comes from imposing conformal flatness on the right equation in (35b). The α has been isolated from (35a). The reason for choosing to substitute α instead of β is that the latter, as can be seen in Figure 8A, changes sign over the compactified domain (the shift β^r is positive at the horizon and negative at ℐ⁺) and thus cannot be easily substituted. After the substitutions, the left equation in (35b) reads

$\begin{array}{ll}{\overline{\omega }}^{\prime }=\frac{\overline{\omega }}{\Omega }\left(\pm \frac{\sqrt{{\beta }^2+c^2{\Omega }^2(1-2M\overline{\omega }\Omega )}}{cr}-{\Omega }^{\prime }\right).& (38)\end{array}$

The sign providing the RHS that coincides with CMC's ${\overline{\Omega }}^{\prime }$ from Figure 5 for substituted CMC data is chosen (the minus one, in this case). The resulting equation from (10) is much longer and has been included in Appendix A as (1). Both equations are formally diverging at the trumpet and at ℐ⁺. The ellipticity of the coupled system of equations has not been studied. For convenience, the system is solved for

$\begin{array}{ll}\delta \beta =\beta -\frac{c{K}_{CMC}r}{3}& (39)\end{array}$

instead of β as the former vanishes at ℐ⁺ (see inset in Figure 8A). The explicit form of the conformal factor (9) is also substituted.

Using Taylor expansions in the radial coordinate around the trumpet is not suitable, as α may not necessarily be proportional to an integer power of the radial coordinate [18, 19]. The integration will thus start from ℐ⁺ toward the trumpet, motivating the change of coordinate x = 1 − r. Guesses for the initial values of δβ and $\overline{\omega }$ at ℐ⁺ (x = 0) required to start the shoot-and-match integration from there are obtained by Taylor expanding the equations around ℐ⁺ up to first order. The values of the variables and their first derivatives at future null infinity,

$\begin{array}{l}\delta \beta {|}_{{ℐ}^{+}}\approx 0\\ +x\frac{3(c-1)\left(c^2{\xi }_{cK}{K}_{CMC}^2-3c{K}_{CMC}^2{m}_{cK}-3m_{cK}\left(9{\xi }_{cK}+K_{CMC}^2\right)\right)}{c{K}_{CMC}\left(2\left(3c^2-2\right)K_{CMC}^2+9(3c-4){\xi }_{cK}\right)}& (\mathrm{40a})\end{array}$

$\begin{array}{l}\overline{\omega }{|}_{{ℐ}^{+}}\approx \frac{\left(3c^2-1\right)K_{CMC}^2+9(c-1){\xi }_{cK}}{2c^2{K}_{CMC}^2}·\\ [1+x\frac{1}{c^2{K}_{CMC}^2\left(2\left(3c^2-2\right)K_{CMC}^2+9(3c-4){\xi }_{cK}\right)}·(6c^4{K}_{CMC}^4+\\ +36c^3{\xi }_{cK}{K}_{CMC}^2-c^2{K}_{CMC}^2\left(45{\xi }_{cK}+4K_{CMC}^2+27m_{cK}\right)\\ -243c{\xi }_{cK}{m}_{cK}+27m_{cK}\left(9{\xi }_{cK}+K_{CMC}^2\right))]& (\mathrm{40b})\end{array}$

are obtained by imposing regularity of the expansions at ℐ⁺ for each power of x. The system of equations is solved using Mathematica's NDSolve function on the integration domain x ∈ [10⁻⁵, 1) and a WorkingPrecision of 50. The starting point x = 10⁻⁵ is chosen to avoid the formal divergence of the equations at ℐ⁺, while the integration is carried out up to almost x = 1 (at least as close as x = 0.9996). The starting point needs to look visually very near to future null infinity and be closer to ℐ⁺ than any of the gridpoints in the evolutions (see 6). Conditions (40) are evaluated at x = 10⁻⁵ for the chosen values of m_cK and ξ_cK, and a value for c is set with high precision to start the integration. The criteria to determine whether the found solution is good enough is for the lapse α to be zero at the trumpet. In practice, when using NDSolve, that translates to obtaining profiles of δβ and $\overline{\omega }$ that look smooth and do not diverge at x = 1 – given the difficulty of the integration, a solution with a very small divergence localized beyond x = 0.999 is taken as valid. Unless c is very close to the required value, the solutions very quickly become infinite as being integrated toward the origin, given the formally divergent character of the equations. On top of that, there are some regions in the potential values of c where the solutions become complex, or they cannot be integrated any closer to the origin than a certain point. That point may correspond to the hyperboloidal equivalent of the “critical point” mentioned in Bruegmann [18] and Baumgarte and de Oliveira [19]. Moreover, the level of fine-tuning required for c, so that the solution is well-behaved up to x = 1 (≡ r = 0), is very high. What is meant with “fine-tuning” for c is the system of equations is solved for a specific value of c, depending on the direction in which the obtained profiles for δβ and $\overline{\omega }$ are diverging, the next value of c is selected to decrease the divergence. This procedure is repeated until a value for c that provides regular profiles for the solutions at the origin (and α close enough to zero there) is found. This makes usual methods to choose the values for c in the shooting-and-matching difficult to employ successfully, so that after a careful study of each setup (for chosen values of m_cK and ξ_cK), a manual tuning of c has been used.

The profiles of β and $\overline{\omega }$ for two different parameter choices, with the CMC equivalent for comparison, are displayed in Figure 8. The two configurations considered were (with M = 1 and K_CMC = −1): m_cK = 1 and ξ_cK = 1, requiring c = 0.9996723791389223, and m_cK = 0.1 and ξ_cK = 0.5, with c = 0.9661803490. These values of the parameters were chosen because they provide a long-term stationary solution in evolutions of the Einstein equations (see 6.3 for further comments). The number of significant digits required for the value of c depends on the method and precision used. As can be seen in the plots, the larger m_cK, the more “CMC-like” the slices look near the origin, whereas it has been found experimentally that ξ_cK has more of an impact in the region near ℐ⁺, allowing solutions to be further from the CMC profile for smaller values of ξ_cK. This is also the behavior seen in evolutions of CMC trumpet initial data with those parameter choices for the slicing condition.

With NDSolve, it is not straightforward to estimate the error of the solution (even using options like AccuracyGoal or PrecisionGoal), which in turn makes the study of its convergence difficult. Using x = ·10⁻⁵ as the starting point for the integration gave larger residuals, as expected, but this is not enough to systematically study convergence. In order to overcome this hurdle, simple explicit integrators were implemented to solve the same system of coupled equations. Those were a 1st order Euler method and a 4th order Runge-Kutta (RK4). The explicit integrator was used for the shooting, while a bisection method was used for the matching part (looking how close to zero α|_r = 0 was for the chosen value of c). An example of the results obtained for m_cK = 1 and ξ_cK = 1 (for a value of c = 1.03903643703151 in the Euler method) is shown in Figures 9A, B with a solid black line, also including the NDSolve solution for comparison (black dashed) – there are obvious differences. Two solutions obtained with the RK4 integrator are also shown in Figure 9 in blue, solid for 200 points (with c = 0.9891442037734391040608175) and dashed for 220 points (c = 0.989160040599287855336) although they are only distinguishable from each other in the noisy part near ℐ⁺. While these curves are closer to the NDSolve solution, there are still differences between them. The main obstacle in the explicit integration methods was that ~3% of the gridpoints closest to ℐ⁺ look very noisy, and there is a jump between the values of δβ and $\overline{\omega }$ at both sides of the noise even if a solution for the whole domain was found. Fine-tuning on the correct value of c was difficult (this is why only two solutions are given for the RK4) as the RHS would change sign several times throughout the domain of c considered (quite possibly due to the presence of the noisy part) and not all of the potentially promising values would provide a solution extending all the way to the origin. An extra difficulty was that sometimes the RHS would become complex halfway through the integration and no full solution was found, which also happened with NDSolve. Convergence for the part of the RK4 solution between the origin and the noise is shown in Figure 9C: the Hamiltonian and r-component of the momentum constraints are evaluated for the 200 and 220-point RK4 solutions and rescaled according to the expected 4th order convergence. Both lines overlap perfectly in most of the interior domain.

FIGURE 9

Figure 9. (A) Solutions for δβ. (B) Solutions for $\overline{\omega }$. (C) Constraints for RK4. (A, B): Comparison between the NDSolve results (black dashed) and those obtained for a simple Euler solver (solid black) and the Runge-Kutta 4 (blue lines) for the m_cK = 1 and ξ_cK = 1 case, specifying the number of gridpoints used for Euler and RK4. The non-smooth part near ℐ⁺ for the explicit integration methods (Euler and RK4) is shown in the insets. Note that the profiles of the solutions show a jump to the left and right of the noisy part. (C) Constraints evaluated for the two solutions obtained with the RK4 integrator (with 200 and 220 points, thus a resolution increase of 1.1). The Hamiltonian and r-component of the momentum constraint for the lower resolution are divided by f^p, where f is the increase in resolution of 1.1, and p is the order of convergence (4 for RK4). The curves coincide well in the interior of the domain (the solutions converge there), but the noisy part near ℐ⁺ shown on (A, C) does not converge.

Apart from the unavoidable fact that the equations are formally singular at the extrema, potential explanations for the delicate tuning of c required and the noisy part in the solutions obtained from the explicit integrators are (i) the Taylor expansion (40) used to start the integration are not suitable and (ii) the slicing condition considered (10) together with the imposition of conformal flatness do not provide a solution (see 6.3 for comments on its stationary solution after evolution).

Initial data (constructed from the NDSolve solution) for the evolution variables for the m_cK = 0.1, ξ_cK = 0.5 case, chosen because its solution differs more from the CMC profile, are presented in Figure 10. The quantities are calculated from β, $\overline{\omega }$ and c using relations (37), (30), (31), and (32). The CMC profiles from Figure 7 have also been included in the plot in a light blue color to facilitate comparison. The main qualitative difference between both sets of data is that $\Delta \tilde{K}$ has an positive dependence on r in the non-CMC case. Note that the full trace of the physical extrinsic curvature, $\tilde{K}=K_{CMC}+\Delta \tilde{K}=-1+\Delta \tilde{K}$, is still negative everywhere for the solved-for case. The spatial conformal factor χ is no longer unity at ℐ⁺ and A_rr reaches the origin (the location of the trumpet) with a steeper slope.

FIGURE 10

Figure 10. Solved-for Schwarzschild trumpet BH initial data for the evolution variables with K_CMC = −1, M = 1, m_cK = 0.1, and ξ_cK = 0.5 in black, and CMC data shown in Figure 7 in light blue for comparison. The vertical lines denote the respective location of the horizons. The A_rr component of the trace-free part of the conformal extrinsic curvature shows a steeper slope at the origin for the solved-for data. The values at ℐ⁺ of the corresponding χ, α, and β^r are different from their CMC values, while the solved-for $\Delta \tilde{K}$ is non-zero in the whole domain. As imposed when constructing both trumpet slices, γ_rr = γ_θθ = 1 and, consequently, Λ^r = 0.

The compactification and slicing for the m_cK = 0.1, ξ_cK = 0.5 solution, including the CMC profiles for comparison, are shown in Figure 11. The two compactification factors in Figure 11A are qualitatively the same. The slope at the origin is different because the slices of the solved-for solution reach further into the horizon – the trumpet is located at R₀ = 1.60M. This can be appreciated in Figure 11B. Full details on the construction of Penrose diagrams for numerical data will be included in Vañó-Viñuales [24].

FIGURE 11

Figure 11. (A) Compactification factor. (B) CMC and solved-for slices. Comparison between CMC trumpet and the solution for m_cK = 0.1 and ξ_cK = 0.5. Whereas in the CMC case the trumpet is located at R₀ = 1.91M of the Schwarzschild radius (as indicated in Figure 3C), for the solved-for slices it is at R₀ = 1.60M, closer to the singularity. Accordingly, the slope (~1/R₀) near r = 0 of the compactification factor is steeper. The slices shown on the Penrose diagram correspond to different values of the hyperboloidal time t for CMC and solved-for data; they have been chosen in such a way that the slices coincide at ℐ⁺.

6. Evolution results

6.1. Implementation

Simulations are performed with a spherically symmetric code that uses the method of lines with a 4th order Runge-Kutta time integrator and 4th order finite differences, adding Kreiss-Oliger dissipation [92]. The grid used is staggered (cell-centered), so it avoids the two points where the equations are formally singular - the origin r = 0, that corresponds to the value of the Schwarzschild radial coordinate R₀ where the trumpet asymptotes to, and r = r_ℐ, which corresponds to future null infinity ℐ⁺. Extrapolating boundary conditions such as the outflow boundary conditions in Calabrese and Gundlach [93] are used at both boundaries. Off-centered finite difference stencils in the advection terms' derivatives is known to improve stability and performance of numerical relativity simulations [94–96]. Thus, advection stencils are up-winded toward larger radii (like on the right in figure 2 in Vañó-Viñuales and Husa [48]), where the radial shift component is negative and down-winded toward smaller radii where the shift is positive.

The chosen values of the parameters for the simulations considered here are (unless stated otherwise): κ₁ = 0.5, κ₂ = 0, κ₀ = 0, dissipation strength 0.1, K_CMC = −1, m_cK = 0.1, ξ_cK = 0.5, λ = 1, ${\xi }_{{\beta }^r}=0$, η = 0, 1. Most simulations use 456 spatial discretization points and a timestep of dt = 6.667·10⁻⁴. The number of spatial gridpoints is enough to resolve the hyperboloidal trumpet initial data considered here, while the timestep is chosen to be below the Courant–Friedrichs–Lewy limit. For some configurations (not in the case of the work presented here), the dt needs to be smaller to account for the presence of stiff terms in the equations. Evolutions have been performed with the generalized BSSN system.

As mentioned in the previous section, convergence of the stationary initial data solutions could so far only be shown on part of the integration domain (see Figure 9C). For those data, noisy features makes the reconstructed initial data for the evolution variables (such as A_rr or $\Delta \tilde{K}$) non-smooth enough to pose problems in the evolutions, namely, that the simulations crash due to the spiky profiles. In this study, the focus will be to understand the phenomenological behavior of the solutions. In any case, any reasonable lack of convergence or smoothness of the solved-for initial data will disappear as the evolution progresses as the gauge conditions will drive the data to the real solution. As initial data, the two options depicted in Figure 10 will be used, namely, the hyperboloidal Schwarzschild CMC trumpet data (24) and (25), as well as the solved-for m_cK = 0.1, ξ_cK = 0.5 solution obtained with NDSolve, which from now on will be called “statio” solution. Two different gauge setups will be considered, both using (10) as slicing condition: the first one will use the integrated Gamma-driver (12), of which neither CMC trumpet data nor the statio solution are a stationary solution, while the second one will involve the modified shift condition (13) without advection terms, whose RHS is zero for both sets of initial data.

6.2. Evolution with shift condition including advection terms

The first test is performed with the Gamma-driver (12) shift condition with advection terms and η term, together with the slicing (10). The expectation is that the trumpet readjustment to happen for the statio data should be smaller than for the CMC data.

When evolved with η = 0, a small oscillatory behavior appears in all evolution variables around the initial profiles of the statio data, no matter if the starting point of the simulation is CMC or statio data. The amplitude of these oscillations grows slowly (up to the final t = 100 of these simulations), which indicates that the simulation will crash at some point later in time. The damping term with η was already introduced in Baiotti and Rezzolla [79] to avoid strong oscillations in the shift, while other works have found a small value of η useful to suppress gauge oscillations, such as those affecting eccentricity measurements in Purrer et al. [97]. A further study of suitable values of η and comparison with non-hyperboloidal simulations is left for future work.

If setting η = 1, the dynamics drives the initial data to a stable solution. The differences between the initial and the final profiles are shown in Figure 12. As expected, those differences are larger for CMC initial data, especially the change in A_rr, which gets up to 2.4 (beyond the range shown in the Figure). The two lines for Λ^r lie on top of each other as both their initial and end states are the same. While the statio initial data has indeed the advantage that it has required less trumpet dynamics in the evolution, the profiles of the A_rr and Λ^r quantities near ℐ⁺ look slightly diverging, which may cause convergence problems in general (see comments in next subsection). The depicted changes in α show that its final state is smaller than the CMC initial profile but slightly larger than the statio one. No further study of the parameter space has been performed because anyway it is not yet clear what gauge conditions are best suited for the hyperboloidal setup.

FIGURE 12

Figure 12. In blue are the differences between the initial CMC trumpet data and the final state of the simulation (evaluated at t = 100), while the black lines show the difference between the m_cK = 0.1, ξ_cK = 0.5 initial profile and the final one. The final state is the same for both sets of initial data, as they were both evolved with the same setup: slicing (10) with m_cK = 0.1 and ξ_cK = 0.5, and shift with advection terms (12) with λ = 1 and η = 1. The vertical line denotes the final location of the horizon. The differences are larger for CMC initial data. See main text for further details.

6.3. Evolution with shift condition without advection terms

Both sets of initial data are now evolved with the shift condition (13). The initial dynamics in the CMC case are driven by the slicing condition (10), while the statio solution remains visually static throughout the evolution (ran up to t = 1, 000). Figure 13A aims to capture how fast the final state is attained in the evolutions, via showing the behavior of the L² norm over the whole Λ^r gridfunction over time, for both sets of initial data (CMC and statio), as well as for two runs with CMC initial data and different parameter configurations. The statio Λ^r is set to be exactly zero initially (implied by the explicit conformal flatness imposed) although this cannot be seen in Figure 13A due to the logarithmic plot used in the vertical axis as the evolution starts its value changes probably because the solved-for solution has some small errors and still needs to relax to its truly stationary state. In any case, the change is much smaller than for CMC initial data. Both sets of initial data relax to the final state at the same rate. The quality of the final solution is estimated looking at the values of the Hamiltonian constraint at different values of the radial coordinate as time passes in Figure 13B. In the middle of the compactified domain (dash-dot line), the value of H rapidly approaches a small value. Closer to ℐ⁺ (dashed line), the final value is larger, meaning that the constraint violation there is more pronounced than in the center (which is expected because the equations are formally singular). At the last gridpoint (half a spatial step from the compactified location of future null infinity, solid line), the effect is even larger. One indication that the statio solution must have some small errors is that in the blue solid line, the leftmost value of H for CMC data (at t = 0) is very small – except for the compactification factor, which is solved numerically, the rest is an explicit analytic solution (24). However, the initial value of the Hamiltonian constraint for statio data on the black solid line is large, indicating that the given initial data does not satisfy the constraints in a satisfactory way, at least very near ℐ⁺.

FIGURE 13

Figure 13. (A) L²-type norm for Λ^r over time. (B) Hamiltonian constraint at several r. Long-term behavior of unperturbed simulations with advection-term-free shift condition. (A) L²-type norm for Λ^r's gridfunctions at given times, $\sqrt{{\displaystyle {\sum }_{i=1}^N{({\Lambda }_i^r)}^2}}$ where i ∈ [1, N] denotes the spatial points. As Λ^r = 0 is the stationary solution of (13), this plot shows how fast initial CMC and statio data get to their final state. Apart from the two solid lines representing [data also appearing in (B)], two simulations are included: CMC initial data, and either i) (dashed line) a different value of parameter κ₁ = 1.5 (taken to be 0.5 otherwise) or ii) (dash-dot line) evolved slicing condition with parameter values m_cK = 0.1 and ξ_cK = 0.5. (B) Values of the Hamiltonian constraint over time at different values of the compactified radial coordinate: in the middle of the domain for dash-dot, near ℐ⁺ for dashed and at the closest point to ℐ⁺ for the solid line. Constraint violations are larger closer to future null infinity.

The dashed and dash-dot lines included in Figure 13A correspond to two simulations with CMC initial data but with a different parameter configuration, respectively, κ₁ = 1.5, and slicing with m_cK = ξ_cK = 1. They have been included to shed light on the loss of self-convergence issues detected in the simulations and shown in Figure 14. The L²-type norm of the stationary states of Λ^r for these two different configurations are different from the solid lines (that for the slicing with m_cK = ξ_cK = 1 gets closer to zero).

FIGURE 14

Figure 14. Convergence order over time calculated using the L²-type norm as ${\log }_f\left(\frac{\sqrt{{\displaystyle {\sum }_{i=1}^N{(X_{low,i}-X_{med,i})}^2}}}{\sqrt{{\displaystyle {\sum }_{i=1}^N{(X_{med,i}-X_{high,i})}^2}}}\right)$ with i ∈ [1, N], f = 1.5, N = 304 and X_{low/med/high,i} denoting the evolution variables for the low, medium, and high resolutions at point i. For the blue curves, all the evolution variables were included (so that X was successively $\chi ,{\gamma }_{rr},A_{rr}\Delta \tilde{K},{\Lambda }^r,\alpha$, and β^r), while for the black lines, the summation was performed only on α. The expected 4th order convergence is lost very early, and only recovered later in the evolution, if at all. Interpretation in the main text. Further understanding the convergence behavior in the evolution for the choice of gauge conditions is beyond the scope of this study.

The convergence order calculated from the evolution variables shown in Figure 14 drops from the expected 4 very early in the simulation. Self-convergence of the lapse (black lines) is recovered around t~250. Other fields behave worse (not shown here): for instance and except for the CMC m_cK = ξ_cK = 1 runs, self-convergence of γ_rr is only 2nd order, while A_rr's is between 2nd and 4th. This is why the convergence order that takes into account all of the evolution variables does not go back to 4 even by t = 900. The exception is the CMC m_cK = ξ_cK = 1 case, where the 4th order convergence is recovered around t~300. The reason is that a larger value of ξ_cK damps the deviations of the lapse from its value at ℐ⁺ more efficiently, together with a larger propagation speed near the trumpet given by a bigger m_cK. With just the simple slicing condition (10) considered here, one can appreciate the enormous effect that gauge conditions have on the evolutions. However, convergence is lost before t~300 even in the best case, which is clearly pointing to the presence of a problem. To the author's best knowledge, the shift condition (13) has not been used before and could be partially to blame of the loss of convergence. Other more sophisticated gauge conditions, such as those covered in Vañó-Viñuales and Husa [48], could be tested for comparison – although the statio solution would not be a stationary solution anymore. During the time when self-convergence of the variables was lost virtually everywhere in the domain (0 ≲ t ≲ 300 or later), the constraints continued to converge in the interior at the domain at all times. However, at late times, they did not asymptote to zero but to ~10⁻⁶. Understanding all these interesting aspects and how they relate to the choice of gauge conditions is left for future research.

To study the robustness of the gauge system and study how constraint convergence evolves in the simulations, a constraint-satisfying Gaussian-like gauge perturbation is included in the initial lapse,

$\begin{array}{ll}\alpha ={\alpha }_0+A{e}^{-\frac{{(r^2-r_c^2)}^2}{4{\sigma }^4}},& (41)\end{array}$

where α₀ denotes either its CMC value as in (24) or its statio value. The chosen values of the parameters are A = 0.05, σ = 0.1, and r_c = 0.25. This initial data are run with 304, 456, and 684 spatial points (and corresponding timestep dt = 10⁻³, 6.667·10⁻⁴, 4.444·10⁻⁴), thereby increasing the resolution by 1.5 between runs.

The initial gauge perturbation extends to all evolution variables and gets propagated away, part into the BH and part out through ℐ⁺, leaving behind what under visual inspection looks like the statio solution. Even if the statio initial data do not appropriately converge, the evolution of the gauge perturbation will after a certain amount of time. Looking at convergence of the constraints can give an estimate of when that happens as well as the general reliability of the simulation. Figure 15 shows convergence of the Hamiltonian constraint at an early time t = 0.1 and a later one t = 4 for both CMC and statio initial data. Except at the boundaries, convergence for the CMC case is good as the blue curves coincide very well in the interior of the domain for both figures. That is not the case for the statio case: initially the lack of convergence of the initial data is seen clearly in the non-coincidence between the black curves in Figure 15A. Later, at t = 4 as shown in Figure 15B, most of the discrepancies have disappeared and convergence looks better – except again at the extrema of the radial coordinate, where the Hamiltonian constraint looks very noisy. This is not necessarily an indicator of a problem, as the constraints are formally divergent at the trumpet and at ℐ⁺.⁷

FIGURE 15

Figure 15. (A) Rescaled H constraint at t = 0.1. (B) Rescaled H constraint at t = 4. Rescaled values of the Hamiltonian constraint as a function of r at two different instants of time. H_low uses 304 points, H_med uses 456, and H_high uses 684. The rescalings are of the form f^p, where f is the increase in resolution between runs (1.5 in this case) and p is the order of convergence (4th in this setup). The same legend is valid for both plots.

7. Conclusion

The hyperboloidal approach allows numerical simulations to reach future null infinity from first principles and without complicated constructions. While it needs to be further understood and developed, progress in the non-linear regime is taking place along several fronts [81, 98, 99] and will also benefit from work in the linear one, e.g., [100]. The focus in the present study has been on spherically symmetric hyperboloidal trumpet initial data for puncture-type evolutions of BHs. More specifically, the construction of CMC trumpet data via the height function approach has been reviewed and adapted to the needs of numerical simulations using the puncture approach together with conformal compactification. Gauge conditions play a crucial role in numerical evolutions both in the stability of the setup and the allowed final states of the system. While understanding them within the hyperboloidal approach is still work in progress, some options providing successful numerical simulations are known.

Availability of stationary initial data that are suitable for the evolution and study of perturbations thereof is very desired, especially as part of the development of the hyperboloidal method. This study sets the infrastructure to pursue those solutions for spherically symmetric BH trumpet data within the conformally compactified domain. A procedure to solve a specific numerically-tested slicing condition together with explicit conformal flatness has been developed and tested in examples. The accuracy of the numerical solutions was not fully satisfactory (convergence could only be checked in part of the domain), but still they could be tested in hyperboloidal evolutions and shed some light on the behavior of some choices of shift conditions. The bottomline is that new initial data that approaches the stationary solution of the gauge conditions much more rapidly than CMC data was constructed despite how challenging the procedure ultimately was.

There are several options that could be tested to improve the quality of the stationary solutions and that have been left for future work. At the numerical level, a more sophisticated solving method could be used, such as an elliptic solver for non-linear equations or a relaxation method. The main requirement is to obtain reliable and smooth solutions for which convergence can be checked in the whole integration domain. At the analytical level, a different slicing condition could be considered. To be suitable, it needs to provide a stable stationary final state in evolution for a hyperboloidal trumpet slice. This is not easy to attain, but progress on the gauge condition front in the near future will provide more options. This progress will also contribute to understand the convergence problems detected during the evolutions, which will be solved elsewhere.

The insight gained on the effects of the two different shift conditions tested here will be used to develop other suitable options that may provide smoother profiles of the evolution variables near ℐ⁺. Within the free evolution setup with a BSSN/Z4-type formulation, the shift condition is closely related to the quantity Λ^a, which may impose some limitations as to which final states are allowed or not. Either modifying the definition of Λ^a, making it more compatible with the hyperboloidal framework, or considering a different formulation of the Einstein equations may be beneficial. In either case, it would be worth attempting to solve the shift condition for the compactification, once the behavior of Λ^a is better understood, instead of imposing conformal flatness as done here.

While CMC trumpet initial data for the RN spacetime has been developed, it is still waiting to be tested in hyperboloidal evolution to the author's best knowledge. Moreover, how far the conformally compactified height function approach can be extended to include a non-vanishing cosmological constant, in a similar way to [101] but with views toward non-linear numerical simulations, is still to be found out. Finally, including a massless scalar field perturbation on stationary trumpet initial data will allow to study the former's behavior without mixing it with trumpet relaxation dynamics, which was one of the problems hit in Vañó-Viñuales [49] and Vañó-Viñuales and Husa [61]. Of special relevance are the decay tails and their convergence at ℐ⁺, as preparation of the GWs to be treated in the full 3D case.

This study succeeded in taking a few steps toward a more thorough understanding of the interplay between gauge conditions and stationary solutions on hyperboloidal trumpet slices for the puncture approach. It has provided a framework to determine those solutions for suitable generic gauge conditions as well as exemplifying how challenging both the initial data calculations and the evolutions are. Insights into how to develop better suited formulations for the hyperboloidal approach have also been gained.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

AVV has performed the derivations and simulations described in this article, agrees to be accountable for all aspects of it, and approves its publication.

Funding

The author thanks the Fundacão para a Ciência e Tecnologia (FCT), Portugal, for the financial support to the Center for Astrophysics and Gravitation (CENTRA/IST/ULisboa) through the Grant Project No. UIDB/00099/2020. This work was also supported through the European Research Council Consolidator Grant 647839.

Acknowledgments

The author thanks Edgar Gasperin, Sascha Husa, and David Hilditch for valuable comments on the manuscript. The author would also like to thank Sascha Husa for important feedback in parts of the research presented here. Discussions with him, Sergio Dain, Niall O'Murchadha, and David Hilditch considerably helped understand and solve challenges in this study. Most of the algebraic derivations were performed using the Mathematica package xAct [102].

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

1. ^Past experiments in spherical symmetry (subsection 8.2.1 in Vañó-Viñuales [49] also mentioned here at the end of 3.1) have shown an instability-inducing drift in the variables, not linked to any specific part of the domain. This was related mainly to gauge conditions and how they deal with the trumpet and ℐ⁺ asymptotics.

2. ^This is why the constant parameter K_CMC introduced in (9) and described in 4.3 is negative for hyperboloidal slices reaching future null infinity. If a positive value is chosen for it, then the hyperboloidal slices intersect past null infinity.

3. ^The time-independent quantity n_cK, a function of the radial coordinate, has here a different expression from that used in Vañó-Viñuales and Husa [48].

4. ^There is another potential drawback to this approach: a change in the mass of the BH (for instance, due to some energy that is accreted by it during evolution) would in principle not be taken into account by the source functions, and the gauge conditions may try to force the system into an inappropriate geometry. There is the possibility, at least in spherical symmetry, to evaluate numerically the new value of the BH's mass “on the fly” during the evolution, use it to calculate the new trumpet geometry, and update the source terms accordingly. An example of this recalculation of the BH's mass and the trumpet is shown for some evolution variables in figures 8.30 and 8.31 in Vañó-Viñuales [49].

5. ^The initial ansatz can be as general as $A(\tilde{r})=1-\frac{2M}{\tilde{r}}+\frac{Q^2}{{\tilde{r}}^2}+\frac{\Lambda }{3}{\tilde{r}}^2$, but in this work only the case with vanishing cosmological constant is considered.

6. ^The profiles of the evolution variables at several instances of time during a collapse simulation are shown in figure 8.14 in Vañó-Viñuales [49].

7. ^Constraint equations, zero for the continuum solution, do not have a scale. Self-convergence of the evolution fields is what needs to be satisfactory there. As described above, this is not the case, but this will be studied elsewhere.

References

1. Abbott R, Abbott TD, Abraham S, Acernese F, Ackley K, Adams C, et al. Open data from the first and second observing runs of advanced LIGO and Advanced Virgo. SoftwareX. (2021) 13:100658. doi: 10.1016/j.softx.2021.100658

CrossRef Full Text | Google Scholar

2. Abbott R, Abe H, Acernese F, Ackley K, Adhicary C, Adhikari N, et al. Open data from the third observing run of LIGO, Virgo, KAGRA, and GEO. arXiv [Preprint]. (2023). arXiv: 2302.03676. doi: 10.3847/1538-4365/acdc9f

CrossRef Full Text | Google Scholar

3. Barack L. Late time dynamics of scalar perturbations outside black holes. 1. A Shell toy model. PhysRev. (1999) D59:044016. doi: 10.1103/PhysRevD.59.044016

CrossRef Full Text | Google Scholar

4. Leaver EW. Solutions to a generalized spheroidal wave equation. J Math Phys. (1986) 27:1238. doi: 10.1063/1.527130

CrossRef Full Text | Google Scholar

5. Leaver EW. Spectral decomposition of the perturbation response of the Schwarzschild geometry. Phys Rev D. (1986) 34:384–408. doi: 10.1103/PhysRevD.34.384

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Seidel E, Suen WM. Towards a singularity proof scheme in numerical relativity. Phys Rev Lett. (1992) 69:1845–8. doi: 10.1103/PhysRevLett.69.1845

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Boyle M, Hemberger D, Iozzo DAB, Lovelace G, Ossokine S, Pfeiffer HP, et al. The SXS collaboration catalog of binary black hole simulations. Class Quant Grav. (2019) 36:195006. doi: 10.1088/1361-6382/ab34e2

CrossRef Full Text | Google Scholar

8. Brill DR, Lindquist RW. Interaction energy in geometrostatics. Phys Rev. (1963) 131:471–6. doi: 10.1103/PhysRev.131.471

CrossRef Full Text | Google Scholar

9. Brandt S, Brügmann B. A simple construction of initial data for multiple black holes. Phys Rev Lett. (1997) 78:3606–9. doi: 10.1103/PhysRevLett.78.3606

CrossRef Full Text | Google Scholar

10. Beig R, Husa S. Initial data for general relativity with toroidal conformal symmetry. Phys Rev D. (1994) 50:R7116–8. doi: 10.1103/PhysRevD.50.R7116

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Campanelli M, Lousto CO, Marronetti P, Zlochower Y. Accurate evolutions of orbiting black-hole binaries without excision. Phys Rev Lett. (2006) 96:111101. doi: 10.1103/PhysRevLett.96.111101

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Baker JG, Centrella J, Choi DI, Koppitz M, van Meter J. Gravitational wave extraction from an inspiraling configuration of merging black holes. Phys Rev Lett. (2006) 96:111102. doi: 10.1103/PhysRevLett.96.111102

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Hannam M, Husa S, Brügmann B, Gonzalez JA, Sperhake U. Where do moving punctures go? J Phys Conf Ser. (2007) 66:012047. doi: 10.1088/1742-6596/66/1/012047

CrossRef Full Text | Google Scholar

14. Hannam M, Husa S, Pollney D, Brügmann B, O'Murchadha N. Geometry and regularity of moving punctures. Phys Rev Lett. (2007) 99:241102. doi: 10.1103/PhysRevLett.99.241102

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hannam M, Husa S, Ohme F, Brügmann B, O'Murchadha N. Wormholes and trumpets: the Schwarzschild spacetime for the moving-puncture generation. Phys Rev. (2008) D78:064020. doi: 10.1103/PhysRevD.78.064020

CrossRef Full Text | Google Scholar

16. Baumgarte TW, Naculich SG. Analytical representation of a black hole puncture solution. Phys Rev. (2007) D75:067502. doi: 10.1103/PhysRevD.75.067502

CrossRef Full Text | Google Scholar

17. Dennison KA, Baumgarte TW. A simple family of analytical trumpet slices of the schwarzschild spacetime. Class Quant Grav. (2014) 31:117001. doi: 10.1088/0264-9381/31/11/117001

CrossRef Full Text | Google Scholar

18. Bruegmann B. Schwarzschild black hole as moving puncture in isotropic coordinates. Gen Rel Grav. (2009) 41:2131–51. doi: 10.1007/s10714-009-0818-6

CrossRef Full Text | Google Scholar

19. Baumgarte TW, de Oliveira HP. Bona-Masso slicing conditions and the lapse close to black-hole punctures. Phys Rev D. (2022) 105:064045. doi: 10.1103/PhysRevD.105.064045

CrossRef Full Text | Google Scholar

20. Li SE, Baumgarte TW, Dennison KA, de Oliveira HP. Dynamical perturbations of black-hole punctures: effects of slicing conditions. Phys Rev D. (2023) 107:064003. doi: 10.1103/PhysRevD.107.064003

CrossRef Full Text | Google Scholar

21. Li SE, Baumgarte TW, Dennison KA, de Oliveira HP. Bona-Massó slices of Reissner-Nordström spacetimes. Phys Rev D. (2022) 106:104059. doi: 10.1103/PhysRevD.106.104059

CrossRef Full Text | Google Scholar

22. Dennison KA, Baumgarte TW, Montero PJ. Trumpet slices in kerr spacetimes. Phys Rev Lett. (2014) 113:261101. doi: 10.1103/PhysRevLett.113.261101

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Heissel G. Towards a Numerical Derivation of Maximal Kerr Trumpet Initial Data. Cardiff: Cardiff University (2017).

Google Scholar

24. Vañó-Viñuales A. Penrose Diagrams for Asymptotically Flat Spherically Symmetric Spacetimes. (2023) In preparation.

Google Scholar

25. Winicour J. Characteristic evolution and matching. Living Rev Relat. (2009) 12:3. doi: 10.12942/lrr-2009-3

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Bishop NT. Numerical relativity: combining the Cauchy and characteristic initial value problems. Class Quantum Grav. (1993) 10:333–41. doi: 10.1088/0264-9381/10/2/015

CrossRef Full Text | Google Scholar

27. Szilagyi B. Cauchy Characteristic Matching in General Relativity. Pittsburgh, PA: University of Pittsburgh (2000).

Google Scholar

28. Bishop NT, Gomez R, Lehner L, Winicour J. Cauchy characteristic extraction in numerical relativity. Phys Rev. (1996) D54:6153–65. doi: 10.1103/PhysRevD.54.6153

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Reisswig C, Bishop NT, Pollney D, Szilagyi B. Unambiguous determination of gravitational waveforms from binary black hole mergers. Phys Rev Lett. (2009) 103:221101. doi: 10.1103/PhysRevLett.103.221101

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Babiuc MC, Winicour J, Zlochower Y. Binary black hole waveform extraction at null infinity. Class Quant Grav. (2011) 28:134006. doi: 10.1088/0264-9381/28/13/134006

CrossRef Full Text | Google Scholar

31. Moxon J, Scheel MA, Teukolsky SA, Deppe N, Fischer N, Hébert F, et al. SpECTRE Cauchy-characteristic evolution system for rapid, precise waveform extraction. Phys Rev D. (2023) 107:064013. doi: 10.1103/PhysRevD.107.064013

CrossRef Full Text | Google Scholar

32. Friedrich H. Cauchy problems for the conformal vacuum field equations in general relativity. Comm Math Phys. (1983) 91:445–72. doi: 10.1007/BF01206015

CrossRef Full Text | Google Scholar

33. Friedrich H. On the existence of n-geodesically complete or future complete solutions of Einstein's field equations with smooth asymptotic structure. Comm Math Phys. (1986) 107:587–609. doi: 10.1007/BF01205488

CrossRef Full Text | Google Scholar

34. Frauendiener J. Conformal infinity. Living Rev Relat. (2004) 7:1. doi: 10.12942/lrr-2004-1

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Friedrich H. Conformal Einstein evolution. In: Frauendiener H, Friedrech H, editors. The Conformal Structure of Space-Time: Geometry, Analysis, Numerics. Berlin; Heidelberg: Springer (2002). p. 1-50. doi: 10.1007/3-540-45818-2_1

CrossRef Full Text | Google Scholar

36. Mitman K, Moxon J, Scheel MA, Teukolsky SA, Boyle M, Deppe N, et al. Computation of displacement and spin gravitational memory in numerical relativity. Phys Rev D. (2020) 102:104007. doi: 10.1103/PhysRevD.102.104007

CrossRef Full Text | Google Scholar

37. Vañó-Viñuales A, Husa S, Hilditch D. Spherical symmetry as a test case for unconstrained hyperboloidal evolution. Class Quant Grav. (2015) 32:175010. doi: 10.1088/0264-9381/32/17/175010

CrossRef Full Text | Google Scholar

38. Penrose R. Asymptotic properties of fields and space-times. Phys Rev Lett. (1963) 10:66–8. doi: 10.1103/PhysRevLett.10.66

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Frauendiener J. Numerical treatment of the hyperboloidal initial value problem for the vacuum Einstein equations. 2. The Evolution equations. Phys Rev. (1998) D58:064003. doi: 10.1103/PhysRevD.58.064003

CrossRef Full Text | Google Scholar

40. Husa S. Problems and successes in the numerical approach to the conformal field equations. LectNotes Phys. (2002) 604:239–60. doi: 10.1007/3-540-45818-2_12

CrossRef Full Text | Google Scholar

41. Hilditch D. Dual Foliation Formulations of General Relativity. (2015). doi: 10.48550/arXiv.1509.02071

CrossRef Full Text | Google Scholar

42. Hilditch D, Harms E, Bugner M, Rüter H, Brügmann B. The evolution of hyperboloidal data with the dual foliation formalism: mathematical analysis and wave equation tests. Class Quant Grav. (2018) 35:055003. doi: 10.1088/1361-6382/aaa4ac

CrossRef Full Text | Google Scholar

43. Scheel MA, Pfeiffer HP, Lindblom L, Kidder LE, Rinne O, Teukolsky SA. Solving Einstein's equations with dual coordinate frames. Phys Rev D. (2006) 74:104006. doi: 10.1103/PhysRevD.74.104006

CrossRef Full Text | Google Scholar

44. Zenginoğlu A. Hyperboloidal foliations and SCRI-fixing. Class Quant Grav. (2008) 25:145002. doi: 10.1088/0264-9381/25/14/145002

CrossRef Full Text | Google Scholar

45. Zenginoğlu A. Hyperboloidal evolution with the Einstein equations. Class Quant Grav. (2008) 25:195025. doi: 10.1088/0264-9381/25/19/195025

CrossRef Full Text | Google Scholar

46. Zenginoğlu A. A hyperboloidal study of tail decay rates for scalar and Yang-Mills fields. Class Quant Grav. (2008) 25:175013. doi: 10.1088/0264-9381/25/17/175013

CrossRef Full Text | Google Scholar

47. Zenginoğlu A. A Conformal Approach to Numerical Calculations of Asymptotically Flat Spacetimes. Potsdam: Max Planck Institute for Gravitational Physics (AEI) and University of Potsdam, Institute of Physics and Astronomy (2007).

Google Scholar

48. Vañó-Viñuales A, Husa S. Spherical symmetry as a test case for unconstrained hyperboloidal evolution II: gauge conditions. Class Quant Grav. (2018) 35:045014. doi: 10.1088/1361-6382/aaa4e2

CrossRef Full Text | Google Scholar

49. Vañó-Viñuales A. Free Evolution of the Hyperboloidal Initial Value Problem in Spherical Symmetry. Palma: Universitat de les Illes Balears (2015). Available online at: http://inspirehep.net/record/1407828/files/arXiv:1512.00776.pdf (accessed July 07, 2023).

Google Scholar

50. Malec E, O'Murchadha N. The General spherically symmetric constant mean curvature foliations of the Schwarzschild solution. Phys Rev. (2009) D80:024017. doi: 10.1103/PhysRevD.80.024017

CrossRef Full Text | Google Scholar

51. Cruz-Osorio A, Gonzalez-Juarez A, Guzman FS, Lora-Clavijo FD. Numerical solution of the wave equation on particular space-times using CMC slices and scri-fixing conformal compactification. Rev Mex Fis. (2010) 56:456–68. doi: 10.1063/1.3473871

CrossRef Full Text | Google Scholar

52. Lee KW, Lee YI. Spacelike Spherically Symmetric CMC Hypersurfaces in Schwarzschild Spacetimes (I): Construction. (2011) 11. doi: 10.48550/arXiv.1806.06638

CrossRef Full Text | Google Scholar

53. Tuite P, O'Murchadha N. Constant Mean Curvature Slices of the Reissner-Nordström Spacetime. (2013). doi: 10.48550/arXiv.1307.4657

CrossRef Full Text | Google Scholar

54. Lee KW. Constant Mean Curvature Foliation Properties in the Extended Reissner-Nordstrom Spacetimes. (2018). doi: 10.48550/arXiv.1111.2679

CrossRef Full Text | Google Scholar

55. Schneemann C. Numerische Berechnung von hyperboloidalen Anfangsdaten für die Einstein-Gleichungen. (2006).

Google Scholar

56. Schinkel D, Panosso Macedo R, Ansorg M. Axisymmetric constant mean curvature slices in the Kerr space-time. Class Quant Grav. (2014) 31:075017. doi: 10.1088/0264-9381/31/7/075017

CrossRef Full Text | Google Scholar

57. Schinkel D, Ansorg M, Panosso Macedo R. Initial data for perturbed Kerr black holes on hyperboloidal slices. ClassQuantGrav. (2014) 31:165001. doi: 10.1088/0264-9381/31/16/165001

CrossRef Full Text | Google Scholar

58. Buchman LT, Pfeiffer HP, Bardeen JM. Black hole initial data on hyperboloidal slices. Phys Rev. (2009) D80:084024. doi: 10.1103/PhysRevD.80.084024

CrossRef Full Text | Google Scholar

59. Bardeen JM, Buchman LT. Bondi-Sachs energy-momentum for the CMC initial value problem. Phys Rev. (2012) D85:064035. doi: 10.1103/PhysRevD.85.064035

CrossRef Full Text | Google Scholar

60. Schinkel D. Anfangsdaten für Schwarze Löcher auf Hyperboloidalen Blättern. (2016).

Google Scholar

61. Vañó-Viñuales A, Husa S. Unconstrained hyperboloidal evolution of black holes in spherical symmetry with GBSSN and Z4c. J Phys Conf Ser. (2015) 600:012061. doi: 10.1088/1742-6596/600/1/012061

CrossRef Full Text | Google Scholar

62. Vañó-Viñuales A, Husa S. Free hyperboloidal evolution in spherical symmetry. In: Proceedings, 14th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics, and Relativistic Field Theories (MG14) (In 4 Volumes): Rome, Italy, July 12-18, 2015, Vol. 2. Rome (2017), p. 2025–30. Available online at: https://inspirehep.net/record/1415732/files/arXiv:1601.04079.pdf (accessed July 07, 2023).

Google Scholar

63. Ohme F, Hannam M, Husa S, O'Murchadha N. Stationary hyperboloidal slicings with evolved gauge conditions. Class Quant Grav. (2009) 26:175014. doi: 10.1088/0264-9381/26/17/175014

CrossRef Full Text | Google Scholar

64. Misner CW, Thorne KS, Wheeler JA. Gravitation. San Francisco, CA: W.H. Freeman and Co. (1973).

Google Scholar

65. Nakamura T, Oohara K, Kojima Y. General relativistic collapse to black holes and gravitational waves from black holes. Prog Theor Phys Suppl. (1987) 90:1–218. doi: 10.1143/PTPS.90.1

CrossRef Full Text | Google Scholar

66. Shibata M, Nakamura T. Evolution of three-dimensional gravitational waves: harmonic slicing case. Phys Rev D. (1995) 52:5428–44. doi: 10.1103/PhysRevD.52.5428

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Baumgarte TW, Shapiro SL. On the numerical integration of Einstein's field equations. Phys Rev. (1999) D59:024007. doi: 10.1103/PhysRevD.59.024007

CrossRef Full Text | Google Scholar

68. Brown JD. BSSN in spherical symmetry. Class Quant Grav. (2008) 25:205004. doi: 10.1088/0264-9381/25/20/205004

CrossRef Full Text | Google Scholar

69. Bona C, Ledvinka T, Palenzuela C, Zacek M. General-covariant evolution formalism for numerical relativity. Phys Rev. (2003) D67:104005. doi: 10.1103/PhysRevD.67.104005

CrossRef Full Text | Google Scholar

70. Alic D, Bona-Casas C, Bona C, Rezzolla L, Palenzuela C. Conformal and covariant formulation of the Z4 system with constraint-violation damping. Phys Rev. (2012) D85:064040. doi: 10.1103/PhysRevD.85.064040

CrossRef Full Text | Google Scholar

71. Sanchis-Gual N, Montero PJ, Font JA, Müller E, Baumgarte TW. Fully covariant and conformal formulation of the Z4 system in a reference-metric approach: comparison with the BSSN formulation in spherical symmetry. Phys Rev. (2014) D89:104033. doi: 10.1103/PhysRevD.89.104033

CrossRef Full Text | Google Scholar

72. Bernuzzi S, Hilditch D. Constraint violation in free evolution schemes: comparing BSSNOK with a conformal decomposition of Z4. Phys Rev. (2010) D81:084003. doi: 10.1103/PhysRevD.81.084003

CrossRef Full Text | Google Scholar

73. Weyhausen A, Bernuzzi S, Hilditch D. Constraint damping for the Z4c formulation of general relativity. Phys Rev. (2012) D85:024038. doi: 10.1103/PhysRevD.85.024038

CrossRef Full Text | Google Scholar

74. Husa D. Numerical relativity with the conformal field equations. In: Fernandez-Jambrina L, Gonzalez-Romero LM, editors. Current Trends in Relativistic Astrophysics: Theoretical, Numerical, Observational, Vol. 617. Berlin; Heidelberg: Springer (2002). p. 159–92. doi: 10.1007/3-540-36973-2_9

CrossRef Full Text | Google Scholar

75. Moncrief V, Rinne O. Regularity of the Einstein equations at future null infinity. Class Quant Grav. (2009) 26:125010. doi: 10.1088/0264-9381/26/12/125010

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Rinne O. An axisymmetric evolution code for the Einstein equations on hyperboloidal slices. ClassQuantGrav. (2010) 27:035014. doi: 10.1088/0264-9381/27/3/035014

CrossRef Full Text | Google Scholar

77. Morales MD, Sarbach O. Evolution of scalar fields surrounding black holes on compactified constant mean curvature hypersurfaces. Phys Rev. (2017) D95:044001. doi: 10.1103/PhysRevD.95.044001

CrossRef Full Text | Google Scholar

78. Bona C, Massó J, Seidel E, Stela J. A new formalism for numerical relativity. Phys Rev Lett. (1995) 75:600–3. doi: 10.1103/PhysRevLett.75.600

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Baiotti L, Rezzolla L. Gauge conditions for long term numerical black hole evolutions without excision. Phys Rev. (2003) D67:084023. doi: 10.1103/PhysRevD.67.084023

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Friedreich H, Rendall AD. The cauchy problem for the Einstein equations. In: Schmidt BG, editor. Einstein's Field Equations and Their Physical Implications: Selected Essays in Honour of Jurgen Ehlers, Vol. 540. Berlin; Heidelberg: Springer (2000). p. 127–224. doi: 10.1007/3-540-46580-4_2

CrossRef Full Text | Google Scholar

81. Vañó-Viñuales A. Free Hyperboloidal Evolution in 3D: Gauge Waves and Constraint-Violating Gravitational Waves. (2023) In preparation.

Google Scholar

82. Duarte M, Feng JC, Gasperín E, Hilditch D. Regularizing dual-frame generalized harmonic gauge at null infinity. Class Quant Grav. (2023) 40:025011. doi: 10.1088/1361-6382/aca383

CrossRef Full Text | Google Scholar

83. Alcubierre M. The Appearance of coordinate shocks in hyperbolic formalisms of general relativity. Phys Rev. (1997) D55:5981–91. doi: 10.1103/PhysRevD.55.5981

CrossRef Full Text | Google Scholar

84. Baumgarte TW, Hilditch D. Shock-avoiding slicing conditions: tests and calibrations. Phys Rev D. (2022) 106:044014. doi: 10.1103/PhysRevD.106.044014

CrossRef Full Text | Google Scholar

85. Gentle AP, Holz DE, Kheyfets A, Laguna P, Miller WA, Shoemaker DM. Constant crunch coordinates for black hole simulations. Phys Rev. (2001) D63:064024. doi: 10.1103/PhysRevD.63.064024

CrossRef Full Text | Google Scholar

86. Malec E, O'Murchadha N. Constant mean curvature slices in the extended Schwarzschild solution and collapse of the lapse. Part I. Phys Rev. (2003) D68:124019. doi: 10.1103/PhysRevD.68.124019

CrossRef Full Text | Google Scholar

87. Estabrook F, Wahlquist H, Christensen S, DeWitt B, Smarr L, Tsiang E. Maximally slicing a black hole. Phys Rev D. (1973) 7:2814–7. doi: 10.1103/PhysRevD.7.2814

CrossRef Full Text | Google Scholar

88. Iriondo M, Malec E, O'Murchadha N. The Constant mean curvature slices of asymptotically flat spherical space-times. Phys Rev. (1996) D54:4792–8. doi: 10.1103/PhysRevD.54.4792

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Malec E, O'Murchadha N. Constant Mean Curvature Slices in the Extended Schwarzschild Solution and Collapse of the Lapse. Part II. (2003). doi: 10.48550/arXiv.gr-qc/0307047

CrossRef Full Text | Google Scholar

90. Panosso Macedo R. Hyperboloidal framework for the Kerr spacetime. Class Quant Grav. (2020) 37:065019. doi: 10.1088/1361-6382/ab6e3e

CrossRef Full Text | Google Scholar

91. Ansorg M, Panosso Macedo R. Spectral decomposition of black-hole perturbations on hyperboloidal slices. Phys Rev. (2016) D93:124016. doi: 10.1103/PhysRevD.93.124016

CrossRef Full Text | Google Scholar

92. Kreiss HO, Oliger J. Methods for the Approximate Solution of Time Dependent Problems. GARP publications series No 10 International Council of Scientific Unions. Geneva: World Meteorological Organization (1973).

Google Scholar

93. Calabrese G, Gundlach C. Discrete boundary treatment for the shifted wave equation. ClassQuantGrav. (2006) 23:S343–68. doi: 10.1088/0264-9381/23/16/S04

CrossRef Full Text | Google Scholar

94. Zlochower Y, Baker JG, Campanelli M, Lousto CO. Accurate black hole evolutions by fourth-order numerical relativity. Phys Rev D. (2005) 72:024021. doi: 10.1103/PhysRevD.72.024021

CrossRef Full Text | Google Scholar

95. Husa S, Gonzalez JA, Hannam M, Brügmann B, Sperhake U. Reducing phase error in long numerical binary black hole evolutions with sixth order finite differencing. Class Quant Grav. (2008) 25:105006. doi: 10.1088/0264-9381/25/10/105006

CrossRef Full Text | Google Scholar

96. Chirvasa M, Husa S. Finite difference methods for second order in space, first order in time hyperbolic systems and the linear shifted wave equation as a model problem in numerical relativity. J Comput Phys. (2010) 229:2675–96. doi: 10.1016/j.jcp.2009.12.016

CrossRef Full Text | Google Scholar

97. Purrer M, Husa S, Hannam M. An Efficient iterative method to reduce eccentricity in numerical-relativity simulations of compact binary inspiral. Phys Rev D. (2012) 85:124051. doi: 10.1103/PhysRevD.85.124051

CrossRef Full Text | Google Scholar

98. Peterson C, Gautam S, Rainho I, Vañó-Viñuales A, Hilditch D. 3D evolution of a semilinear wave model for the Einstein field equations on compactified hyperboloidal slices. arXiv [Preprint]. (2023). arXiv: 2303.16190. doi: 10.1103/PhysRevD.108.024067

CrossRef Full Text | Google Scholar

99. Feng J, Gasperin E. Linearised conformal Einstein field equations. Class Quant Grav. (2023) 40:175001. doi: 10.1088/1361-6382/ace606

CrossRef Full Text | Google Scholar

100. Jaramillo JL, Panosso Macedo R, Al Sheikh L. Pseudospectrum and black hole quasinormal mode instability. Phys Rev X. (2021) 11:031003. doi: 10.1103/PhysRevX.11.031003

PubMed Abstract | CrossRef Full Text | Google Scholar

101. Bizoń P, Chmaj T, Mach P. A toy model of hyperboloidal approach to quasinormal modes. Acta Phys Polon B. (2020) 51:1007. doi: 10.5506/APhysPolB.51.1007

CrossRef Full Text | Google Scholar

102. Martín-García JM. xAct: Efficient Tensor Computer Algebra for Mathematica. Available online at: http://www.xact.es/ (accessed July 07, 2023).

Google Scholar

Appendix

A1. Initial data equation from slicing condition

Equation resulting from substituting relations (37) into (10) and solving for β′:

$\begin{array}{l}{\beta }^{\prime }=\frac{1}{9\alpha cr{\Omega }^2{\overline{\omega }}^2(c^2{K}_{CMC}\Omega (2M\Omega \overline{\omega }-1)+6m_{cK}{r}_{ℐ})}[\alpha c{K}_{CMC}\overline{\omega }(-9{\Omega }^2(c^2(\alpha K_{CMC}-{\xi }_{cK})+\frac{3\alpha {\xi }_{cK}}{\sqrt{K_{CMC}^2{r}^2+9{\Omega }^2}})\\ -18c^{2}M{\Omega }^3\overline{\omega }({\xi }_{cK}-\alpha K_{CMC})-\frac{\alpha K_{CMC}^2{r}^2(3{\xi }_{cK}+K_{CMC}r{\Omega }^{\prime })}{\sqrt{K_{CMC}^2{r}^2+9{\Omega }^2}}+\Omega (\frac{\alpha K_{CMC}^3{r}^2}{\sqrt{K_{CMC}^2{r}^2+9{\Omega }^2}}+54\alpha m_{cK}{r}_{ℐ}))\\ -9{\beta }^3\Omega {\overline{\omega }}^2(c^2{K}_{CMC}\Omega (9M\Omega \overline{\omega }-4)+12m_{cK}{r}_{ℐ})+9\beta c\Omega {\overline{\omega }}^2(10c^3{K}_{CMC}{M}^2{\Omega }^5{\overline{\omega }}^2+2{\Omega }^3(c^3{K}_{CMC}+12cM{m}_{cK}{r}_{ℐ}\overline{\omega })\\ -9c^3{K}_{CMC}M{\Omega }^4\overline{\omega }-\alpha c^2{K}_{CMC}r\Omega {\Omega }^{\prime }+2c{\Omega }^2(\alpha c{K}_{CMC}Mr\overline{\omega }{\Omega }^{\prime }-6m_{cK}{r}_{ℐ})\\ +6\alpha m_{cK}r{r}_{ℐ}{\Omega }^{\prime })+9\alpha {\beta }^{2}c{K}_{CMC}\overline{\omega }({\xi }_{cK}-\alpha K_{CMC})+18{\beta }^5{K}_{CMC}{\overline{\omega }}^2].& (42)\end{array}$

Here, α is kept for readability and is to be substituted using (37). The coordinate location of ℐ⁺ at r_ℐ = 1 is also to be set.