The numerical initial boundary value problem for the generalized conformal field equations

TL;DR

This paper develops a numerical method for solving the initial boundary value problem of generalized conformal field equations, enabling long-term simulations of black hole spacetimes while maintaining stability and constraints.

Contribution

It introduces a well-structured formulation using space-spinor formalism and boundary conditions that preserve constraints, suitable for long-time evolution of black hole perturbations.

Findings

Successfully implemented a numerically well-posed system

Demonstrated stability with Minkowski perturbations

Extended applicability to Schwarzschild black hole perturbations

Abstract

In this paper we study a numerical implementation for the initial boundary value formulation for the generalized conformal field equations. We propose a formulation which is well suited for the study of the long-time behaviour of perturbed exact solutions such as a Schwarzschild or even a Kerr black hole. We describe the derivation of the implemented equations which we give in terms of the space-spinor formalism. We discuss the conformal Gauss gauge, and a slight generalization thereof which seems to be particularly useful in the presence of boundaries. We discuss the structure of the equations at the boundary and propose a method for imposing boundary conditions which allow the correct number of degrees of freedom to be freely specified while still preserving the constraints. We show that this implementation yields a numerically well-posed system by testing it on a simple case of gravitational perturbations of Minkowski space-time and subsequently with gravitational perturbations of Schwarzschild space-time.

Equations

$${\tilde{R}}_{ab}=\lambda {\tilde{g}}_{ab},$$

$${\hat{\nabla }}_a{\tilde{g}}_{bc}$$

$${\hat{\nabla }}_a{g}_{bc}$$

$$f_a=b_a-{\Theta }^{-1}{\nabla }_a\Theta .$$

$$h_a:=\Theta b_a=f_a+{\nabla }_a\Theta .$$

$$({\hat{\nabla }}_a-{\nabla }_a){\omega }_b=-f_{ab}{}^c{\omega }_c,$$

$$({\hat{\nabla }}_a-{\nabla }_a)v^c=f_{ab}{}^c{v}^b,$$

$$f_{ab}{}^c={\delta }^c{}_a{f}_b+{\delta }^c{}_b{f}_a-g_{ab}{g}^{cd}{f}_d.$$

$$c_{\mathbf{a}}^{\mu }:=e_{\mathbf{a}}(x^{\mu })={\hat{\nabla }}_{\mathbf{a}}{x}^{\mu },{\hat{\nabla }}_{\mathbf{a}}{e}_{\mathbf{b}}={\hat{\Gamma }}_{\mathbf{ab}}{}^{\mathbf{c}}{e}_{\mathbf{c}}.$$

$$f_{\mathbf{a}}=\frac{1}{4}{\hat{\Gamma }}_{\mathbf{ab}}{}^{\mathbf{b}}.$$

$$\Big{[}\widehat{\nabla}_{a},\widehat{\nabla}_{b}\Big{]}x^{\mu}=0,$$

$$\Big{[}\widehat{\nabla}_{a},\widehat{\nabla}_{b}\Big{]}e_{\mathbf{c}}^{d}=\widehat{R}_{abc}{}^{d}e_{\mathbf{c}}^{c}.$$

$$\widehat{R}_{abc}{}^{d}=2\Big{(}\delta^{d}{}_{[a}\widehat{P}_{b]c}-\delta^{d}{}_{c}\widehat{P}_{[ab]}-g_{c[a}\widehat{P}_{b]}{}^{d}\Big{)}+C_{abc}{}^{d},$$

$$e_{\mathbf{a}}(c_{\mathbf{b}}^{\mu })-e_{\mathbf{b}}(c_{\mathbf{a}}^{\mu })$$

$$e_{\mathbf{a}}({\hat{\Gamma }}_{\mathbf{bc}}{}^{\mathbf{d}})-e_{\mathbf{b}}({\hat{\Gamma }}_{\mathbf{ac}}{}^{\mathbf{d}})$$

$$+C_{\mathbf{abc}}{}^{\mathbf{d}}-2{\eta }_{\mathbf{c}[\mathbf{a}}{\hat{P}}_{\mathbf{b}]}{}^{\mathbf{d}}+2{\delta }_{[\mathbf{a}}{}^{\mathbf{d}}{\hat{P}}_{\mathbf{b}]\mathbf{c}}-2{\hat{P}}_{[\mathbf{ab}]}{\delta }_{\mathbf{c}}{}^{\mathbf{d}}.$$

$${\tilde{\nabla }}_{[e}{\tilde{R}}_{ab]c}{}^d={\tilde{\nabla }}_{[e}{C}_{ab]c}{}^d=0.$$

$${\nabla }_e{K}_{abc}{}^e=0.$$

$${\hat{\nabla }}_e{K}_{abc}{}^e=f_e{K}_{abc}{}^e.$$

$${\hat{\nabla }}_{[e}{\hat{R}}_{ab]c}{}^d=0,$$

$${\hat{\nabla }}_a{\hat{P}}_{bc}-{\hat{\nabla }}_b{\hat{P}}_{ac}=\left({\hat{\nabla }}_e\Theta +\Theta f_e\right)K_{abc}{}^e=h_e{K}_{abc}{}^e.$$

$$u^b{\nabla }_b{u}^a$$

$$u^b{\nabla }_b{h}_a$$

$${\hat{\nabla }}_a{g}_{bc}={\nabla }_a{g}_{bc}-2f_a{g}_{bc}.$$

$$u^b{\hat{\nabla }}_b{u}^a$$

$$u^a{\hat{P}}_{ab}$$

$$0={\hat{\nabla }}_u{e}_{\mathbf{c}}^c=u^a{\tilde{\nabla }}_a{e}_{\mathbf{c}}^c+(b_a{e}_{\mathbf{c}}^a)u^b+(b_a{u}^a)e_{\mathbf{c}}^b-({\tilde{g}}_{ab}{u}^a{e}_{\mathbf{c}}^b){\tilde{g}}^{ec}{b}_e.$$

$${\hat{\nabla }}_u{e}_{\mathbf{c}}^c={\omega }^{\mathbf{b}}{}_{\mathbf{c}}{e}_{\mathbf{b}}^c,\text{for }\mathbf{c}=1,2,3$$

$$g(e_{\mathbf{a}},e_{\mathbf{b}})={\eta }_{\mathbf{ab}}={\Theta }^2\tilde{g}(e_{\mathbf{a}},e_{\mathbf{b}}).$$

$$\dot{\Theta }=(b_e{u}^e)\Theta ,$$

$$\stackrel{...}{\Theta }=0.$$

$$Z=\Theta b_a{u}^a=\dot{\Theta },H={\Theta }^{-1}({\tilde{g}}^{ab}{b}_a{b}_b-\frac{1}{3}\lambda )=2\ddot{\Theta },$$

$$\Theta (s)=\underline{\Theta }+\underline{Z}s+\frac{1}{4}\underline{H}{s}^2.$$

Full text

I Introduction

In the 1960’s Penrose introduced the idea of studying the asymptotic properties of a space-time from the perspective of its conformal (light cone) structure Penrose (1963, 1962, 1964). The idea is to conformally embed the physical space-time in question into a conformally related space-time, usually called the conformal space-time, where its image becomes an open submanifold. The boundary of this submanifold in the conformal space-time is referred to as the conformal boundary, (denoted $\mathscr{I}$), and represents the points at infinity of the physical space-time. A very nice property of this picture is that from the point of view of the conformal space-time, the conformal boundary, and hence the asymptotic structure of the physical space-time, can be investigated using local differential geometry.

Although Penrose constructed the conformal completion of many known space-times the question as to whether there exist large classes of space-times which admit a conformal completion was not resolved over a long period of time. Over the years Friedrich has worked on this question by developing the conformal field equations, which are the mathematical extension of the Einstein equations to the conformal manifold. The Generalized Conformal Field Equations (GCFE) are the latest form of these equations, and were first written down in 1995 in the important work Friedrich (1995). It was only after the surprising result by Corvino and Schoen Corvino and Schoen (2006); Corvino (2000) on the construction of initial data which agree with exact Schwarzschild or Kerr data outside a compact region that the question of existence raised above could be answered in the affirmative by Chruściel and Delay Chruściel and Delay (2002).

With the advent of gravitational wave detectors and, even more so, after the successful detection of (so far) three gravitational wave events Abbott et al. (2016a, b, 2017), the calculation of asymptotic quantities has become a major task for the numerical relativity community. After all, the waveforms used by the detectors as templates to match to observational data are strictly defined only on part of the conformal boundary.

In the numerical relativity community it is standard practice to calculate the waveforms in the so-called “wave-zone”, where the space-time is almost flat and where a decent approximation can be made, and not on the conformal boundary itself (with the exception of, for example, characteristic extraction Reisswig et al. (2010)). The Cauchy-perturbative approach and extrapolation method are among the standard ways of calculating the waveforms in the wave-zone (see Bishop and Rezzolla (2016) and the references within for a comprehensive review of these methods). These rely on solving the full non-linear Einstein equations in the portion of the space-time that has the most dynamics (for example inside a box surrounding a binary black hole system) generally using the BSSN formulation Shibata and Nakamura (1995); Baumgarte and Shapiro (1998), the Generalized Harmonic Formalism Pretorius (2005) or the Z4 formalulations Bona et al. (2003); Alic et al. (2012). These approaches usually solve the Einstein equations in the physical space-time in the form of an Initial Value Problem (IVP) or an Initial Boundary Value Problem (IBVP) with a time-like outer boundary that approaches future time-like infinity.

The main difference between the conformal field equations over the standard approaches is that they grant access to the entire conformal boundary, while the other methods remain in the physical space-time or, in some cases, extend to a finite portion of null infinity. With the standard approaches it is principally impossible to explore the global structure of a space-time completely since one is always confined to a finite portion of the physical space-time. Thus, one could not, for example, study the behaviour of fields across space-like infinity using Friedrich’s cylinder interpretation Friedrich (1998) (see the summary V for more applications). Yet, even with its advantages, there have been comparatively few investigations into the conformal field equation’s potential for numerically studying global properties of asymptotically simple space-times, perhaps due to their apparent size and complexity.

An early version of these equations — nowadays referred to as the metric conformal field equations — was studied in several papers in Frauendiener (1998a, b, 2000); Frauendiener and Hein (2002) and in Hübner (2001, 1999a, 1999b, 1996a, 1996b). It could be demonstrated that the conformal approach does indeed deliver on its promise: one could evolve the global space-time within a finite computational time on a finite grid and determine the gravitational wave-forms with high accuracy. However the downside of this particular formulation was that the conformal factor is an unknown in the system, being evolved alongside other quantities, and hence the location of the conformal boundary is not known a priori, even though it was possible to devise suitable “$\mathscr{I}$-freezing” gauges.

