Testing outer boundary treatments for the Einstein equations
Testing outer boundary treatments for the Einstein
equations
Oliver Rinne, Lee Lindblom and Mark A. Scheel
Theoretical Astrophysics 130-33, California Institute of Technology, Pasadena,
CA 91125, USA
Abstract. Various methods of treating outer boundaries in numerical relativity
are compared using a simple test problem: a Schwarzschild black hole with
an outgoing gravitational wave perturbation. Numerical solutions computed
using different boundary treatments are compared to a ‘reference’ numerical
solution obtained by placing the outer boundary at a very large radius. For
each boundary treatment, the full solutions including constraint violations and
extracted gravitational waves are compared to those of the reference solution,
thereby assessing the reflections caused by the artificial boundary. These tests are
based on a first-order generalized harmonic formulation of the Einstein equations
and are implemented using a pseudo-spectral collocation method. Constraint-
preserving boundary conditions for this system are reviewed, and an improved
boundary condition on the gauge degrees of freedom is presented. Alternate
boundary conditions evaluated here include freezing the incoming characteristic
fields, Sommerfeld boundary conditions, and the constraint-preserving boundary
conditions of Kreiss and Winicour. Rather different approaches to boundary
treatments, such as sponge layers and spatial compactification, are also tested.
Overall the best treatment found here combines boundary conditions that
preserve the constraints, freeze the Newman-Penrose scalar Ψ0, and control gauge
reflections.
PACS numbers: 04.25.Dm, 02.60.Lj, 02.60.Cb
1. Introduction
A fundamental problem in numerical relativity is the need to solve Einstein’s equations
on spatially unbounded domains with finite computer resources. There are various
ways of addressing this issue. Most often, the spatial domain is truncated at a finite
distance and suitable boundary conditions are imposed at the artificial boundary.
A different approach is to compactify the domain by using spatial coordinates that
bring spatial infinity to a finite location on the computational grid. Another method
often used for wave-like problems (although it is not commonly used in numerical
relativity) includes so-called sponge layers which damp the waves near the outer
boundary of the computational domain. The purpose of this paper is to compare these
various methods by testing their ability to accurately reproduce dynamical solutions
of Einstein’s equations.
An ideal boundary treatment would produce a solution to Einstein’s equations
that is identical (within the computational domain) to the corresponding solution
obtained on an unbounded domain. In particular, no spurious gravitational radiation
or constraint violations should enter the computational domain through the artificial
Testing outer boundary treatments for the Einstein equations 2
boundary. We can use this principle to test the various boundary treatments in the
following way. First we compute a reference solution using a very large computational
domain, large enough that its boundary remains out of causal contact with the interior
spacetime region where comparisons are being made. Next we compute the same
solution using a domain truncated at a smaller distance where one of the boundary
treatments is used: we either impose boundary conditions there, compactify spatial
infinity, or add a sponge layer. Finally we compare the solution on the smaller domain
with the reference solution, measuring the reflections and constraint violations caused
by the boundary treatment. Assessing boundary conditions by comparing with a
reference solution on a much larger domain or a known analytic solution is a common
practice in computational science. For applications to numerical relativity see e.g. [1],
chapter 8 of [2], and [3, 4, 5].
The particular test problem used in this paper is a Schwarzschild black hole
with an outgoing gravitational wave perturbation. The interior of the black hole is
excised; all the characteristic fields propagate into the black hole (and out of the
computational domain) at the inner boundary and hence no boundary conditions
are needed there. Our numerical implementation uses a pseudo-spectral collocation
method. See Appendix A for details on the initial data, the numerical methods, and
the quantities that we compare between the solutions.
We perform all of these tests using a first-order generalized harmonic formulation
of the Einstein equations (see [6] and references therein). In section 2 we discuss
the construction of boundary conditions for this system that prevent the influx of
constraint violations, and that limit the spurious incoming gravitational radiation
by controlling the Newman-Penrose scalar Ψ0 at the boundary. We also improve
the boundary conditions on the gauge degrees of freedom by studying small gauge
perturbations of flat spacetime. We then evaluate the performance of these boundary
conditions on our test problem: measuring the reflections and constraint violations
caused by the computational boundary, and determining how these reflections vary
with the radius of the boundary.
Section 3 evaluates the performance of a variety of other widely used boundary
conditions on our test problem. First we test the simple boundary conditions that
freeze all the incoming characteristic fields at the boundary. We also test the
commonly used variant of this, the Sommerfeld boundary conditions, used in many
binary black hole simulations [7, 8, 9, 10, 11] based on the BSSN [12, 13] formulation
of Einstein’s equations. Finally in section 3 we evaluate the constraint-preserving
boundary conditions proposed by Kreiss and Winicour [14], which differ from those
discussed in section 2 mainly by our use of a physical boundary condition that controls
In section 4 we evaluate two boundary treatments that are alternatives to
imposing local boundary conditions at a finite outer boundary. The first is the spatial
compactification method used e.g. by Pretorius [15, 16, 17] in his groundbreaking
binary black hole evolutions. In this treatment a coordinate transformation maps
spatial infinity to a finite location on the computational grid. As waves travel out,
they become increasingly blue-shifted with respect to the compactified coordinates
and ultimately they fail to be resolved. Hence numerical dissipation is applied, which
damps away these short-wavelength features. We measure the reflections and the
constraint violations generated by the waves in our test problem as they interact with
this boundary treatment. Finally in section 4 we implement and test a sponge layer
method for Einstein’s equations.
Testing outer boundary treatments for the Einstein equations 3
One of the main objectives of current binary black hole simulations is the
computation of reliable waveforms for gravitational wave data analysis. Therefore it is
important to evaluate how the various boundary treatments affect the accuracy of the
extracted waveforms. In section 5, we compute the Newman-Penrose scalar Ψ4 (which
describes the outgoing waves) on an extraction sphere close to the outer boundary (or
compactified region, or sponge layer, respectively) and compare it with the analogous
Ψ4 from the reference solution. We also compare the measured reflections caused
by our Ψ0 controlling boundary condition with the analytical predictions of these
reflections made by Buchman and Sarbach [18, 19].
Finally we discuss the implications of our results in section 6, and we also describe
briefly a number of other boundary treatments which we do not test here.
2. Constraint-preserving boundary conditions
In this section, we briefly review the generalized harmonic form of the Einstein
evolution system used in our tests. The method of constructing constraint-preserving
boundary conditions (CPBCs) for this system is also discussed, and an improved
boundary condition for the gauge degrees of freedom is derived. The numerical
performance of these boundary conditions is evaluated using our test problem, and
the dependence of the spurious reflections as a function of the boundary radius is
measured.
2.1. The generalized harmonic evolution system
The formulation of Einstein’s equations employed here uses generalized harmonic
gauge conditions, in which the coordinates xa obey the wave equation
�xa = Ha(x, ψ), (1)
where � = ψab(∂a∂b − Γcab∂c) is the covariant scalar wave operator, with ψab the
spacetime metric and Γcab the associated metric connection. In this formulation of
the Einstein system the gauge source function Ha may be chosen freely as a function
of the coordinates and of the spacetime metric ψab (but not derivatives of ψab).
As is well known, the Einstein equations reduce to a set of coupled wave equations
when the gauge is specified by equation (1). We write this system in first-order form,
both in time and space, by introducing the additional variables Φiab ≡ ∂iψab and
Πab ≡ −tc∂cψab, where tc is the future directed unit normal to the t = const.
hypersurfaces. Here lower-case Latin indices from the beginning of the alphabet
denote four-dimensional spacetime quantities, whereas lower-case Latin indices from
the middle of the alphabet are spatial. The principal parts of these evolution equations
are given by ‡
∂tψab ' 0,
∂tΠab ' Nk∂kΠab −Ngki∂kΦiab − γ2Nk∂kψab, (2)
∂tΦiab ' Nk∂kΦiab −N∂iΠab +Nγ2∂iψab,
where ' indicates that purely algebraic terms have been omitted, gij is the spatial
metric of the t = const. slices, and N and N i are the lapse function and shift vector,
‡ The parameter γ1 of [6] is chosen to be −1, which ensures that the equations are linearly degenerate.
Testing outer boundary treatments for the Einstein equations 4
respectively. The parameter γ2 was introduced in [6] in order to damp violations of
the three-index constraint
Ciab ≡ ∂iψab − Φiab = 0. (3)
We also include terms of lower derivative order that are designed to damp violations
of the harmonic gauge constraint [20]
Ca ≡ −�xa +Ha = ψbcΓabc +Ha = 0. (4)
The system (2) is symmetric hyperbolic. The characteristic fields in the direction
ni (where nata = 0) are given by
u0ab = ψab, speed 0, (5)
u1±ab = Πab ± Φnab − γ2ψab, speed −N
n ±N, (6)
u2Aab = ΦAab, speed −Nn. (7)
For future reference, we also define
ũ1±ab ≡ Πab ± Φnab. (8)
Here and in the following, an index n denotes contraction with ni, while upper-
case Latin indices A,B, . . . are orthogonal to n, e.g. vA = PAivi where Pab ≡
ψab − nanb + tatb. For further details, we refer the interested reader to [6].
2.2. Construction of boundary conditions
Our construction of boundary conditions for the generalized harmonic evolution
system can be divided into three parts: constraint-preserving, physical, and gauge
boundary conditions.
In order to impose constraint-preserving boundary conditions, we derive the
subsidiary evolution system that the constraints (3) and (4) obey as a consequence of
the main evolution equations (2). The incoming modes of the subsidiary system are
then required to vanish at the boundary (cf. [21, 22, 23, 24, 25, 26, 27, 28, 29]). For
instance, the harmonic gauge constraint (4) obeys a wave equation
�Ca = (lower-order terms homogeneous in the constraints) (9)
and the corresponding incoming fields will involve first derivatives of Ca. In terms of
the incoming modes u1−ab (6) of the main evolution equations, the resulting constraint-
preserving boundary conditions can be written in the form
PC cdab ∂nu
cd ≡ ( 12PabP
cd − 2l(aPb)(ckd) + lalbkckd)∂nu1−cd
= (tangential derivatives), (10)
where PC is a projection operator of rank 4 (cf. [6]). Here ni now refers to the outward-
pointing unit spatial normal to the boundary, la = (ta + na)/
2, ka = (ta − na)/
= denotes equality at the boundary. If the shift vector points towards the exterior
at the boundary (Nn >̇ 0), the fields u2Aab (7) are incoming as well and we obtain a
boundary condition on them by requiring the components CnAab of the four-index
constraint
Cijab ≡ −2∂[iΦj]ab (11)
to vanish at the boundary.
An acceptable physical boundary condition should require that no gravitational
radiation enter the computational domain from the outside (except for backscatter
Testing outer boundary treatments for the Einstein equations 5
off the spacetime curvature, an effect that is a first-order correction in M/R).
Gravitational radiation may be described by the evolution system that the Weyl tensor
obeys by virtue of the Bianchi identities (see e.g. [27]). Our boundary condition
requires the incoming characteristic fields of this system to vanish at the outer
boundary. These incoming fields are proportional to the Newman-Penrose scalar Ψ0
(evaluated for a Newman-Penrose null tetrad containing the vectors la and ka). Hence
the physical boundary condition we use is [27, 22, 30, 29, 31]
= 0, (12)
which can be written in a form similar to (10),
PP cdab ∂nu
cd ≡ (Pa
d − 1
cd)∂nu
= (tangential derivatives). (13)
Here PP is a projection operator of rank 2 that is orthogonal to PC [6]. We remark that
(12) still causes some, albeit very small, spurious reflections of gravitational radiation.
It can be viewed as the lowest level in a hierarchy of perfectly absorbing boundary
conditions for linearized gravity [18, 19].
The constraint-preserving (10) and physical (13) boundary conditions together
constrain six components of the main incoming fields u1−ab . The remaining four
components correspond to gauge degrees of freedom. In the past we chose simply
to freeze those components in time [6],
PG cdab ∂tu
= 0, (14)
where PG ≡ I− PC − PP.
The initial-boundary value problem (IBVP) for the boundary conditions discussed
so far was shown in [32] to be boundary-stable, which is a (rather strong) necessary
condition for well posedness. These boundary conditions have been successfully used
in long-term stable evolutions of single and binary black hole spacetimes [6, 33, 34]. In
the following subsection, we present an improvement to the gauge boundary condition
(14) motivated by the evolution of gauge perturbations about flat spacetime.
2.3. Improved gauge boundary condition
Let us assume that near the outer boundary, the spacetime is close to Minkowski space
in standard coordinates (Ha = 0) so that the Einstein equations may be linearized
about that background. This assumption is reasonable because for the dominant
wavenumber of the outgoing pulse (k = 1.6/M) and the boundary radius we typically
consider (R = 41.9M), we have kR � 1 and R � M . Furthermore, we assume that
the outer boundary is a coordinate sphere of radius r = R.
We begin by noting that harmonic gauge does not fix the coordinates completely:
infinitesimal coordinate transformations
xa → xa + ξa (15)
are still allowed provided the displacement vector satisfies the wave equation,
�ξa = 0. (16)
Under such a coordinate transformation, the metric changes by
δψab = −2∂(aξb). (17)
A closer inspection [32] of the projection operator PG in (14) shows that the gauge
boundary conditions control the components laδψab of the perturbations, where
Testing outer boundary treatments for the Einstein equations 6
la ≡ (ta + na)/
2 is the outgoing null vector normal to the boundary. It is
interesting to observe that these components vanish in the ingoing radiation gauge
[35]. However, imposing radiation gauge on the entire spacetime is not possible in
spacetimes containing strong-field regions, which will always generate perturbations
laδψab that propagate into the far field. A reasonable condition to require then is that
these perturbations pass through the boundary without causing strong reflections.
Each Cartesian component of the vector laδψab obeys the scalar wave equation
�ψ = 0. (18)
Solutions to this equation can be written in the form
Ylm(θ, φ)ψl(t, r), (19)
where the Ylm are the standard spherical harmonics and the ψl are linear combinations
of outgoing (+) and incoming (−) solutions
ψ±l (t, r) = r
F±l (r ∓ t), (20)
F±l (x) being arbitrary functions. A boundary condition is needed on ψ that eliminates
the incoming part of these solutions. In [36], a hierarchy of boundary conditions
is constructed that accomplish this task for all l 6 L. This idea was applied to
the evolution of the Weyl curvature in [18] in order to construct improved physical
boundary conditions. For the gauge boundary conditions considered here, we restrict
ourselves to the L = 0 member of the hierarchy, which corresponds to the Sommerfeld
condition §
(∂t + ∂r + r
= 0. (21)
In contrast, our old gauge boundary condition that froze the incoming characteristic
field, as in (14), is given by
(∂t + ∂r + γ2)ψ
= 0, (22)
where γ2 is the constraint damping parameter.
This Sommerfeld boundary condition (21) is much less reflective than the freezing
condition (22). To see this, we consider a solution of the form
ψl = ψ
l + ρlψ
l (23)
with generating functions
F±l (x) = e
±ikx, (24)
where k ∈ R is the wave number. Substituting this solution into the boundary
conditions (21) resp. (22), we solve for the reflection coefficient ρl. Figure 1 shows
|ρl| for a typical range of wave numbers k and outer boundary radii R used for the
numerical tests in this paper. (The dominant wave number of the outgoing pulse is
k ≈ 1.6/M and in most cases, we place the outer boundary at R = 41.9M .) We
see that |ρl| is much smaller (by about 3 orders of magnitude) for the Sommerfeld
condition than for the freezing condition.
§ To avoid confusion, we remark that in [5, 14], the term ‘Sommerfeld condition’ is used in reference
to a condition of the form (∂t +∂r)u
= 0, i.e. without the extra r−1 term due to our polar coordinates.
Testing outer boundary treatments for the Einstein equations 7
0 0.5 1 1.5 2
l = 1
R = 41.9 M
50 100 150 200
R / M
l = 1
k = 1 / M
Figure 1. Predicted reflection coefficients ρl for freezing (dotted) and
Sommerfeld (solid) boundary conditions as functions of wave number k and outer
boundary radius R. The curves for different l are visually indistinguishable in the
freezing case. Note also that ρ0 = 0 for the Sommerfeld condition.
In the notation of the previous subsection, the improved gauge boundary
condition (21) reads (after taking a time derivative),
PG cdab ∂t[u
cd + (γ2 − r
−1)ψcd]
= 0. (25)
We remark that the extra terms in (25) as compared with the old condition (14) are
of lower derivative order, so that the high-frequency stability result of [32] extends
immediately to these modified gauge boundary conditions.
2.4. Numerical results
The numerical tests of the various boundary conditions performed in this paper are
described in some detail in Appendix A. Figure 2 compares the numerical performance
of our new CPBCs (10), (11), (13), (25) with our old ones (10), (11), (13), (14). The
outer boundary is placed at radius R = 41.9M for these particular tests. Shown are
the discrete L∞ and L2 norms of the difference ∆U between the numerical solution
and the reference solution, and also the violations of the constraints C (see Appendix
A.4 for precise definitions of these quantities). The reference solution has an outer
boundary at radius 961.9M and is computed using our old CPBCs; thus for t < 920M
the outer boundary of the reference solution is out of causal contact with the region
where ∆U and C are computed.
In the difference ∆U we see a reflection that originates when the wave reaches the
boundary at t ≈ R and then amplifies as it moves inward in the spherical geometry,
assuming its maximum at t ≈ 2R. This feature is much more prominent in the L∞
norm than in the L2 norm, which is why we display only the L∞ norm in subsequent
plots. The reflection is much smaller (by a factor of ≈ R/M) for the new boundary
conditions as compared with the old ones. Even at later times, the new boundary
conditions result in a smaller ∆U , which in contrast to the old conditions appears to
decrease as resolution is increased.
We would like to stress that ∆U is a coordinate dependent quantity. Hence a
smaller ∆U does not necessarily mean that the boundary treatment is ‘better’ in a
physically meaningful sense. If however the aim is to produce a solution that is as
Testing outer boundary treatments for the Einstein equations 8
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
Figure 2. Old (solid) vs. new (dotted) CPBCs. Four different resolutions are
shown: (Nr, L) = (21, 8), (31, 10), (41, 12), and (51, 14). The outer boundary is
at R = 41.9M .
close to the reference solution in the same coordinates, the choice of gauge boundary
conditions does become important. Gauge reflections can in principle also impair the
numerical accuracy of gauge-invariant quantities because much numerical resolution
is wasted on resolving the gauge reflections. This is particularly the case when the
gauge excitations in question are high-frequency modes such as those produced along
with the so-called ‘junk radiation’ in binary black hole initial data.
There is no discernible difference between the two sets of boundary conditions as
far as constraint violations are concerned, which is what we expect because both of
them are constraint-preserving.
We close this section by investigating the dependence of the reflections on the
radius of the outer boundary (figure 3). The amplitude of the first peak in ||∆U||∞
decreases as the boundary is moved outward, roughly like 1/R. At late times, there
appears to be a power-law growth of that quantity at a rate that increases slightly
with resolution. Inspection of the constraints (also in figure 3) and Ψ4 (figure 10)
suggests that this is a pure gauge effect. This blow-up is completely dominated by
the innermost domain, which contains a long-wavelength feature that is growing in
time. We speculate that this problem might be cured by a more clever choice of gauge
source function close to the black hole horizon.
Testing outer boundary treatments for the Einstein equations 9
0 200 400 600 800 1000
t / M
R/M = 21.9
161.9
201.9
121.9
 = 51, L = 14)
