Simulations of gravitational collapse in null coordinates: I. Formulation and weak-field tests in generalised Bondi gauges
Carsten Gundlach David Hilditch Thomas W. Baumgarte









Mathematical Sciences, University of Southampton,
  Southampton SO17 1BJ, United Kingdom 

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

Department of Physics and Astronomy, Bowdoin College,
  Brunswick, ME 04011, USA 






  We present a code for numerical simulations of the collapse of
  regular initial data to a black hole in null coordinates. We
  restrict to twist-free axisymmetry with scalar field matter. Our
  coordinates are (u,x,y,φ), where the retarded time u
  labels outgoing null cones emerging from a regular central
  worldline, the angles (θ,φ) label the null generators
  of each null cone, and the radial coordinate x labels points along
  these generators. We focus on a class of generalised Bondi radial
  coordinates x with the twin properties that x=0 is the central
  world line and that the numerical domain (u≥0,
  0≤ x≤ x_max) is a subset of the domain of dependence of
  the initial data on (u=0, 0≤ x≤ x_max). In critical
  collapse, an appropriate choice of these coordinates can be made to
  zoom in on the accumulation point of scale echos of the critical
  solution, without the need for explicit mesh refinement. We
  introduce a novel numerical scheme that in effect reduces the
  angular resolution at small radius, such that the time step
  Δ u for an explicit numerical scheme is limited by the radial
  resolution Δ x, rather than Δ x(Δθ)^2. We
  present convergence tests in the weak-field regime, where we have
  exact solutions to the linearised scalar and gravitational-wave
  equations.


Simulations of gravitational collapse in null coordinates:
  I. Formulation and weak-field tests in generalised Bondi gauges
    Thomas W. Baumgarte
    23 April 2024
============================================================================================================================







§ INTRODUCTION






 §.§ Motivation




Formulations of the Einstein equations on null surfaces are attractive
in both mathematical and numerical relativity for several reasons.

As is generally done in the literature, we will here consider null
coordinates where the surfaces of constant coordinate u are null
hypersurfaces, and the lines of constant (u,θ,φ) are
their null generators <cit.>. Such
coordinates are sometimes called “Bondi-like”, not to be confused
with Bondi coordinates, where in addition the radial coordinate x is
chosen to be the area radius.

As we will review, the Einstein equations in Bondi-like coordinates
can be formulated in a way that makes them maximally constrained, in
the following sense. In affine gauge or Bondi gauge, one solves two
evolution equations, representing two polarisations of gravitational
wave (or one in twist-free axisymmetry). On each null slice of
constant u, the metric is completely determined by solving partial
differential equations (PDEs) with derivatives only in the slice.
Moreover these can be solved by explicit integration along the null
cone generators. In double-null coordinates there is one additional
evolution equation for the area radius R.

This is useful in numerical relativity in two ways. First, the
evolution cannot drift away from a consistent state on each time slice
through numerical error, in the sense that the hypersurface equations
are solved on each time slice, not just the initial one, and so only
free data are evolved. 

Secondly, the equations can easily be discretised in a way that is
compatible with causality. This makes it straightforward to evolve on
the domain of dependence of the initial data, or to impose boundary
conditions on timelike inner and outer boundaries, or to extend the
numerical domain to future null infinity.

A third reason for using Bondi-like coordinates in numerical
relativity is that any results obtained in them have immediate
geometric significance, in contrast to, say, harmonic coordinates or
“puncture” coordinates. For all these reasons, null coordinates are
often the method of choice in spherical symmetry (see
Sec. <ref> for references).

Beyond spherical symmetry, there is an obvious problem with null
coordinates: null surfaces generically form caustics. However, we
expect the expansion of outgoing null cones to prevent caustics in
spacetimes that are sufficiently close to being either flat or
spherically symmetric, with the origin of the null cones near the
centre of approximate spherical symmetry. In a companion paper
<cit.> (from now, Paper II) we will consider axisymmetric
spacetimes with an additional reflection symmetry through the
equatorial plane, so that there is a preferred worldline fixed by this
symmetry.

This motivates us to investigate when null coordinates can be
used to simulate non-spherical gravitational collapse, and how best to
do this. We are not, in fact, aware of any use of null
coordinates in the numerical evolution of regular initial data to a
black hole (gravitational collapse) beyond spherical symmetry. Our
ultimate motivation for this is vacuum critical collapse, see
<cit.> for a general review of critical collapse and
<cit.> for what we consider to be the state of
the art, at the time of writing, in vacuum critical collapse.

The present paper is concerned with general considerations, the
derivation of a number of possible radial gauges adapted to critical
collapse, spacetime diagnostics, the presentation of our numerical
methods, and weak-field convergence tests. In Paper II we apply these
methods to axisymmetric scalar field critical collapse. In
another companion paper <cit.>, we consider issues of
hyperbolicity and well-posedness.





 §.§ Previous numerical work in null coordinates





A spacetime coordinate u is called null if the surfaces of constant
u are null or, in terms of the spacetime metric, g^uu=0. We
speak of null coordinates when one of the coordinates is null, and of
double null coordinates when two of them (usually called u and v)
are null. In our terminology, u will always be an outgoing null
coordinate, also called retarded time. For general review papers on
the use of such (single) null coordinates in numerical relativity see
<cit.>.

One natural choice of null surfaces is the set of null cones emerging
from a regular central worldline. The cones are labelled by the
retarded time u and their null generators by (u,θ,φ).
The fourth coordinate, which we generically call x, then
labels points on each generator.

This formulation has been implemented in twist-free axisymmetry, in
Bondi gauge, where x is the area radius R, compactified at future
null infinity, both in vacuum <cit.> and
with perfect fluid matter
<cit.>. Here, the null cones
emanate from a regular centre. This brings about a severe limitation
of the time step to Δ u∼Δ x (Δθ)^2, see
also Appendix <ref>.

Without symmetry restrictions, Bondi gauge in vacuum has been
implemented in <cit.>, and with
perfect fluid matter in <cit.>,
both using stereographic coordinates on the 2-spheres of constant
(u,x). The vacuum case has also been implemented using angular
coordinates on the 2-spheres in
<cit.>. These papers are focused on
gravitational wave extraction, and so their null cones emanate from a
regular timelike world cylinder, on which boundary data must be
given. This also avoids the time step problem at the origin.

A formulation where the radial coordinate x is the affine parameter
λ along the null generators has been implemented in spherical
symmetry in a cosmological setting with fluid matter in
<cit.>, and in an application to spherical scalar
field critical collapse in <cit.>, and to
vacuum in spherical symmetry with initial data on two intersecting
null cones in <cit.>. Affine
gauge in vacuum without symmetries was formulated in
<cit.>, but not implemented in a code.

Affine gauge has also been used in a number of papers in the context
of asymptotically anti-de Sitter spacetimes, on ingoing null
surfaces emanating from the timelike infinity and terminating inside a
black hole apparent horizon <cit.>. Mentioning only the
two applications most relevant for us, in <cit.> this was
done in 4+1 dimensions with two commuting translation symmetries (and
therefore mathematically similar to the twist-free axisymmetric case
in 3+1 dimensions), and in 3+1 dimensions without symmetries in
<cit.>.

As far as we know, no attempt has been made to simulate the collapse
of regular initial data to a black hole on null cones emanating from a
regular centre, except in spherical symmetry. From among the many
successful applications in spherical symmetry, we review here only the
application to the gravitational collapse of a spherically symmetric
massless scalar field. An early study of scalar field collapse in
Bondi coordinates was <cit.>. Essentially the same
algorithm was used in <cit.> for the study of
power-law tails and quasinormal modes in scalar field collapse.
Double-null coordinates were used for the study of spherical scalar
field critical collapse in <cit.>, Bondi
coordinates compactified at future null infinity in
<cit.>, and (as already mentioned) affine
coordinates in <cit.>.

“Type II” critical collapse is characterised by an arbitrarily large
range of spacetime scales, and hence typically requires adaptive mesh
refinement in numerical simulations. In fact, the pioneering paper
<cit.> was made possible only by the first use of
adaptive mesh refinement in numerical relativity. However, with the
benefit of hindsight, the required mesh refinement zooms in on a
single spacetime point. Once that point has been identified, a much
simpler refinement scheme is possible, such as nested boxes in
Cartesian coordinates. A fixed grid in polar-radial coordinates
centred on this point and spaced logarithmically in radius also
provides the required spatial resolution, but at the cost of a time
step everywhere set by the smallest radial grid spacing.

A fixed grid providing the required mesh refinement in space and time
for spherically symmetric critical collapse was implemented by
Garfinkle <cit.> using double-null coordinates. The
numerical domain is (u≥0,x≤ x_0), with x an ingoing null
coordinate. Hence the numerical domain is exactly the domain of
dependence of the initial data. Its outer boundary is the ingoing null
cone x=x_0, which converges to a point. In the evolution of
near-critical initial data, an appropriate choice of x_0 then puts
the apex of the numerical domain near the accumulation point of scale
echos in near-critical solutions. Whenever half the x-grid points
have fallen into the centre, the resolution is doubled by regridding,
with the time step adjusted accordingly. The numerical grid thus
“zooms in” on the (approximate) critical solution as efficiently as
possible, without the considerable complications of standard adaptive
mesh refinement schemes.

Beyond spherical symmetry, ingoing null cones generically develop
caustics, rather than refocus on a central world line. (This is not a
concern for a sufficiently short time in a setup with initial data on
two null cones u=0 and v=0 that intersect in a spacelike
two-sphere, one often used in mathematical relavity for the study
of black-hole spacetimes, see for example
<cit.>).

However, we can rescue the key idea of <cit.> if we
choose x in such a way that x=0 is the regular centre while
x=x_0 is ingoing null (in spherical symmetry only) or future
spacelike (in general). We first implemented this in spherical
symmetry <cit.>, and found that it provides mesh
refinement for critical collapse as efficiently as the algorithm of
<cit.>.

In hindsight this is similar to using an ingoing radial shift with
spacelike time slices to make the outer boundary ingoing null or
future spacelike. This had already been used for spherical critical
collapse on spacelike time slices in <cit.>.





 §.§ Plan of this paper




A fresh approach to critical collapse using null coordinates looks
promising for the following three reasons, already mentioned
above. First, because of their geometric rigidity, the use of null
cones give us any approximate critical solution we find in coordinates
already adapted to discrete self-similarity. Second, we do not expect
problems with constraint violations. Finally, in a suitable
discretisation the outer boundary (assumed future spacelike) can be
treated exactly like the interior points.

A natural stepping stone from spherical symmetry to vacuum collapse is
non-spherical scalar field collapse, which can be examined with an
arbitrary degree of non-sphericity, whereas vacuum collapse is
necessarily very non-spherical. We shall present our results for
scalar field critical collapse in twist-free axisymmetry in Paper II.

The structure of the present paper is as follows. In order to
highlight the general mathematical structure of the Einstein equations
in null coordinates, in Sec. <ref> we review them in
n+2 spacetime dimensions, without symmetry assumptions.

Starting from Sec. <ref>, we restrict to twist-free
axisymmetry in the usual 3+1 dimensions, and add a massless scalar
field as matter. We review standard gauge choices, and propose several
new ones for critical collapse. We also discuss diagnostic quantities
such as the Hawking mass, and how one can hope to identify apparent
and event horizons. 

Sec. <ref> describes our numerical methods, and
in particular a novel method for completely overcoming the time step
problem mentioned above, such that Δ u∼Δ x.

Sec. <ref> describes convergence tests of the full nonlinear
equations in a small data regime where the linearisation of the
equations about Minkowski is a good approximation.

We conclude with a summary and outlook in Sec. <ref>.
 
A number of appendixes give details of
(<ref>) the time step problem, null
coordinates in (<ref>) Minkowski spacetime and
(<ref>) spherical symmetry,
(<ref>) exact solutions of the linearised equations
that we use as testbeds in the small data regime,
(<ref>) the residual gauge freedom in
our coordinate choice, and (<ref>) regularity
conditions for the metric at the origin and on the symmetry axis.





§ NULL COORDINATES ON GENERIC SPACETIMES OF ARBITRARY
  DIMENSION







 §.§ Metric ansatz




Throughout this paper, a,b,… are abstract tensor indices on the
full spacetime (n+2-dimensional in this section, and
2+2-dimensional in the rest of the paper), and ∇_a is the
covariant derivative with respect to the spacetime metric
g_ab. i,j,k=1...n are angular coordinate indices, and
μ,ν,...=1...n+2 are spacetime coordinate indices, with n ≥ 2 
in this section only, and n=2 in the remainder of the paper. 

In n+2 spacetime dimensions, we define null coordinates
(u,x,θ^i), with i=1...n, by demanding that the hypersurfaces
of constant u are null, in the sense that
g^ab∇_au∇_bu=g^uu=0. The vector field

    U^a:=-∇^au

is obviously null, and obeys

    U^a∇_aU_b=U^a∇_bU_a=1 2∇_b(U_aU^a)=0.

Hence U^a is the tangent vector to the affinely parameterised null
geodesics ruling the null surfaces of constant u.  It is easy to verify
that the most general metric obeying g^uu=0 can be written in n+2
form as

    ds^2   =   -2G du dx-H du^2 
          
      +γ̃_kl
    (dθ̃^k+β̂^k du+β̃^k dx)
      (dθ̃^l+β̂^l du+β̃^l dx).

Note there are

    2+n(n+1) 2+2n=(n+3)(n+2) 2-1

metric coefficients
(G,H,γ̃_kl,β̃^k,β̂^k), which is the full
number of independent metric coefficients in n+2 spacetime
dimensions, minus the one coordinate condition g^uu=0.  At this
point we still have two independent n-dimensional angular “shift
vectors” β̃^k and β̂^k.

We now make the gauge transformation from (u,x,θ̃^k) to
(u,x,θ^i) where θ^i(u,x,θ̃^k) is given
implicitly by a solution of the system

    θ̃^k_,x(u,x,θ^i)=-β̃^k[u,x,θ̃^l(u,x,θ^j)]

of n coupled first-order ordinary differential equations
(ODEs) in x for the functions θ̃^k(u,x,θ^i).
In the new coordinates the metric then takes the form

    ds^2   =   -2G du dx-H du^2 
    
          
    +γ_ij(dθ^i+β^i du)
    (dθ^j+β^j du),

where

    γ_ij   =   γ̃_klθ̃^k_,iθ̃^l_,j,
    
    β^i   =   θ^i_,k(θ̃^k_,u+β̂^k),

and θ^i_,k is the matrix inverse θ̃^k_,i.

With the coordinates in the order x^μ:=(u,x,θ^i), the
metric tensor can be written in matrix form as

    g_μν=([ -H+γ_ijβ^iβ^j            -G       γ_jkβ^k;            -G             0             0;       γ_ikβ^k             0          γ_ij;               ]).

Note that the induced metric on the surfaces of constant u, given by
the bottom right sub-matrix of (<ref>), is
degenerate with signature 0++…+, as one would expect. 

The inverse metric is

    g^μν=([     0  -1 G     0;  -1 G H G^2 β^j G;     0 β^i G  γ^ij;       ]),

where γ^ij is the the inverse of γ_ij. We see that
G=0 would be a coordinate singularity where the metric has no
inverse. Hence we assume that G>0. There is no such restriction on
H.

We see from (<ref>) that in our coordinates the
vector field U^a takes the simple form

    U=1 G∂_x,

and so the outgoing null geodesics that rule each surface of constant
u (its generators) are simply the lines of constant
(u,θ^i). (We have already mentioned that such coordinates are
sometimes called “Bondi-like”). With G>0, x is strictly
increasing with the affine parameter of the null cone generators.

For the following discussions, we denote by N^+_u the
outgoing n+1-dimensional null hypersurfaces of constant u, and by
L^+_u,θ^i the outgoing affinely parameterised null
geodesics that rule each N^+_u [lines of constant
  (u,θ^i)]. We denote by S_u,x the n-dimensional
spacelike surfaces of constant u and x, by N^-_u,x the
ingoing null surface that emerges from each S_u,x, and by
L^-_u,x,θ^i the affinely parameterised null geodesics
that rule it. Note that on N^-_u,x and L^-_u,x,θ^i none of the coordinates are constant: the
coordinate values that label them are only starting values.









 §.§ Standard radial gauge choices





We see from (<ref>) that g_xi=0 for i=1… n, and
from (<ref>) that g^uu=0, both by construction. We
have used up n+1 of the possible n+2 coordinate conditions, with
one remaining to be imposed. 

From an n+2 perspective, this final gauge condition should not
single out any spatial coordinate θ^i, and so should involve
only G, H and γ_ij. If we think of this condition as
fixing, for example, the metric coefficient H, we have

    (n+3)(n+2) 2-(n+2)= 1+n(n+1) 2+n

independent metric coefficients: on the left-hand side the number of
algebraically independent metric coefficients g_μν in n+2
spacetime dimensions without symmetry, minus n+2 coordinate
conditions, and on the right-hand side the number of components in
(G,γ_ij,β^k). We are aware of three such conditions in
the literature.





  §.§.§ Bondi
 



One may be able to assume that the surfaces
S_u,x are n-spheres and that their volume increases
monotonically with x. Then Bondi coordinates are defined by the
coordinate condition γ_ij=x^2γ̅_ij, where
γ̅_ij is the unit round metric on S^n in the coordinates
θ^i. x is then called the area radius and is usually, and in
this paper, denoted by r.





  §.§.§ Double-null
 



Double-null coordinates are defined by
g^xx=0, which, from (<ref>), is equivalent to
H=0. x is then a second null coordinate, and is usually, and in
this paper, denoted by v. The N^-_u,v are now surfaces of
constant v. The affinely parameterised L^-_u,v,θ^i
have the tangent vector

    V^a:=-∇^a v,

which in coordinates takes the form

    V=1 G(∂_u-β^i∂_i).

As already mentioned, we are not aware of a numerical application of
double null coordinates beyond spherical symmetry. This may be because
we expect ingoing null cones to develop caustics generically, rather
than converge to a point.





  §.§.§ Affine
 



Finally, one sees from (<ref>) that x
is an affine parameter along the outgoing null geodesics if and only
if the coordinate condition G_,x=0 holds. x is then often called
λ. λ on each outgoing null geodesic is fixed up to an
additive constant by fixing the function G(u,θ^i). The most
common choice is G=1.





 §.§ The hierarchy of Einstein equations




In order to have nontrivial spacetimes with a regular centre even in
the limit of spherical symmetry, we add a massless minimally coupled
scalar field ψ that obeys the wave equation

    ∇^a∇_aψ=0.

The Einstein equations with scalar field matter in any spacetime
  dimension can be written, in trace-reversed form, as

    E_ab:=R_ab-8π∇_aψ∇_bψ=0.

R_ab is the spacetime Ricci tensor, and we use gravitational units
where c=G=1.

We define the ingoing null vector field

    Ξ:=∂_u-H 2G∂_x-β^i∂_i.

Together Ξ^a and U^a span the normal space to each S_u,x. They are normalised so that Ξ^aU_a=-1. U^a, given
in (<ref>), is tangent to the affinely parameterised
generators of N_u^+, while Ξ^a is tangent to the
generators of N_u,x^- where they emerge from S_u,x.

Given only the metric components γ_ij on a surface of
constant u, and boundary values for G, β^i and
β^i_,x on any past boundary x=x_ min(u,θ^i) of
that surface, we can solve the Einstein equation E_xx=0 for G,
E_x^i=0 for β^i, and E_ij=0 for the null derivatives
Ξγ_ij, in this order, by explicit integration in x. In
the terminology of <cit.> these are the “main
equations”. We will call them the “hierarchy equations”. They
contain H only in the combination Ξ. Explicit expressions in
twist-free axisymmetry in 3+1 spacetime dimensions will be given in
Sec. <ref> below. Given the Ξ g_ij, one gauge
condition is required to fix H and so find the “time” derivatives
γ_ij,u. We can then advance γ_ij in u and repeat.

The remaining n+2 Einstein equations E_ux=0 (the “trivial
equation” <cit.>), E_uu=0 and E_u^i=0 (the
“supplementary conditions” <cit.>) contain higher
u-derivatives than the hierarchy equations and are redundant modulo
the n+2 contracted Bianchi identities. The supplementary conditions
also act as constraints on the boundary data imposed at
x=x_min. We do not discuss them here because
x=x_min=0 will always be a regular central world line here
and in Paper II, so no free data can be imposed there. Finally, the trivial
equation is an algebraic consequence of imposing all other equations.

With the shorthand

    B:=H 2G,

we can write (<ref>) as

    ∂_u=Ξ+B∂_x+β^i∂_i.

This suggests that we consider B and β^i as the x and
θ^i components of a “shift vector” representing the
difference between the coordinate time direction ∂_u and the
null vector Ξ.  However, while the future-pointing unit vector
n^a normal to a spacelike hypersurface Σ is unique, the null
vector Ξ^a depends not only on the null hypersurface N_u^+, but also on its foliation by n-surfaces S_u,x.





§ TWIST-FREE AXISYMMETRY IN 3+1 WITH A SCALAR FIELD







 §.§ Metric ansatz and matter field