The only numerical implementations of the generalized conformal field equations to date are due to Zenginoğlu Zenginoğlu (2007a, b) who reproduced the Schwarzschild and Kerr space-times in spherical and axi-symmetric symmetry, and by Beyer Beyer (2007, 2009a, 2008, 2009b) who considered the associated Initial Value Problem (IVP) for investigating spatially compact space-times with a positive cosmological constant $\lambda$.

In this paper we want to study a numerical implementation of the GCFE in the context of an initial boundary value problem. We consider this approach for several reasons. Any numerical formulation of an evolution problem is based on a finite computational domain so, quite generally, some sort of conditions will have to be imposed on its boundary. Even when the mathematically underlying problem does not need any boundary conditions — this will be the case when the problem is based on spatially compact manifolds or when the boundary is a total characteristic such as when considering the finite initial boundary value problem at space-like infinity Friedrich (1998) — the finite size of the computational domain will lead to a numerical IBVP in the sense that the boundary points will need special consideration.

Our approach complements the more traditional hyperboloidal initial value problem pursued in the context of the metric conformal field equations and the finite initial boundary value problem just mentioned. The hyperboloidal IVP is based on space-like hyper-surfaces which extend to null-infinity. Initial data are provided by solving the constraint equations associated with the conformal field equations. By construction the hyperboloidal IVP is semi-global in its nature in the sense that its solution covers the global future of the space-time but only a part of its past, see Fig. 1.

The finite IBVP at space-like infinity is based on asymptotically Euclidean hyper-surfaces which extend out to space-like infinity represented as a cylinder (see Beyer et al. (2012); Frauendiener and Hennig (2014, 2017); Doulis and Frauendiener (2013, 2017) for some studies of this approach in very simplified situations). Initial data are obtained by solving the conformal constraint equations and the solution is global covering the entire past and future of the space-time evolving from the initial data, see Fig. 1.

Our present approach was chosen for its very physical simplicity: we specify the initial data of a known solution exactly and perturb it by incoming radiation. More precisely, we specify initial data on a space-like hyper-surface with boundary which is not assumed to extend to infinity. Therefore, we need to specify boundary conditions in order to make the problem well-posed. If the problem is indeed well-posed (and in this paper we argue that it is at least numerically well-posed) then this setup allows us to study perturbations of an exact solution in a clear physical way. We choose initial data corresponding to the exact solution and specify the perturbation by boundary conditions. In this way we disentangle the object of the disturbance from the cause of the disturbance. Thus, we avoid the solution of the constraint equations which is known to introduce uncontrolled radiative degrees of freedom. Instead, we can specify in a very controlled way the incoming degrees of freedom on the time-like boundary. The solution obtained in this way is time-global, i.e. global in time but not in space in the sense that there is a region near space-like infinity which is not covered by the evolution of the initial data, see Fig. 1.

The GCFE can also be used to evolve asymptotically anti-de Sitter like space-times which naturally leads to an IBVP with boundary conditions to be specified on the conformal boundary Friedrich (1995). Finally, an IBVP is advantageous for studying perturbations of an exact solution as we will do here: the initial conditions are provided by the exact solution and the perturbation is injected from the boundary. Proceeding in this way avoids the necessity to solve the constraint equation for perturbed initial conditions. Hence the main focus of this work is to conduct the first numerical investigations into the viability of a well-posed IBVP framework for the GCFE where the constraints stay satisfied on at least the numerical level.

The structure of this paper is as follows: In section II we summarize the derivation of the GCFE and introduce the conformal Gauss gauge. We describe how we impose Newman and Penrose’s $\eth$-calculus Newman and Penrose (1966); Goldberg et al. (1967); Eastwood and Tod (1982); Penrose and Rindler (1986) to obtain proper spin-weighted equations and end the section with our derivation of the GCFE in the space-spinor formalism. In section III the field equation for the gravitational spinor is analysed and a numerical procedure for imposing constraint preserving boundary conditions is presented. We then test our system numerically using as applications: the IBVP for gravitational perturbations of Minkowski space-time in section IV.1 and the IBVP for gravitational perturbations of Schwarzschild space-time in section IV.2. The paper concludes in section V with a brief summary and discussion of future applications.

We use the conventions of Penrose and Rindler Penrose and Rindler (1986, 1988) throughout. In particular, we use the metric signature $(+,-,-,-)$ and define associated with any torsion-free connection $\nabla$: the Riemann tensor $[\nabla_{a},\nabla_{b}]\alpha_{c}=-R_{abc}{}^{d}\alpha_{d}$, the Ricci tensor $R_{ab}:=R_{acb}{}^{c}$, the Ricci scalar $R:=R_{a}{}^{a}$ and the Schouten tensor as P_{ab}:=-\frac{1}{4}R_{[ab]}-\frac{1}{2}\Big{(}R_{(ab)}-\frac{1}{6}Rg_{ab}\Big{)}, where $R_{[ab]}=0$ if the connection is Levi-Civita. Then the vacuum Einstein field equations take the form $R_{ab}=\lambda g_{ab}$ with cosmological constant $\lambda$. We will use the convention of denoting quantities associated to the physical metric $\widetilde{g}_{ab}$ with a $\;\widetilde{}\;$ and use bold Latin indices to denote frame indices $e^{a}_{\mathbf{i}}$. Transvections of space-spinors with the spin-frame $\{o,\iota\}$ are written as $\alpha_{0}:=\alpha_{A}o^{A}$ and $\alpha_{1}:=\alpha_{A}\iota^{A}$. If a space-spinor is symmetric in a certain number of indices, then we label their components as the number of contractions with $\iota$. For example if $K_{AB}=K_{(AB)}$ then its components are denoted by $K_{0},\;K_{1}$ and $K_{2}$.

II The analytical background

In this section we provide a brief derivation of the GCFE and the conformal Gauss gauge, a particularly useful set of gauge conditions for coordinates, tetrads and conformal factor. We discuss our imposition of the $\eth$-calculus and then use the space-spinor formalism to split the GCFE into evolution and constraint equation. Finally, we discuss the subsidiary system which governs the violation of the constraints. For a more detailed derivation see for example Friedrich (2002).

II.1 The general conformal field equations

Here we give a brief derivation of the GCFE, beginning with some mathematical preliminaries. The most important part of the GCFE is obtained from the Bianchi identities for the physical and the conformal metric together with Einstein’s field equation which we take to be vacuum here, but admitting a non-vanishing cosmological constant $\lambda$. From these we obtain equations for the rescaled Weyl tensor and the Schouten (or, equivalently, the Ricci tensor) for a conformal Weyl connection. We write these equations in terms of a tetrad and obtain the equations for the tetrad components as well as the connection coefficients from Cartan’s structure equations as usual.

We start with a vacuum space-time $(\widetilde{M},\tilde{g})$ where $\tilde{g}$ is a solution to the Einstein field equation with cosmological constant $\lambda$

[TABLE]

with $\tilde{R}_{ab}$ denoting the Ricci tensor of $\tilde{g}$.

Let $g_{ab}=\Theta^{2}\widetilde{g}_{ab}$ be a conformally related metric with conformal factor $\Theta$. We denote by $\nabla_{a}$ and $\widetilde{\nabla}_{a}$ the Levi-Civita connections of $g_{ab}$ and $\widetilde{g}_{ab}$ respectively. Furthermore, let $\widehat{\nabla}_{a}$ be a Weyl connection, i.e. a connection which is torsion free and compatible with the conformal class of $g_{ab}$ but not necessarily compatible with any metric in the conformal class. Then, there exist smooth 1-forms $f_{a}$ and $b_{a}$ such that

[TABLE]

It follows immediately that the 1-forms are related via

[TABLE]

For reasons that will become apparent later, we also define another 1-form as

[TABLE]

Any two covariant derivative operators differ by a $(2,1)$-tensor. For the difference between $\widehat{\nabla}_{a}$ and $\nabla_{a}$ we denote this tensor by $f_{ab}{}^{c}$. Acting on co-vectors yields

[TABLE]

from which it can be shown that

[TABLE]

where $f_{ab}{}^{c}$ is given by

[TABLE]

We now introduce arbitrary coordinates $x^{\mu}$ and frame field $e^{a}_{\mathbf{a}}$, which is orthonormal with respect to the conformal metric $g_{ab}$ and we define $\eta_{\mathbf{ab}}:=g_{ab}e^{a}_{\mathbf{a}}e^{b}_{\mathbf{b}}=\mathrm{diag}(+1,-1,-1,-1)$. We then have (with the obvious notation $\widehat{\nabla}_{\mathbf{a}}=e^{a}_{\mathbf{a}}\nabla_{a}$)

[TABLE]

By writing (2) in the basis $(e_{\mathbf{a}})_{\mathbf{a}=0:3}$ and contracting with $e^{b}_{\mathbf{b}}e^{c}_{\mathbf{c}}$, we find

[TABLE]

We are now in a position to write down the field equations. The first two are Cartan’s two structure equations Kobayashi and Nomizu (1963), where the torsion-free equation is given by

[TABLE]

and the curvature equation is given by

[TABLE]

Contracting these equations with $e^{a}_{\mathbf{a}}e^{b}_{\mathbf{b}}$ and using the decomposition of the Riemann tensor in terms of the Schouten and Weyl tensors (note that the Schouten tensor is not necessarily symmetric for a general Weyl connection)

[TABLE]

yields the two equations

[TABLE]

The next equation comes from the Bianchi identity for the vacuum metric $\widetilde{g}_{ab}$,

[TABLE]

Rewriting this in terms of the conformal connection and gravitational tensor $K_{abc}{}^{d}:=\Theta^{-1}C_{abc}{}^{d}$ gives the simple equation

[TABLE]

Written in terms of the Weyl connection, we find

[TABLE]

The final equation is obtained from the Bianchi identity for the Weyl connection

[TABLE]

which, after a bit of work using (9), gives

[TABLE]

Then the full set of conformal field equations expressed with a Weyl connection is given by (7),(8),(9) and (10), where in (8) $C_{\mathbf{a}\mathbf{b}\mathbf{c}}{}^{\mathbf{d}}$ is replaced by $\Theta K_{\mathbf{a}\mathbf{b}\mathbf{c}}{}^{\mathbf{d}}$.

