Simulations of gravitational collapse in null coordinates. I. Formulation and weak-field tests in generalized Bondi gauges

I. INTRODUCTION A. 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 [1,2]. 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 polarizations 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 discretized 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. I B 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 center of approximate spherical symmetry. In a companion paper [3] (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 nonspherical 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 [4] for a general review of critical collapse and [5] 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 [6], we consider issues of hyperbolicity and well-posedness.

B. 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 [7,8].

One natural choice of null surfaces is the set of null cones emerging from a regular central worldline. The cones are labeled 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 [9] and with perfect fluid matter [10]. Here, the null cones emanate from a regular center. This brings about a severe limitation of the time step to Δu ∼ ΔxðΔθÞ 2 , see also Appendix A.

Without symmetry restrictions, Bondi gauge in vacuum has been implemented in [11], and with perfect fluid matter in [12], 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 [13]. 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 [14], and in an application to spherical scalar field critical collapse in [15], and to vacuum in spherical symmetry with initial data on two intersecting null cones in [16]. Affine gauge in vacuum without symmetries was formulated in [17], 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 [18]. Mentioning only the two applications most relevant for us, in [19] 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 [20,21].

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 center, 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 [22]. Essentially the same algorithm was used in [23] 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 [24], Bondi coordinates compactified at future null infinity in [25], and (as already mentioned) affine coordinates in [15].

"Type II" critical collapse is characterized by an arbitrarily large range of spacetime scales, and hence typically requires adaptive mesh refinement in numerical simulations. In fact, the pioneering paper [26] 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 centered 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 [27] 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 echoes in near-critical solutions. Whenever half the x-grid points have fallen into the center, 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 [28]).

However, we can rescue the key idea of [27] if we choose x in such a way that x ¼ 0 is the regular center while x ¼ x 0 is ingoing null (possible in spherical symmetry only) or future spacelike (ingoing faster than light). We first implemented this in spherical symmetry [29], and found that it provides mesh refinement for critical collapse as efficiently as the algorithm of [27].

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 [30].

C. 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 discretization 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 nonspherical scalar field collapse, which can be examined with an arbitrary degree of nonsphericity, whereas vacuum collapse is necessarily very nonspherical. 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. II we review them in n þ 2 spacetime dimensions, without symmetry assumptions.

Starting from Sec. III, 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.

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

Section V describes convergence tests of the full nonlinear equations in a small data regime where the linearization of the equations about Minkowski is a good approximation.

We conclude with a summary and outlook in Sec. VI.

A number of appendixes give details of (Appendix A) the time step problem, null coordinates in (Appendix B) Minkowski spacetime and (Appendix C) spherical symmetry, (Appendix D) exact solutions of the linearized equations that we use as test beds in the small data regime, (Appendix E) the residual gauge freedom in our coordinate choice, and (Appendix F) regularity conditions for the metric at the origin and on the symmetry axis.

A. 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 ∇ a u∇ b u ¼ g uu ¼ 0. The vector field

is obviously null, and obeys

Hence U a is the tangent vector to the affinely parametrized 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

Note there are

SIMULATIONS OF GRAVITATIONAL …. I. FORMULATION AND … PHYS. REV. D 110, 024018 (2024) 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;

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

where

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

Note that the induced metric on the surfaces of constant u, given by the bottom right submatrix of (9), is degenerate

The inverse metric is

where γ ij is 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 (10) that in our coordinates the vector field U a takes the simple form

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 parametrized 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 parametrized 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. Our coordinates and basis vectors are sketched in Fig. 1.

B. Standard radial gauge choices

We see from (9) that g xi ¼ 0 for i ¼ 1…n, and from (10) 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 det γ ij . If we think of this condition as fixing, for example, the metric coefficient H, we have

1. Schematic spacetime picture showing our coordinates and basis vectors. Shown are two null cones of constant u, four spacelike closed 2-surfaces of constant ðu; xÞ, the central worldline x ¼ 0, and outgoing null rays of constant ðu; y; φÞ (the angular coordinate y ¼cos θ is suppressed here). Solid arrows represent the outgoing null vector U ∝ ∂ x , dashed arrows the ingoing null vector Ξ, and dotted arrows the timelike vector ∂ u (which are timelike near the origin but may tip inward to become spacelike further out).

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.

C. The hierarchy of Einstein equations

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

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