We now restrict to 3+1 spacetime dimensions in spherical polar null
coordinates (u,x,θ,φ). We assume axisymmetry with Killing
vector K:=∂_φ, meaning that g_μν,φ=0 in
those coordinates. In addition, we assume a reflection symmetry
φ→-φ. This is a consistent truncation, in the sense
that if we impose γ_θφ=0 on the initial data, and
the boundary conditions β^φ=β^φ_,x=0 at
x=x_min to start up the integration, the hierarchy equations
give us β^φ=0 and
γ_θφ,u=0. Geometrically, the twist vector
ϵ^abcdK_a∇_bK_c vanishes, hence the name twist-free
axisymmetry. Physically, this symmetry removes one of the two
gravitational wave degrees of freedom, hence the alternative name
polarized axisymmetry.

We identify a worldline on the symmetry axis world sheet, calling it
the central worldline, or centre for short. (There is a preferred
choice for this if the spacetime has a reflection symmetry z→-z,
or θ→π-θ, but in general the choice is arbitrary.  For
now we stay in the general case.) The null cones of constant u are
assumed to have a regular vertex on the central worldline.

We parameterise the metric under these conditions as

    ds^2   =   -2G du dx-H du^2
          
    +R^2[e^2F(dθ+β du)^2+e^-2Fsin^2θ dφ^2],

where the metric coefficients (G,H,R,F,β) depend on the
coordinates (u,x,θ) only. R is the area radius, in the sense
that γ_ij=R^2sin^2θ, and so the area of S_u,x is 4π R^2.  The central worldline is at R=0.





 §.§ Regularisation of the axis




The field equations can be regularised on the symmetry axis by
reparameterising <cit.>

    β=:sinθ b,     
    F=:sin^2θ f

and replacing the coordinate θ by 

    y:=-cosθ,      -1≤ y≤ 1.

Intuitively, b is the z-component of the shift vector
β∂_θ, and therefore regular on the
symmetry axis. The metric becomes

    ds^2   =   -2G du dx-H du^2 
          +R^2[e^2S fS^-1(dy+S b du)^2
    +e^-2S fS dφ^2],

where we have defined the shorthand 

    S:=1-y^2=sin^2θ.

We also have

    Ξ=∂_u -B∂_x-bS ∂_y.

Even though the metric now has a division by S, the hierarchy
equations for (G,b,H,Ξ R,Ξ f,Ξψ) are regular on the axis in
the sense that they contain neither square roots nor divisions by y
or S, see Eqs. (<ref>-<ref>) below.

We note in passing the identities

    ψ_,θ   =   √(1-y^2)ψ_,y, 
    ψ_,θθ   =   (1-y^2)ψ_,yy-yψ_,y.

From (<ref>) we see that the usual regularity condition
for scalars on the symmetry axis, ψ_,θ=0 at θ=0 and
π, corresponding to y=-1 and y=1, does not impose a condition
on ψ_,y there. Rather, we see from (<ref>)
that ψ_,y=±ψ_,θθ there, which is
unconstrained by regularity.





 §.§ The equation hierarchy




In twist-free axisymmetry, the Einstein equations
have seven algebraically independent components E_μν. It is
convenient to define the two linear combinations

    E_±:=e^-2FSE_yy± e^2FS^-1E_φφ.

 
Geometric free data on an outgoing null cone of constant u consist of f and
ψ as functions of R and y. In double-null gauge, we also need
to specify R as a function of x and y there, considering this as
fixing a gauge freedom in the initial data.

The two Einstein equations E_xx=0 and E_xy=0 do not contain
any u-derivatives. They can, respectively, be written as

    (lnG R_,x)_,x   =   S_G[R,f,ψ] 
    (R^4e^2S fb_,x G)_,x   =   S_b[R,f,ψ,G].

The Einstein equations
E_+=0 and E_-=0 and the scalar wave equation all contain
u-derivatives. They can, respectively, be written as

    (R Ξ R)_,x   =   S_R[R,f,ψ,G,b], 
    (R Ξ f)_,x   =   S_f[R,f,ψ,G,b]-(Ξ R)f_,x,
    (R Ξψ)_,x   =   S_ψ[R,f,ψ,G,b]-(Ξ R)ψ_,x,

where Ξ is the derivative operator defined in (<ref>), and
H only appears as part of Ξ. The right-hand sides S[f,…]
contain the derivatives f_,x, f_,y, f_,xy and f_,yy
(but not f_,xx), and the same derivatives of R, G, b and
ψ, with the exceptions of ψ_,xy and b_,yy.

The remaining algebraically independent Einstein equations are
E_uu, E_ux and E_yy, and contain u-derivatives other
than those already appearing in the hierarchy equations. In
particular, E_uu (only) contains R_,uu; E_uu and E_uy
contain R_,uy, f_,uy and G_,uy; E_uu and E_ux
contain G_,ux; and all three contain G_,u. We do not
investigate here what relevance these have as constraints on the data
at an inner boundary x=x_min.

The full right-hand sides of the hierarchy equations are as follows,
beginning with the principal terms on a separate line (S_G is all
non-principal), and the scalar field stress-energy terms at the end:


    S_G   =   R R_,x(S^2f_,x^2+4πψ_,x^2) ,
    
    
    S_b   =     2 R^2 S f_,xy-R^2 G_,xy/G-2 R
    R_,xy
          
    +f_,x(-4 R^2 S^2 f_,y+8 R^2 y (f S-1)+4 R S R_,y)+R^2
       G_,x G_,y/G^2+2 R G_,y
      R_,x/G+2 R_,x R_,y
    -16 π  R^2 ψ_,xψ _,y ,
    
    
    S_R   =   1/4 R^2 S b_,xy+G S e^-2 f S R_,yy/2 R-1/2 G S^2 e^-2 f
       S f_,yy+1/4 S e^-2 f S G_,yy
          
    +R^4 S b_,x^2 e^2 f S/8 G+b_,x(-R^2 y/2-1/2 R S
       R_,y)+R S b_,y R_,x-2 b R y
       R_,x
           
    +e^-2 f S[f_,y(-4 G S y (f S-1)-G S^2 R_,y/R)+G y (2 f
       S-1) R_,y/R+G S^3 f_,y^2-S G_,y^2/8 G
          
    +G_,y(-1/2 S^2 f_,y+f S
       y-y/2)+1/2 G (2 f (S (5-4 f (S-1))-4)-1)-G S R_,y^2/2
       R^2+2 π  G S ψ _,y^2], 
    
    
    S_f   =   R b_,xy/4+e^-2 f S G_,yy/4 R
          
    -3 b R yf_,x+R^3 b_,x^2 e^2 f S/8 G+b_,x(-1/2
       R S f_,y-f R y)+b_,y(1/2 R S
       f_,x+R_,x/2)
    
          
    -2 b f y R_,x
    +e^-2 f S(-G_,y R_,y/2 R^2-G_,y^2/8 G R
    +2 π  G ψ_,y^2/R),
     
    
    S_ψ   =   G S e^-2 f Sψ _,yy/2 R
          
    +1/2 R S b_,yψ _,x-1/2 R S
       b_,xψ _,y-b R yψ _,x 
    +e^-2 f S(-G S^2 f_,y/R+G y (2 f S-1) /R+S
         G_,y/2 R)ψ _,y.








 §.§ Formulation in terms of γ:=ln (G/R_,x) 
and R-derivatives




Note that S_G given in (<ref>) has a division by R_,x.  If we
assume that R_,x>0 everywhere, we can reparameterise the metric
variable G by the new variable

    g:=G R_,x  ⇒   G=R_,xg.

Then all x-derivatives in the hierarchy equations can be eliminated in
favour of the derivative

    D:=1 R_,x∂_x,

that is, the derivative with respect to R at constant u and
y. Moreover, g is invariant under reparameterising x, and in
particular is simply 1 in flat spacetime. Obviously, D does not
commute with ∂_u and ∂_y, so we have to specify the
order of mixed derivatives. As a convention, we apply D last. If
we further replace g by

    γ:=ln g,

the equations simplify a little further, and γ appears
undifferentiated only in the two combinations
exp±(2Sf-γ). Hence in our final numerical formulation we
treat γ as the primary variable.

The hierarchy equations now take the form 

    Dγ   =   S̅_γ[R,f,ψ] 
     D(R^4e^2S f-γ D b)   =   S̅_b[R,f,ψ,γ], 
     D(R Ξ R)   =   S̅_R[R,f,ψ,γ,b], 
     D(R Ξ f)   =   S̅_f[R,f,ψ,γ,b]-(Ξ
    R) D f,
    
    
     D(R Ξψ)   =   S̅_ψ[R,f,ψ,γ,b]-(Ξ R) Dψ,

where S̅=S/R_,x. With G=gR_,x, third derivatives of R
appear in the Einstein equations when we write them in terms of g or
γ . In the hierarchy equations, only R_,xxy and R_,xyy
appear. The former, which appears only in the equation for b, is
problematic numerically, but it can be eliminated by writing

    S̅_b[R,f,ψ,γ]=- D(R^2 D(R_,y))+S̃_b[R,f,ψ,γ].

The modified source S̃_b no longer contains third derivatives,
and we can explicitly integrate the first term on the
right. 

The source terms are now


    S̅_γ   =   R (S^2  D(f)^2+4 π D(ψ )^2), 
    S̃_b   =   R^2 (-(4  D(f) (S^2 f_,y-2 f S
    y+2y) 
    -2 S  D(f_,y)+ D(γ _,y)
    +16 π D(ψ ) ψ _,y)) 
          
    +2 R (2 S  D(f) R_,y+ D(R_,y)
    +γ_,y)+2 R_,y, 
    S̅_R   =   1/8 e^2 f S-γ R^4 S  D(b)^2 
    -1/4 R (-S (R  D(b_,y)+4
    b_,y) 
    +2  D(b) (SR_,y+R y)+8 b y) 
    +1/8 e^γ -2 f S(32 f^2 S y^2 
          
    -4  D(R_,y) (S^2 f_,y-2 f S y+y) 
    -8R_,y/R(S^2 f_,y-2 f S y+y)
    +8 S^3 f_,y^2-8 f (4S^2 y f_,y+5 y^2-1) 
          
    +2 S γ _,y D(R_,y)
    -S  D(R_,y)^2+2 S D(R_,yy)
    +4 S R_,yy/R-S γ _,y^2+2 S
    (γ _,y^2+γ _,yy) - 4
          
    -4 γ _,y(S^2 f_,y-2 f S y+y)
    +32 S y f_,y-4 S^2 f_,yy
    -4 S R_,y^2/R^2 +16 π Sψ _,y^2
    ), 
    S̅_f   =   1/8e^2 f S-γ R^3  D(b)^2 
    +1/4(R (2 S b_,y D(f)+2  D(b) (-S f_,y-2 f
       y)+ D(b_,y))-4 b y (3 R  D(f)+2 f)+2
       b_,y)
           
    +e^γ -2 f S/8 R^2(-4
       R_,yγ _,y+(2 R γ _,y-4 R_,y)
        D(R_,y)-R  D(R_,y)^2+R (2
        D(R_,yy)+γ _,y^2+2
       γ _,yy+16 πψ _,y^2)), 
    S̅_ψ   =   1/2( D(ψ ) (R (S b_,y-2 b y))-R S  D(b) ψ _,y)+e^γ -2 f S/2 R(ψ _,y(-2 S^2 f_,y+4 f S y+S
        D(R_,y)+S γ _,y-2 y)+S ψ
       _,yy).







 §.§ Regular centre





When the null surfaces are cones emanating from a regular central
world line, null geodesics leaving the central worldline at the same
time must carry the same u, and those leaving at different times are
naturally identified as setting off in the same direction
(y,φ) via parallel transport along the central wordline. The
only remaining gauge freedom is to relabel x by x̃(u,x,y),
and u by ũ(u).

Obviously, all metric coefficients must be single-valued on the
central worldline, that is, at x=0 they must be independent of y
for all u. We stress this by writing G(u,0,y)=G(u, 0), and so
for all other quantities that are single-valued at the centre.

We show in Appendix <ref> that when the centre
R=0 is regular, geodesic, and has coordinate location x=0, in
our gauge we have b(u, 0)=f(u, 0)=0 and

    g(u, 0)   =   U'(u),
    
    H(u, 0)   =   U'(u)^2,

where U(u) is proper time along the central worldline. Choosing u
itself to be proper time, we have g(u, 0)=H(u,
  0)=1. Because g and G are single-valued at the origin,
R_,x at the origin must also be single-valued, and we denote it by
R_,x(u, 0).





 §.§ Standard radial gauge choices




We review here the standard radial gauge choices already outlined for
arbitrary spacetime dimension above.





  §.§.§ Bondi
 



The condition R=x=:r defines Bondi coordinates
(u,r,y,φ). In particular, the definition

    Ξ R:=R_,u-B R_,x-bS R_,y,

reduces in Bondi gauge, where R_,u=R_,y=0, to

    B=B_Bondi:=-Ξ R R_,x

Given ψ and f as functions of (x,y) on a surface of constant
u, such as u=0, and the gauge initial data R(0,x,y)=x, the
hierarchy equations are solved for (γ,b,Ξ R,Ξ f,ΞΨ) by
integration. We then find B from Ξ R and (<ref>).





  §.§.§ Double-null
 



The condition H=0 and hence B=0 defines double-null coordinates
(u,v,y,φ), where x=:v becomes the second null
coordinate. Given ψ, f and R as functions of (x,y) on a
surface of constant u, the hierarchy equations are again solved for
(γ,b,Ξ R,Ξ f,Ξψ) by integration.  Here R(0,x,y) can
be thought of as fixing a gauge freedom. As already mentioned, this
formulation is only almost-maximally constrained because the evolution
equation for R does not relate to a physical degree of freedom but
propagates this initial gauge choice.





  §.§.§ Affine
 



The coordinate condition G_,x=0 defines x=:λ to
be an affine parameter along the null geodesic generators of our time
slices. In this case, the first hierarchy equation
(<ref>,<ref>) becomes

    R_,xx+(S^2f_,x^2+4πψ_,x^2)R=0,

and so does not contain G at all, but instead becomes an ODE in x
(at constant u and y) for R, given f and ψ.  
Ξ R is now used for finding H, rather than for evolving R. By
taking an x-derivative of the hierarchy equation (<ref>) for Ξ R and a
u-derivative of (<ref>), and eliminating R_,uxx
between them, we find an equation of the form

    H_,xx=S_H[R,b,f,ψ,Ξ R,Ξ f,Ξψ],

which can be integrated twice to find H.

We stress that any null gauge in which we solve the Raychaudhuri
equation (<ref>) for G or γ, such as Bondi or double
null gauge, breaks down where R_,x=0. But, as we will see later,
the divergence ρ_+ of the null generators of our N_u^+
is proportional to R_,x, so this will generically happen in strong
gravity. The only alternative is to solve (<ref>) for R, for
example in affine gauge. We will consider this elsewhere.





 §.§ Choices of radial gauge for critical collapse




We now present possible choices of radial gauge adapted to critical
collapse that generalise the ideas of <cit.> beyond
spherical symmetry. In these coordinates, we want x=0 to be the
regular centre R=0 and the outer boundary x=x_max to be
future spacelike or null. We therefore evolve on the domain of
dependence of the initial data, without the need for an explicit outer
boundary condition, and the numerical domain shrinks with time. We can
then hope to control this shrinking in such a way that resolution of
the critical solution is maintained without the need for adaptive mesh
refinement.

The domain of dependence of the data on any u=u_0,
0≤ x≤ x_max is bounded by the ingoing null surface
N^-_u_0,x_max, whose null tangent vector at
S_u_0,x_max is Ξ. Hence at x=x_max, x should
not decrease along Ξ, or

    Ξ x=-H 2G≥0  at  x=x_ max,

and this must hold for all u. An equivalent requirement is that the
surface x=x_max is null or spacelike, that is

    |∇ x|^2=H G^2≤0  at  x=x_ max.

In short, H and therefore B must be non-positive at
x=x_max.  This also means that numerically we can
consistently upwind the advection term B∂_x in the time
evolution equations at the outer boundary, with the one-sided stencil pointing
into the numerical domain.

At the inner boundary x=0, the condition that R=0 remains at x=0
fixes B to be given by (<ref>). This is positive, and so
again we can upwind consistently with the one-sided stencil pointing
into the numerical domain.

For applications to critical collapse we impose a condition at some
0<x_0≤ x_max that makes x=x_0 approximately null:
approximately in the sense that B=0 there in some average
sense, or that B≤ 0 with equality at one or more values of
y. If a spacetime approaches a self-similar critical solution, and
x_0 is chosen appropriately, the past lightcone of the critical
solution will then be at x≃ x_0. This gauge condition should
then also imply B≤ 0 at the outer boundary.

We now present a few possible gauges that obey these two boundary
conditions at x=x_0 and x=0.





  §.§.§ Shifted double null gauge




Consider first the choice

    B=B_sdn:=
    (1-x x_0)B_Bondi(
      0),

for some constant parameter x_0≤ x_max. The Bondi radial
shift was defined in (<ref>). In other words,
B decreases linearly in x from its value in Bondi gauge at
the centre to zero at the some x_0. Then x=0 remains at
R=0, and x=x_0 is an ingoing null surface.

Our convention that u is proper time along the central worldline
implies that Ξ R=-1/2 there, and so (<ref>) gives

    B_sdn=(1-x x_0)1
      2R_,x(u, 0).

We shall call this gauge choice shifted double null (from now on, sdn)
gauge because it simply rescales the ingoing null coordinate v.
More precisely, if we choose x=v at u=0, then at fixed u the new
coordinate x is a linear function of v defined by v=v_0 at x=0
and R=0 at x=0 <cit.>.

Therefore, sdn gauge is a continuous version of the repeated
regridding of the double-null coordinate v in
<cit.>. We have previously implemented it in spherical
symmetry in <cit.>, and found that it works exactly
as well as our earlier implementation, in <cit.>, of the
original Garfinkle regridding algorithm. The resulting numerical
domains are identical for x_0=x_ max, but we found in
<cit.> that choosing x_ max>x_0,
which gives us a spacelike buffer zone, has the advantage of revealing
more of the critical solution spacetime.

On a regular spacetime beyond spherical symmetry, we expect sdn gauge
to fail because ingoing null cones, and in particular x=x_0, do not
reconverge on the central worldline. It is, however, extremely useful
in spherical symmetry. We will now discuss alternatives that work
better beyond spherical symmetry, but reduce to sdn gauge in spherical
symmetry.

We have also not been able to find a stable discretization of the
equations in sdn gauge beyond spherical symmetry, even in the limit of
weak fields. The instability looks like a purely numerical
problem at the origin, but we cannot exclude formation of caustics or
ill-posedness as problems in the continuum. We have not explored this
further.





  §.§.§ Global shifted Bondi gauge




An alternative starting point is to demand that

    R(u,x,y)=s(u) x,

where s(u) is a function to be specified. This simplifies the
hierarchy equations in the same way as Bondi gauge does, and in
particular G=g/s(u). Substituting this into the equation defining
Ξ R, we have

    B=B_gsB:=s'(u)x-Ξ R s(u).

We shall this class of gauges global shifted Bondi (from now on, gsB)
gauge.

To avoid potential numerical instabilities from evolving R as a
dynamical variable, we update it directly with R_,u=s'(u)x, rather
than the generic expression based on Ξ R plus shift terms, and
also use (<ref>) in simplifying other derivatives. This gauge
does not suffer from the same numerical instabilities as sdn gauge.

To fix s(u), we demand that the surface x=x_0 is ingoing
null or future spacelike, in the sense that H(u,x_0,y)≤
0. This gives

    s'(u)=1 x_0min_yΞ R(u,x_0,y).

As we shall see in Paper II, gsB gauge behaves very different from sdn
gauge in strong fields, already in spherical symmetry, and does not
seem to be a good choice for critical collapse.





  §.§.§ Local shifted Bondi gauge




Yet another starting point is the observation that the instability we
observe in our implementation of sdn gauge seems to be connected to
the y-dependence of R, while the choice (<ref>) is too
restricted. Hence we can attempt the more general gauge

    R(u,x,y)=R̅(u,x),

by setting 

    B=B_lsB
    :=R̅_,u-Ξ RR̅_,x,

where R̅_,u(u,x) and hence R̅(u,x) is yet to be specified.
We call this local shifted Bondi (from now on, lsB) gauge. We require

    R̅(u,0)=0

to keep the origin regular at x=0 and 

    R̅_,u(u,x_max)≤min_yΞ R(u,x_max,y)

to make B≤ 0 at the outer boundary. 

We restrict the possible choices of R̅_,u by demanding that in
spherical symmetry we revert to sdn gauge. An obvious possibility is
to start from the Bondi shift (<ref>) with its spherical
part subtracted (“non-spherical Bondi”, from now on nsB),

    B_nsB:=-Ξ
      R-(Ξ R)_l=0R̅_,x,

where the suffix l=0 denotes the spherical part, and add the (purely
spherical) sdn shift:

    B=B_lsB1:=B_nsB+B_sdn.

In this gauge x_0 is a null surface “on average”.  

