Simulations of gravitational collapse in null coordinates: III. Hyperbolicity
Carsten Gundlach









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







We investigate the well-posedness of the characteristic initial-boundary value problem for the Einstein equations in Bondi-like coordinates (including Bondi, double-null and affine). We propose a definition of strong hyperbolicity of a system of partial differential equations of any order, and show that the Einstein equations in Bondi-like coordinates in their second-order form used in numerical relativity do not meet it, in agreement with results of Giannakopoulos et al for specific first-order reductions. In the principal part, frozen coefficient approximation that one uses to examine hyperbolicity, we explicitly construct the general solution to identify the solutions that obstruct strong hyperbolicity. Independently, we present a first-order symmetric hyperbolic formulation of the Einstein equations in Bondi gauge, linearised about Schwarzschild, thus completing work by Frittelli. This establishes an energy norm (L^2 in the metric perturbations and selected first and second derivatives), in which the initial-boundary value problem, with initial data on an outgoing null cone and boundary data on a timelike cylinder or an ingoing null cone, is well-posed, thus verifying a conjecture by Giannakopoulos et al. Unfortunately, our method does not extend to the pure initial-value problem on a null cone with regular vertex. 





Simulations of gravitational collapse in null
  coordinates: III. Hyperbolicity
    Carsten Gundlach
    25 April 2024
===============================================================================







§ INTRODUCTION





Well-posedness of a system of partial differential equations (from now
on, PDEs) is defined as the existence and uniqueness of solutions and
their continuous dependence on the initial and boundary data.  For
formulations of the Einstein equations where surfaces of the time
coordinate u are outgoing null cones, three PDE problems, among
others, are of interest: 1) the double null initial value problem,
with free data posed on an outgoing null cone u=0 and ingoing null
cone v=0 that intersect in a spacelike 2-sphere; 2) the
initial-boundary value problem, with free initial data posed on u=0
and boundary data on a timelike cylinder r=r_0, again intersecting
in a spacelike 2-sphere; and 3) the initial value problem on a null
cone u=0 with regular vertex r=0.

The mathematical literature has focused on the proof of well-posedness
of such problems in some coordinate system and formulation that
is convenient for the proof. This might be called geometric
well-posedness. By contrast, in numerical relativity, well-posedness
of the continuum problem is necessary for the existence of a stable
discretisation, but other considerations are equally important for the
choice of formulation of the Einstein equations, and so one wants a
proof of well-posedness in the formulation of choice.

We give two examples of this difference in emphasis. A famous proof of
geometric well-posedness of the Cauchy problem was given by
Fourès-Bruhat in harmonic coordinates <cit.>, but
these became useful for numerical relativity only through the
breakthrough work of Pretorius <cit.>, which incorporated
two key modifications of lower-order terms: specific choices of the
“gauge source functions” of <cit.> that avoid coordinate
singularities, and the addition of “constraint damping” terms to
suppress violations of the harmonic and Einstein constraints arising
from numerical error <cit.>. In a second example,
much of the mathematical relativity literature uses double-null
coordinates on ingoing and outgoing null cones, but these are not
expected to be useful in numerical relativity beyond spherical
symmetry because the ingoing null cones form caustics. (It is one
purpose of the present series of papers to establish if outgoing null
cones also form caustics in specific strong-field applications.)

Geometric well-posedness for the double initial value problem was
established by Rendall <cit.> for smooth data, covering a
small neighbourhood of the intersection 2-sphere. The proof used
harmonic coordinates, made to coincide with u and v on the initial
surfaces u=0 and v=0.  This result was improved by Luk
<cit.> for a larger region (small finite distances to the
future of the two initial null surfaces) and rougher data. The proof
used double-null coordinates, and estimates in H^1 of a null tetrad
and the Ricci rotation coefficients and curvature components in this
tetrad.  The initial value problem on a regular null cone is even harder
because of the need to characterise data near the tip that will give a
regular solution. Existence was proved by Chruściel
<cit.>, see also <cit.>.

In the present series of papers, we focus on class of coordinates
(u,x,θ,φ) introduced in
<cit.>, and called “Bondi-like” in
<cit.>, where the surfaces of constant (retarded) time u are
outgoing null cones and the lines of constant (u,θ,φ) are
their generators. (To avoid confusion about terminology, we mention
already that “Bondi coordinates” are further defined by the radial
coordinate x being the area radius.)

In addition, we are interested in formulations of the Einstein
equations where a maximum number of equations can be solved as
ordinary differential equations (from now on, ODEs) along the
generators of the coordinate null cones, and only a minimum number,
namely two or three, contain u-derivatives: intuitively, the latter
are evolution equations for the two polarisations of gravitational waves,
plus, in double-null coordinates only, a third one for the area radius
R.

Such formulations are attractive for numerical relativity because they
are “maximally constrained”, meaning that one solves the same
constraint equations on each time slice as on the initial time slice
(so that only free data are evolved), and because these constraints
are not elliptic equations but inhomogeneous first-order linear ODEs
that can be solved by integration. (An exception to this statement are
affine null coordinates, where one of the constraints is a homogeneous
linear second-order ODE.) 

Time evolutions beyond spherical symmetry on null cones emanating from
a regular centre have been carried out, for example, for supernova
core collapse <cit.> and for scalar field critical
collapse <cit.>, and vacuum evolutions using
Cauchy-characteristic matching along a timelike cylinder, for example in
<cit.>.

However, we are not aware of any previous well-posedness result
specifically in Bondi-like coordinates. In an incomplete attempt,
Frittelli <cit.> constructed a first-order symmetric
hyperbolic form of the vacuum Einstein equations in Bondi coordinates,
linearised about Minkowski space. An “energy” estimate, in L^2 of
the reduction variables, then follows. However, the perturbation
Ṽ of the metric coefficient V was omitted from the
system. This is possible in the linearisation about Schwarzschild
because Ṽ couples to the other perturbations, but not vice versa.

The present paper is an attempt to reconcile the results of
<cit.>, and the suggestive
incomplete result of <cit.>, with recent work of
Giannakopoulos and collaborators. Their paper <cit.> found
that a first-order reduction of the null cone formulation is weakly,
but not strongly, hyperbolic (and hence not symmetric hyperbolic),
that the lack of strong hyperbolicity is essentially a gauge problem
in any Bondi-like gauge <cit.>, and that it appears to
indeed break the convergence of numerical solutions with resolution,
both in toy models and in an open-source code for the null
initial-boundary value problem <cit.>.

The structure of the paper is as follows. As in the previous papers in
this series <cit.>, in
Sec. <ref> we restrict to twist-free
axisymmetry, with a minimally coupled massless scalar field ψ as
matter. We briefly restate the metric, the mathematical structure of
the field equations and the gauge choices we consider.  We then
linearise the equations, first around an arbitrary background and then
around Minkowski spacetime, drop lower-order terms, and “freeze” the
background coefficients by treating them as constants. We now have a
homogeneous system of linear PDEs with constant coefficients that we
shall call the “toy model” of the original system. The symbol of the
toy model is also the principal symbol of the original system, and
strong hyperbolicity is an algebraic property of this principal
symbol. 

In contrast to <cit.>, we do not reduce
the PDEs to first order, but leave them in the form in which they are
solved numerically. We generalise the textbook definition of strong
hyperbolicity of first-order systems of PDEs to systems of PDEs of
arbitrary order. We show that, by this criterion, all known Bondi-like
null gauges are only weakly hyperbolic. We also find the general
solution of the toy model itself in closed form, and hence identify
the polynomial solutions that obstruct strong
hyperbolicity. (Appendix <ref> reminds the reader
of a textbook example of this phenomenon).

In Sec. <ref>, we use a completely different
approach. We relax the restriction to twist-free axisymmetry, but
linearise about the Schwarzschild solution, and restrict to Bondi
gauge. We present a first-order symmetric hyperbolic form of the full
linearised equations, thus completing the result of
<cit.>. (We correct some minor errors in
Appendix <ref>.) To include all 6 metric
coefficients in the estimate, we need to include also 10 first and 7
second derivatives of the metric as variables in the system and hence
in the estimate. These are far from all first and second derivatives,
and their choice is crucial for symmetric hyperbolicity. To connect
with Sec. <ref> with
Sec. <ref>, we give the equivalent symmetric
hyperbolic form of the toy model in
Appendix <ref>.

We then present well-posedness estimates for the initial value problem
on two intersecting null cones and the initial-boundary value problem
on a null cone intersecting a timelike world tube, closely following
<cit.>. Unfortunately, these methods do
not allow us to derive an estimate for the pure Cauchy problem on an
outgoing null cone with regular vertex, see
Appendix <ref> for the technical obstructions.

We summarize and conclude in Sec. <ref>.





§ OBSTRUCTIONS TO STRONG HYPERBOLICITY OF THE SECOND-ORDER FORM
OF THE EINSTEIN EQUATIONS IN BONDI-LIKE GAUGES







 §.§ Metric and field equations in twist-free axisymmetry




We begin with a brief review of our setup in this ection, see
<cit.> for full details. We can write the metric
of any twist-free axisymmetric spacetime in the form

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

We use the angular coordinate y:=-cosθ, so that the range
0≤θ≤π corresponds to -1≤ y≤ 1, and the shorthand
S:=1-y^2. The other angular coordinate φ has the usual range
0≤φ<2π. The Killing vector generating the axisymmetry is
∂_φ. (We use the convention of equating vector fields
with derivative operators.) (G,H,R,b,f) and the scalar field ψ
depend on (u,x,y). In <cit.> we assumed that
R=0 occurs at x=0 and is a regular centre, but we do not make this
assumption here.
 
Each surface N^+_u of constant u is an outgoing null cone,
and is ruled by the null geodesics L^+_u,y,φ, which
are also coordinate lines. Each surface S_u,x of constant
u and x is assumed to be spacelike, and has topology S^2 and
area 4π R^2. The
outgoing future-directed null vector field normal to S_u,x
is U:=G^-1∂_x, and is also tangent to the affinely
parameterised generators of N^+_u.  The ingoing null normal
on S_u,x is

    Ξ=∂_u -H 2G∂_x-bS ∂_y.

It is normalised to Ξ^aU_a=-1.

The field equations we want to solve are the Einstein equations

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

and the massless, minimally coupled wave equation

    ∇^a∇_aψ=0.

(We use units where c=G=1.) A subset of the Einstein equations, plus
the wave equation, take the form

    (lnG R_,x)_,x   =   S_G[R,f,ψ], 
    (R^4e^2S fb_,x G)_,x   =   S_b[R,f,ψ,G], 
    (R Ξ R)_,x   =   S_R[R,f,ψ,G,b], 
    (R Ξ f)_,x   =   S_f[R,f,ψ,G,b]-(Ξ R)f_,x,
    (R Ξψ)_,x   =   S_ψ[R,f,ψ,G,b]-(Ξ R)ψ_,x,

