Spherical Evolution of the Generalized Harmonic Gauge Formulation of General Relativity on Compactified Hyperboloidal Slices
Christian Peterson Shalabh Gautam Alex Vañó-Viñuales David Hilditch





§ ABSTRACT 
  We report on the successful numerical evolution of the compactified
  hyperboloidal initial value problem in general relativity using
  generalized harmonic gauge. We work in spherical symmetry, using a
  massless scalar field to drive dynamics. Our treatment is based on
  the dual-foliation approach, proceeding either by using a height
  function or by solving the eikonal equation to map between frames.
  Both are tested here with a naive implementation and with
  hyperboloidal layers. We present a broad suite of numerical
  evolutions, including pure gauge perturbations, constraint violating
  and satisfying data with and without scalar field matter. We present
  calculations of spacetimes with a regular center. For black hole
  spacetimes we use excision to remove part of the black hole
  interior. We demonstrate both pointwise and norm convergence at the
  expected rate of our discretization. We present evolutions in which
  the scalar field collapses to form a black hole. Evolving nonlinear
  scalar field perturbations of the Schwarzschild spacetime, we
  recover the expected quasinormal frequencies and tail decay rates
  from linear theory.








  ^1CENTRA, Departamento de Física, Instituto Superior Técnico IST,
  Universidade de Lisboa UL, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal

  ^2International Centre for Theoretical Sciences (ICTS), Survey No. 151,
  Shivakote, Hesaraghatta Hobli, Bengaluru - 560 089, India


Spherical Evolution of the Generalized Harmonic Gauge
  Formulation of General Relativity on Compactified Hyperboloidal
  Slices
    David Hilditch^10000-0001-9960-5293
    Received xx; accepted xx
================================================================================================================================






§ INTRODUCTION



Asymptotic flatness is the natural assumption under which to model
isolated systems in general relativity (GR) <cit.>. It may be
formalized in a variety of ways, but common to all is the existence of
a region far from the `center' in which the metric becomes ever closer
to that of the Minkowski spacetime. This leads to the definition of
future null infinity, ℐ^+, which can be thought of as the
collection of endpoints of future directed null geodesics within this
asymptotic region. Future null infinity is crucial in various
mathematical definitions and, crucially for astrophysics, is the place
where gravitational waves (GW) can be unambiguously defined. It is
therefore of vital importance to numerical relativity (NR) to have
access to it.

The most common approach used in NR is to solve the Einstein field
equations (EFEs) in a truncated domain with a timelike outer boundary,
evaluating an approximation to outgoing waves on a set of concentric
spheres and then extrapolate this data to infinity at fixed retarded
time. Eventually waves computed in this way will be affected by
artificial boundary conditions. There are however several proposals to
include ℐ^+ within the computational domain directly. One
popular suggestion is to solve the field equations on compactified
outgoing null-slices. Depending on whether or not data from the
characteristic domain couples back to the method used to treat the
central region, this is called either
Cauchy-Characteristic-Matching (CCM) or
Cauchy-Characteristic-Extraction (CCE). For details see the
review <cit.>. Recent numerical work can be found
in <cit.>. Well-posedness and
numerical convergence analysis of CCE and CCM setups in Bondi-like
gauges can be found in <cit.>.

Another proposal, pioneered by Friedrich <cit.>, is
to foliate spacetime via hyperboloidal slices. These are by definition
spacelike hypersurfaces that terminate at null infinity rather than
spatial infinity like Cauchy slices. For numerical applications
hyperboloidal slices can be combined with a compactified radial
coordinate. Following <cit.> this strategy is now completely
standard for perturbative work in a range of applications, see for
instance <cit.>. The essential subtlety
in working with compactified coordinates, which do the work of
bringing ℐ^+ to a finite coordinate distance, is that they
introduce divergent quantities into the problem. Fortunately in the
asymptotically flat setting these can be off-set by the smallness
coming from decay near ℐ^+. The specific rates therefore
matter. For applications to full GR there are two broad approaches to
this regularization. The first is to introduce curvature quantities as
variables and then work with conformally related variables. This
ultimately leads to the conformal Einstein field
equations <cit.>, which have the advantage of complete
regularity at ℐ^+. Numerical work using the conformal EFEs
is reviewed in <cit.>,
see <cit.>. The second broad category is to
work with evolved variables involving at most one derivative of the
metric, which is more standard in NR. Unfortunately there is no known
formulation of this type that is completely regular. Instead we have
to cope with expressions that are formally singular, but which are
expected to take finite limits at ℐ^+. In this sense such
formulations exist as an edge-case for numerical applications. Key
contributions in this setting include those of
Zenginoğlu <cit.>, who uses harmonic coordinates
for full GR, those of Moncrief and
Rinne <cit.>, who offered a partially
constrained formulation with elliptic gauge conditions, of Bardeen,
Sarbach, Buchman and later Morales <cit.>, with
a frame based approach, and of Vañó-Viñuales and
collaborators <cit.>, who employ
variations of the popular moving-puncture gauge. All of these setups
use a conformally related metric but without making curvature an
evolved variable.

We turn now to give an overview of the approach we follow here, which
was proposed in <cit.>. We work with compactified
hyperboloidal coordinates x^μ, but with evolved variables
associated with a coordinate tensor basis _μ
and dX^ν, as one would obtain in the standard solution of
the Cauchy problem. The idea is to take the evolved variables to
include a rescaling that knocks out their leading order decay, and to
then arrive at equations of motion that are as regular as possible. We
work in the second broad category discussed above, introducing at most
first derivatives of the metric as evolved quantities. The generalized
harmonic gauge (GHG) formulation is among the most popular in use in
NR. It is symmetric hyperbolic, possessing a very simple
characteristic structure with speed of light propagation. We thus take
the uppercase coordinates X^α to be generalized
harmonic, so that X^α=F^α, with gauge
source functions F^α which we can choose freely. Working
in the heuristic asymptotic systems setting of
Hörmander <cit.> as applied to great effect in the proof
of nonlinear stability of Minkowski <cit.>, it was found that
specific constraint addition <cit.> and choices for the gauge
source functions <cit.>
can be expected to improve the leading asymptotic decay of solutions
to GHG within a large class of initial data.

In the GHG formulation the metric components satisfy a system of
coupled nonlinear wave equations. Therefore, to assess the numerical
feasibility of our approach, we previously studied model systems of
both linear and nonlinear type, focusing in particular on the GBU and
GBUF systems, which were constructed to capture the asymptotic leading
behavior of GR in GHG. Promising numerical results have been presented
both in spherical symmetry <cit.> and full
3d <cit.>. Therefore, here we move on to give a thorough
treatment of spherical GR.

In section <ref> we give details of our geometric
setup, the formulation of GR that we employ, and our approach to
solving for constraint satisfying initial data. Afterwards, in
section <ref>, we briefly discuss our numerical
implementation and then present a suite of hyperboloidal evolutions of
full GR, placing particular emphasis on convergence tests. Results are
given with `pure' hyperboloids and with hyperboloidal
layers <cit.>, for both a height-function and an eikonal
approach to the construction of the hyperboloidal slices. Our
evolutions include various different types of physical initial data,
including pure gauge waves, constraint violating and satisfying data,
spacetimes with and without scalar field matter, perturbations of the
Minkoswki and Schwarzschild spacetimes. For the latter we examine both
quasinormal modes (QNMs) and late time power-law tail decay. We also
present results for initial data that start from a regular center and
then collapse to form a black hole. Section <ref>
contains our conclusions. Latin indices are abstract. Unless otherwise
stated, underlined Greek indices refer to the X^μ basis,
while standard Greek refer to that of x^μ. The metric is taken to
have mostly + signature. Geometric units are used throughout.




§ GEOMETRIC SETUP AND THE EINSTEIN FIELD EQUATIONS




We work in explicit spherical symmetry with spherical polar
coordinates X^μ = (T,R,θ,ϕ). As we develop
our formalism we implicitly assume that in these coordinates the
metric asymptotes to the standard form of Minkowski in spherical
polars near both spatial and future null infinity. Such coordinates
are guaranteed to exist by any reasonable definition of asymptotic
flatness, but for now we do not try to establish the slowest possible
decay that could be dealt with, and instead focus on the inclusion of
a large class of physical spacetimes. The first two variables, C_+
and C_-, are defined by requiring that the vectors

    ξ^a = _T^a + C_+ _R^a ,    ξ^a = _T^a+C_-_R^a

are null. Next, the function δ is defined by demanding that the
covectors

    σ_a = e^-δξ_a ,    σ_a = e^-δξ_a

are normalized by

    σ_a _R^a = -σ_a _R^a = 1  .

Finally, we define the areal radius

    R≡ e^ϵ/2 R   .

Due to spherical symmetry, all
variables { C_+, C_-, δ, ϵ} are functions of (T,R)
only. With all these elements, the metric takes the form

    (g_μν) = (
    [  2 e^δ C_+ C_-/C_+ - C_- e^δ(C_- + C_+ )/C_- -C_+                        0                        0; e^δ(C_- +C_+ )/C_- - C_+          2 e^δ/C_+ - C_-                        0                        0;                        0                        0                      R^2                        0;                        0                        0                        0              R^2 sin^2 θ;                          ])  .

As usual we denote the Levi-Civita derivative of g_ab
by ∇_a.

The metric is split naturally as

    g_ab = g_ab + g_ab  ,

where g_ab is the {T,R} part of the metric
and g_ab is the metric defined on a sphere of radius R
at time T. These metrics are converted into projection operators
onto their respective subspaces by raising one of their indices by the
inverse metric g^ab.

Our variables admit a simple interpretation. C_± are the local
radial coordinate lightspeeds in coordinates X^μ. The
variable δ determines the determinant of the
two-metric g_ab in these coordinates through

    g = e^2δ ,

and ϵ paramaterizes the difference between the coordinate and
areal radii.

In spherical symmetry the stress-energy tensor reduces to

    ( T_μν ) = (
    [         T_TT         T_TR            0            0;         T_TR         T_RR            0            0;            0            0         T_θθ            0;            0            0            0 T_θθ sin^2 θ;              ])   ,