R ab is the spacetime Ricci tensor, and we use gravitational units where c ¼ G ¼ 1.

We define the ingoing null vector field

Together Ξ a and U a span the normal space to each S u;x . They are normalized so that Ξ a U a ¼ -1. U a , given in (11), is tangent to the affinely parametrized 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 [2] 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. III below. Given the Ξγ 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" [2]), E uu ¼ 0 and E u i ¼ 0 (the "supplementary conditions" [2]) 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

we can write (17) as

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 .

A. 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 ϵ abcd K a ∇ b K 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 center 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 parametrize the metric under these conditions as

where the metric coefficients ðG; H; R; F; βÞ depend on the coordinates ðu; x; θÞ only. R is the area radius, in the sense that det γ ij ¼ R 2 sin 2 θ, and so the area of S u;x is 4πR 2 . The central worldline is at R ¼ 0.

C. 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

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 ln G R ;x

The Einstein equations E þ ¼ 0 and E -¼ 0 and the scalar wave equation all contain u-derivatives. They can, respectively, be written as

where Ξ is the derivative operator defined in (25), 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 nonprincipal), and the scalar field stressenergy terms at the end

Note that S G given in (34) has a division by R ;x . If we assume that R ;x > 0 everywhere, we can reparametrize the metric variable G by the new variable

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

that is, the derivative with respect to R at constant u and y. Moreover, g is invariant under reparametrizing 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

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

The hierarchy equations now take the form

DðRΞψÞ ¼ Sψ ½R; f; ψ; γ; b -ðΞRÞDψ; ð46Þ 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 Sb ½R; f; ψ; γ ¼ -DðR 2 DðR ;y ÞÞ þ Sb ½R; f; ψ; γ:

The modified source Sb no longer contains third derivatives, and we can explicitly integrate the first term on the right. The source terms are now

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 center.

We show in Appendix B that when the center 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 0 ðuÞ;

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Þ.

G. Choices of radial gauge for critical collapse

We now present possible choices of radial gauge adapted to critical collapse that generalize the ideas of [27] beyond spherical symmetry. In these coordinates, we want x ¼ 0 to be the regular center 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

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

In short, H and therefore B must be nonpositive 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 onesided 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 (56). 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 light cone 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

for some constant parameter x 0 ≤ x max . The Bondi radial shift was defined in (56). In other words, B decreases linearly in x from its value in Bondi gauge at the center 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 (56) gives

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

Therefore, sdn gauge is a continuous version of the repeated regridding of the double-null coordinate v in [27].

We have previously implemented it in spherical symmetry in [29], and found that it works exactly as well as our earlier implementation, in [31], of the original Garfinkle regridding algorithm. The resulting numerical domains are identical for x 0 ¼ x max , but we found in [29,31] 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.

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 (63) is too restricted. Hence we can attempt the more general gauge Rðu; x; yÞ ¼ Rðu; xÞ; ð66Þ by setting

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

to keep the origin regular at x ¼ 0 and R;u ðu; x max Þ ≤ min y ΞRðu; x max ; yÞ ð 69Þ 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 (56) with its spherical part subtracted ("nonspherical Bondi", from now on nsB),

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

In this gauge x 0 is a null surface "on average."

We can add a further term that makes the shift nonpositive everywhere at x 0 , giving

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")

and again add the sdn shift,

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 center, to full sdn at the outer boundary, so that the outer boundary is null everywhere. In other words, we consider

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.

Hawking mass and Hawking compactness

Following [32], 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, normalized such that l a n 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

and is unique. We then define the null congruence expansions

It is easy to see that

where the first equality holds when n a is null and the second when l a is an affinely parametrized 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 parametrized geodesics. Then the product

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

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

where A ≔ R 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 outertrapped 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 uu 0 ¼ 0;

x-x 0 ðyÞ ¼ 0:

A basis of tangent vectors to S 0 is given by

The outgoing and ingoing null vectors orthogonal to S 0 , normalized so that n a l a ¼ -1 and l a ∇ a u ¼ 1, are

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

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

for later reference.

The induced metric on S 0 is

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

We therefore have

and

Note that the integral

After an integration by parts to eliminate the x 00 0 term, C becomes

where

If we now further restrict to surfaces S u;x of constant u and x, (91) simplifies to