where Ξ is the derivative operator defined in (<ref>). We
call these the “hierarchy equations”. H and ∂_u appear
only in the combination Ξ. The right-hand sides S[f,…] are
given in full in <cit.>. They contain the
derivatives f_,x, f_,y, f_,xy and f_,yy (but not
f_,xx), and similarly for R, G, b and ψ, with
the exception that ψ_,xy and b_,yy do not appear. 





 §.§ Gauges and solution algorithm




With H given (for example H=0 in double-null gauge), the hierarchy
equations take the form of first-order linear ODEs in x for b,
G, Ξ R, Ξ f, and Ξψ that can be solved by
integration. With R given (for example R=x in Bondi gauge),
(<ref>) is solved for H, given R_,u.  In a third group of
gauges, where G is given (for example G=1 in affine gauge),
(<ref>) becomes a second-order linear ODE in x for R that is
solved first, the other equations are again solved by integration for
b, Ξ R and Ξ f, and (<ref>)_,x and
(<ref>)_,u are combined to find an equation that can be
solved for H by integration.

With f_,u, R_,u and ψ_,u now known, f, ψ, and in
double-null gauge also R, are now advanced in u, and the hierarchy
equations are then solved again. Note that the algorithm is
“maximally constrained” in the sense that the hypersurface equations
are solved at each time step, as they are from the free initial data.

Consider now the characteristic initial-boundary value problem, with
outgoing null boundary (u=0,x>0) and timelike or null inner boundary
(x=0,u>0). (Geometrically, one can think of the problem on two
intersecting null cones as a pure initial-value problem, but in
Bondi-like coordinates it is more natural to think of data on u=0 as
initial data and data on v=0 as boundary data, as we solve the
constraints by integration in v.)

In double null gauge we specify R, f and ψ on
u=0, and Ξ R, Ξ f, Ξψ, b, b_,x and G on
x=0, or nine functions of two variables. In Bondi gauge we specify
f and ψ on u=0, and Ξ f, Ξψ, b, b_,x, G
and H on x=0, or eight functions of two variables. Finally, in
affine gauge, we specify f and ψ on u=0, and R, R_,x,
Ξ f, Ξψ, b, b_,x, H and H_,x on x=0, or ten
functions of two variables.

In each case, the data on u=0 can be specified freely, while the
data on x=0 are constrained by some of the remaining Einstein
equations, which do not contain x-derivatives. We do not discuss
these constraints here, but evaluate the hyperbolicity of the
evolution equations with the constraints relaxed.





 §.§ Definition of strong hyperbolicity




As already noted, a necessary condition for well-posedness of a
nonlinear PDE is well-posedness of its linearisation, in the principal
part, frozen coefficient approximation.  We denote the field equations
linearised around a background solution ϕ by 

    L( x,ϕ,∇)δϕ=0,

We denote the principal part of L by L^p. The principal symbol is
the matrix-valued function L^p( x,ϕ,i k) of the wave number
covector k.  For a quasilinear system (such as the Einstein
equations), the principal symbol of the linearisation is the same as
the principal part of the full equations.

We can use the principal symbol to find plane-wave solutions to the
linearised field equations in the frozen coefficient approximation:
substituting δϕ( x)=e^i k· xδϕ( k)
into L(∇)δϕ( x)=0 gives
L(i k)δϕ( k)=0, which is now a system of linear
algebraic equations. Nontrivial solutions exist only for k such
that L(i k)=0. These k are called characteristic covectors,
and vectors δϕ( k) in the corresponding null space
of L(i k) are called characteristic variables. The linearised
problem with frozen coefficients is well-posed in L^2 if the
plane-wave solutions with real k are complete in the sense that we
can map them smoothly and one-to-one to the boundary and initial data.

Consider now the Cauchy problem for a system of linear PDEs with
constant coefficients. Initial data are imposed on a spacelike surface
t=0, with normal covector n:=dt, and we assume that the initial
data fall off as | x|→ 0, so that we can Fourier-transform the
initial data and the solution in x. Then a sufficient criterion for
well-posedness is given by strong hyperbolicity. We offer the
following definitions:

A system of PDE is defined to be strictly hyperbolic in the time
direction n, in a neighbourhood of the solution ϕ, if, first,
L^p(ϕ,i n) 0, and second, all roots ω of L^p(ϕ,iω n+ik̅)=0 are real and distinct, for all real k̅
that are not zero or proportional to n (see, for example,
<cit.>). (We shall use the notation k for arbitrary
covectors, k̅ for what intuitively are vectors in “space”, and
k_ω:=ω n+k̅.)

If the roots are not all distinct, then the system is strongly
  hyperbolic if the null spaces of L^p(ϕ,iω
  n+ik̅) have dimension corresponding to the multiplicity
of each root ω, and these null spaces depend smoothly on
ω. This generalises the textbook definition for first-order
systems (see, for example, <cit.>).

If the roots are all real, but the system is not strongly hyperbolic,
it is called weakly hyperbolic.

Strict hyperbolicity implies strong hyperbolicity, and for the
linearised equations with frozen coefficients both mean that after
Fourier-transforming in x, initial data on a surface with
normal covector n can be decomposed into plane waves
exp(i k· x)δϕ( k), with δϕ(
  k) in the null spaces mentioned above, each of which propagates
with a speed ω in the direction - k. This in turn
implies local well-posedness of the pure initial-value problem.

If we reduce a strongly hyperbolic system of arbitrary higher order to
first order by introducing suitable derivatives of the original
variables as reduction variables, the general solution translated into
these variables still consists of purely oscillating plane
waves. Hence any such first-order reduction is strongly
hyperbolic. However, the converse is not true: a higher order system
that is only weakly hyperbolic according to our definition may or may
not admit a strongly hyperbolic first-order reduction.





 §.§ The principal symbol





To fix notation, we write our variables in the order
ϕ^†:=(G,b,R,f,H,ψ), the coordinates in the order
x:=(u,x,y), so that ∇:=(∂_u,∂_x,∂_y),
and we use the notation k:=(μ,ξ,η) for the wave number
covector.

For our field equations, the principal terms clearly include the
second derivatives of G, b, R, f and ψ. The only
derivative of H that appears is H_,x, and this should therefore
be considered principal. Finally, because we want to think of
(<ref>) as determining G, we include the first derivative
G_,x as a principal term in this equation, but not in the other
equations, which also contain second derivatives of G. The equations
are quasilinear in this sense. 

Note that until we have fixed a specific null gauge, the principal
symbol will have five rows, corresponding to the linearisations of the
field equations (<ref>-<ref>), but six columns,
corresponding to six metric perturbations. The preliminary principal
symbol and the corresponding vector of perturbation variables are


    L^p( x,ϕ,i k)=(
    [            i Gξ               0       1 R_,xξ^2               0               0               0;       -R^2 Gη ξ -e^2S fR^4 Gξ^2        -
    2R η ξ       2R^2S η ξ               0               0;         R^2  4Y       R^2S 4η ξ        X+GR  2Y      -GR^2S  2Y     -iRR_,x 2Gξ               0;           R 4SY          R 4η ξ               0               X    -iR f_,x 2Gξ               0;               0               0               0               0    -iR ψ_,x 2Gξ        X+GR  2Y;                 ]),
        δϕ( x):=([ δ G; δ b; δ R; δ f; δ H;  δψ ]).



We have defined the shorthands

    X   :=   R ξ(-μ+b S  η+H 2G ξ), 
    
    Y   :=   e^-2S fS R^2η^2.

X is the principal symbol of the derivative operator
∂_xRΞ that appears on the left-hand side of each of the
hierarchy equations.  Y is the symbol of the principal angular
derivative -g^yy∂_y∂_y.

Each of the three classes of Bondi-like gauges that we have already
discussed simply eliminates one of the columns of (<ref>), giving
us a square 5× 5 symbol: this is the δ H column in
generalised double-null gauges, where H is given and δ H=0;
the δ R column in Bondi and related gauges, where R is given
and δ R=0; and the δ G column in affine and related
gauges, where G is given and δ G=0. We do not
write the 5× 5 principal symbols for these gauges out
explicitly, as they can be read off trivially from (<ref>).

In any null gauge, the “time” direction n=du=(1,0,0), that is, using
u as the time coordinate, fails the first criterion for strong
hyperbolicity. That is of course expected, as the surfaces of constant
u are characteristic. Instead we follow <cit.> and analyse
hyperbolicity on spacelike time slices with normal covector

    n:=(1,1 A,0),

with A>0. For now we leave A and k̅=(μ,ξ,η) general, so
that

    k_ω=(ω+μ,ω A+ξ,η).






 §.§ Double-null and related gauges





In double-null gauge H=0 and δ H=0. In related gauges
<cit.>, δ H is given in terms of δ
G, and so the resulting δ H_,x is not principal in the
equations where it appears, because these contain second derivatives
of G. In either case we can delete the δ H column in the
temporary principal symbol to obtain the 5× 5 principal symbol
for these gauges, and we remover δ H from δϕ. 

The first row of L^p(ϕ,i k) now contains first and second order
in k terms, and the other rows are homogeneous in k of second
order. However, the cofactor of the entry ξ^2/R_,x in the
first row, third column, of L^p(ϕ,i k) is zero, and so L^p(ϕ,i k) is homogeneous of order nine, rather than ten, in k.

The nine roots ω of L^p(ϕ,iω n+i k)=0 are

    ω_0:=-Aξ

with multiplicity four,

    ω_±:=ω_c±Δω,

each with multiplicity two, and

    ω_2:=2ω_c+Aξ=-2Zμ̅-Aξ,

with multiplicity one. Here
we have defined the shorthands

    ω_c   :=   -Zμ̅-Aξ, 
    Δω   :=   √(Z^2μ̅^2
    +AGZ Y), 
    μ̅   :=   μ-Aξ-bSη, 
    
    Z   :=   AG 2AG-H.

The roots ω_± with multiplicity two have corresponding null
spaces of dimension two, and the null space of
ω_2 has dimension one, but the null space of ω_0 has
only dimension two, not four,  so two fewer than the multiplicity of the root.

Attention must also be given to special directions in which roots
ω merge. In the special direction given by μ̅=0, the
roots ω_2=ω_0 have merged and ω_0 is now a
quintuple root, but the corresponding null space is still only
two-dimensional, so (exceptionally) three fewer than the multiplicity
of the root.

The characteristic covectors k_ω themselves are independent of
the parameter A of the time slicing, as long as A>0, and to find
all of them it is sufficiently general to give k̅ only two
algebraically independent components. In particular, k_± and
k_2 can be written as

    k_±   =   (GY 2ξ̃+H
      2Gξ̃+bSη,ξ̃,η), 
     k_2   =   (H
      2Gξ̃+bSη,ξ̃,η),

where now ξ̃ and η are arbitrary, except that
ξ̃>0 for k_+ and ξ̃<0 for k_-. Similarly,
k_0 can be written as

    k_0=(μ̃,0,η),