These equations are invariant under several transformations. Since they are obtained from geometric differential equations they are invariant under arbitrary coordinate and tetrad transformations. However, they are also invariant under the simultaneous rescalings $\Theta\mapsto\Theta\phi$ and $g_{ab}\mapsto\phi^{2}g_{ab}$ with some arbitrary scalar field $\phi$ and under the change of Weyl connection $\widehat{\nabla}_{a}$ which amounts to the choice of a 1-form $b_{a}$. Thus, we will have to fix the freedom in these transformations in order to get determined systems. This means that we need to fix coordinates, the tetrad field, the conformal factor $\Theta$ and the 1-form $b_{a}$ or, with $\Theta$ fixed, of $h_{a}$.

A particularly useful gauge is the conformal Gauss gauge (CGG) which is obtained from the conformal structure itself and which we discuss next.

II.2 The conformal Gauss gauge

The conformal Gauss gauge is fixed solely from the conformal structure of the space-time, making use of the associated conformal geodesics. These are curves governed by the equations

[TABLE]

for the tangent vector to the curves $u^{a}$ and a smooth 1-form $h_{a}$ along them. Here, $\nabla_{a}$ is any covariant derivative compatible with the conformal class $[\widetilde{g}_{ab}]$ and $P_{ab}$ is its corresponding Schouten tensor. The metric $g_{ab}$ is any representative of the conformal class. Note, that the covariant derivative is not necessarily compatible with this metric. The equations are covariant under the change of Weyl connection in the following sense: let us suppose that we have a solution $\{u^{a},h_{a}\}$ of the conformal geodesic equations above. Let $f_{a}$ be any smooth 1-form then $\{u^{a},h_{a}-f_{a}\}$ is a solution of the conformal geodesic equation with $\nabla_{a}$ replaced by another Weyl connection $\widehat{\nabla}_{a}$ which is defined by

[TABLE]

We can now introduce the CGG according to the following procedure: we pick an initial space-like hyper-surface $\widetilde{\Sigma}_{0}$ in the physical vacuum space-time $(\widetilde{M},\widetilde{g})$. Then we can write the conformal geodesic equation in terms of the physical metric $\widetilde{g}_{ab}$ and its associated Levi-Civita derivative $\widetilde{\nabla}_{a}$. Note, that the Schouten tensor in these equations will reduce to a term proportional to $\lambda u_{a}$. We now select initial data for these equations on $\widetilde{\Sigma}_{0}$. We denote quantities which are defined only on $\widetilde{\Sigma}_{0}$ by an underline. On $\widetilde{\Sigma}_{0}$, we choose a conformal factor $\underline{\Theta}$, a 1-form $\underline{b}_{a}$, coordinates $(\underline{x}^{1},\underline{x}^{2},\underline{x}^{3})$, and a frame $\underline{e}_{\mathbf{a}}$ which is orthonormal with respect to the conformal metric $\underline{g}_{ab}=\underline{\Theta}^{2}\underline{\widetilde{g}}_{ab}$ and such that $\underline{e}_{\mathbf{0}}$ is normal to $\widetilde{\Sigma}_{0}$. Then we put $u^{a}:=e_{\mathbf{0}}^{a}$ and use $(\underline{u}^{a},\underline{b}_{a})$ as the initial data for the conformal geodesic equations. Their solutions $(u^{a},h_{a})$ provide a unique time-like conformal geodesic through each point of $\widetilde{\Sigma}_{0}$ at least near $\widetilde{\Sigma}_{0}$. Provided the initial data are sufficiently smooth, this congruence is smooth and caustic free in a local neighbourhood $U\subset\widetilde{M}$ of $\widetilde{\Sigma}_{0}$. Thus, at each point in $U$ we have a 1-form $h_{a}$. We denote this 1-form by $b_{a}$ and use it in place of $f_{a}$ to define a Weyl connection $\widehat{\nabla}_{a}$ as shown above. Thus, the 1-form $b_{a}$ mediates between $\widehat{\nabla}_{a}$ and the physical connection $\widetilde{\nabla}_{a}$. Written with respect to the new Weyl connection, the conformal geodesic equations simplify to

[TABLE]

We can now use this congruence of time-like curves to define on $U$ the parameter $s$ as the (time) coordinate $x^{0}$ and the coordinates $(x^{1},x^{2},x^{3})$ by the requirement that they be constant along the curves and agree with $(\underline{x}^{1},\underline{x}^{2},\underline{x}^{3})$ on $\widetilde{\Sigma}_{0}$. A tetrad field $e_{\mathbf{a}}$ is defined on $U$ by parallel transport of the chosen frame on $\widetilde{\Sigma}_{0}$ using the Weyl connection $\widehat{\nabla}_{a}$. Thus, each frame vector satisfies the equation

[TABLE]

We point out here that there is an additional freedom in the transport of the frame along the congruence. Instead of parallel transport we can rotate the spatial legs of the frame in an arbitrary way without affecting the geometry. This amounts to imposing the slightly more general equation

[TABLE]

where $\omega^{\mathbf{b}}{}_{\mathbf{c}}$ is an infinitesimal spatial rotation, i.e. it is characterized by $\omega_{\mathbf{bc}}=-\omega_{\mathbf{cb}}$ and $\omega_{\mathbf{0b}}=0$. It can be arbitrarily prescribed and we will use this freedom partially below.

The tetrad field defines a metric $g_{ab}$ which is necessarily conformal to the physical metric $\widetilde{g}_{ab}$ with a conformal factor $\Theta$ given by the equation

[TABLE]

Note, that these equations imply in particular that $\Theta^{-2}=\widetilde{g}(u,u)$ from which we find using (2) and (12a)

[TABLE]

where we used the over-dot to denote the derivative along $u^{a}$. This equation can be used to propagate the conformal factor. Using this and the conformal geodesic equations, we could also find a propagation equation for the 1-form $h_{a}$ that appears in the GCFE, thus fixing the gauge completely. However, a closer look at these equations reveals a surprising consequence: when looking at successive derivatives of the conformal factor one finds

[TABLE]

Thus, $\Theta$ is given as a quadratic polynomial in $s$. Its derivatives can be expressed in terms of the scalars

[TABLE]

and with these we obtain an explicit formula for the conformal factor, namely

[TABLE]

in terms of initial data defined on $\widetilde{\Sigma}_{0}$.

Furthermore, the conformal geodesic equation and the fact that the frame is transported along the curve using (14) yield an equation for the frame components of the 1-form $b_{a}$. With vanishing infinitesimal rotation $\omega^{\mathbf{b}}{}_{\mathbf{c}}$ and written with respect to the physical geometry it becomes

[TABLE]

Since $u^{a}=e^{a}_{\mathbf{0}}$ and with (15) we can obtain an equation for the frame components of the 1-form $h_{a}=\Theta b_{a}$ appearing in the GCFE

[TABLE]

Noting that $\dot{H}=2\dddot{\Theta}=0$, we find the remarkable result that

[TABLE]

When the more general frame transport equation (14) with non-vanishing $\omega$ is imposed, then the equation for the spatial components of $h_{a}$ changes into

[TABLE]

Then, the equation cannot be solved explicitly anymore, but must be added as an evolution equation to the system.

In summary then, the conformal Gauss gauge is characterized by

the coordinate conditions

[TABLE] 2. 2.

the frame conditions $u^{a}=e_{\mathbf{0}}^{a}$ and

[TABLE] 3. 3.

the explicit form (16) for the conformal factor and 4. 4.

the conditions for the conformal connection expressed either in the explicit form (17) or the transport equation (18) for $h_{\mathbf{a}}$ and the condition on the Schouten tensor (12b) which when expanded in the frame becomes

[TABLE]

Let us point out a particular property of this gauge. The fact that the frame is parallel along the curves implies that, in general, it will not be aligned with the hyper-surfaces defined by constant $s$. This implies that the expansion of the tetrad vectors with respect to the coordinate basis will include a component along $\partial_{s}$ and that the causal character of these hyper-surfaces may change during the evolution.

II.3 Imposing the $\eth$-calculus

The GCFE will allow us to evolve data from the initial surface up to and beyond null-infinity. This is a hyper-surface which is time-like, null or space-like depending on the cosmological constant, but it is always topologically of the form $\mathbb{R}\times S^{2}$. This property suggests to set up numerical methods which are adapted to this spherical topology. However, this is not straightforward since the sphere cannot be covered by a single coordinate chart. Numerically, this leads to the so called pole problem, the fact that the usual polar coordinates become singular on the poles. Another aspect of this problem appears when one uses a frame formalism, since there are no globally defined frames on the sphere.

These problems have been tackled by several researchers using interpolation between two charts Gómez et al. (1997), the cubed sphere Lehner et al. (2005), pseudo-spectral methods Bonazzola et al. (1999) and the eth-calculus Bartnik and Norton (2000). To our knowledge there is no Finite Element treatment of the sphere. Here, we elaborate on the $\eth$-calculus which we have implemented in our code. This has already been discussed extensively in previous papers Beyer et al. (2014, 2015) so we will keep the exposition here reasonably short.

The $\eth$-calculus can be regarded as a generalization of spectral methods on the circle $S^{1}$ based on Fourier series to the sphere $S^{2}$. It is based on harmonic analysis on $S^{3}$ and transferred to $S^{2}$ using the so called Hopf fibration $S^{3}\to S^{2}$. In the context of General Relativity it was first discussed by Newman and Penrose Newman and Penrose (1966), see also Goldberg et al. (1967); Dray (1985, 1986). There are two essential points: the first is the fact that components of tensors of different rank on $S^{2}$ can be expanded in terms of a complete system of “functions” on the sphere, called the spin-weighted spherical harmonics ${}_{s}\!\mathrm{Y}_{lm}$. And the second point is that the covariant derivative in the direction of a certain complex vector is diagonal when written with respect to the ${}_{s}\!\mathrm{Y}_{lm}$.

To briefly discuss the fundamentals of the $\eth$-formalism we consider the unit sphere $S^{2}$ with its metric $q$ which we take to be negative definite here in order to agree with our conventions for the space-time metric. At each point of $S^{2}$ we can introduce an orthonormal frame or, equivalently, a complex tangent vector $\mathbf{M}$ with the properties $q(\mathbf{M},\mathbf{M})=0$ and $q(\mathbf{M},\overline{\mathbf{M}})=-1$. Obviously, this vector is not uniquely defined but only up to a phase transformation

[TABLE]

which corresponds to a frame rotation with the angle $\alpha$.