and (97) simplifies to

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 nondecreasing with x, and non-negative. To prove this, note that for R ¼ Rðu; xÞ, the integral for A is trivial, giving ffiffiffiffiffiffiffiffiffiffi ffi A=4π p ¼ R. We can pull this factor of R into the integral for Cðu; xÞ to obtain

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

Hence a sufficient condition for DM ≔ M ;x =R ;x to be nonnegative is that ρ þ ρ -< 0. (This result is a special case of Eardley's observation [33] 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 DM ≥ 0 and M ¼ 0 along the central worldline we then obtain M ≥ 0.

In lsB gauge, where Mðu; xÞ ¼ Rðu; xÞ 2 Cðu; xÞ ð 106Þ we can therefore either compute M from (106) with C given by (97), or by integrating (105) as Mðu; x; Þ ≔ Z x 0 M ;x ðu; x 0 Þdx 0 ; ð107Þ with Cðu; xÞ ≔ 2 Mð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 nontrivial test of the accuracy of our discretization.

No MOTS on coordinate light cones 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

where SG 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 backward 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 backward null surface cannot in general be a coordinate null cone with regular vertex.

We now ask if one can find a nonmarginally outertrapped surface (OTS) on u ¼ u 0 , one where ρ þ ≤ 0 everywhere. The geometric nongenericity 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 C): 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 [29,31]. 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 [25,27].

How then can we generalize this procedure in spherical symmetry to the nonspherical 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 maximize C½u 0 ; x 0 ðyÞ given by (97). 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. ( 27) 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 (97) 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. (85) 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 fu ¼ u Ã g ∩ Σ 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 ¼ AE1. 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.

IV. NUMERICAL METHODS

A. Legendre pseudospectral method in y 1. 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

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 D that if we represent them as

then in the linearization 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 backward and forward 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 ð0Þ il ≔ 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 ð1Þ il ≔ P 0 l ðy i Þ, and for f we use S ð2Þ il ≔ P 00 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.

Discrete versions of continuum identities and consistent truncation

We define the N × N matrices

where Y denotes the diagonal matrix

and we have defined the diagonal matrices

where

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 (D15), (D33) and (D38) between Legendre polynomials given in Appendix D, which in this notation can be written concisely as

We also have

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. ( 120) is a discrete version of the identity (D29), which we need for the linearized hierarchy equation for f, Eq. (D26), to be discretized exactly in y.

With the definition

we can derive the identity

Except for the presence of the factor I 0 on the left, (123) is the discrete version of the identity (D30), which is needed for the linearized hierarchy equation for b, Eq. (D25).

Without the right factor of I 0 , the left-hand side of (123) would have nonzero 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 linearized 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

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

Note that S with N ¼ 5 or N ¼ 3 is a 3 × 3 matrix, etc. These obey Āð0;2Þ Sð0;2Þ ¼ I;

To understand this, consider again the case N ¼ 3: we can represent P l for l ¼ 0, 2, 4, and P ð2Þ l for l ¼ 2, 4, 6, but P ð1Þ l 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,

We have then the same identities as already discussed, again with the proviso that P

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 pseudospectral method, we have explicitly evaluated the following matrices, which all vanish in the continuum, numerically in the F90 code (in double precision).

The matrices A, S, Δ and Λ are understood as either halfrange 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 -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 shorthand 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

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 0 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Þ .

B. 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 Þ ð148Þ in terms of the intermediate quantities

and the basic variables (metric coefficients and matter field)

These equations can be solved by integration as

int ðu; x; yÞ ¼ ϕ int ðu; 0; yÞ þ Z x 0 Sðu; x 0 ; yÞdx 0 ; ð151Þ for the intermediate quantities, and hence the constrained variables ϕ cons ≔ ðγ; b; ΞR; Ξf; ΞψÞ; ð152Þ in this order. The evolved variables ϕ evol ≔ ðf; R; ψÞ ð 153Þ

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 only through 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 labeled 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 shorthand

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

We then discretize the integrations (151) with the midpoint rule, that is

using the discretizations (156), (157), the discretizations of ϕ ;y and ϕ ;yy using D and D 2 , and the discretization of ϕ ;xy and ϕ ;xyy resulting from combining (156) 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 righthand gridpoint to respect causality.

In the formulation where we use D defined in (40) rather than ∂ x , we differentiate

and integrate

For the evolved variables (153), the definition of Ξ gives