where now μ̃ and η are arbitrary. We see that
k_0 parameterises plane waves that do not depend on x. The
k_± are null covectors, while k_0 and k_2 are
spacelike. They obey U· k_0=0 and Ξ· k_2=0.






 §.§ Bondi and related gauges




In Bondi and related gauges, R is given and δ R=0, and so we
eliminate the column corresponding to it. We now have eight roots
ω: ω_0 with multiplicity four and null space of
dimension two, so two fewer than the multiplicity of the root, and
ω_± with multiplicity two and the same null spaces of
dimension two as in the double null case. There is no problem of
merging roots.





 §.§ Affine and related gauges





In affine gauge, G is given and δ G=0, and so we eliminate
the column corresponding to it. With G=G(u),
  (<ref>) and (<ref>) take the form

    R_,xx   =   (...)R, 
    
    (RR_,u)_,x+1 2G(RR_,xH)_,x   =   (...),

where the dots stand for terms already known at the point where the
equations are solved. (<ref>) is then solved as a
second-order ODE in x for R.  Starting from
(<ref>)_,x, we use (<ref>)_,u to
eliminate R_,uxx and (<ref>) to eliminate R_,ux,
and obtain an equation of the form H_,xx=(...). 

To reflect this, in affine gauge we multiply the third row of
(<ref>), which comes from (<ref>), by iξ. This will
obviously multiply L_p(i k) by a factor of iξ, and hence
will gives us an extra root ω_0, without changing the dimension
of the corresponding null space.

We now have ten roots ω: ω_0 with multiplicity six and
null space of dimension three, so three fewer than the multiplicity of
the root and ω_± with multiplicity two and the same null
spaces of dimension two as in the double null case. There is no
problem of merging roots.





 §.§ Linearisation about Minkowski spacetime and frozen
  coefficient approximation




The principal symbol, and hence the characteristic covectors and
variables are simpler in the linearisation about Minkowski, but the
multiplicities of the roots and dimensions of the corresponding null
spaces are unchanged, so we now restrict to this case.

Linearising about Minkowski spacetime, so that b=f=ψ=0, G=H=1
and R=x in the background, we have the 5× 6 principal symbol


    L^p_Mink( x,i k)=(
    [                 iξ                  0                ξ^2                  0                  0                  0;             -r^2ηξ            -r^4ξ^2              -2rηξ            2r^2Sηξ                  0                  0;             S 4η^2           r^2S 4ηξ S 2rη^2-rμξ+r 2ξ^2         -S^2  2η^2             -ir 2ξ                  0;            1 4rη^2              r 4ηξ                  0         -rμξ+r2ξ^2                  0                  0;                  0                  0                  0                  0                  0 S 2rη^2-rμξ+r 2ξ^2;                    ]).

In the frozen coefficient approximation, we denote the frozen value of
R in the background by r, but the radial argument of the linear
perturbations by x. 

We can simplify (<ref>) further, as follows. We define the
coordinate y̅ and corresponding wave number k as

    k:=√(S) rη,   y̅:=r√(S)y,
      ⇒  
    ky̅=η y,   ∂_y=r√(S)∂_y̅.

Note that in the frozen coefficient approximation both y and y̅ are considered to have infinite range, while S appears only as
a frozen coefficient, ignoring that S=1-y^2=sin^2θ. Similarly,
r stands for a frozen coefficient, while x is considered to have
infinite range. We multiply the rows of (<ref>) by the
constants (1,√(S),4,4S,2) respectively, and renormalise the
perturbations by constants as

    δϕ^†:=(δ G,(r√(S))δ b,(2/r) δ
      R,(2 S)δ f,(2/r)δ H,δψ) 
    :=(δ G,δ b,
    δ R,δ f,δ H,
    δψ),

which corresponds to multiplying the columns of (<ref>) by
constants. With this notation, the preliminary 5× 6 principal
symbol of the linearisation about Minkowski space takes the form

    L^p_Mink,frozen(i k)=(
    [          iξ           0      r 2ξ^2           0           0           0;         -kξ        -ξ^2        - kξ          kξ           0           0;         k^2          kξ k^2-2μξ+ξ^2        -k^2         -iξ           0;         k^2          kξ           0    -2μξ+ξ^2           0           0;           0           0           0           0           0 k^2-2μξ+ξ^2;             ]).


We have been able to absorb all but one of the background-dependent
coefficients into our redefinitions.

We also simplify the calculation of characteristic covectors by
restricting our ansatz to n:=(1,1,0) and k̅=(0,ξ,k), so that

    k_ω=(ω,ω+ξ,k).

In all three classes of Bondi-like gauges, L^p(i k_ω)=0 has roots

    ω_±=±√(ξ^2+k^2),

corresponding to characteristic covectors

    k_±=(±√(ξ^2+k^2),ξ±√(ξ^2+k^2),k),

with multiplicity two and null spaces of
dimension two. In addition, in double null gauge we have
ω_0=-ξ, or

    k_0=(-ξ,0,k),

with multiplicity four but null space of dimension two,
and ω_2=ξ, or

    k_2=(ξ,2ξ,k),
 with multiplicity one.
In Bondi gauge we have ω_0 with multiplicity four but null
space of dimension two, and in affine gauge ω_0 with multiplicity
six but null space of dimension four. So in each gauge the
obstruction to strong hyperbolicity is that the multiplicity of the
characteristic vector k_0 is two more than the dimension of the
corresponding null space.

In the special case ξ=0, we have

    k_±=(± k,± k,k),      k_2= k_0=(0,0,k),

and in the special case k=0, we have

    k_+= k_2=(ξ,2ξ,0),      k_-= k_0=(-ξ,0,0).

For ξ=k=0, there are no characteristic covectors.





 §.§ A toy model for the obstruction to strong hyperbolicity








If a weakly hyperbolic system of homogeneous linear PDEs
with constant coefficients fails to be strongly hyperbolic because
there are not enough plane-wave solutions, the missing solutions must
be a polynomial times a plane wave. (The reader is reminded of a
well-known textbook example in Appendix <ref>.)
 
We now use the principal part, frozen coefficient approximation as a
toy model for the full linearised equations. The fact that it has
constant coefficients allows us to find the general solution in closed
form, and in particular the polynomial solutions that obstruct strong
hyperbolicity.

In the linearisation about Minkowski spacetime the scalar wave
equation decouples from the metric perturbation in any choice of
Bondi-like null gauge, and so we can consider it separately. In the
frozen coefficient approximation, it takes the form

    -2δψ_,ux+δψ_,xx+δψ_,y̅y̅=0.

Note this is the actual wave equation on 2+1-dimensional Minkowski
spacetime with metric

    ds^2=-2 du dx-du^2+dx^2+dy̅^2.

A plane-wave ansatz gives the solution

    δψ(u,x,y̅) =∑_±∫∫δψ_±(ξ,k) e^i k_±· x dξ dk.

The functions δψ_±(ξ,k) map one-to-one to, for
example, characteristic data in L^2 on null cones u=0 and
v:=u+2x=0, or to Cauchy data in L^2 on t:=u+x=0.  Hence this
solution is complete by function counting, with no polynomial
solutions required.

We now derive the general solution for the metric perturbations. We
restrict to vacuum, δψ=0, without loss of generality.
As the missing plane-wave solutions are for k_0=(-ξ,0,k), which
parameterises functions that are independent of x, we expect the
missing solutions to be polynomial in x. The physical gravitational
waves can be expressed in terms of derivatives of a solution Ψ of
the scalar wave equation <cit.>, which can be parameterised as

    Ψ(u,x,y̅) =∑_±∫∫ĉ_±(ξ,k) e^i k_±· x dξ dk,

in complete parallel to (<ref>). In double-null gauge,
we also introduce the shorthand

    Ψ_2(u,x,y̅) = ∫∫ĉ_2(ξ,k) e^i k_2· x dξ dk.

for the general solution in L^2 of the advection equation

    2Ψ_2,u-Ψ_2,x=0

along ingoing null cones.

In double null gauge and its generalisations, we set
δ H=0. The equations in the principal part, frozen coefficient
approximation can be read off from (<ref>) with the
δ H column deleted. They are

    δ G_,x-r 2 δ R_,xx   =   0, 
    δ G_,xy̅+δ b_,xx
    +δ R_,xy̅
    -δ f_,xy̅   =   0, 
    
    -δ G_,y̅y̅
    -δ b_,xy̅
    +2δ R_,ux-δ R_,xx
    -δ R_,y̅y̅
     
    +δ f_,y̅y̅   =   0, 
    
    -δ G_,y̅y̅-δ b_,xy̅
    +2δ f_,ux -δ f_,xx   =   0.

The general solution is

    δ G     =   δ G_0,u + r 2Ψ_2,x,   
    δ b   =   δ b_0 
    + δ b_1,u x -  r 2Ψ_2,y̅ +Ψ_,y̅, 
    δ R   =   δ f_0 
    + 1 2(δ b_1,y̅ + δ G_0,y̅y̅)x 
    + Ψ_2, 
    δ f   =   δ f_0 
    + 1 2(δ b_1,y̅ + δ G_0,y̅y̅)x 
    +Ψ_2  +Ψ_,x.

The five free functions (of two variables) ĉ_±, ĉ_2,
δ b_0, δ f_0 parameterise plane
waves, while the free functions δ G_0 and δ
  b_1 parameterise linear-in-x solutions. These seven free
functions correspond to seven characteristic covectors: k_± and
k_2, all with multiplicity one, and k_0 with multiplicity
four. They also correspond to the freedom to set five functions
(b,b_,x,G,Ξ R,Ξ f) of (u,y̅) at x=0 and two functions
(R,f) of (x,y̅) at u=0. Hence our solution is complete by
function counting.

The linear-in-x solutions are problematic not because they grow (and
so are not in L^2), but because they grow arbitrarily rapidly in x
for boundary data at x=0 that oscillate arbitrarily rapidly in y.

In Bondi gauge and its generalisations, we set δ
  R=0. The equations in the principal part, frozen coefficient
approximation can be read off from (<ref>) with the
δ R column deleted. They are

    δ G_,x   =   0, 
    δ G_,xy̅+δ b_,xx-δ
      f_,xy̅   =   0, 
    
    -δ G_,y̅y̅-δ
      b_,xy̅+δ f_,y̅y̅-δ H_,x   =   0, 
    
    -δ G_,y̅y̅-δ b_,xy̅+2δ
      f_,ux -δ
      f_,xx   =   0.

The general solution is

    δ G     =   δ G_0,   
    δ b   =   δ b_0 
    -δ G_0,y̅ x + Ψ_,y̅, 
    δ f   =   δ f_0  
    +Ψ_,x, 
    δ H   =   δ H_0 
    +δ f_0,y̅y̅ x.