It is well known that the sphere does not admit a globally defined frame field, hence no globally defined vector $\mathbf{M}$. To reach every point on $S^{2}$ one needs to patch several charts with frames defined on them together in order to cover the entire sphere. On a point in the overlap region between two charts the frames respectively, the complex vectors from the two patches differ by a rotation, respectively a phase transformation of the form (19). Components of tensor fields on $S^{2}$ transform under a well defined way under this change of frame on the overlaps. As an example consider a 1-form $v$. Its contraction $\left<v,\mathbf{M}\right>$ with $\mathbf{M}$ transforms under (19) according to $\left<v,\mathbf{M}\right>=\mathrm{e}^{\mathrm{i}\alpha}\left<v,\mathbf{M}\right>$. Components of higher rank tensors transform in similar characteristic ways under frame rotations. This leads to the notion of a spin-weighted “function” with spin-weight $s$ on $S^{2}$, i.e., a quantity $\eta$ which transforms as $\eta\mapsto\mathrm{e}^{\mathrm{i}\alpha s}\eta$ under (19). Such a quantity can be regarded as globally defined on $S^{2}$ since it essentially incorporates the collection of a component of a globally defined tensor field with respect to different charts 111Strictly speaking, spin-weighted functions are global sections of an appropriate line bundle over $S^{2}$ associated with the frame bundle..

A spin-$s$ quantity $\eta$ can be represented as a series

[TABLE]

where the ${}_{s}\!\mathrm{Y}_{lm}$ are the spin-weighted spherical harmonics, a generalization of the well-known spherical harmonics $Y_{lm}={}_{0}\!\mathrm{Y}_{lm}$ (see Penrose and Rindler (1986) for more details).

With spin-weighted quantities being globally defined their derivatives should share this property. The usual coordinate derivatives or directional derivatives along frame vectors do not achieve this. Instead one uses the connection $\overset{q}{\nabla}$ associated with the metric $q$ to define covariant derivatives along $\mathbf{M}$ to obtain well-defined spin-weighted quantities. This leads to the definition of the $\eth$-operator: suppose $\eta=T(\mathbf{M},\ldots,\overline{\mathbf{M}})$ is a spin-weighted quantity obtained from a tensor $T$ by evaluation on $k$ copies of $\mathbf{M}$ and $l$ copies of $\overline{\mathbf{M}}$ so that $\eta$ has spin-weight $s=k-l$. Then

[TABLE]

In similar way, the $\eth^{\prime}$-operator is defined

[TABLE]

Thus, $\eth$ raises the spin-weight by one, while $\eth^{\prime}$ lowers it by one. It follows from these expressions that the action of $\eth$ and $\eth^{\prime}$ on a spin-$s$ quantity $\eta$ can be alternatively expressed by

[TABLE]

where $\alpha$ is the single (complex) connection coefficient of $\overset{q}{\nabla}$ defined in terms of the commutator $[\mathbf{M},\overline{\mathbf{M}}]=\alpha\mathbf{M}-\bar{\alpha}\overline{\mathbf{M}}$.

The action of these operators on the spin-weighted spherical harmonics is particularly simple

[TABLE]

It is this property which makes these operators useful for numerical purposes: to compute their action on any spin-weighted quantity $\eta=\sum_{l,m}\eta_{lm}\,{}_{s}\!\mathrm{Y}_{lm}$ is just a matter of multiplying the coefficients by the appropriate numbers and replacing ${}_{s}\!\mathrm{Y}_{lm}$ with ${}_{s\pm 1}\!\mathrm{Y}_{lm}$.

To set up the framework for the $\eth$-calculus we need to make a few assumptions. As discussed in Sec. II.2, imposing the CGG fixes coordinates and tetrad in the space-time $M$ leaving us only the freedom to fix them on the initial hyper-surface $\Sigma_{0}$. We use this freedom to postulate that the initial hyper-surface $\Sigma_{0}$ is foliated by a family of 2-surfaces with spherical topology, so that it has topology $\mathbb{R}\times S^{2}$. We choose a radial coordinate $\rho$ on $\Sigma_{0}$ labelling the spherical leaves of this foliation. The remaining two coordinates label points on the spheres.

On each sphere $S_{\rho}\subset\Sigma_{0}$ we introduce a (fiducial) negative definite unit-sphere metric $q=-(\mathrm{d}\theta^{2}+\sin^{2}\theta\,\mathrm{d}\phi^{2})$. This also introduces coordinates $\theta$ and $\phi$ which can be interpreted as the usual polar coordinates on each $S_{\rho}$. We assume that $\theta$ and $\phi$ are smooth on $\Sigma_{0}$ (except for the usual problems at the poles of the spheres). On each sphere $S_{\rho}$ we can also introduce a complex null-vector $\mathbf{M}$ as before.

The coordinates $(\rho,\theta,\phi)$ are constant along the conformal geodesics so that every event in the region of $M$ covered by the conformal Gauss gauge is characterized by coordinates $(s,\rho,\theta,\phi)$ with $s$ the parameter along the conformal geodesics. Thus, each event lies on a unique sphere $S_{s,\rho}$ defined by constant $s$ and $\rho$, carrying a negative definite unit-sphere metric $q$ and a complex vector $\mathbf{M}$.

We also need to pick a triad $(e_{1},e_{2},e_{3})$ on $\Sigma_{0}$. We can arrange this to be adapted to the foliation in the sense that at every point $e_{2}$ and $e_{3}$ are tangent to the sphere $S_{\rho}$ through that point. Then $e_{1}$ is necessarily perpendicular to $S_{\rho}$. We can further arrange that $m:=\frac{1}{\sqrt{2}}(e_{2}-\mathrm{i}e_{3})$ is proportional to $\mathbf{M}$ on $\Sigma_{0}$. The ambiguity in the choice of $\mathbf{M}$ is the same as the ambiguity in the choice of an adapted triad, namely a phase transformation, respectively a frame rotation leaving $e_{1}$ fixed.

As described in section II.2, the conformal Gauss gauge provides a space-time frame $(e_{0},e_{1},e_{2},e_{3})$ on $M$ which it is orthonormal with respect to the conformal metric $g$. When the GCFE are referred to this frame they turn into equations for components with respect to the frame in terms of directional derivatives along the frame vector fields. The gauge conditions imply that the frame vectors can be written in terms of coordinate derivatives or, since $\partial_{\theta}$ and $\partial_{\phi}$ are intrinsic to the spheres $S_{s,\rho}$, in terms of $\partial_{s}$, $\partial_{\rho}$, $\mathbf{M}$, and $\overline{\mathbf{M}}$.

Since $\mathbf{M}$ is Lie dragged along the conformal geodesics while the frame (and hence $m$) is parallel with respect to the conformal connection, the vectors $\mathbf{M}$ and $m$ will in general differ at events past $\Sigma_{0}$. It is not even clear that $m$ will remain tangent to the spheres $S_{s,\rho}$. This discrepancy raises an issue: we take components with respect to the space-time frame $(e_{0},e_{1},e_{2},e_{3})$ but we use the derivative operators defined with respect to the fiducial $\mathbf{M}$ and $\overline{\mathbf{M}}$. In order to consistently use the $\eth$ and $\eth^{\prime}$ operator as defined above in terms of $\mathbf{M}$ and its complex conjugate we have to make sure that the components taken with the space-time frame nevertheless are properly spin-weighted when $\mathbf{M}$ is transformed according to (19).

This is indeed the case as the following argument shows. We first note that both frame transports have the property that a frame rotated by an angle $\alpha$ is transported into a frame rotated by the same angle. Thus, a phase transformation of $\mathbf{M}$ by an angle $\alpha$ at an event $P_{s}$ with coordinates $(s,\rho,\theta,\phi)$ corresponds to a phase transformation of $\mathbf{M}$ by $\alpha$ at the point $P_{0}$ with coordinates $(0,\rho,\theta,\phi)$ on $\Sigma_{0}$, and hence to the same phase transformation of $m$ at $P_{0}$, which via parallel transport corresponds to a phase transformation by $\alpha$ for $m$ at $P_{s}$. In this way, the phase transformations of $\mathbf{M}$ and $m$ are tied together at every event and, therefore, the components taken with respect to space-time frame are properly spin-weighted under changes of $\mathbf{M}$, and we may apply the $\eth$-formalism to them. The same argument applies to the “steered transport” described in section II.2 under certain circumstances that we will describe next.

In order to simplify the boundary treatment and the analysis on $\mathscr{I}$ we steer the spatial part of the tetrad using the free infinitesimal rotation in such a way that it is adapted at every point in $M$ to the sphere through that point: we want $e_{2}$ and $e_{3}$ to remain tangent to the sphere. This implies that $e_{2}(\rho)=0=e_{3}(\rho)$ must be satisfied throughout the evolution. Taking a time derivative of $e_{2}(\rho)$ along the vector $u$ and using the frame transport equation (14) yields

[TABLE]

and, similarly, for $e_{3}(\rho)$. If $e_{2}(\rho)=e_{3}(\rho)=0$ everywhere, then a part of the infinitesimal rotation is necessarily determined by

[TABLE]

On the other hand, if these conditions hold, then the equations derived above form a homogeneous linear system for $e_{2}(\rho)$ and $e_{3}(\rho)$ and hence these functions vanish everywhere if they vanish initially. Note, that these gauge conditions still do not fix the frame completely since $\omega^{2}{}_{3}$ remains undetermined, i.e., we are still free to arbitrarily rotate $e_{2}$ and $e_{3}$ in the tangent plane to the spheres. As a consequence the relationship between the complex vectors $m$ and $\mathbf{M}$ is fixed in exactly the same way as in the case without steering on the initial surface $\Sigma_{0}$ and the application of the $\eth$-formalism is consistent.

The conditions (25) have a small flow-on effect on the frame and tetrad equations in the GCFE system, but not on the tensorial equations for the Schouten and Weyl tensors.

At each point in $M$ we can expand the frame vectors $e_{c}$ in the coordinate basis $(\partial_{s},\partial_{\rho},\partial_{\theta},\partial_{\phi})$ or equivalently, in terms of $\mathbf{M}$ and $\overline{\mathbf{M}}$ in the form

[TABLE]

When we assume the frame to be transported using steering, then the coefficients $c_{2}^{1}$ and $c_{3}^{1}$ vanish identically.

The preparation of the GCFE as given in (7–10) now proceeds in a straightforward way. Expand the equation with respect to the frame $(e_{0},e_{1},e_{2},e_{3})$, replace the directional derivatives along the frame vectors by the coordinate derivatives $\partial_{s}$, $\partial_{\rho}$ and the derivatives along $\mathbf{M}$ and finally replace the latter by the appropriate $\eth$ operators. This procedure yields a system of equations which is properly spin-weighted.

In our code we represent the spin-$s$ quantities $\eta$ as grid functions on the sphere sampled at points $(\theta_{i},\phi_{k})$ and we use

[TABLE]

Then the connection coefficient $\alpha$ in (23) becomes

[TABLE]