We can add a further term that makes the shift non-positive everywhere
at x_0, giving

    B_lsB2:=B_lsB1
    +x x_0(min_yΞ R)-(Ξ R)_l=0 R_,x(x_0).

In this gauge, B_,u is discontinuous at such values of u where
the location in y of min_yΞ R(u,x_0,y) changes
discontinuously. In either lsB1 or lsB2, there is no guarantee that
H<0 at the outer boundary to make it future spacelike.

Another possibility is to subtract, at every (u,x), the global
maximum over y of the Bondi radial shift instead of its spherical part
(“non-negative Bondi”)

    B_nnB:=-Ξ R-min_y(Ξ R)R̅_,x,

and again add the sdn shift,

    B_lsB4:=B_nnB+B_sdn.

This version has the property that every surface of constant
x≥ x_0, and so the outer boundary in particular, is null or future
spacelike. It has the disadvantage that B_,x and B_,u are
discontinuous at all values of (u,x) where the location in y of
min_yΞ R(u,x,y) changes discontinuously.

We can also make a transition from lsb gauge, with the spherical part
given by sdn, near the centre, to full sdn at the outer boundary, so
that the outer boundary is null everywhere. In other words, we
consider

    B_lsBtosdn   :=   B_sdn
          +[1-K_01(x-λ x_0 x_0-λ
        x_0)]
    [
    (Ξ R R_,x)_l=0
    -Ξ R R_,x],

where K_01(x) is a sufficiently smooth switching function with
K_01(x≤ 0)=0 and K_01(x≥ 1)=1, and 0<λ≤ 1 is a
parameter.  We then have pure lsb gauge for 0≤ x≤λ x_0,
and a transition to pure sdn over the interval λ x_0≤ x≤
x_0.  For K_01 on the transition range 0≤ x≤ 1 we have
settled on the 9th order polynomial defined by the first four
derivatives vanishing at x=0 and x=1, and the symmetry condition
K_01(1/2)=1/2.





 §.§ Spacetime diagnostics








  §.§.§ Affine parameter
 



In this section we look at a number of ways of extracting geometric
information from our simulations. After fixing the central worldline,
a completely geometric coordinate system is given by
(u,λ,y,φ), where u labels the null cones emanating
from the central worldline, and (y,φ) label the generators of
these null cones. u is fixed to be the proper time along the central
worldline. The same (y,φ) at different u are identified by
parallel transport of the vector ∂_u along the central
worldline. λ is the affine parameter along the generators,
with origin λ=0 at R=0 and normalisation λ≃ R
near the origin.