The four free functions ĉ_±, δ b_0 and
δ H_0 parameterise plane waves, while
δ G_0 and δ f_0 parameterise
linear-in-x solutions. These six free functions correspond to six
characteristic covectors: k_± with multiplicity one and k_0
with multiplicity four. They also correspond to the freedom to set
five functions (b,b_,x,G,Ξ f,H) at x=0 and one function f at
u=0. Hence our solution is complete by function counting.  

Finally, in affine gauge and its generalisations, we set
δ G=0. The equations in the principal part, frozen coefficient
approximation can be read off from (<ref>) with the
δ G column deleted. They are

    δ R_,xx   =   0, 
    δ b_,xx
    +δ R_,xy̅
    -δ f_,xy̅   =   0, 
    
    -δ b_,xy̅
    +2δ R_,ux-δ R_,xx
    -δ R_,y̅y̅
    
    +δ f_,y̅y̅-δ H_,x   =   0, 
    
    -δ b_,xy̅
    +2δ f_,ux -δ f_,xx   =   0.

Following what is done in the full nonlinear equations, we take an
x-derivative of (<ref>), and use derivatives of the
other equations to simplify it. In the principal part, frozen coefficient
approximation the result is simply

    δ H_,xx=0.


The general solution of (<ref>-<ref>,<ref>)
is

    δ b   =   δ b_0 
    +2δ R_1,u x + Ψ_,y̅, 
    δ R   =   δ R_0 + δ R_1,y̅ x,  
    δ f   =   δ f_0 + δ R_1,y̅ x
    +Ψ_,x, 
    δ H   =   δ H_0 
    +δ H_1 x.

The six free functions ĉ_±, δ b_0,
δ H_0, δ R_0 and δ
  f_0 parameterise plane waves, while δ R_1 and
δ H_1 parameterise linear-in-x solutions. These
eight free constants correspond to eight characteristic covectors:
k_± with multiplicity one and k_0 with multiplicity six. They
also correspond to the freedom to set seven functions
(b,b_,x,G,R,R_,x,Ξ f,H,H_,x) at x=0 and one function f
at u=0. Hence our solution is once again complete by function
counting.

If we solved (<ref>) instead of (<ref>) then
δ H_1 would not be free, but given by
δ H_1=δ f_0,y̅y̅-δ
  R_0,y̅y̅. However, this is a feature of the toy
model only, where are able to solve in closed form for
δ R.

Table <ref> gives an overview of free null data
and normal data, the characteristic covectors, and the free functions
in our explicit solution (in the frozen-coefficient approximation, in
the linearisation about Minkowski). In each case, if we add
δψ back in, the number of free function, matching
the free data, becomes 9, 8 and 10 respectively, as discussed above in
Sec. <ref>-<ref>. In all gauges, the
boundary data at x=0 are subject to constraints, which we do not
impose here.





§ SYMMETRIC HYPERBOLIC FIRST-ORDER
  REDUCTION OF THE LINEARISED EINSTEIN EQUATIONS IN BONDI GAUGE









 §.§ Balance laws from symmetric hyperbolic first-order systems





This subsection reviews relevant parts of <cit.> for
completeness and to establish notation. A system of linear first-order PDEs

    C^μϕ_,μ=Dϕ,

where ϕ∈ R^N is a vector of dependent variables, C^μ and
D are real N× N matrices that depend smoothly on
x^μ∈ R^n, and the C^μ are symmetric, implies the balance law

    j^μ_,μ   =   S, 
    
    j^μ   :=   ϕ^† C^μϕ, 
    
    S   :=   ϕ^† Dϕ, 
     D   :=   D+D^†+C^μ_,μ.

The matrix D is real and symmetric, and so can be
diagonalised with real eigenvalues.

If furthermore C^μ k_μ depends smoothly on the covector
k_μ 0, Gauss' law gives

    ∫_∂ Vj^μ dΣ_μ=∫_V S dV

for any “control” spacetime volume V with boundary ∂ V.

If furthermore there exists a time coordinate t (a smooth function
with (dt)_μ everywhere future pointing), such that

    j^t:=j^μ (dt)_μ

is positive definite in ϕ, the system is called symmetric
hyperbolic with respect to t.
 
We now consider two elementary generalisations. If we redefine the
dependent variables as

    ϕ:=Aϕ̃,

where A is invertible and depends smoothly on the
x^μ, (<ref>) becomes

    C̃^μϕ̃_,μ=D̃ϕ̃,

where

    C̃^μ   :=   A^† C^μ A, 
    D̃   :=   A^†(DA-C^μ A_,μ).

The resulting balance law is the same as before, namely

    j̃^μ   :=   ϕ̃^†C̃^μϕ̃=j^μ, 
    S̃   :=   ϕ̃^†D̃ϕ̃=S, 
    D̃   :=   A^† DA,

In particular, the A_,μ term in D̃ in (<ref>)
cancels out of D̃ in (<ref>). However, we
can change the eigenvalues of D̃ and C̃^μ by
choosing an invertible but non-orthogonal A. To see this, note that
because D is symmetric, we have

    D=RΛ R^†,      R^† R=I,

so

    D̃=A^† RΛ R^†A,

Hence D̃ has the same eigenvalues as D if and
only if A is orthogonal.  The same is true for the matrices
the C̃^μ. 

As a second modification, we consider a change of integration weight
dV in (<ref>). The balance law
(<ref>) is equivalent to

    (ωj^μ)_,μ=ω S_ω

with

    S_ω:=S+ω_,μωj^μ
    =ϕ^†( D+ω_,μωC^μ)ϕ,

and so we have

    ∫_∂ Vj^μ ω dΣ_μ=∫_V S_ω ω dV

for any smooth function ω>0.






 §.§ Estimate on an arbitrary control volume





We now review how a symmetric hyperbolic system of
homogeneous first-order PDEs gives rises to “energy” estimates that
demonstrate well-posedness in the corresponding “energy norm”.  We
follow the basic idea of <cit.>, but give an
alternative derivation.

Let n_μ 0 denote an outward-pointing covector field on
∂ V. As the C^μ are real and symmetric, at any point in
∂ V the matrix C^μ n_μ can be diagonalised with real
eigenvalues, and so the space R^N of dependent variables ϕ can
be written as the sum of the positive, negative and zero eigenspaces
of C^μ n_μ, or R^N=V_+⊕ V_-⊕ V_0. Hence the
outward-pointing flux j^μ n_μ=ϕ^† C^μ n_μϕ can
be written as a term that is positive definite on V_+ plus one that
is negative definite on V_-. Integrating over ∂ V, we can
then write (<ref>) schematically as

    ∫_V S dV=:||out||^2-||in||^2.

In applications, we will split ||out||^2 into

    ||out||^2=||out'||^2+||out”||^2,

where ||out”||^2 denotes any part of the outgoing flux that
we do not want to include in our estimate. We trivially obtain

    ||out'||^2≤ ||in||^2+∫_V S dV.


Now assume there is a slicing of V by hypersurfaces of constant t
such that the system is symmetric hyperbolic with respect to t. Let
c>0 be the smallest eigenvalue of C^t:=C^μ(dt)_μ anywhere in
V, and let d be the largest positive eigenvalue of D
anywhere in V, or zero if D is negative definite everywhere
in V. We have S≤ dϕ^†ϕ and j^t≥
cϕ^†ϕ, and so we can bound the source term S of the
balance law as

    S≤d c j^t.


We slice V into surfaces of constant t,

    ∫_V S dV=∫_t_0^t_1 S_t  dt,     
      S_t:=∫_V∩Σ_t S  dΣ_t,

where t_1:=sup_V t and t_0:=inf_V t, and use the bound
(<ref>) to obtain

    S_t≤d c E_t,     
      E_t:=∫_V∩Σ_t j^t  dΣ_t.

We now evaluate (<ref>) on the control volume
V_t:=V∩{t'<t} (the part of V to the past of t'=t) to obtain

    E_t   =   ||in||^2_(∂ V)_t
    -||out||^2_(∂ V)_t
    +∫_t_0^t S_t' dt' 
    
       ≤   ||in||^2+d c∫_t_0^t  E_t' dt',