The code alternates between representing $\eta$ as a grid function or a truncated series (20) with $s$ and $\rho$ dependent coefficients $\eta_{lm}$. The series representation is used to compute the derivatives along the sphere using (24), while the grid representation is used to perform algebraic operations such as multiplication etc. These issues are explained in much more detail in Beyer et al. (2015, 2014, 2016).

II.4 The GCFE in the space-spinor formalism

In this section we present the full GCFE in the space-spinor formalism Sommers (1980). There are two reasons for this: first, the space-spinor formalism allows us to easily perform the time-space split of the equations into constraints and evolution equations by simply computing irreducible parts of the resulting spinor equations. Second, the formalism deals with complex quantities which means that the number of the actual equations after taking components is reduced compared to the tensorial equations by almost half. This leads to less numerical complexity and therefore ultimately also to lower round-off error (even though we have not checked that). Furthermore, it is easier to deal with symmetric trace-free quantities such as the electric and magnetic parts of the rescaled Weyl tensor because the trace-free condition is automatically satisfied in a spinor version compared to the tensorial version. Even though this means more algebraic work we feel that this is worthwhile. The space-spinor formalism is briefly summarized in the appendix A.

The idea is to express the covariant derivative operator $\widehat{\nabla}_{a}$ in terms of $\nabla_{a}$, the Levi-Civita operator of the conformal metric $g_{ab}$ and the mediating 1-form $f_{a}$ (see (5)), then split them into evolution and constraint equations by decomposing the resulting equations. Once this is done the CGG is imposed and we obtain a complete system of evolution and constraint equations. We will not present the full derivations here (see Stevens (2016) for a full account), instead we give a brief summary.

The standard procedure is:

• •

Replace $\widehat{\nabla}$’s with $\nabla$’s using the transformation laws (5) and (6),

• •

convert the tensor fields to space-spinor fields

• •

convert $\nabla_{a}$ into $\partial$ and $\partial_{AB}$ as described in app. A,

• •

decompose the equations into irreducible pieces to obtain evolution and constraint equations.

We apply this procedure to the field equations (8) and (9) for $\widehat{P}_{ab}$ and $K_{abc}{}^{d}$. In doing so we will also change our notation slightly: the Schouten tensor $\widehat{P}_{ab}$ related to the Weyl connection is the only Schouten tensor appearing in the equations, so in order to avoid confusion with complex conjugation of space-spinors we will drop the hat. Thus, from now on $P_{ab}$ will denote the Schouten tensor for the Weyl connection.

To illustrate the procedure, we show how to treat the equation for $K_{abc}{}^{d}$. We first define $\psi_{ABCD}=\Theta^{-1}\Psi_{ABCD}$ as the only irreducible piece of $K_{abcd}$ in the spinor formalism. Skipping the first step as (9) is already in terms of the conformal connection, we convert tensor to spinor indices and then convert primed to unprimed indices. Splitting the covariant derivative into $D$ and $D_{AB}$ derivatives and expressing these in terms of $\partial$ and $\partial_{AB}$ yields the equation

[TABLE]

This equation is symmetric in $ABC$. Taking the totally symmetric part gives an evolution equation for $\psi_{ABCD}$

[TABLE]

while the anti-symmetric part gives us a constraint

[TABLE]

Up to this point the choice of $t^{a}$ was irrelevant. However, to derive the gauge-related equations from the remaining equations (7) and (8), we must fix the time-like vector $t^{a}$ and derive the spinor form of the CGG equations. We fix $t^{a}$ to be proportional to the tangent vector $u^{a}$ of the curves. Then the spin-frame $(o_{A},\iota_{A})$ corresponding to the tetrad must satisfy the usual normalization condition $o_{A}\iota^{A}=1$ and we impose the condition $t_{AA^{\prime}}=o_{A}o_{A^{\prime}}+\iota_{A}\iota_{A^{\prime}}$. Then we have

[TABLE]

and the transport equations (13) or (14) for the frame imply equations for the spin frame. We will explicitly discuss here the equations for the pure CGG without frame rotations, the other case being very similar. Thus, the spin frame must satisfy the equation $t^{a}\widehat{\nabla}_{a}o_{B}=0$ (and its complex conjugate). We define the components of the frame vectors in terms of the coordinate basis by

[TABLE]

Since both $t^{a}$ and the time-like frame vector are parallel to $u^{a}$ we have

[TABLE]

The gauge conditions then imply

[TABLE]

Furthermore, we define the spinor fields

[TABLE]

They satisfy the relations

[TABLE]

These spinor fields encode part of the connection coefficients $\widehat{\Gamma}_{\mathbf{ab}}{}^{\mathbf{c}}$, the other pieces being contained in $K_{AB}$, $K_{ABCD}$ and $f_{AB}$.

The tetrad is parallel along the conformal geodesics so that we have $u^{a}\widehat{\nabla}_{a}g_{bc}=0$ from which we conclude together with eq. (2) that

[TABLE]

Furthermore,

[TABLE]

This and the analogous equation for $\iota_{C}$ yield the equations

[TABLE]

which imply

[TABLE]

The last consequence of the gauge conditions involves the Schouten spinor of the Weyl connection, which we denote in the space-spinor formalism as $P_{ABCD}:=t_{B}{}^{A^{\prime}}t_{D}{}^{B^{\prime}}P_{AA^{\prime}CB^{\prime}}$ (again note, that we dropped the hat here to avoid confusion with complex conjugation as mentioned above). This spinor field satisfies the reality condition

[TABLE]

and the gauge condition (12b) implies

[TABLE]

This concludes the translation of the CGG to the space spinor formalism. We can now proceed to translate the tensorial equations to spinorial form following the procedure outlined above. Incorporating the gauge conditions we arrive at a fully determined system for the unknowns

[TABLE]

in which $\Theta$ and the components of $h_{AB}$ are considered known since they are determined a priori from initial data by the equations (16) and (17).

The evolution equations are

[TABLE]

The constraint equations express the vanishing of the following “zero quantities”

[TABLE]

In order to impose the $\eth$-calculus on our system, we replace the expansion $\partial_{AB}=c^{\mu}_{AB}\partial_{\mu}$ with

[TABLE]

where $X$ and $Y$ are complex functions while $R$ can be chosen to be real as a consequence of the remaining gauge freedom on the spheres. This is essentially a repackaging of the quantities $c^{2}_{AB}$ and $c^{3}_{AB}$. It is now a straightforward but tedious procedure to replace the derivatives along $\mathbf{M}$ and its complex conjugate in terms of $\eth$ and $\eth^{\prime}$ using (23), which ultimately yields properly spin-weighted equations for the spinor components of the fields. The equations for the new functions $R$, $X$ and $Y$ come as evolution equations

[TABLE]

and constraint equations

[TABLE]

where the “algebraic terms” are quite lengthy containing components of $K_{ABCD}$, the spin-weighted components of $\gamma_{ABC}$ and the frame components.

The new proper spin-weighted system then comprises of the evolution equations, constraint equations and explicit expressions for $\Theta$ and $h_{a}$ given in (28), (29) (16) and (17) but now the equations for $c^{2}_{AB}$ and $c^{3}_{AB}$ are replaced by equations for $R,\;X$ and $Y$ given in (30) and (31).

As a partial justification for the correctness of the final system, we expressed the Minkowski and Anti-de Sitter space-times analytically in the CGG adapted to spherical symmetry. The unknowns of the system were then computed from these solutions and were shown to satisfy both the evolution and constraint equations identically.

II.5 The subsidiary system

The constraint equations derived from the GCFE propagate on an analytical level, as shown by Friedrich Friedrich (1995), in the sense that the “zero-quantities” themselves satisfy a semi-linear homogeneous symmetric hyperbolic system of PDE. For the initial value problem this implies that if these quantities vanish initially then they vanish everywhere due to the uniqueness of solutions. Constraint propagation at the numerical level is another story entirely. Clearly, the violation of the constraints, i.e., non-vanishing “zero-quantities” are governed by this propagation system. Therefore, it is important to understand the properties of this system since it has implications on the stability and hence on the well-posedness of the evolution.

We have derived the subsidiary system for all the constraints in our formulation. This is a straightforward but awkward process that we will not elaborate here. Instead, we list only the propagation equation for the constraints arising from the system for the rescaled Weyl spinor since this is the only one needed in the sequel.

[TABLE]

Here, $\widetilde{K}_{ABCD}$, $\widetilde{K}_{AB}$ and $\widetilde{K}$ are the irreducible parts of $K_{ABCD}$. The complete subsidiary system can be found in appendix B.

III The IBVP framework

We now suppose, and will assume for the rest of the document, that the CGG has the coordinate system $\{s,\rho,\theta,\phi\}$. This is to make the distinction between the coordinates used in the CGG gauge and other frequently used coordinates such as $t$ and $r$, which we reserve for exact solutions such as the Schwarzschild space-time.

III.1 Maximally dissipative boundary conditions

We choose the spatial extent of the computational domain to be defined by $\rho\in[\rho_{0},\rho_{1}]=:I$ resulting in the cylindrical domain $I\times S^{2}$. Due to our choice of coordinates, the boundary of the domain is given by the equations $\rho=\rho_{0}$ on the left boundary and $\rho=\rho_{1}$ on the right boundary, defining two hyper-surfaces which (at least for some time) are time-like. Inspection of the evolution equations (28) and (30) shows that the only equation which needs boundary conditions is the evolution equation (28g) for $\psi_{ABCD}$. All other equations propagate along the boundary.

It is well known that the field equation for the gravitational spinor admits a symmetric hyperbolic system of PDEs (see for example Friedrich (2004), Friedrich (1996)). The equation (28g) is already in symmetric hyperbolic form. This is best verified by checking that the symbol of the equation is Hermitian. For any pair $(p,p_{AB})$ of a positive real number $p$ and real symmetric spinor $p_{AB}$ consider the sesquilinear form on totally symmetric spinors of rank 4 defined by

[TABLE]

It is obtained by replacing the derivative operators in the principal part of (28g) by $p$ and $p_{AB}$, respectively. Now it is easy to see that this form is Hermitian, i.e., that

[TABLE]

for any choice of $(p,p_{AB})$ and positive definite for $p_{AB}=0$. Therefore, the system (28g) is symmetric hyperbolic.

To impose boundary conditions it is useful to first analyse the characteristics of the system near the boundaries. We choose the right boundary here, the left being treated analogously. We are interested in the characteristics that cross the time-like hyper-surface defined by $\rho=\rho_{2}$. Since the system is symmetric hyperbolic, at each point $P$ of the boundary there is a family of characteristic cones with vertex $P$ opening towards the past along which the field propagates with characteristic speeds. We want to find the intersection of these cones with the 2-dimensional plane spanned by $\partial_{s}$ and $\partial_{\rho}$. These will be lines with slopes given by the characteristic speeds $\lambda$. They are obtained from a generalized eigen-value problem which arises from the principal part of (28g) by inserting the ansatz $\psi_{ABCD}=X_{ABCD}\,\mathrm{e}^{\mathrm{i}(\rho-\lambda s)}$. Then the characteristic speeds $\lambda$ are obtained as those values for which the algebraic system