The analogy of B with the x-component of a shift vector suggests that, in contrast to the x-derivatives in the righthand side of (148), ϕ evol;x in (161) 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 inner 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 secondorder accurate three-point formulas on grid points, namely

Where R ;x ðu; 0Þ is needed, we use R 0;j ¼ 0 in (162) to evaluate the right derivative as

and then average over y to obtain R ;x ðu; 0Þ.

D. The time step problem, and a resolution

An important, if somewhat vaguely defined, necessary condition for numerical stability of any explicit timeevolution 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 [31] 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 center, that is at r ∼ Δr. In Appendix A, 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 linearized wave equation into spherical harmonics as in (110) and finitedifference 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 [35]. 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 center, 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 time slices in [36]. Here the filtering uses fast Fourier transforms in θ and φ.

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

where we have defined the local time step criteria

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 jR ;x =ΞRj can be understood in two ways: as jR ;x j=jR ;u j 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 center this holds independently of l max .

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 center, 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 þ R x x l …dx 0 . We now obtain the values of those integration constants that correspond to a regular center (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

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 F 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 (171), and truncating this expansion as

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

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 singlevalued at the origin. Moreover, we show in Appendix F that regularity requires R lðlÞ ¼ 0 for all l, so that

Generalizing from the behavior of analytic solutions to the Einstein equations linearized about Minkowski spacetime (see Appendix D 5), and consistent with the regularity requirements in Appendix F we expand

Here b 1ð0Þ ðuÞ is free, corresponding to a gauge choice, see Appendix E. Finally, from

we expect that the expansions of Ξ-derivatives start at one power of x lower, that is ΞR ¼ X ∞ k¼0 X kþ1 l¼0 ðΞRÞ lðkÞ ðuÞx k P l ðyÞ; ð177Þ Ξf ¼ X ∞ k¼1 X kþ1 l¼2 ðΞfÞ lðkÞ ðuÞx k P 00 l ðyÞ; ð178Þ Ξψ ¼ X ∞ k¼0 X kþ1 l¼0 ðΞψÞ 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

Note that P 0 ðyÞ ¼ 1, P 1 ðyÞ ¼ y, P 2 ðyÞ ¼ ð3y 2 -1Þ=2, P 0 1 ðyÞ ¼ 1, P 0 2 ðyÞ ¼ 3y, P 00 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

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:

with all other components initialized to zero. The almostlinear 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 behavior of the error is unchanged.

A. Linearized solutions as test beds

In the small-data regime we can use exact solutions of the linearized field equations as test beds. Small data here means in practice that the difference between the solutions of the linearized 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 linearized 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 linearized 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. (B15)-(B16) in Appendix B.

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 D 3. This means that the linearized 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 linearized equations). By contrast, in sdn gauge δR develops dynamically even if it is set to zero in the initial data.

B. Convergence test method

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

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 (165) scales Δu in proportion to Δx. If the discretization in u is also at least second-order accurate, then we expect that the discretization 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 (165). 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

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

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 discretization error in x, or in u and x, at the reference resolution Δx ref and the fixed resolution 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

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 secondorder error, assumed to be leading. From pairs of numerical evolutions, we then calculate the quantity

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. 234, 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. 5, 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.

Generalized d'Alembert exact solutions

We have tested two kinds of exact solutions. The first are the generalized d'Alembert solutions for a single spherical harmonic derived in Appendix D and given in (D23) for the scalar field δψ l , and in (D48)-(D49) 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 test beds, previously in [9].

We choose the free function χ to be a Gaussian with center 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 factors of two. For these 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. (B14), 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 flavor of gsb is unstable in the ðG; xÞ formulation but stable in the ðγ; RÞ formulation (with and without integration by parts), and the lsB2 flavor 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 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.

dashed, and N x ¼ 1024 dot-dot-dashed. The middle plot in each column enlarges the lower part of the upper plot, and for clarity omits N y ¼ 17. The lower plot additionally omits N x ¼ 64, and is scaled with ðN x =128Þ 2 , to test for (local-in-u) second-order convergence with N x . The horizontal range is always 0 ≤ u ≤ 1.5, but the vertical range (rms error) has been chosen differently in each of the nine plots.

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. 2. However, we get pointwise second-order convergence at fixed x for all resolutions Δx < x=2.

Figure 3 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 .

Figure 4 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 Ny ¼ 17. We now set initial data directly for f l and ψ l , namely a Gaussian again with center 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 Ny ¼ 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 exact solutions