where we have defined (∂ V)_t:=∂ V∩{t'<t}⊆∂ V. In (<ref>) we have used
||in||^2_∂ V_t≤ ||in||^2, which follows
from ∂ V_t⊆∂ V.

We now differentiate to turn the integral inequality
(<ref>) into a
differential one,

    d dt E_t≤d c E_t,      E_t_0≤ ||in||^2,

and solve this to obtain

    E_t≤ e^d c(t-t_0) ||in||^2.

We then have

    ∫_V S dV≤∫_t_0^t_1d c E_t  dt
    ≤(e^d c(t_1-t_0)-1) ||in||^2 ,

and substituting (<ref>) into (<ref>) we obtain the
desired estimate

    ||out'||^2≤  e^d c(t_1-t_0)||in||^2,

as derived in <cit.>, but here for arbitrary t_1-t_0 and V.





 §.§ The scalar wave equation on Schwarzschild





In spherical polar coordinates, symmetry of the matrices C^μ is
less obvious than it is in Cartesian coordinates . This is best
illustrated if we consider the scalar wave equation on the
Schwarzschild background, which decouples to linear order from the
metric perturbations, as already noted. It is

    -2ψ_,ur-2 rψ_,u+ A(ψ_,rr+2
      rψ_,r)+2m r^2ψ_,r +1
      r^2∇^a∇_aψ=0.

Here and in the following, we write all equations in covariant form
with respect to the coordinates x^a on S^2. Following
<cit.>, we denote by q_ab the abstract round unit
metric on S^2, by q^ab its inverse, and by ∇_a the
covariant derivative with respect to q_ab, so
∇_aq_bc=0. Note that ∇_a commutes with partial
∂_u and ∂_x.

Following <cit.>, we define the reduction variables

    P:=(rψ)_,r,      Q_a:=∇_aψ,

and obtain the first-order system

    2P_,u- AP_,r  -1 r∇^aQ_a   =   2m r^2(P-ψ), 
    
    Q_a,r-1 r∇_a P   =   -1 rQ_a, 
    ψ_,r   =   1 r(P-ψ).


Let x^a=(θ,φ) be the usual coordinates on S^2, in terms
of which q_ab=diag(1,sin^2θ). Then

    ∇^a∇_aψ=ψ_,θθ+θ ψ_,θ
    +1sin^2θψ_,φφ.

To make the non-diagonal matrices C^θ and C^φ
symmetric, we need to expand the covector Q_a in components with
respect to the non-coordinate, orthonormal basis

    ∂_θ,   ∂_φ̂:=1sinθ ∂_φ.

The system becomes 

    2P_,u- AP_,r      
    
    -1 r(Q_θ,θ
    +1sinθQ_φ̂,φ)   =   2m r^2(P-ψ)+1 rθ Q_θ, 
    
    
    Q_θ,r-1 rP_,θ   =   -1 rQ_θ, 
    
    Q_φ̂,r-1 rsinθ P_,φ   =   -1 rQ_φ̂, 
    ψ_,r   =   1 r(P-ψ).

The matrices C^u and C^r are diagonal and C^θ and
C^φ are now symmetric. Replacing the coordinate θ by
y:=-cosθ gets rid of the θ term in D, see
<cit.> (s there is -y here).

However, a more elegant approach to both establishing
  the symmetry of the C^a and keeping track of Christoffel terms
  from covariant derivatives is to keep the equations covariant on
S^2. From (<ref>-<ref>) we read off

    j^u   =   2P^2, 
     
    j^r   =   - AP^2+Q^aQ_a+ψ^2, 
     
    j^a   =   -2
      rPQ^a, 
     
    S   =   4m r^2P^2-2 r(Q_aQ^a+ψ^2)+A
      rPψ.

It is easy to check that

    j^u_,u+j^r_,r+∇_aj^a=S

holds if and only if (<ref>-<ref>) hold, where
∇_aj^a is the covariant divergence. This means that our
integration measure dV must contain the covariant measure dΩ
on the round two-sphere, or dΩ=sinθ dθ dφ in
the standard coordinates.





 §.§ The vacuum Einstein equations, linearised in Bondi gauge
  about Schwarzschild





In <cit.>, the metric is written as

    ds^2   =   -e^2βV rdu^2-2e^2β du dr 
          +r^2 h_ab(dx^a-U^a du)(dx^b-U^b du).

Keeping in mind that y=-cosθ and S:=1-y^2=sin^2θ, we
read off the identifications of the metric components of
Secs. <ref> and <ref> given in
Table <ref>. 









We denote the perturbations of V, β, h_ab and U^a
about the vacuum Schwarzschild solution by
Ṽ, β̃, h̃_ab and Ũ^a.  (In
contrast to <cit.>, we have added the tildes on the
perturbations β̃ and Ũ^a to distinguish them
from the full variables.) We also replace
  Ũ_a and Ṽ by 

    ũ_a   :=   rŨ_a, 
    ṽ   :=   Ṽ r.

This makes all variables dimensionless, and therefore all lower-order terms
become proportional to 2m/r^2 or 1/r. We introduce the reduction
variables of <cit.>, which are

    P_ab   :=   (rh̃_ab)_,r, 
    
    Q_a   :=   (rũ_a)_,r-2β̃_,a, 
    
    J_a   :=   ∇^bh̃_ab+(rũ_a)_,r, 
    
    T_a   :=   β̃_,a.

These comprise 8 first derivatives of the 6 metric perturbations.

Note that the linearisation of the Bondi gauge condition h_ab= q_ab is q^abh̃_ab=0, or in coordinates,

    h̃_φφ=-sin^2θ h̃_θθ.

By definition, P_ab is also trace-free, and so

    P_φφ=-sin^2θ  P_θθ.


To complete the system with an evolution equation for ṽ while
maintaining the symmetric hyperbolic form, we introduce the further
variables

    P_a   :=   ∇^bP_ab, 
     Q   :=   ∇^aQ_a, 
    Q̂   :=   ϵ^ab∇_aQ_b, 
     J   :=   ∇^aJ_a, 
    Ĵ   :=   ϵ^ab∇_aJ_b, 
     U   :=   ∇^aũ_a, 
    Û   :=   ϵ^ab∇_aũ_b, 
     T   :=   ∇^aT_a.

These comprise 2 additional first derivatives of the metric (namely
U and Û, for a total of 10) and 7 second
derivatives of the metric.  We move tensor
indices a,b,... on S^2 with q_ab and q^ab, and note that
this commutes with taking derivatives in u and r. ϵ_ab
is the volume form on the unit 2-sphere, with defining properties

    ∇_aϵ_bc=0,     ϵ_acϵ^b_c=q_ab.


We can decompose the vector field Q_a in terms of potentials Q and
Q̂ that are determined (up to a constant) as solutions of the
Poisson equation on the unit 2-sphere, as follows:

    Q_a   =   ∇_aQ-ϵ_a^b∇_bQ̂, 
    ∇^a∇_aQ   =    Q, 
    ∇^a∇_aQ̂   =   Q̂.

We can use (<ref>) to show that

    ∇^b(∇_aQ_b+∇_bQ_a-q_ab∇^cQ_c)
      
       =   ∇^b∇_b Q_a+Q_a 
       =   ∇_a Q-ϵ_a^bQ̂+2Q_a,

where in both lines we have used

    (∇_a∇_b-∇_b∇_a)Q^b=-R_abQ^b=-Q_a,

with R_ab=q_ab the Ricci tensor on the unit round 2-sphere. 
The potentials Q and Q̂ were introduced only to
derive (<ref>), and are not part of our system.

The evolution equations for the reduction variables already introduced
in <cit.> are

    2P_ab,u- A P_ab,r      
    
    +1 r(2∇_(aQ_b)-q_ab∇^cQ_c)    =   2m r^2(P_ab-h̃_ab), 
    
    Q_a,r+1 r∇^bP_ab   =   1 r(J_a-Q_a-2T_a+2ũ_a), 
    
    
    
    J_a,r   =   1 r(2ũ_a-4T_a), 
    ũ_a,r   =   1 r(Q_a+2T_a-u_a). 
    
    
    T_a,r   =   0, 
    h̃_ab,r   =   1 r(P_ab-h̃_ab), 
    β̃_,r   =   0.

Eqs. (<ref>-<ref>) incorporate minor
corrections of their counterparts in <cit.>, see
Appendix <ref>.  The evolution equations
for our additional variables are

    2 P_a,u- A P_a,r      
     +1 r
    (∇_a Q-ϵ_a^b∇_bQ̂)   =   -2 rQ_a
    
          -1cm  +2m r^2( P_a-J_a+Q_a+2T_a),
     Q_,r+1 r∇^a P_a   =   1 r(
      J- Q-2 T+2 U), 
    Q̂_,r+1 rϵ^ab∇_a
      P_b   =   1 r(Ĵ-Q̂+2Û), 
     J_,r   =   1 r(2 U-4 T), 
    Ĵ_,r   =   2 rÛ, 
     U_,r   =   1 r( Q+2 T- U), 
    Û_,r   =   1 r(Q̂-Û), 
     T_,r   =   0, 
    ṽ_,r   =   1 r(1 2 J- T+
      U+2β̃-ṽ).


The explicit matrices C^μ are given in
Appendix <ref>. For the more elegant
covariant-on-S^2 approach to symmetric hyperbolicity, we define

    X_ab^cd:=1
      2(q_a^cq_b^d+q_b^cq_a^d-q_abq^cd),

the projection operator into the space of symmetric tracefree
2-tensors on S^2. We can then write
(<ref>,<ref>) as

    2P_ab,u- AP_ab,r+2
      rX_ab^cd∇_dQ_c   =   l.o., 
    Q_c,r+1 rX_abc^d∇_d P_ab   =   l.o.,

and the conserved current as

    j^u   :=   2 P^† P
    
    j^r   :=   - A P^† P+ Q^† Q
    
    j^d   :=   2 rX^abcdP_abQ_c+1 r(q^da
      Q+ϵ^daQ̂)P_a, 
    
    S   =   1 r(...)+2m r^2(...),

where we have defined the shorthands

    P^† P   :=   1 2P_abP^ab+P_aP^a 
       =   P_θθ^2+P_θφ̂^2+P_θ^2+P_φ̂^2, 
     Q^† Q   :=   Q_aQ^a+...=Q_θ^2+Q_φ̂^2+...,


The factor of 1/2 in front of P_abP^ab compensates for
double-counting of its algebraically independent components, see
(<ref>). The balance
law(<ref>) holds if and only if our
first-order reduction of the linearised Einstein equations holds.

For brevity, we have not written out S in full.
Our introduction of ũ_a and ṽ has the advantage,
relative to <cit.>, that in the Minkowski case m=0
D becomes 1/r times a matrix of integers, so all its
eigenvalues take the form λ_i(r)=λ̅_i/r, and d=d̅/r_0, where d̅ is the largest λ_i. 

Looking at the evolution equations, in order to add ṽ to our
system and then close it as first-order symmetric hyperbolic system we
have effectively introduced a system of variables and equations
that duplicates Frittelli's original system, but at one derivative
higher. Geometrically, this is the level of curvature, rather than of
the connection. However, our estimate includes, beside the 6 metric
perturbations, only 10 first and 7 second derivatives of the metric,
far short of the full set of Ricci rotation coefficients and null
curvature components. These are listed in
Table <ref>.

At both the level of the connection and the level of curvature, we
have in effect two pairs of wave equations, but in
different objects, (P_ab,Q_a), with P_a^a=0, and P_a, Q,Q̂. Using (<ref>) in
(<ref>) was essential for bringing the pair
(<ref>,<ref>) into a form similar to
(<ref>,<ref>), with the same matrices
C^μ.

As the background solution is spherically symmetric, all perturbations
can be split into polar and axial parts, where the axial perturbations
change sign under a reflection of S^2, or reflection in space. All
genuine scalars on S^2, in our case ψ, Ṽ and β̃ are polar. Vectors on S^2 can be split into polar and axial
parts as in (<ref>) for the example of
Q_a. Tracefree symmetric tensors can be split as, for example,

    P_ab=(2∇_(a∇_b)-q_ab∇^c∇_c)P
    +2ϵ_(a^c∇_b)∇_c P̂,

in terms of a scalar P and pseudoscalar P̂.  Axial parts are
hatted. Then axial and polar perturbations decouple from each other.

In twist-free axisymmetry, only the polar
perturbations are present. We note in passing that substituting the solution
β̃=T_a=0 of β̃_,r=0, and restricting the
background solution to the Minkowski spacetime by setting m=0, makes
Eqs. (<ref>,<ref>) equivalent to
Eqs. (D25,D26) of Paper I.

We could also write the linear perturbation equations in terms of
scalars and pseudoscalars only, in order to explicitly decouple polar
and axial perturbations. This form of the system would not be
symmetric hyperbolic, as all angular derivatives would appear in the
form ∇^c∇_c, but it could be made symmetric hyperbolic
again by re-introducing first angular derivatives as reduction
variables, as for the scalar wave equation. However, this would mean
duplicating all vectors, for example defining Q_a=∇_aQ to be a
true vector (polar) and adding Q̂_a=ϵ_ab∇^bQ̂
as a pseudovector (axial).

Looking back, to merely write the second-order Einstein equations
first-order form with a minimum number of variables, one already has
to introduce all the reduction variables (first derivatives of the
metric) (<ref>-<ref>). This is true if ṽ is
included or not. There is one reduction constraint

    ∇^bP_ab=(r∇^bh̃_ab)_,r=(r(J_a-Q_a-2T_a))_,r.

The second derivative (r∇^bh̃_ab)_,r that can be thus
written in two ways appears in only one place, and so there is a
one-parameter family of first-order reductions with inequivalent
principal parts. It turns out this contains the symmetric hyperbolic
first-order form of the equations with the equation for ṽ
excluded found by Frittelli. To bring the full system, with ṽ
included, into a first-order symmetric hyperbolic form, we had to
further add all the variables (<ref>-<ref>).




 




 §.§ Estimates for characteristic initial value and
  initial-boundary value problems





We now use the symmetric hyperbolic form of the linearised Einstein
equations in Bondi gauge to obtain energy estimates on control volumes
of interest. We introduce the Schwarzschild time coordinate

    t:=u+r_*  ⇒   dt=du+dr_*=du+ A^-1dr,

where r_* is the usual tortoise radius, and where we have defined the
shorthand

    A :=1-2m r.

Following <cit.>, we observe that

    C^t:=C^μ(dt)_,μ=C^u+ A^-1C^r,

where C^u and C^r are explicitly given in
(<ref>,<ref>), is positive definite on R^N with
smallest eigenvalue 1 for our system, independently of r, so for
the smallest eigenvalue anywhere in V we have c=1. Equivalently,
with j^u and j^r given by (<ref>) and (<ref>), we have

    j^t:=j^u+ A^-1j^r=
     P^† P+ A^-1 Q^† Q,

where P^† P and Q^† Q were
defined in (<ref>-<ref>). 

The estimates in <cit.> are
on the control volume V_1 shown in Fig. <ref>: the product of S^2
with the triangle bounded by u=u_0, r=r_0, and t=t_1. We have

    -∫_u=u_0j^u-∫_r=r_0j^r
    +∫_t=t_1j^t=∫_V_1 S.

The signs come from du at u=u_0 and dr at r=r_0 pointing into
V, and dt at t=t_1 pointing out. We will ignore the flux out of
r=r_0 (as is standard practice for initial-boundary value problems),
and so

    ||in||^2   =   ∫_u=u_02 P^† P
    +∫_r=r_0 Q^† Q,
    
    ||out'||^2   =   ∫_t=t_1 P^† P+ A^-1 Q^† Q
    
    ||out”||^2   =   ∫_r=r_0 A P^† P.

Hence, from the general formula (<ref>), and with
c=1, our estimate is

    ∫_t=t_1(
     P^† P+ A^-1 Q^† Q) 
       ≤    e^d(r_0) (t_1-t_0)(
     ∫_u=u_02 P^† P
    +∫_r=r_0 Q^† Q).

Recall that d is defined as the largest positive eigenvalue of
D in V. All elements of D are constants times
1/r or 2m/r^2, so d depends only on r_0, the smallest
value of r in V, and we have written d=d(r_0) to emphasise this.
On the Minkowski background, the 2m/r^2 terms are absent, and so
d=d̅/r_0, where d̅ is the largest eigenvalue of r
  D (with m=0). 

A second control volume of interest is shown in Fig. <ref>: the product
of S^2 with the null rectangle triangle bounded by u=u_0 and
v=v_0, and u=u_1 and v=v_1. The null coordinate v on the
Schwarzschild background is defined by

    v:=u+2r_*  ⇒   dv=du+2 A^-1dr,

and hence 

    j^v:=j^u+2 A^-1j^r=2 A^-1 Q^† Q.

The corresponding estimate is

    ∫_u=u_12 P^† P
    +∫_v=v_12 A^-1 Q^† Q
       ≤    e^d(r_0) (t_1-t_0)(
     ∫_u=u_02 P^† P
    +∫_v=v_02 A^-1 Q^† Q).


A third control volume of interest is shown in Fig. <ref>: the
product of S^2 with the null right trapezoid bounded by u=u_0, r=r_0, and
u=u_1 and v=v_1. Again we will ignore the flux out of r=r_0. The
corresponding estimate is

    ∫_u=u_12 P^† P
    +∫_v=v_12 A^-1 Q^† Q
       ≤    e^d(r_0) (t_1-t_0)(
     ∫_u=u_02 P^† P
    +∫_r=r_0 Q^† Q).














A hypothetical fourth estimate would be (<ref>) with m=0
and r_0=0, that is the pure null initial-value problem on an
outgoing cone with regular vertex at r=0. (We need to set m=0 so the
background has a regular centre). This is the PDE problem we have
been considering in Papers I and II in this series. Regularity on the
central worldline r=0 requires that all metric variables take their
Minkowski values. Hence their perturbations, and in particular the
right-moving perturbations Q vanish at r=0. Independently,
r=0 is no longer a boundary. Hence the last term in
(<ref>) would be absent.

However, we recall that on Minkowski spacetime d(r_0)=d̅/r_0,
with d̅≃ 5.0>0 . Hence exp d(r_0) (t_1-t_0) grows
arbitrarily rapidly with t as r_0→ 0, so the limit r_0→ 0 of
(<ref>) as written does not exist. This problem already
arises for the scalar wave equation on Minkowski. In
Appendix <ref> we attempt some simple ways around this
problem, and show that one of them works for the scalar wave equation
on Minkowski, but none work for the metric perturbations. 
In Appendix <ref> we carry out some
numerical tests of the hypothetical estimate on a regular null cones.





§ CONCLUSIONS





In the Introduction, we motivated the desirability of a well-posedness
proof for formulations of the Einstein equations on null cones, not
just geometrically, but for specific formulations of the Einstein
equations that are used in numerical relativity. These formulations
use “Bondi-like” coordinates, where coordinate lines of constant
(u,θ,φ) are outgoing null rays, and where we evolve only
free data and solve ODEs along the null rays to reconstruct the full
metric on each time slice.

One of these formulations, using Bondi coordinates, has been used
successfully for Cauchy-characteristic matching <cit.>
in full generality, see for example <cit.>. This and the
incomplete well-posedness result of <cit.> appeared to
be in tension with recent results <cit.>
that showed that minimal first-order reductions of this and similar
formulations are not strongly hyperbolic, and which we have verified
in Sec. <ref>.

Our tentative resolution is that the characteristic initial-boundary
value problem is well-posed in L^2 of the metric plus selected
first and second derivatives. This is the “skewed” norm whose
existence was conjectured in <cit.>. In
Sec. <ref>, we have proved this for the
linearisation of Bondi gauge about Schwarzschild, for the null
initial-boundary value problem and the double null initial value
problem. 

Based on this proof, we conjecture that 1) the initial-boundary value
problem is well-posed for the linearisation about an arbitrary
background; 2) this holds for any Bondi-like gauge; 3) this holds also
for the initial value-problem on a null cone with regular vertex.