with (T_TT,T_TR,T_RR,T_θθ) being functions
of (T,R) only. In trace reversed from, the field equations are

    R_ab = 8 π( T_ab - 1/2 g_ab T_c^c )  ,

where R_ab is the Ricci tensor of g_ab. Defining D_a as the
covariant derivative associated with g_ab, and
contracting equation (<ref>) with our null-vectors and tracing
in the angular sector, we get 

    - D_σ D_σR
        + D_σR(D_σ - D_σ) C_+/κ
        = 4 π R  T_σσ  , 
        - D_σ  D_σR
        + D_σR  (D_σ - D_σ) C_-/κ
        = 4 πR  T_σσ  , 
       12_2 δ - D_a[e^δ/κ^2(σ^aD_σ C_- 
        - σ^aD_σC_+) ]
        + 2/R^3M_MS
       
        + e^δ/κ^3[ D_σC_+ D_σ C_-
        - D_σ C_+ D_σ C_- ] 
        = 8 π  T_θθ/R^2
        + 8 πe^δ/κT_σσ  , 
       _2R^2 - 2
        + 16 πe^δ/κR^2 T_σσ 
        = 0   ,

where we use the notation D_σ≡σ^aD_a
and D_σ≡σ^aD_a to denote directional
derivatives along the null-vectors σ
and σ. Likewise subscripts σ and σ
denote contraction with these vectors on that slot of the respective
tensor. For the two-dimensional d'Alembert operator in the TR plane
we write _2≡g^abD_aD_b and define the
shorthand κ≡ C_+ - C_-. Finally, the Misner-Sharp
mass <cit.> is given by

    M_MS≡12R(
      2e^δ/κ (D_σR) (D_σR) + 1 ) .

The difference between the projected d'Alembert operator and the
full 3+1 dimensional version, defined
by ≡ g^ab∇_a∇_b, is

    (-_2)φ = - 2/Re^δ/κ( D_σRD_σφ
      + D_σRD_σφ)

when acting on a spherically symmetric function φ. We observe that
the field equations are highly structured when expressed in these
variables. For instance, null directional derivatives D_σ
and D_σ of C_± always appear with a κ^-1
prefactor and, taking this into account, outside of the principal part
the variable δ appears only in the
combination e^δ/κ. Regularity at the origin is discussed
below.

Due to spherical symmetry the two radial null-vectors σ^a
and σ^a must be tangent to outgoing and incoming geodesic
null-curves. They satisfy

    D_σσ^a   =(D_bσ^b)σ^a
      =κ^-1[(D_σ-D_σ)C_+]σ^a ,
    
      D_σσ^a   = (D_bσ^b)σ^a
      = κ^-1[(D_σ-D_σ)C_-]σ^a .


We shall consider a minimally coupled massless scalar field as the
matter model, whose equation of motion is

    ψ≡ g^ab∇_a∇_bψ = 0   .

The null components of the stress energy tensor for this matter
content are given by

    T_σσ = (D_σψ)^2  ,
         T_σσ
      = (D_σψ)^2  ,
      
       T_σσ = 0  ,
         T_θθ = e^δ/κR^2
      D_σψ D_σψ .





 §.§ Generalized Harmonic Gauge




In GHG the coordinates satisfy wave
equations X^α=F^α. In practice this is
imposed by rewriting these wave equations and defining constraints
that measure whether or not they are satisfied. This results in

    C^μ≡Γ^μ + F^μ =
    0 ,

which we will refer to as GHG or harmonic constraints,
where Γ^μ = g^νλ Γ^μ_νλ are the contracted
Christoffels with

    Γ^μ_νλ
      = 1/2 g^μρ
      (_νg_ρλ
      + _λ g_νρ
      - _ρ g_νλ)   ,

and F^μ's are the gauge source functions, and we recall
that X^α=-Γ^α. When the GHG
constraints are satisfied throughout the evolution, the EFEs are
equivalent to the reduced Einstein equations (rEFEs),

    R_ab - ∇_(a C_b) + W_ab = 8π(
      T_ab - 12 g_ab T_c^c ) ,

where the constraint addition tensor W_ab=W_(ab)(C_c) may be any
tensor constructed from the harmonic constraints with the property
that W_ab(0)=0, so that constraint propagation is maintained. As
usual, curved parentheses in subscripts denote the symmetric
part. With this adjustment, metric components satisfy nonlinear
curved-space wave equations. In spherical symmetry, the only free
components of F^μ are F^T and F^R, or equivalently the
null components F^σ≡ F^aσ_a
and F^σ≡ F^aσ_a. The angular components
of these constraints are satisfied identically, provided that we
choose

    F^θ = R^-2θ  ,    F^ϕ = 0  .

The null components of these constraints then read

    C^σ   ≡ C^aσ_a = F^σ
                   +2D_σC_+/κ
                   -2D_σR/R
    
      C^σ   ≡ C^aσ_a = F^σ
                        -2D_σ C_-/κ
                        -2D_σR/R

The gauge source functions F^σ and F^σ will later
be used to help impose asymptotic properties of solutions to the rEFEs
towards ℐ^+, but we can already keep in mind
that F^σ≃ F^σ≃ 2/R, so that the first
and third terms in the right-hand-sides of these equations cancel each
other to leading order near ℐ^+. Demanding regularity of
the rEFEs at the origin will furthermore restrict the limit of these
functions at the origin.

In GHG the metric components would be naively expected to decay like
solutions to the wave equation. Incoming null-derivatives of C_+
and ϵ therefore ought to decay at best
like R^-1 as we head out to ℐ^+. The
harmonic constraints (<ref>), however,
assert that they should in fact be equal to terms that decay
faster. Following <cit.>, who used asymptotic expansions, it
is possible to include constraint additions in the rEFEs so that even
when the constraints are violated, these two specific incoming
derivatives are expected to decay more rapidly. Making such constraint
addition and subsequently redefining the constraint addition
tensor W_ab, we can write the rEFEs as,

    D_σ(2/κR^2D_σC_+ )
        +RD_σ(R F^σ)
        -D_σR^2D_σ C_+/κ
        -R^2 W_σσ 
      = -8πR^2T_σσ , 
        D_σ(2/κR^2D_σ C_-)
        -RD_σ(RF^σ)
        -D_σR^2D_σ C_-/κ
        +R^2 W_σσ 
      = 8πR^2T_σσ  , 
       _2δ + D_a(g^a_bF^b)
        + 2e^δ/κ^3[ D_σC_+ D_σ C_-
        - D_σ C_+ D_σ C_- ]  
        +2/R^2(1-2M_MS/R)
        +2e^δ/κW_σσ
        = 16 π  T_θθ/R^2  , 
       _2R^2 - 2 -2R^2W_θθ
        = -16πe^δ/κR^2 T_σσ  .


When performing a fully first-order reduction of these equations in
terms of null derivatives, a commutator between σ^a
and σ^a derivatives has to be calculated. The action of
this commutator acting on a general spherically symmetric function f
is

    [σ,σ]^a    =1/κ
       (D_σ C_- - D_σC_+)(σ^a - σ^a ) 
         +(D_σδ)σ^a
       -(D_σδ) σ^a   .





 §.§ Regularizing the origin




When describing a spacetime with regular origin in spherical-like
coordinates the metric components have a well-defined parity. Diagonal
components are even functions of R, whereas TR components are
odd. Translating these conditions to our variables we find
that C_+(T,R) + C_-(T,R) is an odd function of R,
whereas C_+(T,R) - C_-(T,R), δ(T,R) and ϵ(T,R) are
even. These parity conditions imply

    C_± (T,-R) = -C_∓ (T,R)   .

Derivatives inherit parity in the obvious way. Regular initial data
moreover require

    lim_R → 0[ ϵ(T,R) - δ(T,R)
      + ln (κ(T,R)/2) ] = O(R^2)   .

Since all the terms in the above limit are even functions of R, this
limit is automatically satisfied by imposing the condition

    ϵ(T,0) = δ(T,0) - ln(κ(T,0)/2)   .

In other words, this condition says that the 0th order term in
Taylor's expansion of ϵ at the origin should satisfy the
above condition.

All these parity conditions assure regularity of the EFEs at the
origin if the initial data are smooth there. To ensure a similar
regularity for the rEFEs, we need to choose F^μ such
that C^μ and its first derivatives remain regular
there. By definition, F^T and F^R are even and odd functions
of R, respectively. From the expressions of the
constraints (<ref>), it can be seen that
the necessary and sufficient conditions on constraints for a regular
origin are that F^T should be regular and F^R must take the
leading limit 2e^-ϵ/R at the origin, with an additional
regular odd function permitted. These conditions translate to the null
components as

    F^σ≃2e^-ϵ/2/R ,
         F^σ≃ -2e^-ϵ/2/R ,

near the origin.

The constraint addition tensor as defined
in (<ref>) also needs to be regular at the
origin. Since the rEFEs are already regular with the previous choices,
and the purpose of constraint addition is to regularize the
asymptotics, we just take a sufficiently rapidly vanishing constraint
addition tensor at the origin. With the adjusted definitions
of W_ab in (<ref>), this corresponds instead
to

    W_σσ   = (e^-δ_R C_+ 
     + D_σ(ϵ +2ln R))C^σ , 
    
     W_σσ   = (e^-δ_R C_- 
     + D_σ(ϵ +2ln R))C^σ , 
    
     W_θθ   = Re^δ/κ R^2(D_σR C^σ + D_σR C^σ ) ,
    
      W_σσ   = 0  .





 §.§ Gauge Sources and Regularization
  at ℐ^+



Below we will change from the generalized harmonic
coordinates X^α to compactified hyperboloidal
coordinates x^α. Pushing the field
Eqs. (<ref>) through this change will only result
in a set of PDEs that can be treated numerically if solutions decay
fast enough. As discussed above, the asymptotic decay of C_+^R
and ϵ can be influenced by adding suitable combinations of
the harmonic constraints to the field equations,
see <cit.>. These have been incorporated
into (<ref>), so that

    W_σσ = W_σσ =
      W_θθ = W_σσ = 0 ,

