Asymptotics of solutions of a hyperbolic formulation of the constraint equations

TL;DR

This paper investigates the asymptotic behavior of solutions to a hyperbolic formulation of Einstein's constraint equations, using numerical methods to analyze perturbations of Schwarzschild data and their asymptotic properties.

Contribution

It introduces a numerical approach to study the asymptotics of hyperbolic constraint solutions and highlights potential issues with asymptotic conditions in generic initial data.

Findings

Most generic initial data may violate fundamental asymptotic conditions

Numerical construction of perturbations of Schwarzschild data

Potential need to exploit free data to ensure proper asymptotics

Abstract

In this paper we consider the hyperbolic formulation of the constraints introduced by R\'acz. Using the numerical framework recently developed by us we construct initial data sets which can be interpreted as nonlinear perturbations of Schwarzschild data in Kerr-Schild coordinates and investigate their asymptotics. Our results suggest that, unless one finds a way to exploit the freedom to pick the free part of the initial data in some suitable way, generic initial data sets obtained by this method may violate fundamental asymptotic conditions.

Equations

$${}^{(3)}R+{\chi }^2-{\chi }_{ij}{\chi }^{ij}=0,D^j{\chi }_{ji}-D_i\chi =0,$$

$$h_{ij}$$

$${\chi }_{ij}$$

$$\varrho :=\sum _{i=1}^3({\tilde{x}}^i{)}^2\to \infty .$$

$$h_{ij}$$

$${\chi }_{ij}$$

$${\hat{n}}_i=\hat{N}{D}_i\rho$$

$${\hat{\gamma }}_{ij}=h_{ij}-{\hat{n}}_i{\hat{n}}_j,$$

$${\hat{\gamma }}_j^i={\delta }_j^i-{\hat{n}}^i{\hat{n}}_j$$

$${\chi }_{ij}=\text{boldmathκ}{\hat{n}}_i{\hat{n}}_j+2{\hat{n}}_{(i}{\mathbf{k}}_{j)}+{\mathbf{K}}_{ij},$$

$$\text{boldmathκ}={\hat{n}}^i{\hat{n}}^j{\chi }_{ij},{\mathbf{k}}_j={\hat{\gamma }}_j^n{\hat{n}}^m{\chi }_{nm},{\mathbf{K}}_{ij}={\hat{\gamma }}_i^n{\hat{\gamma }}_j^m{\chi }_{nm}.$$

$${\mathbf{K}}_{ij}={\mathring{K}}_{ij}+\frac{1}{2}{\hat{\gamma }}_{ij}\mathbf{K},$$

$$\mathbf{K}={\mathbf{K}}_l^l,{\mathring{K}}_{ij}{\hat{\gamma }}^{ij}=0,$$

$$D_i{\hat{n}}_j={\hat{K}}_{ij}+{\hat{n}}_i{\dot{\hat{n}}}_j$$

$${\hat{K}}_{ij}={\hat{\gamma }}_i^n{\hat{\gamma }}_j^m{\hat{K}}_{mn},{\dot{\hat{n}}}_k={\hat{n}}^i{D}_i{\hat{n}}_j.$$

$${\mathcal{L}}_{\hat{n}}\mathbf{K}-{\hat{D}}^l{\mathbf{k}}_l$$

$${\mathcal{L}}_{\hat{n}}{\mathbf{k}}_i-\frac{2}{\mathbf{K}}{\mathbf{k}}^l{\hat{D}}_i{\mathbf{k}}_l+\frac{\text{boldmathκ}}{\mathbf{K}}{\hat{D}}_i\mathbf{K}$$

$$\hat{K}={\widehat{K^i}}_i,{\kappa }_0={}^{(3)}R-{\mathring{K}}_{ij}{\mathring{K}}^{ij},\text{boldmathκ}=\frac{1}{2\mathbf{K}}\left(2{\mathbf{k}}^i{\mathbf{k}}_i-\frac{{\mathbf{K}}^2}{2}-{\kappa }_0\right).$$

$${\rho }^i=\hat{N}{\hat{n}}^i+{\hat{N}}^i,$$

$${\hat{\gamma }}_{ij},\hat{N},{\hat{N}}^i,{\mathring{K}}_{ij}\text{and}{\hat{K}}_{ij},$$

$$\text{boldmathκ}\mathbf{K}<0$$

$$\left(\begin{array}{cc}{\hat{n}}^j& -{\hat{\gamma }}^{jl}\\ \frac{\text{boldmathκ}}{\mathbf{K}}{{\delta }_i}^j& {\hat{n}}^j{{\delta }_i}^l-\frac{2}{\mathbf{K}}{\hat{\gamma }}^{ml}{{\delta }_i}^j{\mathbf{k}}_m\end{array}\right)D_j\left(\begin{array}{c}\mathbf{K}\\ {\mathbf{k}}_l\end{array}\right).$$

$$\left(\begin{array}{cc}-\frac{\text{boldmathκ}}{\mathbf{K}}& 0\\ 0& {\gamma }^{ik}\end{array}\right).$$

$$\left(\begin{array}{cc}-\frac{\text{boldmathκ}}{\mathbf{K}}& 0\\ 0& {\gamma }^{ik}\end{array}\right)\cdot \left(\begin{array}{cc}{\hat{n}}^j& -{\hat{\gamma }}^{jl}\\ \frac{\text{boldmathκ}}{\mathbf{K}}{{\delta }_i}^j& {\hat{n}}^j{{\delta }_i}^l-\frac{2}{\mathbf{K}}{\hat{\gamma }}^{ml}{{\delta }_i}^j{\mathbf{k}}_m\end{array}\right){\hat{n}}_j$$

$${\left(\text{boldmathκ}\mathbf{K}\right)|}_{\rho ={\rho }_0}<0.$$

$$m^i:=\frac{1}{\sqrt{2}}\left({\partial }_{\theta }^i-\frac{\text{i}}{\sin \theta }{\partial }_{\phi }^i\right);$$

$${\omega }_k:=\frac{1}{\sqrt{2}}\left(\mathrm{d}{\theta }_k+\text{i}\sin \theta \mathrm{d}{\phi }_k\right).$$

$${\hat{\gamma }}_{ij}=2{\hat{\gamma }}_{12}{\omega }_{(i}{\overline{\omega }}_{j)}+{\hat{\gamma }}_{11}{\omega }_i{\omega }_j+{\hat{\gamma }}_{22}{\overline{\omega }}_i{\overline{\omega }}_j.$$

$$f(\theta ,\phi )=\sum _{l=|s|}^{\infty }\sum _{m=-l}^l{}_s{a}_{lm}{}_s{Y}_{lm}(\theta ,\phi ),$$

$$⟨{}_s{Y}_{l_1{m}_1}|{}_s{Y}_{l_2{m}_2}⟩:=\underset{{\mathbb{S}}^2}{\int }{}_s{Y}_{l_1{m}_1}(\theta ,\phi ){}_s{\overline{Y}}_{l_2{m}_2}(\theta ,\phi )\mathrm{d}\Omega ={\delta }_{l_1{l}_2}{\delta }_{m_1{m}_2},$$

$${}_s{a}_{lm}=⟨f|{}_s{Y}_{l_2{m}_2}⟩.$$

$$ðf:={\partial }_{\theta }f-\frac{\text{i}}{\sin \theta }{\partial }_{\phi }f-s\cot \theta f,\overline{ð}f:={\partial }_{\theta }f+\frac{\text{i}}{\sin \theta }{\partial }_{\phi }f+s\cot \theta f.$$

