Spectral Cauchy-characteristic extraction of the gravitational wave news function
Kevin Barkett,
1
Jordan Moxon,
1
Mark A. Scheel,
1
and B ́ela Szil ́agyi
2, 1
1
TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA
2
Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91106, USA
(Dated: October 23, 2019)
We present an improved spectral algorithm for Cauchy characteristic extraction (CCE) and char-
acteristic evolution of gravitational waves in numerical relativity. The new algorithms improve
spectral convergence both at the poles of the spherical-polar grid and at future null infinity, as well
as increasing the temporal resolution of the code. Key to the success of these algorithms is a new
set of high-accuracy tests, which we present here. We demonstrate the accuracy of the code and
compare with the existing Pittnull implementation.
I. INTRODUCTION
The discovery of GW150914 [1] heralded the begin-
ning of gravitational wave astronomy. In the subsequent
years that detection has been followed up by a number
of other signals observed from binary black hole (BBH)
mergers [2–5], as well as from the merger of a binary
neutron star (BNS) system [6]. As the aLIGO [7] and
Virgo [8] detectors push to ever greater sensitivities, the
number of expected observations will continue to grow.
Extracting the signals from the noise involves match-
ing the incoming data against a template bank of the-
oretically expected waveforms generated across possible
binary configurations. The efficacy of extracting the con-
figuration parameters (for instance, masses and spins of
the binary components) from a given signal depends on
the fidelity of the computed waveforms comprising the
template bank; this is because errors in the template
bank will bias the estimated parameters. The only
ab
initio
method of generating accurate theoretical wave-
forms for merging BBH systems is via numerical relativ-
ity: the numerical solution of the full Einstein equations
on a computer. Other methods of generating theoretical
BBH waveforms, such as Effective One-Body solutions [9]
and phenomenological models [10, 11], are calibrated to
numerical relativity.
One limitation of numerical relativity simulations is
that they all rely on a Cauchy approach in which the
spacetime is decomposed into a foliation of spacelike
slices, and the solution marches from one slice to the next.
Such an approach can compute the solution to Einstein’s
equations only in a region of spacetime with finite spatial
and temporal extents bounded around the compact ob-
jects, whereas the gravitational radiation is defined at fu-
ture null infinity
I
+
. While some work has gone into hy-
perboloidal compactification methods for simulating the
propagation of gravitational waves to
I
+
[12–14], these
methods have never been fully implemented in the non-
linear regime. Without them, extracting the waveform
signal from the simulations with these finite extents re-
quires additional work.
The most common method of extracting the grav-
itational radiation from a numerical relativity simula-
tion is to compute quantities such as the Newman-
Penrose scalar Ψ
4
[15] or the Regge-Wheeler and Zerilli
scalars [16] at some large but finite distance from the
near zone (perhaps 100–1000
M
, where
M
is the total
mass of the system), typically on coordinate spheres of
constant surface area coordinate coordinate
r
. Because
these quantities or the methods of computing them in-
clude finite-radius effects, these quantities are computed
on a series of shells at different radii
r
, fit to a polynomial
in 1
/r
, and then extrapolated to infinity by reading off
the 1
/r
coefficient of the polynomial [17]. As the extrac-
tion surfaces are shells of constant coordinate radii, the
choice of gauge implemented in the simulation can con-
taminate the resulting waveforms. Furthermore, if the
shells are too close to the orbiting binary, the extrapo-
lation procedure might not remove all of the near-zone
effects.
An alternative method for computing gravitational ra-
diation in numerical relativity is to solve the full Einstein
equations in a domain that extends all the way to
I
+
,
where gravitational waves can be measured. This can be
done by rewriting Einstein’s equations using a character-
istic formalism [18–20], in which the equations are solved
on outgoing null surfaces that extend to
I
+
. This for-
malism chooses coordinates that correspond to distinct
outward propagating null rays, so it fails in the dynam-
ical, strong field regime at any location where outgoing
null rays intersect (i.e., caustics). Because of this, charac-
teristic evolution is unable to evolve the near-field region
of a merging binary system, so it cannot accomplish a
BBH simulation on its own. However, it is possible to
combine an interior numerical relativity code that solves
the equations on Cauchy slices with an exterior charac-
teristic code that solves them on null slices; the determi-
nation of characteristic quantities from Cauchy data is
known as Cauchy-Characteristic Extraction (CCE) (See
Fig. 1.), and the subsequent numerical evolution of those
quantities is known as characteristic evolution.
Specifically, CCE uses the metric and its derivatives
computed from a Cauchy evolution (red region in Fig. 1)
and evaluated on a world tube Γ (thick red line) that
lies on or inside the boundary of the Cauchy region.
These quantities on the world tube are then used as in-
ner boundary data for a characteristic evolution (blue
region) based on outgoing null slices (blue curves). Be-
arXiv:1910.09677v1 [gr-qc] 21 Oct 2019
2
FIG. 1. Penrose diagram showing a typical CCE setup. The
metric is evolved using 3+1 methods in the Cauchy region
(shaded red) and with null methods in the characteristic re-
gion (shaded blue). The Cauchy and characteristic regions
overlap. Curves of constant
̆
t
or ̆
r
, the Cauchy coordinates,
are shown in red, and are shown as dashed curves outside the
Cauchy region, where they extend to spatial infinity
i
0
or fu-
ture temporal infinity
i
+
. Null curves of constant
u
are shown
in blue. Given data on an worldtube Γ (thick red curve) and
on an initial null slice (thick blue curve), the characteristic
evolution computes the full metric in the characteristic re-
gion. In Section III we describe the interface from Cauchy to
Bondi coordinates on Γ. In Section IV we describe the char-
acteristic evolution. In Section V we discuss computing the
news function at
I
+
(thick green curve) and transforming it
to coordinates corresponding to a free-falling observer.
cause the combined CCE system uses the full Einstein
equations for both the Cauchy and characteristic evolu-
tions, it produces the correct solution at
I
+
, with the
characteristic evolution properly resolving near-zone ef-
fects. The gravitational radiation is computed according
to a particular inertial observer at
I
+
(green curve).
This observer is related to any other inertial observer by
a single BMS transformation [19] (the group of Lorentz
boosts, rotations, and supertranslations [21]), so up to
this BMS transformation the waveform is independent of
the gauge chosen by the Cauchy evolution.
The first code to implement CCE and characteristic
evolution was the PittNull code [22–24]. Since its ini-
tial implementation there have been a number of im-
provements made, and the current iteration of that code
utilizes stereographic angular coordinate patches, finite-
differencing, and a null parallelogram scheme with fixed
time steps for integrating in the null and time directions.
Overall the code is second-order convergent with reso-
lution [25, 26] (although a fourth-order implementation
also exists, see [27]). Compared to waveforms computed
from a Cauchy code by evaluating Ψ
4
at finite radii and
extrapolating to
r
→ ∞
as described above, waveforms
extracted via CCE using PittNull were shown to better
remove gauge effects and to better resolve the
m
= 0
memory modes [28–30].
Currently, PittNull requires thousands of CPU-hours
to compute a waveform at
I
+
given worldtube output
from a typical Cauchy BBH simulation at multiple res-
olutions [31]. While that cost is smaller than the com-
putational expense of the Cauchy simulation, it is still
unwieldy, and is likely one reason that most Cauchy
numerical-relativity codes do not use CCE and charac-
teristic evolution despite the availability of PittNull. Be-
cause the metric in the characteristic region is smooth,
the computational cost of characteristic evolution should
be greatly reduced by using spectral methods instead of
finite differencing. Such a spectral implementation of
characteristic evolution has been introduced in the SpEC
framework [31–33]. Their tests showed improved speed
and accuracy over the finite-difference implementation of
PittNull [31, 32].
Our work here describes improvements in accuracy, ef-
ficiency, and robustness to the code described in [31–33].
In particular, we discuss an improved handling of the in-
tegration along the null slices, we clarify issues related to
the particular choice of coordinates along the null slice,
and we implement better handling of the inertial coor-
dinates at
I
+
. We demonstrate through a series of an-
alytic tests that our version of CCE and characteristic
evolution can compute waveforms with much lower com-
putational cost than PittNull. An earlier version of our
implementation has been used to probe the near-field re-
gion of a binary black hole ringdown [34].
We start with a brief summary of the Bondi metric and
the null formulation of the Einstein equations in Sec II. A
detailed explanation for how CCE and characteristic evo-
lution works can be broken up into three distinct parts:
the inner boundary formalism, the volume characteris-
tic evolution, and the
I
+
extraction, which we describe
in subsequent sections. Sec III describes the means by
which the metric known on a world tube is converted
into Bondi form to serve as the inner boundary values for
the characteristic evolution system. Sec IV discusses the
process of evolving Einstein’s equations from the inner
boundary to
I
+
. Sec V explains how to take the metric
computed on
I
+
and extract the Bondi news function
in the frame of an inertial observer at
I
+
. In Sec VI, we
describe code tests and performance.
Throughout this paper, indices with Greek letters
(
μ,ν,...
) correspond to spacetime coordinates, lowercase
Roman letters (
i,j,...
) to spatial coordinates, and capi-
talized Roman letters (
A,B,...
) to angular coordinates,
and we choose a system of geometrized units (
c
=
G
= 1).
For convenience, we have included a definitions key in ap-
pendix C.
3
II. SUMMARY OF CHARACTERISTIC
FORMULATION
In the characteristic region (see Fig. 1), we adopt a
coordinate system
x
μ
= (
u,r,x
A
), where
u
is the coordi-
nate labeling the outgoing null cones,
r
is an areal radial
coordinate, and
x
A
are the angular coordinates. Note
that a curve of constant (
u,x
A
) is an outgoing null ray
parameterized by
r
; for this reason we sometimes call
r
a
“radinull” coordinate. The metric can then be expressed
in the Bondi-Sachs form [18, 19],
ds
2
=
−
(
e
2
β
(1 +
rW
)
−
r
2
h
AB
U
A
U
B
)
du
2
−
2
e
2
β
dudr
−
2
r
2
h
AB
U
B
dudx
A
+
r
2
h
AB
dx
A
dx
B
,
(1)
where
W
corresponds to the mass aspect,
U
A
to the
shift,
β
to the lapse, and
h
AB
to the spherical 2-metric.
The quantity
h
AB
has the same determinant as the unit
sphere metric
q
AB
,
|
h
AB
|
=
|
q
AB
|
. Note that the met-
ric Eq. (1) is not constrained to be asymptotically flat,
as required by Bondi-Sachs coordinates. Instead, we im-
pose the weaker constraint that all metric components
of Eq. (1) are asymptotically finite at
I
+
. To empha-
size this subtle difference with Bondi-Sachs coordinates,
we refer to the spacetime metric as having the “Bondi-
Sachs form” rather than being expressed in Bondi-Sachs
coordinates. An additional intermediate quantity,
Q
A
, is
defined to reduce the evolution equations to a series of
1st order PDEs,
Q
A
=
r
2
e
−
2
β
h
AB
U
B
,r
.
(2)
Instead of expressing the metric in terms of tensorial
objects, we employ a complex dyad so that the metric
components can be computed as spin-weighted scalars,
and each of these scalars can be expanded in terms of
Spin-Weighted Spherical Harmonics (SWSHes) of the ap-
propriate spin weight; see Appendix A for details about
SWSHes. The dyad
q
A
has the following properties:
q
A
q
A
=0
,
(3)
q
A
̄
q
A
=2
.
(4)
If we define
q
AB
and
q
AB
such that
q
AB
=
1
2
(
q
A
̄
q
B
+ ̄
q
A
q
B
)
,
(5)
q
AC
q
CB
=
δ
A
B
,
(6)
then
q
A
=
q
AB
q
B
.
(7)
We express the metric coefficients and the quantity
Q
A
in terms of spin-weighted scalars
J
,
K
,
U
, and
Q
, defined
by
J
=
1
2
h
AB
q
A
q
B
,
(8)
K
=
1
2
h
AB
q
A
̄
q
B
,
(9)
U
=
q
A
U
A
,
(10)
Q
=
Q
A
q
A
.
(11)
The determinant condition on
h
AB
defines a relationship
between
J
and
K
as
K
=
√
1 +
J
̄
J.
(12)
We introduce one more intermediate variable
H
, the time
derivative of
J
along slices of constant
r
,
H
=
J
,u
|
x
A
,r
=const
(13)
The quantities
J,β,
and
Q
are all dimensionless while
U,W,
and
H
have units of 1
/R
(identically, units of 1
/u
in the case of
H
).
Evaluating the components of the Einstein equation
G
μν
= 0 provides a system of equations for the quantities
β
,
Q
,
U
,
W
, and
H
:
β
,r
=
N
β
,
(14)
(
r
2
Q
)
,r
=
−
r
2
(
̄
ð
J
+
ð
K
)
,r
+ 2
r
4
ð
(
r
−
2
β
)
,r
+
N
Q
,
(15)
U
,r
=
r
−
2
e
2
β
Q
+
N
U
,
(16)
(
r
2
W
)
,r
=
1
2
e
2
β
R−
1
−
e
β
ð
̄
ð
e
β
+
1
4
r
−
2
(
r
4
(
ð
̄
U
+
̄
ð
U
))
,r
+
N
W
,
(17)
2(
rH
)
,r
=((1 +
rW
)(
rJ
)
,r
)
,r
−
r
−
1
(
r
2
ð
U
)
,r
+2
r
−
1
e
β
ð
2
e
β
−
(
rW
)
,r
J
+
N
J
,
(18)
where
R
=2
K
−
ð
̄
ð
K
+
1
2
(
̄
ð
2
J
+
ð
2
̄
J
) +
1
4
K
(
̄
ð
̄
J
ð
J
−
̄
ð
J
ð
̄
J
)
,
(19)
and
N
β
,
N
W
,
N
Q
,
N
W
,
and
N
J
are the terms nonlinear
in
J
and its derivatives, as according to [22]. Appendix B
provides the full expressions for these equations.
These equations correspond to different components
of the Einstein equations, namely,
R
rr
= 0 gives the
equation for
β
,r
,
R
rA
q
A
= 0 gives the equation for
U
,r
,
R
AB
h
AB
=0 gives the equation for
W
,r
, and
R
AB
q
A
q
B
=
0 gives the equation for
H
,r
. These cover six of the ten
independent components of Einstein’s equations. As [23]
discusses in more detail, of the four remaining compo-
nents of the Einstein equations, one of these is identi-
cally zero (
R
r
r
= 0) while the other three (
R
r
u
= 0 and
R
r
A
q
A
= 0) serve as constraint conditions for the evolu-
tion on each of the null slices.
However, computing these constraint conditions in-
volve lengthy expressions that include the
u
−
derivatives
of evolution quantities other than
J
,u
. It is not straight-
forward to compute these derivatives to the same accu-
racy achieved by the rest of the code. We leave to future
4
implementations the ability to accurately compute these
constraints as a monitor of how well we obey the full
Einstein equations during the evolution.
The equations are presented in a useful hierarchical or-
der: the right-hand side of the
β
equation involves only
J
and its hypersurface derivatives, the right-hand side of
the
Q
equation involves only
J
and
β
and their hypersur-
face derivatives, and so on for the other equations. There-
fore, given data for all quantities on the inner boundary
as well as
J
on an initial
u
= constant null slice, we
can integrate the series of equations in Eqs. (14)–(18) on
that slice from the inner boundary to
r
=
∞
to obtain
β,Q,U,W
and then
H
in sequence on that slice. Then,
given
H
=
J
,u
|
r
=const
on that slice, we can integrate for-
ward in time to obtain
J
on the next null slice.
III. INNER BOUNDARY FORMALISM
The coordinates used to evolve Einstein’s equations in
the Cauchy region (red area of Fig. 1) are generally dif-
ferent from the coordinates discussed in Section II. The
Cauchy coordinates are chosen to make the interior evo-
lution proceed without encountering coordinate singu-
larities; the procedure for choosing these coordinates is
complicated and typically involves coordinates that are
evolved along with the solution [35–40]. Therefore, for
CCE we must transform from arbitrary Cauchy coordi-
nates to coordinates such that the spacetime metric takes
the Bondi-Sachs form (Eq. (1)) at the worldtube.
Here, in the Cauchy region, for simplicity we assume
Cartesian coordinates (
̆
t,
̆
x
̆ı
) in which the world tube hy-
persurface Γ (which is chosen by the Cauchy code) is a
surface of constant ̆
r
, where ̆
r
=
√
̆
x
2
+ ̆
y
2
+ ̆
z
2
.
We also define angular coordinates ̆
x
̆
A
= (
̆
θ,
̆
φ
) in the
usual way from the Cartesian coordinates ̆
x
̆ı
.
The world tube serves as the inner boundary of the
characteristic domain (see Fig. 1). On this boundary, we
assume that the interior Cauchy code provides the spa-
tial 3-metric
g
̆ı ̆
, the shift
β
̆ı
, and the lapse ̆
α
, along with
the ̆
r
and
̆
t
derivatives of each of these quantities. Angu-
lar derivatives of these quantities are necessary as well;
however, we can compute those numerically within the
worldtube itself, so they need not be provided
a priori
.
Ref. [24] describes how to take the data provided by
the interior Cauchy code and covert it into Bondi form
(Eq. (1)) to extract the inner boundary values of the
evolution quantities (
J
|
Γ
,β
|
Γ
,...
). This section is primar-
ily a summary of their results; however, we use different
notation than Ref. [24]. Additionally, as noted above,
the SpEC CCE treatment takes the inner boundary of
the domain to be the worldtube provided by the Cauchy
code, which is generally not a surface of constant
r
. The
PittNull treatment, on the other hand, uses a surface
of constant
r
as the inner boundary of the domain, and
performs a Taylor expansion in the affine radial coor-
dinate in order to determine inner boundary data on
this surface. Avoiding the Taylor expansion simplifies
the boundary computation, and may provide marginal
precision improvements by avoiding finite Taylor series
truncation error.
A. Affine Null Coordinates
Our goal is to transform from the coordinates (
̆
t,
̆
x
̆ı
) to
coordinates such that the metric takes the Bondi-Sachs
form (Eq. (1)). It is simplest to proceed in two steps: the
first step, described in this subsection, is to construct co-
ordinates foliated by outgoing null geodesics. The second
step, described in Section III B, will be to transform from
these affine coordinates to Bondi coordinates.
We begin by constructing a choice null generator
`
̆
μ
,
which involves the unit outward spatial vector normal to
the world tube’s surface,
s
̆
μ
, and the unit timelike vector
normal to a slice of constant
̆
t
,
n
̆
μ
:
s
̆
μ
=
{
0
,
g
̆ı ̆
̆
x
̆
√
g
̆ı ̆
̆
x
̆ı
̆
x
̆
}
,
(20)
n
̆
μ
=
1
̆
α
{
1
,
−
β
̆ı
}
.
(21)
Eq. (20) depends on our simplifying assumption that the
world tube is spherical in Cauchy coordinates, and can
be generalized. From these equations, the null generator
is
`
̆
μ
=
n
̆
μ
+
s
̆
μ
̆
α
−
g
̆ı ̆
β
̆ı
s
̆
.
(22)
The time derivatives of these vectors are
s
̆
μ
,
̆
t
=
{
0
,
(
−
g
̆ı ̆
+
s
̆ı
s
̆
/
2)
s
̆
k
g
̆
̆
k,
̆
t
}
,
(23)
n
̆
μ
,
̆
t
=
1
̆
α
2
{
−
̆
α
,
̆
t
,
̆
α
,
̆
t
β
̆ı
−
̆
αβ
̆ı
,
̆
t
}
,
(24)
`
̆
μ
,
̆
t
=
n
̆
μ
,
̆
t
+
s
̆
μ
,
̆
t
+
`
̆
μ
(
−
̆
α
,
̆
t
+
g
̆ı ̆
,
̆
t
β
̆ı
s
̆
+
g
̆ı ̆
β
̆ı
,
̆
t
s
̆
+
g
̆ı ̆
β
̆ı
s
̆
,
̆
t
)
̆
α
−
g
̆ı ̆
β
̆ı
s
̆
.
(25)
We will now construct a null coordinate system based
on outgoing null geodesics generated by
`
̆
μ
. Let
̄
λ
be
an affine parameter along these geodesics such that the
value of
̄
λ
on the world tube Γ is
̄
λ
|
Γ
= 0. We also define
a null coordinate ̄
u
and angular coordinates ̄
x
̄
A
= (
̄
θ,
̄
φ
)
that obey ̄
u
=
̆
t
and ̄
x
̄
A
= ̆
x
̆
A
on the world tube, and
are constant along the outgoing null geodesic generated
by
`
̆
μ
. Thus we have defined a new intermediate, affine
coordinate system, ̄
x
̄
μ
= ( ̄
u,
̄
λ,
̄
θ,
̄
φ
) and we will express
the metric
g
̄
μ
̄
ν
in these affine coordinates.
To do this, we will need to write down the coordinate
transformation from ̆
x
̆
μ
to ̄
x
̄
μ
in a neighborhood of the
world tube, not just on the world tube itself, because we
need derivatives of this transformation. In particular, we
will need derivatives with respect to
̄
λ
. The derivative
5
of the metric components
g
̆
μ
̆
ν
along the null direction
simply is
g
̆
μ
̆
ν,
̄
λ
=
`
̆
γ
g
̆
μ
̆
ν,
̆
γ
.
(26)
The evolution of the coordinates ̆
x
̆
μ
along null geodesics
implies that in a neighborhood of the world tube
̆
x
̆
μ
,
̄
λ
=
`
̆
ν
∂
̆
ν
̆
x
̆
μ
=
`
̆
μ
.
(27)
Given the new coordinates ̄
x
̄
μ
, the metric components
in these coordinates are
g
̄
μ
̄
ν
=
∂
̆
x
̆
α
∂
̄
x
̄
μ
∂
̆
x
̆
β
∂
̄
x
̄
ν
g
̆
α
̆
β
.
(28)
On the world tube,
∂
̆
t
∂
̄
x
̄
A
= 0
∂
̆
x
̆
i
∂
̄
x
̄
A
=
∂
̆
x
̆
i
∂
̆
x
̆
A
∂
̆
t
∂
̄
u
= 1
,
∂
̆
x
̆ı
∂
̄
u
= 0
,
(29)
where the term
∂
̆
x
̆ı
/∂
̄
x
̄
A
is the standard Cartesian to
spherical Jacobian. The above values of the Jacobians
hold only on the world tube. In addition to the metric
itself, we will also need first derivatives of the metric,
including the derivative with respect to
̄
λ
. This requires
the
̄
λ
derivatives of the Jacobians evaluated on the world
tube, which we represent here as
∂
2
̆
x
̆
μ
∂
̄
x
̄
A
∂
̄
λ
=
∂`
̆
μ
∂
̄
x
̄
A
=
`
̆
μ
,
̄
A
,
∂
2
̆
x
̆
μ
∂
̄
u∂
̄
λ
=
∂`
̆
μ
∂
̄
u
=
`
̆
μ
,
̄
u
,
(30)
where we have made use of Eq. (27).
We are now ready to write out the metric in these inter-
mediate coordinates by taking the expression in Eq. (28)
and taking the appropriate derivatives,
g
̄
u
̄
λ
=
−
1
,
g
̄
λ
̄
λ
=
g
̄
λ
̄
A
= 0
,
g
̄
u
̄
u
=
g
̆
t
̆
t
,
g
̄
u
̄
A
=
∂
̆
x
̆ı
∂
̄
x
̄
A
g
̆ı
̆
t
,
g
̄
A
̄
B
=
∂
̆
x
̆ı
∂
̄
x
̄
A
∂
̆
x
̆
∂
̄
x
̄
B
g
̆ı ̆
,
g
̄
A
̄
B,
̄
λ
=
∂
̆
x
̆ı
∂
̄
x
̄
A
∂
̆
x
̆
∂
̄
x
̄
B
g
̆ı ̆
,
̄
λ
+
(
`
̆
μ
,
̄
A
∂
̆
x
̆ı
∂
̄
x
̄
B
+
`
̆
μ
,
̄
B
∂
̆
x
̆ı
∂
̄
x
̄
A
)
g
̆
μ
̆ı
,
g
̄
A
̄
B,
̄
u
=
∂
̆
x
̆ı
∂
̄
x
̄
A
∂
̆
x
̆
∂
̄
x
̄
B
g
̆ı ̆
,
̆
t
,
g
̄
u
̄
A,
̄
λ
=
`
̆
μ
,
̄
A
g
̆
t
̆
μ
+
∂
̆
x
̆ı
∂
̄
x
̄
A
(
g
̆ı
̆
t,
̄
λ
+
`
̆
μ
,
̄
u
g
̆ı ̆
μ
)
,
(31)
and
g
̄
u
̄
u
=
g
̄
u
̄
A
= 0
,
g
̄
u
̄
λ
=
−
1
,
g
̄
A
̄
B
g
̄
B
̄
C
=
δ
̄
A
̄
C
,
g
̄
λ
̄
A
=
g
̄
A
̄
B
g
̄
u
̄
B
,
g
̄
λ
̄
λ
=
−
g
̄
u
̄
u
+
g
̄
λ
̄
A
g
̄
u
̄
A
,
g
̄
A
̄
B
,
̄
λ
=
−
g
̄
A
̄
C
g
̄
B
̄
D
g
̄
C
̄
D,
̄
λ
,
g
̄
λ
̄
A
,
̄
λ
=
g
̄
A
̄
B
(
g
̄
u
̄
B,
̄
λ
−
g
̄
λ
̄
C
g
̄
B
̄
C,
̄
λ
)
.
(32)
B. Bondi Form of Metric
Given the intermediate null coordinates and the metric
in that coordinate system, we apply one last coordinate
transformation to put the spacetime metric in Bondi-
Sachs form (Eq. (1)). We define coordinates (
u,r,θ,φ
),
where
r
is a surface area coordinate,
u
= ̄
u
,
θ
=
̄
θ
, and
φ
=
̄
φ
. The surface area coordinate
r
is defined by
r
=
(
|
g
AB
|
|
q
AB
|
)
1
4
=
(
|
g
̄
A
̄
B
|
|
q
̄
A
̄
B
|
)
1
4
,
(33)
where
q
̄
A
̄
B
is the unit sphere metric.
The components of the metric in Bondi coordinates are
then
g
μν
=
∂x
μ
∂
̄
x
̄
α
∂x
ν
∂
̄
x
̄
β
g
̄
α
̄
β
.
(34)
The Jacobians include the derivatives of the surface area
coordinate
r
. We compute
r
,
̄
α
=
r
4
(
g
̄
A
̄
B
g
̄
A
̄
B,
̄
α
−
|
q
̄
A
̄
B
|
,
̄
α
|
q
̄
A
̄
B
|
)
.
(35)
Since the only difference between the final boundary
coordinates (
u,r,θ,φ
) and intermediate coordinates is
the choice of radinull coordinates, the Jacobians for the
u
,
θ
and
φ
directions are trivial. Eq. (32) gives us
g
uu
=
g
̄
u
̄
u
= 0
,
g
uA
=
g
̄
u
̄
A
= 0
,
g
AB
=
g
̄
A
̄
B
.
(36)
The other metric components are
g
ur
=
∂r
∂
̄
x
̄
μ
g
̄
u
̄
μ
=
−
r
,
̄
λ
,
g
rr
=
∂r
∂
̄
x
̄
μ
∂r
∂
̄
x
̄
ν
g
̄
μ
̄
ν
=
(
r
,
̄
λ
)
2
g
̄
λ
̄
λ
6
+2
r
,
̄
λ
(
r
,
̄
A
g
̄
λ
̄
A
−
r
,
̄
u
)
+
r
,
̄
A
r
,
̄
B
g
̄
A
̄
B
,
g
rA
=
∂r
∂
̄
x
̄
μ
g
̄
A
̄
μ
=
r
,
̄
λ
g
̄
λ
̄
A
+
r
,
̄
B
g
̄
A
̄
B
.
(37)
From this we can also construct the inverse Jacobian
elements. The elements of that Jacobian we shall need
are
∂
̄
u
∂u
= 1
,
∂
̄
u
∂x
i
= 0
,
∂
̄
λ
∂u
=
−
r
,
̄
u
r
,
̄
λ
,
∂
̄
x
̄
A
∂x
A
=
δ
̄
A
A
.
∂
̄
x
̄
A
∂r
=
∂
̄
x
̄
A
∂u
= 0
.
(38)
The final metric element we shall want is
g
AB
which we
can compute as
g
AB
=
∂
̄
x
̄
α
∂x
A
∂
̄
x
̄
β
∂x
B
g
̄
α
̄
β
=
g
̄
A
̄
B
+
∂
̄
λ
∂x
B
g
̄
λ
̄
A
+
∂
̄
λ
∂x
A
g
̄
λ
̄
B
+
∂
̄
λ
∂x
A
∂
̄
λ
∂x
B
g
̄
λ
̄
λ
=
g
̄
A
̄
B
(39)
where we made use of the fact that
g
̄
λ
̄
λ
=
g
̄
λ
̄
A
= 0.
Because
u
and
x
A
are equal to
̆
t
and ̆
x
̆
A
on the world
tube and are constant along outgoing null geodesics, the
time and angular coordinates (
̆
t,
̆
x
̆
A
) on the world tube
determine the coordinates
u
and
x
A
throughout the char-
acteristic region, including on
I
+
. Thus, the coordinates
at
I
+
will be gauge-dependent, since
̆
t
and ̄
x
̄
A
are de-
pendent upon the gauge choices made in the 3+1 Cauchy
evolution. We will later eliminate this gauge dependence
by evolving and transforming to the coordinates of free-
falling observers on
I
+
, as described below in Sec. V B.
C. Inner Boundary Values of Characteristic
Variables
Now that we have the full metric in Bondi-Sachs
form (Eq. (1)), we assemble the inner boundary values
for the various evolution variables used in the volume,
J,β,Q,U,W,
and
H
. We write out the complex dyads as
q
A
=
{−
1
,
−
i
sin
θ
}
,
q
A
=
{
−
1
,
−
i
sin
θ
}
.
(40)
Because of the identification between the intermediate
angular coordinates ̄
x
̄
A
and the Characteristic coordi-
nates
x
A
, the dyads are identified,
q
A
=
q
̄
A
and
q
A
=
q
̄
A
.
Then, as a consequence of Eq. (40),
q
A,
̄
λ
=
q
A
,
̄
λ
= 0 and
q
A,
̄
u
=
q
A
,
̄
u
= 0.
Inverting the metric in Eq. (1),
g
μν
=
0
−
e
−
2
β
0
A
−
e
−
2
β
(1 +
rW
)
e
−
2
β
−
e
−
2
β
U
A
0
B
−
e
−
2
β
U
B
r
−
2
h
AB
,
(41)
where
h
AB
h
BC
=
δ
C
A
and
|
h
AB
|
=
|
q
AB
|
.
In the PittNull code, the quantities
J
,
β
,
Q
,
U
, and
W
and their
̄
λ
derivatives are computed using an expan-
sion in affine coordinates to compute their values along a
surface of constant surface area coordinate
r
[24]. Pitt-
Null then chooses its internal compactified radinull co-
ordinates in the characteristic region to be surfaces of
constant
r
. However, in Ref. [31] and here, we choose
our inner boundary to be the worldtube. The value of
the surface area coordinate
r
at the world tube we define
as
R
(
u,x
A
),
R
=
r
|
Γ
,
(42)
R
,
̄
λ
=
r
,
̄
λ
|
Γ
,
(43)
R
,
̄
u
=
r
,
̄
u
|
Γ
.
(44)
The consequences of this change in inner boundary hy-
persurface are discussed in more detail within Sec. IV A.
We can now write down the inner boundary values of
the characteristic variables in terms of the metric coef-
ficients that we have computed at the inner boundary.
Going back to the definition of
J
=
1
2
q
A
q
B
h
AB
, we get
the expressions
J
|
Γ
=
1
2
R
2
q
A
q
B
g
AB
=
1
2
R
2
q
̄
A
q
̄
B
g
̄
A
̄
B
,
(45)
K
|
Γ
=
√
1 +
J
|
Γ
̄
J
|
Γ
,
(46)
J
,
̄
λ
|
Γ
=
1
2
R
2
q
̄
A
q
̄
B
g
̄
A
̄
B,
̄
λ
−
2
R
,
̄
λ
R
J
|
Γ
,
(47)
J
,
̄
u
|
Γ
=
1
2
R
2
q
̄
A
q
̄
B
g
̄
A
̄
B,
̄
u
−
2
R
,
̄
u
R
J
|
Γ
.
(48)
To get the inner boundary value of
H
, we expand
J
,u
as
J
,u
=
∂
̄
u
∂u
J
,
̄
u
+
∂
̄
λ
∂u
J
,
̄
λ
(49)
so then we find after substituting and simplifying that,
H
|
Γ
=
1
2
R
2
q
̄
A
q
̄
B
(
g
̄
A
̄
B,
̄
u
−
R
,
̄
u
R
,
̄
λ
g
̄
A
̄
B,
̄
λ
)
.
(50)
We can read off the value for
g
ur
to compute
β
,
β
|
Γ
=
−
1
2
ln
(
R
,
̄
λ
)
.
(51)
We will also need
β
,
̄
λ
|
Γ
in order to compute
Q
|
Γ
. Directly
differentiating Eq. (51) yields
β
,
̄
λ
|
Γ
=
−
R
,
̄
λ
̄
λ
2
R
,
̄
λ
,
(52)
7
but this involves the quantity
R
,
̄
λ
̄
λ
, which appears to
depend on second derivatives of the metric. So we instead
compute
β
,
̄
λ
|
Γ
using
β
’s evolution equation, Eq. (B1):
β
,
̄
λ
|
Γ
=
R
8
R
,
̄
λ
(
J
,
̄
λ
|
Γ
̄
J
,
̄
λ
|
Γ
−
(
K
,
̄
λ
|
Γ
)
2
)
,
(53)
which involves only first derivatives.
The quantities
U
and
W
can similarly be read off from
the metric:
U
|
Γ
=
g
rA
g
ur
q
A
,
(54)
W
|
Γ
=
1
R
(
−
g
rr
g
ur
−
1
)
.
(55)
To get
Q
|
Γ
, we will also need
U
,
̄
λ
|
Γ
, which we com-
pute by differentiating the expression for
U
|
Γ
and using
Eq. (52) to eliminate
R
,
̄
λ
̄
λ
in favor of
β
,
̄
λ
|
Γ
:
U
,
̄
λ
|
Γ
=
−
(
g
̄
λ
̄
A
,
̄
λ
+
R
,
̄
λ
̄
B
R
,
̄
λ
g
̄
A
̄
B
+
R
,
̄
B
R
,
̄
λ
g
̄
A
̄
B
,
̄
λ
)
q
̄
A
+2
β
,
̄
λ
|
Γ
(
U
|
Γ
+
g
̄
λ
̄
A
q
̄
A
)
,
(56)
where it is understood that
β
,
̄
λ
|
Γ
is to be evaluated using
Eq. (53). Now that we have an expression for
U
,
̄
λ
|
Γ
, the
inner boundary value of
Q
is given by
Q
|
Γ
=
R
2
(
J
|
Γ
̄
U
,
̄
λ
|
Γ
+
K
|
Γ
U
,
̄
λ
|
Γ
)
.
(57)
D. Computational Domain
We implement angular basis functions through the
use of the external code packages
Spherepack
[41, 42],
which can handle standard spherical harmonics, and
Spinsfast
[43], which is capable of handling SWSHes.
The world tube metric and most of the intermediate
quantities of the inner boundary formalism are real, ten-
sorial metric quantities (i.e. representable by the typi-
cal spherical harmonics), so we use
Spherepack
. Once
all of the inner boundary values of the Bondi evolution
quantities are computed, they are then projected onto
the basis utilized by
Spinsfast
for use during the vol-
ume evolution. Because Cauchy codes evaluate the world
tube data at discrete time slices, we use cubic interpola-
tion to evaluate each of the metric quantities at arbitrary
time values.
IV. VOLUME EVOLUTION
A. Computational Domain
Because the domain of characteristic evolution extends
all of the way out to
I
+
where the surface area coordi-
nate
r
is infinite, to express
I
+
on a finite computational
domain, we define a compactified coordinate,
ρ
,
ρ
=
r
R
+
r
(58)
where
R
is the surface area coordinate of the world tube
given in Eq. (42) so that
ρ
runs from
ρ
|
Γ
= 1
/
2 to
ρ
|
I
+
= 1. This choice of compactification is subtly dif-
ferent from that which is used in PittNull [27]. Because
they expand in affine coordinates to obtain a hypersur-
face of constant Bondi radius, their compactification pa-
rameter is constant and unchanging during their entire
evolution. By tying our compactification parameter to a
fixed Cauchy coordinate radius ̆
r
and allowing the surface
area coordinate
r
to change freely, we must be careful in
how we define our derivatives.
One consequence of utilizing
ρ
is that angular deriva-
tives computed numerically on our grid,
ð
|
ρ
, are evalu-
ated at a constant value of
ρ
, so these are not the same
as angular derivatives defined on surfaces of constant
r
,
which we denote as
ð
. Since Eqs. (14)–(18) involve
ð
and
not
ð
|
ρ
, we must apply a correction factor to compute
ð
from
ð
|
ρ
:
ð
F
=
ð
|
ρ
F
−
F
,ρ
ð
|
ρ
ρ
=
ð
|
ρ
F
−
F
,ρ
ρ
(1
−
ρ
)
R
ð
|
ρ
R,
(59)
for an arbitrary spin-weighted scalar quantity
F
. Similar
correction factors are needed for second derivatives that
appear in the evolution equations:
(
ð
F
)
,ρ
=
ð
|
ρ
F
,ρ
−
F
,ρ
1
−
2
ρ
R
ð
|
ρ
R
−
F
,ρρ
ρ
(1
−
ρ
)
R
ð
|
ρ
R,
(60)
̄
ðð
F
=
̄
ð
|
ρ
ð
|
ρ
F
+
F
,ρ
(
ρ
(1
−
ρ
)
R
2
)
(
2(1
−
ρ
)
̄
ð
|
ρ
R
ð
|
ρ
R
−
R
̄
ð
|
ρ
ð
|
ρ
R
)
−
ð
|
ρ
F
,ρ
(
ρ
(1
−
ρ
)
R
̄
ð
|
ρ
R
)
−
̄
ð
|
ρ
F
,ρ
(
ρ
(1
−
ρ
)
R
ð
|
ρ
R
)
+
F
,ρρ
(
ρ
(1
−
ρ
)
R
)
2
̄
ð
|
ρ
R
ð
|
ρ
R.
(61)
Correction factors for
̄
ð
F
,
̄
ð
F
,ρ
,
ðð
F
,
ð
̄
ð
F
, and
̄
ð
̄
ð
F
are obtained by appropriately interchanging
ð
and
̄
ð
in
8
Eqs. (59)–(61).
Numerical derivatives with respect to
t
and
u
are also
taken at constant
ρ
on our grid, but at constant
r
in the
equations, so similar correction factors are required there
as well, as discussed below in Sec. IV E.
We employ computational grid meshes suitable for
spectral methods, Chebyshev-Gauss-Lobatto for the rad-
inull direction and
Spinsfast
mesh for the angular di-
rections with a uniform
φ
and
θ
grids.
B. Spectral representability
Spectral techniques represent functions over a finite
numerical domain as a series of polynomial functions.
Such representations are of greatest use when the nu-
merical evolution gives rise to smooth solutions, which
converge exponentially with resolution in the spectral ex-
pansion. However, any defect in the solution, such as dis-
continuities, corners, cusps, or the presence of logarith-
mic dependence, will spoil the exponential convergence
of a spectral method, and potentially introduce spurious
oscillatory contributions to the numerical result. For this
reason, it is of great importance to the characteristic evo-
lution code in SpEC to minimize or eliminate sources of
such non-regular contributions to the hypersurface equa-
tions.
The nature of the characteristic hypersurface equations
permits terms proportional to log(
r
) to develop in the
solution of the characteristic evolution system. These
terms are not representable by polynomial expansions
in 1
/r
or by polynomial expansions in
ρ
, so if present
they spoil exponential convergence. Such terms creep
into the evolved solutions by three principal avenues: 1)
via the initial data choice, which if constructed naively,
can excite logarithmic modes, 2) via poorly-chosen co-
ordinates of the metric on the
u
=const hypersurfaces,
and 3) via incomplete numerical cancellation in the equa-
tions, which possess nontrivial pole structure. Points 1)
and 2) arise from the use of the asymptotically non-flat
Bondi form of the spacetime metric, Eq. (1). In that
form, even mathematically faithful solutions to the hy-
persurface equations for generic worldtube data possess
logarithmic dependence. These logarithmic terms would
vanish in an asymptotically flat coordinate system, so are
a
pure gauge
contribution.
In sections IV C and IV D, we explain our methods
for minimizing the logarithmic contributions in the char-
acteristic evolution system. As part of the discussion
in Sec. IV C, we describe the choice of initial data that
eliminates logarithmic dependence from the first hyper-
surface of the characteristic evolution system, address-
ing point 1) above. In Sec. IV D, we describe improve-
ments to the integration techniques that address point
3) above. These methods reduce logarithmic dependence
to the point where it is not noticable in the tests pre-
sented here. However, the full remedy for point 2) re-
quires careful re-examination of the characteristic evo-
lution equations and a set of coordinate transformations
for the evolution system that will be considered for future
development of spectral characteristic techniques, but is
beyond the scope of this paper.
C. Initial Data Slice
The characteristic evolution equations require bound-
ary data on two boundaries: the worldtube (thick red
curve in Fig. 1) and an initial slice
u
=
u
0
(thick blue
curve in Fig. 1). Boundary values on the worldtube were
treated in Sec. III above; here we discuss values on the
initial slice. Given the hierarchical nature of the evolu-
tion equations, the only piece of the metric we need to
specify on the initial slice is
J
, as we can compute all
of the other evolution quantities from
J
using Eqs. (14)–
(18). The main mathematical consideration for choosing
J
for the initial slice is ensuring the regularity of
J
at
I
+
;
the main physical consideration in typical applications is
choosing a
J
that corresponds to no incoming radiation,
either by a linearized approximation [26], or by matching
to a linearized solution [44]. Finally, there is the numeri-
cal consideration mentioned in Sec. IV B that we wish to
minimize the excitation of pure-gauge logarithmic depen-
dence and keep the initial data
C
∞
over the numerical
domain.
When choosing
J
on the initial
u
=
u
0
slice, we wish to
match the worldtube data provided by the Cauchy code
as closely as possible. The worldtube data that we take
as input (see Sec. III) consist of the full spacetime metric
and its first radial and time derivatives, which are suffi-
cient to constrain the value of
J
and the value of
∂
r
J
on
the worldtube. By careful analysis of the characteristic
evolution equations, one can show that the initial
u
=
u
0
hypersurface is free of logarithmic dependence if [45]
∂
2
`
J
−
J
(
(
∂
`
K
)
2
−
∂
`
J∂
`
̄
J
)
= 0 at
I
+
. This condition
is satisfied by the simpler conditions
J
=
J
,``
= 0 at
I
+
,
so we construct an initial
J
that satisfies
J
=
J
,``
= 0 at
I
+
and matches the world tube data. This construction
is consistent with the input Cauchy data in the overlap
region of (Fig. 1) to linear order in a radial expansion.
Our initial choice of
J
, determined by the functions
J
|
Γ
and
∂
r
J
|
Γ
, is
J
initial
=
R
2
r
(3
J
|
Γ
+
R∂
r
J
|
Γ
)
−
R
3
2
r
3
(
J
|
Γ
+
R∂
r
J
|
Γ
)
=
R
2
(
1
ρ
−
1
)
(3
J
|
Γ
+
R∂
r
J
|
Γ
)
−
R
3
2
(
1
ρ
−
1
)
3
(
J
|
Γ
+
R∂
r
J
|
Γ
)
.
(62)
D. Radinull Integration
The characteristic equations Eqs. (14)–(18) can be
solved in sequence by integration in
r
from the worldtube
9
to
I
+
. We use a numerical radinull grid in the compacti-
fied variable
ρ
, and we re-express the characteristic equa-
tions in terms of
ρ
derivatives; see Eqs. (B1)–(B6). The
grid points in
ρ
are chosen at Chebyshev-Gauss-Lobatto
quadrature points. The radinull equations for
β
,ρ
and
U
,ρ
(Eq. (B1), (B3)) both lend themselves to straight-
forward Chebyshev-Gauss-Lobatto quadrature. Starting
at the inner boundary values of
β
|
Γ
(Eq. (51)) and
U
|
Γ
(Eq. (55)), these evolution variables are integrated out
to
I
+
.
A quick examination of the radinull equations for the
evolution quantities
Q
,ρ
,W
,ρ
,
and
H
,ρ
(Eq. (B2), (B5),
(B6)) reveals powers of (
ρ
−
1) in denominators, so regu-
larity at
I
+
(
ρ
= 1) is not guaranteed by the form of
the equations. A previous version of this same spec-
tral characteristic evolution method [31] utilized inte-
gration by parts in order to rewrite the equations in a
form without poles, allowing them to be integrated di-
rectly via Chebyshev-Gauss-Lobatto quadrature. How-
ever, integration by parts introduced logarithmic terms
like log(1
−
ρ
) which canceled analytically in the final
results of gauge-invariants such as the Bondi news, but
which were not well represented by a Chebyshev-Gauss-
Lobatto spectral expansion in
ρ
. These logarithmic terms
spoiled exponential convergence and led to a large noise
floor, limiting the accuracy of the method. We choose an
alternative approach here.
The evolution equation for
Q
, Eq. (B2), can be written
in the form
(
r
2
Q
)
,ρ
=
Q
C
(1
−
ρ
)
2
+
Q
D
(1
−
ρ
)
3
,
(63)
where
Q
C
corresponds to the 1
/
(1
−
ρ
)
2
term and
Q
D
is
the 1
/
(1
−
ρ
)
3
term in Eq. (B2), and all factors of (1
−
ρ
)
in denominators have been written explicitly.
To better characterize the asymptotic behavior of this
equation, we rewrite the system in terms of the inverse
radinull coordinate
x
=
R/r
= 1
/ρ
−
1. Then Eq. (63)
becomes
(
Q
x
2
)
,x
=
C
x
2
+
D
x
3
,
(64)
where
C
=
−
Q
C
+
Q
D
R
2
,
(65)
D
=
−
Q
D
R
2
.
(66)
We know the right-hand side quantities
C
and
D
are
regular at
x
= 0, and we seek a solution
Q
that is also
regular there. So we introduce new variables, motivated
by Taylor series expansions of
Q
,
C
, and
D
about
I
+
(
x
= 0),
Q
=
Q
−
Q
|
I
+
−
xQ
,x
|
I
+
,
(67)
C
=
C
−
C
|
I
+
−
xC
,x
|
I
+
,
(68)
D
=
D
−
D
|
I
+
−
xD
,x
|
I
+
−
x
2
2
D
,xx
|
I
+
.
(69)
Thus, by construction,
Q
and
C
are both guaranteed to
behave like
x
2
near
x
= 0 while
D
behaves as
x
3
. Sub-
stituting these variables into Eq. (64) and gathering like
terms, we find the differential equation
(
Q
x
2
)
,x
=
C
x
2
+
D
x
3
+
2
C
,x
|
I
+
+
D
,xx
|
I
+
2
x
+
Q
,x
|
I
+
+
C
|
I
+
+
D
,x
|
I
+
x
2
+
2
Q
|
I
+
+
D
|
I
+
x
3
.
(70)
Because of how we have defined
Q
,
C
and
D
, any potential
singularity issues are confined to the last three terms. To
satisfy Eq. (70) for all
x
, the numerators of each of these
terms must identically vanish, providing constraints and
boundary conditions on the asymptotic values of
Q
,
C
,
and
D
,
Q
|
I
+
=
−
D
|
I
+
2
,
(71)
Q
,x
|
I
+
=
−
C
|
I
+
−
D
,x
|
I
+
,
(72)
0 =
−
C
,x
|
I
+
−
1
2
D
,xx
|
I
+
.
(73)
The last equation, Eq. (73) is a regularity condition on
C
and
D
. If satisfied, it ensures no logarithmic depen-
dence in the solution to the
Q
equation. A careful anal-
ysis of the differential equations, which will be presented
in complete detail in future work, shows that the lead-
ing violation of Eq. (73) is
∝
̄
ð
∂
2
x
J
|
I
+
, and that Eq.
(73) is entirely satisfied if
J
= 0 and
J
,xx
= 0 at
I
+
.
The leading violation of the conditions on
J
can be de-
termined through further analysis to have the leading
contribution of
U
(
∂
x
J
)
2
|
I
+
. These pure-gauge regular-
ity violations are important to note for precision studies
and for unusual regimes for characteristic evolution, but
for the practical evolutions, the scales we observe do not
typically exceed
U
∼
10
−
6
,
J
∼
10
−
3
. So, even for long
evolutions, the logarithmic dependence does not grow to
a significant fraction the main contribution.
We now integrate the equation
(
Q
x
2
)
,x
=
C
x
2
+
D
x
3
(74)
with inner boundary value
Q
|
Γ
=
Q
|
Γ
+
D
|
I
+
2
+
(
C
|
I
+
+
D
,x
|
I
+
)
(75)
to obtain
Q
at all radinull points. Then we reconstruct
Q
by adding back in its asymptotic values,
Q
=
Q−
D
|
I
+
2
−
x
(
C
|
I
+
+
D
,x
|
I
+
)
.
(76)
Because the equation for
Q
does not mix the real and
imaginary parts of
Q
, we follow [31] and solve for real
and imaginary parts of
Q
separately.
10
Examining the evolution equation for
W
, Eq. (B5), we
recognize that it has the same form as the equation for
Q
,
Eq. (B2). Therefore, in order to solve for
W
, we use the
same procedure as we do for
Q
, following from Eq. (63)
through Eq. (76) but replacing all of quantities specific
to
Q
with their
W
equivalents.
The radinull equation for
H
, Eq. (B6) can be written
as
(
rH
)
,ρ
−
rJ
2
(
H
̄
T
+
̄
HT
) =
H
A
+
H
B
1
−
ρ
+
H
C
(1
−
ρ
)
2
,
(77)
where
H
B
= Σ
i
H
Bi
. The form of this equation is very
similar to that of Eq. (63) that governs the
Q
(and
W
)
radinull evolution. However, there is now the additional
complication that
H
,ρ
has a term proportional to not only
H
, but also to
̄
H
. This couples the real and imaginary
parts of the equation.
The previous version of this code employed the Magnus
expansion in order to handle this difficulty [31]. While
the Magnus expansion might be useful for systems where
the terms in its expansion are rapidly shrinking, there
is no guarantee that will hold in general. Instead, we
will write the system as a matrix differential equation,
expressing
H
(and
H
A
,H
B
,
and
H
C
) as column vectors
like
H
=
(
<
(
H
)
=
(
H
)
)
,
(78)
and defining the quantity
M
as
M
≡
(
<
(
J
)
<
(
T
)
<
(
J
)
=
(
T
)
=
(
J
)
<
(
T
)
=
(
J
)
=
(
T
)
)
,
(79)
so that
MH
here represents matrix multiplication. Then
Eq.(77) becomes the matrix equation,
(
rH
)
,ρ
−
rMH
=
H
A
+
H
B
1
−
ρ
+
H
C
(1
−
ρ
)
2
,
(80)
As before, we convert from
ρ
into the inverse radinull
coordinate
x
=
R/r
= 1
/ρ
−
1 to better characterize its
behavior near
I
+
,
(
H
x
)
,x
+
M
H
x
=
A
+
B
x
+
C
x
2
(81)
where
M
=
M
(1 +
x
)
2
,
(82)
A
=
−
H
A
R
(1 +
x
)
2
,
(83)
B
=
−
H
B
R
(1 +
x
)
,
(84)
C
=
−
H
C
R
.
(85)
As we did with the
Q
equation, we shall introduce one
final set of variables, motivated by Taylor series expan-
sions of
H,B,
and
C
about
x
= 0.
H
=
H
−
H
|
I
+
,
(86)
B
=
B
−
B
|
I
+
−M
H
|
I
+
+
M
|
I
+
H
|
I
+
,
(87)
C
=
C
−
C
|
I
+
−
xC
,x
|
I
+
.
(88)
Once again, these variables are constructed so that
H
and
B
behave as
x
and
C
behaves as
x
2
in a neighborhood
about
x
= 0. Substituting these into Eq. (81), we get
(
H
x
)
,x
+
M
H
x
=
A
+
B
x
+
C
x
2
+
H
|
I
+
+
C
|
I
+
x
2
+
B
|
I
+
+
C
,x
|
I
+
−M
|
I
+
H
|
I
+
x
(89)
As before, the numerators of the last two terms must
vanish, which gives us a boundary condition on
H
at
I
+
,
H
|
I
+
=
−
C
|
I
+
,
(90)
and a boundary constraint on
B
,
C
, and
M
,
0 =
B
|
I
+
+
C
,x
|
I
+
+
M
|
I
+
C
|
I
+
.
(91)
The last constraint is a regularity condition that is guar-
anteed to be satisfied provided the input spin-weighted
scalars
β,Q,U,
and
W
themselves are regular [45]. Of
course, the small violation that arises from the
Q
and
W
equations will lead to a similarly small violation in the
regularity of
H
. In principle, a carefully chosen coordi-
nate transformation could fully address all of these small
violations.
We then integrate the equation
(
H
x
)
,x
+
M
H
x
=
A
+
B
x
+
C
x
2
(92)
from the worldtube to
I
+
, with boundary value
H
|
Γ
=
H
|
Γ
+
C
|
I
+
, to obtain
H
on the entire null slice. We
reconstruct
H
by computing
H
=
H−
C
|
I
+
.
(93)
To help ensure the stability of the system, we perform
spectral filtering for each of the evolution quantities
J
,
β
,
Q
,
U
,
W
, and
H
after every time we compute them,
similar to [31]. For the angular filtering, we set to 0
the highest two
`
−
modes in the spectral decomposition
on each shell of constant
ρ
. Thus, resolving the system
up through
`
max
modes requires storing and evolving the
evolution quantities in the volume up through
`
=
`
max
+
2 modes. We filter along the radinull direction by taking
the spectral expansion of the evolution quantities along
each null ray and scaling the
i
−
th coefficient by
e
−
108(
i/
(
n
ρ
−
1))
16
,
(94)