[TABLE]

has non-trivial solutions. Define the spinor $T_{AB}=\frac{1}{2}T\varepsilon_{AB}+S_{AB}$ with $T=\sqrt{2}\lambda$ and $S_{AB}=c^{1}_{AB}-\lambda c^{0}_{AB}$. Near the initial hyper-surface we can write $S_{AB}=s\alpha_{(A}\beta_{B)}$ in terms of its principal null spinors with $\alpha_{A}\beta^{A}=1$. Then we find that the system has non-trivial solutions only if

[TABLE]

and we find the resulting speeds

[TABLE]

Here, we have defined $f^{ij}=c_{AB}^{i}c^{jAB}$ where $f^{ii}<0$. The corresponding solutions $X_{ABCD}^{k}$ are simply the totally symmetric outer products of the principal spinors $\alpha_{A}$ and $\beta_{A}$ with $k$ giving the number of $\beta_{A}$ in the product. Thus, $X_{ABCD}^{1}=\alpha_{(A}\alpha_{B}\alpha_{C}\beta_{D)}$ etc. It is not difficult to see that the lines defined by $\lambda_{0}$ and $\lambda_{4}$ correspond to the intersection of the light-cone with the plane spanned by $\partial_{s}$ and $\partial_{\rho}$, while $\lambda_{1}$ and $\lambda_{3}$ correspond to a time-like cone inside the light-cone. Thus, the components $X^{0}_{ABCD}$ and $X^{4}_{ABCD}$ propagate with the speed of light.

We decompose the spinor $\psi_{ABCD}$ with respect to the spin-frame $(\alpha_{A},\beta_{A})$, writing $\psi_{ABCD}=\sum_{k=0}^{4}\widetilde{\psi}_{k}X^{k}_{ABCD}$. Since $\lambda_{k}<0$ for $k=0,1$, the corresponding components $\widetilde{\psi}_{k}$ propagate towards decreasing $\rho$ while $\widetilde{\psi}_{3}$ and $\widetilde{\psi}_{4}$ travel towards increasing $\rho$. The component $\widetilde{\psi}_{2}$ propagates along the boundary. Thus, on the right boundary we need to provide boundary conditions for $\widetilde{\psi}_{0}$ and $\widetilde{\psi}_{4}$.

A very general class of boundary conditions are the maximally dissipative boundary conditions, which express the incoming components in terms of the outgoing ones and freely specifiable functions on the boundary. In our case on the right boundary we write

[TABLE]

where the $q_{i}$ are free boundary data and $H$ is a $2\times 2$ matrix defined on the boundary satisfying

[TABLE]

where $I_{2}$ is the identity matrix. The matrix $H$ encodes the reflective properties of the boundary and the inequality ensures that the reflected ‘energy’ does not exceed the ‘energy’ hitting the boundary. We will use the simplest case $H=0$, which corresponds to the boundary being completely transparent to the outgoing modes.

In order to implement these conditions we need to run through the following procedure at every point on the boundary: determine the principal null directions of the spinor $S_{AB}$, compute the characteristic speeds, project out the outgoing characteristic modes from $\psi_{ABCD}$, use the boundary condition to compute the ingoing characteristic modes, and then reassemble the spinor $\psi_{ABCD}$ on the boundary. This yields a system of the form

[TABLE]

Here, $M$ is a $2\times 5$-matrix made up from the components of the spin-frame transformation between $(o_{A},\iota_{A})$ and $(\alpha_{A},\beta_{A})$. Then, imposing the boundary conditions means solving this system for $\psi_{0}$ and $\psi_{1}$, regarding the other components as given with free functions $q_{k}$.

Note, that the expressions for $\lambda_{0}$ and $\lambda_{4}$ become singular when $1+f^{00}=0$. Since $g(\mathrm{d}s,\mathrm{d}s)=g^{ss}=1+f^{00}$ this will happen, when the hyper-surface $s=\mathrm{const}$ becomes null at the boundary. Due to our choice of gauge we cannot prevent this from happening. However when the frame components $c^{0}_{AB}$ vanish for all time, as discussed in section IV.1, we are in a similar situation to the Friedrich-Nagy gauge Friedrich and Nagy (1999) where the spatial frame vector normal to the boundary remains normal for all time. We will have more to say about these issues later.

The characteristics for this particular choice of evolution system are simple, there will only ever be two ingoing modes to provide boundary conditions for. However in general, this will not be the case. Adding combinations of the constraints (29) to the evolution equations will alter their characteristics, and hence the number of ingoing modes. The characteristic speeds are functions of the frame components, and hence the number of ingoing modes may change over time.

III.2 Constraint preserving boundary conditions

In the previous section we employed maximally dissipative boundary conditions to obtain a stable evolution system, however this did not take into account the propagation of the constraint system. As the constraint equations coming from the field equation for $\psi_{ABCD}$ are PDEs, if they are violated on the boundary these violations may propagate into the interior. As one can see in the subsidiary system given in section II.5, many of the constraint propagation equations in the subsidiary system are written in terms of the gravitational spinor’s constraints. Thus if they are violated, they will cause other constraints to be violated too.

At each boundary, the maximal dissipative boundary formalism yields two degrees of freedom for the boundary data for our evolution system. We, however, expect only one physical degree of freedom. This is most likely the ingoing mode travelling on the light-cone, i.e., $\widetilde{\psi}_{0}$ or $\widetilde{\psi}_{4}$ depending on the boundary. The discrepancy in the number of degrees of freedom is related to the violation of the constraints and we will use the additional ingoing mode to kill ingoing constraint violating modes. In order to define the required constraint preserving boundary conditions for our evolution system, we need to analyse the subsidiary system.

The principal part of the subsidiary equation for the gravitational spinor constraint $G_{AB}$ (44c) is

[TABLE]

This system is again symmetric hyperbolic and we can go through the same analysis as we did on the $\psi_{ABCD}$ evolution system in the previous section. We find that there are three characteristic modes propagating with three different characteristic speeds. The remarkable fact is that these three speeds agree with the characteristic speeds $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ of the $\psi_{ABCD}$ system, see (33). Furthermore, the characteristic modes are again given in terms of the principal null spinors $\alpha_{A}$ and $\beta_{A}$, being proportional to $\alpha_{A}\alpha_{B}$, $\alpha_{(A}\beta_{B)}$, and $\beta_{A}\beta_{B}$, respectively. Thus, the constraints propagate on characteristics which agree with the time-like characteristics of the Weyl system. In order to kill the ingoing constraint violating mode on the right boundary we need to ensure that the ingoing component proportional to $\alpha_{A}\alpha_{B}$ vanishes on the boundary.

The problem is then to find a way to prescribe boundary data for the $\psi_{ABCD}$ system on the boundary which gives interesting behaviour as well as having no ingoing mode in the subsidiary system. Since the physically relevant quantities at the boundary are $\psi_{0}$ resp. $\psi_{4}$ one may ask whether it is possible to leave their prescription free, while fixing the remaining freedom in choosing $\psi_{1}$ resp. $\psi_{3}$ in such a way that the ingoing mode of the subsidiary system is zero. This is in fact possible, and is the main result of the paper.

Our proposed resolution on the right boundary (the left is analogous) for the case of $\psi_{0},\psi_{1}$ ingoing, $\psi_{3},\psi_{4}$ outgoing and $\psi_{2}$ propagating along the boundary is to derive an ODE for the free datum $q_{1}$ as follows:

• •

Regard $\psi_{0}$ as a free function, i.e., the boundary data $q_{0}$ is a free spin-2 function.

• •

Take a time derivative of the boundary equations (36) for $\psi_{0}$ and $\psi_{1}$ and replace all the time derivatives with their evolution equations. This will result in two equations containing the terms $\partial_{s}q_{1},\partial_{\rho}\psi_{0}$ and $\partial_{\rho}\psi_{1}$.

• •

Solve these equations simultaneously for $\partial_{\rho}\psi_{0}$ and $\partial_{\rho}\psi_{1}$.

• •

Now consider the single equation which requires the vanishing of the ingoing constraint violating mode. This equation also involves $\partial_{\rho}\psi_{0}$ and $\partial_{\rho}\psi_{1}$. Replace them with the expressions from the previous step.

• •

Finally, solve this equation for $\partial_{s}q_{1}$.

The resulting equation is free of $\partial_{\rho}\psi_{0}$ and $\partial_{\rho}\psi_{1}$ and gives us an ODE for $q_{1}$ on the boundary. Thus, evolving this equation along the boundary yields $q_{2}$ which is used as the additional “unphysical” degree of freedom for the evolution that should kill the ingoing mode from the subsidiary system.

The final remark here is how to extend this approach to incorporate the changing propagation directions of the $\psi_{ABCD}$ components if another evolution system was chosen with different characteristics. We note that there is a relationship between the characteristics of the evolution system and the subsidiary system, see for example Gundlach and Martin-Garcia (2004). This relationship is that there are always $n-1$ ingoing modes for the subsidiary system and $n$ ingoing modes for the evolution system. Hence there is always a degree of freedom to choose ingoing gravitational radiation for example, while just the right number of ingoing modes remaining to kill those of the subsidiary system.

At this point it might be useful to take count of the degrees of freedom in the system. As in every covariant geometric system of equations we have the usual gauge freedom of the choice of coordinates and basis of the tangent vectors. Since we formulate our equations in a conformally covariant way we have further gauge freedom, namely the choice of a conformal factor and of the conformal Weyl connection. All these choices have to be fixed in order to get a well-posed evolution system. In our approach they are all satisfied by imposing the Gauss gauge and choosing appropriate initial data for coordinates, frame, conformal factor and Weyl 1-form. Given our choice of evolution system, this reduces the freedom to two complex functions on each boundary, describing the incoming modes. However, as discussed above, one of these functions is determined by the requirement that the ingoing constraint violating mode should vanish. This leaves us with exactly one free function on each boundary, the ingoing spin-2 component of the gravitational field. This is a complex-valued function of spin-weight $\pm 2$, so it encodes exactly the two polarization degrees of freedom in a gravitational wave. Note that when a different evolution system is chosen, the relationship between the evolution and subsidiary systems described in the previous paragraph ensures that there is still only ever one complex valued degree of freedom on each boundary.

IV Numerics