Full text

1 Introduction

The initial value problem plays a hugely important role on the seemingly eternal quest through the largely unmapped terrain of solutions of Einstein’s equations. In its most common formulation, see [21, 1, 36], one performs a $3+1$-split and thereby decomposes the spacetime and all corresponding geometric tensor fields into spatial and timelike components. As a consequence Einstein’s equations naturally split into constraints and evolution equations. Thanks to the ground-breaking work by Choquet-Bruhat et al. [19, 14] we know that each solution of the constraints determines a unique solution of the Einstein’s equations – a so-called maximal globally hyperbolic development – in which the initial data set arises as the induced geometry of some spacelike surface. In all of what follows, any triple $(\Sigma,h_{ij},\chi_{ij})$ of a $3$-dimensional differentiable manifold $\Sigma$, a Riemannian metric $h_{ij}$ and a smooth symmetric tensor field $\chi_{ij}$ on $\Sigma$ is called an initial data set if it satisfies the (vacuum) constraint equations

[TABLE]

everywhere on $\Sigma$. The quantities ${}^{(3)}R$ and $D_{i}$ are the Ricci scalar and the covariant derivative operator associated with $h_{ij}$, respectively, and $\chi=\chi^{j}_{\ j}$. For this whole paper we agree that operations involving abstract spatial indices $i,j,\ldots$ are performed with $h_{ij}$.

One easily confirms that the constraints are under-determined in that there are far more degrees of freedom in the two tensor fields $h_{ij}$ and $\chi_{ij}$ than there are constraint equations. The constraints therefore leave us with many different ways to separate the free part of the initial data set from the part which we want to determine by solving the constraints. There are thus many quite distinct formulations of the same equations. It turns out that the particular PDE character and hence the PDE problem naturally associated with it (boundary value problem, initial (boundary) value problem etc.) largely depends on such choices. The most common formulations of the constraints in the literature are based on the Lichnerowicz-York conformal method, see for instance [4, 5, 26]. There are several variations, but all these formulations have in common that they lead to elliptic equations which one then attempts to solve as boundary value problems.

While this is mathematically convenient, and also seemingly natural as a procedure to construct initial data sets, it was realized that the resulting initial data sets can contain “junk radiation” and other non-physical features. In numerical relativity one commonly works with initial data sets for which the spatial metric, which is part of the free data in this formulation up to a conformal factor, is conformally flat. However, it was shown [22] that rotating black hole spacetimes do not admit conformally flat slices. This discrepancy is one reason for the observed undesired spurious gravitational wave content (see, for instance [1, 16]). In the literature one can find several attempts to go beyond conformal flatness within the conformal formulation of the constraints, for example, by Matzner, Huq and Shoemaker [27] where the free part of the data is determined from a superposition of boosted Kerr–Schild metrics [16].

In general, however, it is unclear how to address this obvious discrepancy between physical meaningfulness and mathematical convenience given that geometric and gauge degrees of freedom are intermixed in a highly nonlinear and largely obscure manner. One possible approach is to study several different formulations of the constraints, construct solutions and then compare the results in order to see if there is one method which stands out both mathematically and physically. One would therefore like to consider approaches which differ significantly from the standard one above.

In fact, one sees easily that it is not necessary to formulate the constraints as purely elliptic PDEs. Breaking the restrictive rigidity of elliptic PDEs was for instance a crucial step in the work by Corvino [17] who casted the constraints partly into a different form. Early works where the constraints were formulated as fully non-elliptic evolution problems – i.e., initial value problems as opposed to boundary value problems – are [12, 3]. Recently, Rácz has brought this idea to a new level in a series of papers [31, 32, 33]. He proposes two different formulations of the constraints, one where the equations are strongly hyperbolic, and one where they take the form of a parabolic-hyperbolic system (at least under certain restrictions). In the present paper we focus on the former whose details are discussed in Section 2.1. The fact that the constraints here are solved as an initial value problem breaks with standard methodologies and interpretations (and maybe even ideologies) which had been developed for the constraints over many decades. Similar to the $3+1$-decomposition of spacetime and associated evolutions with respect to some foliation with spacelike $3$-dimensional surfaces, it is the initial hypersurface $\Sigma$ which is decomposed in a $2+1$-fashion here. The constraints are then solved as an evolution problem with respect to a foliation with $2$-dimensional surfaces. For most of this paper we have in mind that these $2$-dimensional hypersurfaces are spheres, that the initial one of these spheres is finite and that the constraints are then evolved outwards radially towards infinity. This sounds convenient, for example, because numerical codes, in particular [7, 6, 8, 9], which had been developed to solve the evolution equations of Einstein’s equations, can be used without major changes here as well. However, it also has an immediate draw-back: In contrast to the boundary value problem we only have little control over the behavior of the solutions at infinity. In particular, it is, as a matter of principle, impossible to impose that the resulting initial data set should be asymptotically flat, asymptotically hyperbolic or whatever asymptotic behavior is of interest. As a consequence of this lack of control, one may therefore end up with initial data sets which are physically irrelevant. Here we focus on exactly this problem and analyze the asymptotics of non-trivial initial data set families obtained by this method. We restrict to the asymptotically flat case.

There are several notions of asymptotic flatness in the literature most of which are motivated by either mathematical or physical considerations (or by both, see [24, 37]). Asymptotically flat initial data sets are supposed to represent a moment of time for isolated systems whose gravitational field decays towards infinity in an appropriate manner. The strongest notion of asymptotic flatness, which we shall consider here, is taken from [18]. We say that an initial data set $(\Sigma,h_{ij},\chi_{ij})$ is strongly asymptotically flat if $\Sigma$ is diffeomorphic to $\mathbb{R}^{3}$ (possibly minus a ball of finite size) and if there exist coordinates $\{\tilde{x}^{i}\}$ on $\Sigma$ such that the components of $h_{ij}$ and $\chi_{ij}$ with respect to these coordinates satisfy444We use the $O$ symbol rather informally in the usual sense $f=O(g)\iff f/g<\infty$ in the relevant limit. Similarly, $f=o(g)\iff f/g\to 0$. In this paper, we avoid the technicalities of norms which must be used to make the definition of the $O$-symbols precise. In fact, in order to make the above notions of asymptotic flatness precise and physically meaningful, the $O$-symbol must be defined with respect to a norm which does not only control the decay of the fields themselves, but also guarantees a sufficient decay speed of an appropriate number of derivatives. If this is the case, “asymptotic flatness” in the sense above can be shown to imply that the curvature tensor also decays at infinity with some known rate. The interested reader can find the details in the references above.

[TABLE]

in the limit

[TABLE]

Here $\delta_{ij}$ denotes the Euclidean metric in Cartesian coordinates and $M\geq 0$ is the ADM mass. This notion of asymptotic flatness is motivated by the fact that it allows to conformally compactify the initial data set in a way where $\varrho=\infty$ is represented by a single point $i^{0}$ on the compactified manifold and the fields $h_{ij}$ and $\chi_{ij}$ extend smoothly through $i^{0}$ after certain conformal rescalings. This is necessary to apply the extraordinarily detailed analysis by Friedrich [20] which yields conditions for which the corresponding solution of the (full) Einstein equations has a smooth conformal compactification. This is a crucial physical property of spacetimes describing isolated systems. In fact, the entire theory of gravitational radiation in general relativity (beyond the linearized level) relies on the existence of such structures.