0 200 400 600 800 1000
t / M
 = 51, L = 14)
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
R = 121.9 M
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
R = 121.9 M
Figure 3. New CPBCs at different radii. Top half: all radii at the highest
resolution, bottom half: R = 121.9M at all resolutions. In the top right panel,
curves for all outer boundary radii coincide.
3. Alternate boundary conditions
In this section, we consider several alternate boundary conditions that are often used
in numerical relativity. All of these are local conditions imposed at a finite boundary
radius, then in section 4 we consider some additional non-local boundary treatments.
We run the alternate boundary conditions on our test problem and compare the results
with the CPBCs (using the new gauge boundary condition (25)).
3.1. Freezing the incoming fields
A very simple boundary condition is obtained by freezing in time all the incoming
fields at the boundary, i.e.,
= 0 (and ∂tu
= 0 if Nn >̇ 0). (26)
This boundary condition is attractive from a mathematical point of view because it
is of maximally dissipative type and hence, together with the symmetric hyperbolic
evolution equations (2), yields a strongly well-posed IBVP [37, 38, 39]. However, in
general this boundary condition is not compatible with the constraints.
Testing outer boundary treatments for the Einstein equations 10
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
Figure 4. Freezing (solid) vs. new CPBCs (dotted). Four different resolutions
are shown: (Nr, L) = (21, 8), (31, 10), (41, 12), and (51, 14). For freezing
boundary conditions, both ||∆U|| and C converge to a nonzero function with
increasing resolution. The outer boundary is at R = 41.9M .
The left side of figure 4 demonstrates that freezing boundary conditions cause a
significantly larger (by ≈ 3 orders of magnitude) initial reflection than our CPBCs.
The difference with respect to the reference solution remains large in the subsequent
evolution and unlike for the CPBCs does not decrease with increasing resolution.
Furthermore, the violations of the constraints (right side of figure 4) do not converge
away. This means that a solution to the Einstein equations is not obtained in the
continuum limit.
3.2. Sommerfeld boundary conditions
A boundary condition that is often imposed in conjunction with the BSSN [12, 13]
formulation of the Einstein equations is a Sommerfeld condition on all the components
of the spatial metric gij and extrinsic curvature Kij ,
(∂t + ∂r + r
gij − δij
= 0. (27)
This condition has been used for example in many recent binary black hole
simulations [7, 8, 9, 10, 11]. We cannot impose precisely the conditions (27) in our
simulations because there is no one-to-one relationship between gij and Kij , and the
incoming characteristic fields of our generalized harmonic formulation of Einstein’s
equations. Instead we consider the similar condition
(∂t + ∂r + r
−1)(ψab − ηab)
= 0 (28)
on all the components of the spacetime metric (ηab being the Minkowski metric). A
very similar boundary condition (without the r−1 term) has recently been used in the
generalized harmonic evolutions of [40].
In our formulation, boundary conditions are required not on the spacetime metric
itself but only on certain combinations of its derivatives. By taking a time derivative
of (28) and rewriting in terms of incoming characteristic fields, we obtain
ab + (γ2 − r
−1)ψab]
= 0. (29)
Testing outer boundary treatments for the Einstein equations 11
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
Figure 5. Sommerfeld (solid) vs. new CPBCs (dotted). Four different
resolutions are shown: (Nr, L) = (21, 8), (31, 10), (41, 12), and (51, 14). The
outer boundary is at R = 41.9M .
This then is our version of the Sommerfeld boundary condition (cf. (25)), to be imposed
on a spherical boundary in the far field (where linearized theory is assumed to be valid).
Because the BSSN formulations using (27) are usually second-order in space,
there is no analogue of our three-index constraint (3) in that system. To mimic this
situation in our tests of equation (29), we also impose a CPBC on u2Aab as discussed in
section 2.2, which together with our constraint damping terms ensures that violations
of the three-index constraint (3) are exponentially damped.
Our version of Sommerfeld boundary conditions performs similarly on our test
problem (figure 5) to the freezing boundary conditions (26) (figure 4). The initial pulse
of reflections is smaller by ≈ 2 orders of magnitude, but later ||∆U|| grows to a similar
level as for freezing boundary conditions. Again the constraints do not converge away,
although this non-convergence appears only at somewhat higher resolutions than in
the freezing case.
3.3. Kreiss-Winicour boundary conditions
Recently, Kreiss and Winicour [14] proposed a set of ‘Sommerfeld-like’ CPBCs for the
harmonic Einstein equations and showed that they result in an IBVP that is well-posed
in the generalized sense. Their boundary conditions were implemented and tested in
[5]; here we compare their performance with the various other boundary treatments.
The Kreiss-Winicour boundary conditions are obtained by requiring the harmonic
constraint to vanish at the boundary,
= 0. (30)
In our notation, this can be written as an algebraic condition on part of the incoming
fields u1−,
= Fa, (31)
where
[2k(cδad) − kaψcd],
lbu1+ab −
bcu1+bc + P
iju2ija − 12Pa
iψbcu2ibc (32)
Testing outer boundary treatments for the Einstein equations 12
− γ2ta +Ha.
The range of the projection operator PC
is identical with that of PC defined in (10).
For the unconstrained incoming fields ũ1− (i.e. u1− without the γ2 term, equation
(8)), Kreiss and Winicour [14] specify certain free boundary data qPab and q
ab. In our
notation,
PP cdab ũ
cd = q
ab, P
ab ũ
cd = q
ab. (33)
In the linearized wave and gauge wave tests of [5], these boundary data are obtained
from the known exact solutions. In the absence of an exact solution, it is suggested
that the data could be obtained from an exterior Cauchy-characteristic or Cauchy-
perturbative code. However, since we do not have such an exterior code, we compute
the boundary data from the background solution, i.e. Schwarzschild spacetime. As in
the Sommerfeld case (section 3.2), we use a constraint-preserving boundary condition
on u2Aab to emulate the second-order formulation of [5, 14], and this value of u
Aab is
then used to compute Fa in (32).
Figure 6 shows the numerical results for our test problem. The magnitude of the
initial reflections lies between that of freezing and Sommerfeld boundary conditions
and is somewhat smaller at later times, though still larger than for our CPBCs at the
higher resolutions. The constraints converge away with increasing resolution, as they
should for a boundary condition that is consistent with the constraints. In a numerical
simulation, violations of the constraints are in general present in the interior of the
computational domain. These propagate as described by the constraint evolution
system (9) and some may hit the outer boundary. The Dirichlet boundary conditions
(30) might be expected to cause more reflections of constraint violations than our
no-incoming-field conditions (10), however, no indications of this are seen in figure
6. Probably the constraint damping we use is sufficiently effective in eliminating the
source of these reflections.
We shall see in section 5.1 that the Kreiss-Winicour boundary conditions also
cause larger errors in the physical degrees of freedom than our CPBCs. Since the
main difference between the two sets of boundary conditions is our use of a physical
boundary condition ∂tΨ0
= 0, we conclude that such a condition is crucial in reducing
the reflections from the outer boundary.
4. Alternate approaches
So far we have only considered boundary conditions that are local algebraic or
differential conditions imposed at the boundary of some finite computational domain.
There are of course many ways of treating the outer boundary that do not fall into that
category. In this section, we evaluate two such approaches: spatial compactification
and sponge layers.
4.1. Spatial compactification
Spatial compactification is a method that has been widely used in numerical relativity,
for instance in [41, 42] or more recently in the generalized harmonic binary black hole
simulations of Pretorius [15, 16, 17].
The basic idea is to introduce spatial coordinates that map spacelike infinity to a
finite location. Here we consider mappings that are functions of coordinate radius only
(whereas Pretorius applies the mapping to each Cartesian coordinate separately). We
Testing outer boundary treatments for the Einstein equations 13
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
Figure 6. Kreiss-Winicour (solid) vs. new CPBCs (dotted). Four different
resolutions are shown: (Nr, L) = (21, 8), (31, 10), (41, 12), and (51, 14). The
outer boundary is at R = 41.9M .
have used two such mappings, named Tan and Inverse, as detailed in Appendix B.1.
Each map has a scaleR across which the mapping is (essentially) linear. The outermost
grid point is placed at a very large but finite uncompactified radius (r = 1017M). With
respect to the compactified radial coordinate, the characteristic speeds are below
numerical roundoff there and hence no boundary condition should be needed. The
following results were produced using constraint-preserving boundary conditions; we
have checked for one simulation that using no boundary condition at all yields results
that are visually indistinguishable from the ones presented here on the scales of figures
7, 8, and 10.
As the waves travel outward, they become more and more blue-shifted with
respect to the computational grid and are eventually no longer properly resolved.
However, some form of artificial numerical dissipation is applied that acts as a low-
pass filter and causes the waves to be damped as they become increasingly distorted.
We have experimented with various such filters; see Appendix B.1 for details. One of
them (referred to as number 2 in the following) is designed to emulate as closely as
possible the fourth-order Kreiss-Oliger dissipation used by Pretorius.
In the following numerical comparisons, we evaluate the differences with respect
to the reference solution only in the part of the domain where the compatification map
is essentially linear, i.e. for r 6 R. First we compare the various filtering methods at
a fixed resolution, using the Tan compactification mapping (figure 7). The filters that
are applied to the right side of the evolution equations (numbers 1 and 3, cf. table B1)
do somewhat better than those applied to the solution itself (numbers 2 and 4), and
the Exponential filters (numbers 3 and 4) are slightly better than the Kreiss-Oliger
filters (numbers 1 and 2). All of them are outperformed by the CPBCs (imposed at
r = R). For our closest approximation to the dissipation used by Pretorius (number
2), ||∆U|| is comparable to constraint-preserving boundary conditions at the peak of
reflections (at t ≈ 2R) but becomes larger by about 2 orders of magnitude at later
times. The compactification methods also generate considerable constraint violations.
Next we focus on the best filter (number 4) of the previous test but vary the
resolution (figure 8). We do see convergence of ||∆U|| initially but the convergence
degrades at later times. This is surprising at first because with increasing resolution,
Testing outer boundary treatments for the Einstein equations 14
0 200 400 600 800 1000
t / M
New CPBC
TAN, Filter 1
TAN, Filter 2
TAN, Filter 3
TAN, Filter 4
0 200 400 600 800 1000
t / M
New CPBC
TAN, Filter 1
TAN, Filter 2
TAN, Filter 3
TAN, Filter 4
Figure 7. Tan compactification with various filters vs. new CPBCs. Only the
highest resolution (Nr, L) = (51, 14) is shown. The compactification scale (and
the radius of the outer boundary in the CPBC case) is R = 41.9M .
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
Figure 8. Tan compactification with filter 4 (solid) vs. new CPBCs (dotted).
Four different resolutions are shown: (Nr, L) = (21, 8), (31, 10), (41, 12), and
(51, 14). The compactification scale (and the radius of the outer boundary in the
CPBC case) is R = 41.9M .
the waves travel a longer distance before they fail to be resolved. Note however that
the high-frequency filter is applied at each time step, as is done in the simulations
of Pretorius. For higher resolutions, the time steps are smaller because of the CFL
condition and the filter is applied more often, thus leading to a stronger damping of
the waves. This may well lead to the observed loss of convergence with increasing
resolution. The constraints appear to converge away in this test, although from figure
8 it appears that this will not persist for even higher resolutions.
We have also evaluated the Inverse mapping described in Appendix B.1. The
results are similar, but somewhat worse than the Tan mapping results shown here.
Testing outer boundary treatments for the Einstein equations 15
4.2. Sponge layers
A method that has been used for a long time in computational science, in particular
for spectral methods (see e.g. section 17.2.3 of [43] and references therein), involves
so-called sponge layers. A sponge layer is introduced by modifying the evolution
equations according to
∂tu = . . .− γ(r)(u− u0), (34)
where u0 refers to the unperturbed background solution (Schwarzschild spacetime in
our case) and the smooth sponge function γ(r) > 0 is large only close to the outer
boundary of the computational domain. (Here we use uncompactified coordinates as
in sections 2 and 3.) In this way, the waves are damped exponentially as they approach
the outer boundary. Details on our particular choice of γ(r) can be found in Appendix
We compare the sponge layer method with our CPBCs in figure 9. For the
CPBCs, the boundary is either placed at R = 41.9M (the outer edge of the sponge-
free region) or at R = 121.9M (the outer edge of the sponge). At early times (t . 2R),
the ||∆U||∞ of the sponge layer method lies between that of the CPBCs for the two
choices of outer boundary radius, whereas at later times, it is much larger than both
versions of CPBCs. The constraint violations in the sponge runs do not converge
away.
5. Physical gravitational waves
Perhaps the most important predictions of numerical relativity simulations at the
present time are the gravitational waveforms produced by astrophysical systems like
binary black holes. It is important therefore to understand how the accuracy of these
waveforms is affected by the choice of boundary treatment. Physical gravitational
radiation can be described by the Newman-Penrose scalars Ψ4 and Ψ0. The scalar Ψ4
is dominated by the outgoing radiation (its ingoing part is suppressed by a factor of
(kr)4, where k is the wavenumber), whereas Ψ0 is dominated by the ingoing radiation
(its outgoing part is suppressed by a factor of (kr)4). In this section we compare
the gravitational waves extracted from the various boundary treatment solutions on
a sphere of radius r = Rex, using the methods described in Appendix A.5.
We note that Ψ4 (Ψ0) has a coordinate-invariant meaning only in the limit as
future (past) null infinity is approached. The quantities computed at finite radius r will
differ from those observed at infinity by terms of the order O(1/r). In the particular
case of perturbed Schwarzschild spacetime considered here, a gauge-invariant wave
extraction method does exist even at finite radius (see e.g. [44] and references therein)
but we do not adopt it here. Our purpose in this paper is merely to measure the
effects on Ψ4 caused by the various boundary treatments.
5.1. Difference of Ψ4 with respect to the reference solution
We begin by evaluating ∆Ψ4 ≡ Ψ4 − Ψref4 , where Ψ4 is the Newman-Penrose scalar
computed using one of the various boundary methods and Ψref4 is the same quantity
computed from the reference solution at the same extraction radius. The curves shown
in figure 10 plot the maximum value of |∆Ψ4| over time intervals of length 20M (this
time filtering averages over the high frequency quasi-normal oscillations of the black
hole), normalized by the maximum value of |∆Ψ4| over the entire evolution. The
Testing outer boundary treatments for the Einstein equations 16
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
CPBC at R = 41.9 M
0 200 400 600 800 1000
t / M
∞ (Nr,L)=(21,8)
(31,10)
(41,12)
(51,14)
CPBC at R = 41.9 M
0 200 400 600 800 1000
t / M
,L)=(21,8)
(31,10)
(41,12)
(51,14)
CPBC at R = 121.9 M
0 200 400 600 800 1000
t / M
∞ (Nr,L)=(21,8)
(31,10)
(41,12)
(51,14)
CPBC at R = 121.9 M
Figure 9. Sponge layer method (solid) vs. new CPBCs at two different radii
(dotted). Four different resolutions are shown: (Nr, L) = (21, 8), (31, 10), (41, 12),
and (51, 14). The size of the sponge-free region is R = 41.9M and ||∆U||∞ is only
computed for r 6 R.
radius of the outer boundary (or the compactification scale, or the size of the sponge-
free region, respectively) used for these comparisons is R = 41.9M , and the radiation
is extracted nearby at Rex = 40M .
The first peak in |∆Ψ4| seen in figure 10 arises as the wave in our test problem
passes outward through the extraction sphere at t ≈ Rex. This peak is caused by
a presently unknown (probably gauge) interaction between the outer boundary (or
compactified region etc.) and the spacetime near the extraction sphere. We have
verified that this interaction and its influence on the peak in ∆Ψ4 goes away if we
move the outer boundary (or the extraction surface) so that they are not in causal
contact as the outgoing wave pulse passes the extraction surface.
Some of the outgoing radiation is reflected off the boundary. Most of this reflected
radiation is subsequently absorbed by the black hole, but some of it excites the hole,
which then emits quasi-normal mode radiation of exponentially decaying amplitude.
This exponential decay can be clearly seen for most of the boundary treatments.
In the case of freezing boundary conditions, nearly all of the outgoing quasi-
normal mode radiation is reflected from the boundary because the reflection coefficient
is nearly 1 for the wave number of the dominant mode, k = 0.376/M (cf. figure 1).
Testing outer boundary treatments for the Einstein equations 17
It then re-excites the black hole, which again radiates and so forth. On average the
amplitude of the reflections remains roughly constant in time for this case. This
behaviour is consistent with the result shown in figure 3 of [6] for a similar perturbed
black hole simulation.
For the Sommerfeld and Kreiss-Winicour boundary conditions, the reflections are
much smaller but still considerably larger (by 2 to 3 orders of magnitude) than for
our CPBCs. We attribute this difference largely to our use of the physical boundary
condition (12).
The spatial compactification method has the largest difference |∆Ψ4|, particularly
at early times t ∼ R (about 4 orders of magnitude larger than for the CPBCs). We
suspect that this may be a consequence of the use of artificial dissipation, as discussed
in section 4.1.
The sponge layer method has the smallest errors at early times. This is not
surprising because the outer boundary of the sponge layer is much further out at
R = 121.9M . However at later times when the waves begin to interact with the
sponge layer, this method causes reflections comparable in amplitude to those using
Sommerfeld boundary conditions.
We also note that at late times the level of |∆Ψ4| decreases significantly with
resolution for the CPBCs, but not generally for the other boundary treatments.
We think it is remarkable that the maximum relative error in the extracted
physical radiation is quite small (10−5 to 10−3) in these tests, even for the less
sophisticated boundary treatments such as the freezing or Sommerfeld boundary
conditions. This success is due in part to the fact that the extraction radius,
Rex = 40M , for this test problem is about ten wavelengths (of the initial radiation
pulse) away from the central black hole. Our results are likely to be more accurate
than those from typical binary black hole simulations, which place the outer boundary
at two or three wavelengths. This suggests that current binary black hole codes
using, for instance, Sommerfeld boundary conditions, can still produce waveforms
that are useful for some aspects of gravitational wave data analysis provided the
outer boundary is placed sufficiently far out. Data analysis applications needing
high precision waveforms, however, such as source parameter measurement or high-
amplitude supermassive binary black hole signal subtraction for LISA, will need to
use a more sophisticated boundary treatment that produces smaller errors in Ψ4.
5.2. Comparison with the predicted reflection coefficient
Buchman and Sarbach [18, 19] have recently developed a hierarchy of increasingly
absorbing physical boundary conditions for the Einstein equations by analyzing
the equations describing the evolution of the Weyl curvature on both a flat and
a Schwarzschild background spacetime. Their analysis predicts, in particular, the
reflection coefficient ρ (defined as the ratio of the ingoing to the outgoing parts of the
solution) that arises from the ∂tΨ0
= 0 physical boundary condition that we use.
For quadrupolar radiation (as in our numerical tests), this reflection coefficient is
given by equation (89) of [18],
ρ(kR) = 3
(kR)−4 +O(kR)−5, (35)
where k is the wave number of the gravitational radiation and R is the boundary
radius. (As explained at the beginning of section 2.3, we assume the background
spacetime to be flat; effects due to the backscattering would only enter at O(M/R).)
Testing outer boundary treatments for the Einstein equations 18
0 200 400 600 800 1000
t / M
,L)=(31,10)
(51,14)
Freezing
0 200 400 600 800 1000
t / M
,L)=(31,10)
(51,14)
Sommerfeld
0 200 400 600 800 1000
t / M
,L)=(31,10)
(51,14)
Kreiss-Winicour
0 200 400 600 800 1000
t / M
,L)=(31,10)
(51,14)
TAN compactification (filter 4)
0 200 400 600 800 1000
t / M
,L)=(31,10)
(51,14)
Sponge layer
Figure 10. Difference of Ψ4 for the various alternate methods (solid) vs. the
new CPBCs (dotted). Two resolutions are shown: (Nr, L) = (31, 10) and (51, 14).
The radius of the outer boundary (or the compactification scale, or the size of
the sponge-free region, respectively) is R = 41.9M and the waves are extracted
at Rex = 40M .
Testing outer boundary treatments for the Einstein equations 19
0 1 2 3 4 5 6 7
(k R)
R = 21.9 M
0 1 2 3 4 5 6 7
(k R)
R = 121.9 M
Figure 11. Comparison of the time Fourier transform of the measured Ψ0(t)
with 3
(kR)−4Ψ4, which is the predicted value using the reflection coefficient of
[18].
By evaluating Ψ0 and Ψ4 at the extraction radius of our test, we find that the ratio
Ψ0/Ψ4 agrees with their predicted ρ to leading order in 1/(kR). We note that the
tetrad we use for wave extraction (Appendix A.5) does not agree exactly with that
of [18]. However, the tetrads do agree for the unperturbed Schwarzschild solution,
so that the errors introduced into Ψ0 and Ψ4 due to our different choice of tetrad
are second-order small in perturbation theory and hence the comparison with [18] is
consistent.
For a numerical solution using our new CPBCs, we evaluate the Newman-Penrose
scalars Ψ0(t) and Ψ4(t) on extraction spheres located 1.9M inside the outer boundary.
In figure 11 we plot the time Fourier transforms of these quantities. We also plot
(kR)−4Ψ4, which by the above argument should agree with Ψ0 to leading order in
1/(kR). Figure 11 shows that the numerical agreement is reasonably good: roughly
at the expected level of accuracy. The overall dependence of the predicted reflection
coefficient ρ on k and R is captured very well. We surmise that the levelling off of
our numerical Ψ0 for k & 3 is due to numerical roundoff effects. (Note the magnitude
of Ψ0 at those frequencies.) For radii R & 200M , Ψ0 is at the roundoff level for all
frequencies.
6. Discussion
The purpose of this paper is to compare various methods of treating the outer
boundary of the computational domain. We evaluate the performance of several often-
used boundary treatments in numerical relativity by measuring the amount of spurious
reflections and constraint violations they generate. To this end, we consider as a test
problem an outgoing gravitational wave superimposed on a Schwarzschild black hole
spacetime. First we compute this numerical solution on a reference domain, large
enough that the influence of the outer boundary can be neglected. Then we repeat the
evolution on smaller domains using one of the boundary treatments, either imposing
local boundary conditions, compactifying the domain using a radial coordinate map,
or installing a sponge layer. We use a first-order generalized harmonic formulation
of the Einstein equations, although these boundary methods can be applied to other
Testing outer boundary treatments for the Einstein equations 20
formulations as well. We believe our results are fairly independent of the particular
formulation used.
Our main conclusion is that our version of constraint-preserving boundary
conditions performs better than any of the alternate treatments that we tested. Our
boundary conditions include a limitation on the influx of spurious gravitational waves
by freezing the Newman-Penrose scalar Ψ0 at the boundary. We also introduce and
test an improved boundary condition for the gauge degrees of freedom.
For some of the simple boundary conditions, such as freezing or Sommerfeld
conditions, we find constraint violations that do not converge away with increasing
resolution. The continuum limit does not satisfy Einstein’s equations in these cases.
Most of the alternate boundary conditions also generate considerable reflections as
measured by ∆U , the norm of the difference with respect to the reference solution. In
many cases, these reflections do not decrease significantly with increasing resolution.
The difference norm ∆U that we use to measure boundary reflections includes the
entire spacetime metric, not just the physical degrees of freedom. It is important then
to evaluate separately the effects of the various boundary treatments on the physical
degrees of freedom. We use the extracted outgoing radiation as approximated by
the Newman-Penrose scalar Ψ4 for this purpose. Here our conclusions are somewhat
different. Rather surprisingly, most of the boundary methods we consider generate
relatively small errors in Ψ4. This suggests that if gravitational waveforms are only
needed to an accuracy of, say, 1% (which is comparable to the discrepancies between
recent binary black hole simulations [45]) then even the simple Sommerfeld conditions
might be good enough. (For those, we find relative errors ∼ 10−5.) The largest
relative errors in Ψ4 we find (∼ 10−2) occur with our implementation of the spatial
compactification method used by Pretorius [15, 16, 17]. We attribute these largely to
the use of artificial dissipation. Undesirable effects of dissipation might be somewhat
less severe in binary black hole evolutions, which have much larger wavelengths
(λ ∼ 20 − 100M) than ours (λ ∼ 4M). Our tests suggest that the errors in Ψ4 can
be made to decrease significantly with resolution only by using more sophisticated
constraint preserving and physical boundary conditions. The importance of using a
physical boundary condition on Ψ0 is illustrated in particular by the difference between
the performance of our boundary conditions and those of Kreiss and Winicour [14].
Some caveats regarding the interpretation of our results must be stated. First,
the ratio of the dominant wavelength to the radius of the outer boundary is typically
much larger for binary black hole evolutions (where λ/R & 0.5) than for the simple
test problem considered here (where λ/R ∼ 0.1). Boundary treatments generally work
better for smaller λ/R, i.e. when the boundary is well out in the wave zone. Hence the
results presented here are likely to be more accurate than those from typical binary
black hole simulations. Second, we use spectral methods rather than finite-difference
methods, which are more commonly used in numerical relativity at this time. This
complicates the implementation of the kind of numerical dissipation that is crucial for
the spatial compactification method to work. While we have attempted to construct
a filter that mimics the finite-difference dissipation as closely as possible, a direct
comparison is clearly impossible. In finite-difference methods, the error introduced by
the type of numerical dissipation considered here is below the truncation error. Hence
tests similar to ours but performed with a finite-difference method would not be able
to detect the effect of dissipation.
There are several directions in which the present work could be extended. For
large values of the outer boundary radius, we observe a non-convergent power-
Testing outer boundary treatments for the Einstein equations 21
law growth of the error in our test problem when constraint-preserving boundary
conditions are used; the origin of this growth should be investigated further. It would
be interesting to implement and test the hierarchy of physical boundary conditions that
are perfectly absorbing for linearized gravity (including leading-order corrections due
to the curvature and backscatter) found recently by Buchman and Sarbach [18, 19].
Our boundary conditions could also be tested using known exact solutions such as
gauge waves, and comparisons could be made with the results found in [5].
For completeness we also mention a number of additional outer boundary
approaches that were not addressed in this paper, but would also be interesting future
extensions of this research. In [46, 47], boundary conditions for the full nonlinear
Einstein equations on a finite domain are obtained by matching to exact solutions of
the linearized field equations at the boundary. Alternatively, the interior code could
be matched to an ‘outer module’ that solves the linearized field equations numerically
[48, 49, 50, 51]. Other approaches involve matching the interior nonlinear Cauchy code
to an outer characteristic code (see [52] for a review) or using hyperboloidal spacetime
slices that can be compactified towards null infinity (see [53] for a review).
Appendix A. Details on the numerical test problem
Appendix A.1. Initial data
The initial data used for our numerical tests are the same as in [27]. The background
solution is a Schwarzschild black hole in Kerr-Schild coordinates,
ds2 = −dt2 + 2M
(dt+ dr)2 + dr2 + r2dΩ2. (A.1)
Throughout the paper, M refers to the bare black hole mass of the unperturbed
background. We superpose an odd-parity outgoing quadrupolar wave perturbation
constructed using Teukolsky’s method [54]. Its generating function is taken to be a
Gaussian G(r) = A exp[−(r − r0)2/w2] with A = 4× 10−3, r0 = 5M , and w = 1.5M .
The full non-linear initial value equations in the conformal thin sandwich formulation
are then solved to obtain initial data that satisfy the constraints [55]. This yields
initial values for the spatial metric, extrinsic curvature, lapse function, and shift
vector. We note that after the superposition, the resulting solution is still nearly
but not completely outgoing.
Our generalized harmonic formulation of Einstein’s equations requires initial data
for the full spacetime metric and its first time derivative. These can be computed from
the 3+1 quantities obtained above, provided we also choose initial values for the time
derivatives of the lapse function and shift vector. These initial time derivatives are
freely specifiable and are equivalent to the initial choice of the gauge source function
Ha; we choose ∂tN = 0 and ∂tN i = 0 at t = 0.
Appendix A.2. Numerical method
We use a pseudospectral collocation method as described for example in [27].
The computational domain for the test problem considered here is taken to
be a spherical shell extending from r = 1.9M (just inside the horizon) out to
some r = R. This domain is subdivided into spherical-shell subdomains of extent
∆r = 10M . On each subdomain, the numerical solution is expanded in Chebyshev
polynomials in the radial direction and in spherical harmonics in the angular directions
Testing outer boundary treatments for the Einstein equations 22
(where each Cartesian tensor component is expanded in the standard scalar spherical
harmonics). Typical resolutions are Nr ∈ {21, 31, 41, 51} coefficients per subdomain
for the Chebyshev series and l 6 L with L ∈ {8, 10, 12, 14} for the spherical harmonics.
We change the outer boundary radius R by changing the number of subdomains
while keeping the width ∆r of each subdomain fixed; this facilitates direct comparisons
between runs with different values of R. For example, the innermost four subdomains
of the reference solution (which has a total of 96 subdomains and R = 961.9M) are
identical to the four subdomains used to compute the solution with R = 41.9M .
The evolution equations are integrated in time using a fourth-order Runge-Kutta
scheme, with a Courant factor ∆t/∆xmin of at most 2.25, where ∆xmin is the smallest
distance between two neighbouring collocation points. As described in [27], the top
four coefficients in the tensor spherical harmonic expansion of each of our evolved
quantities is set to zero after each time step; this eliminates an instability associated
with the inconsistent mixing of tensor spherical harmonics in our approach.
We use two methods of numerically implementing boundary conditions; the choice
of method depends on the type of boundary conditions. Boundary conditions that can
be expressed as algebraic relations involving the characteristic fields are implemented
using a penalty method (see [56] and references therein; in the context of finite-
difference methods see also [57] and references therein). In particular, we use a penalty
method to implement the Kreiss-Winicour boundary conditions (cf. section 3.3)
and to impose boundary conditions at the internal boundaries between neighbouring
subdomains. Boundary conditions that are expressed in terms of the time derivatives
of the characteristic fields are implemented using the method of Bjørhus [58], where
the time derivatives of the incoming characteristic fields are replaced at the boundary
with the relevant boundary condition. All boundary conditions in this paper besides
those mentioned above are implemented using the Bjørhus method.
Appendix A.3. Gauge source functions
Our generalized harmonic formulation [6] of Einstein’s equations allows for gauge
source functions that depend arbitrarily on the coordinates and the spacetime metric:
Ha = Ha(t, x, ψ). The generalized harmonic evolution equations are equivalent to
Einstein’s equations only if the constraint (4) remains satisfied.
We choose the time derivatives of lapse and shift to be zero at the beginning of
the simulation; this determines the initial value of Ha via the constraint (4). For the
subsequent evolution, we hold this Ha fixed in time.
Appendix A.4. Error quantities
We use two different measurements of the errors in our solutions, which we monitor
during our numerical evolutions. First, given a numerical solution (ψab,Πab,Φiab),
the difference between that solution and the reference solution (ψ(ref)ab ,Π
(ref)
ab ,Φ
(ref)
iab ) is
computed with the following norm at each point in space,
δabδcd(M−2∆ψac∆ψbd + ∆Πac∆Πbd
+gij∆Φiac∆Φjbd)
, (A.2)
where ∆ψab means ψab−ψ(ref)ab , and similarly for ∆Πab and ∆Φiab. Second, we define
a quantity C that measures the violations in all of the constraints of our system,
δab(FaFb + gij(CiaCjb + gklδcdCikacCjlbd)
Testing outer boundary treatments for the Einstein equations 23
+M−2(CaCb + gijδcdCiacCjbd))
, (A.3)
where Fa and Cia are first derivatives of Ca defined in [6]. To compute global error
measures, a spatial norm ||·||, either the L∞ norm or the L2 norm, is applied separately
to ∆U and C.
The question often arises as to the significance of particular values of ||∆U|| and
||C||. For example, is a simulation with ||C|| = 10−2 good to one percent accuracy? To
make it easier to answer such questions, we normalize both ||∆U|| and ||C|| as follows,
and we always plot normalized quantities.
We divide ||∆U|| by a normalization factor ||∆U0||, defined as the difference
between a given solution at t = 0 and the unperturbed Schwarzschild background;
i.e., the quantity ||∆U0|| is computed from (A.2) using the unperturbed Schwarzschild
solution instead of the reference solution. Since ||∆U0|| is evaluated at t = 0, it depends
only on the initial data used in the simulation, and is a measure of the amplitude
of the superposed gravitational wave perturbation. For the initial data used here,
||∆U0||∞ = 6× 10−3 and ||∆U0||2 = 1.4× 10−4. The quantity ||∆U||/||∆U0|| is more
easily interpreted than ||∆U||; for example, ||∆U||/||∆U0|| is unity when the difference
from the reference solution is of the same size as the initial perturbation.
Similarly, the constraint energy norm ||C|| is divided by the norm of the first
derivatives ||∂U|| (at the respective time),
gijδabδcd(M−2∂iψac∂lψbd + ∂iΠac∂jΠbd
+gkl∂iΦkac∂jΦlbd)
. (A.4)
The constraints for our system are linear combinations of the first derivatives of the
fields, hence ||C||/||∂U|| ∼ 1 corresponds to a complete violation of the constraints.
Appendix A.5. Wave extraction
For evaluating gravitational waveforms, we compute the Newman-Penrose scalars
Ψ0 = −Cabcdlamblcmd, Ψ4 = −Cabcdkam̄bkcm̄d, (A.5)
where Cabcd is the Weyl tensor, la and ka are outgoing and ingoing null vectors
normalized according to laka = −1, ma is a complex unit null spatial vector orthogonal
to la and ka, and m̄a is the complex conjugate of ma. For perturbations of flat
spacetime, there is a standard choice for the vectors la, ka, and ma. In general curved
spacetimes, however, no such prescription for the tetrad exists that would produce
coordinate-independent quantities Ψ0 and Ψ4 at finite radius. We choose the null
vectors according to
la = 1√
(ta + na) , ka = 1√
(ta − na) , (A.6)
where ta is the future-pointing unit timelike normal to the t = const. slices and na is
the unit spacelike normal to the extraction sphere. Finally, we choose
sin θ
, (A.7)
where (r, θ, φ) are spherical coordinates on the r = Rex = const. extraction sphere.
Note that our choice of ma is not exactly null nor of unit magnitude at finite extraction
radius. However, the tetrad is orthonormal for the unperturbed Schwarzschild
solution, so that the errors introduced into Ψ0 and Ψ4 because of the lack of tetrad
orthonormality will be second-order small in perturbation theory.
Testing outer boundary treatments for the Einstein equations 24
The quantity Ψ4 corresponds to outgoing radiation in the limit of r → ∞,
t − r = const., i.e. as future null infinity is approached. Similarly Ψ0 corresponds
to ingoing radiation as past null infinity is approached. At finite extraction radius,
Ψ4 and Ψ0 will disagree with the waveforms observed at infinity by terms of the order
O(Rex)−1.
We decompose the quantities Ψ4 and Ψ0 in terms of spin-weighted spherical
harmonics of spin-weight −2 on the extraction surface. Since our perturbation is
an odd-parity quadrupole wave, the imaginary part of the (l = 2, m = 0) spherical
harmonic is by far the dominant contribution to Ψ4, and we only display that mode
in our plots. We normalize the curves in our graphs by the maximum (in time) value
of |Ψ4| at the extraction radius Rex, which for Rex = 40M is max |Ψ4| = 6× 10−4.
Appendix B. Details of the alternate approaches
In this appendix, we provide some more details on the alternate boundary treatments
discussed in section 4: spatial compactification and sponge layers.
Appendix B.1. Spatial compactification
We implement spatial compactification by introducing a radial coordinate
transformation x → r(x) that maps a compact ball on the computational grid with
x ∈ [0, xmax] to the full unbounded physical slice with r ∈ [0,∞]. We consider two
such mappings. The Tan mapping is similar to the one used by Pretorius [15, 16, 17]
and is given by
rTan(x) = R tan
, 0 6 x < 2R. (B.1)
The scale R determines the range in physical radius r across which the map is
essentially linear (see figure B1). When comparing compactification with other
boundary treatments, we compare quantities only in the region r < R. (The scale
R is equal to unity in the work of Pretorius. He uses mesh refinement to obtain the
appropriate resolution close to the origin, while we fix the resolution and choose the
scale R appropriately.) We also tested an Inverse map defined by
rInverse(x) =
x, 0 6 x 6 R ,
2R− x, R < x < 2R,
(B.2)
see figure B1. This map is only C1 at x = R, but we maintain spectral accuracy in
our tests by placing this surface at the boundary between spectral subdomains.
Dissipation is needed to remove the short wavelength components of the
waves as they travel outward on the compactified computational grid and become
unresolved. We apply this dissipation only in the radial direction, but everywhere
in the computational domain. In spectral methods, dissipation can be conveniently
implemented in the form of a spectral filter. This filter is applied by multiplying each
spectral expansion coefficient of index k by a function f(k). (See Appendix A.2 for
details on the pseudospectral method we use.) Higher values of k correspond to shorter
wavelengths in the numerical approximation; let kmax be the highest index used in the
spectral expansion. The first filter function we consider is the closest analogue in the
context of our spectral methods to Kreiss-Oliger [59] dissipation,
fKreiss-Oliger(k) = 1− � sin4
2kmax
, 0 6 � 6 1. (B.3)
Testing outer boundary treatments for the Einstein equations 25
0 1 2
x / R
INVERSE
0 0.5 1
k / k
Filters 1 and 2
Filters 3 and 4
Figure B1. Compactification mappings (left) and filter functions (right). The
dashed line indicates the boundary of the region in where the compactification
mapping is (essentially) linear.
Typical values of the parameter � used by Pretorius are � ∈ [0.2, 0.5]; we use � = 0.25.
This filter was derived via a comparison with finite-difference methods as follows.
In the finite-difference approach, a numerical solution u is represented on a set of
equidistant grid points xj . (It suffices to consider the one-dimensional case here.)
Some form of numerical dissipation is usually required for the finite-difference method
to be stable. The one that is most often used for second-order accurate methods is
fourth-order Kreiss-Oliger dissipation [59]. One possible implementation of this, used
e.g. by Pretorius, amounts to replacing
u→ F [u] ≡
u (B.4)
at each time step, where h is the grid spacing and D4 is the second-order accurate
centred finite difference operator approximating the fourth derivative,
D4ui = h
−4(uj−2 − 4uj−1 + 6uj − 4uj+1 + uj+2). (B.5)
Taking u to be a Fourier mode u(k)j = exp(ikxj), it follows that the mode is damped
by a frequency-dependent factor,
u(k) → F [u(k)] ≡
1− � sin4
2kmax
u(k), (B.6)
where kmax = π/(2h) is the Nyquist frequency. Thus we obtain the filter function
(B.3). Strictly speaking, the above analysis only applies to Fourier expansions and
not to the Chebyshev expansions we use. Nevertheless, we apply the filter in the
form (B.3) to our Chebyshev expansion coefficients. Note that in (B.6), each spectral
coefficient u(k) is filtered separately; this is not true for the analogous calculation for
a Chebyshev expansion.
We also use a different filter function, which we call the Exponential filter, that
is often used in spectral methods (see [60] and references therein),
fExponential(k) = exp
σkmax
. (B.7)
Testing outer boundary treatments for the Einstein equations 26
No. Type Parameters Applied to
1 Kreiss-Oliger � = 0.25 right side
2 Kreiss-Oliger � = 0.25 solution
3 Exponential σ = 0.76, p = 13 right side
4 Exponential σ = 0.76, p = 13 solution
Table B1. Details of the filtering methods
Typical values of the parameters are σ = 0.76 and p = 13. This choice of parameters
gives less dissipation at small values of k than the Kreiss-Oliger filter, and also ensures
that f(kmax) ≈ 10−16 is at the level of the numerical roundoff error.
There are various ways the filters can be applied in a numerical evolution. We
have experimented with two different methods. In the first method, the filter is applied
to the right side of the equations, i.e. the evolution equations ∂tu = S are modified
according to ∂tu = F [S], where F [S] is the filtered right side. In the second method,
the filter is instead applied to the solution itself, i.e. after each substep of the time
integrator (cf. Appendix A.2), the numerical solution u is replaced with its filtered
version F [u]. This second method is closest to how the Kreiss-Oliger filter is applied
by Pretorius.
For our numerical tests, we have used four different combinations of the various
options described above. They are summarized in table B1.
Appendix B.2. Sponge layers
For sponge layers we must specify a sponge profile function γ(r), as defined in (34). We
choose γ(r) to be nonzero only outside some sponge-free region of radius R, and when
comparing sponge layers with other boundary treatments, we compare quantities only
in the sponge-free region r < R.
The sponge profile function γ(r) we use is a Gaussian centred at the outer
boundary, which we choose to place at r = 3R,
γ(r) = γ0 exp
r − 3R
. (B.8)
The amplitude of the Gaussian is taken to be γ0 = 1. The width σ is chosen so that
γ(r) 6 10−16 (the numerical roundoff error) for r 6 R, which requires σ . R/3. In
our numerical example, we take R = 41.9M and σ = 13.3M . Hence σ is considerably
larger than the wavelength λ ≈ 4M of the gravitational wave, which is required in
order to avoid reflections from the sponge layer (cf. section 17.2.3 of [43]). Figure B2
shows a plot of this sponge profile.
Acknowledgments
We thank Luisa Buchman, Jan Hesthaven, Larry Kidder, Harald Pfeiffer, Olivier
Sarbach, and Jeff Winicour for helpful discussions concerning this work. The
numerical simulations presented here were performed using the Spectral Einstein
Code (SpEC) developed at Caltech and Cornell primarily by Larry Kidder, Mark
Scheel and Harald Pfeiffer. This work was supported in part by grants from the
Sherman Fairchild Foundation, and from the Brinson Foundation; by NSF grants
Testing outer boundary treatments for the Einstein equations 27
0 1 2 3
r / R
Figure B2. The sponge profile function γ(r). The dashed line indicates the
boundary of the region where γ is below the numerical roundoff error.
PHY-0099568, PHY-0244906, PHY-0601459, DMS-0553302 and NASA grants NAG5-
12834, NNG05GG52G.
References
[1] Novak J and Bonazzola S 2004 Absorbing boundary conditions for simulation of gravitational
waves with spectral methods in spherical coordinates J. Comput. Phys. 197 86–196
[2] Rinne O 2005 Axisymmetric Numerical Relativity Ph.D. thesis Univ. of Cambridge Preprint
http://www.arxiv.org/abs/gr-qc/0601064
[3] Lau S L 2004 Rapid evaluation of radiation boundarz kernels for time-domain wave propagation
on black holes: implementation and numerical tests Class. Quantum Grav. 21 4147–4192
[4] Babiuc M C, Szilágyi B and Winicour J 2006 Harmonic initial-boundary evolution in general
relativity Phys. Rev. D 73 064017
[5] Babiuc M C, Kreiss H O and Winicour J 2007 Constraint-preserving Sommerfeld conditions for
the harmonic Einstein equations Phys. Rev. D 75 044002
[6] Lindblom L, Scheel M A, Kidder L E, Owen R and Rinne O 2006 A new generalized harmonic
evolution system Class. Quantum Grav. 23 S447–S462
[7] Brügmann B, Tichy W and Jansen N 2004 Numerical simulation of orbiting black holes Phys.
Rev. Lett. 92 211101
[8] Campanelli M, Lousto C O, Marronetti P and Zlochower Y 2006 Accurate evolutions of orbiting
black-hole binaries without excision Phys. Rev. Lett. 96 111101
[9] Baker J G, Centrella J, Choi D I, Koppitz M and van Meter J 2006 Gravitational-wave extraction
from an inspiraling configuration of merging black holes Phys. Rev. Lett. 96 111102
[10] Diener P, Herrmann F, Pollney D, Schnetter E, Seidel E, Takahashi R, Thornburg J and Ventrella
J 2006 Accurate evolution of orbiting binary black holes Phys. Rev. Lett. 96 121101
[11] Herrmann F, Hinder I, Shoemaker D and Laguna P 2007 Unequal-mass binary black hole plunges
and gravitational recoil Class. Quantum Grav. 24 S33–S42
[12] Shibata M and Nakamura T 1995 Evolution of three-dimensional gravitational waves: Harmonic
slicing case Phys. Rev. D 52 5428
[13] Baumgarte T W and Shapiro S L 1998 Numerical integration of Einstein’s field equations Phys.
Rev. D 59 024007
[14] Kreiss H O and Winicour J 2006 Problems which are well-posed in a generalized sense with
applications to the Einstein equations Class. Quantum Grav. 16 S405–S420
[15] Pretorius F 2005 Numerical relativity using a generalized harmonic decomposition Class.
Quantum Grav. 22 425–452
[16] Pretorius F 2005 Evolution of binary black hole spacetimes Phys. Rev. Lett. 95 121101
[17] Pretorius F 2006 Simulation of binary black hole spacetimes with a harmonic evolution scheme
Class. Quantum Grav. 23 S529–S552
[18] Buchman L T and Sarbach O C A 2006 Towards absorbing outer boundaries in general relativity
Class. Quantum Grav. 23 6709–6744
Testing outer boundary treatments for the Einstein equations 28
[19] Buchman L T and Sarbach O C A 2007 Improved outer boundary conditions for Einstein’s field
equations Class. Quantum Grav. 24 S307–S326
[20] Gundlach C, Calabrese G, Hinder I and Mart́ın-Garćıa J M 2005 Constraint damping in the Z4
formulation and harmonic gauge Class. Quantum Grav. 22 3767–3774
[21] Stewart J M 1998 The Cauchy problem and the initial boundary value problem in numerical
relativity Class. Quantum Grav. 15 2865–2889
[22] Friedrich H and Nagy G 1999 The initial boundary value problem for Einstein’s vacuum field
equations Comm. Math. Phys. 201 619–655
[23] Iriondo M S and Reula O A 2002 Free evolution of self-gravitating, spherically symmetric waves
Phys. Rev. D 65 044024
[24] Calabrese G, Lehner L and Tiglio M 2002 Constraint-preserving boundary conditions in
numerical relativity Phys. Rev. D 65 104031
[25] Calabrese G and Sarbach O 2003 Detecting ill-posed boundary conditions in general relativity
J. Math. Phys. 44 3888–3899
[26] Calabrese G, Pullin J, Reula O, Sarbach O and Tiglio M 2003 Well posed constraint-preserving
boundary conditions for the linearized Einstein equations Comm. Math. Phys. 240 377–395
[27] Kidder L E, Lindblom L, Scheel M A, Buchman L T and Pfeiffer H P 2005 Boundary conditions
for the Einstein evolution system Phys. Rev. D 71 064020
[28] Bona C, Ledvinka T, Palenzuela-Luque C and Žáček M 2005 Constraint-preserving boundary
conditions in the Z4 numerical relativity formalism Class. Quantum Grav. 22 2615–2634
[29] Sarbach O and Tiglio M 2005 Boundary conditions for Einstein’s field equations: Analytical and
numerical analysis J. Hyp. Diff. Eq. 2 839–883
[30] Bardeen J M and Buchman L T 2002 Numerical tests of evolution systems, gauge conditions,
and boundary conditions for 1D colliding gravitational plane waves Phys. Rev. D 65 064037
[31] Nagy G and Sarbach O 2006 A minimization problem for the lapse and the initial-boundary
value problem for Einstein’s field equations Class. Quantum Grav. 16 S477–S504
[32] Rinne O 2006 Stable radiation-controlling boundary conditions for the generalized harmonic
Einstein equations Class. Quantum Grav. 23 6275–6300
[33] Scheel M A, Pfeiffer H P, Lindblom L, Kidder L E, Rinne O and Teukolsky S A 2006 Solving
Einstein’s equations with dual coordinate frames Phys. Rev. D 74 104006
[34] Pfeiffer H P, Brown D A, Kidder L E, Lindblom L, Lovelace G and Scheel M A 2007 Reducing
orbital eccentricity in binary black hole simulations Class. Quantum Grav. 24 S59–S81
[35] Chrzanowski P L 1975 Vector potential and metric perturbations of a rotating black hole Phys.
Rev. D 11 2042–2062
[36] Bayliss A and Turkel E 1980 Radiation boundary conditions for wave-like equations Comm. Pure
Appl. Math. 33 707–725
[37] Rauch J 1985 Symmetric positive systems with boundary characteristics of constant multiplicity
Trans. Am. Math. Soc. 291 167–187
[38] Secchi P 1996 The initial boundary value problem for linear symmetric hyperbolic systems with
characteristic boundary of constant multiplicity Diff. Int. Eq. 9 671–700
[39] Secchi P 1996 Well-posedness of characteristic symmetric hyperbolic systems
Arch. Rat. Mech. Anal. 134 155–197
[40] Szilágyi B, Pollney D, Rezzolla L, Thornburg J and Winicour J 2007 An explicit harmonic code
for black-hole evolution using excision Class. Quantum Grav. 24 S275–S293
[41] Garfinkle D and Duncan G 2001 Numerical evolution of Brill waves Phys. Rev. D 63 044011
[42] Choptuik M, Lehner L, Olabarrieta I, Petryk R, Pretorius F and Villegas H 2003 Towards the
final fate of an unstable black string Phys. Rev. D 68 044001
[43] Boyd J P 2001 Chebyshev and Fourier Spectral Methods 2nd ed (Dover publications)
[44] Pazos E, Dorband E N, Nagar A amd Palenzuela C, Schnetter E and Tiglio M 2007 How far
away is far enough for extracting numerical waveforms, and how much do they depend on the
extraction method? Class. Quantum Grav. 24 S341–S368
[45] Baker J G, Campanelli M, Pretorius F and Zlochower Y 2007 Comparisons of binary black hole
merger waveforms Class. Quantum Grav. 24 S25–S31
[46] Abrahams A M and Evans C R 1988 Reading off the gravitational radiation waveforms in
numerical relativity calculations: Matching to linearized gravity Phys. Rev. D 37 318
[47] Abrahams A M and Evans C R 1990 Gauge-invariant treatment of gravitational radiation near
the source: Analysis and numerical simulations Phys. Rev. D 42 2585
[48] Abrahams A M et al. 1998 Gravitational wave extraction and outer boundary conditions by
perturbative matching Phys. Rev. Lett. 80 1812–1815
[49] Rupright M E, Abrahams A M and Rezzolla L 1998 Cauchy-perturbative matching and outer
boundary conditions I: Methods and tests Phys. Rev. D 58 044005
Testing outer boundary treatments for the Einstein equations 29
[50] Rezzolla L, Abrahams A M, Matzner R A, Rupright M E and Shapiro S L 1999 Cauchy-
perturbative matching and outer boundary conditions: Computational studies Phys. Rev.
D 59 064001
[51] Zink B, Pazos E, Diener P and Tiglio M 2006 Cauchy-perturbative matching reexamined: Tests
in spherical symmetry Phys. Rev. D 73 084011
[52] Winicour J 2005 Characteristic evolution and matching Living Rev. Relativity 8(10)
[53] Frauendiener J 2004 Conformal infinity Living Rev. Relativity 7(1)
[54] Teukolsky S A 1982 Linearized quadrupole waves in general relativity and the motion of test
particles Phys. Rev. D 26 745–750
[55] Pfeiffer H P, Kidder L E, Scheel M A and Shoemaker D 2005 Initial data for Einstein’s equations
with superposed gravitational waves Phys. Rev. D 71 024020
[56] Hesthaven J S 2000 Spectral penalty methods Appl. Numer. Math. 33 23–41
[57] Schnetter E, Diener P, Dorband E N and Tiglio M 2006 A multi-block infrastructure for three-
dimensional time-dependent numerical relativity Class. Quantum Grav. 23 S553–S578
[58] Bjørhus M 1995 The ODE formulation of hyperbolic PDEs discretized by the spectral collocation
method SIAM J. Sci. Comput. 16 542–557
[59] Kreiss H O and Oliger J 1973 Methods for the approximate solution of time dependent problems
Global Atmospheric Research Programme (Publication Series No. 10)
[60] Gottlieb D and Hesthaven J S 2001 Spectral methods for hyperbolic problems J. Comput. Appl.
Math. 128 83–131
	Introduction
	Constraint-preserving boundary conditions
	The generalized harmonic evolution system
	Construction of boundary conditions
	Improved gauge boundary condition
	Numerical results
	Alternate boundary conditions
	Freezing the incoming fields
	Sommerfeld boundary conditions
	Kreiss-Winicour boundary conditions
	Alternate approaches
	Spatial compactification
	Sponge layers
	Physical gravitational waves
	Difference of 4 with respect to the reference solution
	Comparison with the predicted reflection coefficient
	Discussion
	Details on the numerical test problem
	Initial data
	Numerical method
	Gauge source functions
	Error quantities
	Wave extraction
	Details of the alternate approaches
	Spatial compactification
	Sponge layers