already includes damping terms in the field equations that, at least
within a large class of initial data, should result in decay
like D_σC_+=D_σϵ=O(R^-2)
near ℐ^+. Ideally, we would like to obtain similar
improved decay (beyond that expected for the wave equation) for the
variables C_- and δ. The remaining tool we have to achieve
this is to use the gauge source functions F^σ
and F^σ. Following the approach of <cit.>, we take,

    F^σ   = 2/R 
      + 2p/R(e^δ -1)  ,
    
      F^σ   = -2/R 
      -p/R(1+C_-)
      + 1/Rf_D .

In spherical vacuum the choice p=0 and f_D=0 grants δ
and C_- asymptotic decay like solutions to the wave equation,
whereas p=1 should give improved decay on incoming null derivatives
like that of C_+ and ϵ. Once we introduce the scalar field
the situation is more complicated because, without care, the slow
decay of T_σσ results in poor decay for the
incoming lightspeed
like C_-∼-1+(lnR)/R. In 3+1 dimensions
without symmetry gravitational waves induce similar behavior.
Fortunately, this shortcoming of plain harmonic gauge can be overcome
by using the gauge driver function f_D to absorb the logarithmic
terms. Model problems for this have been studied both in spherical
symmetry and in full 3d <cit.>. In particular, we take the equation of motion

    f_D -2/χ(R)∂_T f_D 
     - 32 π (∂_T ψ)^2 = 0  ,

for the gauge driver f_D. The basic idea is that by insisting on a
wave-equation principal part hyperbolicity is guaranteed whilst
simultaneously the second term suppresses the natural radiation field
associated with the wave operator, and the third forces it equal to a
desired value that eradicates the slowest decay in the worst behaved
of the Einstein equations (the wave equation for C_-
in (<ref>)). Details can be found in the
references above. Here we have defined χ≡√(1+R^2), so
that χ(R) is an even function of R, χ∼ 1 near the
origin and χ∼ R near ℐ^+. When starting from black
hole initial data we adjust slightly the
choice (<ref>) based on compatibility
with the Schwarzschild solution in a reference coordinate system. The
specifics are explained below.

We have described separately how the choice of gauge and constraint
addition have to be taken in order to have regularity at the two
potentially problematic ends, the origin and ℐ^+. To
transition smoothly from the origin choice (<ref>) to the
asymptotic choice (<ref>) in applications we multiply the
former (origin choice) by a function ι(R) that is
identically 1 in an open region containing the origin and decays as
a Gaussian asymptotically, namely 

    ι(R) = 
     
      1 ,      R<R_0 
    
      e^-((R-R_0)/σ_0)^4 ,     R≥ R_0  .

We multiply the latter (asymptotic choice) by 1-ι(R).




 §.§ First order reduction and rescaling




In our numerical implementation we use a first order reduction of the
field equations. For this we introduce the first order reduction (FOR)
variables

    θ^±≡D_σ C_±/κ  ,
        θ^±≡D_σ C_±/κ ,
        ζ^+ ≡ D_σ ζ  ,
        ζ^- ≡ D_σ ζ ,

where ζ stands for either δ, ϵ, ψ or f_D.

Parity conditions for the FOR variables are easily obtained from the
definitions (<ref>) and the parity conditions of
the metric components described above, plus the fact that the scalar
field ψ and the gauge driver f_D are even functions
of R. These conditions are

    θ^± (T,-R) = - θ^∓ (T,R)   , 
      θ^± (T,-R) = - θ^∓ (T,R)   , 
       ζ^+(T,-R) = ζ^-(T,R)   ,   ζ^-(T,-R) = ζ^+(T,R)   .

Time derivatives of these variables satisfy the same parity
conditions, whereas radial derivatives flip the signs on the right
hand sides.

Introducing these variables creates new constraints. According to our
definitions, they read

    C_C_±≡_R C_± 
        - e^δ(θ^± - θ^±)   , 
        C_ζ≡_R ζ - e^δζ^+ - ζ^-/κ  .


As discussed earlier, at least within a large class of initial data, a
particular behavior of the variables is expected as ℐ^+ is
approached. Therefore, in order to obtain O(1) variables throughout
the whole domain, we rescale the evolved functions according to their
expected asymptotic decay. We define

    C̃_±≡χ (C_±∓ 1)  , 
          Θ^±≡χ^2θ^± , 
          Θ^±≡χθ^± , 
       Z ≡χζ , 
          Z^+ ≡χ^2 ζ^+  , 
          Z^- ≡χζ^-  ,

where the ∓ 1 in the C̃_± is their Minkowski value, and
so it is only this difference that decays. In words, the rule for the
definitions is that rescaling by the given powers of χ maps from
the “ζ” variables and their null derivatives to the
capitalized variables “Z”. In our numerical implementation all of
the evolution equations and constraints are written in terms of these
rescaled variables.

For brevity we do not state the full symmetric hyperbolic equations of
motion for the rescaled reduction variables, but they are
straightforwardly derived, and can be found in the Mathematica
notebooks <cit.> that accompanies this
paper. Instead, to illustrate the procedure and the basic shape of the
expressions obtained in the reduction, in particular
at ℐ^+, let us consider a variation of the `ugly' model
equation employed in earlier work, namely

    φ + 2 p R^-1(D^aR)D_aφ   =S ,

with p a non-negative integer and where the source term S may be
thought of as decaying like R^-3
near ℐ^+. (The constant p
in (<ref>) corresponds to the value
of p in (<ref>) for the variables δ and C_-). As
elsewhere, we assume sphericity, but now with an arbitrary given
asymptotically flat metric using the same notation as above. We
introduce the rescaled variable =Rφ, and
define rescaled reduction variables

    ^+=R^2D_ξφ ,  ^-=RD_ξφ .

This gives rise to the reduction constraint

    C_φ   =_Rφ-κ^-1e^δR^-1
                 [R^-1^+-^-] .

The equations of motion for the reduction variables can then be
written as

    D_σ^- + pR^-1(D_σR)^-
        +(1+p)R^-2(D_σR)^+
       -1κ(D_σC_+-D_σC_-)^-
        + κ2RS = 0 ,
       D_σ^+
        - (1-p)R^-1(D_σR)^+
        +(1+p)DR^-
       -1κ(D_σC_+-D_σC_-)^+
        +κ2R^2S = 0 .

Introducing the change to hyperboloidal coordinates we effectively
have to multiply the first of these equations by a power R'∼ R^n,
with n∈(1,2]. If the fields have wave equation asymptotics (p=0)
and S=O(R^-3), this creates no formally singular
terms. If instead p is taken to be a positive integer, suitable
initial data allows for ^- to decay faster
than O(R^-1), but this still renders the second term in
the first of these equations formally singular. With our first order
reduction, the EFEs are similar to this system with either p=0 or
with p=1, so that all formally singular terms appear with an
identical structure and can be straightforwardly treated by
L'Hôpital's rule.




 §.§ Hyperboloidal Foliations and the
  Dual-Frame Formalism




The purpose of the present work is to perform the numerical evolutions
within hyperboloidal slicings of spacetime. To do so, we need to
introduce the appropriate changes of coordinates that define them.
Following the dual-foliation strategy of <cit.>,
this change is made without transforming the tensor basis we use. In
our setting, this gives greater freedom in the choice of coordinates
but without interfering with hyperbolicity, and ultimately permits us
to include ℐ^+ within the computational domain.

The first step in our construction is to introduce a compactification
so that we can bring R→∞ to a finite coordinate
distance. This is accomplished by defining a new coordinate r
through

    R(r)    = r_m + r-r_m/Ω(r)^1/n-1Θ(r-r_m)
               , 
    Ω(r)    =1-(r-r_m)^2/(r_ℐ-r_m)^2  ,
                    1< n ≤ 2   ,

where Θ is the Heaviside function, so
that R→∞ corresponds
to r→ r_ℐ. For simplicity we
take r_ℐ = 1. The compactification is the identity in
the range r∈ [0,r_m]. Note that R is a monotonically increasing
function of r, and so it is invertible. The derivative diverges
asymptotically at the rate dR/dr≡ R'∼ R^n, so that the
parameter n serves to control the rate of
compactification <cit.>. By definition a hyperboloidal
slice is one which remains spacelike everywhere but terminates at null
infinity. We construct our hyperboloidal time coordinate according to
two strategies, which we explain in the following subsections.




  
*The height function approach: The first method we use to
construct hyperboloidal time employs a height function. Here, level
sets of the time coordinate are explicitly lifted up by a given
function in a spacetime diagram relative to those of the harmonic time
coordinate T, which in contrast are taken to terminate at spatial
infinity. For this we define a time function that asymptotes to
retarded time through

    t = T - H(R)  ,

where H(R) is called the height function because it encodes how the
new slices are lifted. Recalling the expected
asymptotics D_σ C_+= D_σ C_+ = O(R^-2)
of the outgoing coordinate lightspeed, it follows that the associated
`mass-term' at ℐ^+

    C̃_+|_ℐ^+ = m_C_+ ,

is constant. Imposing that the outgoing coordinate lightspeed c_+ in
the hyperboloidal coordinates is bounded <cit.> requires knowledge of this mass-term. We therefore
choose

    H(R) = R - m_C_+ln R - r  .

Observe that the term involving m_C_+ mimics the tortoise
coordinate of Schwarzschild spacetime asymptotically. The appearance
of m_C_+ here is the reason that we need improved asymptotic decay
in C_+, since if the mass-term were time dependent the
height-function ansatz as defined in (<ref>) would
fail.

If we take initial data such that C_+ decays sufficiently fast, in
particular so that m_C_+≡ 0, we recover the definition used
in earlier numerical works to treat systems of linear and nonlinear
wave equations in the Minkowski spacetime <cit.>. The inclusion of the term m_C_+ here
is important both to guarantee that the outgoing radial coordinate
lightspeed in the compactified hyperboloidal coordinates is O(1),
but also, as illustrated in Figure <ref>, to obtain the
desired global structure of the slices.
  



 
With this construction, the Jacobians to change from (T,R)
to (t,r) derivatives are given by

    _R  = 1/R'(r)_r 
      + ( 1 - 1/R'(r) 
      - m_C_+/ R(r) ) _t    , 
       _T = _t  ,