We refer to such initial data as strongly asymptotically flat because there are also weaker notions in the literature. According to [2] the ADM mass of an initial data set is defined if the initial data set is weakly asymptotically flat, i.e., if

[TABLE]

Spacetime developments determined by such weakly asymptotically flat initial data sets close to the trivial data set were studied in quite some detail in [11, 10]. The assumption of weak asymptotic flatness in these works was a significant improvement over the hypotheses in the original ground-breaking work by Christoudoulou and Klainerman [15] regarding the nonlinear stability of Minkowski spacetime. In any case, while the extremely difficult and subtle achievements in these works allow to control the decay of fields at infinity just enough to prove stability, only weak conclusions can be drawn regarding the global conformal structure and therefore regarding gravitational radiation.

In summary it is fair to say that it depends on the particular situation, application and interest whether one wants to consider the strong notion of asymptotic flatness or rather the the weak version above. We are therefore motivated to investigate both here.

The main result of our paper is evidence that, at least in the particular setting we consider here, asymptotically flat initial data sets are unstable (in both the weak and the strong sense) within the family of solutions of the hyperbolic formulation of the constraints. It is therefore likely to be difficult to construct asymptotically flat initial data sets by means of the hyperbolic formulation of the constraints by Rácz in practice because both the free data and the initial data of the initial value problem of the constraints may need to be fine-tuned in some way. We argue that this conjectured instability has its origin in the existence of spherically symmetric slices in the Schwarzschild spacetime with “the wrong mass” which are not asymptotically flat. All this is explained in detail in Section 4. In any case, our results are not conclusive enough to claim that Rácz’s formulation of the constraints is “useless in practice”. After all, we only study a particular setting here, and it may be possible that the freedom in choosing the free part of the data can be exploited in a clever way to solve this issue. Our results do confirm however that, without extra care, the lack of a-priori control of the asymptotics of solutions for this formulations of the constraints can have negative consequences and can in particular lead to non-physical initial data sets.

We emphasize that our paper is mostly concerned with the asymptotics of solutions of the constraints at spatial infinity. The interesting recent findings in [39] therefore have no direct consequences for us.

Our paper is organized as follows. In Section 2.1 we introduce the $2+1$-decomposition and the hyperbolic formulation of the constraint equations. A particularly interesting setting is obtained when the leaves of the foliations are chosen as $2$-dimensional spheres. This is interesting physically, but also turns out to be convenient technically because one can work with the so-called eth-formalism. All this is discussed in Section 2.2. In Section 2.3 we consider the Kerr and Schwarzschild solutions in Kerr-Schild coordinates. These yield natural families of exact solutions of the constraints which will be a useful foundation for the analysis in the following sections. Section 3 is devoted to numerical tests. The main part is Section 4 where we analyze the asymptotics of a certain class of solutions of the constraints by means of a combination of analytic and numerical methods and thereby derive the main conclusions of our paper.

2 The vacuum constraint equations as a hyperbolic system

2.1 Hyperbolic formulation of the constraints

Let us now discuss the particular hyperbolic formulation of the constraints that we consider in this paper using more or less the original notation introduced in [33]. Given any initial data set $(\Sigma,h_{ij},\chi_{ij})$ as described above, one starts from a $2+1$-decomposition of $\Sigma$ (in full analogy with standard $3+1$-decompositions of spacetimes). One picks a smooth function $\rho:\Sigma\rightarrow\mathbb{R}$ such that all the $\rho=const$-level sets are smooth hypersurfaces $S_{\rho}$ in $\Sigma$ and such that the family of all these hypersurfaces is a foliation of $\Sigma$. There is therefore a smooth lapse function $\hat{N}$ on $\Sigma$ so that

[TABLE]

is a unit normal to the surfaces $S_{\rho}$. Observe that this lapse $\hat{N}$ has nothing to do with the lapse of the $3+1$-decomposition of spacetime and, in contrast to the $3+1$-decomposition, $\hat{n}_{i}$ is clearly always spacelike here (because the metric $h_{ij}$ is Riemannian). Since $\hat{N}$ does not vanish as a consequence of the above assumptions, we assume, without loss of generality, that $\hat{N}>0$ everywhere. The induced metric on $S_{\rho}$ is

[TABLE]

and

[TABLE]

is therefore the map which projects any tensor field on $\Sigma$ orthogonally to a tensor field which is intrinsic to the surfaces $S_{\rho}$. Any tensor field, which is annihilated by contractions with $\hat{n}_{i}$ or $\hat{n}^{i}$ on each of its indices, is referred to as intrinsic (to the surfaces $S_{\rho}$) in all of what follows.

Let us now perform a standard decomposition of tensor fields on $\Sigma$ with respect to the normal and tangential directions of $S_{\rho}$. As a result, all of the relevant geometric quantities will be either scalars or intrinsic tensor fields. We stress that such decompositions are always unique. For $\chi_{ij}$ this gives

[TABLE]

with

[TABLE]

Clearly, the fields $\mathbf{k}_{j}$ and $\mathbf{K}_{ij}$ are intrinsic and $\mathbf{K}_{ij}$ is symmetric. The quantity $\mathbf{K}_{ij}$ can be decomposed further as

[TABLE]

where

[TABLE]

and where $\mathring{K}_{ij}$ is symmetric. Next we decompose $D_{i}\hat{n}_{j}$:

[TABLE]

where

[TABLE]

Both $\hat{K}_{ij}$ and $\dot{\hat{n}}_{k}$ are purely intrinsic and $\hat{K}_{ij}$ is symmetric. By means of standard calculations, see for example [33], one demonstrates that the constraint equations (1) can be expressed in terms of these quantities as follows

[TABLE]

where $\hat{D}$ is the covariant derivative associated with $\hat{\gamma}_{ij}$, and,

[TABLE]

In order to interpret Eqs. (12) – (13) with Eq. (14) as a system of evolution equations now, we pick, in addition to the above, a smooth vector field $\rho^{i}$ normalized as $\rho^{i}D_{i}\rho=1$. Effectively, solving this system as an evolution problem means that we shall integrate along this vector field; see below. We find

[TABLE]

where $\hat{N}^{i}$ is an intrinsic shift vector field, i.e., $\hat{N}^{i}=\hat{\gamma}^{i}_{\ j}\rho^{j}$; recall the definition of the lapse $\hat{N}$ in Eq. (7). Let us now introduce the following terminology:

**Free data: **

The fields

[TABLE]

introduced above are considered as freely specifiable data on $\Sigma$, i.e., free data, for the system Eqs. (12) – (14). Observe that the full $3$-dimensional metric $h_{ij}$ is determined by these free data alone. In particular, one can therefore calculate ${}^{(3)}R$ in Eq. (14) and $\dot{\hat{n}}_{k}$ according to Eq. (11) before the system is solved. Given $\hat{N}$ and $\hat{N}^{i}$, the Lie derivatives with respect to $\hat{n}^{i}$ in Eqs. (12) and (13) can be written as derivatives along $\rho^{i}$ (“evolution derivatives”) and derivatives tangential to $S_{\rho}$. This justifies the intuition of considering our system as an evolution problem where one integrates along the field $\rho^{i}$.

**Unknowns: **

The fields $\mathbf{K}$ and $\mathbf{k}_{i}$ are considered as the unknowns of the system Eqs. (12) – (14). Once a solution has been found, one can calculate $\kappa$ from Eq. (14), and therefore the complete initial data set $h_{ij}$ and $\chi_{ij}$.

**Hyperbolicity condition: **

The inequality

[TABLE]

is referred to as the hyperbolicity condition.