We show in Appendix D 6 that any solution of the scalar wave equation gives rise to a solution of the linearized Einstein equations, without the need for a spherical harmonic decomposition. In particular, the scalar wave equation admits plane-wave exact solutions, and in Appendix D 6 b we give these in Eq. (D61) for δψ and in (D62)-(D63) 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 test bed.

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

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 center 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, Ny ¼ 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 ¼ AE -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. 5.

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 color.

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 discretization 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 discretization 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 negligible, 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-discretization at early times and small N y and by x-discretization at late times and large N y . Plotting single-l components of the error against u and x further shows that the early discretization error in y arises mostly at x ≳ 1.5 and the late error discretization 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 ffiffi ffi 2 p , as there are now two identical waves of error instead of one in the same numerical domain. We therefore do not present plots here.

E. Tests of the computation of the Hawking mass

As a different indication of numerical error, we have implemented the expression for M given by ( 84) with (103), a centered finite differencing of this to compute the resulting M ;x , and, independently, the direct expression (105) 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 nontrivial test of the correctness of (84) with (103), (105), the hierarchy equations, and their discretizations. 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 center 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

The error with N y ¼ 5, 9, 17, 33, 65, 129 is the same, and so for Ny ¼ 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 nonvanishing 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 center 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.

A. 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 center, already identified in [9], 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 [27] 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 [30].

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.

C. Outlook

Two key questions remain open. The first is whether outgoing null cones emanating from a regular center 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 nonsphericity, rather than black hole formation. One reason to be optimistic is that in [5] we have constructed outgoing null cones emanating from a regular vertex in postprocessing 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 generalized affine parameter gauge. We leave this to future work.

The CFL condition ds 2 > 0 translates into

The choice s ¼ 0, q ¼ AE1 of points in the stencil, with the choice p ¼ 1 to locate the stencil next to the center, gives

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 (A3) 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

We let du ¼ -Δu, dx ¼ -sΔx, dθ ¼ qΔθ, and R ¼ pΔx.

At the center 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

With s ¼ 1, p ¼ 1, q ¼ AE1, and assuming Δθ ≪ 1, the equivalent of (A3) is now

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

Recall that G > 0. Choosing s ¼ AE1 with sign opposite to that of H we obtain

where B ≔ H=ð2GÞ as defined previously. This is just the CFL condition for the x-advection term in

APPENDIX B: 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

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

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

and obtain the metric

where

and hence

Therefore the general solution of (B1), (B2) with a regular center is (B7) and (B9): note this has only one free function UðuÞ. (B7) can be written as

In Bondi coordinates, where R ;u ¼ 0, we have

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 (B6) 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 (B3) and make the specific coordinate transformation

where x 0 > 0 is a parameter, to obtain (B6) with

It follows that

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

APPENDIX C: 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

The hierarchy equations become

where

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

and, with j∇Rj 2 ¼ -2ΞR=g, the Hawking mass in spherical symmetry is equal to the well-known Misner-Sharp mass

(Unlike the Hawking mass in general, the Misner-Sharp mass in spherical symmetry derives from a conserved stress-energy current [37].) The mass aspect (105) is

and the redshift defined in (77) is

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 center to infinity diverging. Assuming that the time step is limited by the Bondi shift (56), as in Eq. (165), we have

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

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 (C11) and (C1), we can write the resulting three equations as

and from this we can derive

From (C12) and (C17) we find

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

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 ¼ α.

APPENDIX D: AXISYMMETRIC LINEAR PERTURBATIONS OF MINKOWSKI SPACETIME 1. Linearized field equations in any radial gauge

As a test of our formulation and numerical methods beyond spherical symmetry, we linearize 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 linearized equations, as they appear in the linearization process.

The linearized hierarchy equations are

for δG, then

for δb and δf,

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

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

With the background solution G 0 ¼ R 0;x , and the boundary condition δG ¼ δR ;x at the origin, which follows from linearizing the gauge condition G ¼ R ;x that u is proper time at the origin, the solution of (D1) is

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

Linearized field equations in linearized 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 linearized 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 linearized Bondi gauge, defined by δR ¼ δG ¼ 0, δf and δb obey the coupled equations