where we see that the mass-term places a correction that dominates the
term coming from the compactification itself.




  
*The eikonal approach: Our second approach to construct a
hyperboloidal time coordinate is by defining

    t = u + r

where u satisfies the eikonal equation

    ∇^cu∇_cu = 0  .

As explained in <cit.>, demanding
that ∇^a u ∝ξ^a leads to having outgoing coordinate
lightspeeds identically one in the compactified hyperboloidal
coordinates. This means that hyperboloidal slices built this way adapt
dynamically so that we can control outgoing speeds and should help
avoid undesirable coordinate red and blue-shift effects on outgoing
signals.

The idea in the eikonal approach is to derive
equation (<ref>) and project it
along σ^a to get an evolution equation
for U^-≡∇_σu, which can then be solved
alongside the rEFEs, while the
condition -∇^a u ∝σ^a gives a fixed functional form
to U^+≡ D_σ u in terms of U^- and the metric
variables. The equation we get for U^- is of advection-type, with
principal part decoupled from the principal part of the rEFEs, so
symmetric-hyperbolicity of the composite system is trivially
preserved.

Similarly to the previous sections, the choice we take for the
function u that satisfies the eikonal equation can only be made
asymptotically, as parity of u at the origin is complicated and in
any case is incompatible with black hole excision. The reason for this
is that since the eikonal coordinates forces the outgoing lightspeeds
to be identically one, we necessarily need boundary conditions at the
black hole boundary. To overcome these obstacles, we place a
source S in the right-hand-side of equation (<ref>),
instead solving

    ∇^cu∇_cu = S  ,

and choosing S so that u ≃ T-R near the origin/horizon and
satisfies the eikonal equation identically only
near ℐ^+. This reduces the Jacobians to the identity
when R is small whilst allowing them to take the desired form
near ℐ^+. (See <cit.> for full details).

Concretely, if we decompose the vector u^a≡ -∇^a u with
the identity

    u^a = e^δ/κ ( U^- σ^a + U^+ σ^a )

the evolution equation for U^- reads

    u^a D_a U^- = ( u^a D_a σ^b )u_b 
     + 12D_σS .

When the expressions on the right-hand-side are expanded out we see
that they in fact contain incoming null derivatives
of C_+. Rewriting this equation in the compactified hyperboloidal
coordinates requires multiplying by R^n, so the improved asymptotic
decay of C_+ is needed just as in the height function setting. The
sourced eikonal equation (<ref>) places a constraint
for U^+ of the form

    U^+ = -κ/2U^-e^-δS  .


The Jacobians in this setting that change from (T,R) derivatives
to (t,r) ones are

    _R    = 1/R'(r)_r + ( e^δ/κ (U^+ - U^-)
       + 1/R'(r)) _t   , 
    _T    = e^δ/κ
       (C_+ U^- - C_- U^+)_t   .





 §.§ Bondi Mass




In the spherical setting, for each hyperboloidal slice, the Bondi
  mass <cit.> can be found by taking the limit of the
  Misner-Sharp mass as one tends to ℐ^+. From our
  expression for M_MS, eq. (<ref>), this
  leads to the expression

    M_B   ≡lim_R→∞ M_MS
       = lim_R→∞1/4 RE^-
                     +1/8( 2C̃^- -2C̃_+
                     + 4Δ -2E^+ 
          +(C̃_+ + C̃_-)E^-  + E^-E^+
                     - 4E +3E^-E  )

Recalling that ϵ has improved asymptotic decay, the first
term is formally singular, meaning it attains a finite limit by a
product of a term that diverges, R, and a term that decays at least
as O(R^-1), namely E^-. As opposed to previous formally singular
terms appearing in the equations of motion, the evaluation of this
term by use of L'Hôpital's rule is cumbersome. However, since at the
continuum level all physical quantities are defined up to constraint
addition, and noting that precisely this reduction variable appears in
the σ^a component of the GHG constraints
(eq. (<ref>)), the expression
for M_B can be regularized by a constraint addition that
asymptotes at leading order to -R^2 C^σ/4.  With
this particular choice the regularized expression for M_B
reads instead

    M_B = 1/4( -C̃_+ -C̃_- +F_D
      - E^+ -2Θ^- )

where clearly all terms are now regular O(1) at ℐ^+. It
is this expression that will be
used when we evaluate M_B for our numerical simulations.

The Bondi mass has to satisfy two requirements for a physical solution
of the EFEs. First it has to be non-negative, and second it should be
monotonically decreasing as radiation leaves the spacetime
through ℐ^+. This last property can be deduced at the
continuum level by the Bondi mass-loss formula, which in terms of our
matter model, variables, our particular constraint addition and
choice of gauge reads

    Ṁ_B = -π (Ψ^-)^2  .





 §.§ Constraint Satisfying Initial Data




The evolution Eqs. (<ref>) are equivalent to
the EFEs only when the GHG constraints are identically satisfied over
the entire domain. As these constraints satisfy a system of coupled
second-order nonlinear homogeneous hyperbolic PDEs <cit.>, they remain satisfied throughout the evolution if they
are satisfied in the initial data up to their first time
derivative. To obtain this formal evolution system satisfied by the
harmonic constraints we just have to take a divergence of the
rEFEs (<ref>). Being a free evolution scheme,
the GHG formulation of the EFEs provides us with no other way to
ensure that these constraints are satisfied throughout the evolution.

Recall that we employ two sets of coordinates,
namely X^α, which are used to define the tensor basis
for the GHG formulation, and the compactified hyperboloidal
coordinates x^α. To avoid confusion in our discussion of the
constraints it is therefore most convenient to rely on the rEFEs in
abstract index form as
in (<ref>). Trace-reversing the rEFEs and
contracting once with n^a, the future pointing unit normal to our
hyperboloidal (constant t) slices gives

    ∇_n C_a    = 2 ℳ_a
      + (n_a γ^bc - n^c γ_a^b) ∇_b C_c
         + 2 W_abn^b-n_aW  ,

with ℳ_a=(G_ab-8π T_ab)n^b similar to the
expression given in <cit.>, W the trace of the
constraint addition tensor W_ab and γ_ab the spatial
metric. Here C_a is the one form of the GHG constraints defined in
the X^α tensor basis by (<ref>),
and ℳ_a = 0 encodes the Hamiltonian and momentum
constraints associated with the hyperboloidal foliation. From this
expression we observe that it is sufficient to choose initial data
that satisfy the harmonic constraints together with data that satisfy
the standard Hamiltonian and momentum constraints as usual. From this
we then need to reconstruct the evolution variables employed in our
formulation.

Given the spatial metric and extrinsic curvature associated with t
in the initial data, a general strategy for the latter step can be
formulated within the language of <cit.>. In broad strokes, this
would involve constructing appropriate projections of the Jacobians
that map between the two coordinate tensor bases, choosing the lapse
function of the uppercase T foliation, together with its
Lie-derivative along the uppercase normal vector to the T-foliation,
and then combining these quantities to build the uppercase spatial
metric and extrinsic curvature. From there we could choose the
uppercase shift vector and switch to our desired choice for the
evolved variables. This warrants a careful treatment without symmetry,
but for now we settle on a bespoke spherical approach sufficient for
our present needs.

We begin with a collection of useful expressions. First we express the
spatial metric and extrinsic curvature associated with the time
coordinate t in terms of our variables. To accomplish this, we
treat t as a general time coordinate that defines a foliation
and r to be a radial coordinate on those slices. Later, we shall
take the special case of hyperboloidal coordinates. The uppercase
coordinates (T,R) are then taken as

    T ≡ T(t,r)   ,    R ≡ R(r)   .

The metric in the lowercase coordinates can then be defined in terms
of the Jacobians J_μ^μ≡_μ X^μ as

    g_μν = J_μ^μ J_ν^ν  g_μν  ,

with g_μν given
in (<ref>). Denoting t and r derivatives by
dot () and prime ('), respectively, and the future directed
unit normal to the constant t hypersurfaces by n^a, the
standard ADM variables lapse, shift and spatial metric, denoted in the
standard notation, are expressed in terms of the GHG ones as,

    α = e^δ/2  R'  Ṫ √(κ)/√(2) √(R' - C_- T') √(R' - C_+ T')  ,
        
       β^r =  - Ṫ( (C_+ + C_-)
          R' - 2 C_+ C_- T' )/2   (R' - C_- T')
          (R' - C_+ T')  , 
       γ_ij = (
    [              γ_rr                 0                 0;                 0         R^2   e^ϵ                 0;                 0                 0 R^2   e^ϵsin ^2 θ;                   ])   ,

with

    γ_rr = 1/γ^rr
        = 2 e^δ/κ(R'-C_- T')
      (R'-C_+ T')   .

Similarly, the extrinsic curvature, given
by K_ij = - _n γ_ij/2, takes the form

    K_ij = (
    [          K_rr             0             0;             0          K_θθ             0;             0             0 K_θθ sin ^2 θ;               ])  ,

with

    K_rr≡β^r  γ_rr' + 2   (β^r)'  γ_rr - γ̇_rr/2 α  , 
        K_θθ≡R^2   e^ϵ/α( R'  β^r/R + (β^r  ϵ' - ϵ̇)/2)   .

Plugging (<ref>) into the former of these
expressions results in an evolution equation involving our variables.
We also define the normal and radial derivatives for the massless
scalar field ψ described above as

    ψ_n ≡_n ψ = n^μ_μψ  ,
        ψ_r ≡_r ψ  ,

and the corresponding charge and current densities by

    ρ_ψ = 1/2( ψ_n^2 + γ^rr ψ_r^2 )   ,
         j^r_ψ = γ^rr ψ_n  ψ_r   .


The nontrivial components of the vector ℳ_a are then

    ℋ≡ 2  ℳ_n ≡
      2 n^a n^b ( R_ab - 1/2 g_ab R - 8 π T_ab)   ,

and

    ℳ_r    ≡ n^a γ_r^b ( R_ab 
      - 1/2 g_ab R - 8 π T_ab)   .

