Spherically symmetric black hole spacetimes on hyperboloidal slices

Gravitational radiation and some global properties of spacetimes can only be unambiguously measured at future null inﬁnity ( I + ). This motivates the interest in reaching it within simulations of coalescing compact objects, whose waveforms are extracted for gravitational wave modelling purposes. One promising method to include future null inﬁnity in the numerical domain is the evolution on hyperboloidal slices: smooth spacelike slices that reach future null inﬁnity. The main challenge in this approach is the treatment of the compactiﬁed asymptotic region at I + . Evolution on a hyperboloidal slice of a spacetime including a BH entails an extra layer of diﬃculty, in part due to the ﬁnite coordinate distance between the BH and future null inﬁnity. Spherical symmetry is considered here as 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¨om BH 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 speciﬁc sets of parameters. Finally, results of testing these initial data in evolution are presented.


I. 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 I + , 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][4][5], such as GW interferometers, so that is where GWs signals should ideally be extracted from simulations.
Two main description 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.Still, 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][9][10].In evolutions 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][14][15], where the proper distance becomes infinite while reaching towards the symmetric point to future timelike infinity i + .For embedding diagrams of the wormhole and trumpet geometries see e.g.figures 1 and 2 in [15].Construction of maximal trumpet slices has been tackled for Schwarzschild [16][17][18][19][20], for Reissner-Nordström (RN) [21], for Kerr [22,23].Asymptotically, the slices considered in those works are spacelike Cauchy, and thus reach spatial infinity i 0 .The trumpet puncture approach is also chosen in the present work, 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 behaviour of slices inside of the horizon, useful e.g. 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][26][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][29][30][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][33][34][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 travelling out is slowed down and becomes underresolved, the outward propagation speed of signals on compactified hyperboloidal slices is of order unity and they reach I + at a finite coordinate time without any loss of resolution.Figure 5 in [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 I + with a finite value of the coordinates.The core idea is that instead of working with the physical metric gab 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 I + at the appropriate rate ḡab = Ω 2 gab . (1) One of the most difficult aspects of this approach to the hyperboloidal initial value problem [35,39] is the regularisation of the resulting formally singular equations (see (2) in section II) in a way that works numerically and avoids instabilities arising from the continuum equations, in particular for hyperbolic free evolution schemes as considered here.At the analytical level, the equations were shown to be manifestly regular at I + [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 [43], aims to minimise the divergent terms in the equations, making them as regular as possible.The present implementation follows Zenginoglu's approach [44][45][46][47] to conformal compactification, using free evolution and a timeindependent conformal factor Ω. Stable evolutions in spherical symmetry of regular initial data that do not form BHs were presented in [37], while [48] covers experiments with suitable hyperbolic gauge conditions.Evolving a hyperboloidal slice of spacetime including BHs is particularly challenging1 , especially in the puncture approach where both the regions inside of the horizons and the asymptotic far field are compactified.Constant-meancurvature (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][51][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 a generalization 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 works have considered hyperboloidal initial data including BHs.Configurations in spherical and in axial symmetry were presented in [55], while [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 [58], whereas properties such as the Bondi-Sachs energy and momentum of the above setups were presented in [59].The binary BH scenario was also studied in [60].However, these works were designed with the aim to treat the BHS via excision, and thus not a lot of effort was put into regularising 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 previous work [49,61,62], 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 [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 [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 work are to review hyperbolic CMC trumpet BH initial data suitable for numerical evolutions with the puncture approach in mind (section IV), and to set the basic infrastructure in terms of initial data and gauge conditions to calculate stationary trumpet slices (section V).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 (section VI).For this purpose, the already nontrivial hyperboloidal initial value problem in spherical symmetry is considered, as it still contains the critical part of the radial treatment.
This paper is organised as follows: in section II the used formulation of the conformally compactified Einstein equations is briefly described, and the gauge conditions considered here are covered in section III.Initial data including a BH is treated in the following two sections: as constant-mean-curvature (CMC) in section IV, while an example of solving hyperboloidal-adapted gauge conditions is provided in section V. Section VI presents basic evolution results, and final thoughts on this work are gathered in the conclusions.The appendix collects an equation used in section V.Sections II, III and IV cover previously treated material, while sections V and VI 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 negative2 value means expansion of the normals.Notation for the metrics is the same one as used in [37]: the 4-dimensional physical spacetime metric is denoted as g, the 4D conformal metric as ḡ, the 3D conformal spatial metric (induced by ḡ) as γ, the 3D twice conformal metric γ and the 3D twice conformal background metric γ.

II. FORMULATION
The emphasis in this work 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 The well-posed formulations considered are either the Generalized BSSN [65][66][67][68] or a similar conformal version of the Z4 [69][70][71], the Z4c equations [72,73].The full derivation of these equations in terms of the conformally rescaled metric is described in [37] and in Chapter 2 of [49].The equations used in the simulations are those included in appendix C of [37] (or appendix A in [48]) and again in Chapter 2 of [49].There is a modification related to the evolution of BH spacetimes: a constraint damping term of the form where Z r is a Z4 variable and κ 0 a freely specifiable parameter, can be added to Λ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 where γab is the spatial metric induced from ḡab , and χ is the spatial conformal factor.The conformal extrinsic curvature tensor Kab is decomposed into its conformal trace-free part and (in this formulation) its physical trace, mixed with the physical Z4 variable Θ, Evolved are A ab , and K's variation with respect to its initial value ∆ K = K − K0 = K − K CM C (this last parameter will be explained in subsection IV C).The quantity Θ is evolved as well if using the Z4 formulation.The Z4 variable Z a is absorbed into the vector where Γ a bc are the Christoffel symbols calculated from γ ab and Γa bc the ones built from a time-independent background metric γ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 β i .

A. Spherically symmetric reduction variables
The following spherically symmetric ansatz is used for the spherically symmetric line element in the conformally compactified domain (with σ 2 ≡ dθ 2 + sin 2 θdϕ 2 ) where the freedom introduced by the spatial conformal factor χ is fixed by eliminating γ θθ = γ −1/2 rr .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 .Also 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 , ∆ K, Λ r , α, β r and Θ.
The conformal factor Ω is set to be a time-independent function of the compactified radial coordinate r as with r I the coordinate location of future null infinity (set to r I = 1 in the implementation without restricting generality) and K CM C a negative parameter described in subsection IV C.This expression satisfies that Ω(r) is a regular function that becomes zero at I , with non-vanishing derivative there (compare e.g.[55,74]).The origin of this expression is explained in subsection IV E.

III. GAUGE CONDITIONS
Hyperboloidal constrained evolutions [75][76][77] have used suitable gauges imposed via the resolution of elliptic constraint equations.In this work 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 β i , control the behaviour 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 [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 the outside, making it a natural boundary for the numerical integration domain, where no boundary conditions need to be imposedradiation just needs to be allowed to leave the spacetime.It is possible to fix I + to a specific coordinate location in the numerical grid (r I (9) in the present setup) for compactified hyperboloidal slices.This procedure is called scri-fixing [39,44].
The background behaviour of hyperboloidal slices differs from Cauchy ones in that the trace of the physical extrinsic curvature K is non-zero asymptotically.This requires a modification of the usual slicing conditions commonly used in numerical simulations.See e.g. 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 [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, here by including a BH in the spacetime, and elsewhere by being tested in the full 3D case [81].Work towards finding suitable gauge conditions [82] is also being tackled from the dual foliation approach.

A. Hyperbolic gauge conditions tested with BH spacetimes
When applying the gauge conditions discussed in [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 neighbourhood of a BH if excision is not used.However, near I + 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 I + and something different close to the BH.A condition that has provided successful evolutions using trumpet initial data is, with ˙≡ ∂ t and where ξ cK is a parameter used to damp the behaviour of the lapse at I + .This equation is equivalent to (20) in [48] with ξ 1 α = ξ cK , ξ 2 = 0 and α 2 f = n cK + α 2 , later setting n cK3 to be proportional to Ω.Note that the coefficient in front of ∆ K is similar to the shock-avoiding slicing condition [20,83,84].This form was chosen for the following considerations.The α 2 part provides physical propagation speeds for the gauge modes (the first three lines listed in figure 1), as mentioned above.This is desired at I + , 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.
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: α = 0 ↔ α = α, β r = βr , ∆ 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 Γa bc in (7), are They are obtained (subsection IV C) from ( 21) or (24) setting A( r Ω ) = 1, M = 0, C CM C = 0 and Ω = Ω.For the shift condition, two different options are considered.One is a variant of the integrated Gamma-driver [79] adapted to hyperboloidal slices mostly the same as (26) in [48].The coefficient in front of Λ r is chosen in such a way that the associated propagation speeds will be the physical ones at I + .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 section V).However, dropping the advection terms modifies the characteristic propagation associated to the shift condition.In order to ensure that the related ingoing speed at I + is still zero, the coefficient in front of Λ r is modified as The resulting outgoing propagation speed is also modified: it is smaller (although still positive) at I + (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 behaviour in numerical evolutions, both at the origin and at I + .Hyperboloidal CMC trumpet initial data (( 24) and (25) as derived in section IV, or any initial data satisfying the relations in subsection V B) are a stationary solution of the Einstein equations as described in section II.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 e.g. the Zero speed c 0 = −β r and incoming (−) and outgoing (+) lightspeeds c ± = −β r ± α χ γrr plotted for a CMC slice of a Schwarzschild BH with M = 1, using K CM C = −1 and r I = 1.The vertical line indicates the radial position of the horizon, where c 0 , c − < 0 and c + = 0.At I + we have c 0 , c + > 0 and c − = 0. Except for c 0 , 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 subsection V C. While the speeds go to zero at the location of the trumpet (r = 0) and coincide whith the lightspeeds at I + , 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 I + , but that forced the outgoing one to be equal to c 0 (instead of c + ) there.Compare to a Minkowski equivalent (with different parameter choices) in figure 1 in [48].plot on the right of figure 2 in [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 [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 is a stationary but not a stable solution for the gauge conditions with trumpet source terms (more on this in section V).Still, the growth in these simulations is slow enough to study small scalar field perturbations, as presented in [61].Whether a different choice of trumpet slice or form of the gauge conditions would not cause the growth is an open question 4 .Meanwhile, an attempt to combine stability and stationarity together is described in 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 [49] section V, where a solution for the gauge conditions with hyperboloidal Minkowski source functions is determined for a specific setup.

IV. CONSTANT-MEAN-CURVATURE INITIAL DATA A. 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 t via the height function h(r) [85,86] as The height function satisfies dh/dr < 1 everywhere except asymptotically, where dh/dr| I = 1 holds, thus characterizing the hyperboloidal slices as spacelike but extending to I + .In order to reach future null infinity with a finite value of the spatial coordinates, the radial coordinate r on a hyperboloidal slice is compactified into a new r using a compactification factor Ω(r) Following (1), the line element is conformally rescaled by the conformal factor Ω, to provide regular metric components all the way to I + The compactification factor Ω 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) behaviour near I + , their behaviour 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.

B. 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, expressed first in terms of the usual time 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 where A and h ′ are functions of r Ω , while Ω and Ω depend on r.

C. 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 ( K) is a constant.A special case of CMC slices are maximal slices [87], where K = 0 and spatial infinity is reached asymptotically.Maximal Schwarzschild trumpet slices are analytically described in [16].Generalizations for slices with a non-vanishing K are considered in [85,86,88,89] and including the critical case for trumpets in [50,58], while a study of CMC slices in the RN geometry has been performed in [53].
This derivation of a height function providing CMC slices follows [85,86].See e.g.subsection 3.2.2 in [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 Setting it equal to a constant value of K = K CM C < 0 and introducing C CM C as an integration constant, the first derivative of the height function is isolated to give The expression for the flat spacetime case is obtained by setting A(r) = 1 and C CM C = 0, and the height function can be integrated to 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 0 ≡ X(t = 0) is used: A height function determined by imposing CMC is not the only suitable choice.It could also be only asymptotically CMC [57].Other possible options are e.g.[42], where h ′ (r) is chosen to provide unit outgoing radial coordinate lightspeed, and [90], where h(r) is introduced as part of the minimal gauge [91].

D. Hyperboloidal CMC trumpet slices
The choice of the integration constant C CM C 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 CM C depends on M , Q and K CM C and is calculated [16,87] by setting to zero the discriminant of the denominator (6th order polynomial) of5 γ rr 0 (r) = 1 With this critical choice of C CM C (see (3.42) in [49] for the explicit expression for RN), the denominator now has a double real root at r = R 0 .This finite value of the radial coordinate r is where the slice (reaching I + 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 CM C = 0 the double root is R 0 = 3M/2 [16,87], whereas for K CM C → −∞ it tends to R 0 → 2M .The dependence of the double root R 0 on the charge Q and the value of The effect of C CM C 's value on the CMC slices is depicted in figure 1 in [58] and in [53], and illustrated in figure 3 in the form of Carter-Penrose diagrams of the Schwarzschild spacetime.For a value of C CM C smaller than the critical one, the denominator of (22) has two different real roots (R 1 , R 2 ) for r, the outer one R 2 corresponding to the location r2 , but in this work only the case with vanishing cosmological constant is considered.

FIG. 2:
Innermost value of the Schwarzschild-like radial coordinate (the double root R 0 ) reached by the outer CMC trumpet slice, as a function of K CM C and Q (for zero cosmological constant).Taken from [49].
of the minimal surfaces mentioned in [58].If (such as the example shown in figure 3a), the slices reach inside of the white hole, while for a larger value of C CM C (figure 3b) they enter the BH.Quantities become complex for r ∈ (R 1 , R 2 ); the corresponding part of the diagrams is left in white.As mentioned above, for the critical value of C CM C a double root appears and complete CMC trumpet slices (figure 3c) exist, joining either I + to the symmetric point to future timelike infinity (i + ) (the outer slices) or the singularity to i + (the inner ones).For C CM C 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 CM C for different values of K CM C is given in figure 3.5 in [49], showing the maximal case K = 0 corresponding to the usual trumpet slices [16] in the first subfigure.For a positive value of K CM C (in the current sign convention), the hyperboloidal slices reach past null infinity I − instead of I + .Equivalent Penrose diagrams depicting CMC slices for the RN spacetime (with ) are shown in figures 3.11, 3.12, 3.13 and 3.14 in [49].The non-extremal case is illustrated by the choice Q = 0.9M .The trumpet case (with critical C CM C ) should compare to panel 2 in figure 2 in [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 [49] and it has the feature that all slices are trumpet ones with R 0 = M always, as can also be seen in figure 2. For the over-extreme case (Q > M ) no critical value of C CM C is found, thus no trumpet slices can be constructed (at least with this method).
The effect of C CM C can also be seen in the Schwarzschild examples of the height function h(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 (r = 2M ): it diverges downward for C CM C = 2 because it enters the white hole, and upwards for the cases crossing the BH horizon.The region between (R 1 , R 2 ) is complex for the subcritical values of C CM C (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 r = R 0 = 1.905.For the supercritical C CM C = 4, the height function attains a finite value at the singularity r = 0.The integration constants for each case have been set in such a way that as r → ∞ the Schwarzschild height functions approach the flat spacetime one Initial data developed in [56][57][58] aims for its evolution using excision.For that purpose, any CMC slices with (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 that the most suitable choice is the CMC trumpet slice for the critical value of C CM C .This is the case that will be considered from now on.
(a) Incomplete slices entering the white hole.[49].In subfigures 3a and 3b, the height function is imaginary between R 1 and R 2 .For the critical C CM C (subfigure 3c), the slices that reach the singularity and those that reach I + are separated by a thick line located at r = R 0 that corresponds to the double root.In the C CM C = 4 case on subfigure 3d, the hyperboloidal slices go all the way from I + into the singularity.This last diagram is qualitatively the same (with different values of the parameters) as the Penrose diagram in figure 10 in [44].The present numerical experiments use the outer slices (those reaching I + ) of subfigure 3c, as the critical value of C CM C allows to map the trumpet slice to the whole of the radial coordinate range r ∈ (0, r I ].

E. Compactification Ω and conformal rescaling Ω for CMC slices
It is convenient to determine the compactification factor Ω by imposing a conformally flat initial spatial metric, which is in a way equivalent to transforming to the isotropic radial coordinate, as this is also a simple choice compatible with an initial zero (32).The factor Ω is expected to vanish at the same rate as the conformal factor Ω at I + , but it is allowed to have a different behaviour elsewhere.Expression ( 23) is solved numerically for Ω; a suitable procedure is described in subsection 6.6.1 of [49].The result for the Schwarzschild case is shown as the solid line in figure 5, with the dotted line representing the linear behaviour of Ω near r = 0, corresponding to the asymptotic behaviour near the trumpet r = R 0 (see (5) in [16]).The resulting Ω will only extend to r = 0 if the critical value of C CM C is used.If C CM C is larger than the critical value, Ω will diverge as r → 0, while using a value smaller than the critical C CM C will provide a Ω that does not reach the origin (and does not vanish at the smaller r it gets to).
C CMC = 2 In flat spacetime (A( r Ω ) = 1), condition ( 23) can be solved analytically resulting in expression (9), where Ω is substituted by Ω.This result is commonly seen in the literature [55,74], and has been used in preceding work [37,48] as conformal and compactification factors for the flat spacetime case.In order to ensure that the slices reach I + (instead of ending at spacelike infinity i 0 ), (9) satisfies that Ω| I = 0 and ∇ a Ω| I ̸ = 0.Then, its behaviour at I + is unaffected by the choices of the parameters M and C CM C 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 work dealing with BH spacetimes, as already mentioned towards the end of subsection II A. The profile of this choice of Ω appears in figure 5  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 Ω 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 recognised 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 [49] for a representation of the lapse) compare to figure 5 in [13] and figure 2 in [16], with the difference that here the data are compactified on a hyperboloidal slice instead of a spacelike one (with vanishing K).

F. Vacuum initial data: Schwarzschild spacetime
From now on, the described CMC trumpet BH initial data, suitable for evolutions using the formalism in section II, will be restricted to the Schwarzschild (A(r) = 1 − 2M r ) spacetime.The Reissner-Nordström case works equivalently.The consideration of a non-vanishing cosmological constant is left for future work.After imposing conformal flatness (23), Ω′ can be isolated from there and introduced into (21), yielding for the metric components and for the extrinsic curvature ( K = K CM C ) and rest of the variables (calculated from (30), (31) and ( 32) using( 24)) The background metric γab is chosen to be the initial value of the evolved one, that is, conformally flat γrr = γθθ = 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 subsection III A, CMC initial data (24) and ( 25) are a stationary solution of the Einstein equations, but not of the gauge conditions described in section III.-4. -3. -2. -1. 0. -0.25

V. STABLE STATIONARY INITIAL DATA FOR GIVEN GAUGE CONDITIONS
Here "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 section IV an example was given describing how CMC initial data were a stationary solution of the gauge conditions with CMCconstructed 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 section IV, 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 [62]), where the latter also coincides with the final state of a collapsed scalar field creating a BH with the same total mass6 (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 CM C ) is to impose the conditions (listed in subsection V B) that satisfy the Einstein equations in the non-dynamical regime, and then solve the gauge conditions for stationarity.The latter means setting α = 0 and β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 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 I + 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 I + (see fourth panels in figure 8.31 in [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 I + ).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 subsection V C, and solved for two example configurations.

A. 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σ 2 ≡ dθ 2 + sin 2 θdϕ 2 ) 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 t instead and relating to (17a)).
For the conformally compactified version, as will be used in initial data calculations in subsection V C and relates to the physical one as in (16), set 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 [49]) The quantity χ 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 B. Relations holding in the stationary regime Subsection V C will tackle the derivation of stationary initial data in relation to the gauge conditions, from the same starting point as [63].However, here 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 timeindependent expressions for the trace of the extrinsic curvature, given in (30) and (31).Setting the RHS of (2.82a) in [49] to zero together with the evolved Z4 constraint Θ = 0, gives the following expression for the quantity ∆ K = K − K CM C (the variation of the physical trace of the extrinsic curvature with respect to the background value K CM C ) in terms of the conformally compactified metric components In essence, the relation above is equivalent to (7) in [63], only here different variables are used and the relation holds in the conformally compactified domain.Setting now the RHS of (2.82b) in [49] to zero provides the following expression for A rr to hold in the stationary regime 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 [49] is After the second equality above γrr = γθθ = 1 have been set and the substitution γ θθ = γ −1/2 rr 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 not 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 subsection V C.
The Schwarzschild spacetime is a static solution of the Einstein equations, given by ( 17) with A(r) = 1 − 2M r .The following relations between metric components in the uncompactified physical domain hold: The first condition is e.g. ( 9) in [63] and corresponds to using the Killing lapse and shift.The second one is used in (8) also in [63], introducing 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 The compactification factor Ω cannot be chosen freely (the conformal factor Ω can), but it is to be substituted from the second relation in (35b), 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 → c t, in the same way as in the transformation to hyperboloidal time (55) in [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 CM C (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.

C. Solving the slicing condition and imposing an explicitly conformally flat metric
A delicate evolution variable in the BSSN/Z4 formulations is Λ i , 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 γ θθ = γ −1/2 rr holds and γrr = γθθ = 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 section IV, 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 , 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 I + .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 subsection V B are used to substitute all quantities in those two equations in terms of β (defined in (34)), the following rescaling of the compactification factor 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 ω for the CMC case, with c = 1, are depicted in figure 8.The explicit substitutions to be performed are 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 I + ) and thus cannot be easily substituted.After the substitutions, the left equation in (35b) reads The sign providing the RHS that coincides with CMC's Ω′ 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 (A1).Both equations are formally diverging at the trumpet and at I + .The ellipticity of the coupled system of equations has not been studied.For convenience, the system is solved for instead of β, as the former vanishes at I + (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 I + towards the trumpet, motivating the change of coordinate x = 1 − r.Guesses for the initial values of δβ and ω at I + (x = 0) required to start the shoot-and-match integration from there are obtained by Taylor-expanding the equations around I + up to first order.The values of the variables and their first derivatives at future null infinity, are obtained by imposing regularity of the expansions at I + for each power of x.The system of equations is solved using Mathematica's NDSolve function on the integration domain x ∈ [10 −5 , 1) and a WorkingPrecision of 50.The starting point x = 10 −5 is chosen to avoid the formal divergence of the equations at I + , 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 I + than any of the gridpoints in the evolutions (see section VI).Conditions ( 40) are evaluated at x = 10 −5 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 ω that look smooth and do not diverge at x = 1 -given the difficulty of the integration, a solution with a very small divergence localised 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 towards 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 [18,19].Besides, 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 and, depending on the direction in which the obtained profiles for δβ and ω 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 ω 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 CM C = −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 subsection VI C 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 I + , allowing solutions to be further from the CMC profile for smaller values of ξ cK .This is also the behaviour 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 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 −5 as 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 subfigures 9a and 9b 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 I + .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 the ∼ 3% of the gridpoints closest to I + look very noisy and there is a jump between the values of δβ and ω 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 subfigure 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.
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, ii) the slicing condition considered (10) together with the imposition of conformal flatness do not provide a solution (see subsection VI C 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 β, ω 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 ∆ K has an positive dependence on r in the non-CMC case.Note that the full trace of the physical extrinsic curvature, K = K CM C + ∆ K = −1 + ∆ K, is still negative everywhere for the solved-for case.The spatial conformal factor χ is no longer unity at I + and A rr reaches the origin (the location of the trumpet) with a steeper slope.
The compactification and slicing for the m cK = 0.1, ξ cK = 0.5 solution, including the CMC profiles for comparison, is 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 0 = 1.60M .This can be appreciated in figure 11b.Full details on the construction of Penrose diagrams for numerical data will be included in [24].).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 I + shown on subfigures 9a and 9b does not converge.

A. 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 that 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 0 where the trumpet asymptotes to, and r = r I , which corresponds to future null infinity I + .Extrapolating boundary conditions like the outflow boundary conditions in [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][95][96].Thus, advection stencils are upwinded towards larger radii (like on the right in figure 2 in [48]) where the radial shift component is negative, and down-winded towards smaller radii where the shift is positive.
The chosen values of the parameters for the simulations considered here are (unless stated otherwise): κ 1 = 0.5, κ 2 = 0, κ 0 = 0, dissipation strength 0.1, K CM C = −1, m cK = 0.1, ξ cK = 0.5, λ = 1, ξ β r = 0, η = 0, 1.Most simulations use 456 spatial discretization points and a timestep of dt = 6.667 • 10 −4 .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 subfigure 9c).For those data, noisy features makes the reconstructed initial data for the evolution variables (such as A rr or ∆ K) non-smooth enough to pose problems in the evolutions, namely that the simulations crash due to the spiky profiles.Here the focus will be to understand the phenomenological behaviour 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)  -4. -3. -2. -1. 0.
-0.25 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.

B. 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 behaviour 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 [79] to avoid strong oscillations in the shift, while other works have found a small value of η useful to suppress gauge oscillations, like those affecting eccentricity measurements in [97].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 FIG.11: 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 0 = 1.91M of the Schwarzschild radius (as indicated in figure 3c), for the solved-for slices it is at R 0 = 1.60M , closer to the singularity.Accordingly, the slope (∼ 1/R 0 ) 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 I + .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 I + 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.

C. 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 = 1000).Figure 13a aims to capture how fast the final state is attained in the evolutions, via showing the behaviour of the L 2 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 veritcal 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 I + (dashed line), the final value is larger, meaning the constraint violation there is more pronounced than in the centre (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 I + .
The dashed and dash-dot lines included in subfigure 13a correspond to two simulations with CMC initial data, but with a different parameter configuration: respectively κ 1 = 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 FIG.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.
L 2 -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).
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 eolution variables does not go back to 4 even by t = 900.The exception is the CMC m cK = ξ cK = 1 case, where 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 I + 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.Still, 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 [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 −6 .Understanding all these interesting aspects and how they 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 subfigure 13b), two simulations are included: CMC initial data, and either i) (dashed line) a different value of parameter κ 1 = 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 I + for dashed and at the closest point to I + for the solid line.Constraint violations are larger closer to future null infinity.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, where here α 0 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 is run with 304, 456 and 684 spatial points (and corresponding timestep dt = 10 −3 , 6.667 • 10 −4 , 4.444 • 10 −4 ), so 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 I + , leaving behind what under visual inspection looks like the statio solution.Even if the statio initial data do not appropriately converge, 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 I +7 .
All: m cK =0.1, ξ cK =0.5 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 χ, γ rr , A rr ∆ K, Λ r , α 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 behaviour in the evolution for the choice of gauge conditions is beyond the scope of this work.

VII. CONCLUSIONS
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, like e.g.[100].The focus in the present work 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 work 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 behaviour 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

FIG. 15:
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.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 I + .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 behaviour 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.Also, how far the conformally compactified height function approach can be extended to include a non-vanishing cosmological constant, in a similar way to e.g.[101] but with views towards 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 behaviour without mixing it with trumpet relaxation dynamics, which was one of the problems hit in [49,61].Of special relevance are the decay tails and their convergence at I + , as preparation of the GWs to be treated in the full 3D case.This work succeeded in taking a few steps towards 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.
(d) Slices reaching the singularity.

FIG. 3 :
FIG.3: Carter-Penrose diagrams representing hyperboloidal CMC Schwarzschild foliations with M = 1 and K CM C = −1 and the the choices of C CM C of the height functions depicted in figure4.They correct the versions included in figures 3.3 and 3.4 in[49].In subfigures 3a and 3b, the height function is imaginary between R 1 and R 2 .For the critical C CM C (subfigure 3c), the slices that reach the singularity and those that reach I + are separated by a thick line located at r = R 0 that corresponds to the double root.In the C CM C = 4 case on subfigure 3d, the hyperboloidal slices go all the way from I + into the singularity.This last diagram is qualitatively the same (with different values of the parameters) as the Penrose diagram in figure10in[44].The present numerical experiments use the outer slices (those reaching I + ) of subfigure 3c, as the critical value of C CM C allows to map the trumpet slice to the whole of the radial coordinate range r ∈ (0, r I ].

FIG. 4 :
FIG. 4:Example of numerically integrated height functions for the Schwarzschild spacetime with M = 1 andK CM C = −1.The values of C CM C correspond to those used in figure3.The integration constant has been set in such a way that as r → ∞ the Schwarzschild height function approaches the flat spacetime one (h(r) = (3/K CM C ) 2 + r2 ).In blue the parts of the height function that are imaginary (for C CM C smaller than the critical one).

FIG. 5 :
FIG.5: Conformal (Ω) and compactification ( Ω) factors, as well as the behaviour of Ω ∼ r/R 0 near the origin (with R 0 the location of the trumpet in uncompactified Schwarzschild radial coordinate), for K CM C = −1 and critical C CM C .This figure is also included in[61].

FIG. 7 :
FIG. 7: CMC Schwarzschild trumpet BH initial data for the evolution variables with K CM C = −1, M = 1 and critical C CM C .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 I + .Compare to solved-for stationary data shown in figure 10.

CMCmFIG. 8 :
FIG.8: 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.
Solutions for ω.

H 200 / 1 .1 4 M- 3 . × 10 - 7 - 2 . × 10 - 7 - 1 . × 10 - 7 0 1 . × 10 - 7 2. × 10 - 7 3. × 10 - 7 r
FIG. 9: (A) and (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 I + 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 I + shown on subfigures 9a and 9b does not converge.

5 FIG. 10 :
FIG.10: Solved-for Schwarzschild trumpet BH initial data for the evolution variables with K CM C = −1, M = 1, m cK = 0.1 and ξ cK = 0.5 in black, and CMC data shown in figure7in 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 I + of the corresponding χ, α and β r are different from their CMC values, while the solved-for ∆ K is non-zero in the whole domain.As imposed when constructing both trumpet slices, γ rr = γ θθ = 1 and consequently Λ r = 0.