These are just the first lines of (D2), (D3) above. The linearized wave equation (D5) does not simplify further. Using G 0 ¼ R 0;x , as well as δG ¼ δR ¼ 0, the linearized hierarchy equation (D4) becomes

Given a solution ðδf; δbÞ of (D7), (D8), (D9) 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 linearized Bondi gauge (D8), (D7) are given in Bondi background coordinates as Eqs. (D25), (D26) below.

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

Replacing the retarded time coordinate u with the Minkowski time coordinate

and writing δψðu; r; yÞ ≕ δϕðt; r; yÞ ¼ δϕðu þ r; r; yÞ; ðD12Þ this takes the more familiar form -δϕ ;tt þ δϕ ;rr þ 2 r δϕ ;r þ 1 r 2 ½ð1y 2 Þδϕ ;y ;y ¼ 0: ðD13Þ We initially work in coordinates ðt; r; yÞ and switch back to ðu; r; yÞ later. a. Separating off the y-dependence. Spherical symmetry of the background allows us to separate the angle y with the ansatz δϕðt; r; yÞ ≕ X ∞ l¼0 δϕ l ðt; rÞP l ðyÞ;

where P l ðyÞ is the Legendre polynomial of order l, obeying

We then have the one-dimensional wave equation

for the l-th partial wave, or equivalently We can find the general solution of (D16) that is regular at r ¼ 0 in terms of a single free function χ l of one variable, in the form of a generalized d'Alembert solution [38], as

where

and χ ðnÞ denotes the nth derivative of χ.

Expanding χðt AE rÞ into its Taylor series about r ¼ 0, we obtain the double sum the scalar field δψ governing the linearized 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.

d. The cases l = 0 and l = 1

Our separation of variables ansatz gives b 0 ¼ f 0 ¼ f 1 ¼ 0 because P 0 0 ¼ P 00 0 ¼ P 00 1 ¼ 0. In addition the potential ansatz (D48) gives δb 1 ¼ 0 because λ ¼ 0. However, as P 1 ¼ y and P 0 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 linearized equations in Bondi gauge without making the potential ansatz (D48), (D49). With δb ¼ δb 1 ðu; rÞ and δf ¼ 0, (D25) reduces to

while (D26) is obeyed trivially. The solution of (D52) that is finite at the center is constant in r. The resulting regular perturbation of H is given by (D9). Putting both together, we have

We show in Appendix E 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.

f. Consistent boundary conditions for the linearized Einstein equations decomposed into spherical harmonics

We now return to the case of a nontrivial r l ðuÞ > 0. From (D48) and (D49) we read off that δb l ðu; r l ðuÞÞ ¼ 0; ðD55Þ δf l ðu; r l ðuÞÞ ¼ 0; ðD56Þ δb l;r ðu; r l ðuÞÞ ¼ 2λ r l ðuÞ 2 δψ l ðu; r l ðuÞÞ; ðD57Þ δ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 linearized 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 linearized 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 linearized about Minkowski and split into spherical harmonics.

6. Solution without spherical harmonic decomposition a. 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 (D18), (D49) and (D50), we obtain the unseparated expression δfðu; r; yÞ ¼ δψ ;yy ðu; r; yÞδψ ;yy ðu; 0; yÞ þ 2 Z r 0 δψ ;yy ðu; r; yÞ r dr ðD59Þ for δf in terms of δψ. Substituting (D51) into (D25) and integrating twice, we obtain the expression δbðu; r; yÞ ¼ Z r 1 r4 Z r 0 r2 À 2ð1y 2 Þδf ;ry ðu; r; yÞ -8yδf ;r ðu; r; yÞ Á dr ðD60Þ for δb in terms of δf. b. Plane-wave solutions The scalar wave equation admits the plane-wave solution δψ ¼ χðt AE zÞ, that is δψðu; r; yÞ ¼ χðu AE Þ; u AE ≔ u þ rð1 AE yÞ; ðD61Þ for arbitrary functions χ. Substituting this into (D59), we can carry out the integration explicitly to obtain the plane wave solution δfðu; r; yÞ ¼ r 2 χ 00 ðu AE Þ þ 2 rχ 0 ðu AE Þ 1 AE y -2 χðu AE Þ -χðuÞ ð1 AE yÞ 2 ; ðD62Þ and substituting this into (D60), we can again carry out the integrations in closed form to obtain δbðu; r; yÞ ¼ AE2 rð1 ∓ yÞχ 00 ðu AE Þ -1 AE 3y 1 AE y ðχ 0 ðu AE Þχ 0 ðuÞÞ ! : ðD63Þ 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 AE yÞ. We can then write δψðu; r; yÞ ¼ χðu þ ϵÞ; ðD64Þ δfðu; r; yÞ ¼ r 2 χ 00 ðu þ ϵÞ þ 2r 2 d dϵ χðu þ ϵÞ -χðuÞ ϵ ; δbðu; r; yÞ ¼ AE2r ð1 ∓ yÞχ 00 ðu þ ϵÞ -ð1 AE 3yÞ χðu þ ϵÞ -χðuÞ ϵ ! ;

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