In terms of the above variables,
and ^γ K_rr≡ K_rrγ^rr
and ^γ K_T ≡ K_θθγ^θθ,
the Hamiltonian and momentum constraint equations are expressed as

    2  ^γ K_T^2 + 4  ^γ K_rr ^γ K_T
        - 16 πρ - 2 (R'   (γ^rr)'+ 2   R” γ^rr)/R
        + 2 (e^-ϵ - R'^2  γ^rr)/R^2
        - ((γ^rr)'
        + 6   R'  γ^rr/R) ϵ'
        - 3/2γ^rr ϵ'^2 
        - 2  γ^rr ϵ” = 0  ,

and

    - 8 π  j_r + 2   R'/R
        (^γ K_T - ^γ K_rr)
        + (^γ K_T - ^γ K_rr)  ϵ' 
        + 2  ^γ K_T' = 0   ,

respectively. Along with the GHG constraints these are the equations
that need to be satisfied by initial data to ensure we obtain
solutions of GR. Our bespoke procedure for spherical initial data is
divided into three steps:


  Step 1: Reformulate the Hamiltonian and momentum constraints,
  eqs. (<ref>) and (<ref>). Here, we work
  assuming that the initial matter distribution ψ_n and ψ_r
  and the Jacobian are given.

  Step 2: Formulate the GHG
  constraints (<ref>). At this stage one
  can choose initial data for the gauge driver f_D defined
  in (<ref>).

  Step 3: Make a concrete choice for the foliation and solve the
  constraints. Here, we can specify T(t,r) at t = 0 to construct
  the initial slice of interest, which is hyperboloidal in our
  case. This way, we keep the entire procedure general and keep the
  first two steps agnostic to the choice of foliation.


We see that both Hamiltonian and momentum constraint equations are
formally singular at the origin and at ℐ^+. To regularize
these equations, we extensively use the fact that spherically
symmetric vacuum solutions, like Minkowski and Schwarzschild, are the
solutions to these equations. To facilitate this, we shall take the
general solution to a nontrivial matter distribution to be corrections
over these vacuum solutions. This idea has been extensively explored
in a series of works about formulations of the constraint
equations <cit.> with different PDE character.
  
The GHG constraint equations are much easier to solve due the fact
that they contain time derivatives of our evolved fields C_± and
which can simply be set. Trivial data for f_D leads to undesired
asymptotics for Ċ_-. Nonetheless, the latter can be fully
controlled by carefully choosing f_D in the initial data, which we
shall do while solving for the GHG constraints.

A significant advantage of this approach is that all these results
apply to any spacetime foliation. For our specific interests, we shall
restrict ourselves to the hyperboloidal foliations constructed via the
height function and eikonal approach described above, and study the
nonlinear perturbations of Minkowski and Schwarzschild backgrounds.




  §.§.§ Minkowski perturbations
 



For the Minkowski spacetime in global inertial time and standard
spherical polar coordinates, C_+ = - C_- = 1
and δ = ϵ = 0, all with vanishing time derivatives. This
gives a simplified form to the ADM variables defined in
Eqs. (<ref>)-(<ref>), that
we shall denote with the subscript “MK", for Minkowski, for
example α_MK. We will utilize these values to simplify
the general constraint equations.

As mentioned above, our first step is to solve the Hamiltonian and
momentum constraint equations. We choose to solve them
for γ^rr and ^γ K_T, respectively, as these
equations are linear in these variables. Assigning the Minkowski
values obtained above to ^γ K_rr and ϵ, and taking
the substitution

    ^γ K_T = ^γ K^MK_T
      + K_T^(1)/R  ,
         j^r = j^r_ψ  ,

reduces the momentum constraint equation to

    (K_T^(1))' = 4 π  R  ψ_n  ψ_r   .

To set ^γ K^MK_rr and ^γ K^MK_T
and so forth, we choose the Minkowski values for our evolved variables
just stated, then apply the Jacobian transformations given above. This
equation is regular whenever ψ_n  ψ_r ∼ 1/R in the
initial data. Thereby, the entire impact of the current j_r is now
contained in the correction term K_T^(1). Setting
moreover ψ_n = 0 and imposing the boundary
condition K_T^(1) = 0 at the origin, we obtain the trivial
solution for K_T^(1) in the initial data. Further, taking

    γ^rr = (γ_MK)^rr
      + γ^(1)/R'^2   R  ,
        ρ = ρ_ψ  ,

gives the simple form

    (γ^(1))'    =    - 4 π R^2 R'  ψ_r^2 
                        ( (γ_MK)^rr
                        + γ^(1)/R R')   ,

to the Hamiltonian constraint equation, which is regular
whenever ψ_r ∼ 1/R in the initial data. For this class of
initial data, all information of the initial matter
distribution ψ_r is now encoded in γ^(1). This is the
final equation we solve for generating the Hamiltonian and momentum
constraint satisfying initial data in the Minkowski case.

The second step involves solving the GHG
constraints (<ref>). We will solve them
for Ċ_+ and Ċ_- using the gauge source functions given
in (<ref>), along with all the
simplifying assumptions and solutions to the Hamiltonian and momentum
constraint equations given above. We further take p = 0 in the gauge
source functions and C_+ = - C_- = 1 in the initial data, which
translates to the ADM variables as

    α = ( α_MK √(γ_rr  (γ_MK)^rr)) 
       ,   β^μ = β^μ_MK  .

Now we observe that the trivial data for the gauge driver f_D gives
data for Ċ_- which is nonzero and O(1)
at ℐ^+. This inconsistency contradicts our regularization
scheme adopted from <cit.>. As previously noted, the data for f_D is
sufficient to establish the asymptotic behavior of Ċ_- in the
initial data to any desired level. For our convenience, we opt for a
choice that causes it to decay more rapidly than any inverse power
of R. One such choice is

    f_D = -    ( γ^(1)/R'^2 R
                - β_MK^r/Ṫ R'(1 - √(1 + (γ_MK)_rr γ^(1)/R R'^2)) )  ×
       2 χ(R) R' (R' + T')/R  ,

and Ḟ_D = 0. Interestingly, this data depends on the
foliation. This exhausts all our requirements to build the data. To
sum up, we have assigned Minkowski values to C_±, ϵ
and ^γ K_rr, which translates to δ̇. The
Minkowski values are determined by transforming the global inertial
values of our evolved variables through the appropriate Jacobians. We
must then solve the constraints to obtain Ċ_±, γ^rr
and ^γ K_T. The last two translate to δ
and ϵ̇ in the data. We have also established data for
the matter and the gauge driver while ensuring that the asymptotic
properties of the ADM variables are satisfied.

The final step is choosing the Jacobians. We consider both the height
function and eikonal approaches to define the foliation T(t,r) and
keep the same compactification R(r) defined above in both cases. For
Minkowski perturbations, we set m_C_+ = 0 in the height
function (<ref>). In our numerical
implementation we simply integrate out the resulting ODEs.

The Misner-Sharp mass corresponding to this class of initial data is

    M_MS^MK = -1/2γ^(1)  .

Interestingly, this expression for the mass and the rescalings in the
correction terms in (<ref>) and (<ref>) are
independent of the foliation. As we conclude from
Eq. (<ref>), and observe in
Figure <ref>, M_MS remains positive
semi-definite everywhere and O(1) whenever ψ_r falls-off at
least like 1/R towards ℐ^+. We will consider the more
general case with ψ_n  0 in the future.








  §.§.§ Schwarzschild perturbations




Our procedure for perturbed black hole initial data is very similar.
The Schwarzschild metric in Kerr-Schild coordinates is given by

    (g_SS)_μν = (
    [   -1 + 2 M/R        2 M/R            0            0;        2 M/R    1 + 2 M/R            0            0;            0            0          R^2            0;            0            0            0 R^2 sin ^2 θ;              ])   .

Here, M represents the mass of the black hole and the
subscript “SS" stands for Schwarzschild. Comparing this with
the metric in (<ref>), we observe that for the
Schwarzschild metric in Kerr-Schild coordinates, we have

    C_+ = 1 - 2M/R/1 + 2M/R  , 
         C_- = -1   , 
        δ = ϵ = 0   ,

all with vanishing time derivatives. This gives the associated ADM
variables that we denote by the subscript “SS".

Like before, we set the Schwarzschild values for ^γ K_rr
and ϵ in the initial data and take

    ^γ K_T = ^γ K^SS_T + (K_SS)_T^(1)/R

to get the momentum constraint equation similar
to (<ref>). Once again, we set ψ_n = 0 to get
the trivial solution (K_SS)_T^(1) = 0 in the initial
data. Likewise, taking

    γ^rr = (γ_SS)^rr 
     + (γ_SS)^(1)/R'^2   R  ,

gives an equation for (γ_SS)^(1) similar
to (<ref>). This equation is regular, as before, with
a solution that is regular throughout the domain. Thus we again embed
all information of the initial matter distribution
in (γ_SS)^(1) and solve it to generate the
Hamiltonian and momentum constraints satisfying initial data.

We again take the simplifying conditions

    α = ( α_SS √(γ_rr  (γ_SS)^rr)) 
        ,   β^μ = β^μ_SS  ,

to get the Schwarzschild values for C_± and solve the GHG
constraints for Ċ_± with the gauge source functions

    F^σ_SS = F^σ + C_+ - 1/R  ,

and F^σ given
in (<ref>) with p = 0. We,
additionally, take the following initial data for the gauge driver

    F^SS_D =     - ( γ^(1)/R'^2 R
      - β_SS^r/Ṫ R'(1 - √(1 + (γ_SS)_rr  (γ_SS)^(1)/R R'^2)) ) ×
          2 R' (R' + T') 
       ,

and ḟ_D = 0 to correct the asymptotics of Ċ_-. This
specific choice of f_D gives trivial data for Ċ_- and,
interestingly, reduces to (<ref>) for
vanishing M and χ(R) = R. This exhausts all our requirements to
generate the initial data.

As before, we consider both constructions of hyperboloidal foliation
defined above. We take m_C_+ = -4M in the height function.

The Misner-Sharp mass in this case is given by,

    M_MS^SS = M -1/2 (γ_SS)^(1)  ,

with (γ_SS)^(1) having similar properties as in the
Minkowski case, and effectively increasing the mass of the initial
data.