It is shown in [33] that for any smooth free data Eq. (16) on $\Sigma$, the system Eqs. (12) – (13) with Eq. (14) is a quasilinear strongly (in fact, strictly) hyperbolic system with respect to the foliation of surfaces $S_{\rho}$ for the unknowns $\mathbf{K}$ and $\mathbf{k}_{i}$ provided that (17) holds. Therefore, it is natural to refer to such an inequality as the hyperbolicity condition as above. Observe that the principal part of the system can be written as

[TABLE]

Eqs. (12) – (13) can therefore always be symmetrized [33] by multiplying it from the left with the matrix

[TABLE]

Since the matrix

[TABLE]

is always positive definite provided Eq. (17) holds, Eqs. (12) – (13) with Eq. (14) can indeed be transformed into a quasilinear symmetric hyperbolic system with respect to the foliation of surfaces $S_{\rho}$. In particular, this implies that given any smooth free data as in Eq. (16) on $\Sigma$, the initial value problem of Eqs. (12) – (13) with Eq. (14) is well-posed for the class of all smooth initial data $\left.\mathbf{K}\right|_{\rho=\rho_{0}}$ and $\left.\mathbf{k}_{i}\right|_{\rho=\rho_{0}}$ prescribed on any smooth initial surface $S_{\rho_{0}}$ for which

[TABLE]

This means that for any such data (free data and initial data) there exists a unique smooth solution $(\mathbf{K},\mathbf{k}_{i})$ of Eqs. (12) – (13) with Eq. (14), at least on some (possibly small) interval $[\rho_{0},\rho_{1})$.

Suppose for the moment that $\Sigma$ can be identified with $\mathbb{R}^{3}$ (possibly outside some compact subset) and that the limit $\rho\rightarrow\infty$ describes the geometric asymptotics at spacelike infinity. The well-posedness of the initial value problem of our system, which we have just established under the conditions above, is certainly crucial. However, there are two further questions: (i) Do any of the solutions extend to $\rho\rightarrow\infty$? After all, the hyperbolicity condition (17) involves the unknowns and could therefore break down at any point of the evolution. (ii) What is the asymptotic geometric behavior of the resulting initial data set? The answers to these questions are completely unclear and will (at least partly) be addressed in this article.

2.2 Foliation by spheres and the eth-formalism

For the rest of this paper we shall restrict to the case that $\Sigma$ can be identified with $\mathbb{R}^{3}$ (possibly minus some compact subset) and that the surfaces $S_{\rho}$ introduced above are diffeomorphic to $2$-spheres for all $\rho$ larger than some $\rho_{0}>0$. The initial value problem of Eqs. (12) – (13) with Eq. (14) can then be understood as an evolution problem with “time” $\rho$ and where the “spatial topology” is $\mathbb{S}^{2}$. Some of the technical complications in solving this evolution problem both analytically and numerically arise from the fact that $\mathbb{S}^{2}$ cannot be covered by a single regular coordinate chart. The description of tensor fields on $\mathbb{S}^{2}$ in terms of any such single coordinate chart breaks down at some “spatial” point. This problem is often called the pole problem because in standard polar coordinates for $\mathbb{S}^{2}$ these issues appear at the poles.

This problem is resolved in the so-called $\eth$-formalism (or eth-formalism) where all formally singular terms associated with the coordinate singularity at the poles are absorbed into geometrically motivated regular differential operators. Without going into the details (see [30, 7, 6, 8]), here is a brief summary. One first introduces coordinates $(\rho,\theta,\varphi)$ on $\Sigma$ for $\rho>\rho_{0}$ which are adapted to the foliation, i.e., the field $\rho^{i}$ in Eq. (15) is identified with the coordinate field $\partial_{\rho}^{i}$ and $(\theta,\varphi)$ are standard polar coordinates on each $S_{\rho}$ (diffeomorphic to $\mathbb{S}^{2}$). Given these, one introduces the frame $(m^{i},\overline{m}^{i})$ (defined almost everywhere on $S_{\rho}$ for any $\rho$) where

[TABLE]

complex conjugation is denoted by a bar and i is the imaginary unit. The corresponding dual frame is $(\omega_{k},\overline{\omega}_{k})$ where

[TABLE]

Any tensor field which is intrinsic to $S_{\rho}$ can be expressed in terms of this frame at each value of $\rho>\rho_{0}$. Our convention is that the index $1$ refers to the projection onto the $m^{i}$ or $\omega_{i}$ direction, while the index $2$ is associated analogously with $\overline{m}^{i}$ or $\overline{\omega}_{i}$. For example, we therefore have (at any fixed value of $\rho$)

[TABLE]

Observe that $\hat{\gamma}_{12}$ must therefore be a real function while $\hat{\gamma}_{11}=\overline{\hat{\gamma}_{22}}$ are complex functions on $\mathbb{S}^{2}$. Note, that the metric of the unit-sphere is obtained with $\hat{\gamma}_{11}=0$ and $\hat{\gamma}_{12}=1$. Note further, that the metric remains unchanged under frame rotations $m^{i}\mapsto\text{e}^{\text{i}\alpha}\,m^{i}$ for arbitrary functions $\alpha:\Sigma\to\mathbb{R}$.

Given this frame, the spin-weight formalism (see the references above) assigns to any component of a smooth intrinsic tensor field with respect to this frame a well-defined integer spin-weight. The spin-weight describes the transformation behavior of the components of tensor fields under pointwise frame rotations of the frame $(m^{i},\overline{m}^{i})$ (and correspondingly of $(\omega_{i},\overline{\omega}_{i})$) as defined above. By convention, $m^{i}$ has the spin-weight $1$, $\overline{m}^{i}$ the spin-weight $-1$, $\omega_{i}$ the spin-weight $-1$ and $\overline{\omega}_{i}$ the spin-weight $1$. Correspondingly, the function $\hat{\gamma}_{12}=\hat{\gamma}_{ij}m^{i}\overline{m}^{j}$ has the spin-weight [math], $\hat{\gamma}_{11}=\hat{\gamma}_{ij}m^{i}m^{j}$ the spin-weight $2$ and $\hat{\gamma}_{22}=\hat{\gamma}_{ij}\overline{m}^{i}\overline{m}^{j}$ the spin-weight $-2$ etc. Any sufficiently regular function $f$ on $\mathbb{S}^{2}$ (and therefore on $S_{\rho}$ for any fixed $\rho>\rho_{0}$) with spin-weight $s$ can be expanded in terms of spin-weighted spherical harmonics (SWSH) ${}_{s}Y_{lm}$ as follows

[TABLE]

for complex coefficients ${}_{s}a_{lm}$. Thus, a spin-weighted quantity can be represented in terms of local coordinates $(\theta,\varphi)$ (on the left) of in terms of its expansion coefficients ${}_{s}a_{lm}$ (on the right).

Due to the orthonormality relation

[TABLE]

the coefficients ${}_{s}a_{lm}$ in Eq. (19) can be calculated as

[TABLE]

If $f$ is a function of spin-weight $s$ on $\mathbb{S}^{2}$, the eth-operators $\eth$ and $\bar{\eth}$ are defined as

[TABLE]

Using Eqs. (22) and (18), we can therefore express directional derivatives along the frame vectors $(m^{i},\overline{m}^{i})$ in terms of the eth-operators as

[TABLE]

The operator $\eth$ raises the spin-weight by one while $\bar{\eth}$ lowers it by one. Of particular importance are the identities

[TABLE]