We expect that generalising the symmetric hyperbolic form of the
linearised Einstein equations to an arbitrary background solution can
be done because we expect all additional terms to be of lower order.
The existence of a symmetric hyperbolic form of the linearised
Einstein equations in other gauges seems plausible because they can be
written as two wave equations coupled to transport equations along the
outgoing null cones, coupled only through lower-order terms.

The third conjecture appears to be the most challenging.
Unfortunately, the methods of <cit.> that
we have used for the estimates above do not allow an estimate for the
initial value problem with initial data on a regular null cone.
Chrusciel has given an existence proof for solutions of the Einstein
equations with initial data specified on a null cone with regular
vertex <cit.>. This uses harmonic coordinates, relying
on results of Dossa <cit.> for quasilinear wave
equations. Those proofs suggest that we may need to split off the
lowest powers of r and corresponding lowest spherical harmonics
Y_lm as an approximate solution (resulting in a polynomial in x
and y), and control only the remainder in an energy norm.

For a proof of well-posedness of Cauchy-characteristic matching
along the lines of the present paper to succeed, the set of variables
that are passed between the Cauchy code and the null code probably has
to coincide with the set of variables that appear at the matching
surface in the separate well-posedness estimates on both sides.

Setting aside the difficulties with proving well-posedness, one may
wonder if the linearised Einstein equations in Bondi coordinates on
null cones emanating from a regular centre are actually already
well-posed in the norm ∫ P^† P dS on those
null cones, and only our estimates for the lower-order terms are not
sharp enough to see this. Numerical experiments with the code of
<cit.> are described in
Appendix <ref>. While these cannot of
course prove stability, they seem to be compatible with it.





The author is grateful to Jonathan Luk, David Hilditch, Thanasis
Giannakopoulos and Piotr Chrusciel for helpful discussions, and
particularly to TG for detailed comments on the manuscript.











§ A TEXTBOOK EXAMPLE OF A WEAKLY HYPERBOLIC SYSTEM





A textbook example of a PDE system that is weakly but not
strongly hyperbolic is given in <cit.>, pp. 29-30: 

    ϕ_,t=([ λ 1; 0 λ;   ])ϕ_,x,     ϕ:=([ ϕ_1; ϕ_2;     ]).

It has the general solution

    ϕ(t,x)=∫_-∞^∞[
    a_(k)([ 1; 0;   ])
    +b_(k)([ ikt;   1;     ])]e^ik(x+λ t) dk.

Clearly this is not well-posed in the L^2 norm
∫(ϕ_1^2+ϕ_2^2) dx because the ikt term grows arbitrarily
rapidly in t for k arbitrarily large. The very simplest system
with this property occurs for λ=0, so that ϕ_2,x=0 and
ϕ_1,t=ϕ_2,x.

Such systems are called “weakly well-posed”, which means they are
not well-posed in L^2 but can be made well-posed if a higher
derivative norm is used for the initial data than for the solutions,
see <cit.>, p. 39. In this example, we would include
ϕ_2,x^2 in the norm at t=0 only, or we could introduce an
additional variable ϕ_3:=ϕ_2,x. They can become completely
ill-posed if lower-order terms are included, see
<cit.>, p. 40 for an example.

We give this example to make two points: 1) when a plane-wave ansatz
does not give all the solutions of a linear PDE system with constant
coefficients, we should look for additional solutions that are
polynomial; 2) the problem with polynomial solutions is not that they
grow but that they can grow arbitrarily rapidly in time for initial
data that oscillate arbitrarily rapidly in space. (In our null toy
problem, read x for “time” and y̅ for “space”.)





§ NOTES ON <CIT.>





Eq. (<ref>) corrects
the right-hand side of (15a) of <cit.> by a factor of
-1/2. This error occurs already between the nonlinear field equation
(8) of <cit.> and its linearisation
(11a). 

Eq. (<ref>) corrects the right-hand side of (15b) of
<cit.> by the addition of
Q_a+2T_a. Eq. (<ref>) corrects the right-hand side of
(15e) of <cit.>, by the subtraction of
Q_a+2T_a. These errors occurs between (11c) and (15b) and (15e),
respectively. Note there is no error between (6b) and its
linearisation (11c).

The last two errors are related because the definitions (<ref>)
and (<ref>) give rise to the integrability condition

    ∇^bP_ab=(r∇^bh̃_ab)_,r=(r(J_a-Q_a-2T_a))_,r,