§ NUMERICAL EVOLUTIONS



Our numerical implementation lies within the infrastructure used in
earlier works, the most similar systems being <cit.>. The method itself is
entirely standard, so we give just a brief overview. Evolution is made
under the method of lines with fourth order Runge-Kutta. We use second
order finite differences in space. Our first order reduction makes the
treatment of the origin quite subtle. Following our earlier work we
have adapted Evans method <cit.> (see also <cit.>
for variations) in the obvious manner for the metric components (and
their reduction variables) as well as the scalar field. We define in
particular two second order accurate finite differencing operators as
acting on some given grid-function f,

    D_0f    = 1/hf_i+1-f_i-1/2 ,
    D̃f     =1/h (p+1) r^p_i+1f_i+1-r^p_i-1
                   f_i-1/r^p+1_i+1-r^p+1_i-1 ,

with h the grid-spacing. The parameter p here is not directly
related to those that appeared
in (<ref>)
and (<ref>)). Since we are focused in this work on
proof-of-principle numerics, we have not tried to extend the code
beyond second order accuracy. This remains an important task for the
future, but there is no particular reason to expect difficulties in
doing so. Terms in the evolution equations
like _rψ + prψ are treated with D̃,
whereas plain derivatives are approximated by D_0. As mentioned
above, ghostzones to the left of the origin are populated by
parity. With this said, the EFEs still contain formally singular terms
at the origin. These are managed by application of L'Hôpital's
rule. At ℐ^+ we do not require continuum boundary
conditions, and so it is permissible simply to shift the finite
differencing stencils to the left. To minimize reflections from the
outer boundary we use truncation error matching,

    (D f)_N = 1/4h(f_N-4-5 f_N-3+10 f_N-2-11 f_N-1+5 f_N)   .

The idea is that by matching the form of the finite differencing error
at the boundary with that of the interior, high frequency errors are
reduced. Formally singular terms at ℐ^+, which appear in
our first order reduction in exactly the same shape as those of
the GBUF model we studied earlier <cit.>, are likewise
treated with L'Hôpital's rule. As is
standard when using the GHG formulation, we use excision to treat the
black hole, should there be one. The strategy is to keep the apparent
horizon, the position where the expansion of outgoing null-geodesics
vanishes, within the computational domain with the boundary itself
remaining outflow, so that no boundary continuum conditions are
needed. Due to the relatively simple dynamics in the strong-field
region in our experiments, this can be done without the full control
system machinery <cit.> used in binary spacetimes, and
without any careful imposition of the outflow
condition <cit.>.  Instead we simply monitor the position
of the apparent horizon, which is very simple in spherical symmetry
(see <cit.> for a textbook discussion), and monitor the
coordinate lightspeeds at the boundary to check that the outflow
condition is satisfied. Typically we take the domain to extend a small
number of points inside the apparent horizon. To compute derivatives
at the excision boundary we populate ghostzones by fourth order
extrapolation. When evolving initial data that collapse to form a
black hole we must use our two finite differencing operators. If
instead we start from data that already contains a black hole we can
do away with the D̃ operator. We have implemented both
possibilities. We study numerical evolution of various different
choices of initial data. They are described in turn in the following
subsections.








 §.§ Gauge Perturbations




As a first test we perform numerical evolutions of gauge perturbations
of Minkowski spacetime. The simplest procedure to input this type of
initial data is to construct our slices with hyperboloidal
layers, as explained in <cit.>, so the nontrivial field
content lies inside the region where the change of coordinates, in
this case from Minkowski in global inertial time and spherical polars,
to compactified hyperboloidal coordinates, is simply the
identity. Concretely, we use the same
expression (<ref>) for the height function
with m_C_+=0 and we take r_m=0.4. In
the 3+1 language we place a Gaussian perturbation in the lapse with
amplitude A and zero shift. In terms of our variables this
corresponds to

    C_+(0,r) = 1 + A e^-R(r)^2/σ_0^2  , 
        C_-(0,r) = -1 - A e^-R(r)^2/σ_0^2  ,  
       δ(0,r) = ln ( 1 + A e^-R(r)^2/σ_0^2 )  ,  
        ϵ(0,r) = 0   , 
       ψ(0,r) = 0   ,  
         f_D(0,r) = _Tf_D(0,r) =  0  .

where we take σ_0^2=0.02.

The form of the perturbation differs from that in <cit.>
only in that the present perturbation is centered at the
origin. Observe that the choice of GHG implies an evolution equation
for the lapse and shift. In our variables this corresponds to a
non-zero time derivative of the shift which is equivalent to

    _T C_+(0,r) =  -2A e^-2R^2/σ_0^2/σ_0^2( A + e^R^2/σ_0^2)R         , 
       _T C_-(0,r) =    -2A e^-2R^2/σ_0^2/σ_0^2( A + e^R^2/σ_0^2)R        ,  
       _T δ(0,r) = 0  ,    _Tϵ(0,r) = 0   .

With these expressions we construct the FOR fields to perform the
numerical evolutions. By construction, this initial data satisfies the
constraints up to numerical error.

Evolving the data, the dynamics look similar in both the height
function and eikonal setup. The gauge pulses initially propagate
outwards and through ℐ^+ with no sign of reflection. They
leave behind small perturbations in the evolved fields which later
gradually decay. We chose to work with the gauge parameters p=0
in (<ref>), meaning that the
variables C_+ and ϵ should have improved asymptotic decay
relative to the wave equation, whilst C_- and δ should behave
like solutions to the wave equation. Since the scalar field vanishes
identically in this case, the gauge driver field f_D, which we feed
trivial initial data, does too. Even at moderate resolution the method
successfully runs for long-times, at least until t=10^4 in code
units.  The differences in asymptotics of the variables is well
captured by the numerical evolution, and the harmonic constraints
remain small throughout. Evolving for instance with 200 radial
grid-points, at t=10 the harmonic constraints are of
order O(10^-8). In Figure <ref> we
present snapshots of the evolved fields for an evolution with eikonal
Jacobians with compactification parameter n=1.5 at this
resolution. These evolutions exhibit both pointwise and norm
convergence. To demonstrate this, in
Figure <ref> we plot snapshots
of rescaled differences at three resolutions for the same setup as in
the snapshots of Figure <ref>. To avoid
overpopulating the figure we have taken the sum of absolute value of
these differences over all four of the fields plotted in the first
figure. The data are perfectly compatible with second order
convergence as desired. These results are in excellent qualitative
agreement with those of <cit.>, with a different
formulation of GR.
  
Interestingly, not all compactification parameters n behave in the
same manner numerically at the moderate resolution we employ. Despite
the continuum freedom to choose it at the analytical level, only for
values n≤ 1.5 we get the appropriate norm-convergence with
increasing resolution. This can be understood from the empirical
observation that for bigger values the reduction fields get sharp
features at the grid-point at ℐ^+, therefore affecting the
precise numerical cancellation that needs to happen in order for this
scheme to work. We believe that by adjusting the compactification we
will be able to achieve convergent results for the entire range
of n, but since this fact is seen in all of our subsequent numerical
evolutions, in practice we work in the range n∈ [1.25,1.5] here.

Pure gauge wave evolutions on the Schwarzschild background behave
similarly, but since physical dynamics, which we discuss thoroughly
below, necessarily excite gauge waves we do not present them in detail
here.








 §.§ Constraint Violating Initial Data




We now move on to our hardest set of numerical tests, which comprise
of constraint violating initial data with non-vanishing scalar
field. It is important to consider constraint violating data because
in general numerical error violates the constraints in any
free-evolution setup, and so we must be confident that at least
reasonably small finite errors of this type will not cause a
catastrophic failure of the method. Since these initial data excite
also gauge pulses, all aspects of the solution space including gauge,
constraint violations and physical are probed. In all of the following
tests we see similar results for height function and eikonal Jacobian,
so we present a selection of representative plots from each.




  
*Minkowski perturbations: In the following we performed
tests both with and without layers, with r_m = 0.4 and r_m=0. Both
behave similarly, so we present results only for the case r_m=0. We
begin by perturbing the Minkowski metric, setting Gaussian data on all
variables according to

    C_+(0,r) = 1 + C_0 e^-R(r)^2  , 
         C_-(0,r) = -1 - C_0 e^-R(r)^2  , 
       δ(0,r) = δ_0 e^-R(r)^2  , 
        ϵ(0,r) = (δ_0 - ln(1 + C_0)) e^-R(r)^2  , 
       ψ(0,r) = ψ_0 e^-R(r)^2  , 
       f_D(0,r) = 0  ,

with vanishing time derivatives. The initial data for the rescaled FOR
variables are then calculated accordingly. Therefore in this first
test the reduction constraints should remain satisfied at the
continuum level.  When evolving with the eikonal Jacobians we choose
data appropriate for the Minkowski spacetime in global inertial
coordinates, namely

    U^-= e^-δ(0,r)( 2 + C_0e^-R(r)^2)  ,

with the variable U^+ taken from the sourced eikonal equation
itself (<ref>).

The coefficients of the Gaussian in the initial data for C_± are
chosen to give the correct parity at the origin. Similarly, the
coefficient of the Gaussian in the data for ϵ is chosen from
the regularity
condition (<ref>) at the origin
at T = 0. Unsurprisingly, if we take data that violate parity the
code is observed to crash near the origin in finite time.

As a first test we kept the initial perturbations small enough to
avoid complete gravitational collapse. For simplicity, we
took C_0 = δ_0 = ψ_0 = 10^-3. This data do not satisfy
the GHG, Hamiltonian or Momentum constraints. The magnitude of the GHG
constraint violation in the initial data is 10^-4, comparable with
the scalar field itself. We performed numerical evolutions with both
the height-function and eikonal Jacobians, with several different
values of the compactification parameter n within the range
specified at the end of the previous subsection. In all cases the time
evolution is comparable to the gauge wave discussed above in
section <ref>. The initial constraint
violation initially grows and eventually decays, so that by t=10 it
has reduced by 1 order of magnitude. Similar comments apply to
the L^2 norm of the constraints. These data are clearly physically
wrong from the outset since the Bondi mass initially vanishes despite
the fact that the scalar field is non-trivial, and even takes negative
values during the evolution. We observe this behavior both in the
original (<ref>) and
regularized (<ref>) expressions for the Bondi
mass, with the former not even varying monotonically due to constraint
violations.