which can also be used to define the spin-weighted spherical harmonics ${}_{s}Y_{lm}$ for any integer $s$ from the standard spherical harmonics $Y_{lm}={}_{0}Y_{lm}$. It is worthwhile to point out that the power of the $\eth$-formalism arises from the fact that the SWSH and the $\eth$ and $\bar{\eth}$ operator are globally defined on $\mathbb{S}^{2}$. We have introduced them above in terms of their local representations in the usual polar coordinates with respect to the frame (18). This globality and, in particular the raising and lowering properties (24), allow us to compute derivatives and integrals of spin-weighted quantities very easily using fast transformations between their representations with respect to the SWSH and the local coordinates.

In order to apply this framework to Eqs. (12) – (14) one proceeds as follows. First all intrinsic fields in these equations are expressed in terms of the frame above. Each quantity in the equations is then described in terms of functions on $(\rho_{0},\infty)\times\mathbb{S}^{2}$ (i.e., they depend on $\rho$, $\theta$ and $\varphi$) each of which has a well-defined spin-weight at each fixed $\rho$ which is constant in $\rho$. Next, all directional derivatives along $m^{i}$ and $\overline{m}^{i}$ are expressed in terms of $\eth$ and $\bar{\eth}$ according to Eq. (23). After simplifying the resulting equations algebraically, one obtains a system of equations, where all terms in each equation have the same spin-weight and all terms are explicitly regular (i.e., the pole problem disappears). Because the full set of equations which one obtains are quite lengthy, we will not write them down now; see [35]. For brevity we only write them down in the special case below in Eq. (29).

2.3 Some explicit solutions of the constraints

This subsection is devoted to some explicit solutions of the constraints Eqs. (12) – (14). On the one hand, we will use these to test our numerical code in Section 3. On the other hand, however, they will also be the foundation for our analysis of non-trivial nonlinearly perturbed initial data in Section 4.

Kerr-Schild initial data.

We start with the Kerr spacetime [25]. As noted before, Eq. (17) is incompatible with time-symmetric slices. The initial data sets obtained by $t=const$-slices of the Kerr metric in Boyer-Lindquist coordinates are therefore not compatible with our formulation of the constraints. This problem can be reconciled by considering instead $t=const$-surfaces of the Kerr metric in Kerr-Schild coordinates. As these slices are not time-symmetric, but asymptotically flat, they shall turn out to be useful for our studies here.

In Kerr-Schild coordinates $(t,x,y,z)$, the Kerr metric with mass $M$ and angular momentum parameter $a$ takes the form [28, 34]

[TABLE]

where

[TABLE]

the field $\eta_{\mu\nu}$ is the Minkowski metric in standard Cartesian coordinates $(t,x,y,z)$, and $r$ is the unique positive real solution of the equation

[TABLE]

The limit $r\rightarrow\infty$ corresponds to spacelike infinity. Observe that $l_{\mu}$ is null both with respect to $g_{\mu\nu}$ and $\eta_{\mu\nu}$.

The induced metric $h_{ij}$ and second fundamental form $\chi_{ij}$ on any $t=const$-surface $\Sigma$ in this spacetime is an initial data set. It can therefore be expressed as in Sections 2.1 and 2.2. The $\rho=r=const$-surfaces on $\Sigma$ are diffeomorphic to spheres and provide a foliation of $\Sigma$ for any $r>0$. Therefore, in all of what follows we will write $r$ instead of $\rho$.

The components of the intrinsic metric on $\mathbb{S}^{2}$ are

[TABLE]

and $\hat{\gamma}_{22}=\overline{\hat{\gamma}_{11}}$. The lapse function takes the form

[TABLE]

and the components of the shift vector are

[TABLE]

and $\hat{N}^{2}=\overline{\hat{N}^{1}}$. The components of $\mathbf{K}_{ij}$ are

[TABLE]

with

[TABLE]

and $\mathbf{K}_{11}=\overline{\mathbf{K}_{22}}$. Via Eq. (9), this allows us to calculate the components of $\mathring{K}_{ij}$ and of $\mathbf{K}$. The latter function and the components of $\mathbf{k}_{i}$ are given by

[TABLE]

where

[TABLE]

and $\mathbf{k}_{2}=\overline{\mathbf{k}_{1}}$.

The Schwarzschild case and nonlinear perturbations.

The special case of Kerr-Schild Schwarzschild initial data is obtained from the above by setting $a=0$ for any $M\geq 0$. In this case we have

[TABLE]

Observe that the $S_{r}$ level sets are therefore round spheres with radius $r$. We calculate from Eq. (14) that

[TABLE]

Moreover, we have

[TABLE]

In order to define nonlinear perturbations of these exact Kerr-Schild Schwarzschild data now we proceed as in [34]: We keep the exact Schwarzschild quantities Eq. (26) as the free data, but then consider $\mathbf{K}$ and $\mathbf{k}_{1}$ as the general unknowns of Eqs. (12) – (14), which thereby take the form

[TABLE]

with $\mathbf{k}_{2}=\overline{\mathbf{k}_{1}}$ and Eq. (27); in applying the eth-operators in these equations, we must make use of the spin-weights $1$ for $\mathbf{k}_{1}$ and $-1$ for $\mathbf{k}_{2}$ consistently. One can show easily [34] that the hyperbolicity condition is satisfied for Eq. (29) for all initial data sufficiently close to data determined by Eq. (28) on any initial sphere $r=r_{0}>0$. The initial value problem is therefore well-posed. The solutions of Eq. (29) given by different initial data for $\mathbf{K}$ and $\mathbf{k}_{1}$ are interpreted as different nonlinear perturbations of the Kerr-Schild Schwarzschild initial data set.

One can easily find an explicit family of such nonlinear perturbations by assuming that $\mathbf{k}_{1}=\mathbf{k}_{2}=0$ and that $\mathbf{K}$ only depends on $r$. In this case, Eq. (29) reduces to the single ordinary differential equation

[TABLE]

Using Eq. (27) we obtain the following family of solutions

[TABLE]

where $C$ is an integration constant. We impose the restriction $C\geq 0$ so that $\mathbf{K}$ is real for all large values of $r$. The case of exact Kerr-Schild Schwarzschild data corresponds to the case $C=0$, and, by choice of the time orientation, to the overall negative sign. In Section 4.1 we discuss the significantly different properties of the initial data sets obtained for $C=0$ and $C>0$.

3 Numerical implementation and tests

Since the general hyperbolic constraint equations can be formulated completely in the eth-formalism under the assumptions in Section 2 (this is not restricted to the special case in Eq. (29)), the pseudo-spectral numerical framework developed in [7, 6, 8, 9] applies naturally. The idea is that the unknowns, all of which have a well-defined spin-weight in this formalism, are expanded in terms of spin-weighted spherical harmonics as in Eq. (19) at any value $r=\rho$. The complex coefficients in the expansions hence only depend on $r$. The “spatial derivatives” in the equations can be calculated explicitly using Eqs. (24). In this paper we restrict to the axially symmetric case where we assume that all the unknowns are constant along $\varphi$. In this case, we can use the highly efficient and accurate code developed in [8]. We refer the reader to these references for details.

The main purpose of this section is to discuss numerical tests using the exact Kerr-Schild Kerr data solution of the constraints for $M=1$ and $a=0.5$ considered in Section 2.3. For comparing the numerical and the exact solution, we define

[TABLE]

