Testing outer boundary treatments for the Einstein equations

TL;DR

This paper compares various outer boundary treatments in numerical relativity, assessing their effectiveness in minimizing reflections and constraint violations using a Schwarzschild black hole test case.

Contribution

It introduces an improved boundary condition for gauge degrees of freedom and evaluates multiple boundary treatments within a unified framework.

Findings

Constraint-preserving boundary conditions perform well in reducing reflections.

Freezing the Newman-Penrose scalar Psi_0 minimizes boundary reflections.

Combining constraint preservation with gauge control yields the best results.

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 use a first-order generalized harmonic formulation of the Einstein equations. 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 Psi_0, and control gauge reflections.

Figures

Tables

Table 1. Table 1: Details of the filtering methods

No.	Type	Parameters	Applied to
1	Kreiss-Oliger	$ϵ=0.25$	right side
2	Kreiss-Oliger	$ϵ=0.25$	solution
3	Exponential	$\sigma =0.76$, $p=13$	right side
4	Exponential	$\sigma =0.76$, $p=13$	solution

Equations

$$\square x^a=H^a(x,\psi ),$$

$${\partial }_t{\psi }_{ab}\simeq 0,$$

$${\partial }_t{\Pi }_{ab}\simeq N^k{\partial }_k{\Pi }_{ab}-N{g}^{ki}{\partial }_k{\Phi }_{iab}-{\gamma }_2{N}^k{\partial }_k{\psi }_{ab},$$

$${\partial }_t{\Phi }_{iab}\simeq N^k{\partial }_k{\Phi }_{iab}-N{\partial }_i{\Pi }_{ab}+N{\gamma }_2{\partial }_i{\psi }_{ab},$$

$${\mathcal{C}}_{iab}\equiv {\partial }_i{\psi }_{ab}-{\Phi }_{iab}=0.$$

$${\mathcal{C}}_a\equiv -\square x_a+H_a={\psi }^{bc}{\Gamma }_{abc}+H_a=0.$$

$$u_{ab}^0={\psi }_{ab},$$

$$u_{ab}^{1\pm }={\Pi }_{ab}\pm {\Phi }_{nab}-{\gamma }_2{\psi }_{ab},$$

$$u_{Aab}^2={\Phi }_{Aab},$$

$${\tilde{u}}_{ab}^{1\pm }\equiv {\Pi }_{ab}\pm {\Phi }_{nab}.$$

$$\square {\mathcal{C}}_a=\text{(lower-order terms homogeneous in the constraints})$$

$$P_{ab}^{\mathrm{C}cd}{\partial }_n{u}_{cd}^{1-}$$

$${\mathcal{C}}_{ijab}\equiv -2{\partial }_{[i}{\Phi }_{j]ab}$$

$${\partial }_t{\Psi }_0\doteq 0,$$

$$P_{ab}^{\mathrm{P}cd}{\partial }_n{u}_{cd}^{1-}$$

$$P_{ab}^{\mathrm{G}\mathrm{cd}}{\partial }_t{u}_{cd}^{1-}\doteq 0,$$

$$x^a\to x^a+{\xi }^a$$

$$\square {\xi }^a=0.$$

$$\delta {\psi }_{ab}=-2{\partial }_{(a}{\xi }_{b)}.$$

$$\square \psi =0.$$

$$\psi =\sum _{l=1}^{\infty }\sum _{m=-l}^l{Y}_{lm}(\theta ,\varphi ){\psi }_l(t,r),$$

$${\psi }_l^{\pm }(t,r)=r^{l-1}{\left(\frac{\partial }{\partial r}\frac{1}{r}\right)}^l{F}_l^{\pm }(r\mp t),$$

$$({\partial }_t+{\partial }_r+r^{-1})\psi \doteq 0.$$

$$({\partial }_t+{\partial }_r+{\gamma }_2)\psi \doteq 0,$$

$${\psi }_l={\psi }_l^{+}+{\rho }_l{\psi }_l^{-}$$

$$F_l^{\pm }(x)={\text{rme}}^{\pm ikx},$$

$$P_{ab}^{\mathrm{G}\mathrm{cd}}{\partial }_t[u_{cd}^{1-}+({\gamma }_2-r^{-1}){\psi }_{cd}]\doteq 0.$$

$${\partial }_t{u}_{ab}^{1-}\doteq 0(\mathrm{and}{\partial }_t{u}_{Aab}^2\doteq 0\mathrm{if}{N}^n\dot{>}0).$$