Another important difference between these experiments and the pure
gauge waves is that the gauge driver variable f_D now actually
varies. Recall that the purpose of the gauge driver is to prevent
log-terms, which are known to afflict plain harmonic gauge, from
appearing in the variable C_-. Using the gauge driver
condition (<ref>) we see no evidence that these
log-terms are present. Given the rescaling in the definition of our
variables (<ref>) such log-terms would result
in divergence of the solution, and so should be very obvious. To check
this, we evolved the constraint violating
data (<ref>) without using the gauge
driver, with identical initial data, and find that the growth
in C̃_- is indeed evident, as can be seen in
Figure <ref>. In fact, evolutions without the gauge driver
can only be performed in a staggered grid, as having a grid-point
at ℐ^+ leads to an explosion of the simulation after the
first time step. All of our simulations except this were performed
with a with a grid-point at future null infinity. Earlier work has
been presented with a combination of these two setups. For instance
the evolutions of <cit.> were performed on a staggered
grid, whereas in <cit.> the grid setup was identical to that
in the majority of our simulations.
 

 

  

 


Next we switched the gauge driver back on. We performed four numerical
evolutions, doubling resolution each time, so that we could examine
multiple curves. Pointwise convergence of the fields at early times
looks very much like that presented in
Figure <ref>. In
Figure <ref> we plot the norm self-convergence
rate from these long experiments. For the test itself we apply the
standard technique, injecting the higher resolution data on to the
coarse resolution grid, taking the norm of the differences between the
very high (V) and high (H) resolutions, the high and medium (M)
resolutions, and finally medium and low (L) resolutions. We then
plot

    Q_1   =log_2(||M-L||/||H-M||) ,  
           Q_2=log_2(||H-M||/||V-H||) ,