as our error measure. The quantities $\mathbf{K}^{(num)}$ and $\mathbf{k}_{1}^{(num)}$ represent the numerical solutions while $\mathbf{K}$ and $\mathbf{k}_{1}$ represent the exact solution given in Section 2.3. The norm555For simplicity we are not using the correct $\mathbb{S}^{2}$-volume element to approximate this “norm”; we find that the inclusion of the correct volume element here does not change the results significantly. $\|\cdot\|_{L^{2}(\mathbb{S}^{2})}$ is numerically approximated by

[TABLE]

for any function $f$ sampled on the grid of $N$ points in the $\theta$ direction; recall that due to axial symmetry we do not need to sample in the $\varphi$ direction.

In Fig. 2 we carry out a convergence test with RK4 (the 4th-order Runge-Kutta scheme) as the numerical integrator and $dr$ as the size of the “time” step. We see clearly that the numerical solutions converge in $4$th order which suggests that the numerical error is dominated by the “time” discretization. The “spatial” discretization error is expected to be negligible in this regime because the spectral discretization in space is expected to be highly accurate. All this is in full agreement with other numerical evolutions we have performed with the same code in previous investigations [8, 9].

Notice that Fig. 2 restricts to small values of $r$. Depending on how we implemented the equations practically in our code, we were surprised to find the following numerical instability in our attempts to extend the evolutions to large $r$ values. We first attributed this to some sort of aliasing instability caused by non-linear terms in the equations. However, after a judicious study, we recognized terms of the form $\kappa_{0}/\mathbf{K}$ in the equation as the source. In fact both quantities $\kappa_{0}$ and $\mathbf{K}$ decay in the limit $r\rightarrow\infty$; Eq. (27) yields the decay rate of $r^{-4}$ for $\kappa_{0}$ while $\mathbf{K}$ decays like $r^{-2}$ in the case of exact Kerr-Schild Kerr data according to Eq. (25). Later we will see that $\mathbf{K}$ may decay slower than that for more general classes of solutions. But within the particular family of solutions considered in this paper it will always approach zero at infinity. As a consequence of the finite accuracy of numbers in computers, any expression which involves $\kappa_{0}/\mathbf{K}$ is therefore potentially problematic.

In Fig. 2 we investigate this issue and demonstrate that the numerical error depends strongly on the way we numerically evaluate the $\kappa_{0}/\mathbf{K}$-terms. We perform three different runs each of which has a different way to estimate these problematic terms. All these runs were performed with the Runge-Kutta-Fehlberg (RKF) method as the numerical integrator with absolute tolerance $10^{-12}$ and $N=16$. Recall that in general $\kappa_{0}$ is given by Eq. (14) which, in the case of Kerr-Schild Kerr data leads to an exact expression which is too long to be printed out here. Since we designed our code to be as general as possible, we wrote a routine that calculates the term ${}^{(3)}R$ in Eq. (14) from any given free data $\hat{\gamma}_{ij}$, $\hat{N}$ and $\hat{N}^{i}$ fully numerically. While this makes it in principle possible to perform runs with arbitrary free data without changing the code, it turns out that the associated numerical errors are the source of the above mentioned instability. In fact, this routine is used in the run which yields the continuous (second) curve in Fig. 2, for which the numerical error is clearly unstable. If we instead evaluate the terms $\kappa_{0}/\mathbf{K}$ in the equations by calculating ${}^{(3)}R$ with exact computer algebra for Kerr-Schild Kerr data the numerical instability disappears as demonstrated by the dashed (third) curve in Fig. 2. In comparison, we also perform a numerical run where we set $\kappa_{0}/\mathbf{K}$ as identically zero which should be a good approximation for large $r$ according to the decay rates given above. The associated numerical error yields the dotted (first) curve in Fig. 2.

These tests clearly show that terms of the form $\kappa_{0}/\mathbf{K}$ may be problematic numerically if they are not treated with care. We have found, however, that this problem disappears when we pick the method corresponding to the third curve in Fig. 2. This is what we will do in all of what follows. As a matter of fact, there are indeed several standard techniques in the literature for dealing with such issues which we could have employed as well (for instance, we could have rescaled the unknowns with appropriate powers of $r$). Without going into the details, we refer the interested reader to [23].

We conclude this section with a quick remark about the CFL condition. Given that we use explicit integrators, in particular the Runge-Kutta-Fehlberg method, one could think that the CFL condition may be violated in our numerical runs, for example for large values of $r$. If this was the case, our conclusions about the decay of the unknown quantities in the limit $r\rightarrow\infty$ below could not be trusted. Observe that while for other codes, $r$ can range from small to large values at each time step, in particular when large spherical grids are used, it plays the role of the “time parameter” for us. In our case, the CFL condition therefore restricts the “time” step depending on the current single value of $r$. Since the Runge-Kutta-Fehlberg method adapts the time step dynamically during the evolution according to certain prescribed tolerances and truncation errors and the time step is therefore always compatible with the stability region of the integrator, the CFL condition is always satisfied – even for large values of $r$. For an extended discussion of related issues, see for instance [13].

4 Asymptotics of nonlinearly perturbed Schwarzschild initial data

4.1 Spherically symmetric data

Now we want to present the main results of this paper in detail. Let us start by recalling the family of exact solutions of the constraint equations given by Eqs. (26) and (31) for any $C\geq 0$. The exact Kerr-Schild Schwarzschild initial data set corresponds to $C=0$. The first main question we address is whether all members of this family of spherically symmetric initial data sets are asymptotically flat.

The metric $h_{ij}$ given by the free data Eq. (26) takes the form

[TABLE]

When we introduce the new radial coordinate

[TABLE]

it takes the form

[TABLE]

The metric is therefore consistent with both the weak and the strong notions of asymptotic flatness when $R$ is identified with $\varrho$ (and when Cartesian coordinates are introduced in the standard way from polar coordinates $(R,\theta,\varphi)$), cf. Eqs. (4), (2) and (5).

Regarding the second fundamental form, we find

[TABLE]

which follows directly from Eq. (8). A necessary condition for the validity of the limit Eq. (3) (or (6)) is therefore that $\chi=O(R^{-2})$ (or $\chi=o(R^{-3/2})$) in the limit $R\rightarrow\infty$, which is equivalent to $\chi=O(r^{-2})$ (or $\chi=o(r^{-3/2})$) in the limit $r\rightarrow\infty$. Given that

[TABLE]

as a consequence of Eq. (31) it follows that

[TABLE]

using that $\kappa$ can be calculated from Eq. (14) and the quantities before. In conclusion, if $C>0$, then both limits Eqs. (3) and (6) are violated. The initial data set is therefore not asymptotically flat, neither in the strong nor in the weak sense. Therefore, only in the case $C=0$, the initial data set is asymptotically flat (both in the strong and the weak sense). In the latter case, one can that show that also all derivatives of these fields decay with appropriate rates. In particular all curvature quantities therefore decay at infinity. This is the case only for $C=0$.

Owing to the fact that the two cases $C=0$ and $C>0$ are therefore significantly different we shall spend the remainder of this subsection on the following geometric interpretation of these initial data sets. We first observe that since these initial data sets are all spherically symmetric, it must be possible to (locally) embed them isometrically into the Schwarzschild or Minkowski spacetime as a consequence of Birkhoff’s theorem. In fact, we can construct these isometric embeddings explicitly. They are depicted in Fig. 3 and can be described as follows:

Case $C=0$:

Let $(\Sigma,h_{ij})$ be the $3$-dimensional Riemannian manifold where the metric is given by Eq. (32) for any $M\geq 0$. As before, we use coordinates $(r,\theta,\varphi)$ on $\Sigma$. The two possible families of isometric embeddings of $(\Sigma,h_{ij})$ into the Schwarzschild spacetime with mass $m=M$ are