The aim of this section is to formulate a numerical implementation of the IBVP for the GCFE system and apply it in a variety of situations. The idea is to choose an initial space-like hyper-surface and prescribe on it data for the unknowns (27) so that the constraints (29) are satisfied. Then we evolve the initial data using the evolution equations (28) using boundary conditions (where applicable) according to process proposed in the previous section so the constraints remain satisfied during the evolution.

We must first discretize the continuous version of the equations in order to evolve the system on the computer using various numerical methods. We will discretize our system using the method of lines. In full generality, we need to discretize three spatial dimensions, two of which are tangent to unit 2-spheres 222This is because we introduced the unit-sphere derivative operators $\eth$, $\eth^{\prime}$. while the third dimension is transverse to them. The action of the spherical operators $\eth$, $\eth^{\prime}$ can be approximated using pseudo-spectral methods using the algorithm given in Beyer et al. (2016). The radial direction is approximated using a straightforward fourth-order finite difference method with the summation-by-parts property as described in Strand (1994). Discretizing in this way yields a semi-discrete system of ODEs in time, which we solve using a standard fourth-order explicit Runge-Kutta scheme. In order to impose stable boundary conditions, we utilize the simultaneous-approximation-term (SAT) method that is described in Carpenter et al. (1999).

Before carrying out a numerical investigation into the validity of this scheme, we first tested the system and numerical methods by solving several IVP’s for the Kottler space-time. This was done using different values and signs of the cosmological constant and by using the setup detailed in section IV.2. As these space-times are spherically symmetric, we found that our evolution and subsidiary systems reduced to ODEs in time, hence boundary conditions were not required. We discovered that we could reproduce the conformal structure of Schwarzschild, Schwarzschild-de Sitter (for the case of distinct cosmological and event horizons) and Schwarzschild-Anti-de Sitter space-times and the constraints were shown to propagate.

We also conducted noise tests for the evolution system to test for numerical stability by checking that high frequency perturbations did not trigger non-linear growth. To this end we followed the apples with apples robustness test Alcubierre et al. (2003) by using randomly perturbed Minkowskian initial and boundary data. We ran multiple simulations in axi-symmetry using $\rho$ and $\theta$ resolutions of $200$ and $32$ respectively with a timestep of $0.0025$ for $2000$ iterations. We calculated a “global” $l_{2}$-norm on each timeslice by summing the squares of each variable on each point, dividing by the number of points and then taking the square root. As can be seen in FIG. 2 this did not grow substantially over time, and for all practical purposes is negligible. This gives a good indication of numerical stability.

IV.1 The IBVP for non-linear gravitational perturbations of Minkowski space-time

The next test concerns the boundary implementation. In particular, we need to check to what extent the method we discussed in section III.2 actually kills ingoing constraint violating modes. We choose Minkowski initial data and imposing appropriate boundary conditions to shoot in gravitational waves, represented by $\psi_{0}$ and $\psi_{4}$ (propagating in the $\rho$-direction from right to left and left to right respectively). The boundary conditions are chosen to be axi-symmetric, i.e. independent of $\phi$. The conformal Gauss gauge preserves the $\phi$-independence throughout the evolution. Hence we have a $2+1$ problem.

We will consider the physical representation of a part of Minkowski space-time with metric

[TABLE]

and $\rho\in[\rho_{1},\rho_{2}]$ for some values $\rho_{1}<\rho_{2}$ and then perturb it with incoming gravitational waves. To choose the gauge conditions we take the conformal factor and the 1-form $h_{a}$ from the exact Minkowski space-time and impose them also in the perturbed space-time. Since we look at the physical representation of the Minkowski metric we have $\Theta=1$ and $h_{a}=0$. These are also the initial conditions for $\Theta$ and $h_{a}$. This choice says that the gauge is adapted to metric time-like geodesics and it has the consequence the frame components $c^{0}_{AB}$ vanish for all time, which can be seen from the evolution equations. Therefore, there is no $\partial_{s}$ contribution from the expansion of $\partial_{AB}$ during the evolution, so that the spatial frame vectors will stay tangent to the $s=\text{constant}$ hyper-surfaces.

The initial metric is

[TABLE]

and with the choice of initial triad as indicated in section II.3 we find the only non-vanishing system components to be

[TABLE]

The characteristic speeds of the components of $\psi_{ABCD}$ used in the evolution are given by

[TABLE]

for $\psi_{0},\cdots,\psi_{4}$ respectively.

We discretize the spatial directions by choosing equi-distant points in the 2-dimensional interval $[0.25,1.25]\times[0,\pi]$. As boundary conditions, we shoot in a gravitational quadrupole $(l=2)$ wave from each boundary by choosing the free data $q_{0}$ for $\psi_{0}$ on the right boundary and $q_{4}$ for $\psi_{4}$ on the left boundary as

[TABLE]

with the spin-weighted spherical harmonic ${}_{2}Y_{20}(\theta)\propto\sin^{2}{\theta}$ when written in the usual polar coordinates. We adopt our boundary treatment outlined in section III which fixes the boundary conditions for $\psi_{1}$ on the right boundary and $\psi_{3}$ on the left boundary. Hence we have fixed all the free data on the boundaries. We evolve up until $s=1$ with a $\theta$-resolution of $32$, $\rho$-resolutions of $\{25,50,100,200,400\}$ and use time-steps of $0.5\,\Delta\rho$. We find that this $\theta$-resolution is enough to represent all the functions in our system on the spheres $s=\text{constant}$, $\rho=\text{constant}$ at machine precision.

First we look at the case without employing the subsidiary-mode-killing boundary treatment, to confirm that constraint violating modes are indeed propagated in from the boundary. We use the maximally dissipative boundary conditions but choose the free data (i.e. the “$q$”) for $\psi_{1}$ on the right and $\psi_{3}$ on the left to be zero. Looking at the system variables on the slices $s=\text{constant},\;\theta=\frac{\pi}{2}$ in FIG. 3 we see that they converge everywhere with increasing resolution at the correct order of four. This means that the IBVP for the evolution system alone is well-posed. However, the constraint does not converge to zero everywhere. The constraint violating modes should propagate in a time-like manner i.e., with less than the speed of light and this is clearly the case: with each snapshot the region where the constraint converges to zero becomes smaller while near the boundaries the constraint mode seems to converge to a non-zero limit. This indicates that the constraint is not satisfied.

Next, we test the boundary treatment which supposedly kills the incoming constraint violating mode. In FIG. 4 we display the same constraint as in FIG. 3 except with our constraint preserving boundary conditions. One can see that there is no longer a constraint violation propagating inward from the boundary and the constraints converge across the entire grid at the correct order.

Analogous plots are seen in all the other constraints and also for different choices of $\theta=$constant. Although this was a very simple case, the premise of our boundary treatment method has been verified.

IV.2 The IBVP for non-linear gravitational perturbations of Schwarzschild space-time

In the previous section we used simple initial data and fixed the gauge freedom appropriately so that we could test our framework for the simple case of $c^{0}_{AB}(s,\rho,\theta)=0$, i.e. the spatial frame vectors remain tangential to the $s=\text{constant}$ hyper-surfaces, in analogy to the Friedrich-Nagy gauge. However, this is a very special case and we want to explore the more general case of the non-vanishing $c^{0}_{AB}$ as well. We choose Schwarzschild space-time in isotropic coordinates as the test case and solve the associated IBVP with axi-symmetric boundary conditions as in the previous section. The IBVP setup is shown in FIG. 5.

It has been shown by Friedrich Friedrich (2003) that there exists a specific choice of initial data for the CGG that globally covers the Schwarzschild-Kruskal space-time smoothly and without degeneracy. This involves writing the Schwarzschild metric in isotropic coordinates and choosing appropriate initial data for the conformal geodesics. In this section we present how we (following Friedrich (2003)) set up the GCFE system to investigate null and time-like infinity of Schwarzschild space-time. We then numerically evolve the resulting initial data and see how the system behaves near the conformal boundary and the singularity.

The Schwarzschild metric written in isotropic coordinates is

[TABLE]

We now take the hyper-surface $t=0$ as our initial hyper-surface $\Sigma_{0}$ and compute the corresponding initial data induced on it by the Schwarzschild metric. We also need to choose initial data $\underline{H}$, $\underline{Z}$ and $\underline{\Theta}$ so that we can compute the 1-form $h_{a}$ and conformal factor $\Theta$. Friedrich makes the choice

[TABLE]

which does not compactify the initial surface in the $\rho$ direction. Next, the 1-form is fixed initially by setting $\underline{f}_{a}=0$ so that $\underline{h}_{a}=\nabla_{a}\underline{\Theta}$. This also determines $\underline{H}$ and we choose $\underline{Z}=0$ to fix the remaining freedom. Using (16) and (17), this gives the explicit expression for $\Theta$ as

[TABLE]

and $h_{a}$, whose spinor representation we write in decomposed form as $h_{AB}+\frac{1}{2}\varepsilon_{AB}h$

[TABLE]

The final step is to prescribe the extrinsic curvature $K_{ABCD}$ which vanishes initially since the initial hyper-surface is time-symmetric.

From the initial metric chosen as (38) evaluated at $t=0$, expressions for $\Theta$ (39), $h_{AB}$ (40) and $\underline{K}_{ABCD}=0$ we find the remaining non-vanishing initial data to be

[TABLE]

A convenient function to know during the evolution is the original Schwarzschild radius $r(s,\rho)$. In spherical symmetry the radius of the 2-spheres in the GCFE system is $R$ and hence we can relate it the Schwarzschild radius by

[TABLE]

This is very useful as it will tell us where the event horizon is located and if or when we end up at the curvature singularity. This is exactly true in spherical symmetry, i.e., for the unperturbed space-time but it will also be approximately correct in the perturbed case.

We choose $m=0.5$ and choose the spatial extent of the computational domain as the 2-dimensional interval $[0.25,1.25]\times[0,\pi]$ as in the previous section. Finally, we need to specify boundary conditions. The left boundary starts at the cross-over surface at $\rho=m/2=\frac{1}{4}$. The Kruskal extension has a reflection isometry $\rho\mapsto m^{2}/4\rho$ which fixes the hyper-surface $\rho=m/2$. We impose this reflection symmetry as boundary condition on the left boundary. This implies the following conditions for the ingoing components of the rescaled Weyl spinor

[TABLE]

There will be no constraint violation propagating in from the inner boundary as we have not violated the constraints there, it will remain as Schwarzschild space-time until the gravitational wave coming in from the right boundary reaches it, which does not happen in our simulations.

On the right boundary we implement our constraint preserving boundary treatment. This is done by first choosing the free wave profile $q_{0}$ for $\psi_{0}$ to be