or, using (<ref>),

    rQ_a,r+∇_bP_ab=rJ_a,r+J_a-Q_a-2T_a.

Using the correct expression (<ref>) for rJ_a,r
then gives the correct equation (<ref>) for Q_a,r.

(<ref>) is the covariant form of Eq. (15c) of
<cit.>. A factor of 2 is missing from the last two terms
of Eq. (12) of <cit.>, but is restored in Eq. (15c).

A second set of minor corrections concerns the explicit matrices
C^μ.  With ϕ^†:=(P_ab,Q_a,...), Frittelli states
that C^u=diag(1,1,0...,0) and
C^r=diag(- A,1,1,...,1), but this cannot be true as stated
for any choice of normalisation of the equations. 

The matrices C^a are not given explicitly, but become symmetric,
without terms sin^2θ appearing in C^u and C^r, only if one
introduces frame components on the 2-sphere as in (<ref>). With the
variables in the order

    ϕ^†=(P_θθ,P_θφ̂,Q_θ,Q_φ̂,
     P_θ, P_φ̂, Q,Q̂,...),

where the dots stand for all other variables, in any order, then with
our normalisation of the equations the matrices are

    C^u   =   diag(C^u_4,C^u_4,0,...,0), 
    
    C^r   =   diag(C^r_4,C^r_4,1,...,1), 
    
    C^θ   =   diag(C^θ_4,C^θ_4,0,...,0), 
    
    C^φ   =   1sinθdiag(C^φ̂_4,C^φ̂_4,0,...,0),

where we have defined

    C^u_4   =   diag(2,2,0,0), 
    
    C^r_4   =   diag(- A,- A,1,1), 
    
    C^θ_4   :=   1 r([ 0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0;   ]) 
    
    C^φ̂   :=   1 r([  0  0  0 -1;  0  0  1  0;  0  1  0  0; -1  0  0  0;    ]).






§ SYMMETRIC HYPERBOLIC FORM OF THE TOY MODEL





To link the considerations of Sec. <ref> and
Sec. <ref>, in this Appendix we construct a
symmetric hyperbolic reduction of the toy model of
Sec. <ref> in Bondi gauge, using the methods of
Sec. <ref>.

We introduce the reduction variables 

    P   :=   δ f_,x, 
    
    Q   :=   δ b_,x+δ G_,y̅, 
    
    J   :=   δ f_,y̅-δ b_,x, 
    
    T   :=   δ G_,y̅, 
    
    L   :=   J_,y-T_,y

Here P, Q, J and T are closely related to P_ab, Q_a,
J_a and T_a of Sec. <ref>. They do not carry
indices because of the restriction to twist-free axisymmmetry. 

A symmetric hyperbolic first-order reduction of the system
(<ref>-<ref>) in terms of these variables is

    2P_,u-P_,x-Q_,y̅   =   0, 
    
    Q_,x-P_,y̅   =   0,
    δ b_,x   =   Q-T, 
    
    J_,x   =    0, 
    
    T_,x   =    0, 
    δ f_,x   =    P, 
    δ G_,x   =    0, 
    
    L_,x   =   0, 
    δ H_,x   =   L,

in the variables

    ϕ^†:=(P,Q,δ
      b,J,T,δ f,δ G,L,δ H),

The three matrices C^μ are symmetric and C^t:=C^u+C^x is
positive definite. 

The bad solution δ b=-δ G_0,y̅ x of
the toy model becomes δ b=-T_0x in the
symmetric-hyperbolic form, and the bad solution δ
  H=δ f_0,y̅y̅ x becomes δ
  H=L_0x. To see that the linear growth in x is by itself in
conflict with our estimates, consider a toy model of the toy model,
namely

    T_,x   =   0, 
    
    b_,x   =   -T,

where for simplicity T and b depend only on u and x. The
general solution is

    T   =   T_0(u), 
    
    b   =   b_0(u)-T_0(u)(x-x_0).

Comparing to our general framework for estimates, we see that c=1,
d=1, Q:=(T,b) and P is absent. Hence the estimate
(<ref>) reduces to 

    ∫_t=t_1 Q^† Q≤ e^t_1-t_0∫_x=x_0 Q^† Q.

We can use u=t-x, u_1=t_1-x_0, u_0=t_0-x_0 to write both
sides as integrals over u, namely 

    ∫_t=t_1 Q^† Q   =   ∫_u_0^u_1[T_0(u)^2+b_0(u)^2] du, 
    ∫_x=x_0 Q^† Q   =   ∫_u_0^u_1[T_0(u)^2 
          +(b_0(u)-T_0(u)(u_1-u))^2 ] du.

One can show that

    ∫_t=t_1 Q^† Q∫_x=x_0 Q^† Q   ≤    1+Δ^2 2(1+√(1+Δ^2 4)), 
       ≤    e^t_1-t_0 for  t_1≥ t_0,

where Δ:=u_1-u_0=t_1-t_0, and the first inequality is sharp.
We should stress again that this is only a toy model, similar to the
one considered in Sec. 6.4 of <cit.>.





§ ATTEMPTS TO EXTEND THE ESTIMATES TO R_0=0





We attempt to overcome the problem that on the Minkowski background
d=d̅/r with d̅>0, which means that we cannot take r_0→
0 in any of our estimates. In this Appendix we try three simple ways of
overcoming this, first by a change of integration measure dV, as in
Eq. (<ref>).  The most general weight compatible
with the time translation and rotation symmetries of the Schwarzschild
or Minkowski background is

    dV=du ω(r) dr dΩ,

where dΩ=sinθ dθ dφ is the integration
weight on S^2 induced by the round unit metric q_ab. The
choice ω(r)=r^2 gives dV induced by the background
Schwarzschild metric g_μν. 

From (<ref>) we see that both d and c aquire
a factor of ω, so ω cancels out of the ratio d/c that
appears in our estimates. This leaves us with the addition of C^r to
play with. For simplicity, we make the ansatz ω(r)=r^p, so that
(<ref>) becomes

    S_ω=S+p r j^r=ϕ^†( D+p rC^r)ϕ.

For the wave equation on Minkowski, we find that d̅=0 for p=1
only, while it is positive for all other p. It is also positive on
Schwarzschild for all r<∞, for all p. So for the wave
equation on Minkowski only, the unique choice ω=r gives us
d=0. 

For the linearised Einstein equations on the Minkowski background, we
find that the largest eigenvalue d̅ of r D+p C^r has a
minimum of d̅≃ 3.4 at p≃ -2.2, that is, we cannot
make d̅ non-positive with any choice of p. Hence we cannot
make d̅ non-positive in this way.

For a second attempt, we observe that, with β̃_,r=0 and
β̃=0 at regular centre, β and its angular derivatives T_a
and T vanish identically, as they must vanish at the
centre. If we experimentally take them out of the system, then d̅≃ 3.1. If we additionally bring ω(r)=r^p into play, d̅
has a minimum of d̅≃ 2.6 at p≃ -1.2. Again, this is
not enough to make d̅ non-positive.

For a third attempt, we attempt to change the ratio d/c by a linear
recombination of variables as in Eq. (<ref>). We have already seen
that the eigenvalues of D and C^t can change if A is not
orthogonal. For the full system we cannot try out all possible
matrices A by brute force, but we can for the scalar wave equation
on Minkowski, as we shall see now. 

From (<ref>-<ref>) we read off
that on Minkowski C^t:=C^u+C^r is the unit matrix, while

    D=D+D^†=1 r([  0  0  0  1;  0 -2  0  0;  0  0 -2  0;  1  0  0 -2;    ]).

We have d̅=√(2)-1≃ 0.4. We have also just seen that
with p=1 we can make this d̅=0. However, we now focus on
linear recombinations of the variables other than by an overall
r-dependent factor.

Geometrically, it makes no sense to mix the components of the vector
Q_a either with each other or with the scalars P and ψ,
and we can fix an overall factor in A by leaving them completely
unchanged, so A acts non-trivially only on the pair
(P,ψ). Hence we can assume

    A=([ α 0 0 β; 0 1 0 0; 0 0 1 0; γ 0 0 δ;   ]).

Two of the eigenvalues D̃:=A^† DA are then
always -2/r for any A. As they are negative, this is not a
problem. Furthermore, if we write A=A̅ R where R is
orthogonal, R does not change the eigenvalues of C̃^t and
D̃. Hence we can choose the rotation R between P and
ψ such that
γ=0, and without loss of generality we can then rescale ψ to
set δ=1.  The non-trivial eigenvalues of D̃ are now
(β-1)±√((β-1)^2+α^2). We have already used up the rotation
between P and ψ that would allows us to set α=0, and hence
one of the eigenvalues of D̃ is always positive.

In summary, none of three relatively trivial ways of changing d/c in
the estimate manage allows us to make it non-positive, and as d/c is always
proportional to 1/r, we cannot let r_0→ 0 in estimates based on
our symmetric hyperbolic reduction.





§ NUMERICAL EXPERIMENTS





We have not been able to obtain an estimate (for the linearised
Einstein equations) for the PDE problems we have been solving
numerically in <cit.>, where
our time slices are outgoing null cones emanating from a regular
  centre.  The reason for this was that the factor exp
d(r_0) (t_1-t_0) in our estimates diverges as r_0→ 0. However,
we did not actually notice an instability in our code. In this
appendix, we test the code a little harder by running it with
small amplitude but random initial data, and see if there is a
discrete norm that remains bounded in by its initial value. The
control volume is as in Fig. <ref>, but with r_0=0,
understood as a regular centre.

We now construct a hypothesis to test by modifying the estimate
(<ref>) as follows. r=0 is not actually a boundary, but
an interior world line, and no free data can be imposed there, so the
term at r=r_0 disappears. In the setup of our code we are not
interested in controlling the energy leaving the control volume
through v=v_1, so we drop that term from the estimate
(<ref>). The hypothesis we want to test is that the function
norms on u=0 and u=u_1 are the ones that we derived but that an
unknown function K(u_1-u_0) takes the place of exp
d(r_0) (u_1-u_0). Hence our hypothesis is

    ∫_u=u_1 P^† P≤ K(u_1-u_0) ∫_u=u_0 P^† P,

where P^† P was defined above in
(<ref>).  

Restricting to vacuum in twist-free axisymmetry,
and translated into the notation of our code, the two terms in
P^† P are

    1 2P^abP_ab   :=   1 2q^acq^bd
    (rh̃_ab)_,r(rh̃_cd)_,r
       =   (rh̃_θθ)_,r^2 
       =   4S^2(rf)_,r^2

and

    P^aP_a   :=   q^ab(r∇^ch̃_ac)_,r
    (r∇^dh̃_bd)_,r
       =   [r(h̃_θθ,θ+2θ h̃_θθ)]_,r^2
       =   4S[S(rf)_,ry-8y(rf)_,r]^2.