The tangent vector to the affinely parameterised generators of our
null cones is U given by (<ref>). We have
dx/dλ:=U^a∇_ax=G^-1 and λ=0 on the central
worldline R=x=0. Integration then gives us

    λ(u,x,y)=∫_0^x G(u,x',y) dx'=∫_0^R g(u,x',y) dR'.

for the affine parameter. Recall that in our convention g=1 at the
origin, so we have the required normalisation.





  §.§.§ Redshift
 



Let τ be the proper time measured by a timelike observer at
coordinate location (u,x,y,φ) and with 4-velocity u^a
(normalised to u^au_a=-1). The redshift of photons emitted from the
central world line at (u, 0), measured by this observer, is

    Z:=dτ du=1 u^a∇_au

where we have used our convention that dτ/du=1 along the central
world line.

One natural choice of a family of timelike observers puts them at constant
R, y and φ. The ansatz 

    u^a=Z^-1((∂_u)^a-R_,u R_,x(∂_x)^a)

for the corresponding u^a gives
u^a∇_aR=u^a∇_ay=u^a∇_aφ=0 and u^a∇_a
u=Z^-1 as required. From u_au^a=-1 we then find that Z is given
by

    Z=(-2g(Ξ R+S bR_,y)-R^2e^2S fS
    b^2)^1 2.






  §.§.§ Hawking mass and Hawking compactness




Following <cit.>, let S be a smooth closed
spacelike 2-surface. Let l^a and n^a be a pair of future-pointing
null vectors normal to S, normalised such that
l^an_a=-1. Let n^a be outgoing and l_a be ingoing. This is
unique up to multiplying n^a by a positive function e^Λ on
S and multiplying l_a by e^-Λ. The projection
operator onto the tangent space of S is

    π_ab:=g_ab+l_an_b+n_al_b,

and is unique. We then define the null congruence expansions

    ρ_+:=1 2π^ab∇_an_b,     ρ_-:=1
      2π^ab∇_al_b,

It is easy to see that

    π^ab∇_a(e^Λ n_b)=e^Λ(∇_an^a+n^bl^a∇_a
    l_b)=e^Λ∇_an^a,

where the first equality holds when n^a is null and the second when
l^a is an affinely parameterised geodesic. A similar result then
holds for e^-Λl_a. Without loss of generality we now define
n^a and l^a to be continued off S as affinely
parameterised geodesics. Then the product

    ρ_+ρ_-=1 4(∇_al^a)(∇_bn^b)

is independent of the normalisation e^Λ, and therefore
uniquely determined by the spacetime geometry and S.

From ρ_+ρ_-, the Hawking mass M of S is now defined
by

    M( S)   :=   1 2√(A( S) 4π) C( S), 
    
    C( S)   :=   1+1 2π∫_ Sρ_+ρ_- dS,

where A:=∫_ S dS is the area of S. We have defined
the “Hawking compactness” C as an intermediate quantity that is of
interest in its own right. In particular, a marginally outer-trapped
surface, defined by ρ_+=0 and ρ_-<0 at each point, has
Hawking compactness C=1, and an outer-trapped surface, defined by
ρ_+<0 and ρ_-<0 at each point, has Hawking compactness
C>1, but the converses are not true.

We now restrict attention to spacelike 2-surfaces S_0 that
lie within a single coordinate null cone. We define such surfaces by

    u-u_0=0,      x-x_0(y)=0.

A basis of tangent vectors to S_0 is given by

    t_(y):=∂_y+x_0' ∂_x,     t_(φ):=∂_φ.

The outgoing and ingoing null vectors orthogonal to S_0,
normalised so that n^al_a=-1 and l^a∇_a u=1, are

    n   =   U, 
     l   =   Ξ+e^-2S fS gR_,x R^2(x_0' ∂_y+x_0'^2 2∂_x).

n^a=U^a holds because the outgoing null surface that emanates from
S_0 is simply the part of N^+_u that lies to the
future of S_0. If and only if x_0(y) is constant, we also
have l^a=Ξ^a.

The corresponding null geodesic expansions are

    Rρ_+   =   1 g, 
    
    2Rρ_-   =   K_0+K_1 x_0'+K_2 x_0'^2+K_3 x_0”.

The K_i are evaluated at (u_0,x_0(y),y), but do not contain
derivatives of x_0(y). Here we only give

    K_0=2Ξ R-R(S b)_,y

for later reference.

The induced metric on S_0 is

    h_ij=x^μ_,ix^ν_,jh_μν=x^μ_,ix^ν_,jg_μν=g_ij,

where x^i:=(y,φ) are the coordinates on S_0. Its
determinant is therefore simply R^4, independent of x_0(y), so the
volume element on S_0 is

    dS=R^2 dy dφ.

We therefore have

    √(A( S_0) 4π)=(1 2∫_-1^1
    R^2 dy)^1/2

and 

    C( S_0)=1+1 2∫_-1^1 2ρ_+ρ_-R^2 dy.

Note that the integral 12∫_-1^1 ...dy is equal to the
l=0 component of its integrand, and that 2R^2ρ_+ρ_-=-1 for
light cones in flat spacetime.

After an integration by parts to eliminate the x_0” term,
C becomes

    C( S_0)=1+1 2∫_-1^1 
    L(x_0,x_0') dy,

where

    L   =   L_0+L_1 x_0'+L_2 x_0'^2, 
    
    L_0   :=   K_0 g=2Ξ R -R(S b)_,y g, 
    
    L_1   :=   S(e^-2S f
    R_,x(R g)_,y-R^3b_,x) gR^2, 
    
    L_2   :=   e^-2S fS g_,xR_,x gR.


If we now further restrict to surfaces S_u,x of
constant u and x, (<ref>) simplifies to 

    2Rρ_-=K_0,

and (<ref>) simplifies to

    C(u,x):=C( S_u,x)=1+1 2∫_-1^1 L_0 dy.


Beyond spherical symmetry, M( S_u,x) need not be monotonic
in x or non-negative. However, in the lsB class of gauges, one can
show that it is non-decreasing with x, and non-negative. To prove
this, note that for R=R(u,x), the integral for A is trivial,
giving √(A/4π)=R. We can pull this factor of R into the
integral for C(u,x) to obtain

    M(u,x):=M( S_u,x)=1 4∫_-1^1 R(1 +2ρ_+ρ_-R^2) dy.

To evaluate M_,x( S_u,x), one can pull ∂_x into
the integral and use integration by parts in y to express
it as 

    M(u,x)_,x   =   R_,x 4∫_-1^1{e^-2S fS 4g^2(g_,y-R^2e^2S f D b)^2 
          
    -2R^4ρ_+ρ_-[S^2( D f)^2+4π(
        Dψ)^2] 
          
    +4π e^-2S fSψ_,y^2} dy.

Hence a sufficient condition for DM:=M_,x/R_,x to be
non-negative is that ρ_+ρ_-<0. (This result is a special case of
Eardley's observation <cit.> that the Hawking mass
increases on a foliation of an outgoing null surface by luminosity
distance, as long as ρ_+ρ_-<0 and the matter obeys the dominant
energy condition.)  From D M≥ 0 and M=0 along the central
worldline we then obtain M≥ 0.

In lsB gauge, where 

    M(u,x)=R(u,x) 2C(u,x)

we can therefore either compute M from (<ref>) with C given
by (<ref>), or by integrating (<ref>) as

    M̃(u,x,):=∫_0^xM_,x(u,x') dx',

with

    C̃(u,x):=2M̃(u,x) R(u,x)

then defined from M̃. In the continuum, C=C̃ and M=M̃, but their discretizations are very different, to the extent that
their agreement is a highly non-trivial test of the accuracy of our
discretization.





  §.§.§ No MOTS on coordinate lightcones u=const





A marginally outer-trapped surface (MOTS) in 4-dimensional spacetime
is a smooth closed 2-dimensional spacelike surface in the spacetime,
such that the generators of its outgoing null cone have zero
divergence. We now ask if a MOTS S_0 can exist that lies in a
single null time slice, defined by ρ_+=0 on the surface u=u_0,
x=x_0(y). 

A first problem with this is that in any gauge where R is specified
as free data and evolved, including sdn, gsB and lsB gauge, the
Raychaudhuri equation will necessarily involve division by
R_,x. To see this, write it as

    (lnG R_,x)_,x=R R_,xS̃_G,
        S̃_G:=S^2f_,x^2+4πψ_,x^2,

where S̃_G is now regular even when R_,x=0. It is now easy
to check that if we write this as a first-order ODE in x for G,
g:=G/R_,x or ρ_+=R_,x/(RG), or as a second-order ODE in x
for the affine parameter λ with λ_,x=G, these ODEs
involve a division by R_,x, and so become singular where the null
expansion ρ_+ changes sign. Where R_,x>0, we can absorb it
into D:=d/dR (as we have done), but D is not defined
where R_,x=0. This division by R_,x is purely a gauge problem:
in any gauge where we specify G a priori and solve the Raychaudhuri
equation for R, for example in affine gauge, R_,x and therefore
ρ_+ can change sign without a problem.

A second, purely geometrical, problem is that when we trace the
future-outgoing null rays that emerge from the MOTS backwards in time
we generically do not expect them to converge to a point, but rather
to form caustics, except in spherical symmetry. In particular this
means that this backwards null surface cannot in general be a coordinate null cone with regular vertex.

We now ask if one can find a non-marginally outer-trapped
surface (OTS) on u=u_0, one where ρ_+≤ 0 everywhere. The
geometric non-genericity argument then does not apply, so we expect to
be able to find OTS (but not MOTS) in affine gauge.

But in twist-free vacuum axisymmetry even this is not possible because
of a third, again purely geometrical, problem: the shear along the
symmetry axis is zero, (as we see from S=0 there), so in vacuum the
Raychaudhuri equation on the symmetry axis on the symmetry axis
reduces to g_,x=0 or R_,xx=0, and so ρ_+=1/(gR)=R_,x/G
cannot become zero or change sign. There may of course be OTS and MOTS in
such a spacetime, but they cannot be embedded in an outgoing null cone
with regular vertex.





  §.§.§ Surfaces of maximal Hawking compactness





In spite of the obstacles set out in the last subsection, in spherical
symmetry we and other authors have happily identified MOTS, and used
their Hawking mass as an estimate of the initial black hole mass
while working in double-null, Bondi or sdn gauge. How did this work?

Quite naively, we were led by what one does in spherical symmetry on
Cauchy surfaces, for example in polar-radial coordinates (t,r) (see
Appendix <ref>): these coordinates become
singular on a MOTS, but one simply identifies the
first appearance of a local maximum in r of C(t,r) with
C≥ 0.99, say, as an approximate MOTS, and estimates M≃ r/2
at its location.

Similarly, in null coordinates we called a surface S_u,x an
approximate MOTS when it was a local maximum in x of C(u,x) with
C≥ 0.99, say <cit.>. Unlike
ρ_+(u,x), which is often a decreasing function of x (while
remaining strictly positive), C(u,x) typically does have one or more
local maxima in x. This singled out a specific S_u,x as
our MOTS candidate. The same approach was also taken by the authors of
<cit.>.

How then can we generalise this procedure in spherical symmetry to the
non-spherical case? The obvious difference is that we no longer have a
preferred foliation of our coordinate null cones into 2-surfaces.

Ideally, we would look for surfaces x=x_0(y) in u=u_0 that
maximise C[u_0,x_0(y)] given by (<ref>). The resulting
Euler-Lagrange equation is a quasilinear second-order ODE for
x_0(y). The boundary term usually obtained in deriving the
Euler-Lagrange equation vanishes, and so the variation is
unconstrained, consistent with our earlier observation following
Eq. (<ref>) that a regular scalar need not have
vanishing y-derivative on the symmetry axis.

In Paper II, we shall for simplicity limit ourselves to finding the
local maxima in x of C(u,x) given by (<ref>) for the
compactness of the coordinate 2-spheres S_u,x, as we
did in spherical symmetry. It turns out that with a smaller threshold
value, say C≥ 0.8, our heuristic criterion consistently
distinguishes collapsing and dispersing solutions. Recall, however,
our earlier observation following Eq. (<ref>) that beyond
spherical symmetry C>1 is necessary but not sufficient for a
spacelike 2-surface to be outer-trapped.





  §.§.§ Event horizons and coordinate null cones




We now show that if a spacetime admits an event horizon, this has at
least one generator in common with each of at least two coordinate
null cones.

The intersection H∩Σ of an event horizon H
with a Cauchy surface Σ with topology R^3 (which
excludes external black holes with two spatial infinities) at
sufficiently late time is a 2-sphere (or consists of disconnected
2-spheres). Hence on H∩Σ the coordinate u attains a
global minimum and global maximum. If the intersection is at least
C^2 in our coordinate system, they are also local extrema. At any
such local extremum u_*, H∩Σ and
{u=u_*}∩Σ have the same (2-dimensional, spacelike) tangent
space V. It follows that the unique future outgoing null geodesic
through that point and normal to V is a generator of both H
and u=u_*.

Independently, in axisymmetry there will be (at least) two local
extrema u_* of u on H∩Σ located at the poles
y=± 1.  If the spacetime has an additional y→-y symmetry
(reflection through the equatorial plane), and H∩Σ
(or one of its components) straddles the equatorial plane, there is a
third local extremum on the equator y=0, while the other two are
related by the reflection symmetry.

It seems unlikely that we can identify any horizon generators, even if
they are also coordinate null cone generators.





§ NUMERICAL METHODS









 §.§ Legendre pseudospectral method in y








  §.§.§ Choice of basis functions, and synthesis matrices




In the physical scenarios we are interested in, it is natural to
maintain constant angular resolution, and to expect the solution to be
smooth. We therefore use a pseudospectral method in the angle
θ, or y. In this subsection, we write N for N_y, the
number of grid points in y.

We represent ψ as

    ψ=∑_l=0^N-1ψ_l(u,x)P_l(y),

where P_l is the Legendre polynomial of order l, proportional to
the spherical harmonic Y_l0. We represent any other quantity that
transforms as a scalar under coordinate changes on the 2-spheres
S_u,x, including the metric components R, G and H, in
the same way.

By contrast, b and f are not scalars but components of a vector
and symmetric 2-tensor on S^2, respectively. We show in
Appendix <ref> that if we represent them as

    b   =   ∑_l=1^N b_l(u,x)P_l'(y), 
    
    f   =   ∑_l=2^N+1 f_l(u,x)P_l”(y),

then in the linearisation of the Einstein and wave equations about
spherical symmetry the different l decouple. We choose this spectral
representation for b, f and the scalars also in the nonlinear
case.

In any pseudospectral method, one goes backwards and forwards between a
finite number of coefficients in Fourier space and an equal number of
carefully chosen collocation points in real space, carrying out
differentiation in Fourier space and nonlinear algebraic operations in
real space.

Here we transform between collocation points y_i and spherical
harmonic components l by full matrix multiplication. This is clearly
inefficient for large N, in contrast to the fast Fourier transform
available for a Fourier series or Chebyshev polynomials.

We call the matrices that take variables from Fourier space (with
index l) to real space (with index i) synthesis matrices. For all
variables that transform as a scalar under coordinate changes on the
coordinate 2-spheres, we use the synthesis matrix
S_il^(0):=P_l(y_i), where y_i are the collocation points (to
be determined below) and l is the spectral index. In other words,
each column of S^(0) represents one P_l. For b we use
S_il^(1):=P_l'(y_i), and for f we use
S_il^(2):=P_l”(y_i). Note this means that with i=1,… N,
we can represent l=0,… N-1 for scalars, l=1,… N for b
and l=2,… N+1 for f.





  §.§.§ Differentiation in y and choice of collocation points




Rather than transforming back to Fourier space in order to
differentiate in y, we differentiate directly in real space. Either
method requires full-matrix multiplication. Making differentiation
exact for certain functions of y then fixes the collocation points.

We implement first y-derivatives through the Legendre-Gauss-Lobatto
differentiation matrix <cit.>

    D_ij=P_N-1(y_i) (y_i-y_j)P_N-1(y_j),     i j, 
    
    -N(N-1) 4,     i=j=1, 
    N(N-1) 4,     i=j=N, 
    
    0,    otherwise,

where the grid points y_2,…,y_N-1 are the zeros of
P_N-1'(y) in increasing order, y_1=-1 and y_N=1. For second
derivatives, we use the matrix square of D. This discretisation is
known to be exact for polynomials up to order 2N-3.  We choose
Legendre-Gauss-Lobatto because we want the north and south poles
y=± 1 to be on the grid. In order to also have the equator on the
grid, we then choose N to be odd, typically 2^K+1.






  §.§.§ Analysis matrices




An analysis matrix takes grid functions from physical to Fourier
space. One might compute the analysis matrix A_li^(0) using the
formulas for Legendre-Gauss-Lobatto quadrature and the fact that
∫ P_mP_n dy=(2n+1)δ_mn/2. However, with A^(0) defined
this way, the product A^(0)S^(0) differs from the expected unit
matrix in that its bottom right element is 2+1/(N-1). This is due to
the fact that Legendre-Gauss-Lobatto quadrature is exact for
polynomials in y up to order 2N-3, whereas this bottom right
element requires integration of P_N-1P_N-1, which is a
polynomial of order 2N-2.

Therefore we use as our analysis matrix A_li^(0) the matrix
inverse of S^(0)_il, which differs from the naive analysis matrix
only in its last row. (Having to make this choice could have been
avoided by using Gauss-Legendre quadrature, which is exact for
polynomials up to order 2N-1.) We similarly define A_li^(1)
and A_li^(2) as the matrix inverses of S_il^(1)
and S_il^(2).





  §.§.§ Discrete versions of continuum identities and
  consistent truncation




We define the N× N matrices

    Δ^(s):=(1-Y^2)D^2-2(s+1)YD,

where Y denotes the diagonal matrix

    Y:= diag{y_i}_i=1^N,

and we have defined the diagonal matrices

    Λ^(s):= diag{λ^(s)_l}_l=s^N-1+s,

where

    λ^(s)_l:=-(l+s+1)(l-s),    l≥ s,

are the eigenvalues of the Laplace operator on S^2 for the “spins”
s=0,1,2.

Our synthesis and differentiation matrices obey discrete equivalents
of the continuum identities (<ref>), (<ref>) and (<ref>)
between Legendre polynomials given in
Appendix <ref>, which in this notation can be
written concisely as

    Δ^(s)S^(s)=S^(s)Λ^(s),

for s=0,1,2.
We also have 

    DS^(0)   =   S^(1)I_+, 
    
    DS^(1)   =   S^(2)I_+, 
    ⇒   D^2S^(0)   =   S^(2)I_+^2,

relating the different spins. Here I_+ denotes the matrix with ones
in the super-diagonal. We also define I_- as the matrix with ones in
the sub-diagonal, and I_0 as the matrix with ones in the diagonal,
except for a zero in the bottom right element. These obey
I_+I_-=I_0. In particular, Eq. (<ref>) is a discrete version
of the identity (<ref>), which we need for the linearised
hierarchy equation for f, Eq. (<ref>), to be discretised
exactly in y. 

With the definition

    Δ^(1):=(1-Y^2)D-4Y   ⇒  Δ^(1)D=Δ^(1),

we can derive the identity 

    Δ^(1)S^(2)I_0= S^(1)Λ^(1)I_-.

Except for the presence of the factor I_0 on the left,
(<ref>) is the discrete version of the identity
(<ref>), which is needed for the linearised hierarchy equation
for b, Eq. (<ref>). Without the right factor of I_0,
the left-hand side of (<ref>) would have non-zero entries
for l=N-1 and N-3 in its last column. With N grid points we can
represent f_N+1 but not b_N+1, so one can think of this as an
aliasing error arising from the spectral truncation. To avoid this
error, which would lead to an instability, we need to suppress the
unpartnered f_N+1.

In the linearised field equations on a general spherically symmetric
background, or in a gauge other than Bondi gauge, f_l and b_l
couple not only to each other but also to ψ_l, G_l, R_l and
H_l. Therefore we must suppress f_N+1, f_N and b_N, as all
three have no scalar partners.





  §.§.§ Half-range spectral method




Most papers on axisymmetric gravitational collapse assume an
additional reflection symmetry, z→-z in cylindrical coordinates,
or y→-y in our spherical coordinates. In such situations, we can
save computing time by only representing the even l in Fourier
space, or the values -1≤ y≤ 0 in real space. Let N̅:=(N+1)/2 denote the number of grid points on the half-range
equivalent to N grid points on the full range. We define
reduced analysis and synthesis matrices as

    A̅^(0,2)   :=   Q_+A^(0,2)V_+^T, 
    S̅^(0,2)   :=   XS^(0,2)Q_+^T, 
    A̅^(1)   :=   Q_-A^(1)V_-^T, 
    S̅^(1)   :=   XS^(1)Q_-^T,

where we have defined (with N̅=3⇔ N=5 serving as a
prototype)

    Q_+:   =   ([ 1 0      ;   0 1 0  ;       0 1;   ]), 
    
    Q_-:   =   ([ 0 1 0    ;     0 1 0;         0;   ]), 
    
    X:   =   ([ 1        ;   1      ;     1 0 0;   ]), 
    
    V_+:   =   ([ 1       1;   1   1  ;     1    ;   ]), 
    
    V_-:   =   ([  1          -1;     1    -1   ;        0      ;    ]).

Note that S̅ with N=5 or N̅=3 is a 3×3 matrix,
etc.  These obey

    A̅^(0,2)S̅^(0,2)   =   I, 
    A̅^(1)S̅^(1)   =   I_0.

To understand this, consider again the case N̅=3: we can
represent P_l for l=0,2,4, and P_l^(2) for l=2,4,6, but
P_l^(1) only for l=2,4.  Q_- has row rank N̅-1, and
therefore S^(1) and A^(1) have rank N̅-1.

We also define separate derivative matrices acting on even and odd
functions,

    D_±   :=   XDV_±^T, 
    
    D_±^2   :=   D_∓ D_±.

We have then the same identities as already discussed, again with the
proviso that P^(2)_2N̅=P^(2)_N+1 can be represented but
must be suppressed for consistency.





  §.§.§ Tests at finite numerical precision





We numerically calculate the synthesis, analysis and differentiation
matrices for selected values of N in Mathematica. The collocation
points (zeros of the Legendre polynomials) need to be determined
numerically, and we do this with precision 10^-50. The matrix
calculations are then carried out with the same precision. We then
save the resulting expressions for the collocation points y_i and
the synthesis, analysis and differentiation matrices result into ascii
files in slightly more than double precision, and read this into the
F90 code at runtime. 

As as test of the numerical error of the pseudo-spectral method, we
have explicitly evaluated the following matrices, which all vanish in
the continuum, numerically in the F90 code (in double precision).

    T_0   :=    A^(0)S^(0)-I, 
    
    T_1   :=    A^(1)S^(1)-I
    I_0
      , 
    
    T_2   :=    A^(2)S^(2)-I, 
    
    T_3   :=   Δ^(0)_+S^(0)-S^(0)Λ^(0), 
    
    T_4   :=   Δ^(1)_-S^(1)-S^(1)Λ^(1), 
    
    T_5   :=   Δ^(2)_+S^(2)-S^(2)Λ^(2), 
    
    T_6   :=    D_+S^(0)-S^(1)I_+, 
    
    T_7   :=    D_-S^(1)-S^(2)I_+
    I_0
      , 
    
    T_8   :=    D_+^2S^(0)-S^(2)I_+^2
    I_+
      , 
    
    T_9   :=   Δ_+^(1)S^(2)I_+
    I_0
       -S^(1)Λ^(1).

The matrices A, S, Δ and Λ are understood as either
half-range or full-range. Where there is a case distinction, the upper
case applies to the full-range setup and the lower case to the
half-range setup. In the full-range case D_+=D_-:=D, and similarly
for Δ and Δ.

To avoid duplication of code, in the numerical code we use the
half-range notation throughout, and simply set both D_+ and D_- to
D in the full-range case, and similarly D_+^2 and D_-^2 to
D^2. By contrast, the combinations Δ and Δ
do not explicitly appear in the nonlinear equations and so are not
stored as matrices, but are used here as shorthands for their
definitions.

We also check that each row of D or D_+ adds up to zero, as the
differencing of a constant grid function should give zero. However, no
such test applies to D_-, which acts only on odd functions of
y. Hence we define another test

    T_10:=D_+ 1,

where 1 is the column vector with 1 in every row.

At N=3 or N̅=2, the error in all tests is zero or at machine
precision. At all higher resolutions, the largest error occurs in
T_5. All errors at the two equivalent resolutions N and
N̅=(N+1)/2 are approximately the same, as we would expect.  At
N=5, N̅=3, this error is ≃ 10^-13. At each doubling
of resolution, it increases by a factor of ∼ 100, up to
4· 10^-5 at N=65, N̅=33, and 8· 10^-3 at
N=129, N̅=65.

We believe that the observed error in these tests is essentially
round-off error in double precision computations in the F90 code, and
that it increases so rapidly with N because the synthesis and
analysis matrices become increasingly ill-conditioned with both
resolution N and spin s. 

We note that while the error in T_0 at N=65, evaluated within
Mathematica, is 10^-49 as expected, in T_1 and T_2 it is
already 10^-23. We can avoid this by noting that P'_l is a
linear combination with integer coefficients of P_l-1, P_l-3
and so on. Hence we can write S^(1)=S^(0)T^(01) and
S^(2)=S^(0)T^(02), where T^(01) and T^(02) are lower
diagonal matrices with integer coefficients, and so can be inverted
exactly. T_1 and T_2 then have a much smaller internal error of
10^-44. The maximum difference between A^(1) or A^(2)
obtained in the two ways is still only 10^-24 at N=65, so they
are identical when reduced to double precision in Fortran
code. However, for N=129, Mathematica cannot find A^(1) and
A^(2) by direct matrix inversion, but we can still find them from
A^(0) and the inverses of T^(01) and T^(02).





  §.§.§ High-frequency filtering





As in any pseudo-spectral code, nonlinear terms spuriously excite
high-l modes through aliasing, and without care this eventually
leads to numerical instabilities, even if the linearised code is
stable. For high-frequency filtering of a grid function in real space,
we therefore multiply by F:=SF̂A, where F̂ is a diagonal
matrix where only the entries corresponding to l≤ l_max
diagonal entries are one, and those for l>l_max are zero. (A
smooth transition from one to zero would also be possible, but we have
not tried this).

When the Einstein equations are linearised around any spherically
symmetric solution, ψ_l, f_l, b_l, R_l and γ_l
couple only for the same l. This suggests that for consistency we
truncate at the same l_max for all variables, rather than the
same degree of polynomial in y, using the appropriate analysis and
synthesis matrices. Define F̂(k) to be the
diagonal matrix with 1 in its first 0≤ k≤ N entries, and 0 in
the remaining ones. We then filter the scalars with F̅^(0):=S^(0)F̂(l_max+1)A^(0), the variable b with F̅^(1):=S^(1)F̂(l_max)A^(1), and the variable f with F̅^(2):=S^(2)F̂(l_max-1)A^(2). These matrices have therefore different rank.





 §.§ Finite differencing in x of the hierarchy equations




In this subsection, we write N for N_x. As we have seen, the
hierarchy equations take the form

    ϕ_ int,x=S(y,ϕ,ϕ_,x,ϕ_,y,ϕ_,xy,ϕ_,yy)

in terms of the intermediate quantities

    ϕ_ int:=(γ, R^4e^2S f-γ D
    b,b,RΞ R,RΞ f,RΞψ),

and the basic variables (metric coefficients and matter field)

    ϕ:=(γ,R, b,f,ψ).

These equations can be solved by integration as

    ϕ_ int(u,x,y)=ϕ_ int(u,0,y)+∫_0^xS(u,x',y) dx',

for the intermediate quantities, and hence the constrained variables

    ϕ_ cons:=(γ,b,Ξ R,Ξ f,Ξψ),

in this order. The evolved variables

    ϕ_ evol:=(f,R,ψ)

are specified freely at u=0 and evolved in u using their
Ξ-derivatives.

In double-null gauge no influence propagates from larger to smaller
x. In a general gauge, this happens onlythrough the shift term
B∂_x. This suggests to us that the hierarchy equations should
be solved by an integration scheme in x that does not evaluate the
integrand to the right of the point where the integral is being
approximated. Then for H≤ 0, strictly no information flows to
larger x. In the following, we present a simple second-order
accurate scheme that has these properties.

We use a grid equally spaced in x, x_i=iΔ x, with x_0=0 and
x_N=x_max, so Δ x=x_max/N. In the code, array
storage for fields at the mid-point x_i+1/2 is labelled by the
array index i, that is by the grid point to its left. We do not
store field values at the first grid point x_0=0 and first mid-point
x_1/2=Δ x/2. Therefore arrays representing grid points have
array index ranging from 1 to N, and arrays representing
mid-points from 1 to N-1, corresponding to i=3/2 to
N-1/2. With this convention in mind, we define the shorthands

    Δ_i ϕ   :=   ϕ_i+1-ϕ_i, 
    ϕ̅_i   :=   ϕ_i+1+ϕ_i 2.

Here we write only the index i corresponding to the coordinate x,
but not the index j corresponding to the coordinate y.

For sufficiently smooth functions ϕ we then have

    (ϕ_,x)_i+1/2   =   Δ_i ϕΔ x+O(Δ x^2), 
    ϕ_i+1/2   =   ϕ̅_i+O(Δ x^2).

We then discretize the integrations (<ref>) with the
midpoint rule, that is

    ϕ_ int,i+1=ϕ_ int,i+S_i+1 2Δ x,

using the discretisations (<ref>,<ref>), the
discretizations of ϕ_,y and ϕ_,yy using D and D^2,
and the discretization of ϕ_,xy and  ϕ_,xyy resulting from
combining (<ref>) with D. The resulting ϕ_
    int,i+1 is second-order accurate but depends only on
ϕ_i+1 and ϕ_i, thus respecting causality. By contrast, the
trapezoid rule would require left derivatives at the right-hand
gridpoint to respect causality. 

In the formulation where we use D defined in (<ref>)
rather than ∂_x, we differentiate

    ( Dϕ)_i+1/2≃(ϕ_,x)_i+1/2
      (R_,x)_i+1/2 = Δ_i ϕΔ_i R

and integrate 

    ϕ_ int,i+1=ϕ_ int,i+S̅_i+1 2Δ_i
    R.






 §.§ Shift terms




For the evolved variables (<ref>), the definition of Ξ
gives

    ϕ_ evol,u=Ξϕ_ evol+Bϕ_ evol,x+Sbϕ_ evol,y

The analogy of B with the x-component of a shift vector suggests
that, in contrast to the x-derivatives in the right-hand side of
(<ref>), ϕ_ evol ,x in (<ref>) should be
upwinded: as right derivatives ϕ_,x^+ for B>0, and left
derivatives ϕ_,x^- for B<0. In particular, no numerical
boundary condition is then required at an outer boundary x=x_
  max as long as B≤ 0 there, or at the innner boundary where
B=B_Bondi>0. We do not upwind the shift term in the y
direction, as we take all y-derivatives spectrally. 

The upwind x-derivatives are evaluated with the second-order accurate
three-point formulas on grid points, namely

    (ϕ_,x)^+_i   :=   2ϕ_i+1-3 2ϕ_i-1 2ϕ_i+2Δ
        x, 
    
    (ϕ_,x)^-_i   :=   -2ϕ_i-1-3 2ϕ_i-1 2ϕ_i-2Δ
        x.


Where R_,x(u, 0) is needed, we use R_0,j=0 in
(<ref>) to evaluate the right derivative as

    (R_,x)^+_0,j=2R_1,j-1 2R_2,jΔ x,

and then average over y to obtain R_,x(u, 0). 





 §.§ The time step problem, and a resolution




An important, if somewhat vaguely defined, necessary condition for
numerical stability of any explicit time-evolution scheme is the
Courant-Friedrichs-Lewy (from now, CFL) condition. This says that the
outermost characteristic cone, in our case the light cone, is
contained in the spacetime numerical stencil. We consider this as a
heuristic guide to a stability limit on the time step, without
carrying out an actual discrete stability analysis. Unusually, in
spherical symmetry, in double-null coordinates (u,v), this causality
condition does not impose any restriction on the time step Δ
u. However, and not surprisingly, we found in <cit.> that
even then a limit Δ u∼Δ x is required for
stability. Beyond spherical symmetry, however, the combination of
spherical polar coordinates and null coordinates imposes a severe
restriction, and we need to address this problem in order to make our
code efficient.

Explicit methods for hyperbolic problems using spacelike time slices
and Cartesian coordinates typically have a time step condition Δ
t∼Δ x, where Δ x is the grid spacing of the Cartesian
spatial coordinates. By contrast, polar spatial coordinates on
spacelike slices give rise to a time step condition of Δ t∼Δ rΔθ, which comes from evaluating the CFL condition
in the tangential direction, Δ t∼ rΔθ, near the
centre, that is at r∼Δ r. In Appendix
<ref>, we show that in polar coordinates on
null cones, the stability criterion is an even worse Δ u∼Δ rΔθ^2.

This is a problem not only when we use finite differencing in the
angle θ. When we split the linearised wave equation into
spherical harmonics as in (<ref>) and finite-difference
in r for given l, we find empirically that the time step limit for
evolving ψ_l(u,r) is Δ u∼ l^-2Δ r. There is no
differentiation in y to give rise to a CFL condition in the
tangential direction, but instead the wave equation decomposed into
spherical harmonics now contains an l(l+1)/r^2 potential. It turns
out that this requires the same restriction on the time stesp as if we
had used finite differencing in θ, with Δθ∼
1/l_max.

In the context of the wave equation on Minkowski spacetime, the
l(l+1)/r^2 barrier can be transformed into a 2(l+1)/r barrier for
a first-order reduction of the wave equation, and this was made
stable for Δ t∼Δ r by the introduction of a suitable
summation by parts (SBP) differencing scheme in r <cit.>. We
could not see how to do this in null coordinates, even for the
flat-space scalar wave equation. 

However, the spectral method in y that we use suggests a possible
remedy to the time step restrictions in polar coordinates, both on
null slices and on the usual spacelike slices: we simply filter out
all spatial frequencies above l∼ i, where l is the spherical
harmonic index, and i the grid index in the radial coordinate
x. In other words, at radius x only spherical harmonics up to
l∼ x/(Δ x) are represented. We discuss how this is done in
the code in the next subsection.

An intuitive understanding of why this boundary condition removes the
extra restrictions on the time step due to spherical polar coordinates
and null coordinates is that, with Δθ∼
1/l_max∼ 1/i∼Δ x/x near
the centre, we now have an effective angular resolution such that
xΔθ≃Δ x, so that the effective “grid cells”
have roughly equal sides, giving us the stability benefit of a
Cartesian grid. 

Outside the central region, for i≳ l_ max, the angular
resolution is constant, giving us the physical benefits of a spherical
grid, namely efficient resolution of ingoing and outgoing waves of
finite l and an outer boundary with spherical topology.

A similar filtering approach to overcoming the stricter CFL limit in
spherical polar coordinates has been presented on spacelike times
lices in <cit.>. Here the filtering uses fast Fourier
transforms in θ and φ.

These adjustments allow us to run the axisymmetric code with the time
step

    Δ u = min_i,jmin(C_1Δ u_1, C_2Δ u_2),

where we have defined the local time step criteria

    Δ u_1(u,x,y)   :=   Δ x
    |R_,xΞ R|, 
    Δ u_2(u,x,y)   :=   Δ x|B|^-1.

The parameters C_1 and C_2 are independent of Δ x and
Δ y (or l_max), and are of order one. As B is the
shift in the x-direction, Δ u≤Δ u_2 is the standard
limit on the time step for any explicit finite differencing of an
advection equation. In double-null gauge, B=0, and so this criterion
is empty, but empirically the limit Δ u≤Δ u_1 on the
time step is required even then. The quantity |R_,x/Ξ R| can be
understood in two ways: as |R_,x|/|R_,u| in double-null gauge,
implying that R always changes less per time step than per grid
point, or as the value taken by B in Bondi gauge. Empirically, we
find that this term guarantees stability in double null and Bondi
gauges and their variants. Because of our gradual suppression of high
angular frequencies near the centre this holds independently of
l_max.





 §.§ Treatment of the central region






  §.§.§ High-frequency filtering after each time step




In order to allow for a time step Δ u∼Δ x in the way
just discussed, in the initial data and after each full time step we
apply a filter to the evolved variables R, f and ψ as
discussed in Sec. <ref>, but with a local l_max.
This sets

    ϕ_evol,l=0  for  l>l_max,local(i)

where 

    l_max,local(i):=min(max(2,2i-2),l_max,global),

which equals 2,2,4,6,... for i=1,2,3,4...
This filtering means that at i=1,2, only l=0,1,2 are present. At i=3,
l=0,1,2,3,4 are present, at i=4, l=0...6, and so on.
l_max,global is even to treat the full and half-range
discretisation equally. We use the relevant analysis and synthesis
matrices for the scalars, b and f, respectively.

The filter always removes at least the top two l-modes of f on
the full-range grid, or the top mode on the half-range grid, as these
do not have a counterpart in γ, H, R, ψ and b. This
means that l_max,global≤ N_y-1 on the full y-range and
2(N̅_y-1) on the half-range. (We always choose N_y and
N̅_y to be odd.)





  §.§.§ Expansion of the field equations about the origin





The filtering near the origin that gets round the bad CFL condition de
facto imposes unphysical numerical boundary conditions, for example
ψ_l(u,x_l)=0, at some x_l∼ lΔ x. This is compatible
with second (or higher) order numerical accuracy in Δ x for
large l, but not for the smallest l. For example, in a regular
continuum solution ψ_l(u,x)∼ x^l near the centre, so imposing
ψ_l(u,x_l)=0 imposes an error of O(Δ x^l). If we want the
code to be second-order accurate, this is acceptable for l≥ 2, but
for l=0 and l=1 we need to impose a more accurate, nonzero,
boundary condition.

In the code, we impose boundary conditions at
an inner boundary (including an unphysical one) as integration
constants when we solve the hierarchy equations by integration in x,
for example the constant c_l in
(RΞψ)_l(u,x)=c_l+∫_x_l^x...dx'. We now obtain the values
of those integration constants that correspond to a regular centre (to
a given order of accuracy) by expanding the full hierarchy equations
in powers of x.

We shall see that for second-order accuracy we need nonzero
integration constants only for the spherical harmonic components
b_1,2, (Ξ R)_0,1,2, (Ξ f)_2 and (Ξψ)_0,1,2
(where the suffix denotes the value of l). The general argument
above applies also to γ_0,1,2 but their integration constants
turn out to be zero.

We assume that in a regular axisymmetric solution of the
Einstein-scalar equations, the scalar field ψ(u,x,y) admits a
convergent expansion around the origin x=0 of the form

    ψ   =   ∑_l=0^∞∑_k=l^∞ψ_l(k)(u) x^k P_l(y) 
       =   ∑_k=0^∞∑_l=0^k …,

where the re-ordering obviously requires convergence. The first suffix
in ψ_l(k) denotes the spherical harmonic and the second one the
power of x.  Appendix <ref> shows that this
[together with analyticity of the ψ_l(k)(u)] corresponds to
analyticity in suitable Cartesian coordinates (t,ξ,η,z).

Expanding any hierarchy equation F=0 as
(<ref>), and truncating this expansion as

    F≃∑_k=0^k_max∑_l=0^k F_l(k)(u) x^k P_l(y)

the result is a polynomial of finite order k_max in both x
and y. This observation guarantees that when, order by order in x,
we set the coefficients of all nonvanishing powers of y to zero
separately to obtain a system of algebraic equations for the
ψ_l(k), the expansion remains exact in y.

A priori, we expand γ and R in the same way as
ψ. However, we impose γ(u, 0)=γ_0(0)=0 in
order to make u proper time at the origin, and so

    γ= ∑_k=1^∞∑_l=0^k γ_l(k)(u) x^k P_l(y).

We additionally impose R(u, 0)=R_0(0)=0 to locate the origin
R=0 at x=0 and R_1(1)=0 to make R_,x single-valued at the
origin. Moreover, we show in Appendix <ref> that
regularity requires R_l(l)=0 for all l, so that

    R =∑_k=1^∞∑_l=0^k-1R_l(k)(u) x^k P_l(y).

Generalising from the behaviour of analytic solutions to the Einstein
equations linearised about Minkowski spacetime (see
Appendix <ref>), and consistent with the
regularity requirements in Appendix <ref> we
expand

    f   =   ∑_k=2^∞∑_l=2^k f_l(k)(u) x^k P”_l(y), 
    
    b   =   ∑_k=0^∞∑_l=1^k+1 b_l(k)(u) x^k P'_l(y).

Here b_1(0)(u) is free, corresponding to a gauge choice, see
Appendix <ref>.  Finally, from
Ξ=∂_u-B∂ x+... we expect that the expansions of
Ξ-derivatives start at one power of x lower, that is

    Ξ R   =   ∑_k=0^∞∑_l=0^k+1(Ξ
    R)_l(k)(u) x^k P_l(y),
    Ξ f   =   ∑_k=1^∞∑_l=2^k+1(Ξ
    f)_l(k)(u) x^k P”_l(y), 
    Ξψ   =   ∑_k=0^∞∑_l=0^k+1(Ξψ)_l(k)(u) x^k P_l(y).


As a test of consistency, we have explicitly expanded all fields to
O(x^5). This allows us to consistently expand the hierarchy equation
for γ to O(x^2), for b to O(x^5) and for Ξ R, Ξ f
and Ξψ to O(x^3), and the resulting coefficient equations can
be solved for (γ,b,Ξ R,Ξ f,Ξψ) to O(x^3).

In the code, we only need the expansions to O(x), as the error of
O(x^2) corresponds to O(Δ x)^2 for the innermost few grid
points. The nonvanishing terms then involve only spherical harmonics
up to l=2, and are

    γ   =    O(x^2), 
    
    b   =   [b_1(0)+( -4π
    R_0(1)^-1ψ_0(1)ψ_1(1)-R_0(1)^-2R_1(3)
           +2R_0(1)^-3R_0(2)R_1(2)) x]P'_1(y) 
          +[-4R_0(1)^-1(f_2(2)+π
      3ψ_1(1)^2) -R_0(1)^-2R_2(3)
           +23R_0(1)^-3R_1(2)^2] x P'_2(y) + O(x^2), 
    Ξ R   =   -1 2P_0(y) 
          +(-R_0(1)b_1(0)-R_0(1)^-1R_1(2))
     x  P_1(y) 
          +O(x^2), 
    Ξ f   =   -R_0(1)^-1f_2(2) x P_2”(y)+O(x^2), 
    Ξψ   =   [1 2R_0(1)^-1ψ_0(1) 
    +1 2(R_0(1)^-1ψ_0(2)
          +R_0(1)^-2(R_1(2)ψ_1(1)-R_0(2)ψ_0(1))) x
    ]P_0(y)
          -1 2R_0(1)^-1ψ_1(1) P_1(y) 
          -R_0(1)^-1ψ_2(2) x P_2(y)+O(x^2).

Note that P_0(y)=1, P_1(y)=y, P_2(y)=(3y^2-1)/2, P'_1(y)=1,
P'_2(y)=3y, P”_2(y)=3 but we have not substituted these values
for clarity of exposition.

Assuming both b_10=0 (the origin is geodesic) and R_,y=0 (lsB
gauge), as we do here and in Paper II, these equations simplify to 

    γ   =    O(x^2), 
    
    b   =   -4π
    R_0(1)^-1ψ_0(1)ψ_1(1) x P'_1(y) 
          -4R_0(1)^-1(f_2(2)+π
      3ψ_1(1)^2)  x P'_2(y) 
          + O(x^2), 
    Ξ R   =   -1 2P_0(y)+O(x^2), 
    Ξ f   =   -R_0(1)^-1f_2(2) x P_2”(y)+O(x^2), 
    Ξψ   =   [1 2R_0(1)^-1ψ_0(1) 
    +1 2(R_0(1)^-1ψ_0(2)
          -R_0(1)^-2R_0(2)ψ_0(1)) x
    ]P_0(y)
          -1 2R_0(1)^-1ψ_1(1) P_1(y) 
          -R_0(1)^-1ψ_2(2) x P_2(y)+O(x^2).


This means that we need to fit only the following coefficients from
the evolved variables R, f and ψ: R_0(1), R_0(2),
f_2(2), ψ_0(0), ψ_0(1), ψ_0(2), ψ_1(1)
and ψ_2(2). ψ_0(0) must be fitted for consistency but is
not used.

We have implemented both a least-squares fit of, say, ψ_l to
ax^l+bx^l+1 (direct fit), and a least-squares fit of ψ_l/x^l
to a+bx (linear fit), and similarly for f_l and R_l. The direct
fit weights grid points by x^l relative to the linear fit. In each
case we fit to the first n_fit points. We obtain R_0(1)
as the value of the left difference at x=0, as this term will have
to cancel the equivalent transport term, and we then fit to
R-R_0(1)(u)x with the general method. We choose to also fit
R_0(3), f_2(3), ψ_1(2) and ψ_2(3), which are not
used, in order to fit two powers of x to each function. The exception
is that we fit three powers to ψ_0, which then requires n_fit≥ 3.

If we restrict to linear perturbations about Minkowski spacetime in
Bondi gauge, with R=R_0(1)(u) x exactly in the background
solution, we are dropping all other R_l(k) and all products of
expansion coefficients. In an early version of the code, we did that,
and also truncated the expansions at different orders from the above,
resulting in the following expansion:

    γ_0,1,2   =    O(x^2), 
    
    b_1   =    O(x^2), 
    R^4b_2,x G   =   -4R_0(1)^2(f_2(2)x^4+6
        5f_2(3)x^5), 
    
    
    b_2   =   -4R_0(1)^-1(f_2(2)x+3
        5f_2(3)x^2), 
    
    
    (RΞ R)_0   =   -1 2R, 
    
    (RΞ f)_2   =   -f_2(2)x^2-3 10f_2(3)x^3, 
    
    (RΞψ)_0   =   1 2ψ_0(1)x, 
    
    (RΞψ)_1   =   -1 2ψ_1(1)x, 
    
    (RΞψ)_2   =   -ψ_2(2)x^2,

with all other components initialised to zero. The almost-linear
simulations presented here were carried out with this expansion, but
we have checked since that the fully nonlinear expansion makes only a
very small difference to the error we have measured in
convergence tests. In particular, the magnitude and qualitative
behaviour of the error is unchanged.





  §.§.§ Integration of the hierarchy equations




We initialise the integrals for all hierarchy equations up to
i=i_expand≥ 1, and integrate from there. We truncate the
integrand to l_max,local(i) when integrating from x_i-1
to x_i. We then use Ξ R to find the x-shift B. This allows
us to calculate the upwinded x-derivatives of the evolved variables
R, f and ψ, as these depend on the local sign of B.

In the integrals for RΞ f and RΞψ we truncate not the
integrand but the whole integral to l≤ l_max,local(i). We
also set

    (RΞψ)_l=-(R Bψ_,x)_upwind,l,
         l>l_max,local(i).

With (<ref>), this gives (Rψ_,u)_l=0 for
l>l_max,local(i). As R does not depend on y in lsB and
gsB gauges, this then also sets (ψ_,u)_l=0, consistent with the
boundary condition ψ_l=0 that we impose on
l>l_max,local(i). We treat RΞ f similarly.





 §.§ Time evolution




We set initial data for f and ψ on u=0,
0≤ x≤ x_max.  In gauges other than affine gauge we must
also initialise R, and in these gauges we think of this
initialisation as pure gauge. We choose

    R(0,x,y)=x/2,

in analogy with R=(v+u)/2 in the standard double null coordinates
on Minkowski spacetime. 

After solving all hierarchy equations, and applying the shift terms in
Ξ, the resulting “time” derivatives ϕ_ evol,u are
discretized using the second-order Runge-Kutta method, with all
hierarchy equations and gauge conditions evaluated at each Runge-Kutta
sub-step, so that our time update can be characterised as the “method
of lines”.

In double-null or sdn gauges, R(u,x,y) is genuinely evolved, but in
Bondi gauge R=x, and only f and ψ are evolved. In gsB gauge
we have R(u,x,y)=s(u)x. Numerically, we evolve only s(u) (as an
auxiliary variable). In lsB gauge we have R(u,x,y)=R̅(u,x). Numerically, we evolve R(u,x,y) but filter out the l>0
components that are created by numerical error after each full time
step.





§ TESTS IN THE ALMOST-LINEAR REGIME





We test convergence of the full nonlinear code in a regime of small
deviations from Minkowski spacetime (with ψ=0). We evolve
in several of the nonlinear gauge choices we have discussed
above. Specifically, we compare sdn gauge (<ref>), gsB gauge
(<ref>) with s(u) given by (<ref>), lsB2 gauge
(<ref>), and lsBtosdn gauge (<ref>), each
using the (G,x) formulation and the (γ,R) formulation, the
latter with and without the integration by parts (<ref>).
All tests use direct fits near the origin with n_fit=2 and
n_expand=1 and the expansion
(<ref>-<ref>). 





 §.§ Linearised solutions as testbeds




In the small-data regime we can use exact solutions of the linearised field equations as testbeds. Small data here means in
practice that the difference between the solutions of the linearised
and nonlinear field equations can be neglected in comparison with the
numerical error in the nonlinear evolution.

For clarity, in this subsection we denote the exact solutions of the
linearised equations by δψ, δ f, δ b,
δ R. These will then be good approximations to nonlinear but
small ψ, f and b, and a small nonspherical part of R.

In the linearised equations, δψ on the one hand, and
δ f, δ b and δ R on the other, evolve
independently, while different spherical harmonic components l also
decouple from each other.

We write the Minkowski background f=b=ψ=0 in a gauge which agrees
with all of our gauge choices. The spherical background metric
coefficients G, R and H in this gauge are given by
Eqs. (<ref>-<ref>) in Appendix <ref>.

The quantities δ b and δ R for the same physical
solution disagree in different gauges. By contrast δψ is
linearly gauge-invariant, and δ f is linearly gauge-invariant
within the class of lsB and gsB gauges; see also
Appendix <ref>. This means that the
linearised solution in Bondi gauge also gives us δψ in any
other gauge, and δ f in any lsB or gsB gauge. In these gauges,
δ R=0 (in the linearised equations). By contrast, in sdn gauge
δ R develops dynamically even if it is set to zero in the
initial data.







 §.§ Convergence test method




To look for second-order self-convergence of a variable
  ϕ with respect to Δ x, we
assume that the Richardson expansion

    ϕ_Δ x(x)=ϕ_0(x)+ϕ_2(x) Δ x^2+O(Δ x^3)

holds, where we have suppressed the other arguments of
ϕ. ϕ_0(x) is the solution in the continuum limit in x,
and ϕ_2(x) is the second-order error, assumed to be leading.

We now distinguish two cases. If ϕ(u,x,y) obeys a hypersurface
equation (PDE in x and y only), then u is just a parameter for
the purposes of the convergence test. At fixed finite resolution in
y we can then think of the PDE as a (large) system of ODEs in x.

If, on the other hand, ϕ obeys an evolution equation (PDE in u,
x and y), then at fixed finite resolution in y we can consider
it as a (large) system of PDEs in u and x. Our time step criterion
(<ref>) scales Δ u in proportion to Δ x. If the
discretisation in u is also at least second-order accurate, then we
expect that the discretisation errors in both u and x are
proportional to Δ x^2, so we are testing convergence in u and
x together.

In halving Δ x exactly, Δ u is halved approximately, but
not exactly, by the application of the time step condition
(<ref>). To compensate for this, we align the output times at
both resolutions exactly by adjusting the last time step coming up to
the scheduled output time.

From pairs of numerical evolutions, we calculate the self-convergence
error estimate

    E_ϕ,Δ x   :=   ϕ_Δ x-ϕ_Δ x/2Δ x^2-(Δ x/2)^2 Δ x_ref^2 
       =   3 4(Δ x_refΔ x)^2
    [ϕ_Δ x-ϕ_Δ x/2] 
       ≃   Δ x_ref^2 ϕ_2.

Any pair of resolutions could be used to estimate the error, but for
simplicity we use Δ x and Δ x/2, as the coarse grid is
then aligned with the fine grid, so we can evaluate the difference on
the coarse grid without interpolation.
If we know the continuum solution ϕ_0, we can also compute
the alternative error estimate

    E_ϕ,Δ x   :=   (Δ x_refΔ x)^2
    [ϕ_Δ x-ϕ_0] 
       ≃   Δ x_ref^2 ϕ_2.


Note that E_ϕ,Δ x depends on the reference
resolution Δ x_ref, the fixed resolution Δ y (or
l_max), and (u,x,y), but for brevity we do not write these
arguments. If E_ϕ,Δ x calculated for two or more
pairs of resolutions is similar at all Δ x below some
threshold, we have pointwise second-order convergence in Δ x,
and E_ϕ,Δ x itself is approximately equal,
pointwise, to the discretisation error in x, or in u and x, at
the reference resolution Δ x_ref and the fixed
resolution Δ y. The error at any other (smaller) Δ x can
be estimated by scaling with (Δ x/Δ
x_ref)^2.

In this paper we do not yet carry out systematic convergence testing
in y, as we expect different spherical harmonics to decouple in
almost-linear evolutions. Hence the error in ϕ_l(u,x) becomes
negligible once l_max>l, and otherwise ϕ_l(u,x) cannot
be represented at all. However, for completeness, looking ahead to
Paper II, we discuss testing convergence in y already here.

On a smooth nonlinear solution, a spectral method should converge
exponentially, but we shall see in Paper II that our code converges
only to second order in Δ y. Hence, to look for second-order
convergence with respect to Δ y∝ 1/l_max∝
1/N_y, we assume that the Richardson expansion

    ϕ_l_max(y)=ϕ_0(y)+ϕ_2(y) l_max^-2 
    +O(l_max^-3)

holds, where we have again suppressed the other arguments of ϕ
and are keeping the resolution in them fixed. We can then consider the
PDE as a large system of ODEs in y. ϕ_0(y) is the solution
in the continuum limit in y (but at fixed finite resolution in x
and u), and ϕ_2(y) is the second-order error, assumed to be
leading. From pairs of numerical evolutions, we then calculate the
quantity

    E_ϕ,l_max1(x,y)   :=   ϕ_l_max1-ϕ_l_max2 l_max1^-2-l_max2^-2 l_maxref^-2
       ≃    l_maxref^-2 ϕ_2.

We can compare ϕ_l(u,x) at different l_max, and so do
not need to align y-grids at different resolutions.

For code checks, it is often useful to plot single Fourier components
E_ϕ_l,Δ x(x,u) of the error against x, and
animate these plots with time u. In
Figs. <ref>-<ref>, for data
dominated by a single spherical harmonic, we show such single-l
errors against x, at a representative moment of time u. In
Fig. <ref>, for data containing all
spherical harmonics, we instead take the root-mean-square (from now,
rms) norm of E_ϕ,Δ x(x,u,y) over
0≤ x≤ x_max and -1≤ y≤ 1, and plot this against
u. We also evaluate the maximum norm of
E_ϕ_l,Δ x(x,u) over x, and the maximum norm of
E_ϕ,Δ x(x,u,y) over x and y, but we do not
present plots here.





 §.§ Single-l tests









  §.§.§ Generalised d'Alembert exact solutions




We have tested two kinds of exact solutions. The first are the
generalised d'Alembert solutions for a single spherical harmonic
derived in Appendix <ref> and given in
(<ref>) for the scalar field δψ_l, and in
(<ref>-<ref>) for the coupled metric perturbations
(δ b_l,δ f_l), both in terms of an arbitrary function
χ of one variable.  These were derived, and used as testbeds,
previously in <cit.>.

We choose the free function χ to be a Gaussian with centre at
0.8 and width 0.2. As the amplitude for χ we choose
10^-11-l, which results in a maximum value ∼ 10^-11 for
ψ and f in the initial data. The numerical domain is
0≤ x≤ 3 with x_0=2 and 0≤ u≤ 1.9. We set N_y=5, and
N_x=64...8192, increasing by fSactors of two. For these tests, and in
contrast to critical collapse applications, the fact that the grid
shrinks to a point is irrelevant, and this is why we stop at u=1.9,
when the grid has shrunk to 0.05 of its original size, and the
waves have essentially left the grid.

At high resolution, our numerical evaluation of the d'Alembert exact
solution at the innermost grid points suffers from large round-off
error, resulting from the division by powers of r up to
r^l+1. To get round this, at small r we implement a truncated
power-series expansion of the d'Alembert solution.

In initial data in lsB and gsB gauge, R is set to the spherical
background value given in Eq. (<ref>), that is
R_0(0,x)=x/2. In sdn gauge, we add to this a Gaussian δ R, so
that we do not have to wait for δ R to develop dynamically.

To start with a summary of the numerical results that follow, sdn and
lsBtosdn gauge are unstable, our flavour of gsb is unstable in the
(G,x) formulation but stable in the (γ,R) formulation (with
and without integration by parts), and the lsB2 flavour of lsB gauge
is stable in all three formulations. We then find pointwise
second-order convergence (both self-convergence and convergence
against the exact solution) for all l-components of all variables.





  §.§.§ Results in sdn gauge
 



In sdn gauge, for l=3 d'Alembert data and l=3 d'Alembert plus
Gaussian-in-R data, we have only gone to 1024 grid points to see
that in all three formulations the error (against the exact solution)
has a constant in x, oscillating in u, part that does not decrease
with resolution. In other words, in sdn gauge there is an instability
at the origin.





  §.§.§ Results in lsB tosdn gauge




In lsBtosdn gauge, f is unstable at the centre and not
converging from the start for l=2,3. (ψ still converges,
presumably because it sees essentially flat spacetime). It is
sufficient to go to 1024 grid points to see this. We have only tried
the (γ,R) formulation with integration by parts.

A mild instability at the centre is still present when the blending is
between 0.3 x_0 and x_0, with pure lsB gauge for x≤
0.3 x_0. This could be because we are then not projecting R
down to l=0 even at those x. We also see something like a linear
gauge shock at x∼ 1.45, in the middle of the blending zone.





  §.§.§ Results in gsB gauge




In the (G,x) formulation, l=2 is unstable in our flavour of gsb
gauge, with again an oscillating instability at the centre giving rise
to an error constant in x that increases with better resolution. We
needed to go to 4096 points to see this. In the (γ,R)
formulation with and without integration by parts, it is stable and
perfectly second order convergent to 8192 points.

By contrast, l=3 in gsb gauge is stable up to 8192 points
already in the (G,x) formulation.





  §.§.§ Results in lsB2 gauge




We have tested lsB2 gauge in the (G,x) formulation only for l=3 and
only up to 1024 gridpoints. The error at 1024 grid points is visually
identical with that in the (γ,R) formulation with integration
by parts. We have tested lsB2 gauge in the (γ,R) formulation
with integration by parts in more detail. All error plots in the
following are in this gauge and formulation.

In particular, we have tested the d'Alembert solution for l=0,1
(ψ only), l=2,3,4 (ψ, f and b) and l=5 (f and b
only, not suppressing this highest frequency for once).
We find second-order pointwise convergence of ψ_0, ψ_1
and, for l≥ 2, of ψ_l, f_l and b_f, from N_x=64 to
N_x=8092 radial grid points.

For l=2, the above statement needs to be qualified. The error in
ψ_2 is not smooth near the origin, but has a blip at the second
grid point (rather than a specific value of x). This is demonstrated
in Fig. <ref>. However, we get pointwise
second-order convergence at fixed x for all resolutions Δ x<x/2.

Fig. <ref> shows the scaled error in f_2. Here,
it is not clear if we have reached convergence even at N_x=4096, and
it appears again that the error is not smooth at the origin, although
in a different manner to that in ψ_2.

Fig. <ref> shows the scaled error in b_2. Near
the origin, the scaled error is E_b_2,Δ x(x)∼
(Δ x)^2/x. At the innermost grid point this evaluates to ∼Δ x. Hence b_2 converges pointwise to second order at all
points, including near the origin, but converges only to first order
in the maximum norm, which at high resolution is dominated by that
innermost grid point. In other norms, such as the root-mean-square
norms, the convergence would be at an intermediate order.

For l>2, we suspect that the error is still technically unsmooth
near the origin, but this is hard to see as the errors, like the
variables themselves, are suppressed at the origin by powers of x^l
for f, ψ and the other scalars, and x^l-1 for b. Hence
the second-order pointwise convergence looks fine. It is possible that
our methods can be improved to make the error better behaved at the
origin, but we leave this to future work.

We note finally that, where we have the exact solution, the
self-convergence estimate of the error is pointwise approximately
equal to the true error (difference from the exact solution).













To test convergence at l>5, we change from the d'Alembert solution (which
becomes increasingly complicated) to simple Gaussian data
for a single l. As we need larger N_y to represent larger l, we
go to the half-range, with N̅_y=17. We now set initial data
directly for f_l and ψ_l, namely a Gaussian again with centre
at x=0.8, width 0.2, and amplitude 10^-11. The numerical
domain is again 0≤ x≤ 3 with x_0=2 and 0≤ u≤ 1.9, with
N̅_y=17, and N_x=64...8192, increasing by factors of two. At
l=2, we find the same as with the d'Alembert l=2 data. At l=4,8,16,32
we find apparent perfect second-order pointwise convergence (as for
the l=4 d'Alembert data). 





 §.§ Plane-wave tests









  §.§.§ Plane-wave exact solutions




We show in Appendix <ref> that any
solution of the scalar wave equation gives rise to a solution of the
linearised Einstein equations, without the need for a spherical
harmonic decomposition. In particular, the scalar wave equation admits
plane-wave exact solutions, and in Appendix <ref>
we give these in Eq. (<ref>) for δψ and in
(<ref>-<ref>) for the related (δ
b_l,δ f_l). Of course, these are just the translation into our
coordinates of the linear limit of the well-known plane-symmetric
gravitational waves. This is our second exact solution testbed.

A scalar field plane wave moving in the negative z-direction on flat
spacetime takes the form

    ψ=χ(t+z)=χ(u+r+z)=χ(u+r(1+y)).

Its contours at constant z are parabolas in the (ρ,z) plane
(where ρ^2:=r^2+z^2), as
one expects of the intersection of a null plane with a null cone. When
this intersection reaches the vertex of the cone it closes up and
disappears.

We choose χ to be a Gaussian with centre 1.0 and
width 0.1, and with amplitude 10^-11 for ψ, which means that
the maximum of ψ is exactly 10^-11, and amplitude
10^-13 for f, which means that the maximum of f is
approximately 4· 10^-11. The numerical domain is again 0≤
x≤ 3 with x_0=2, now with 0≤ u≤ 1.5. This means that the
plane wave crosses the origin and then disappears out of the numerical
domain during the evolution.

The (single) plane wave solutions contain both even and odd spherical
harmonics. In order to test plane waves with our half-range
formulation, we add a copy of the same wave moving in the opposite
direction to make ψ and f even functions of y (at constant
u and x), and b an odd function.





  §.§.§ Results in lsB2 gauge




We evolve with all combinations of N_y = 17,33,65,129, N̅_y=9,17,33,65
(for the double plane wave only), and N_x = 64...1024. 

This double plane wave solution contains only even spherical
harmonics, and so can be run on the full or half range in y. We have
verified that the errors in the double plane wave test are identical
in the equivalent resolution pairs.  We remove the top two frequencies
in f and top frequency in b in the exact solution before the
comparison with the numerical solution, as the numerical solution is
similarly truncated.

Beginning with the single plane wave, we evaluate the maximum and rms
norms of the error against u. The norms are taken over all grid
values of x and y, with y=±-1 weighted half in the rms norm in
order to make it the same on the full and half-grid. The rms error in
ψ, f and b for the single plane wave at the different
resolutions is shown in Fig. <ref>.

The figure consists of nine plots laid out in a square. The three
columns from left to right show the variables ψ, f and
b. Focus initially on the middle column, showing f. Each plot shows
four differences of angular resolution and five differences of radial
resolution. Different radial resolutions are distinguished by line
type, and different angular resolution by line colour.

In the top plots, the rms error is not scaled. We see that the
error decreases quickly with angular resolution N_y, but is almost
independent of radial resolution N_x. This indicates that for
N_y=17,33 and N_x=64...1024, the total error budget is
dominated by angular discretisation error.

The middle plots zoom in on the smallest errors, which arise from
the highest angular resolutions. They demonstrate that at N_y=65,129
the error does decrease with increasing radial resolution N_x. In
the lowest plot these errors are scaled with (N_x/128)^2. The
scaled curves with N_y=129 (omitting N_x=64) lie on top of each
other, indicating second-order convergence with N_x. This indicates
that at N_y=129, for N_x=128...1024 the error budget is dominated
by the radial discretisation error.

Comparing now the three variables ψ, f and b, we note that as
b is computed from f at each time step, it shows a numerical error
already in the initial data, whereas ψ and f are evolved from
analytic initial data, and so their error is zero at u=0. For
N_y=129 only, the initial error in b is neglible, and the
subsequent error in b is similar to the error in f or ψ.

We see perfect second-order convergence with N_x in ψ (but not
f and b) already at N_y=65. Moreover, the errors in ψ with
N_y=65 and 129 (purple and red curves) are indistinguishable in
our plots.  We conclude that a plane scalar wave requires less angular
resolution than a plane gravitational wave of the same shape χ
(where χ is a solution of the scalar wave equation). This may be
because the exact solution for f and b involves derivatives of χ.

A closer look at the middle and bottom plots shows a transition from
an error dominated by y-discretisation at early times and small
N_y and by x-discretisation at late times and large
N_y. Plotting single-l components of the error against u and x
further shows that the early discretisation error in y arises mostly
at x≳ 1.5 and the late error discretisation error in x
mostly at x≲ 1.5.

The rms errors in the double plane wave test are very similar as
functions of u, but larger by roughly √(2), as there are now
two identical waves of error instead of one in the same numerical
domain. We therefore do not present plots here.









 §.§ Tests of the computation of the Hawking mass




As a different indication of numerical error, we have implemented the
expression for M given by (<ref>) with (<ref>), a
centred finite differencing of this to compute the resulting M_,x,
and, independently, the direct expression (<ref>) for
M_,x. The two expressions for M_,x are affected in different
ways by finite differencing error in x and spectral error in y,
and so their agreement in the continuum limit is a non-trivial test of
the correctness of (<ref>) with (<ref>), (<ref>),
the hierarchy equations, and their discretisations. The cleanest test
is one where the solution is close to Minkowski, so that the hierarchy
equations are approximately linear in ψ, b and f, and the
expressions for γ and M are approximately quadratic.

We have tested single-l Gaussian initial data in f or ψ, and
have varied the amplitude, N_x, N_y, and the numerical method. We
have not evolved these data.

We first take a Gaussian in f only with centre at x=1, width
0.25, x_ max=5, l=2, and amplitude 10^-4. Then
M∼ 9· 10^-8 at the outer boundary, and M_,x has a
maximum of ∼2.5· 10^-7. Our baseline resolution is
N_x=1000 and N_y=5. At this baseline, the difference between the
two versions of M_,x, an estimate of the error in either, has
maximum ∼ 1.1· 10^-10, a relative error of
∼ 4· 10^-4.

Decreasing the amplitude by a factor of 10 reduces M and the error
in M_,x by a factor of 100, as expected. If we decrease the
amplitude much more, the error is dominated by noise. This would be
expected as round-off-error in the cancellation between 1 and
2R^2ρ_+ρ_- in C.

The error with N_y=5,9,17,33,65,129 is the same, and so for N̅_y=3,5,9,17,33,65 on the half-range. Doubling N_x from the
baseline reduces the error by a factor of 4. Hence the error in M in
almost-linear evolutions is dominated by finite differencing in x,
while spectral error in y is negligible.

Initial data in which only ψ is non-vanishing test other parts of
the two expressions for M_,x. Note that b 0 when hierarchy
equations are solved for these data even with f=0. We use a Gaussian
with the same centre and width as for f above. An amplitude of
10^-4 gives M∼ 5.5· 10^-8. M_,x has a maximum of
∼ 1.5· 10^-7, and the error at baseline has a maximum of
∼ 4.5· 10^-11, a relative error of ∼ 3·
10^-4. All other comments for pure f data just above also apply
here.





§ CONCLUSIONS





We have begun an investigation of the use of null coordinates in
numerical relativity, applied to gravitational collapse. In this first
paper we have both made progress and identified new difficulties.





 §.§ Progress




On the purely mathematical side, we have shown that the ingoing null
derivative Ξ normal to the surfaces of constant (u,x) plays a
role similar to the Lie derivative L_n normal to spacelike
time slices. Specifically, and oversimplifying a bit, the Einstein
equations give us the geometric time derivative L_n
  L_n g_ij on the usual spacelike time slices, but Ξ g_ij
on null slices. Writing the hypersurface equations in terms of Ξ
removes any explicit appearance of the radial shift B, just as
writing the 3+1 equations in terms of L_n removes any
explicit appearance of the lapse and shift.

On the numerical side, we have removed a practical obstacle to using
null cones with a regular centre, already identified in
<cit.>, namely that the time step
Δ u∼Δ x(Δθ)^2 is unreasonably small at high
angular resolution. We have been able to replace this by
Δ u∼Δ x, similar to Cartesian coordinates, by reducing
the angular resolution at small radius.

For applications to critical collapse in spherical polar coordinates,
we have found an equivalent of the method of repeated radial
regridding of <cit.> that works beyond spherical
symmetry. Essentially this is done by adding to a standard gauge
choice, such as Bondi gauge, an ingoing radial shift that shrinks the
grid continuously, with the outer boundary becoming future spacelike
<cit.>.

In the present paper, we have demonstrated convergence of these
numerical methods in the evolution of weak data, and in Paper II we
successfully apply them axisymmmetric scalar field critical
collapse.





 §.§ Problems




As we have discussed, one cannot find marginally outer-trapped
surfaces embedded in null cones with a regular centre, and in vacuum
axisymmetry not even any (non-marginally) outer-trapped surfaces. We
propose using closed 2-surfaces of large Hawking compactness (Hawking
mass/area) as an alternative diagnostic of black hole formation. In
Paper II, this allows us to find the threshold of black-hole formation
by bisection, and demonstrate critical scaling of the black hole mass.

It becomes clear in Paper II that in sufficiently non-spherical
spacetimes the divergence of the congruence of the generators of our
null cones becomes negative in some directions even when no black hole
is formed. Bondi coordinates and double-null coordinates, and the
generalised Bondi coordinatees of Paper II, break down when this happens. 





 §.§ Outlook




Two key questions remain open. The first is wheter outgoing null cones
emanating from a regular centre remain regular in strong gravity
further away from spherical symmetry, and in particular in
electromagnetic or vacuum critical collapse. This is a purely
geometric question, independent of gauge or numerical methods. Note
also that the potential problem is non-sphericity, rather than black
hole formation. One reason to be optimistic is that in
<cit.> we have constructed outgoing null cones
emanating from a regular vertex in post-processing of 3+1
near-critical vacuum evolutions, as a way of comparing evolutions in
different coordinate systems, without finding caustics.

A second question is how much useful information about any newly
formed black hole we can find, given that outgoing null coordinates
cannot penetrate into the horizon, or at least not very deeply, and
that we cannot find MOTS on a single coordinate null cone.

As already discussed, going further will probably require a change to
a generalised affine parameter gauge. We leave this to future work.







We would like to thank Luis Lehner, Nigel Bishop, David Garfinkle,
Bernd Brügmann, Tomáš Ledvinka, Anton Khirnov and Jonathan Luk
for helpful discussions, and the Mathematical Research Institute
Oberwolfach for supporting this work through its “Oberwolfach
Research Fellows” scheme. Bastien Roy and Bernardo Porto-Veronese
contributed through M1 internships in 2021 funded by École
Polytechnique. DH was supported in part by FCT (Portugal) Project
No. UIDB/00099/2020. TWB was supported in part by National Science
Foundation (NSF) grants PHY-2010394 and PHY-2341984 to Bowdoin
College.













§ THE TIME STEP PROBLEM IN NULL COORDINATES





We show here that the time step for explicit finite differencing
schemes is even more severely limited in the combination of polar
spatial coordinates with a regular centre and a null “time”
coordinate u than for polar coordinates and the usual (spacelike)
time coordinate t. For simplicity we restrict to axisymmetry in 3+1
dimensions. We focus on the causal geometry, rather than giving a
rigorous argument. Our presentation closely follows
<cit.>.

Consider at first the Minkowski metric in standard spherical polar
coordinates,

    ds^2 = -dt^2+dr^2+r^2(dθ^2 +
    sin^2θ dφ^2).

In an abuse of notation, we now let dx^μ stand for the small
finite coordinate distance from the grid point that is being updated
to any one of the other grid points from which the update is
calculated. These will be integer multiples of the grid spacings
Δ x^μ.  For an explicit numerical method, the update is found
from only a few neigbouring grid points to the past. The set of these
points is called the “stencil” of the point being updated. Causality
requires that the numerical stencil be wider than the physical light
cone, or, in our notation, ds^2>0 for the outermost grid points in
the stencil. Empirically, this is also a necessary criterion for the
stability of any explicit finite difference scheme, often referred to
as the Courant-Friedrichs-Lewy (CFL) condition.

For simplicity we assume that the stencil involves, besides the point
being updated, only grid points on the previous time level. We
parameterise such a stencil by dt=-Δ t, dr=sΔ r,
dθ=qΔθ, with Δ t,Δ r,Δθ>0 by
definition, and s,q=-1,0,1 as a simple example of a 3× 3
point stencil in r and θ. We assume axisymmetry, so nothing
depends on φ.  φ. We also have r=pΔ r, where,
for example, p=1,2,3,… on a centred equally spaced grid. What
matters here is only that s and q are of order unity and so is the
smallest possible value of p, which occurs next to the centre.

The CFL condition ds^2>0 translates into

    Δ t<Δ r√(s^2+p^2q^2Δθ^2).

The choice s=0, q=± 1 of points in the stencil, with the choice
p=1 to locate the stencil next to the centre, gives

    Δ t≲Δ r Δθ,

where we have replaced < by ≲ to allow for more general
stencils and for other O(1) factors in the argument. As is
well-known, this is worse than the time step restriction
Δ t≲min(Δ x^i) in Cartesian coordinates by the
factor Δθ≪ 1. Put differently, the right-hand side of
(<ref>) is quadratic in small quantities.
Halving the grid spacing in all spatial directions halves Δ t
in Cartesian coordinates, but reduces it to a quarter in polar
coordinates.

Consider now the Minkowski metric in the null coordinate form

    ds^2 = -2G du dx-Hdu^2+R^2(dθ^2 +
    sin^2θ dφ^2).

We let du=-Δ u, dx=-sΔ x, dθ=qΔθ, and
R=pΔ x. At the centre R=0, any radial gauge must approach
Bondi gauge to keep it at x=0, and setting H=1 and R=x there without loss
of generality, we also have G=1.  The CFL
condition(<ref>) is now

    Δ u≲Δ x(√(s^2+p^2q^2(Δθ)^2)-s).

With s=1, p=1, q=± 1, and assuming Δθ≪ 1,
the equivalent of (<ref>) is now

    Δ u≲1 2Δ x (Δθ)^2,

worse by a second factor of Δθ≪ 1 than for polar
coordinates on spacelike time slices. Except in spherical
symmetry, this is a problem of any explicit numerical time evolution
scheme on null cones with a regular vertex.

At large R, the sharpest CFL condition arises from stencil points with
q=0 (and p then drops out). We now consider arbitrary G and H.
The CFL condition becomes

    Δ u<-2sG H Δ x.

Recall that G>0.  Choosing s=± 1 with sign opposite to that of
H we obtain

    Δ u<Δ x B,

where B:=H/(2G) as defined previously. This is just the CFL
condition for the x-advection term in
∂_u=Ξ+B∂_x+Sb∂_y.





§ MINKOWSKI SPACETIME





In Minkowski spacetime, we can choose coordinates where f=b=0 and
the remaining metric coefficients R, G and H depend only on
u and x. We denote them by R_0, G_0 and H_0.  With these
assumptions, there are only two nontrivial hierarchy equations, namely

    (lnG_0 R_0,x)_,x   =   0, 
    (R_0Ξ_0R_0)   _,x=   -1 2G_0.

Clearly, these can be integrated in terms of two arbitrary functions of u.

To clarify what these two functions are in a regular spacetime, we
note that in the standard double null coordinates (U,V), the
Minkowski metric is

    ds_0^2=-dU dV+(V-U 2)^2 dΩ^2.

We now change to the most general null coordinates u and x adapted
to the spherical symmetry, defined by

    U   =   U(u), 
    
    V   =   2R_0(u,x)+U(u),

and obtain the metric

    ds^2=-2G_0 du dx-H_0 du^2+R_0^2 dΩ^2,

where

    G_0   =   U'R_0,x, 
    
    H_0   =   U'(U'+2R_0,u),

and hence 

    Ξ_0R_0=-1 2U'.

Therefore the general solution of (<ref>,
<ref>) with a regular centre is (<ref>) and
(<ref>): note this has only one free function
U(u). (<ref>) can be written as

    g_0=U'(u).

In Bondi coordinates, where R_,u=0, we have

    H_0,Bondi=U'(u)^2.


If we choose x=0 to be a geodesic and the coordinate basis vectors
to be parallely transported along it, the metric at x=0 is given by the
Minkowski metric (<ref>) at x=0 for all times. 

In flat spacetime, shifted double null coordinates, shifted global
Bondi coordinates and a natural choice of local shifted Bondi
coordinates are all identical. To find the metric in these coordinates, we
start again from the metric (<ref>) and make the
specific coordinate transformation

    U   =   u-x_0, 
    
    V   =   (u-x_0)(1-x x_0),

where x_0>0 is a parameter,
to obtain (<ref>) with 

    R_0   =   1 2(1-u x_0)x, 
    
    G_0   =   1 2(1-u x_0), 
    
    H_0   =    1-x x_0.

It follows that

    g_0   =   1, 
    
    B_0:=H_0 2G_0   =   x_0-x x_0-u, 
    Ξ_0R_0   =   -12.

We call this the shifted Minkowski (from now on, sM) background gauge.





§ SPHERICAL SYMMETRY





We restrict to spherical symmetry, but bring in a spherical scalar
field as matter, by setting f=b=0 and making R, G, H and
ψ functions of u and x only. The metric is

    ds^2=-2G du dx-H du^2
    +R^2 dΩ^2.

The hierarchy equations become

    Dln g   =   4π
    R( Dψ)^2 
     D(RΞ R)   =   -g2, 
     D(RΞψ)   =   -(Ξ R) Dψ,
 
where

    Ξ=∂_u-H g D.


In spherical symmetry, diagnosing collapse and estimating the horizon
mass is straightforward. The Hawking compactness C and mass M
are given by

    M=R 2 C,      C=1+2Ξ R g,

and, with |∇ R|^2=-2Ξ R/g, the Hawking mass in spherical
symmetry is equal to the well-known Misner-Sharp mass

    M=R 2(1-|∇ R|^2).

(Unlike the Hawking mass in general, the Misner-Sharp mass in
spherical symmetry derives from a conserved stress-energy current
<cit.>.) The mass aspect (<ref>) is

    M_,x=-4π R^2R_,xΞ R g( Dψ)^2,

and the redshift defined in (<ref>) is

    Z=√(-2gΞ R),

We see that, as long as Ξ R remains finite, g→∞ gives
both Z→∞ and C→ 1, so this is the obvious criterion for
black hole formation, both from the compactness reaching one and the
red shift from the centre to infinity diverging.

Assuming that the time step is limited by the Bondi shift
(<ref>), as in Eq. (<ref>), we have

    Δ uΔ x≲|R_,xΞ R|
    = G |g||Ξ R|.

At finite G and Ξ R, this goes to zero again as g→∞, or
equivalently R_,x→ 0.

A commonly used non-null coordinate system in spherical symmetry is
the polar-radial one, defined by

    ds^2=-α^2 dt^2+a^2 dR^2+R^2 dΩ^2,

where R is now a coordinate. It is instructive to relate this to our
null coordinates.

Expressing (t,R) in terms of (u,x), and comparing coefficients of
(<ref>) and (<ref>), we can write
the resulting three equations as

    t_,x   =   aαR_,x, 
    
    g   =   -a^2(R_,u-α at_,u), 
    
    H   =   g(R_,u+α at_,u),

and from this we can derive

    Ξ R   =   -g 2a^2
    =1 2(R_,u-α at_,u), 
    
    C   =   1-1 a^2=1-|∇ R^2|, 
    
    Z   =   g a=α t_,u-aR_,u.

From (<ref>) and (<ref>) we find

    Z du=α dt-a dR,

and so the redshift of the centre with respect to other observers at
constant R is

    Z=α.dt du|_R
    =.dt_0 du|_ R,

where t_0 is proper time along worldlines of constant R. In a
static spacetime only, we can set t_,u=1 without loss of
generality and R_,u=0, and we obtain Z=α.





§ AXISYMMETRIC LINEAR PERTURBATIONS OF MINKOWSKI SPACETIME









 §.§ Linearised field equations in any radial gauge




As a test of our formulation and numerical methods beyond spherical
symmetry, we linearise about flat spacetime. We denote the background
values of all fields by a subscript 0, and so ψ_0=0. We adapt
the background coordinates to spherical symmetry by making G_0,
H_0 and R_0 functions of u and x only. Although G_0=R_0,x
holds in the background, for clarity we write either G_0 or
R_0,x in the linearised equations, as they appear in the
linearisation process.

The linearised hierarchy equations are

    (δ G G_0-δ R_,x
      R_0,x)_,x=0

for δ G, then

    (R_0^4 δ b_,x G_0)_,x   =    
    R_0^2
    (2S δ f_,xy-8yδ f_,x) 
          -R_0^2 G_0δ G_,xy -2R_0δ R_,xy+2R_0,xδ
    R_,y
          +(R_0^2G_0,x G_0^2+2R_0R_0,x G_0)
    δ G_,y,  
    (R_0Ξ_0 δ f)_,x   =   R_0 4xδ b_,xy+R_0,x 2δ b_,y
    -(Ξ_0R_0)δ f_,x
          +1 4R_0δ G_,yy

for δ b and δ f,

    (R_0Ξ_0δ R+(Ξ_0R_0)δ R
    -R_0R_0,xδ B 
    )_,x=δ S_R

for either δ R or δ H (we have not written out δ
S_R as it is long), and 

    (R_0Ξ_0δψ)_,x=
    G_0S 2R_0δψ_,yy-G_0y R_0δψ_,y
    -(Ξ_0R_0)δψ_,x

for δψ. We see that δψ decouples from the metric equations in
any radial gauge, and is independent of the choice of linearised
  radial gauge.

With the background solution G_0=R_0,x, and the boundary
condition δ G=δ R_,x at the origin, which follows
from linearising the gauge condition G=R_,x that u is proper
time at the origin, the solution of
(<ref>) is

    δ G=δ R_,x

everywhere. This is far as we can go without choosing a (linearised)
radial gauge. 





 §.§ Linearised field equations in linearised Bondi gauge




The perturbation equations take their simplest form in Bondi
gauge. This is defined by R=x=:r in the full equations, and hence by
R_0=x, G_0=H_0=1 in the background, and δ R=δ G=0 for
the perturbations. 

However, the choices of background gauge and linear perturbation gauge
are in principle independent.  In particular, we can choose to use sM
gauge in the background, but linearised Bondi gauge δ R=δ
G=0 for the perturbations.  Under a change of background gauge, for
example from Bondi to sM, the linear perturbations δ b, δ
f and δψ change only their argument, as if they were
scalars.

In linearised Bondi gauge, defined by δ R=δ G=0, δ
f and δ b obey the coupled equations

    (R_0^4 b_,x G_0)_,x   =    
    R_0^2
    (2S δ f_,xy-8yδ f_,x), 
    (R_0Ξ_0 δ f)_,x   =   R_0 4δ b_,xy+R_0,x 2δ b_,y
    -(Ξ_0R_0)δ f_,x.

These are just the first lines of (<ref>), (<ref>) above.
The linearised wave equation (<ref>) does not simplify further.
Using G_0=R_0,x, as well as δ G=δ R=0, the linearised
hierarchy equation (<ref>) becomes

    (-R_0 2δ H)_,x   =   R_0^2 4(Sδ b_,xy-2yδ b_,x) 
          +R_0,x[(6S-4)δ f
    -2yR_0δ b 
          +S 2(2R_0δ b_,y+8yδ f_,y
    -Sδ f_,yy)].

Given a solution (δ f,δ b) of (<ref>,<ref>),
(<ref>) can be solved for δ H by
integration, but as δ H does not couple back
into the equations for δ b and δ f, we can ignore it. 

The perturbation equations in linearised Bondi gauge
(<ref>,<ref>) are given in Bondi background
coordinates as Eqs. (<ref>), (<ref>) below.





 §.§ Other linearised radial gauges





Consider now the equations linearised about Minkowski in any linear
gauge obtained by linearising a nonlinear gauge. The different
spherical harmonics l still evolve independently.

If R_,y=0 in the nonlinear gauge, such as in gsB and lsB gauge,
the l 0 components of δ R are absent, and hence, from
(<ref>), also the l 0 components of δ
G. (<ref>,<ref>) then still reduce to
(<ref>,<ref>), and exact solutions derived in
linear Bondi gauge are still relevant, after changing only their
argument.

By contrast, in any gauge where R_,y 0, such as sdn gauge,
δ R and δ G couple back to δ f and δ b.
However, δ f transforms as a scalar under changes of radial gauge,
even nonlinearly, and because f is already a perturbation the change
due to the change of argument is quadratically small. By contrast,
δ b changes already to linear order.  δ R does transform
as a scalar, but R_,x 0 in the background solution, so the
change of argument changes δ R to linear order. In summary, we
can still use the exact solution for δ f as a testbed, but
not δ b.





 §.§ Solution of the scalar wave equation in Bondi gauge





We now find explicit solutions of the scalar wave equation for
δψ, and of the coupled equations for δ f and δ
b in linear Bondi gauge, using separation of variables. We use Bondi
gauge R=x in the background, and to indicate this we write r for
x. At the end we translate the results back into sM background
gauge. We begin in this subsection with the wave equation.

The linear wave equation on flat spacetime in Bondi coordinates is

    -2δψ_,ur+δψ_,rr+2 r(δψ_,r-δψ_,u)
    +1 r^2[(1-y^2)δψ_,y]_,y=0.

Replacing the retarded time coordinate u with the Minkowski time coordinate

    t:=u+r,

and writing

    δψ(u,r,y)=:δϕ(t,r,y)=δϕ(u+r,r,y),

this takes the more familiar form 

    -δϕ_,tt+δϕ_,rr+2
      rδϕ_,r+1 r^2[(1-y^2)δϕ_,y]_,y=0.

We initially work in coordinates (t,r,y) and switch back to
(u,r,y) later.





  §.§.§ Separating off the y-dependence




Spherical symmetry of the background allows us to separate the angle
y with the ansatz

    δϕ(t,r,y)=:∑_l=0^∞δϕ_l(t,r)P_l(y),

where P_l(y) is the Legendre polynomial of order l, obeying

    L^2P_l:=(1-y^2)P_l”-2yP_l'=-l(l+1)P_l.

We then have the one-dimensional wave equation

    -δϕ_l,tt+δϕ_l,rr+2 rδϕ_l,r-l(l+1) r^2δϕ_l=0,

for the l-th partial wave, or equivalently

    -2δψ_l,ur+δψ_l,rr+2 r(δψ_l,r-δψ_l,u)
    -l(l+1) r^2δψ_l=0,

where

    δψ(u,r,y)=:∑_l=0^∞δψ_l(u,r)P_l(y),






  §.§.§ General solution of l-th partial wave equation





We can find the general solution of (<ref>) that is regular at
r=0 in terms of a single free function χ_l of
one variable, in the form of a generalised d'Alembert solution
<cit.>, as

    δϕ_l(t,r)   =   ∑_p=0^l A_l^p r^-p-1[χ_l^(l-p)(t-r) 
              -(-1)^l-pχ_l^(l-p)(t+r)],

where 

    A_l^p:=(l+p)! 2^p p! (l-p)!,

and χ^(n) denotes the n-th derivative of χ.

Expanding χ(t± r) into its Taylor series about r=0, we obtain the
double sum

    δϕ_l(t,r)=2∑_p=0^l∑_q=0 
     l-p+q  odd^∞(-1)^q A_l^p q! r^q-p-1χ_l^(l-p+q)(t).

Noting that l-p+q≥ 0 and odd in all terms, we parameterise q as
q=2k+1-l+p. We then have

    δϕ_l(t,r)   =   2∑_k=0^∞(∑_p= max
    (0,l-2k-1)^l(-1)^l-pA_l^p (2k+1-l+p)!)
    
           r^2k-lχ_l^(k)(t) .

The inner sum (in round brackets) vanishes for k<l, and so the outer
sum starts up only at k=l. Hence, the expression
(<ref>) admits a formal expansion in even(odd) powers of
r at constant t, for l even(odd), starting at r^l. In
particular, δψ_l(t,-r)=(-1)^lδψ_l(t,r), and
δψ_l∼ r^l at the origin.

Switching back from (t,r,y) to (u,r,y), we trivially find
(<ref>) with

    δψ_l(t,r)   =   ∑_p=0^l A_l^p r^-p-1[χ_l^(l-p)(u) 
              -(-1)^l-pχ_l^(l-p)(u+2r)],

Similarly, by setting t=u+r, in (<ref>), we
trivially obtain

    δψ_l(t,r)   =   2∑_k=0^∞(∑_p= max
    (0,l-2k-1)^l(-1)^l-pA_l^p (2k+1-l+p)!)
    
           r^2k-lχ_l^(k)(u+r) .






 §.§ Relating the perturbed Einstein equations in Bondi gauge
  to the scalar wave equation, in spherical harmonics





In the vacuum Einstein equations in Bondi gauge, linearised about flat
spacetime, the coupled equations for δ b and δ f
(<ref>,<ref>) become

    r^2δ b_,rr+4r δ b_,r-2(1-y^2)δ f_,r
      y+8yδ f_,r   =   0, 
    
    
    r(4δ f_,ur-2δ f_,rr)+4(δ f_,u-δ
      f_,r)       
    
    -r δ b_,r y-2δ b_,y   =   0,

where we have used x=R_0=:r and Ξ_0 R_0=-1/2. These equations
are equivalent to Eqs. (6,7) of
<cit.>. In that paper solutions of these
coupled equations were related to solutions of the scalar wave
equation. Here we give a self-contained derivation of that relation.





  §.§.§ Separating off the y-dependence




Spherical symmetry of the background means we can make the separation
of variables ansatz

    δ b   =   δ b_l(u,r)P^(1)_l(y), 
    δ f   =   δ f_l(u,r)P^(2)_l(y),

with the P^(1,2)_l(y) to be found below.
We obtain the four separated equations

    P_l^(1)'   =   P_l^(2),
    
    (1-y^2)P_l^(2)'-4yP_l^(2)   =   λ P_l^(1), 
    
    (r^4δ b_l,r)_,r   =   2λ r^2δ f_l,r, 
    
    2r(r δ f_l,u)_,r-(r^2δ f_l,r)_,r   =   1 2(r^2δ b_l)_,r,

where, without loss of generality, we have fixed a separation
constant in (<ref>) to one. Clearly (<ref>) and (<ref>)
together give

    (1-y^2)P_l^(1)”-4yP_l^(1)'=λ P_l^(1).







  §.§.§ Tensor spherical harmonics




We would like to relate the unfamiliar spherical harmonic-like
functions P^(1)_l and P^(2)_l to the scalar spherical harmonics
P_l. To motivate this, we note, following <cit.>,
that a geometrically natural ansatz for a perturbation δ
g_ij of the round unit metric g_ij^0 on S^2 with covariant
derivative ∇̅_i is

    δ g_ij=g_1Y_lm g_ij^0+g_2 ∇̅_i∇̅_jY_lm

where we have gone beyond axisymmetry but have restricted to polar
perturbations, and where Y_lm is a scalar spherical
harmonic. Similarly, the perturbation of any vector on S^2 such as
the shift must take the form

    δβ^i=g_3 g^ij_0∇̅_jY_lm.

Here g_1,2,3 are functions of u and r.  In axisymmetry, the
spherical harmonic Y_lm reduces to the Legendre polynomial
Y_l0=P_l. Comparing (<ref>,<ref>) with our metric ansatz
(<ref>) and separation of variables ansatz
(<ref>,<ref>) suggests

    P^(1)_l= P_l',    P^(2)_l=P_l”,

consistent with (<ref>). Substituting (<ref>) into
(<ref>), and using (<ref>) and its y-derivative, we find
that (<ref>) is indeed obeyed with the separation constant given
by

    λ=-l(l+1)+2 = -(l+2)(l-1).

We note in passing that P_l^(2) obeys the ODE

    (1-y^2)P_l^(2)”-6P_l^(2)'=λ̃P_l^(2),

where

    λ̃:=-l(l+1)+6=-(l+3)(l-2).

We are now left with the two coupled PDEs (<ref>,<ref>) for
δ b_l and δ f_l, with λ given by (<ref>).




  §.§.§ Potential ansatz




We can solve (<ref>) identically by introducing the potential
Ψ_l, in terms of which

    δ b_l(u,r)   =    2λ∫_r_l(u)^r Ψ_l(u,r̅)r̅^4dr̅, 
    δ f_l(u,r)   =   ∫_r_l(u)^rΨ_l,r(u,r̅)r̅^2dr̅.

r_l(u) is arbitrary, but will later be set to zero.
Substituting this into (<ref>), dividing by r and
differentiating with respect to r in order to eliminate the
integrals, we obtain the third-order PDE E_3=0, where

    E_3:=-2r^3Ψ_l,urr+r^3Ψ_l,rrr-r^2Ψ_l,rr+λ(rΨ_l,r-Ψ).

We can write

    E_3=r^2(E_2 r)_,r,

where we have defined

    E_2:=-2r^2Ψ_l,ur+r^2Ψ_l,rr+2rΨ_l,u-2rΨ_l,r+λΨ_l.

Hence the general solution of E_3=0 is 

    E_2=Cr,

where C is an arbitrary constant, and the general solution of this
is in turn

    Ψ_l^(3)(u,r)=Ψ_l^(2)(u,r)+Cr  l(l+1),

where Ψ_l^(2) is the general solution of the second-order ODE
E_2=0. With the substitution

    Ψ_l=:r^2δψ_l,

E_2=0 for Ψ_l becomes the scalar wave equation (<ref>)
for δψ_l.

Substituting (<ref>) into (<ref>,<ref>), we have

    δ b_l(u,r)   =   2λ∫_r_l(u)^r δψ_l(u,r̅)r̅^2dr̅, 
    δ f_l(u,r)   =   δψ_l(u,r)-δψ_l(u,r_l(u))
    
          
    +2∫_r_l(u)^rδψ_l(u,r̅)r̅ dr̅.


Note that, in addition to an arbitrary solution δψ_l of the scalar
wave equation, the gravitational wave (δ f_l,δ b_l) is parameterised
also by the constant C and the function r_l(u). The particular
integral parameterised by the constant C is both non-dynamical and
gives rise to δ b and δ f that are singular at r=0. We therefore set
C=0, and so obtain E_2=0 as the master equation governing
(axisymmetric, polar) linear gravitational waves. We believe that
r_l(u) parameterises some kind of gauge freedom, but have not tried
to show this. We now set it to zero. 

Using the fact that δψ_l∼ r^l at constant t, we see from
(<ref>,<ref>) that near the centre, each spherical harmonic
component of the linear   gravitational wave has the spatial dependence
δ b∼ P_l'r^l-1 and δ f∼ P_l” r^l, for l≥ 2. In
particular, near the origin r=0 we therefore have δ b∼ ry and
δ f∼ r^2 for generic regular solutions with non-vanishing l=2
components.

Substituting the general regular solution (<ref>) for
δψ_l(u,r) with a sufficiently simple function χ_l(u) into
(<ref>,<ref>), and with r_l(u)=0, we can carry out the
integrations and so obtain δ b_l(u,r) and δ f_l(u,r) in closed
form. In particular, this works when χ_l(u) is a Gaussian, and
gives us a class of exact solutions for testing our code.

The general solution of (<ref>,<ref>) is
obtained by summing

    δ f(u,r,y)   =   ∑_l=2^∞δ f_l(u,r)P_l”(y), 
    δ b(u,r,y)   =   ∑_l=2^∞δ b_l(u,r)P_l'(y).

If we now substitute (<ref>) and (<ref>), and then
(<ref>), we see that at constant u, δ f_l admits a
formal expansion in even (odd) powers of r when l is is even
(odd), starting at r^l. This is just as for δψ_l. Similarly,
δ b_l admits a formal expansion in odd (even) powers of r when l
is is even (odd), starting at r^l-1.

Obviously, the scalar field δψ governing the linearised
gravitational waves is completely independent from the matter scalar
field δψ. We have used the same notation merely to stress that
they obey identical wave equations. 





  §.§.§ The cases l=0 and l=1




Our separation of variables ansatz gives b_0=f_0=f_1=0 because
P_0'=P_0”=P_1”=0. In addition the potential ansatz (<ref>)
gives δ b_1=0 because λ=0. However, as P_1=y and
P_1'=1, there can be nonvanishing perturbations δ b_1, δ
H_1, even though there is no f_1.

We therefore look at the case l=1 of the linearised equations in
Bondi gauge without making the potential ansatz
(<ref>,<ref>). With δ b=δ b_1(u,r) and δ f=0,
(<ref>) reduces to

    (r^4δ b_1,r)_,r=0,

while (<ref>) is obeyed trivially. The solution of
(<ref>) that is finite at the centre is constant in r. 
The resulting regular perturbation of H is given by
(<ref>). Putting both together, we have 

    δ b_1(u,r)=δ b_1(u),     δ H_1=2rδ b_1(u).

We show in Appendix <ref>
that this represents a linear gauge transformation, which
physically corresponds to an acceleration -δ b_1(u) of the
origin of our coordinate system along the symmetry axis.





  §.§.§ Transformation to shifted Minkowski background gauge




Having found δ b(u,r,y), δ f(u,r,y) and
δψ(u,r,y) in Bondi background gauge, we transform them to
shifted Minkowski background gauge, with the coordinate change from
(u,r,y) to (u,x,y) given by

    r=x 2( 1-u x_0),

compare (<ref>). Under this coordinate transformation, only
δ H changes nontrivially. The metric coefficients δ b,
δ f and R and the scalar field δψ change only as if
they were scalar fields, and δ G remains zero. We are only
interested in δ b, δ f and δψ, and we only
need to transform their argument r to x using
(<ref>).





  §.§.§ Consistent boundary conditions 
  for the linearised Einstein equations decomposed into spherical
  harmonics




We now return to the case of a non-trivial r_l(u)>0.
From (<ref>) and (<ref>) we read off that 

    δ b_l(u,r_l(u))   =   0, 
    δ f_l(u,r_l(u))   =   0, 
    δ b_l,r(u,r_l(u))   =   2λ r_l(u)^2δψ_l(u,r_l(u)), 
    δ f_l,r(u,r_l(u))   =   δψ_l,r(u,r_l(u))+2
      r_l(u)δψ_l(u,r_l(u)).

The first three equations gives us a class of startup conditions for
integrating the linearised hierarchy equations for δ b_l and
δ f_l that have a consistent continuum limit, in terms of a
boundary condition at r=r_l(u) for δ b_l,r, which is
equivalent to a boundary condition on the underlying scalar wave
δψ_l.

The simplest choice is to set δ b_l,r=0 at r=r_l(u). As
this corresponds to the homogeneous Dirichlet boundary condition
δψ_l=0 at r=r_l(u) on the underlying scalar field, this
boundary condition gives rise to a well-posed PDE problem for the
spherical harmonic component ψ_l(u,r). We have therefore shown
that, for the equations linearised about flat spacetime in linear
Bondi gauge, it is consistent to impose the three boundary conditions
δ f_l=δ b_l=δ b_l,r=0 at r=r_l(u). 

This is of course what we do, as a numerical trick, for the full
nonlinear Einstein equations, with r_l∼ lΔ x. What we have
shown here is that this has a continuum limit for the Einstein
equations linearised about Minkowski and split into spherical harmonics.





 §.§ Solution without spherical harmonic decomposition









  §.§.§ Relation between solutions for δψ and for (δ f,δ b)




We can remove the spherical harmonic decomposition, and obtain
expressions for δ b(u,r,y) and δ f(u,r,y) in terms of a
solution δψ(u,r,y) of the scalar wave equation.  Combining
(<ref>), (<ref>) and (<ref>), we obtain the
unseparated expression

    δ f(u,r,y)   =   δψ_,yy(u,r,y)-δψ_,yy(u,0,y) 
          
    +2∫_0^rδψ_,yy(u,r̅,y)r̅ dr̅

for δ f in terms of δψ.  Substituting
(<ref>) into (<ref>) and integrating twice, we
obtain the expression

    δ b(u,r,y)   =   ∫_0^r1r̅^4∫_0^r̅r̃^2(
    2(1-y^2)δ f_,ry(u,r̃,y)
          -8yδ f_,r(u,r̃,y)) dr̃

for δ b in terms of δ f. 





  §.§.§ Plane-wave solutions





The scalar wave equation admits the plane-wave solution
δψ=χ(t± z), that is 

    δψ(u,r,y)=χ(u_±),      u_±:=u+r(1± y),

for arbitrary functions χ. Substituting this into
(<ref>), we can carry out the integration explicitly
to obtain the plane wave solution 

    δ f(u,r,y)=r^2χ”(u_±)+2rχ'(u_±) 1± y
    -2χ(u_±)-χ(u) (1± y)^2,

and substituting this into (<ref>), we can again carry
out the integrations in closed form to obtain 

    δ b(u,r,y)   =   ± 2[r(1∓ y)χ”(u_±) 
           -1± 3y
      1± y(χ'(u_±)-χ'(u))].

Note that ψ(u,r,-y)=ψ(u,r,y), f(u,r,-y)=f(u,r,y) and b(u,r,-y)=-b(u,r,y).

To see that these expressions are regular as y→∓ 1, we define
the auxiliary variable ϵ:=r(1± y). We can then write

    δψ(u,r,y)   =   χ(u+ϵ), 
    δ f(u,r,y)   =   r^2χ”(u+ϵ)+2r^2d
      dϵ(χ(u+ϵ)-χ(u)ϵ), 
    δ b(u,r,y)   =   ± 2r[(1∓ y)χ”(u+ϵ) 
          -(1± 3y)(χ(u+ϵ)-χ(u)ϵ)],

where we can now see that the fraction in large round brackets is
regular as ϵ→ 0 if χ(u) is regular.





§ RESIDUAL GAUGE FREEDOM





A redefinition of the coordinates (u,x,y) leaves the form
(<ref>) of the metric invariant if g_xx=g_xy=0 in the
new, as in the old, coordinates. A nonlinear gauge transformation that
does this and that can be written explicitly in terms of free
functions is

    u   =   u(ũ), 
    
    y   =   y(ũ,ỹ), 
    
    x   =   x(ũ,x̃,ỹ).

This means we separately relabel the outgoing null cones N^+_u, the outgoing null rays L^+_u,y,φ on each
outgoing cone, and the coordinate x along each outgoing
ray. However, we cannot express the most general gauge freedom, which
also moves the central worldline, in explicit form.

Instead, we look for the most general infinitesimal gauge
transformation that leaves the form of the metric invariant.  We make
the ansatz δ x^μ=:ξ^μ, and hence

    δ g_μν=2∇_(μξ_ν),

where

    ξ^μ:=(ξ^u,ξ^x,ξ^y,0),

are functions of (u,x,y). The general solution of δ
g_xx=δ g_xy=0 is

    ξ^u   =   A(u,y), 
    ξ^x   =   ξ^x(u,x,y), 
    ξ^y   =   B(u,y)+
    A_,y(u,y) S∫_0^x G e^2S fR^2 dx'.

For A=A(u) this reduces to the infinitesimal version of
(<ref>)-(<ref>). 

For simplicity, we now restrict the background metric to flat
spacetime in Bondi gauge, with R=x, G=H=1 and b=f=0. (<ref>)
then simplifies to 

    ξ^y=B(u,y)-SA_,y(u,y) x.

The resulting metric perturbation is 

    δ G   =   A_,u+ξ^x_,x,
    δ H   =   2(A_,u+ξ^x_,u), 
    δ R   =   ξ^x+xB_,y 2-SA_,yy-2yA_,y 2, 
    δ f   =   -A_,yy 2x+B_,y+2yB 2S, 
    δ b   =   -A_,y x^2+B_,u S-xA_,uy+ξ^x_,y x^2.


For analysis, we relate the coordinates (u,x,y) to the
cylindrical Minkowski coordinates defined by

    t:=u+x,      z:=-yx,     ρ:=√(S)x,

with φ completing both coordinate systems. The inverse
transformation is

    u   =   t-√(ρ^2+z^2), 
    
    x   =   √(ρ^2+z^2),
    
    y   =   -z√(ρ^2+z^2),

and hence the coordinate basis vectors transform as

    ∂_t   =   ∂_u, 
    ∂_z   =   y∂_u-y∂_x-S x∂_y, 
    ∂_ρ   =   -√(S)∂_u+√(S)∂_x-y√(S) x∂_y,

with ∂_φ in both coordinate systems. 

From (<ref>-<ref>) and (<ref>-<ref>), the
residual gauge transformations expressed in this basis are

    ξ   =   (A+ξ^x) ∂_t+(√(S)(ξ^x+yA_,y)
    -xyB√(S))∂_ρ
          +(SA_,y-xB-yξ^x) ∂_z.

ξ is a regular vector field if and only if the coefficients of
∂_t,∂_z,∂_ρ are regular functions of
(t,z,ρ). In particular, they must be finite and single-valued at
the origin x=0 and on the axis S=0.

We now focus on gauge transformations that map the tz-plane to
itself, as these include all that move the central worldline
along the symmetry axis. Setting the component
ξ^ρ to zero is equivalent to

    ξ^x=xyB S-yA_,y,

and substituting this back into (<ref>), we have 

    ξ=(A-yA_,y+xyB S)∂_t
    +(A_,y-xB S)∂_z.

The resulting metric perturbations are regular at the origin and on
the axis if and only if

    A=T(u)+Z(u) y,      B=-SZ'(u),

where T and Z are arbitrary functions. The corresponding regular
gauge vector field is

    ξ=[T(u)-z Z'(u)] ∂_t
    +[Z(u)+x Z'(u)] ∂_z,

which explains why we have named the free coefficients T and
Z. The corresponding pure gauge metric perturbations are

    δ H   =   2(T'(u)+zZ”(u)), 
    δ G   =   T'(u), 
    δ b   =   -Z”(u),

with δ R=δ f=0. This family includes the three
Killing vector fields ∂_t, ∂_z and
z ∂_t+t ∂_z of Minkowski spacetime.

Restricting to T(u)=0 leaves us with

    ξ=Z'(u) z ∂_t+[Z(u)+x Z'(u)] ∂_z.

This is equal to Z(u)∂_z at the origin, and so is the gauge
transformation that moves the origin along the symmetry axis. The
resulting metric perturbation is

    δ H=-2Z”(u) x y,     δ b=-Z”(u),

with the other metric coefficients unchanged. Hence an infinitesimal
non-inertial motion of the origin along the symmetry axis gives
exactly the l=1 regular metric perturbation (<ref>)
that we already derived from the linearised Einstein equations.

From elementary flatness at the origin, the same statement must be
true for an arbitrary regular, twistfree axisymmetric
spacetime. Hence we have shown that b(u, 0)=0 is precisely the
gauge condition that forces the origin of our coordinate system to
move on a geodesic along the z-axis. We are imposing this as a
boundary condition on (<ref>) or (<ref>). Conversely,
by imposing an arbitrary value of b(u, 0), we could accelerate
the origin during the evolution.





§ REGULARITY OF THE METRIC AT THE ORIGIN





Following common practice, Rinne <cit.> defines a scalar
field as regular if it is analytic in Cartesian coordinates
(t,ξ,η,z) (which in turn we assume to give rise to a chart
without coordinate singularities). Transforming to the cylindrical
coordinates (t,ρ,φ,z) defined by

    ξ:=ρcosφ,     η:=ρsinφ,

with t and z unchanged, Rinne shows that a scalar field is regular
and axisymmetric if and only if it is an analytic function of
(t,ρ^2,z), that is, analytic in (t,ρ,z), with only even
powers of ρ.

Rinne then shows that an axisymmetric metric is regular, in the sense
that its components in the Cartesian coordinate basis are analytic
functions of the Cartesian coordinates, if and only if in the
cylindrical coordinate basis it takes the form

    ([             A             B           ρ D         ρ^2 F;             B             C           ρ E         ρ^2 G;           ρ D           ρ E       H+ρ^2 J         ρ^3 K;         ρ^2 F         ρ^2 G         ρ^3 K ρ^2( H-ρ^2 J);               ]),
 with the coordinates in order (t,z,ρ,φ),
and the “Rinne coefficients” A, B, C,
D, E, F, G, K, H, J regular, that is analytic functions of
(t,z,ρ^2). Clearly, the metric is twist-free if and only if F= G= K=0, and we now restrict to this case.

In the following, we define a scalar field and metric as
“Rinne-regular” if they obey these definitions. We shall now ask
what a Rinne-regular metric and scalar field look like, first in
spherical polar coordinates retaining the regular time slicing t,
and secondly replacing t by a retarded time u. 

In the first step, we change from the cylindrical coordinates
(t,ρ,z,φ) to the spherical coordinates (t,x,y,φ)
defined by

    ρ:=x√(1-y^2),      z:=-xy,

or equivalently

    x=√(ρ^2+z^2),    y=-z/√(ρ^2+z^2),

with t and φ unchanged.  Note that x and y are our
radial and angular coordinates, not the Cartesian coordinates
(ξ,η,z) defined above.

With ρ^2=x^2-z^2, a function is axisymmetric and regular if and
only if it is an analytic function of (t,x^2-z^2,z), but this is the
case if and only if it is an analytic function of (t,x^2,z) or,
expressed entirely in spherical coordinates, of (t,x^2,-xy). In
particular, this means it is even in x at constant z, not constant
y: this will be important below.

By transforming to the coordinate basis with respect to
(t,x,y,φ), and reparameterising the Rinne coefficients
as

    B   =   z D- X_1, 
     E   =   2z J- X_2-z X_3,
     C   =    H+(2z^2-ρ^2) J-2z
      X_2+(ρ^2-z^2) X_3,

we can show that the Rinne-regular twist-free axisymmetric metric is
spherically symmetric if and only if X_1= X_2=
  X_3=0.

Moreover, because the coefficients of the redefinitions
(<ref>-<ref>) themselves are Rinne-regular
functions, and because we cannot absorb them into redefinitions of
X_1,2,3 without making those singular, the reparameterised
twist-free axisymmetric is Rinne-regular if and only if the
“modified Rinne coefficients” A, D, H,
J, X_1,2,3 are Rinne-regular.

In the second step, we change from the spherical coordinates
(t,x,y,φ) to the retarded coordinates (u,x,y,φ),
where the retarded time
u is defined as

    u(t,x,y):=t-h(t,x,y).

Without loss of generality, we set h(t,0,y)=0, and we split h
into its odd and even in x (at constant z) parts as

    h(t,x,y):=x h_o(t,x^2,-xy)+x^2 h_e(t,x^2,-xy).

For u to be a non-trivial retarded time, we assume that h_o
0. Hence h is not Rinne-regular even if h_o and h_e are. The
simplest choice would be u=t-x, the standard retarded null
coordinate on Minkowski spacetime. In this case, we can invert the
definition of u to give t=u+x, and so a Rinne-regular axisymmetric
function would be an analytic function of (u+x,x^2,-xy). 

We now impose that the lines of constant (u,y,φ) are null,
which is equivalent to g_xx=g_xy=0. These are linear equations
relating the modified Rinne coefficients, but their coefficients are
neither even nor odd in x (at constant z), so they have
Rinne-regular solutions only if we impose their even and odd parts
separately, giving us four equations. Hence we see that beyond
spherical symmetry we need the two generic Rinne-regular functions
h_e and h_o in order to be left with five regular functions to
match to our metric coefficients (G,H,R,f,b).

Which four of A, D, H, J,
X_1,2,3, h_e and h_o should we solve for? In spherical
symmetry the simple height function h=x is already sufficiently
general. This suggests that we should always solve for four of the
seven modified Rinne coefficients, rather than the height function,
and that this solution should be regular also in the special case
h=x. These requirements force us to solve for D,
H, X_1 and X_2.

We can now express our five metric coefficients in terms of the five
remaining arbitrary regular functions A, J, X_3, h_e and h_o. To avoid square roots and logarithms in the
solution, we solve not for R and f but for R^4 and
e^4Sf. The resulting expressions are not immediately
helpful, as they are explicit functions of (t,x,y), not (u,x,y),
and we do not have an explicit expression for t(u,x,y).

However, we can expand (G,H,R,f,b) in powers of x at
constant (u,y) by using the derivative operator

    .∂∂ x|_u,y=
    .∂∂ x|_t,y
    -u_,x(t,x,y) u_,t(t,x,y).∂∂ t|_x,y

to find the coefficients of the Taylor series. We are then in effect
inverting the height function order by order. The Taylor coefficients are
given in terms of A, J, X_3, h_e and
h_o and their partial derivatives with respect to (t,x^2,z), 
evaluated at x=z=0, which are finite and independent of each other. 

We do not give the series expressions here, but it is clear from their
construction that G, H, R/x, f/x^2 and b are analytic
functions of (x^2,-xy), at constant u, with the coefficients of
the Taylor series given explicitly in terms of the coefficients of the
Taylor series in (x^2,-xy), at constant t=u, of A,
J, X_3, h_e. This confirms the ansatz that we
found to be consistent for expanding the field equations in powers of
x at constant u in Sec. <ref>.



99



BondiBurgMetzner1962 H. Bondi, M. G. J. van der Burg, and
  A. W. K. Metzner, Gravitational Waves in General
  Relativity. VII. Waves from axi-symmetric isolated
  systems. https://www.jstor.org/stable/2414436?seq=9Proc.
    Roy. Soc. A, 269, 21 (1962).

Sachs1962 R. K. Sachs, Gravitational waves in general
  relativity VIII. Waves in asymptotically flat spacetime,
  https://royalsocietypublishing.org/doi/10.1098/rspa.1962.0206Proc. 
Roy. Soc. A 270, 103 (1962).

paper2 C. Gundlach, T. W. Baumgarte and D. Hilditch,
  Simulations of gravitational collapse in null coordinates:
  II. Critical collapse of an axisymmetric scalar field.

GundlachLRR C. Gundlach and J. M. Martín-García,
  Critical phenomena in gravitational collapse,
  https://link.springer.com/article/10.12942/lrr-2007-5Living
    Rev. Relativ. 07, 05 (2007).

Vacuum_Collapse T. W. Baumgarte, B. Brügmann, D. Cors
  Agulló, C. Gundlach, D. Hilditch, A. Khirnov, T. Ledvinka,
  S. Renkhoff, and I. Suárez Fernández, Critical phenomena in the
  collapse of gravitational waves,
  https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.131.181401Phys. Rev. Lett. 131, 181401 (2023).

paper3 C. Gundlach, Simulations of gravitational collapse
  in null coordinates: III. Hyperbolicity.

Winicour2005 J. Winicour, Characteristic evolution and
  matching, http://www.livingreviews.org/lrr-2005-10Living
    Rev. Relativ. 5, 10 (2005).

BishopRezzolla2016 N. T. Bishop and L. Rezzolla, Extraction of
  gravitational waves in numerical relativity, https://link.springer.com/article/10.1007/s41114-016-0001-9Living
  Rev. Relativ. 19, 2 (2016).
  
GomezPapadopoulosWinicour1994 R. Gómez, P. Papadopoulos,
  and J. Winicour, Null cone evolution of axisymmetric vacuum
  spacetimes, https://aip.scitation.org/doi/abs/10.1063/1.530848J. Math. Phys. 35, 4184 (1994).

SiebelFontMuellerPapaodopoulos2002 F. Siebel, J.A. Font,
  E. Müller and P. Papadopoulos, Simulating the dynamics of relativistic
  stars via a light-cone approach, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.65.064038Phys. Rev. D 65, 064038 (2002).

BishopGomezLehnerMaharajWinicour1997 N. T. Bishop,
  R. Gómez, L. Lehner, M. Maharaj and J. Winicour, High-powered
  gravitational news,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.56.6298Phys. Rev. D
    56, 6298 (1997).

BishopGomezLehnerMaharajWinicour1999 N. T. Bishop,
  R. Gómez, L. Lehner, M. Maharaj and J. Winicour, Incorporation of
  matter into characteristic numerical relativity, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.60.024005Phys.  Rev. D 60, 024005
  (1999). 

ReisswigBishopPollney2013 C. Reisswig, N. T. Bishop and
  D. Pollney, General relativistic null-cone evolutions with a
  high-order scheme, 
  https://link.springer.com/article/10.1007/s10714-013-1513-1Gen. Rel. Gravit. 45, 1069
  (2013).  

VanderWaltBishop2012 P. J. van der Walt and N. T. Bishop,
  Observational cosmology using characteristic numerical relativity:
  characteristic formalism on null geodesics, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.85.044016Phys. Rev. D 85,
  044016
  (2012). 

CrespoOliveiraWinicour2019 J. A. Crespo, H. P. de Oliveira and
  J. Winicour, The affine-null formulation of the gravitational
  equations: spherical case, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.100.104017Phys. Rev. D 100, 104017 (2019). 

Maedler2019 T. Mädler, Affine-null metric formulation of
  general relativity at two intersecting null hypersurfaces,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.99.104048Phys. Rev. D
    99, 104048 (2019).

Winicour2013 J. Winicour, The affine-null metric formulation
  of the gravitational equations, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.87.124027Phys. Rev. D 87, 124027
  (2013).

CheslerYaffe2014 P. M. Chesler and L. G. Yaffe, Numerical
  solution of gravitational dynamics in asymptotically anti-de Sitter
  spacetimes,
  https://link.springer.com/article/10.1007/JHEP07(2014)086JHEP
    2014, 86 (2014).

CheslerYaffe2011 P. M. Chesler and L. G. Yaffe, Holography
  and colliding gravitational shock saves in asymptotically AdS5
  spacetime,
  https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.106.021601Phys. Rev. Lett. 106, 021601 (2011).

CheslerLowe2019 P. M. Chesler and D. A. Lowe, Nonlinear
  evolution of the AdS_4 black hole bomb, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.122.181101Rev. Lett. 122,
  181101
  (2019). 

Chesler2022 P. M. Chesler, Hairy black resonators and the
  AdS_4 superradiant instability, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.105.024026Phys. Rev. D 105, 024026
  (2022). 

GoldwirthPiran1987 D. S. Goldwirth and T. Piran,
  Gravitational collapse of massless scalar field and cosmic
  censorship, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.36.3575Phys. Rev. D 36, 3575 (1987).

GundlachPricePullin1994 C. Gundlach, R. Price and J. Pullin,
  Late-time behavior of stellar collapse and explosions: II. Nonlinear
  evolution, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.49.890Phys. Rev. D 49, 890 (1994).

HamadeStewart1996 R. S. Hamade and J. M. Stewart, The
  spherically symmetric collapse of a massless scalar field,
  https://iopscience.iop.org/article/10.1088/0264-9381/13/3/014Class. Quant. Grav. 13, 497 (1996). 

PuerrerHusaAichelburg M. Pürrer, S. Husa and
  P. C. Aichelburg, News from critical collapse: Bondi mass, tails,
  and quasinormal modes,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.71.104005Phys. Rev. D
    71, 104005 (2005).

Choptuik1993 M. W. Choptuik, Universality and scaling in
  gravitational collapse of a massless scalar field,
  https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.70.9Phys. Rev. Lett. 10, 1103 (1993).

Garfinkle1995 D. Garfinkle, Choptuik scaling in null
  coordinates,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.51.5558Phys. Rev. D
    51, 5558 (1995).

DafermosHolzegelRodnianskiTaylor2021 M. Dafermos,
  G. Holzegel, I. Rodnianski and M. Taylor, The non-linear stability
  of the Schwarzschild family of black holes (2021), https://arxiv.org/abs/2104.08222arXiv:2104.08222.

PortoGundlach2022 B. Porto-Veronese and C. Gundlach,
  Critical phenomena in gravitational collapse with competing scalar
  field and gravitational waves in 4+1 dimensions,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.106.104044Phys. Rev. D
    106, 104044 (2022).

Rinne2020 O. Rinne, Type II critical collapse on a single
  fixed grid: a gauge-driven ingoing boundary method, https://link.springer.com/article/10.1007/s10714-020-02768-xGen. Rel. and Grav. 52, 117
  (2020). 

ymscalar C. Gundlach, T. W. Baumgarte and D. Hilditch, Critical
  phenomena in gravitational collapse with two competing massless
  matter fields, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.100.104010Phys. Rev. D 100, 104010 (2019).

SzabadosLRR L. B. Szabados, Quasi-Local Energy-Momentum and
  Angular Momentum in GR,
  https://link.springer.com/article/10.12942/lrr-2004-4Living
    Rev. in Relativ.  4, 4 (2004).

EardleyBattelle D. M. Eardley, in Sources of
  Gravitational Radiation (Proceedings of the Battelle Seattle
  Workshop, July 24 - August 8, 1978), edited by L. L. Smarr,
  Cambridge University Press, 1979, 
  https://www.jstor.org/stable/43852846JSTOR.

Boyd J. P. Boyd, Chebyshev and Fourier Spectral Methods
  (Second Edition),
  Dover Publications, New York
  2000, https://depts.washington.edu/ph506/Boyd.pdfaccessed online.

SBP C. Gundlach, J. M. Martín-García,
  D. Garfinkle, Summation by parts methods for the spherical harmonic
  decomposition of the wave equation in arbitrary dimensions,
  https://iopscience.iop.org/article/10.1088/0264-9381/30/14/145003?pageTitle=IOPscienceClass. Quantum Grav. 30, 145003 (2013).

Jietal2023 L. Ji, V. Mewes, Y. Zlochower, L. Ennoggi,
  F. G. Lopez Armengol, M. Campanelli, F. Cipolletta, and
  Z. B. Etienne, Ameliorating the Courant-Friedrichs-Lewy condition in
  spherical coordinates: A double FFT filter method for general
  relativistic MHD in dynamical spacetimes,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.108.104005Phys. Rev. D
    108, 104005 (2023).

Price R. H. Price, Nonspherical perturbations of
  relativistic gravitational collapse. I. Scalar and gravitational
  perturbations, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.5.2439Phys. Rev. D 5, 2419 (1972).

GerlachSengupta U. H. Gerlach and U. K. Sengupta,
  Relativistic equations for aspherical gravitational collapse,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.18.1789Phys. Rev. D 18, 1789 (1978).

RinnePhD O. Rinne, Axisymmetric general relativity, PhD
  thesis, Cambridge 2005,
  https://arxiv.org/abs/gr-qc/0601064gr-qc/0601064.

Kinoshita S. Kinoshita, On geometrical origin of Kodama
  vector (2024), https://arxiv.org/abs/2402.164842402.16484.