[TABLE]

for any $t_{0}\in\mathbb{R}$ and $r>2M$, where $(t_{\mathrm{SS}},r_{\mathrm{SS}},\theta_{\mathrm{SS}},\varphi_{\mathrm{SS}})$ are standard Schwarzschild coordinates. The pull-back of the second fundamental form induced from the Schwarzschild metric with mass $m$ on any such slice agrees with the field $\chi_{ij}$ determined as in Section 2.1 from the fields given by Eqs. (26), (27) and (31) for $C=0$, $M=m$ and the sign $\mp$ in Eq. (31). Notice that Eq. (37) is the well-known formula for (in- and outgoing) Kerr-Schild constant-time slices with respect to standard Schwarzschild coordinates.

Case $C>0$:

Let $(\Sigma,h_{ij})$ be the $3$-dimensional Riemannian manifold where the metric is given by Eq. (32) for any $M\in\mathbb{R}$. As before, we use coordinates $(r,\theta,\varphi)$ on $\Sigma$. The two possible families of isometric embeddings of $(\Sigma,h_{ij})$ into any Schwarzschild spacetime with mass $m>\max\{0,M\}$ (i.e., $m\not=M$) are

[TABLE]

for any $t_{0}\in\mathbb{R}$ and $r>2\max\{m,-M\}$, where, as above, $(t_{\mathrm{SS}},r_{\mathrm{SS}},\theta_{\mathrm{SS}},\varphi_{\mathrm{SS}})$ are standard Schwarzschild coordinates. The pull-back of the second fundamental form induced from the Schwarzschild metric with mass $m$ on any such slice agrees with the field $\chi_{ij}$ determined as in Section 2.1 from the fields given by Eqs. (26), (27) and (31) for

[TABLE]

and the sign $\mp$ in Eq. (31). Observe that this finding is in full consistency with a result by Ó Murchadha and Roszkowski in [29] (even though the authors do not provide the explicit embedding formulas).

Loosely speaking, the case $C>0$ yields initial data sets with the “wrong mass” $m\not=M$; it is therefore not surprising that such initial data sets are not asymptotically flat. We have found that within the family of spherically symmetric slices in the Schwarzschild spacetime, this subfamily of slices with the “wrong mass” is generic. The main practical consequence for the initial value problem of the hyperbolic formulation of the constraint equations in Section 2 is that the initial data must therefore be fine-tuned in order to obtain an asymptotically flat initial data set (i.e., one with the “right mass” $m=M$). While in the spherically symmetric case, we know explicitly how to do this, the problem of finding asymptotically flat initial data sets in less symmetric situations using the approach in Section 2 is therefore likely to be very hard. In order to obtain further insights we shall now study perturbations of the above data sets in the axially symmetric case.

4.2 Axially symmetric data

In this subsection we keep the “exact Kerr-Schild Schwarzschild free data” in Eq. (26), but then allow in principle arbitrary axially symmetric solutions $\mathbf{K}$ and $\mathbf{k}_{1}$ of the fully nonlinear system Eq. (29) (as long as the hyperbolicity condition Eq. (17) holds). The resulting initial data sets are therefore axially symmetric nonlinear perturbations of the Kerr-Schild Schwarzschild initial data set in the sense discussed in Section 2.3. The objective is then to study the asymptotics of these data sets in the limit $r\rightarrow\infty$ in the light of our results of the spherically symmetric case in Section 4.1. In a first step we restrict to heuristic and approximate studies in order to develop a first feeling of what might happen. In Section 4.3 we shall then tackle the fully nonlinear setting numerically.

To this end, we suppose that we have a one-parameter family of axially symmetric solutions of Eq. (29) of the form

[TABLE]

where $\mathbf{K}^{(0)}$ only depends on $r$ and where $\epsilon$ is a small parameter. When we plug this ansatz into Eq. (29) and evaluate it for $\epsilon=0$, the equations become that of the spherically symmetric case Eq. (30), i.e.,

[TABLE]

for $\kappa_{0}$ given by Eq. (27). As before, the general solution, which is real for all large values of $r$, is

[TABLE]

for any $C\geq 0$.

Next we average the first equation in Eq. (29) over any surface $S_{r}$ with $r>0$, i.e., we integrate with respect to $\theta$ (and $\varphi$, but recall that we restrict to the axially symmetric case) using the standard volume element of the round unit $2$-sphere and then normalize this by dividing by $4\pi$. When we now (i) use the fact that such an average of the term $\eth\mathbf{k}_{2}+\bar{\eth}\mathbf{k}_{1}$ in this equation is zero (because $\mathbf{k}_{1}$ and $\mathbf{k}_{2}$ have spin-weight $1$ and $-1$, respectively), (ii) consider the term $4|\mathbf{k}_{1}|^{2}/(r^{3}\mathbf{K})$ in the equation as negligible (because it is $O(\epsilon^{2})$), and, (iii) observe that

[TABLE]

we find that we can identify $\mathbf{K}^{(0)}(r)$ with the average

[TABLE]

in leading-order in $\epsilon$. The average of $\mathbf{K}^{(1)}$ in Eq. (39) must therefore be identically zero at every $r$. Because of this we shall now always write $\underline{\mathbf{K}}$ instead of $\mathbf{K}^{(0)}$ in all of what follows. From Eq. (40), it thus follows that

[TABLE]

in the limit $r\rightarrow\infty$.

Observe that in the language of spin-weighted spherical harmonics, see Section 2.2, $\underline{\mathbf{K}}=\mathbf{K}^{(0)}$ is proportional to the $l=0$-mode of $\mathbf{K}$, while $\mathbf{K}^{(1)}$ is given by the collective contributions of all $l>0$-modes and hence describes the variations of $\mathbf{K}$ along the spheres at each value of $r$.

Now let us derive certain bounds for $\mathbf{K}^{(1)}$ and $\mathbf{k}_{1}$ according to the ansatz Eq. (39) which shall consequently allow us to deduce the asymptotics of these quantities. To this end, we assume that $\epsilon$ is so small such that

[TABLE]

holds for every $r$, $\theta$ and $\varphi$. If this is true it implies that $|\mathbf{K}|<2|\underline{\mathbf{K}}|$ everywhere. We can thus conclude that $\mathbf{K}$ decays (at least) as fast as $\underline{\mathbf{K}}$ in Eq. (42). If we also assume that the hyperbolicity condition (17) holds everywhere, we can deduce that

[TABLE]

where we used Eq. (14) and previous estimates. Since $\kappa_{0}={O}(r^{-4})$ (see Eq. (27)), Eq. (42) implies

[TABLE]

Notice that this may not at all be a sharp decay estimate. In fact, the better the hyperbolicity condition is satisfied, i.e., the more negative the function \text{\boldmath\kappa}\ \mathbf{K} is, the less optimal we expect Eq. (43) to be.

In any case, given all these estimates, Eq. (14) now implies that

[TABLE]

and Eq. (34) yields

[TABLE]