(Recall y:=-cosθ and S:=1-y^2.) Note, however, that δ H forms
part of the variables Q, that these are not coming in through
r=0, and we do not control their leaving through v=v_1. Hence
instead of our full estimate we can use Frittelli's truncated
estimate, with ony the term (<ref>) in the integrand.

In the linearised equations, the scalar matter field decouples from
the metric perturbations, and we conjecture that it obeys a similar
estimate with integrand

    P^2:=(rψ)_,r^2.


Our integration measure in axisymmetry is 

    ∫ dS=∫ dΩ dR=2π∫_-1^1dy ∫_0^x_maxR_,x dx,

and we now drop the factor 2π. Our code requires the outer
boundary to be null or future spacelike, so we cannot literally run in
Bondi gauge R=x. However, runnning in lsB gauge and setting
x_0=x_max will turn the outer boundary x=x_max into
an ingoing null In weak gravity, this will result in R=R(u,x)≃
c(u)x, with c(u) determined so that x=x_0 is ingoing null. We can
then identify x=x_0 with v=v_1. See <cit.> for
details.

To test our conjecture, we therefore evaluate the norms (at time u)
defined by

    ||(rψ)_,r||^2   :=   ∫_-1^1dy ∫_0^x_max[(Rψ)_,x]^2 R_,x dx, 
    
    ||(rf)_,r||_1^2   :=   ∫_-1^1 (1-y^2) dy ∫...f...

and their discretizations

    ||(rψ)_,r||_1^2   ≃    2 ∑_j=1^N_yA^(0)_l=0,j∑_i=1^N_x(R_iψ_i,j-R_i-1ψ_i-1,j)^2
      R_i-R_i-1,  
    
    
    ||(rf)_,r||_1^2   ≃    2∑ A^(0)_l=0,j(1-y_j^2)∑...f...,

where R_i has only one index as R=R(u,x) in lsB gauge as R does
not depend on y. We have also used that 1/2∫ψ dy is the
l=0 component of ψ, and have used the analysis matrix A^(0)
of our pseudospectral framework to determine it. Our radial grid
starts at i=1, but in the formulas we use R=0 at i=0. See again
<cit.> for details.

We have evolved with noise of amplitude 10^-10 in f and ψ,
at N_x=64...1024 and N_y=17...65. We set x_max=1.1 and
x_0=1, and evolve to u=0.95, when the range of R has gone down
from its initial value of 0.5 by a factor of 1-0.95=0.05. We do
not filter out high frequencies, other than the components of f with
l=N_x+1 and N_x (see <cit.> for why we do
this.) At each x_i, we also set all spherical harmonic components
with l>2i to zero (see again <cit.> for why we do
this.)

The discrete L^2 norms of ψ and f themselves do not decrease
with u, but the discretised L^2 norms ||(rψ)_,r|| and
||(rf)_,r||_1 given in (<ref>,<ref>)
do. From about u=0.01, both decrease
approximately as u^-0.35. Near the end, ||(rψ)_,r||
decreases much more rapidly, while ||(rf)_,r|| increases a bit
in an apparently random manor, then decreases again. In summary, our
numerical experiments are easily compatible with the estimate
(<ref>) for K(u_1-u_0)=1 and including only P^abP_ab,
that is

    ∫_u=u_1(rf)_,r^2 dr dy 
    ≤∫_u=u_0(rf)_,r^2 dr dy,

and a similar estimate with ψ instead of f. 

Looking at individual spherical harmonic components f_l(u,x) and
ψ_l(u,x) of f(u,x,y) and ψ(u,x,y), we see that the random
initial dominated by the grid frequency data quickly smooth out into
well-resolved data on frequencies lower than the grid frequency. Hence
our method seems to be quite dissipative. In a second phase we then
see that these appear to be “stretched” in x, as the grid in fact
zooms in on them, with R_max shrinking linearly in u.

For definiteness, we have implemented a specific continuum norm and
its discretisation, but we should stress that this was just an
informed guess. In particular, we have arbitrarily chosen
dV=dΩ dR, rather than, for example, dV=dΩ R^p dR. The
“correct” choice would depend on what (if any) estimate can be
proved. We have also not gone to particularly large values of N_x
and N_y, but only ones that we have also used in physics
simulations. (A soft upper limit on N_x is set by computing time,
while a hard limit N_y≤ 128 is set by the accuracy of our spectral
method in y.) Hence our results should only be considered as a
slightly more challenging stability test of our code motivated by the
results of this paper, not a numerical test of well-posedness.



99



FouresBruhat1952 Y. Fourès-Bruhat, Théorème
  d’existence pour certains systèmes d'èquations aux dérivées
  partielles non linéaires,
  https://gallica.bnf.fr/ark:/12148/bpt6k3182n/f618.imageActa
    Math. 88, 141 (1952), see also
  https://link.springer.com/article/10.1007/s10714-022-02917-4
Gen. Relativ. Gravit. 54, 35 (2022) for an English translation.

Pretorius2005 F. Pretorius, Evolution of binary black-hole
  spacetimes,
  https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.95.121101
       Phys. Rev. Lett. 95 121101 (2005).

Friedrich1985 H. Friedrich, On the hyperbolicity of Einstein's
  and other gauge field equations, https://link.springer.com/article/10.1007/BF0121772Comm. Math. Phys. 100, 525 (1985).

constraintdamping C. Gundlach, G. Calabrese, I. Hinder and
  J. M. Martín-García, Constraint damping in the Z4
  formulation and harmonic gauge,
  https://iopscience.iop.org/article/10.1088/0264-9381/22/17/025Class. Quant. Grav. 22, 3767 (2005).

Rendall1990 A. D. Rendall, Reduction of the characteristic
  initial value problem to the Cauchy problem and its applications to
  the Einstein equations, https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1990.0009Proc. Roy. Soc. A 427, 1872 (1990).

Luk2011 J. Luk, On the local existence for the
  characteristic initial value Problem in general Relativity,
  https://academic.oup.com/imrn/article/2012/20/4625/658337Int. Mat. Res. Notices, 20,4625-4678 (2012), accessed online as
https://arxiv.org/abs/1107.0898arXiv:1107.0898. 

Chrusciel2013 P. T. Chruściel, The existence theorem for the
  general relativistic cauchy problem on the light-cone,
  https://www.cambridge.org/core/journals/forum-of-mathematics-sigma/article/existence-theorem-for-the-general-relativistic-cauchy-problem-on-the-lightcone/AD55028BB47B55110EA9349AFC021817Forum
    of Mathematics, Sigma, 2:e10 (2014).

ChruscielJezierski2011 P. T. Chruściel and J. Jezierski,
  On free general relativistic initial data on the light cone,
  https://doi.org/10.1016/j.geomphys.2011.04.015J. of Geometry
    and Physics 62, 578 (2012).

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

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

GHZ2020 T. Giannakopoulos, D. Hilditch, and M. Zilhão,
  Hyperbolicity of general relativity in Bondi-like gauges, 
https://journals.aps.org/prd/abstract/10.1103/PhysRevD.102.064035Phys. Rev. D 102, 064035 (2020).

Siebeletal2003 F. Siebel, J. A. Font, E. Müller and
  P. Papadopoulos, Axisymmetric core collapse simulations using
  characteristic numerical relativity,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.67.124018Phys. 
Rev. D 67, 124018 (2003).

axinull_critscalar C. Gundlach, T. W. Baumgarte and
  D. Hilditch, Simulations of gravitational collapse in null
  coordinates: II. Critical collapse of an axisymmetric scalar field, https://arxiv.org/abs/2404.15839arXiv:2404.15839.

Maetal2023 S. Ma, J. Moxon, M. A. Scheel, K. C. Nelli,
  N. Deppe, M. S. Bonilla, L. E. Kidder, P. Kumar, G. Lovelace,
  W. Throwe and N. L. Vu, Fully relativistic three-dimensional
  Cauchy-characteristic matching, 
  https://arxiv.org/abs/2308.10361arXiv:2308.10361. 

Frittelli2005a S. Frittelli, Well-posed first-order reduction
  of the characteristic problem of the linearized Einstein equations,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.71.024021Phys. Rev. D
    71, 024021 (2005).

GBHPZ2022 T. Giannakopoulos, N. T. Bishop, D. Hilditch,
  D. Pollney, and M. Zilhão, Gauge structure of the Einstein field equations in Bondi-like coordinates,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.105.084055Phys.
    Rev. D 105, 084055 (2022).

GBHPZ2023 T. Giannakopoulos, N. T. Bishop, D. Hilditch,
  D. Pollney, and M. Zilhão, Numerical convergence of model
  Cauchy-characteristic extraction and matching,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.108.104033Phys.
    Rev. D 108, 104033 (2023).

axinull_formulation C. Gundlach, D. Hilditch and
  T. W. Baumgarte, Simulations of gravitational collapse in null
  coordinates: I. Formulation and weak field tests in Bondi-like
  gauges, https://arxiv.org/abs/2404.15105arXiv:2404.15105.

Frittelli2004 S. Frittelli, Estimates for the characteristic
  problem of the first-order reduction of the wave equation, https://iopscience.iop.org/article/10.1088/0305-4470/37/36/005J.
    Phys. A 37, 8639 (2004).

Frittelli2005b S. Frittelli, Estimates for first-order homogeneous linear characteristic problems,
  https://iopscience.iop.org/article/10.1088/0305-4470/38/19/010J.
    Phys. A 38, 4209 (2005).

RenardyRogers M. Renardy and R. C. Rogers, An
  Introduction to Partial Differential Equations, 2nd edition,
  Springer (New York), 2004.

GKO B. Gustafsson, H.-O. Kreiss and J. Oliger,
  Time-dependent problems and difference methods, Wiley (New
  York), 1995.

GomezPapadopoulosWinicour1994 R. Gómez, P. Papadopoulos,
  and J. Winicour, Null cone evolution of axisymmetric vacuum
  spacetimes,
  https://aip.scitation.org/doi/abs/10.1063/1.530848J. Math. Phys. 35, 4184 (1994).

WinicourLRR2012 J. Winicour, Characteristic evolution and
  matching,
  https://link.springer.com/article/10.12942/lrr-2012-2Living
    Rev. Rel. 15, 2 (2012).

KreissLorenz2004 H.-O. Kreiss and J. Lorenz,
  Initial-boundary value problems and the Navier-Stokes equations,
  SIAM (Philadelphia) 2004.

Dossa1997 M. Dossa, Espaces de Sobolev non isotropes, à
  poids et problèmes de Cauchy quasi-linéaires sur un conoïde
  caractéristique,
  http://www.numdam.org/item/AIHPA_1997__66_1_37_0/Ann. de
    l'Inst. Henri Poincaré Physique théorique, 66, 37 (1997).

Giannakopoulosthesis T. Giannakopoulos, Characteristic
  formulations of general relativity and applications, PhD thesis,
  Instituto Superior Tecnico, University of Lisbon, 2022, 
  https://arxiv.org/abs/2308.16001arXiv:2308.16001.