$$(\partial_{t}+\partial_{r}+r^{-1})\left(\begin{array}[]{c}g_{ij}-\delta_{ij}\\ K_{ij}\end{array}\right)\doteq 0.$$

$$({\partial }_t+{\partial }_r+r^{-1})({\psi }_{ab}-{\eta }_{ab})\doteq 0$$

$${\partial }_t[u_{ab}^{1-}+({\gamma }_2-r^{-1}){\psi }_{ab}]\doteq 0.$$

$${\mathcal{C}}_a\doteq 0.$$

Full text

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 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 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 $\Psi_{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 $\Psi_{0}$.

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.

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 $\Psi_{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 $\Psi_{4}$ from the reference solution. We also compare the measured reflections caused by our $\Psi_{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 $x^{a}$ obey the wave equation

[TABLE]

where $\Box=\psi^{ab}(\partial_{a}\partial_{b}-\Gamma^{c}{}_{ab}\partial_{c})$ is the covariant scalar wave operator, with $\psi_{ab}$ the spacetime metric and $\Gamma^{c}{}_{ab}$ the associated metric connection. In this formulation of the Einstein system the gauge source function $H^{a}$ may be chosen freely as a function of the coordinates and of the spacetime metric $\psi_{ab}$ (but not derivatives of $\psi_{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 $\Phi_{iab}\equiv\partial_{i}\psi_{ab}$ and $\Pi_{ab}\equiv-t^{c}\partial_{c}\psi_{ab}$, where $t^{c}$ is the future directed unit normal to the $t=\mathrm{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 111The parameter $\gamma_{1}$ of [6] is chosen to be $-1$, which ensures that the equations are linearly degenerate.

[TABLE]

where $\simeq$ indicates that purely algebraic terms have been omitted, $g_{ij}$ is the spatial metric of the $t=\mathrm{const.}$ slices, and $N$ and $N^{i}$ are the lapse function and shift vector, respectively. The parameter $\gamma_{2}$ was introduced in [6] in order to damp violations of the three-index constraint

[TABLE]

We also include terms of lower derivative order that are designed to damp violations of the harmonic gauge constraint [20]

[TABLE]

The system (2) is symmetric hyperbolic. The characteristic fields in the direction $n_{i}$ (where $n_{a}t^{a}=0$) are given by

[TABLE]

For future reference, we also define

[TABLE]

Here and in the following, an index $n$ denotes contraction with $n_{i}$, while upper-case Latin indices $A,B,\ldots$ are orthogonal to $n$, e.g. $v_{A}=P_{Ai}v^{i}$ where $P_{ab}\equiv\psi_{ab}-n_{a}n_{b}+t_{a}t_{b}$. 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

[TABLE]

and the corresponding incoming fields will involve first derivatives of $\mathcal{C}_{a}$. In terms of the incoming modes $u^{1-}_{ab}$ (6) of the main evolution equations, the resulting constraint-preserving boundary conditions can be written in the form

[TABLE]

where $P^{\mathrm{C}}$ is a projection operator of rank $4$ (cf. [6]). Here $n_{i}$ now refers to the outward-pointing unit spatial normal to the boundary, $l^{a}=(t^{a}+n^{a})/\sqrt{2}$, $k^{a}=(t^{a}-n^{a})/\sqrt{2}$, and $\doteq$ denotes equality at the boundary. If the shift vector points towards the exterior at the boundary ($N^{n}\,\dot{>}\,0$), the fields $u^{2}_{Aab}$ (7) are incoming as well and we obtain a boundary condition on them by requiring the components $\mathcal{C}_{nAab}$ of the four-index constraint

[TABLE]

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 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 $\Psi_{0}$ (evaluated for a Newman-Penrose null tetrad containing the vectors $l^{a}$ and $k^{a}$). Hence the physical boundary condition we use is [27, 22, 30, 29, 31]

[TABLE]

which can be written in a form similar to (10),

[TABLE]

Here $P^{\mathrm{P}}$ is a projection operator of rank $2$ that is orthogonal to $P^{\mathrm{C}}$ [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 $u^{1-}_{ab}$. The remaining four components correspond to gauge degrees of freedom. In the past we chose simply to freeze those components in time [6],

[TABLE]

where $P^{\mathrm{G}}\equiv\mathbb{I}-P^{\mathrm{C}}-P^{\mathrm{P}}$.

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 ($H_{a}=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\gg 1$ and $R\gg 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

[TABLE]

are still allowed provided the displacement vector satisfies the wave equation,

[TABLE]

Under such a coordinate transformation, the metric changes by

[TABLE]

A closer inspection [32] of the projection operator $P^{\mathrm{G}}$ in (14) shows that the gauge boundary conditions control the components $l^{a}\delta\psi_{ab}$ of the perturbations, where $l^{a}\equiv(t^{a}+n^{a})/\sqrt{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 $l^{a}\delta\psi_{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 $l^{a}\delta\psi_{ab}$ obeys the scalar wave equation

[TABLE]

Solutions to this equation can be written in the form

[TABLE]

where the $Y_{lm}$ are the standard spherical harmonics and the $\psi_{l}$ are linear combinations of outgoing ($+$) and incoming ($-$) solutions

[TABLE]

$F^{\pm}_{l}(x)$ being arbitrary functions. A boundary condition is needed on $\psi$ that eliminates the incoming part of these solutions. In [36], a hierarchy of boundary conditions is constructed that accomplish this task for all $l\leqslant 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 222 To avoid confusion, we remark that in [5, 14], the term ‘Sommerfeld condition’ is used in reference to a condition of the form $(\partial_{t}+\partial_{r})u\doteq 0$, i.e. without the extra $r^{-1}$ term due to our polar coordinates.

[TABLE]

In contrast, our old gauge boundary condition that froze the incoming characteristic field, as in (14), is given by

[TABLE]

where $\gamma_{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

[TABLE]

with generating functions

[TABLE]

where $k\in\mathbb{R}$ is the wave number. Substituting this solution into the boundary conditions (21) resp. (22), we solve for the reflection coefficient $\rho_{l}$. Figure 1 shows $|\rho_{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\approx 1.6/M$ and in most cases, we place the outer boundary at $R=41.9M$.) We see that $|\rho_{l}|$ is much smaller (by about $3$ orders of magnitude) for the Sommerfeld condition than for the freezing condition.

In the notation of the previous subsection, the improved gauge boundary condition (21) reads (after taking a time derivative),

[TABLE]

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 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^{\infty}$ and $L^{2}$ norms of the difference $\Delta\mathcal{U}$ between the numerical solution and the reference solution, and also the violations of the constraints $\mathcal{C}$ (see 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 $\Delta\mathcal{U}$ and $\mathcal{C}$ are computed.

In the difference $\Delta\mathcal{U}$ we see a reflection that originates when the wave reaches the boundary at $t\approx R$ and then amplifies as it moves inward in the spherical geometry, assuming its maximum at $t\approx 2R$. This feature is much more prominent in the $L^{\infty}$ norm than in the $L^{2}$ norm, which is why we display only the $L^{\infty}$ norm in subsequent plots. The reflection is much smaller (by a factor of $\approx 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 $\Delta\mathcal{U}$, which in contrast to the old conditions appears to decrease as resolution is increased.

We would like to stress that $\Delta\mathcal{U}$ is a coordinate dependent quantity. Hence a smaller $\Delta\mathcal{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 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 $||\Delta\mathcal{U}||_{\infty}$ 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 $\Psi_{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.

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

[TABLE]

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.

The left side of figure 4 demonstrates that freezing boundary conditions cause a significantly larger (by $\approx 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 $g_{ij}$ and extrinsic curvature $K_{ij}$,

[TABLE]

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 $g_{ij}$ and $K_{ij}$, and the incoming characteristic fields of our generalized harmonic formulation of Einstein’s equations. Instead we consider the similar condition

[TABLE]

on all the components of the spacetime metric ($\eta_{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

[TABLE]

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 $u^{2}_{Aab}$ 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 $\approx 2$ orders of magnitude, but later $||\Delta\mathcal{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,

[TABLE]

In our notation, this can be written as an algebraic condition on part of the incoming fields $u^{1-}$,

[TABLE]

where

[TABLE]

The range of the projection operator $P^{\mathrm{C^{\prime}}}$ is identical with that of $P^{\mathrm{C}}$ defined in (10). For the unconstrained incoming fields $\tilde{u}^{1-}$ (i.e. $u^{1-}$ without the $\gamma_{2}$ term, equation (8)), Kreiss and Winicour [14] specify certain free boundary data $q^{\mathrm{P}}_{ab}$ and $q^{\mathrm{G}}_{ab}$. In our notation,

[TABLE]

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 $u^{2}_{Aab}$ to emulate the second-order formulation of [5, 14], and this value of $u^{2}_{Aab}$ is then used to compute $F_{a}$ in (3.3).

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 $\partial_{t}\Psi_{0}\doteq 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 have used two such mappings, named Tan and Inverse, as detailed in B.1. Each map has a scale $R$ across which the mapping is (essentially) linear. The outermost grid point is placed at a very large but finite uncompactified radius ($r=10^{17}M$). 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 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\leqslant 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 1) 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), $||\Delta\mathcal{U}||$ is comparable to constraint-preserving boundary conditions at the peak of reflections (at $t\approx 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 $||\Delta\mathcal{U}||$ initially but the convergence degrades at later times. This is surprising at first because with increasing resolution, 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 B.1. The results are similar, but somewhat worse than the Tan mapping results shown here.

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

[TABLE]

where $u_{0}$ refers to the unperturbed background solution (Schwarzschild spacetime in our case) and the smooth sponge function $\gamma(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 $\gamma(r)$ can be found in B.2.

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\lesssim 2R$), the $||\Delta\mathcal{U}||_{\infty}$ 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 $\Psi_{4}$ and $\Psi_{0}$. The scalar $\Psi_{4}$ is dominated by the outgoing radiation (its ingoing part is suppressed by a factor of $(kr)^{4}$, where $k$ is the wavenumber), whereas $\Psi_{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=R_{\mathrm{ex}}$, using the methods described in A.5.

We note that $\Psi_{4}$ ($\Psi_{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 $\Or(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 $\Psi_{4}$ caused by the various boundary treatments.

5.1 Difference of $\Psi_{4}$ with respect to the reference solution

We begin by evaluating $\Delta\Psi_{4}\equiv\Psi_{4}-\Psi_{4}^{\mathrm{ref}}$, where $\Psi_{4}$ is the Newman-Penrose scalar computed using one of the various boundary methods and $\Psi_{4}^{\mathrm{ref}}$ 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 $|\Delta\Psi_{4}|$ over time intervals of length 20$M$ (this time filtering averages over the high frequency quasi-normal oscillations of the black hole), normalized by the maximum value of $|\Delta\Psi_{4}|$ over the entire evolution. The 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 $R_{\mathrm{ex}}=40M$.

The first peak in $|\Delta\Psi_{4}|$ seen in figure 10 arises as the wave in our test problem passes outward through the extraction sphere at $t\approx R_{\mathrm{ex}}$. 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 $\Delta\Psi_{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). 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 $|\Delta\Psi_{4}|$, particularly at early times $t\sim 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 $|\Delta\Psi_{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, $R_{\mathrm{ex}}=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 $\Psi_{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 $\rho$ (defined as the ratio of the ingoing to the outgoing parts of the solution) that arises from the $\partial_{t}\Psi_{0}\doteq 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],

[TABLE]

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 $\Or(M/R)$.) By evaluating $\Psi_{0}$ and $\Psi_{4}$ at the extraction radius of our test, we find that the ratio $\Psi_{0}/\Psi_{4}$ agrees with their predicted $\rho$ to leading order in $1/(kR)$. We note that the tetrad we use for wave extraction (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 $\Psi_{0}$ and $\Psi_{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 $\Psi_{0}(t)$ and $\Psi_{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 $\frac{3}{2}(kR)^{-4}\Psi_{4}$, which by the above argument should agree with $\Psi_{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 $\rho$ on $k$ and $R$ is captured very well. We surmise that the levelling off of our numerical $\Psi_{0}$ for $k\gtrsim 3$ is due to numerical roundoff effects. (Note the magnitude of $\Psi_{0}$ at those frequencies.) For radii $R\gtrsim 200M$, $\Psi_{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 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 $\Psi_{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 $\Delta\mathcal{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 $\Delta\mathcal{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 $\Psi_{4}$ for this purpose. Here our conclusions are somewhat different. Rather surprisingly, most of the boundary methods we consider generate relatively small errors in $\Psi_{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 $\sim 10^{-5}$.) The largest relative errors in $\Psi_{4}$ we find ($\sim 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 ($\lambda\sim 20-100M$) than ours ($\lambda\sim 4M$). Our tests suggest that the errors in $\Psi_{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 $\Psi_{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 $\lambda/R\gtrsim 0.5$) than for the simple test problem considered here (where $\lambda/R\sim 0.1$). Boundary treatments generally work better for smaller $\lambda/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-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

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,

[TABLE]

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-r_{0})^{2}/w^{2}]$ with $A=4\times 10^{-3}$, $r_{0}=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 $H_{a}$; we choose $\partial_{t}N=0$ and $\partial_{t}N^{i}=0$ at $t=0$.

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 $\Delta 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 (where each Cartesian tensor component is expanded in the standard scalar spherical harmonics). Typical resolutions are $N_{r}\in\{21,31,41,51\}$ coefficients per subdomain for the Chebyshev series and $l\leqslant L$ with $L\in\{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 $\Delta 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 $\Delta t/\Delta x_{\mathrm{min}}$ of at most $2.25$, where $\Delta x_{\mathrm{min}}$ 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.

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: $H_{a}=H_{a}(t,x,\psi)$. 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 $H_{a}$ via the constraint (4). For the subsequent evolution, we hold this $H_{a}$ fixed in time.

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 $(\psi_{ab},\Pi_{ab},\Phi_{iab})$, the difference between that solution and the reference solution $(\psi_{ab}^{\mathrm{(ref)}},\Pi_{ab}^{\mathrm{(ref)}},\Phi_{iab}^{\mathrm{(ref)}})$ is computed with the following norm at each point in space,

[TABLE]

where $\Delta\psi_{ab}$ means $\psi_{ab}-\psi_{ab}^{\mathrm{(ref)}}$, and similarly for $\Delta\Pi_{ab}$ and $\Delta\Phi_{iab}$. Second, we define a quantity $\mathcal{C}$ that measures the violations in all of the constraints of our system,

[TABLE]

where $\mathcal{F}_{a}$ and $\mathcal{C}_{ia}$ are first derivatives of $\mathcal{C}_{a}$ defined in [6]. To compute global error measures, a spatial norm $||\cdot||$, either the $L^{\infty}$ norm or the $L^{2}$ norm, is applied separately to $\Delta\mathcal{U}$ and $\mathcal{C}$.

The question often arises as to the significance of particular values of $||\Delta\mathcal{U}||$ and $||\mathcal{C}||$. For example, is a simulation with $||\mathcal{C}||=10^{-2}$ good to one percent accuracy? To make it easier to answer such questions, we normalize both $||\Delta\mathcal{U}||$ and $||\mathcal{C}||$ as follows, and we always plot normalized quantities.

We divide $||\Delta\mathcal{U}||$ by a normalization factor $||\Delta\mathcal{U}_{0}||$, defined as the difference between a given solution at $t=0$ and the unperturbed Schwarzschild background; i.e., the quantity $||\Delta\mathcal{U}_{0}||$ is computed from (37) using the unperturbed Schwarzschild solution instead of the reference solution. Since $||\Delta\mathcal{U}_{0}||$ 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, $||\Delta\mathcal{U}_{0}||_{\infty}=6\times 10^{-3}$ and $||\Delta\mathcal{U}_{0}||_{2}=1.4\times 10^{-4}$. The quantity $||\Delta\mathcal{U}||/||\Delta\mathcal{U}_{0}||$ is more easily interpreted than $||\Delta\mathcal{U}||$; for example, $||\Delta\mathcal{U}||/||\Delta\mathcal{U}_{0}||$ is unity when the difference from the reference solution is of the same size as the initial perturbation.

Similarly, the constraint energy norm $||\mathcal{C}||$ is divided by the norm of the first derivatives $||\partial\mathcal{U}||$ (at the respective time),

[TABLE]

The constraints for our system are linear combinations of the first derivatives of the fields, hence $||\mathcal{C}||/||\partial\mathcal{U}||\sim 1$ corresponds to a complete violation of the constraints.

A.5 Wave extraction

For evaluating gravitational waveforms, we compute the Newman-Penrose scalars

[TABLE]

where $C_{abcd}$ is the Weyl tensor, $l^{a}$ and $k^{a}$ are outgoing and ingoing null vectors normalized according to $l^{a}k_{a}=-1$, $m^{a}$ is a complex unit null spatial vector orthogonal to $l^{a}$ and $k^{a}$, and $\bar{m}^{a}$ is the complex conjugate of $m^{a}$. For perturbations of flat spacetime, there is a standard choice for the vectors $l^{a}$, $k^{a}$, and $m^{a}$. In general curved spacetimes, however, no such prescription for the tetrad exists that would produce coordinate-independent quantities $\Psi_{0}$ and $\Psi_{4}$ at finite radius. We choose the null vectors according to

[TABLE]

where $t^{a}$ is the future-pointing unit timelike normal to the $t=\mathrm{const.}$ slices and $n^{a}$ is the unit spacelike normal to the extraction sphere. Finally, we choose

[TABLE]

where $(r,\theta,\phi)$ are spherical coordinates on the $r=R_{\mathrm{ex}}=\mathrm{const.}$ extraction sphere. Note that our choice of $m^{a}$ 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 $\Psi_{0}$ and $\Psi_{4}$ because of the lack of tetrad orthonormality will be second-order small in perturbation theory.

The quantity $\Psi_{4}$ corresponds to outgoing radiation in the limit of $r\rightarrow\infty$, $t-r=\mathrm{const.}$, i.e. as future null infinity is approached. Similarly $\Psi_{0}$ corresponds to ingoing radiation as past null infinity is approached. At finite extraction radius, $\Psi_{4}$ and $\Psi_{0}$ will disagree with the waveforms observed at infinity by terms of the order $\Or(R_{\mathrm{ex}})^{-1}$.

We decompose the quantities $\Psi_{4}$ and $\Psi_{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 $\Psi_{4}$, and we only display that mode in our plots. We normalize the curves in our graphs by the maximum (in time) value of $|\Psi_{4}|$ at the extraction radius $R_{\mathrm{ex}}$, which for $R_{\mathrm{ex}}=40M$ is $\max|\Psi_{4}|=6\times 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.

B.1 Spatial compactification

We implement spatial compactification by introducing a radial coordinate transformation $x\rightarrow r(x)$ that maps a compact ball on the computational grid with $x\in[0,x_{\mathrm{max}}]$ to the full unbounded physical slice with $r\in[0,\infty]$. We consider two such mappings. The Tan mapping is similar to the one used by Pretorius [15, 16, 17] and is given by

[TABLE]

The scale $R$ determines the range in physical radius $r$ across which the map is essentially linear (see figure 12). 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

[TABLE]

see figure 12. This map is only $C^{1}$ 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 A.2 for details on the pseudospectral method we use.) Higher values of $k$ correspond to shorter wavelengths in the numerical approximation; let $k_{\mathrm{max}}$ 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,

[TABLE]

Typical values of the parameter $\epsilon$ used by Pretorius are $\epsilon\in[0.2,\,0.5]$; we use $\epsilon=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 $x_{j}$. (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

[TABLE]

at each time step, where $h$ is the grid spacing and $D^{4}$ is the second-order accurate centred finite difference operator approximating the fourth derivative,

[TABLE]

Taking $u$ to be a Fourier mode $u^{(k)}_{j}=\exp(\rmi kx_{j})$, it follows that the mode is damped by a frequency-dependent factor,

[TABLE]

where $k_{\mathrm{max}}=\pi/(2h)$ is the Nyquist frequency. Thus we obtain the filter function (45). 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 (45) to our Chebyshev expansion coefficients. Note that in (48), 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),

[TABLE]

Typical values of the parameters are $\sigma=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(k_{\mathrm{max}})\approx 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 $\partial_{t}u=S$ are modified according to $\partial_{t}u=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. 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 1.

B.2 Sponge layers

For sponge layers we must specify a sponge profile function $\gamma(r)$, as defined in (34). We choose $\gamma(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 $\gamma(r)$ we use is a Gaussian centred at the outer boundary, which we choose to place at $r=3R$,

[TABLE]

The amplitude of the Gaussian is taken to be $\gamma_{0}=1$. The width $\sigma$ is chosen so that $\gamma(r)\leqslant 10^{-16}$ (the numerical roundoff error) for $r\leqslant R$, which requires $\sigma\lesssim R/3$. In our numerical example, we take $R=41.9M$ and $\sigma=13.3M$. Hence $\sigma$ is considerably larger than the wavelength $\lambda\approx 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 13 shows a plot of this sponge profile.

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 PHY-0099568, PHY-0244906, PHY-0601459, DMS-0553302 and NASA grants NAG5-12834, NNG05GG52G.