[TABLE]

where $a$ is a fixed constant representing the amplitude of the wave. The other ingoing mode is $\psi_{1}$. Its boundary value is determined to kill the ingoing constraint violating mode via our resolution.

Now we have all that is needed to start evolving the system. We note that we use a $\theta$-resolution of $64$ so that even at late times of the simulation our fields are well represented in the spin-weighted spherical harmonic basis. We use $\rho$-resolutions of $\{25,50,100,200,400\}$ which gives us ample data to check the propagation of the constraints. We evolve up to $s=1.22$ which incorporates $\mathscr{I}^{+}$ into the simulation, which is reached on the right boundary at around $s=0.83$. Simulations are also performed without our boundary treatment where the appropriate $q_{i},\;i\neq 0$ are set to zero. This will allow us to contrast the before and after of our boundary treatment and emphasize the problem that we resolve.

It is worth noting that we could change the boundary treatment when the right boundary has passed beyond $\mathscr{I}^{+}$. As the constraint violating modes are time-like, they can never propagate through $\mathscr{I}^{+}$ from outside and hence there is no need to kill them any longer. However, while this is true analytically we must be cautious with this kind of argument since the numerical propagation of constraint violation can be faster. This is a purely numerical artefact which will shrink as the $\rho$-resolution is increased.

We first present a convergence plot for the case of the simple, but non-constraint-preserving, choice of setting the free boundary data (except $q_{0}$) to zero. FIG. 6 displays the convergence plots for a component of the $\psi_{ABCD}$ constraint $G_{AB}$ at $s=\text{constant},\;\theta=\frac{\pi}{2}$ slices. One clearly sees that there is a mode propagating in from the right boundary that stops the constraints from converging to machine precision. This happens not just to this constraint, but to all in the constraints in our system.

Now we contrast these plots to the analogous ones that implement our boundary treatment, shown in FIG. 7. Immediately one sees that these convergence plots are exceedingly better than the previous ones. We get convergence toward machine precision at the correct order and in the process have moved the constraints around $1\times 10^{8}$ closer to this. Looking at other constraints in the system we see that the problem has been overcome in all of them and this is seen to be the case for different choices of $\theta$.

Thus we have shown that without our boundary treatment, constraint violating modes of the subsidiary system propagate into the computational domain and prevent the constraints from converging to zero. In the constraint violating case we see again that the evolution system is well-posed since the constraint violation converges to non-zero value. Thus the most general case of our boundary treatment has been shown to be numerically viable.

To end this section, it is interesting to see how the “approximate” Schwarzschild radius $r=R/\Theta$ in the perturbed case differs to that of the unperturbed case, see Figure 8(b). One can see that there is in fact a deviation that gets larger as the simulation progresses. This deviation is large enough to be discerned by the eye alone, which tells us that the location of the event horizon in the perturbed case may have actually moved a significant distance. These plots also showcase the anisotropic nature of the gravitational radiation, with spheres being stretched in one direction and squeezed in the other. This can be deduced by noticing the “perturbed” Schwarzschild radius is larger or smaller than the corresponding unperturbed one, depending on which part of the sphere is considered.

V Summary

In this paper we discussed a numerical implementation of the GCFE in the space-spinor formalism. Subsequently the manifold topology was restricted to be of the form $M_{2}\times\mathbb{S}^{2}$ so that the $\eth$-calculus could be utilized. The resulting system was checked for correctness by an analytic comparison to exact solutions as well as a range of numerical tests.

Most of the evolution equations in the system are ODEs in the sense that they involve only the time derivative $\partial_{s}$ with the notable exception of the subsystem for the gravitational spinor $\psi_{ABCD}$, which is symmetric hyperbolic. We discussed the characteristics of both that system and the part of the subsidiary system concerned with the propagation of its constraints. Using maximally dissipative boundary conditions we proposed a practical way of imposing boundary conditions which are physically reasonable in the sense that there is one freely specifiable degree of freedom and no incoming constraint violating modes. This was found to be numerically stable in several different settings, indicating that the IBVP for this formulation of the conformal Einstein equations is well posed.

This newly developed framework will give us the ability to investigate a range of problems that involve the global structure and global properties of space-times. In a forthcoming paper we will report on studies of the characteristic ringing behaviour of the Schwarzschild space-time under gravitational perturbations and discuss the global issues of wave signal readout on $\mathscr{I}^{+}$ and the mass-loss due to gravitational radiation. After this problem, the most obvious next step to take is to generalize this setup to the Kerr space-time. As it is still not clear whether this space-time is stable under general perturbations, we could investigate this question from a global perspective. Investigating this problem does not require any fundamental changes to the setup for the Schwarzschild space-time, only a new initial data set is needed.

There is also the possibility of investigating the conditions required on an asymptotically flat initial data set so that the resulting vacuum solution has a regular null infinity. Friedrich Friedrich (2007) has restricted the problem to how initial data is chosen on the blowup of the point $i^{0}$ to a 2-sphere. He has conjectured that the necessary condition for a regular null infinity is that the initial data near null infinity are those induced by asymptotically conformally stationary space-times. This still remains as just a conjecture and hence it would be intriguing to probe this question numerically by evolving sets of initial data that do and do not satisfy the necessary conditions of the conjecture. There have already been numerical studies of linearly perturbed space-times which incorporate space-like infinity, see for example Beyer et al. (2012); Doulis and Frauendiener (2017, 2013), while in Frauendiener and Hennig (2014, 2017) we have studied simpler systems which show similar behaviours near space-like infinity.

Another open problem within the scope of this approach is that of the stability of Anti-de Sitter space-time. First brought to attention by Bizoń and Rostworowski in 2012 Bizoń and Rostworowski (2011), it was found that the space-time is non-linearly unstable under a particular class of perturbations. Since then time-stable periodic solutions have been discovered Maliborski and Rostworowski (2013); Balasubramanian et al. (2014); Green et al. (2015). The main issue to address in this context is the fact that a conformal geodesic on the Einstein cylinder “stalls” in the sense that it does only cover a finite interval of conformal time even for an infinite range of its parameter. This means that we need to implement a reparametrization in order to “reset” its parameter and to continue the evolution.

On the more mathematical side of things it would be desirable to have a rigorous proof that our empirical boundary treatment does in fact lead to a well-posed IBVP.

Appendix A The space-spinor formalism

Due to its compact form, the calculus of space-spinors Sommers (1980) is a very enticing formalism in which to write the GCFE. We will use definitions similar to that used in the series of papers by Frauendiener Frauendiener (1998a, b, 2000) throughout.

First, we define a time-like vector field $t^{a}$, normalized with respect to the conformal metric $g$ by $t^{a}t_{a}=2$. Writing $t^{a}$ in terms of spinors, one immediately finds the relationship

[TABLE]

which justifies the choice of normalization constant. One can now use this vector field as a map from the complex spin-space $\bar{S}^{A^{\prime}}$ onto the spin-space $S^{A}$, i.e. we can convert primed indices to unprimed ones. This map is given by

[TABLE]

For example, a spinor $\alpha_{a}=\alpha_{AA^{\prime}}$ can be mapped to $\alpha_{AB}:=t_{B}{}^{A^{\prime}}\alpha_{AA^{\prime}}$. The result can be decomposed into two terms

[TABLE]

This shows us that the trace term corresponds to the part of the spinor that has values in the direction of $t^{AA^{\prime}}$. Thus finding irreducible decompositions of space-spinors is the same as performing a 3+1 splitting. This is incredibly useful for deriving evolution and constraint equations.

We also define a complex conjugation map on the unprimed spin-space via

[TABLE]

This map has the property that for a spinor of rank $n$ we obtain

[TABLE]

We define a rank-2 spinor as real iff it is equal to the negative of its complex conjugate, in accordance with the reality of the SL(2,$\mathbb{C}$) spinors. Since for any two real rank-2 spinors $\alpha_{AB}$ and $\beta_{CD}$ their outer product should also be real, a rank-4 spinor is real if it is equal to its complex conjugate.

We can now split the covariant derivative $\nabla$ into spatial and temporal parts using the mapping $t_{B}{}^{A^{\prime}}\nabla_{AA^{\prime}}$ and its subsequent decomposition. This gives us two new derivative operators

[TABLE]

Two fundamental spinor fields can now be defined as the derivatives of $t^{a}$ with respect to these new derivative operators. We have

[TABLE]

Geometrically, the spinor field $K_{AB}$ corresponds to the acceleration vector of $t^{a}$ while $K_{ABCD}$ is related to the geometry of the distribution defined by vectors $V^{a}$ that satisfy $V^{a}t_{a}=0$ and for a time-like vector field $t^{a}$ satisfying the hyper-surface orthogonal property, it corresponds to the extrinsic curvature. They have the reality properties

[TABLE]

Note that these new derivatives operators are real in the sense that they map real spinors to real spinors, but they do not commute with our definition of complex conjugation, i.e.

[TABLE]

A similar equation holds for $K_{ABCD}$. Hence we introduce new derivative operators

[TABLE]

adjusted to commute with complex conjugation. We now have

[TABLE]

Note that the correction term $K_{AB}$ or $K_{ABCD}$ corrects only one index. So replacing the action of $D$ or $D_{AB}$ on a rank-$n$ spinor with $\partial$ or $\partial_{AB}$ respectively, result in $n$ correction terms. For example,

[TABLE]

We will denote the spin-frame spinors by $o^{A}$ and $\iota^{A}$ along with their primed counterparts.

Appendix B The complete subsidiary system

It is convenient in the calculation of the subsidiary equations to decompose the fields $K_{ABCD}$, $P_{ABCD}$ and $\gamma_{ABC}$ into irreducible pieces. The irreducible decompositions of these fields can be written as

[TABLE]

where the new spinor quantities are totally symmetric. It is also useful to decompose the constraint $U_{ABCD}=\widetilde{U}_{ABCD}+\frac{1}{2}\widetilde{U}_{AB}\varepsilon_{CD}$, where $\widetilde{U}_{ABCD}=\widetilde{U}_{(AB)(CD)}$ and $\widetilde{U}_{AB}=\widetilde{U}_{(AB)}$.

Applying the above procedure to the most general form of the constraints (29) we find the subsidiary system to be

[TABLE]

Acknowledgements.

The authors are grateful to H. Friedrich for discussions. JF wishes to thank the Department of Mathematics at the University of Oslo for hospitality while this article was written. Funding for this visit was provided by the European Research Council through the FP7-IDEAS-ERC Starting Grant scheme, Project No. 278011 STUCCOFIELDS. The authors wish to acknowledge the contribution of NeSI high-performance computing facilities to the results of this research.