as a function of time, from which it can be seen that the simulations
converge at second order as we increase resolution, as expected with
our discretization. Concretely, given a field Z(t,r) and its two
associated first order reduction fields Z^+ , Z^-, see the
discussion around (<ref>), the
continuum limit of the norm we use is,

    ∫[ r^2Z^2
      + ( R'R^2/χ^2)
      ( 2R'-1/2R´χ^2(Z^+)^2
      + 1/2R'(Z^-)^2 )  ] dr .

 




No attempt has been made to tune to the threshold of collapse, but
increasing the amplitude of the data to ψ_0 = 0.355 leads to
apparent horizon formation. As explained above, we also implemented
excision, so that once an apparent horizon is formed, the black hole
interior can be taken out of the evolved region. As a visualization
technique, we put all the evolved variables to zero inside the excised
region so that we can identify where the apparent horizon was found.
Examples of this are shown in
Figure <ref>. Performing longer evolutions we
see that the data appear to settle down, although slow dynamics
continue. At time t=50 substantial constraint violation remains, and
pointwise is in fact substantially larger than the scalar field and
its reduction variables.




  
*Schwarzschild perturbations: Next we take initial data of
a similar type to (<ref>), but now
built as a perturbation on top of the Schwarzschild spacetime. To do
so we adjust several details in our evolution setup, focusing
henceforth exclusively on the case r_m=0.

In order to do excision as previously described, we need to start from
horizon penetrating coordinates. We take the Schwarzschild solution in
Kerr-Schild form, which we recall was given above in our variables in
Eq. (<ref>).

Since this is an exact solution of the EFEs, it is desirable to have a
choice of gauge for which (<ref>) is static in local
coordinates, so that dynamics in the numerical evolution come from the
departure from it. To achieve this we adjust the gauge
sources (<ref>), in which, as
mentioned above (<ref>), the only modification is

    F^σ_SS = F^σ + C_+ - 1/R .

It is easy to see that that expressions (<ref>) are an exact
solution of the rEFEs with the previous choice of gauge. This fact is
reflected numerically in the sense that numerical evolutions with this
exact initial data remain unchanged up to numerical error for long
times, up to t∼ 10^3M, even at very modest resolution. Due to the
way gauge conditions are constructed, this has not yet been achieved
with the approach followed in <cit.>. We observe that the new
term does not alter the asymptotics of the evolved fields, and that we
still need the gauge driver f_D in order to regularize the C_-
field in the presence of a scalar field.

For our next numerical test we put constraint violating Gaussian
perturbations on top of the solution (<ref>), where we take
the BH mass M=1 as defining our units here.  Analogous to the
previous constraint-violating case, the data we start our simulations
with is

    C_+(0,r) = 1-2/R/1-2/R
        + C_p0 e^-(R-3)^2  , 
        C_-(0,r) = -1 + C_m0 e^-(R-3)^2  , 
       δ(0,r) = δ_0 e^-(R-3)^2  , 
          ϵ(0,r) = ϵ_0 e^-(R-3)^2  , 
       ψ(0,r) = ψ_0 e^-(R-3)^2  , 
           f_D(0,r) = 0  , 
        U^-= e^-δ(0,r)
        ( 1 - C_-(0,r) )

with the first-order-reduction fields computed assuming vanishing time
derivatives. Note that the Gaussians are now centered at R=3 and
that in this case there is no relation between ϵ and the
amplitude of the other fields, or between C_+ and C_-, since
we do not have a regular center in this case. We have performed
successful evolutions of this data using again both the
height-function and eikonal Jacobians with several values of n. In
broad terms these data develop in a manner similar to the previous
setups, in so much that part of the the initial pulses still propagate
out to infinity with O(1) speeds. Of course in this case part of the
field content also accretes on to the black hole. At infinity the
scalar field moreover clearly exhibits the direct signal, ringing and
tail phases. These data are particularly important as our first test
in which several of the evolved fields have non-trivial `mass-terms'
at infinity. We observe no particular difficulty in their numerical
treatment, finding both pointwise and norm convergence just as
convincingly as in the previous cases with a regular center. As an
example of this, in
Figure <ref> we show the
norm self-convergence rates with this setup.

To check the effect of the logarithmic mass term with the height
function change of coordinates (<ref>), we
performed also tests with that term omitted. We find that the signal
tends to accumulate at large R without ever leaving the domain,
leading for the evolution to crash in finite time, compatible with our
understanding from above in section <ref> that
without the mass-term included these slices have the `wrong' global
structure, terminating instead at spatial infinity, as depicted in
Figure <ref>. It is clear that the inclusion of the
`mass-terms' is a fundamental ingredient in the method.


 

 
Despite the success in this suite of configurations, we do expect that
if we were to take the initial constraint violations, of whatever
type, sufficiently large then we could cause our numerical method to
fail. We have not attempted to do so, however, since, first, the same
statement would be true even in standard Cauchy evolutions and second,
the task of the hyperboloidal region is to cope with a combination of
stationary features and outgoing waves, with the metric variables
decaying out to ℐ^+. If we were to face a situation in
application in which large errors in the wavezone induced a failure of
the method, either more resolution is needed, or else the wavezone,
which only loosely defined, ought to be taken to `start' further out
and therefore the parameters for the hyperboloidal layer adjusted.




 §.§ Constraint Satisfying Initial Data




The rEFEs are a set of wave equations for all the metric variables, so
up to this point we have successfully tested our regularization
techniques, allowing us to numerically extract the wave signal
at ℐ^+, the ultimate goal of this project. However, in
order to have a physical spacetime that simultaneously solves the EFEs
we need to satisfy all the constraints, namely, GHG, Hamiltonian,
Momentum and FOR constraints, as explained in
section <ref>. Therefore, as a final test
of our scheme, we move on to evolve constraint satisfying initial data
representing perturbations of the Minkowski and Schwarzschild
spacetimes.








  
*Minkowski perturbations: We begin by constructing initial
data for a spacetime that can be thought as a perturbation of the
Minkowski spacetime. We first take the height function approach for
constructing the initial hyperboloidal slice, with m_C_+≡ 0,
as demanded by our method explained in
section <ref>. We
choose ψ(0,r) = ψ_0 e^-R^2 and ψ_n=0,
with ψ_0 = 10^-3 in order to avoid complete gravitational
collapse. With these choices we numerically generate the
solution γ^(1).

A simple way to generate constraint satisfying initial data for
hyperboloidal slices built with the eikonal approach is to
choose U^+ and U^- initial data so that the
Jacobians (<ref>) match the height function
ones (<ref>). Observe that the choice C_+≡ 1
automatically implies that the U^+ constraint,
Eq. (<ref>), is satisfied.

With these details taken care of, we evolve the initial data with both
the eikonal and height-function equations of motion. The basic
dynamics qualitatively resemble those of the constraint violating
case, with regular fields both at the origin and ℐ^+. For
this reason we do not present snapshots in space. Instead, in the
first panel of Figure <ref> we plot the scalar field
waveform at ℐ^+, where we clearly see the field decaying
at late times.

One of the most stringent test on the physics in the present case is
the evaluation of the Bondi mass, Eq. (<ref>) both
for the initial data and the time development. As previously
mentioned, this should be a non-negative and monotonically decreasing
function of time. Both properties can be seen from the middle panel of
Figure <ref>, from which we see that we initially
start with a positive constant value until the time radiation leaves
the domain through ℐ^+. The left and middle panels indeed
indicate that the spacetime asymptotes to the Minkowski spacetime as
the scalar field leaves the numerical domain.

Within our setup the only physical radiation comes from the scalar
field, so our methods can only be claimed successful if this signal is
well-captured numerically. As for the constraint violating initial
data in section <ref>, we performed norm
self-convergence test for the constraint satisfying data, obtaining
similar results to those shown in
Figure <ref>. Focusing instead on the radiation
field, our proxy for gravitational waves, in the third panel of
Figure <ref>, we plot the absolute value of the
rescaled differences of the scalar field at the grid-point
at ℐ^+ as a function of time.  The overlapping of the
three curves in this case shows that the errors of this radiation
signal decrease at the expected rate with increasing resolution, thus
implying that in the limit of infinite resolution we tend to the real
physical solution.

We proceeded to modify the initial data for the scalar field
to ψ=ψ_0 e^-R^2/σ_0^2, with ψ_0=0.8
and σ_0=0.6, in order to generate an apparent horizon
dynamically, which we see at time t∼ 0.5 in code units.
Qualitatively, the evolved fields look much like those presented in
Figure <ref>.  Interestingly our evolved fields
appear somewhat more regular near ℐ^+ than those
of <cit.> in the same physical setup. In contrast to the
constraint violating collapse, the Bondi mass remains positive and
monotonically decreasing for all times settling to a non-zero value
for late times. The apparent horizon mass is non-decreasing. After
black hole formation the code continues to run without problems for at
least t∼ 10^3M at moderate resolutions, where M is the Bondi
mass at late times.




  
*Schwarzschild perturbations: In order to generate
constraint satisfying initial data for Schwarzschild spacetime
perturbations we follow again the steps mentioned in
section <ref>. We started by generating
initial data slice with the height function Jacobians. In the present
case we no longer take m_C_+=0. However, in order for the scheme
to work we need to know the constant m_C_+ exactly, and therefore
we take C_+ identical to the Schwarzschild solution previously
mentioned. We have experimented with various different initial data
choices for the scalar field compatible with our present procedure for
the constraints.

In order to evolve using eikonal hyperboloidal slices in the present
case we again generated initial data for U^+ and U^- so that
eikonal Jacobians match initially the height function
ones. Importantly, the matching of the Jacobians gives a unique
solution for the initial data for U^+ and U^-, so the U^+
constraint (eq. (<ref>)) will not be
satisfied for a generic given S. To overcome this issue we rather
take eq. (<ref>) as defining the
function S. This choice of S does vanish asymptotically, so with
this approach the outgoing radial coordinate lightspeed in lowercase
coordinates still goes to unity, which is the desired property when we
use the eikonal Jacobians.

The basic dynamics proceed as expected, with part of the scalar field
accreting on to the black hole, and the rest gradually propagating out
to null infinity. As a specific example, we take a Gaussian profile
for the scalar field centered at R≃2.1 M, with M=1 from the
reference solution which we perturb,
with ψ(0,r) = ψ_0 e^-(R-2.1)^2/σ^2, σ=0.2
and ψ_0=10^-4. The outcome of a long evolution of this data is
shown in Figure <ref>, where we plot the scalar field
value at ℐ^+ as a function of time for the full non-linear
evolutions with the height function Jacobians. From this we see that
we recover the expected behavior from linear scalar perturbations on
top of Schwarzschild, where we see that the spherically-symmetric
quasi-normal mode of ψ is in good agreement with the first part
of the data, while late time evolution decays as t^-2. This is in
qualitative agreement with the earlier free-evolution results
of <cit.> under different gauges on a staggered grid
(see Figure 3 of <cit.> and Figure 8.24 in <cit.>). For
comparison we also performed evolutions in the Cowling approximation,
which corresponds to taking Schwarzschild spacetime as the fixed
background and evolving the scalar field on top. As expected,
decreasing the amplitude ψ_0 in the initial data for the
non-linear evolutions makes the fitting of the frequencies and tail of
linear theory every time more accurate.

We performed long evolutions of this data with the eikonal Jacobians
as well. They also show quasi-normal mode ringing and tail decay, but
fitting for the known frequencies and tail rate is more involved since
the time at ℐ^+ used in our evolutions then needs to be
post-processed to make a fair comparison, in particular to the Bondi
time coordinate. This is also an issue in the approach
of <cit.>. We postpone a detailed comparison both of these two
cases and of the effect of nonlinearities on the QNM frequencies to
future work.

The Bondi mass of these evolutions has qualitatively the same behavior
as the one displayed in the center panel of
Figure <ref>, namely, it is a strictly positive and
monotonically decaying function of time, where most of its decay
happens when the scalar field leaves the domain
through ℐ^+. Importantly, M_B takes a
value ≃1 at late times, which is close the value we started
with for the constant M of the background perturbed by the scalar
field. Finally, pointwise convergence at ℐ^+ as a function
of time looks qualitatively the same as in the third panel of
Figure <ref>.








§ CONCLUSIONS




Continuing our research program towards the inclusion of future null
infinity in the computational domain in full 3d numerical relativity,
here we presented an implementation of spherical GR in GHG that uses
the dual-foliation formalism to get all the way out. The strategy is
to take the evolved variables to be equivalent to those that would be
solved for in the standard Cauchy problem, subjected to a rescaling to
obtain non-trivial O(1) quantities, and then to change to
compactified hyperboloidal coordinates. In this way, the hope is to
extend the computational domain to null infinity in a manner that
leaves the numerical treatment in the strong-field region absolutely
unchanged, and in the future without symmetry. We examined a broad
suite of initial data, including gauge waves, constraint violating and
satisfying configurations. We considered spacetimes with a regular
center and dynamical black holes, in which case we use the excision
method to remove the interior. In all cases our initial data were
posed on hyperboloidal slices. The coordinate transformation was
managed either by the use of a height-function or by solving the
eikonal equation, both with a given radial compactification containing
a parameter 1<n≤2 that controls how fast the transformation is
made. We have examples in which the compactification takes effect
immediately from the origin, and others which employ hyperboloidal
layers, where it takes effect only further out.

To construct constraint satisfying data from regular equations with
scalar field matter, we treated the overall solution to be a
perturbation of the Minkowski or Schwarzschild spacetimes, taking
suitable initial data for the gauge drivers. In both cases, those
corrections turned out to possess a geometrical and physical
interpretation in terms of the Misner-Sharp mass. We believe that this
technique can be generalized, first to drop simplifying assumptions
within spherical symmetry, but also to full 3d, both of which are
kept for future work.

We find convincing evidence for numerical convergence across the
entire suite of spacetimes we considered. Although we did not make any
push for precision here, for physical initial data we did find good
compatibility vis-à-vis expected frequencies and rates, for instance
in QNMs and Price decay.

It is gratifying to see the line of reasoning developed across the
direct precursors to this study work bear fruit in full GR, even if
only in the spherical setting. In brief, in <cit.> it was
observed that to obtain equations of motion regular enough to treat on
compactified hyperboloidal slices, the coordinate light-speed
variable C_+ needs to display decay beyond that expected for
solutions of the wave-equation. It was then argued in <cit.>
that this could be achieved in the GHG formulation, even when the
constraints are violated, by appropriate constraint addition to the
field equations. In plain harmonic gauge it is known that either slow
decay of the stress-energy tensor or the presence of gravitational
waves serve as an obstruction to decay of the metric
components. In <cit.> it was observed that this
shortcoming of the gauge can be overcome by the use of carefully
chosen gauge source functions. In parallel, numerical studies were
performed with model problems with the same asymptotics as in GR in
GHG. These taught us first <cit.> a convenient choice of
reduction variables, second <cit.> the importance of using
truncation error matching at null infinity, and
third <cit.> that the general strategy to suppress
log-terms is indeed viable in practice.

The interplay between the mathematical and numerical works has also
been important in getting to this stage. For instance,
in <cit.> the suggestion was to force improved asymptotics
for all metric components except for those associated with
gravitational waves. In practice however, with the GBUF model problem,
it was found that this approach would make convergence in the numerics
difficult. (For this reason we have worked here with p=0
in (<ref>) for the C_-
variable). Evidently, all of these ingredients played an important
role in treating spherical GR.

Although the results presented are an important milestone in our
research program, open questions remain both in spherical symmetry and
more generally. So far we have done nothing to chase down sharp
conditions in the required asymptotics of our method. In our current
setup we insist, for instance, on hyperboloidal initial data such that
the radiation field from our scalar matter is O(1)
at ℐ^+. But it is known even in the Minkowski spacetime
that `reasonable' Cauchy data for the wave equation can result in
solutions with logarithmically growing radiation
fields <cit.>. In the future we wish
to understand more clearly the class of data that allows for the
inclusion of ℐ^+ in the computational domain, and how that
class sits within the broader choice that allows analogous growth in
the radiation fields. It would be good also to formulate conditions,
ideally necessary and sufficient, on the choice of gauge that would
allow for the inclusion of ℐ^+ in the computational
domain, even within the better class of data. On a more technical
level we wish to achieve successful numerical evolutions with the most
aggressive compactification parameter n=2, to work without the first
order reduction, and to switch to higher order and pseudospectral
methods. Despite these questions, in view of the 3d toy-model results
of <cit.> and those presented here for spherical GR, we
believe that the essential pieces are now in place to achieve, in the
near-term, the goal of 3d numerical evolutions of full GR in GHG on
compactified hyperboloidal slices.




§ ACKNOWLEDGEMENTS



The authors thank Sukanta Bose, Miguel Duarte, Justin Feng, Edgar
Gasperin, Prayush Kumar and Anil Zenginoğlu for helpful
discussions and or comments on the manuscript.

The Mathematica notebooks associated with this work can be found
at <cit.>.

The authors thank FCT for financial support through Project
No. UIDB/00099/2020 and for funding with DOI
10.54499/DL57/2016/CP1384/CT0090, as well as IST-ID through Project
No. 1801P.00970.1.01.01. This work was partially supported by the ICTS
Knowledge Exchange Grants owned by the ICTS Director and Prayush
Kumar, Ashok and Gita Vaish Early Career Faculty Fellowship owned by
Prayush Kumar at the ICTS. Part of the computational work was
performed on the Sonic cluster at ICTS. SG's research was supported by
the Department of Atomic Energy, Government of India, under project
no. RTI4001, Infosys-TIFR Leading Edge Travel Grant, Ref. No.:
TFR/Efund/44/Leading Edge TG (R-2)/8/, and the University Grants
Commission (UGC), India Senior Research Fellowship. Part of this work
was done at the University of the Balearic Islands (UIB), Spain, the
Department of Mathematics at the University of Valencia and the
Astrophysical and Cosmological Relativity department at the Albert
Einstein Institute, aka the Max-Planck Institute for Gravitational
Physics in Potsdam-Golm. SG thanks Sascha Husa, Isabel Cordero
Carrión and Alessandra Buonanno for local hospitality and travel
support at these institutions.