Since the free data, and therefore the metric $h_{ij}$, are the same as in the spherically symmetric case in Section 4.1, the same analysis regarding asymptotic flatness in Section 4.1 can be performed here. It suggests that generic resulting initial data sets (i.e., those corresponding to $C>0$) are not asymptotically flat in the strong sense, see in particular Eq. (3). Regarding the weak notion of asymptotic flatness, see in particular Eq. (6), the situation is more subtle. In fact, if the decay exponent of $\chi$ is $3/2+\epsilon$ for an arbitrarily small $\epsilon>0$, the initial data set may be asymptotically flat in this weak sense. Arguably it appears to be unlikely that the second fundamental form of perturbed data sets, which we consider here, would decay better than the second fundamental form of generic unperturbed (spherically symmetric) data sets given by $C>0$. In any case, the next subsection is devoted to the fully nonlinear study of such issues.

Before we do this however let us close this subsection with the following remark regarding the case $C=0$ for which our analysis above suggests asymptotic flatness (both in the weak and the strong sense). Recall that our analysis here was based on neglecting terms of higher order in $\epsilon$ and therefore on neglecting nonlinear couplings between the $l=0$-mode and all modes $l>0$. As a consequence of this the average $\underline{\mathbf{K}}$ was a solution of the “spherically symmetric equations” which have the exact solution Eq. (40). Nonlinear interactions of the modes may however produce significant deviations in particular in the limit $r\rightarrow\infty$. It is therefore conceivable that, even if we pick initial data for the hyperbolic constraints which correspond to $C=0$, the asymptotics will nevertheless be described by the generic case $C>0$.

4.3 Numerical analysis

The purpose of this section is to support the largely heuristic arguments presented in the previous subsection by numerically solving the fully nonlinear system of constraints. This will allow us to confirm the decay rates derived before and thereby to provide support for our claims regarding asymptotic flatness of the resulting initial data sets. To this end we choose the same free data in Eq. (26) as before, but now we solve the fully nonlinear system Eq. (29) numerically. As the initial hypersurface we pick $S_{r_{0}}$ with $r_{0}=10^{3}$ and as initial data we pick the family

[TABLE]

for parameters $\beta\geq 0$ and $\epsilon\in\mathbb{R}$ which implies that $\underline{\mathbf{K}}|_{r=r_{0}}$ is given by Eq. (40) with $r=r_{0}$ and

[TABLE]

In the following we always make sure that $\epsilon$ is sufficiently small in comparison to $\beta$ so that $\mathbf{K}|_{r=r_{0}}$ never vanishes anywhere. Observe that we use the “axial symmetric notation” here introduced in [8], that is, ${Y_{l}}(\theta)={{}_{0}Y_{l0}}(\theta,\varphi)$ and ${{}_{1}Y_{l}}(\theta)={{}_{1}Y_{l0}}(\theta,\varphi)$. Moreover, we write

[TABLE]

In all the numerical evolutions we choose $N=16$ as the number of $\theta$-grid points. The inequality (17) is always satisfied during all numerical evolutions. Hence, the system (29) always remains hyperbolic. In the following plots we show the behavior of the modes $a_{l}(r)$ and $b_{l}(r)$ up to order four. In all these figures we have picked constants $\alpha$ and $\tilde{\alpha}$ in order to place certain reference curves optimally. Furthermore, for a better appreciation of the different orders of magnitude, we have used a logarithmic scale for the horizontal axis and a “symmetric logarithmic scale” for the vertical axis (the “symlog” option in matplotlib in python, for details see for instance [38]).

Case 1 ($\beta=10^{-5}$ and $\epsilon=10^{-6}$).

According to Eq. (46) this corresponds to the case $C>0$. From the discussion in Section 4.2, we expect that $\mathbf{K}={O}(r^{-3/2})$ and $\mathbf{k}_{1}={O}(r^{-1/2})$. This is confirmed in Figs. 6 and 6. In particular, we see that $a_{0}$ (recall Eq. (47)) definitely decays slower than $r^{-2}$ and hence (as one can check) the strong notion of asymptotic flatness is violated. There is no evidence that it decays faster than $r^{-3/2}$ (which would be necessary for the weak notion). In fact, the dotted curve in Fig. 6, which shows $r^{3/2}a_{0}$, suggests that $a_{0}$ decays exactly like $r^{-3/2}$, therefore violating the weak notion of asymptotic flatness.

Case 2 ($\beta=0$ and $\epsilon=10^{-6}$).

This is an example of the (presumably) unstable case corresponding to $C=0$. It is very interesting that the numerical results in Figs. 8 and 8 seem to be largely compatible with the decay rates in Section 4.2 despite our argument in the last paragraph there. The solid curve in Fig. 6, however, clearly shows that while the decay rate of $a_{0}$ is faster for some period of the evolution (in fact, there, the decay is compatible with the predictions in Section 4.2), eventually, the decay rate approaches $-3/2$ as expected. The strong notion of asymptotic flatness is therefore definitely violated; most likely also the weak notion.

Case 3 ($\beta=10^{-1}$ and $\epsilon=10^{-2}$).

Now we choose a relatively large value of $\epsilon$. The heuristic analysis in Section 4.2, for which we assume that $\epsilon$ is small, may therefore not yield an accurate description. Our numerical results are shown in Figs. 10 and 10. It is evident in these plots that $a_{0}$ is certainly not compatible with the decay rate $r^{-2}$ associated with the strong notion of asymptotic flatness. The dashed curve in Fig. 6 strongly suggests that as Case 1 before, $a_{0}$ decays exactly like $r^{-3/2}$ for large values of $r$ and therefore probably also violates the weak notion of asymptotic flatness.

5 Conclusions

We have studied the recently proposed hyperbolic/algebraic formulation of the vacuum constraint equations in General Relativity using both heuristic and numerical tools. In particular, we studied the asymptotic behavior of solutions of this system which can be viewed as (axisymmetric) perturbations of the Kerr-Schild data for the Schwarzschild space-time. We found that these solutions are generically not asymptotically flat in the strong sense of Eqns. (2)-(3) and we also see indications that they violate the weak asymptotic flatness conditions of Eqns. (5-6).

It is fair to say that up to now no method to construct initial data is known which yields full control of the actual physics of the initial data set in a clean and direct manner. People have developed a certain level of intuition and experience for more traditional approaches over many years. But even there, many issues are unsolved (for example, regarding “junk radiation”). For this relatively new framework considered here, we are only at the beginning of developing this experience. The goal of our work here was therefore clearly to construct “simple”, rather than “physically realistic” initial data sets. Moreover, we address only one of many crucial aspects which characterize the physical properties of initial data sets, namely asymptotic flatness. The fact that there is no obvious way to guarantee this property for the approach considered here could be seen as a major drawback. Building on these first steps, a next natural idea could be to extend our analysis to some other, possible more physically relevant scenario, like Schwarzschild initial data perturbed by a Teukolsky wave.

One problem with our present approach is that we start from an arbitrary sphere; there are no geometrically distinguished conditions that one could impose on the data. This might be different if one could start the solution process at a single point (a sphere with zero radius) where the conditions of regularity might be enough to fix a unique solution. Then the initial data set would be determined essentially by the free data. Another possibility would be to start the solution at infinity and integrate inwards. The conditions of asymptotic flatness might be enough to guarantee a unique solution. The big problem in both cases is, however, that one cannot enforce things at the other end: when evolving outwards the solution may violate asymptotic flatness and when evolving inwards we may end up with a singular initial data set.

It would be interesting to see whether there are conditions on the free data that could guarantee the correct behavior at the end of the evolution. It would also be interesting to see how things change when one uses the hyper/parabolic formulation of the constraints instead of the hyperbolic/algebraic formulations that we used in the present article.

Acknowledgments

L.E. was partly funded by the University of Otago Research Grant “Dynamical dark energy in the young universe and its consequences for the present and future history” in 2016.