APPENDIX E: RESIDUAL GAUGE FREEDOM

A redefinition of the coordinates ðu; x; yÞ leaves the form (6) 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ð ũÞ; ðE1Þ y ¼ yð ũ; ỹÞ; ðE2Þ

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

μν ¼ 2∇ ðμ ξ νÞ ; ðE4Þ where ξ μ ≔ ðξ u ; ξ x ; ξ y ; 0Þ; ðE5Þ are functions of ðu; x; yÞ. The general solution of δg xx ¼ δg xy ¼ 0 is ξ u ¼ Aðu; yÞ; ðE6Þ ξ x ¼ ξ x ðu; x; yÞ; ðE7Þ ξ y ¼ Bðu; yÞ þ A ;y ðu; yÞS Z x 0 G e 2Sf R 2 dx 0 : ðE8Þ For A ¼ AðuÞ this reduces to the infinitesimal version of (E1)-(E3). For simplicity, we now restrict the background metric to flat spacetime in Bondi gauge, with R ¼ x, G ¼ H ¼ 1 and b ¼ f ¼ 0. (E8) then simplifies to ξ y ¼ Bðu; yÞ -SA ;y ðu; yÞ x :

The resulting metric perturbation is

xB ;y 2 -SA ;yy -2yA ;y 2 ; ðE12Þ δf ¼ -A ;yy 2x þ B ;y þ 2yB 2S ; ðE13Þ δb ¼ -A ;y x 2 þ B ;u S -xA ;uy þ ξ x ;y x 2 : ðE14Þ analysis, we relate the coordinates ðu; x; yÞ to the cylindrical Minkowski coordinates defined by t ≔ u þ x; z ≔ -yx; ρ ≔ ffiffiffi S p x;

with φ completing both coordinate systems. The inverse transformation is

and hence the coordinate basis vectors transform as

with ∂ φ in both coordinate systems. From (E6)-(E7) and (E19)-(E21), the residual gauge transformations expressed in this basis are

ξ 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

and substituting this back into (E22), we have

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

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

which explains why we have named the free coefficients T and Z. The corresponding pure gauge metric perturbations are δH ¼ 2ðT 0 ðuÞ þ zZ 00 ðuÞÞ; ðE27Þ δG ¼ T 0 ðuÞ; ðE28Þ δb ¼ -Z 00 ðuÞ; ðE29Þ 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

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 00 ðuÞxy; δb ¼ -Z 00 ðuÞ; ðE31Þ

with the other metric coefficients unchanged. Hence an infinitesimal noninertial motion of the origin along the symmetry axis gives exactly the l ¼ 1 regular metric perturbation (D53) that we already derived from the linearized 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 (30) or (43). Conversely, by imposing an arbitrary value of bðu; 0Þ, we could accelerate the origin during the evolution.

APPENDIX F: REGULARITY OF THE METRIC AT THE ORIGIN

Following common practice, Rinne [40] 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 φ; ðF1Þ 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 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

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 twistfree 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

or equivalently

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 2z 2 , a function is axisymmetric and regular if and only if it is an analytic function of ðt; x 2z 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 reparametrizing the Rinne coefficients as

we can show that the Rinne-regular twist-free axisymmetric metric is spherically symmetric if and only if

Moreover, because the coefficients of the redefinitions (F5)-(F7) themselves are Rinne-regular functions, and because we cannot absorb them into redefinitions of X 1;2;3 without making those singular, the reparametrized 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Þ ≔ thð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Þ ≔ xh o ðt; x 2 ; -xyÞ þ x 2 h e ðt; x 2 ; -xyÞ: ðF9Þ

For u to be a nontrivial 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 ¼ tx, 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.

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. IV E 2.