Numerical construction of initial data sets of binary black hole type using a parabolic-hyperbolic formulation of the vacuum constraint equations

TL;DR

This paper explores a parabolic-hyperbolic approach to construct initial data for binary black holes, revealing solutions typically have cone geometry, and proposes an iterative scheme to approximate asymptotically Euclidean data.

Contribution

It introduces a novel numerical method for constructing binary black hole initial data using a parabolic-hyperbolic formulation and analyzes the asymptotic properties of these solutions.

Findings

Solutions generally exhibit cone geometry rather than asymptotic flatness

An iterative numerical scheme can approximate asymptotically Euclidean data

The approach successfully models binary black hole initial data without spin

Abstract

In this paper we investigate the parabolic-hyperbolic formulation of the vacuum constraint equations introduced by R{\'a}cz with a view to constructing multiple black hole initial data sets without spin. In order to respect the natural properties of this configuration, we foliate the spatial domain with 2-spheres. It is then a consequence of these equations that they must be solved as an initial value problem evolving outwards towards spacelike infinity. Choosing the free data and the "strong field boundary conditions" for these equations in a way which mimics asymptotically flat and asymptotically spherical binary black hole initial data sets, our focus in this paper is on the analysis of the asymptotics of the solutions. In agreement with our earlier results, our combination of analytical and numerical tools reveals that these solutions are in general not asymptotically flat, but have a cone geometry instead. In order to remedy this and approximate asymptotically Euclidean data sets, we then propose and test an iterative numerical scheme.

Figures

Equations

$$R-K_{ab}{K}^{ab}+K^2=0,{\nabla }_a{K^a}_c-{\nabla }_{c}K=0$$

$${\gamma }_{ij}$$

$$K_{ij}$$

$$R={\left(\sum _{a=1}^3({\tilde{x}}^i{)}^2\right)}^{1/2}\to \infty ,$$

$$N_a=A{\nabla }_a\rho$$

$$h_{ab}={\gamma }_{ab}-N_a{N}_b,$$

$$k_{ab}=-\frac{1}{2}{\mathcal{L}}_N{h}_{ab},$$

$${h^a}_b={{\delta }^a}_b-N^a{N}_b.$$

$$K_{ab}=\kappa N_a{N}_b+N_a{p}_b+N_b{p}_a+q_{ab},$$

$$q_{ab}=Q_{ab}+\frac{1}{2}q{h}_{ab},Q_{ab}{h}^{ab}=0.$$

$${\rho }^a{\nabla }_a\rho =1.$$

$${\rho }^a=A{N}^a+B^a.$$

$$k_{ab}=-A^{-1}\left(\frac{1}{2}{\mathcal{L}}_{\rho }{h}_{ab}-D_{(a}{B}_{b)}\right)=:A^{-1}{\stackrel{\star }{k}}_{ab}.$$

$$\stackrel{\star }{k}=h^{ab}{\stackrel{\star }{k}}_{ab}.$$

$$v_a:=N^b{\nabla }_b{N}_a=-A^{-1}{D}_{a}A.$$

$$\stackrel{\star }{k}{\nabla }_{\rho }A+A^2{D}^a{D}_{a}A-\stackrel{\star }{k}{B}^a{D}_{a}A=$$

$${\nabla }_{\rho }q-B^a{D}_{a}q-A{D}_a{p}^a-2p^a{D}_{a}A=$$

$${h^b}_c{\nabla }_{\rho }{p}_b-B^a{D}_a{p}_c-\frac{1}{2}A{D}_{c}q=$$

$$+\stackrel{\star }{k}{p}_c+p^b{\stackrel{\star }{k}}_{bc}+\kappa D_{c}A,$$

$$E={}^{(2)}R+2\kappa q-2p^a{p}_a-Q_{ab}{Q}^{ab}+\frac{1}{2}{q}^2,F=2({\nabla }_{\rho }\stackrel{\star }{k}-B^a{D}_a\stackrel{\star }{k})-{\stackrel{\star }{k}}_{ab}\stackrel{\star }{k}{}^{ab}-\stackrel{\star }{k}{}^2.$$

$$\begin{array}{rl}{\mathcal{L}}_{\rho }{p}_b-B^a{D}_a{p}_c-\frac{1}{2}A{D}_{c}q=& p_a{D}_b{B}^a-A{D}_a{Q^a}_c+\kappa D_{c}A-{Q^a}_c{D}_{a}A-\frac{1}{2}q{D}_{c}A\\ & +\stackrel{\star }{k}{p}_c+A{D}_c\kappa .\end{array}$$

$$\stackrel{\star }{k}<0.$$

$$g_{\alpha \beta }={\eta }_{\alpha \beta }-V{l}_{\alpha }{l}_{\beta }.$$

$${\gamma }_{ab}$$

$$\alpha$$

$$K_{ab}$$

$${\tilde{l}}^a=({\delta }^{-1}{)}^{ab}{l}_b,$$

$${\tilde{l}}^a{l}_a=({\delta }^{-1}{)}^{ab}{l}_a{l}_b=1.$$

$${\gamma }^{ab}=({\delta }^{-1}{)}^{ab}+\frac{V}{1-V}{\tilde{l}}^a{\tilde{l}}^b$$

$$l^a=\frac{1}{1-V}{\tilde{l}}^a,l^a{l}_a=\frac{1}{1-V}.$$

$$l_a=\pm f{\nabla }_a\rho ,$$

$$f=\frac{1}{\sqrt{({\delta }^{-1}{)}^{ab}{\nabla }_a\rho {\nabla }_b\rho }}.$$

$$N_a=\sqrt{1-V}{l}_a,$$

Full text

1 Introduction

Thanks to the ground-breaking work by Choquet-Bruhat and Geroch [23, 16] we know that each solution of the Einstein vacuum constraint equations (an initial data set, see below) determines a unique solution of the full vacuum Einstein’s equations – the so-called maximal globally hyperbolic development – in which the initial data set arises as the induced geometry of some spacelike hypersurface. A major concern of modern general relativity research is therefore how to construct solutions of these constraints. In all of what follows, a triple $(\Sigma,\gamma_{ab},K_{ab})$ of a $3$-dimensional differentiable manifold $\Sigma$, a Riemannian metric $\gamma_{ab}$ and a smooth symmetric tensor field $K_{ab}$ on $\Sigma$ is called an initial data set if it satisfies the vacuum constraint equations1.11.11.1All physical units in this paper are geometric ($G=c=1$), i.e., distances, time and mass are measured with the same physical unit.

[TABLE]

everywhere on $\Sigma$. Here, $\nabla_{a}$ is the Levi-Civita covariant derivative determined by $\gamma_{ab}$, $R$ the Ricci scalar and $K={K^{a}}_{a}$. For this whole paper we agree that spatial abstract indices $a,b,\ldots$ are raised and lowered with the metric $\gamma_{ab}$.

In general, the physical scenario of interest determines which particular properties one is looking for in solving Eq. (1.1). In this paper here for example we are interested in isolated systems of vacuum black holes. However, the constraints form an under-determined system and, to the best of our knowledge, it is not currently known if there exists a geometrically or physically preferable way to decide what part of the data should be specified and what part should be solved for. One of the most widely used approaches, both analytically and numerically, is that by Lichnerowicz and York (see [5, 6] and references therein) which involves specifying the conformally covariant ingredients of the data and solving a non-linear elliptic boundary value problem equivalent to the constraints. There are several well-established numerical techniques [6] for solving this boundary value problem.

While this approach by Lichnerowicz and York is very successful both analytically and numerically, it has also been criticised on several grounds in the literature. There are analytical problems when solutions are sought with highly non-constant mean curvature (see [21, 2] for an overview and references). There are physical problems related to “spurious radiation” (see e.g. [18], but also [25]). While discussions of such issues are ongoing in the literature, they have motivated several authors [15, 26, 15, 28, 14] to strive for alternative approaches for solving the constraints. The particular approach of interest for this paper here was introduced by Rácz in a series of papers [32, 33, 36, 34] following earlier work by Bartnik [4]. As we discuss below in more detail, the constraint equations are turned into either a hyperbolic-algebraic system or a parabolic-hyperbolic system depending on which components of the data are regarded as free. These systems are to be solved as an initial value problem.

First steps in investigating whether this approach has any advantages over more established methods have been carried out in [35] for the constraints of the Maxwell equations and in [30, 22] for the Einstein constraint equations. One obvious issue with solving the constraints as an initial value problem is the lack of control of the properties of the solutions either in the strong field regime or at spacelike infinity – two regimes which crucially characterise the physical properties of the resulting data set. Winicour considered the hyperbolic-algebraic formulation to obtain linearised perturbations of the Minkowski spacetime [40]. Somewhat disappointingly, he found that there are no suitable Cauchy hypersurfaces in the Minkowski spacetime on which the linearised algebraic-hyperbolic formulation of the constraints are well posed. This paper here now is a continuation of the investigations, which we started in [11] for the hyperbolic-algebraic formulation, of the asymptotics at spacelike infinity for the parabolic-hyperbolic one. In order to understand which questions are of interest to us, let us note that any initial data set should have some reasonable asymptotic behaviour at spacelike infinity if it was to represent black holes. Without this fundamental notions associated with black holes, like event horizons and masses, and with gravitational radiation may not even be defined (see, for example, [38]). The most prominent requirement, which we shall consider in this paper, is asymptotic flatness.

Definition 1

According to [20] we say that an initial data set $(\Sigma,\gamma_{ab},K_{ab})$ is 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 $\gamma_{ab}$ and $K_{ab}$ with respect to these coordinates satisfy1.21.21.2 We use the $O$ symbol rather informally in the usual sense $f=O(g)\iff|f|\leq C|g|$ for some constant $C>0$ in the relevant limit. In the case of tensors, we require that each coordinate component of the tensor satisfies such an estimate with respect to some natural coordinate system (like Cartesian coordinates in Definition 1). In this paper, we avoid the technicalities of choosing appropriate norms for the definition of the $O$-symbol. 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 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 decays at infinity. The interested reader can find details in the references [38, 20].

[TABLE]

in the limit

[TABLE]

Here $\delta_{ab}$ denotes the Euclidean metric and $M\geq 0$ is the ADM mass.

If Eq. (1.3) fails, but the other conditions are satisfied, the initial data is called asymptotically Euclidean. Regarding other notions of asymptotic flatness, which are useful in different contexts, see [13, 17].

In [11] the asymptotics of a class of initial data sets constructed by means of Rácz’s algebraic-hyperbolic formulation of the constraints was investigated for perturbations of a single Schwarzschild black hole. It was found that, in general, these initial data sets violate asymptotic flatness. In this paper now, we continue these explorations by investigating the parabolic-hyperbolic formulation in a multiple black hole setting. One of the key ingredients of Rácz’s approach is that space is foliated in a $2+1$-fashion. While the topology of the $2$-dimensional leaves of the foliation is largely free, a particularly natural choice to capture an isolated black hole configuration, which we adopt in all of what follows, are distorted $2$-spheres centred in the strong field regime around the black holes (in contrast to the choice of $2$-planes in [30, 22]). Because of this choice, we can also reuse the numerical pseudo-spectral methods developed in [12, 9, 10, 7] which were also employed in [11]. As we discuss in detail, it is then a consequence of the parabolic-hyperbolic equations by Rácz that they must be solved as an initial value problem starting from a boundary in the strong field regime evolving outwards towards spacelike infinity. Similar to the findings in [11] the resulting data sets turn out to violate asymptotically flatness; in fact, they are in general not even asymptotically Euclidean. In order to remedy this, we then introduce an iterative procedure which we demonstrate to allow us to approximate asymptotically Euclidean initial data sets. This paper exclusively focusses on the asymptotics at spacelike infinity. We have not performed any searches for apparent horizons or performed any other investigations of the properties of our data sets in the strong field regime.

Our paper is organised as follows. In Section 2.1 we briefly summarise Rácz’s parabolic-hyperbolic formulation of the constraint equations while Section 2.2 introduces the Kerr-Schild formalism which turns out to be useful for us later. In Section 3 we consider some explicit solutions of the parabolic-hyperbolic system and analyse a family of exact spherically symmetric solutions in order to set the scene for our analysis in the following sections. We present a brief discussion of Bishop’s superposition method [15] and then adapt it to the parabolic-hyperbolic formulation of the constraints. We discuss the numerical methods involved in solving the system of equations. In Section 4 we analyse the asymptotic behaviour of the numerical solutions and introduce an iterative procedure for constructing asymptotically Euclidean initial data sets. Finally, in Section 5, we summarise our main findings.

2 Preliminaries

2.1 The parabolic-hyperbolic formulation of the constraints

In this section we briefly summarise Rácz’s parabolic-hyperbolic formulation of the vacuum constraint equations. Further details can be found in [32, 33, 34, 36]. Since our conventions and notation partly deviate from those in these earlier references we provide a table in Appendix A comparing the different conventions.

Consider any data set $(\Sigma,\gamma_{ab},K_{ab})$ where $(\Sigma,\gamma_{ab})$ is a $3$-dimensional Riemannian manifold and $K_{ab}$ is a smooth symmetric tensor field on $\Sigma$. The Levi-Civita covariant derivative associated with $\gamma_{ab}$ is denoted by $\nabla_{a}$. Suppose in addition that there is a smooth function $\rho:\Sigma\rightarrow\mathbb{R}$ whose level sets $S_{\rho}$ are smooth $2$-dimensional surfaces in $\Sigma$ and that the collection $S$ of all these surfaces is a foliation of $\Sigma$. The quantity $t^{a}\nabla_{a}\rho$ vanishes for all vectors $t^{a}$ tangent to the surfaces $S_{\rho}$ and

[TABLE]

is the unit co-normal where $A$ is a strictly positive smooth function called the lapse (function).

Now we decompose the data set $(\Sigma,\gamma_{ab},K_{ab})$ with respect to the foliation defined by the function $\rho$ in full analogy to the standard $3+1$-decomposition of spacetimes, see for example [6].

The first and second fundamental forms of the surfaces $S_{\rho}$ are

[TABLE]

and2.12.12.1Observe our sign convention for second fundamental forms. As outlined in the table in Appendix A this choice of sign is different in most of the earlier literature.

[TABLE]

respectively. The orthogonal projector onto the surfaces is therefore

[TABLE]

The covariant derivative associated with $h_{ab}$ is referred to as $D_{a}$.

We say that a tensor field on $\Sigma$ is intrinsic (to the surfaces $S_{\rho}$) if any contraction with $N_{a}$ or $N^{a}$ vanishes. In particular, $h_{ab}$ and $k_{ab}$ are both intrinsic. Contracting all indices of a tensor field with ${h^{a}}_{b}$ always yields an intrinsic tensor field. In fact, any tensor field can be decomposed uniquely into intrinsic and normal parts as follows

[TABLE]

where we require that $p_{a}$ and $q_{ab}$ are intrinsic. The symmetric intrinsic tensor field $q_{ab}$ can then be decomposed into its trace and trace-free part (with respect to $h_{ab}$)

[TABLE]

Next we pick an arbitrary smooth vector field $\rho^{a}$ normalised by the condition

[TABLE]

Due to Eq. (2.1), this means that there must exist an intrinsic vector field $B^{a}$, the shift (vector), such that

[TABLE]

Given $\rho^{a}$, we can write $k_{ab}$ in Eq. (2.3) as

[TABLE]

We set

[TABLE]

As found in [34] and as we will discuss below, the crucial property of the quantity $\overset{\star}{k}_{ab}$ is that it only depends on $\rho^{a}$, $h_{ab}$ and $B_{a}$ (the “free data” for the constraint equations) and not on $A$ (one of the “unknowns”). Finally we use Eq. (2.1) to write the intrinsic acceleration as

[TABLE]

According to [34] the constraint equations (1.1) now take the following form by decomposing the momentum constraint covector into its normal and intrinsic parts and performing appropriate manipulations of the Hamiltonian constraint2.22.22.2These equations agree with [32], when the matter terms have been set to zero. :

[TABLE]

where $\nabla_{\rho}=\rho^{a}\nabla_{a}$, and,

[TABLE]

The Ricci scalar associated with $h_{ab}$ is called ${}^{(2)}R$. Notice that $\nabla$ can be eliminated in favour of Lie-derivatives. In Eqs. (2.11) and (2.12) we can write $\mathcal{L}_{\rho}A$ instead of $\nabla_{\rho}A$ and $\mathcal{L}_{\rho}q$ instead of $\nabla_{\rho}q$, and Eq. (2.13) can be written as

[TABLE]

These equations suggest to group the various fields introduced above as follows:

Free data

The fields $B_{a}$, $Q_{ab}$, $h_{ab}$ and $\kappa$ are considered as freely specifiable in (2.11)–(2.13) everywhere on $\Sigma$. Notice from the above that $\overset{\star}{k}$, $D_{a}$, ${}^{(2)}R$, $Q_{ab}$ and $F$ (and all index versions of these intrinsic fields) can be calculated from the free data everywhere on $\Sigma$.

Unknowns

The fields $A$, $q$ and $p_{a}$ are considered as unknowns which one attempts to determine as solutions of (2.11)–(2.13). Notice that all coefficients in these equations can be calculated from the free data everywhere on $\Sigma$.

According to [32], it can be shown that given any smooth initial data2.32.32.3The initial datum for $A$ is assumed to be strictly positive. for $A$, $q$ and $p_{a}$ on any $\rho=\rho_{0}$-leaf of the $2+1$-decomposition of $\Sigma$ in addition to smooth free data everywhere $\Sigma$, the initial value problem of Eqs. (2.11) – (2.13) in the increasing $\rho$-direction is well-posed, i.e., the equations have a unique smooth solution $A$, $q$ and $p_{a}$ at least in a neighbourhood of the initial leaf, provided the parabolicity condition2.42.42.4The table in Appendix A explains the sign discrepancy with [34]. holds everywhere on $\Sigma$:

[TABLE]

We remark that if $\overset{\star}{k}$ is positive instead, then the initial value problem in the decreasing $\rho$-direction is well-posed instead. In either case, Eqs. (2.11) – (2.13) is a quasilinear parabolic-hyperbolic system. It is important to notice that $\overset{\star}{k}$ is fully determined by the free data. Eq. (2.15) can therefore be checked before one attempts to solve Eqs. (2.11) – (2.13). We remark that in the context of black hole data sets, the initial data for (2.11)–(2.13) imposed on the initial $\rho=\rho_{0}$-leaf of the $2+1$-decomposition is often referred to as boundary conditions in the strong field regime in the literature.

As in [11] we restrict to the case $\Sigma=\mathbb{R}^{3}\backslash\overline{B}$ where $B$ is some finite ball in $\mathbb{R}^{3}$ in all of what follows. Moreover, we assume that the level sets of $\rho$ are diffeomorphic to $2$-spheres. Following [31, 8, 9, 10, 12, 11], all intrinsic tensor fields can therefore be written in terms of quantities with well-defined spin-weights (see Section B in the appendix for a quick summary). We can also express the intrinsic covariant derivative operator $D_{a}$ (defined with respect to the intrinsic metric $h_{ab}$) in terms of the covariant operator $\hat{D}_{a}$ defined with respect to the round unit-sphere metric $\Omega_{ab}$; recall that $D_{a}-\hat{D}_{a}$ can be expressed by some smooth intrinsic tensor field. Using Section B, we can then express the covariant derivative operator $\hat{D}_{a}$ in terms of the $\eth$- and $\eth^{\prime}$-operators [31]. Once all of this has been completed for all terms in Eqs. (2.11) – (2.13), each of these equation and each term ends up with a consistent well-defined spin-weight. Most importantly, however, all terms are explicitly regular: Standard polar coordinate issues at the poles of the $2$-sphere disappear when all quantities are expanded in terms of spin-weighted spherical harmonics and Eqs. (B.5) and (B.6) are used to calculate the intrinsic derivatives. From the numerical point of view this gives rise to a (pseudo-)spectral scheme. Further details related to our implementation can be found in [11].

2.2 Data sets of Kerr-Schild form

In this subsection we introduce general Kerr-Schild-like data sets. These will play an important role in the remainder of this paper. We stress that we do not yet impose the constraint equations in this subsection. Sometimes we call data sets unphysical or preliminary if the constraints are not imposed.

Let us start our discussion with a $4$-dimensional Lorentzian manifold2.52.52.5Note that we use Greek indices for denoting spacetime abstract indices. Spacetime indices are consistently raised and lowered with $g_{\alpha\beta}$. $(M,g_{\alpha\beta})$ where the metric $g_{\alpha\beta}$ is of Kerr-Schild form

[TABLE]

Here $\eta_{\alpha\beta}$ is the Minkowski metric, $l^{\alpha}$ is a null vector field with respect to $g_{\alpha\beta}$ and $V$ is a smooth spacetime function. Notice that $l_{\alpha}=g_{\alpha\beta}l^{\beta}=\eta_{\alpha\beta}l^{\beta}$ and that $l^{\alpha}$ is therefore automatically also null with respect to $\eta_{\alpha\beta}$. We also assume coordinates $(t,x,y,z)$ such that $\eta_{\alpha\beta}$ takes the form $\mathrm{diag}(-1,1,1,1)$ and normalise $l^{\alpha}$ such that its $t$-component is $1$.

Given this it turns out that the $t=\mathrm{const}$-surfaces $\Sigma_{t}$ are spacelike if $V<1$, as we shall always assume. The future directed (with respect to the $\partial_{t}$) unit normal is2.62.62.6Coordinate one-forms, e.g. $dt$, are either written in index-free notation, or, in abstract index notation $dt_{\alpha}$. Notice that $dt_{\alpha}=\overline{\nabla}_{\alpha}t$ where $\overline{\nabla}_{\alpha}$ is the covariant derivative associated with $g_{\alpha\beta}$. Similarly we write either $\partial_{t}$ or $\partial^{\alpha}_{t}$ for coordinate vector fields. $n_{\mu}=-dt_{\mu}/\sqrt{1-V}$. The spacetime lapse function is therefore $\alpha=1/\sqrt{1-V}$, the induced metric $\gamma_{\mu\nu}=g_{\mu\nu}+dt_{\mu}dt_{\nu}/(1-V)$, the shift is $\beta^{\mu}=\partial_{t}^{\mu}-\alpha n^{\mu}$ and the second fundamental form is $K_{\mu\nu}=-\frac{1}{2}\mathcal{L}_{n}\gamma_{\mu\nu}$. For any fixed $t\in\mathbb{R}$, consider now the embedding $\Phi:\Sigma\rightarrow\Sigma_{t}$, $(x,y,z)\mapsto(t,x,y,z)$ where $\Sigma$ is a smooth $3$-dimensional manifold diffeomorphic to $\Sigma_{t}$. Adopting the index notation introduced in Section 2.1 (in particular, all spatial index operations are performed with $\gamma_{ab}$ as before), we denote the pull-backs of the above quantities from $\Sigma_{t}$ to $\Sigma$ as $\gamma_{ab}$, $K_{ab}$, $l_{a}$, $\beta_{a}$ and $\alpha$, respectively. We find straightforwardly

[TABLE]

where $\delta_{ab}$ is the Euclidean metric which in Cartesian coordinates $(x,y,z)$ takes the form $\text{diag}(1,1,1)$. The field $l_{a}$ has magnitude $1$ with respect to $\delta_{ab}$, i.e., if we set2.72.72.7Observe carefully that $\delta^{ab}$ and $(\delta^{-1})^{ab}$ are different fields: The first one is defined by raising both indices of $\delta_{ab}$ with $\gamma^{ab}$, while the second one is the uniquely determined inverse of $\delta_{ab}$.

[TABLE]

then

[TABLE]

In particular, it follows that

[TABLE]

and

[TABLE]

Now pick a smooth function $\rho$ on $\Sigma$ with the properties discussed in Section 2.1 giving rise to a foliation $S$ in terms of level sets $S_{\rho}$. We restrict to the case where $l_{a}$ is normal to $S_{\rho}$, i.e.,

[TABLE]

with

[TABLE]

From Eqs. (2.1), (2.23) and (2.22) we find that

[TABLE]

which means that the lapse defined in Eq. (2.1) is

[TABLE]

Given Eqs. (2.2) and (2.17), it follows that

[TABLE]

Plugging all this into Eq. (2.19) yields

[TABLE]

where

[TABLE]

and $t^{\mu}=\partial_{t}^{\mu}$ has the coordinate representation $(1,0,0,0)$. From Eq. (2.4) we get

[TABLE]

where $v_{a}$ can be calculated from Eq. (2.10) and $k_{ab}$ from Eq. (2.8). The quantities $q$ and $Q_{ab}$ are given by Eq. (2.5).

Once we have picked intrinsic coordinate systems $(y^{1},y^{2})$ of the $\rho=\mathrm{const}$-surfaces and thereby an “adapted” coordinate system $(\rho,y^{1},y^{2})$ of $\Sigma$, the vector field $\rho^{a}=\partial_{\rho}^{a}$ is determined. The shift $B_{a}$ in Eq. (2.7) is then given as

[TABLE]

and $k_{ab}$, $\overset{\star}{k}_{ab}$ and $\overset{\star}{k}$ can be calculated from Eqs. (2.8) and (2.9).

A particularly important example is the Kerr-Schild Schwarzschild initial data set with mass $M\in\mathbb{R}$. It can be written in Kerr-Schild form as (see [1, 27])2.82.82.8In general one could have $l_{\mu}=-dt_{\mu}\pm dr_{\mu}$ which distinguishes ingoing Kerr-Schild coordinates from outgoing ones. We use the positive sign exclusively here.

[TABLE]

where $r=\sqrt{x^{2}+y^{2}+z^{2}}$. Using the formalism discussed above, we find that

[TABLE]

It is therefore consistent with the assumption Eq. (2.23) when we identify the function $\rho$ with the coordinate $r$, whose level sets are $2$-spheres. Straightforward calculations yield

[TABLE]

where $(\theta,\phi)$ are standard polar coordinates on each $r=\mathrm{const}$-sphere and where we defined

[TABLE]

Observe that the parabolicity condition Eq. (2.15) is satisfied for all $r>0$.

3 Multiple black hole initial data sets

3.1 Outline of our approach

In order to construct multiple black hole initial data sets using the formalism introduced in Section 2.1 we proceed now in two steps. The first step is to produce, without imposing the constraints yet, data sets which can somehow be interpreted as multiple black holes. As discussed in Section 3.3.1 we approach this in a largely ad hoc and certainly not unique way. Because the constraints are not imposed yet, such data sets are referred to as (preliminary) data sets. Only in the second step, see Section 3.3.2, such a preliminary data set is used to obtain the free data and the initial data for solving the constraints as the initial value problem of Eqs. (2.11) – (2.13). The solution is therefore a physical initial data set. Since we are particularly interested in the asymptotics of these, we start our discussion by analysing a family of spherically symmetric (single black hole) data sets, which we show to be limits at spatial infinity of more general initial data sets in Section 4 – in full analogy to the findings in [11].

3.2 Kerr-Schild-like spherically symmetric solutions of the constraints

In this section, we discuss a family of spherically symmetric solutions of the constraints encompassing the single Kerr-Schild Schwarzschild black hole data set introduced at the end of Section 2.2. This family is obtained by finding the general spherically symmetric solution of the parabolic-hyperbolic formulation of the constraints where only the free data are determined by Eq. (2.34). We shall discuss that while all these correspond to slices in a Schwarzschild spacetime with a particular mass (as a consequence of the Birkhoff theorem), almost all of these are not asymptotically flat.

In order to construct this family of spherically symmetric initial data sets now we pick any $M\in\mathbb{R}$ and choose the free data for Eqs. (2.11) – (2.13) by Eq. (2.34). Given these free data, we look for the general spherically symmetric solution of Eqs. (2.11) – (2.13). As a consequence of spherical symmetry we impose that $p_{a}=0$ and that $q$ and $A$ only depend on $r$. Under these assumptions, Eqs. (2.11) – (2.13) become

[TABLE]

Since the first equation is independent of the second one, we can easily express its general solution as

[TABLE]

where $\mathcal{C}$ is an integration constant, and we’ve defined $\lambda(r)$ as in Eq. (2.36). Substituting Eq. (3.3) into Eq. (3.2) leads to

[TABLE]

with $m\in\mathbb{R}$ being another integration constant. For later reference we note that

[TABLE]

The corresponding initial data set is then given by

[TABLE]

where $A$ is given by Eq. (3.4), $q$ by Eq. (3.3) and $h_{ab}$ and $\kappa$ by Eq. (2.34). We see that (3.4) and (3.3) agree with the particular solution (2.35) if and only if $\mathcal{C}=0$ and $m=M$, in which case, $M$ is the ADM mass.

In order to analyse the asymptotics of this family of data sets, we try to bring $\gamma_{ab}$ to the form Eq. (1.2) by introducing a new radial coordinate $R$. It turns out that this is not possible in general. It is always possible however to bring the metric asymptotically to the cone form

[TABLE]

with

[TABLE]

where

[TABLE]

From the asymptotic form (3.8) of the metric one can see that the area of the 2-spheres of constant $R$ is $4\pi R^{2}(1+\mathcal{C}^{2})\geq 4\pi R^{2}$. Therefore, $\mathcal{C}$ determines the excess area beyond the Euclidean value $4\pi R^{2}$. If this was an deficit area then this would be the metric of a cone embedded in a 4-dimensional Euclidean space. Nevertheless, we call $\mathcal{C}$ the cone parameter. The metric is asymptotically Euclidean, and hence is of the form Eq. (1.2), if and only if $\mathcal{C}=0$. Since

[TABLE]

the data set is asymptotically flat with ADM mass $m$ if and only if $\mathcal{C}=0$. Surprisingly we notice that the mass parameter $M$, which determines the free data via Eq. (2.34), has nothing to do with the actual mass of the data set.

The Hawking mass (see [1], Appendix A) of any $r=const$-sphere

[TABLE]

where $|{S}_{r}|$ is the surface area of the $r=\mathrm{const}$-sphere $S_{r}$, and, where

[TABLE]

are the in- and outgoing null expansion scalars defined with respect to suitably normalised future-pointing null normals of ${S}_{r}$, turns out to be

[TABLE]

where $m$ is the integration constant found in (3.4). So, even if $\mathcal{C}\not=0$ and the initial data set is therefore not asymptotically flat, we can still associate the mass $m$ with this data set.

Thus we find the result that all the spherically symmetric data sets, asymptotically flat or not, have a well-defined and finite Hawking mass limit. This is consistent with the following observation made in [11]: It is a consequence of the Birkhoff theorem that for any spherically symmetric data set $(\Sigma,\gamma_{ab},K_{ab})$ (which satisfies the constraints) there is a hypersurface in a Schwarzschild spacetime and a diffeomorphism from (a subset of) $\Sigma$ to (a subset of) that hypersurface such that the pull-backs of the first and second fundamental forms induced on that hypersurface by the Schwarzschild metric agree with $\gamma_{ab}$ and $K_{ab}$ (see also [29]). Indeed, this applies to the asymptotic region of all the spherically symmetric data sets above. Irrespective of the value of $\mathcal{C}$, the mass of that “target Schwarzschild spacetime” turns out to be $m$ as in Eq. (3.13). This embedding into the target Schwarzschild spacetime is described by the formula

[TABLE]

where $t$ is the time coordinate of the target Schwarzschild spacetime with mass $m$ given in Kerr-Schild form by Eqs. (2.16) and (2.33) (where $M$ must be replaced by $m$) and where the radial coordinate of $\Sigma$ is mapped to the target radial Schwarzschild Kerr-Schild coordinate (therefore bearing the same name). The quantity $t_{0}$ is an arbitrary constant.

3.3 Multiple black hole initial data sets

3.3.1 Step 1: Multiple black hole data sets

As outlined above, the idea of our initial data construction procedure is to first construct preliminary multiple black hole data sets in a first step without imposing the constraints. The idea is to make ‘natural’ choices for the fields $l_{a}$, $V$ and $\dot{\gamma}_{ab}$ above such that the data set obtained from the formulas in Section 2.2 resembles a multiple black hole system at a moment of time. Our particular approach for this is inspired by the work in [15].

To this end, we pick $n$ black hole mass parameters $M_{1},\ldots,M_{n}$ and $n$ Cartesian coordinate position vectors $(x_{1},y_{1},z_{1}),\ldots,(x_{n},y_{n},z_{n})$ and set

[TABLE]

with

[TABLE]

We set

[TABLE]

and notice that $\rho$ is a radial coordinate with

[TABLE]

in the limit $r\rightarrow\infty$.

From now on we restrict to the case $n=2$. We write $M_{+}=M_{1}$, $M_{-}=M_{2}$ and set

[TABLE]

for some fixed $Z\geq 0$ and $M_{+},M_{-}\geq 0$. This yields

[TABLE]

with

[TABLE]

In Fig. 1 we show the level-surface $\rho=\rho_{0}$ for two different values of $\rho_{0}$. Obviously, these surfaces undergo a topology change as $\rho_{0}$ varies. Determining the critical value shows that when3.13.13.1The existence of this lower bound suggests that our current approach to construct data sets is somewhat limited as we cannot put the initial $2$-surface arbitrarily close to the black holes. It is clear however that only small modifications of the geometry of these peanuts is necessary to overcome this problem.

[TABLE]

the $\rho=\rho_{0}$-surface is diffeomorphic to a single $2$-sphere; because of the particular shape of these $2$-spheres we sometimes refer to them as peanuts. The set of all $\rho=\rho_{0}$-surfaces satisfying Eq. (3.21) indeed yields a foliation of the exterior region of $\mathbb{R}^{3}$.

With $\rho$ as defined by Eqs. (3.17), (3.19) and (3.20), we define $l_{a}$ as in Eq. (2.23) with the positive sign and $f$ given by Eq. (2.24), i.e.,

[TABLE]

This means that $l^{a}$ is outward pointing. In order to determine a preliminary Kerr-Schild data set now using the formulas in Section 2.2, we pick

[TABLE]

and

[TABLE]

for an, at this point, arbitrary function $\mathcal{F}$. It is clear that the resulting data set equals the standard single Kerr-Schild Schwarzschild black hole initial data set discussed in Section 2.2 if $M_{+}=0$ (or $M_{-}=0$) and $\mathcal{F}=0$ (with a coordinate shift given by $Z$), or, if $Z=0$ and $\mathcal{F}=0$. In general notice that by the particular choice of $\dot{\gamma}_{ab}$ in Eq. (3.24), which we justify in Section 4.1 and Section 4.3, only the quantity $q$ is affected by the function $\mathcal{F}$, see Eqs. (2.30), (2.31) and (2.32).

For arbitrary fixed values of $M_{+},M_{-}$ and $Z$ we can calculate $\overset{\star}{k}$ as

[TABLE]

It follows that there exists a $\rho_{1}>0$ such that $\overset{\star}{k}$ is strictly negative (and hence the parabolicity condition Eq. (2.15) holds) for all $\rho>\rho_{1}$. The fact that we have not been able find an explicit formula for $\rho_{1}$ is not problem in practice because we can check numerically that $\overset{\star}{k}$ is negative on the whole numerical domain.

It is a consequence of Eq. (3.18) that

[TABLE]

and hence that each preliminary data set approaches the physical single Kerr-Schild Schwarzschild black hole data set with mass $M_{+}+M_{-}$ in leading order (provided $\mathcal{F}$ vanishes zero sufficiently fast in the limit $r\rightarrow\infty$).

From now on we shall decorate the quantities associated with the preliminary data sets with [P] in order to distinguish them from corresponding quantities of the final physical initial data set discussed in the next subsection.

3.3.2 Step 2: Solving the parabolic-hyperbolic constraints system

Adapted coordinates.

Given a data set as constructed in Section 3.3.1, we next attempt to impose the constraints. To this end, we can read off the free data and the initial data for solving Eqs. (2.11) – (2.13) as an initial value problem in the increasing $\rho$-direction starting from some $\rho=\rho_{0}$-surface compatible with (3.21). For this we introduce adapted coordinates $(\rho,\vartheta,\varphi)$ where $\rho$ is the radial coordinate along which we perform the evolutions and where each “peanut” $\rho=\text{const}$ is endowed with intrinsic polar coordinates $(\vartheta,\varphi)$. The coordinate transformation from the “original” spherical coordinates $(r,\theta,\phi)$ to these adapted coordinates $(\rho,\vartheta,\varphi)$ is taken to be of the form

[TABLE]

with $u(r,\theta,\phi)$ given by Eqs. (3.19) and (3.20) using

[TABLE]

It is clear that while Eqs. (2.11) – (2.13) are solved in the adapted $(\rho,\vartheta,\varphi)$-coordinate system, the coefficients of these equations must be calculated from the preliminary data set which has been found in the original $(r,\theta,\phi)$-coordinate system before. At each $\rho$-step of the numerical evolution, i.e., on any $\rho=\rho_{0}=\text{const}$-surface, we must therefore compute $r$ as a function of $\vartheta=\theta$. Since there is no explicit formula for this we find this function3.23.23.2We shall now mostly suppress the coordinate $\varphi=\phi$ because all examples considered in this paper are axisymmetric and therefore independent of this coordinate. $\hat{r}(\vartheta)=r(\rho,\vartheta)_{|\rho=\rho_{0}}$ numerically as follows3.33.33.3This is certainly only one way to determine this function numerically. We have not compared this to any other method (for example, the Newton method) yet.:

Calculate

[TABLE]

which follows by inverting the equation $\rho=\rho_{0}$ in the special case $\theta=\vartheta=0$. 2. 2.

Then numerically solve the ODE

[TABLE]

as an initial value problem on the interval $\vartheta\in[0,\pi]$ with the initial datum $\hat{r}(0)$ given by Eq. (3.27).

With the function $\hat{r}(\vartheta)$ determined in this way for each $\rho_{0}$ we proceed as follows. The tensor components of all quantities calculated in the $(r,\theta,\phi)$-coordinate system are transformed to the $(\rho,\vartheta,\varphi)$-coordinates using the Jacobi matrix of the coordinate transformation (3.26) which reads

[TABLE]

Observe that its inverse is

[TABLE]

With $\hat{r}(\vartheta)=r(\rho,\vartheta)_{|\rho=\rho_{0}}$ determined on any $\rho=\rho_{0}=\text{const}$-surface as discussed above, we can express all tensor components completely in terms of the adapted $(\rho,\vartheta,\varphi)$-coordinates as needed to numerically solve the constraint equations.

Numerical method, errors and tests.

At the end of Section 2.1, we have discussed that we can use the spin-weight formalism (see also Section B in the appendix3.43.43.4The parabolic hyperbolic equations written in terms of the SWSH can be found in [37]) to express all fields intrinsic to the $2$-spheres $\rho=\text{const}$, and thereby the constraint equations (2.11) – (2.13), using spin-weighted spherical harmonics; the details of our numerical implementation can be found in [8, 9, 10, 12, 11]. Thus, the “spatial” discretisation used in our code is of (pseudo-)spectral nature. All examples considered so far assume axial symmetry, so there is no dependence on $\varphi$. We can therefore employ uniform $\vartheta$-grids with $N$ points (see in particular [10] regarding details) and exploit axial symmetry. As the “time”-stepping method we choose the adaptive SciPy ODE solver odeint3.53.53.5See https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html.. We denote its absolute error tolerance parameter by $\tilde{\mathcal{E}}$ and the corresponding magnitude by $\varepsilon=-\log(\tilde{\mathcal{E}})$3.63.63.6For this whole paper, $\log$ is the logarithm to base $10$.. The parameter $\varepsilon$ therefore controls the local $\rho$-step size of the numerical evolutions.

We can expect that the parameters $N$ and $\varepsilon$ can be used to control the error that is numerically generated through the time- and space-discretisation. However, in certain settings we find that the error depends neither on $\varepsilon$ nor on $N$. Then we conclude that the error is dominated by the accumulated finite number representation errors in our code. In such a situation, our code provides only limited means to improve the numerical accuracy. We provide evidence for this in the following sections.

We now discuss how $\varepsilon$ and $N$ may be used to control the errors. Let $E[f](\varepsilon,N,\rho,\vartheta)$ be the absolute error of some particular unknown $f$ at $(\rho,\vartheta)$ calculated numerically with discretisation parameters $\varepsilon$ and $N$. In principle, this error can of course only be calculated if the exact solution is known. In practice, when the exact solution is not known, we shall follow the common practice to determine $E[f](\varepsilon,N,\rho,\vartheta)$ by comparing the numerical solution to another numerical solution obtained with some sufficiently high resolution (instead of the exact solution).

On the one hand, we expect that if $N$ is sufficiently large so that the grid resolves all the spatial features of the solution, the numerical error is dominated by the time discretisation. In such a setting, the numerical error should not become smaller when we increase $N$ (in fact, oversampling may be a significant error source). The error should decrease monotonically with $\varepsilon$. Unless stated otherwise we always pick $\varepsilon=12$. On the other hand, the error can be expected to be dominated by the spatial discretisation if $\varepsilon$ is sufficiently large. In this setting, the error should be roughly independent of $\varepsilon$, but should decrease monotonically with $N$. Unless stated otherwise we always take $N=11$.

In regards to the calculation of the error, we set for any fixed $N$ and $\varepsilon$

[TABLE]

which can be used to study the error in both the time and space discretisations.

We now discuss two tests of our code. First we consider a non-trivial test of our numerical implementation by choosing3.73.73.7In all examples in this paper we will choose $M_{+}+M_{-}=1$. Any quantity which carries either a distance, time or mass unit is therefore expressed in units of $M_{+}+M_{-}$ in the geometric physical unit system adopted in this paper. $M_{+}=Z=1$ and $M_{-}=\mathcal{F}=0$. Picking $Z>0$ ‘shifts’ the single black hole so that it is no longer centred at the origin of the $(r,\theta,\phi)$-coordinate system. The solutions obtained in the adapted coordinates $(\rho,\vartheta,\varphi)$ must however agree with the standard single black hole solution after undoing the shift. We pick $\varepsilon=8,10,12$ and calculate the quantity $\mathcal{E}_{\rho}$ using the exact single black hole solution as the reference solution. The results are shown in the left column of Fig. 2. The convergence of the functions $q$ and $A$ is consistent with that of our numerical scheme. The error associated with $p$ and $\bar{p}$ defined as3.83.83.8Notice that $p$ has spin-weight $1$ and its complex conjugate $\bar{p}$ the spin-weight $-1$.

[TABLE]

are below the numerical round-off error which occurs at $10^{-15}$.

The second test case is given by the choice $M_{+}=M_{-}=1/2$, $Z=1$ and $\mathcal{F}=0$. No exact solution is known here and the numerical solution obtained with some higher resolution is therefore chosen as the reference solution to calculate the error. The right most column of Fig. 2 shows convergence plots for the spatial discretisation. When a numerical solution is calculated with a spatial resolution $N_{(1)}$, the reference solution is calculated with a spatial resolution $N_{(2)}=2N_{(1)}-1$ to ensure that both numerical resolutions share grid points. The calculated errors shown in both the middle and right column of Fig. 2 have the expected dependence on $\varepsilon$ and $N$.

4 Asymptotic properties

4.1 Formal asymptotic expansions and their numerical justifications

In Section 3 we have proposed a construction procedure for a class of initial data sets and discussed some numerical test cases. The hope is that these initial data sets describe binary black hole systems. Whether this is really the case remains to be seen. In this section we make an important step towards the physical understanding of these data sets by analysing the asymptotics $\rho\rightarrow\infty$.

To this end first recall the results in the spherically symmetric case in Section 3.2. We found there that every spherically symmetric data set is uniquely determined by the parameters $\mathcal{C}$, $m$ and $M$, and is asymptotically flat with ADM mass $m$ if and only if the cone parameter $\mathcal{C}$ vanishes. In any case, the limit of the Hawking mass at $\rho\rightarrow\infty$ is always $m$ irrespective of the value of $\mathcal{C}$. Intuitively, we would expect that any initial data set obtained as in Section 3 should become spherically symmetric asymptotically. The asymptotics found in the spherically symmetric case should therefore apply to this larger family of initial data sets as well.

In order to provide evidence for this claim, we perform a formal power series analysis. To this end, we first calculate power series expansions (with respect to $\rho$ in the limit $\rho\rightarrow\infty$) of the general class of binary black hole preliminary data sets in Section 3.3.1. Recall that this yields the free data for the constraint equations (2.11) – (2.13). Second, we make the following power series ansatz for the unknowns $q$, $A$ and $p_{a}$ of Eqs. (2.11) – (2.13)4.14.14.1Consistently with footnote 1.2, the $O$-symbol is defined component-wise with respect to a natural coordinate system in the case of tensor quantities (here, the $(\rho,\vartheta,\varphi)$-coordinates).:

[TABLE]

The structure of the equations and the expressions of the free data suggest that it is sufficient to consider integer powers of $\rho$ only4.24.24.2This is significantly different for the results in [11] obtained with the hyperbolic formulation of the constraints where also half-integer powers of $\rho$ are present.. Shortly we will present numerical support for this ansatz. If these expansion do not fairly represent the asymptotic behaviour of the solutions then we would see evidence of this in our numerics.

Before we proceed we notice that all the following expansions are independent of the choice of $\mathcal{F}$ in Eq. (3.24). In fact it follows from the discussion in Section 3.3.1 and Eqs. (2.29) – (2.32) for the particular choice Eq. (3.24) that only the function $\prescript{[P]}{}{q}$, i.e., the function $q$ associated with the preliminary data set in Step 1 of our construction procedure in Section 3.3, is affected by the choice of $\mathcal{F}$. This consequently only affects the initial data of the resulting function $q$ in Step 2 of our construction procedure. While these initial data certainly affect the resulting values of the expansion coefficients in the following, the expansions themselves and the particular relationships between the various expansion coefficients, which we uncover now, hold irrespectively.

It turns out that the leading term of the expansion of Eq. (2.13), i.e., the term of order $\rho^{-1}$, yields the equation

[TABLE]

Recall that $\hat{D}_{a}$ is the covariant derivative associated with the metric $\Omega_{ab}$ of the standard round unit sphere. Motivated by the spherically symmetric case, see Eqs. (3.3) and (3.5), we write

[TABLE]

for some (possibly new) constant $\mathcal{C}$. Then, the leading term of the expansion of Eq. (2.11), i.e., the term of order $\rho^{-2}$, turns out to yield the nonlinear elliptic equation

[TABLE]

where $\hat{\Delta}=\Omega^{ab}\hat{D}_{a}\hat{D}_{b}$ is the Laplacian on the round unit sphere. First we notice that if $A^{(0)}$ is any smooth solution of this equation, then $A^{(0)}$ does not change sign because $A^{(0)}(\vartheta,\varphi)=0$ at any point $(\vartheta,\varphi)$ would clearly violate Eq. (4.5) at that point. We therefore conclude that any smooth solution of this equation is either strictly positive or strictly negative. Suppose we have a smooth strictly negative solution $A^{(0)}$. Then $-A^{(0)}$ is a smooth strictly positive solution. Assuming that $A^{(0)}$ is strictly positive from now on is therefore no loss of generality. A particular smooth strictly positive solution is the constant solution

[TABLE]

compare this to Eq. (3.5). In fact, Eq. (4.6) is the only smooth strictly positive (constant or not) solution of Eq. (4.5): Suppose there were two different smooth strictly positive solutions $A^{(0)}$ and $\tilde{A}^{(0)}$ of Eq. (4.5). Then a standard integration by parts argument implies

[TABLE]

where the norm and the scalar product here are the standard $L^{2}$-norm and $L^{2}$-scalar product on the $2$-sphere with respect to $\Omega_{ab}$. One can easily check that

[TABLE]

Since the fraction on the right-hand side is strictly positive if $A^{(0)}$ and $\tilde{A}^{(0)}$ are strictly positive, the right-hand side of Eq. (4.7) is therefore non-negative. But since the left-hand side is non-positive, $A^{(0)}$ and $\tilde{A}^{(0)}$ must be equal. Thus, we conclude that $A^{(0)}$ given by Eq. (4.6) is the only strictly positive smooth solution of Eq. (4.5). Due to Eq. (4.4) we therefore have

[TABLE]

in agreement with the spherically symmetric case Eq. (3.3). Due to Eqs. (4.6), (4.1) and the leading behaviour of $h_{ab}$ we can now interpret $\mathcal{C}$ again as a cone parameter just as in the spherically symmetric case. In particular it follows that the resulting initial data set is not asymptotically flat if $\mathcal{C}\not=0$.

Again in agreement with the spherically symmetric case, the Hawking mass Eq. (3.11) always has a well-defined limit $m$ as a function of $\rho$ if Eqs. (4.1) – (4.3) hold irrespective of whether $\mathcal{C}$ is zero or not. In fact,

[TABLE]

where $\underline{A^{(1)}}$ is the average of the function $\underline{A^{(1)}}$ over the round unit sphere4.34.34.3According to Eq. (B.9) and Eq. (B.11), we have $\underline{f}=f_{0}/\sqrt{4\pi}$ if $f$ is any axially symmetric smooth function on $\mathbb{S}^{2}$., see Eq. (B.8), and where we have used that both the $\rho^{-3}$-order of Eq. (2.11) and the $\rho^{-3}$-order of Eq. (2.12) imply that

[TABLE]

We have therefore found that if Eqs. (4.1) – (4.3) hold then the general solution of the constraint equations (2.11) – (2.13) with free data determined by the preliminary data sets in Section 3.3.1 is spherically symmetric in leading order and the limit of the Hawking mass is consistent with the spherically symmetric case.

In general, however, the solutions are not spherically symmetric beyond the leading order. We have seen that the initial data set is not asymptotically flat and the $3$-metric is not asymptotically Euclidean if $\mathcal{C}\not=0$. In the case $\mathcal{C}=0$, the next order of the expansion of the constraint equations takes the form

[TABLE]

This implies that

[TABLE]

and the general axially symmetric solution for $q^{(2)}$ is

[TABLE]

where $q^{(2)}_{1}$ is any real constant and where ${}_{0}Y_{1}$ is a spherical harmonic (see Section B). In particular, we have the important equivalence

[TABLE]

We can therefore conclude that the resulting initial data set is asymptotically flat (Definition 1) if and only if $\mathcal{C}=0$ and $q^{(2)}_{1}=0$. While the condition $\mathcal{C}=0$ alone implies that the $3$-metric is asymptotically Euclidean, the initial data set is only asymptotically flat if in addition $q^{(2)}_{1}=0$. Notice that according to [3] the ADM mass is therefore only uniquely defined if both $\mathcal{C}=0$ and $q^{(2)}_{1}=0$. It is remarkable however that the limit of the Hawking mass Eq. (4.9) with respect to our foliation of $2$-spheres always exists. Whether this limit is unique or whether it depends on our particular foliation of $2$-spheres remains open.

In order to support this formal asymptotic analysis we show numerical results for $M_{-}=2/3$, $M_{+}=1/3$ and $Z=1$ in Fig. 3. In order to check our prediction that $A\rightarrow\underline{A}$ in leading-order in the limit $\rho\rightarrow\infty$, the first plot in Fig. 3 shows the quantity

[TABLE]

where we have used Eqs. (B.9), (B.10) and (B.11). Our results suggest that this quantity should be $O(\rho^{-2})$ for large $\rho$. Our numerical findings are certainly consistent with this. Next we check our prediction that $q\rightarrow\underline{q}$ in leading-order in the limit $\rho\rightarrow\infty$ by plotting the quantity

[TABLE]

in Fig. 3. Again we find agreement with this prediction.

Finally, we also provide numerical evidence to support Eq. (4.3) in4.44.44.4This plot uses the ‘symlog’ option in matplotlib in Python, which for a function $f(x)$ essentially plots $\text{sign}(f(x))\log(|f(x)|)$, for details see for instance [39]. Fig. 4 where we plot $\rho p$ as a function of $\rho$ demonstrating that it remains finite for large values of $\rho$. Recall that based on Eqs. (4.10) and (4.11) we do not expect $p_{a}^{(1)}$, and therefore $q^{(2)}_{1}$, to vanish even in the asymptotically Euclidean case $\mathcal{C}=0$.

4.2 Numerical determination of the asymptotic parameters

Given any initial data set as in Section 3, how would we calculate the asymptotic parameters identified in Section 4.1 numerically?

The results in Section 4.1 suggest that

[TABLE]

and $\underline{q}$ therefore agrees with the spherically symmetric case in the two leading orders. Calculating $\underline{q}(\rho)$ for some sufficiently large $\rho$ therefore allows us to numerically estimate the cone parameter $\mathcal{C}$ as follows. We write

[TABLE]

Defining

[TABLE]

it follows that $\mathcal{C}_{N}(\rho)\rightarrow\mathcal{C}$ in the limit $\rho\rightarrow\infty$, in fact,

[TABLE]

So $\mathcal{C}_{N}(\rho)$ in Eq. (4.13) can be understood as a numerical approximation of the asymptotic parameter $\mathcal{C}$ which converges to $\mathcal{C}$ in the limit $\rho\rightarrow\infty$ with the rate $O(\rho^{-2})$. This is confirmed in the first plot of Fig. 5.

In a similar way the above results yield

[TABLE]

Hence,

[TABLE]

In analogy to the above, the quantity

[TABLE]

is therefore a numerical approximation of $\underline{A^{(1)}}$, in fact,

[TABLE]

using that

[TABLE]

as a consequence of Eq. (4.14). Eq. (4.16) is confirmed in the second plot of Fig. 5. This can be used to estimate the limit $m$ of the Hawking mass in Eq. (4.9). We have

[TABLE]

The quantity

[TABLE]

is therefore a numerical approximation of $m$ and

[TABLE]

This decay is confirmed in the third plot in Fig. 5

Even though the estimates $\mathcal{C}_{N}(\rho)$, $m_{N}(\rho)$ for $\mathcal{C}$, $m$ etc. are better the larger the value of $\rho$ is at which we calculate $\mathcal{C}_{N}(\rho)$, $m_{N}(\rho)$, we find that numerical errors in numerically solving the constraint equations become significant when we go further than $\rho\sim 10^{3}$ (for $\epsilon=12$ and $N=11$; see the end of Section 3.3.2). The question is therefore how good the values $\mathcal{C}_{N}(\rho)$ and $m_{N}(\rho)$ at, say, $\rho\sim 10^{3}$ are as approximations for the actual asymptotic parameters $\mathcal{C}$ and $m$. For any of these quantities $\mu$ and the corresponding function $\mu_{N}(\rho)$, we consider

[TABLE]

calculated at some large value of $\rho$ as a measure for the absolute error. Consider $M_{+}=M_{-}=\frac{1}{2}$ and $Z=1$ as an example. Our numerical evolution is carried out for $\rho_{0}=3$. Our estimation of the asymptotic parameters at $\rho=10^{3}$ yields

[TABLE]

Notice that the relative error is of order $\sim 10^{-6}$. As mentioned above, this is likely due to the error of measuring $\mathcal{C}$ and $m$ at a finite point and not at infinity. However, due to the errors in numerically solving the constraints we find that we cannot go further than $\rho\sim 10^{3}$. So we need to accept whatever error we have at that point in the measurement of the matching parameters.

4.3 Iterative construction of asymptotically Euclidean binary black hole data

Consistent with our previous work [11] we have seen here also that solving the constraints as an initial value problem has the major drawback of giving no control over the asymptotics of the solutions. With incorrect or unphysical asymptotics however the resulting initial data sets may not have any reasonable physical interpretation. The last example in Section 4.2 for instance is not asymptotically flat (neither in the weak nor the strong sense) – a property one would expect for any compact isolated gravitationally bound astrophysical system. In fact it is not even asymptotically Euclidean. The problem is that the relationship between the freely specifiable quantities $M_{+}$, $M_{-}$, $Z$ and $\mathcal{F}$ and the resulting corresponding asymptotic quantities $\mathcal{C}$, $m$ and $q^{(2)}_{1}$ is highly nonlinear and nonlocal, and, therefore hard to analyse. As described earlier, this involves first the construction of a preliminary data set for given $M_{+}$, $M_{-}$, $Z$ and $\mathcal{F}$ (i.e., Step 1 in Section 3.3.1), second, the numerical evolution (Step 2 in Section 3.3.2) for as large values of $\rho$ as possible, in order to determine the asymptotic parameters as in Section 4.2.

Now we want to address the following question: Suppose that for some choice of $M_{+}$, $M_{-}$, $Z$ and $\mathcal{F}$ the corresponding asymptotic parameters found as described above are not “favourable”, in particular consider the case that $\mathcal{C}$ does not vanish. How can we change the free quantities $M_{+}$, $M_{-}$, $Z$ and $\mathcal{F}$ in order to “improve the data set”? In fact we shall now discuss an iterative numerical procedure which allows us to decrease $|\mathcal{C}|$ step by step and thereby, in principle, make $|\mathcal{C}|$ as small as we like. We shall see that the function $\mathcal{F}$ introduced in Eq. (3.24) is crucial for this. Recall our discussion in Section 4.1 which yields that $\mathcal{C}=0$ is not sufficient for asymptotic flatness because $q^{(2)}_{1}$ must vanish as well. In this section now we nevertheless restrict our attention to the cone parameter $\mathcal{C}$. We believe that similar ideas apply to deal with $q^{(2)}_{1}$ as well.

The method is based on the following observation. It is a consequence of Eqs. (2.30), (2.31), (2.32), (3.24) and the procedure in Section 3.3.1 that the only quantity of the preliminary data set constructed in Step 1 of our procedure which is affected by the function $\mathcal{F}$ is ${}^{[A]}q$. In fact, if we consider the parameters $M_{+}$, $M_{-}$ and $Z$ as fixed now and consider two preliminary data sets given by two different choices $\hat{\gamma}$ and $\tilde{\gamma}$ of a real parameter $\gamma$ introduced by

[TABLE]

Eq. (2.32) yields that

[TABLE]

It is clear that such two different preliminary data sets lead to two different initial data sets in Step 2 of our method in Section 3.3.2. According to Eqs. (4.2) and (4.8), it can be expected that both resulting data sets have two different cone parameters $\hat{\mathcal{C}}$ and $\tilde{\mathcal{C}}$ and

[TABLE]

At least in certain regimes we could therefore expect that

[TABLE]

for some positive approximately constant quantity4.54.54.5The heuristic argument here might suggest that $\nu\approx 1$. The following iteration scheme turns out to converge faster if we pick $\nu$ as in Eq. (4.21) below. $\nu$. Based on this, we propose the following iterative scheme4.64.64.6We have not yet tried to use any classical root finding method like the bisection or the secant method. It is possible that some of these converge faster or are more reliable than the one discussed here.. As above, fix the free quantities $M_{+}$, $M_{-}$, $Z$, and set

[TABLE]

First calculate the full initial data set as in Section 3.3 for $\gamma_{0}=0$, and determine the corresponding cone parameter $\mathcal{C}_{0}$ as in Section 4.2. Then pick $\gamma_{1}=\nu\mathcal{C}_{0}$, and, again calculate the full initial data set and determine the cone parameter $\mathcal{C}_{1}$. If the (purely experimental) choice of $\nu$ in Eq. (4.21) was “correct” and Eq. (4.20) was therefore exact, we would find $\mathcal{C}_{1}=0$. We would have therefore achieved the goal of “improving” the cone parameter and find an initial data set with the optimal value $\mathcal{C}_{1}=0$. In general, however, the resulting value of $\mathcal{C}_{1}$ will in general not be zero. Numerical evidence suggests that $|\mathcal{C}_{1}|<|\mathcal{C}_{0}|$. We therefore repeat the iteration as often as necessary until eventually the resulting cone parameter is sufficiently close to zero. In summary, our proposed iterative procedure is defined as follows. Fix $M_{+}$, $M_{-}$, $Z$ and set $\nu$ as in Eq. (4.21). Pick an accuracy goal parameter $\mu>0$.

Start condition:

Pick $\gamma_{0}=0$. Calculate the corresponding full initial data set as in Section 3.3 and determine the corresponding cone parameter $\mathcal{C}_{0}$ as in Section 4.2.

Iterative step:

Let $n\geq 0$. Suppose that we have determined $\gamma_{n}$, the corresponding full initial data set and the corresponding cone parameter $\mathcal{C}_{n}$. If $|\mathcal{C}_{n}|<\mu$, stop here. Otherwise, set

[TABLE]

and repeat the iterative step with $n$ replaced by $n+1$.

Again notice that the motivation for Eq. (4.22) is that if Eq. (4.20) was exact with $\nu$ as in Eq. (4.21), then $\mathcal{C}_{n+1}=0$. The scheme is illustrated in Fig. 6.

We shall now apply this to three cases as summarised below. For all these three cases we choose $\rho_{0}=3$ and estimate the asymptotic parameters at $\rho=10^{3}$. We always pick $\epsilon=12$, $N=11$ and $\mu=10^{-9}$. Notice that due to the fixed choice $\rho_{0}=3$, Eq. (3.21) gives a restriction of how much we can vary the parameters $M_{+}$, $M_{-}$ and $Z$.

1. $M_{+}=1/2$, $M_{-}=1/2$, $Z=1$:

For this case, only two iterations are needed until $\mathcal{C}$ has a sufficiently small magnitude4.74.74.7We write $\mathcal{C}=\mathcal{C}_{n}$ and $\gamma_{f}=\gamma_{n}$ where $n$ is the last iteration step.:

[TABLE]

The corresponding value of $\gamma_{f}$ is

[TABLE]

and the mass is

[TABLE] 2. $M_{+}=2/3$, $M_{-}=1/3$, $Z=1$:

A total of eight iterations are needed here to give a sufficiently small cone parameter

[TABLE]

with

[TABLE]

and mass

[TABLE] 3. $M_{+}=1/2$, $M_{-}=1/2$, $Z=3/2$:

Fourteen iterations are used to get the values:

[TABLE]

Here

[TABLE]

and

[TABLE]

Fig. 7 shows colour maps of the solutions, for each case, in “Cartesian” coordinates.

[TABLE]

We wish to point out several things. First, it turns out that in none of the three cases above the resulting quantity $q^{(2)}_{1}$ is zero; the initial data sets are therefore not asymptotically flat. Therefore, one should not regard them as physical. Second, we have noticed that the resulting masses $m$ (the well-defined limits of the Hawking mass) only depend weakly on the mass parameters $M_{+}$ and $M_{-}$ as long as $M_{+}+M_{-}$ is kept constant (here we choose $M_{+}+M_{-}=1$). The masses $m$, however, depend strongly on $Z$. In order to further study the dependence of $m$ and $\gamma_{f}$ on $Z$ we fix $M_{+}=M_{-}=1/2$ and solve the equations for various separation distances $Z\in[0.1,1.5]$. Both $\gamma_{f}$ and $m$ are shown as functions of $Z$ in Fig. 8.

In the left plot in Fig. 8 we see that $\gamma_{f}$ is an increasing function of $Z$. By fitting a fourth order polynomial to the numerical values of $\gamma_{f}$ we are able to interpolate this function. The right plot shows the difference of the values of $\gamma_{f}$ determined by the iteration scheme ($\prescript{(N)}{}{\gamma}_{f}$) and the value given by the fitted polynomial ( $\prescript{(P)}{}{\gamma}_{f}$). Regarding the second plot in Fig. 8, it is interesting to notice that $m$ is a decreasing function of $Z$. For $Z=0$, i.e., the single black hole case, we get $m=1$ as expected. When the separation parameter $Z$ is larger, the mass decreases. This is counter-intuitive, particularly when compared to the Newtonian case, where the gravitational binding energy should become small as $Z$ increases. In fact, even within GR it is expected [19] that the interaction energy of a binary black hole system behaves like

[TABLE]

which should increase the total energy $m=M_{+}+M_{-}+E$ with increasing $Z$. There are several possibilities to interpret this observation: one could simply dismiss this consideration on the basis that the data sets we compute are only asymptotically Euclidean and not asymptotically flat. So they do not provide data for a physical vacuum system but for some system that may have a source at infinity. A second possibility for the behaviour could be that the initial data do not correspond to “two black holes” since we have not demonstrated that there are apparent horizons in the data set. Finally, if we accept that these data sets do possess some physical relevance, then it could still be argued that we are not in the regime in which the asymptotic formula holds. Whatever the true answer is, it is clear that this phenomenon requires a better understanding and more numerical work.

We end with a note about the convergence rate of our iterative scheme. For Case 1 (defined by $2M_{+}=2M_{-}=Z=1$) it took $126$ iterations to get the cone parameter to a magnitude of $10^{-12}$. This is a marked increase from the two iterations to get it to $\sim 10^{-9}$. Use of more advanced root finding methods such as the Newton method or the bisection method would likely increase the rate of convergence.

5 Conclusions

We construct initial data sets using a particular adaptation of [15] to the parabolic-hyperbolic formulation of the constraint equations by Rácz in a multiple black hole setting. Similar to [11], where we constructed perturbed single black hole initial data sets with the algebraic-hyperbolic formulation of the constraint equations, we find that our initial data sets here are in general not asymptotically flat. While the initial data sets in our previous work are at least asymptotically Euclidean (because the full $3$-metric can be prescribed freely), the initial data sets here are asymptotically cone-like unless the asymptotic quantity $\mathcal{C}$ defined above vanishes. This and the other relevant asymptotic quantity $q^{(2)}_{1}$ are determined by the free data and initial data in a complicated non-linear and non-local fashion which is far from being understood. It is therefore impossible to predict which choice of initial data and free data implies $\mathcal{C}=q^{(2)}_{1}=0$ which is necessary for asymptotic flatness. Our results suggest that is only the case for a set of initial data sets of measure zero.

Given this we then demonstrate by incorporating some further degrees of freedom (in particular non-vanishing values for the function $\mathcal{F}$) that asymptotically Euclidean initial data sets can be approximated by an iterative scheme. We have not incorporated the quantity $q^{(2)}_{1}$ into this scheme yet, but are hopeful that similar ideas can be applied to iteratively approximate asymptotically flat initial data sets in a similar way. This could be done by taking advantage of further degrees of freedom present in the choice of $\dot{\gamma}_{ab}$ in Eq. (3.24).

Now, even if one of our initial data set is asymptotically flat, it is not clear whether it represents black holes. The potential for demonstrating the existence of apparent horizons, which would indicate the presence of black holes, is however restricted by the limitation Eq. (2.15) imposed by Rácz’s formalism for foliations based on $2$-spheres. According to this, the constraints must be evolved outwards and can therefore not be used to explore the interior geometry. In contrast to this, the restriction implied by Eq. (3.21) as a consequence of the particular “peanut shape” of our $2$-spheres, which prevents us from putting the “initial” inner boundary arbitrarily close to the “black holes”, can be overcome straightforwardly. In any case, for this paper we have not made any efforts to search for (pieces of) apparent horizons in our data sets. The focus of this paper is purely on the asymptotics at spacelike infinity. We have not investigated other potential remedies for this. One possibility could be to choose the initial value for the lapse $A$ as negative which may allow us to integrate inwards. Another solution could be to do the inwards-integration from the initial $2$-surface with the algebraic-hyperbolic formulation of the constraints, while the outwards-integration is performed with the parabolic-hyperbolic formulation as before.

Acknowledgements

JR was supported by a Master scholarship awarded by the University of Otago. Part of this research was funded by a grant to JF from the Division of Sciences of the University of Otago.

\appendixpage

Appendix A Notation and conventions comparison

In this work we have chosen our sign conventions as in the texts [1, 24]. We have also attempted to simplify the notation, which differs from both Rácz’s work and from our previous work [11]. The following table provides a comparison of the different conventions:

$K_{ab}$ $\kappa$ $p_{a}$ $Q_{ab}$ $q$ $\gamma_{ab}$ $h_{ab}$ $A$ $B_{a}$ $N_{a}$ $v_{a}$ $k_{ab}$ $\nabla_{a}$ $D_{a}$

[11] $-\chi_{ab}$ $-\boldsymbol{\kappa}$ $-\boldsymbol{k}_{a}$ $-\mathring{K}_{ab}$ $-\boldsymbol{K}$ $h_{ab}$ $\hat{\gamma}_{ab}$ $\hat{N}$ $\hat{N}_{a}$ $\hat{n}_{a}$ $\dot{\hat{n}}_{a}$ $-\hat{K}_{ab}$ $D_{a}$ $\hat{D}_{a}$

Rácz $-\chi_{ab}$ $-\boldsymbol{\kappa}$ $-\boldsymbol{k}_{a}$ $-\mathring{K}_{ab}$ $-\boldsymbol{K}$ $h_{ab}$ $\hat{\gamma}_{ab}$ $\hat{N}$ $\hat{N}_{a}$ $\hat{n}_{a}$ $\dot{\hat{n}}_{a}$ $-\hat{K}_{ab}$ $D_{a}$ $\hat{D}_{a}$

Appendix B Spin-weight and spin-weighted spherical harmonics

We say that a function $f$ defined on $\mathbb{S}^{2}$ has spin-weight $s$ if it transforms as $f\to e^{\text{i}s\xi}f$ under a local rotation by an angle $\xi$ in the tangent plane at any point in $\mathbb{S}^{2}$. Let $(\vartheta,\varphi)$ be standard polar coordinates on $\mathbb{S}^{2}$. If $f$ has spin-weight $s$ and is sufficiently smooth, it can be written as

[TABLE]

where ${}_{s}Y_{lm}(\vartheta,\varphi)$ are the spin-weighted spherical harmonics (SWSH) and where $f_{lm}$ are complex numbers. Using the conventions in [31, 8, 9, 10, 12, 11], these functions satisfy

[TABLE]

where $\delta_{lm}$ is the Kronecker delta and $d\Omega$ is the area element of the metric of the round unit sphere. Using this we find that the coefficients $f_{lm}$ in Eq. (B.1) can be calculated as

[TABLE]

The eth-operators $\eth$ and $\eth^{\prime}$ are defined by

[TABLE]

for any function $f$ on $\mathbb{S}^{2}$ with spin-weight $s$. We have

[TABLE]

Thus, using the properties above it is easy to see that $\eth$ raises the spin-weight by one while $\eth^{\prime}$ lowers it by one.

In our discussion we are often interested in the average of a function $f$ with spin-weight [math] on $\mathbb{S}^{2}$ defined by

[TABLE]

Expressing $f$ in terms of SWSH and using Eq. (B.2) it follows

[TABLE]

where we have used the fact that ${}_{0}Y_{00}(\vartheta,\varphi)=(4\pi)^{-1/2}$. Another quantity of interest is the $L^{2}$-norm with respect to the standard round metric on $S^{2}$. The Parseval identity states that

[TABLE]

Finally we notice that all functions considered in this paper are axially symmetric and therefore do not depend on the angle $\varphi$. For such functions, all coefficients with $f_{lm}$ with $m\not=0$ vanish and we use the following short-hand notation to write Eq. (B.1) as

[TABLE]