Numerical initial data deformation exploiting a gluing construction: I. Exterior asymptotic Schwarzschild

TL;DR

This paper introduces a new spectral numerical method for constructing initial data in Einstein's equations by gluing binary black hole solutions to exterior Schwarzschild data without using the conformal approach.

Contribution

It presents a novel spectral scheme inspired by exterior asymptotic gluing that constructs initial data with binary black holes attached to Schwarzschild exterior, bypassing the conformal method.

Findings

Successfully constructs initial data with binary black holes and Schwarzschild exterior.

Demonstrates the method's effectiveness through ADM mass comparisons.

Shows the approach's potential for more general initial data configurations.

Abstract

In this work a new numerical technique to prepare Cauchy data for the initial value problem (IVP) formulation of Einstein's field equations is presented. Directly inspired by the exterior asymptotic gluing (EAG) result of Corvino (2000) our (pseudo)-spectral scheme is demonstrated under the assumption of axisymmetry so as to fashion composite Hamiltonian constraint satisfying initial data featuring internal binary black holes (BBH) as glued to exterior Schwarzschild initial data in isotropic form. The generality of the method is illustrated in a comparison of the ADM mass of EAG initial data sets featuring internal BBHs as modelled by Brill-Lindquist and Misner data.In contrast to the recent work of Doulis and Rinne (2016), and Pook-Kolb and Giulini (2018) we do not make use of the York-Lichnerowicz conformal framework to reformulate the constraints.

Figures

Equations

$$g_S={\left(1+\frac{M_{\mathrm{ADM}}}{2|\mathbf{x}-\mathbf{C}|}\right)}^4{\delta }_{\mathrm{Euc}},$$

$${\overline{g}}_{{\Omega }_{\rho }}:=\chi g_E+(1-\chi )g_S.$$

$$\mathcal{R}[{\overline{g}}_{{\Omega }_{\rho }}+h]\simeq \mathcal{R}[{\overline{g}}_{{\Omega }_{\rho }}]+L_{{\overline{g}}_{{\Omega }_{\rho }}}[h],$$

$$L_{\overline{g}}[h]=-h^{ij}{\overline{\mathrm{Ric}}}_{ij}-{\overline{\nabla }}^2\left[h^k{}_k\right]+{\overline{\nabla }}_i{\overline{\nabla }}_j{h}^{ij}.$$

$${\left(L_{\overline{g}}^{\ast }[f]\right)}_{ij}=-{\overline{\mathrm{Ric}}}_{ij}f-{\overline{g}}_{ij}{\overline{\nabla }}^2\left[f\right]+{\overline{\nabla }}_{(i}{\overline{\nabla }}_{j)}[f],$$

$$\mathcal{V}[u]={\int }_{\Omega }\left(\frac{1}{2}\Vert L_{\overline{g}}^{\ast }u{\Vert }^2-\tilde{\delta }u\right)\omega \mathrm{d}{\mu }_{\overline{g}},$$

$$0={\frac{d}{dt}\left[\mathcal{V}[u+t\eta ]\right]|}_{t=0}⟹{\int }_{\Omega }{\left[L_{\overline{g}}^{\ast }\eta \right]}^{ij}{\left[L_{\overline{g}}^{\ast }u\right]}_{ij}\omega \mathrm{d}{\mu }_{\overline{g}}={\int }_{\Omega }\eta \tilde{\delta }\omega \mathrm{d}{\mu }_{\overline{g}},$$

$$L_{\overline{g}}\left[\omega L_{\overline{g}}^{\ast }\left[u\right]\right]=\tilde{\delta }\omega =\delta ,\left(\delta \in C_{{\omega }^{-1}}^{k,\alpha }(\Omega )\right),$$

$$u$$

$${}^{(k)}\delta =\mathcal{R}[g_{ij}]-\mathcal{R}\left[{\overline{g}}_{ij}+\omega \sum _{m=0}^{k-1}{}^{(m)}{\tilde{h}}_{ij}\right],$$

$$L_{{\delta }_{\mathrm{Euc}}}^{\ast }[f]{}^i{}_j=-\delta {}^i{}_j{\partial }^k[{\partial }_k[f]]+\partial {}^i\partial {}_j[f],$$

$$g_{ij}=(1-N^{-2}{N}_K{N}^K)n_i{n}_j+2N^{-1}n{}_{(i}{N}_{j)}+\delta {}_i^I\delta {}_j^J{\gamma }_{IJ},$$

$${}_{s}T:=T(\underbrace{\boldsymbol{\omega},\,\cdots\boldsymbol{\omega}}_{\mbox{\small$q_{1}$ times}};\,\underbrace{\overline{\boldsymbol{\omega}},\,\cdots,\,\overline{\boldsymbol{\omega}}}_{\mbox{\small$q_{2}$ times}};\,\underbrace{\mathbf{m},\,\cdots\mathbf{m}}_{\mbox{\small$r_{1}$ times}};\,\underbrace{\overline{\mathbf{m}},\,\cdots,\,\overline{\mathbf{m}}}_{\mbox{\small$r_{2}$ times}}),$$

$$\stackrel{\circ }{\gamma }$$

$$m^I=\frac{1}{\sqrt{2}\text{lower}0.4pt\text{vrule}height=6.44444pt,depth=-5.15558pt}\left({\partial }_{\vartheta }^I-i\csc \vartheta {\partial }_{\phi }^I\right)⟹{\omega }_I=\frac{1}{\sqrt{2}\text{lower}0.4pt\text{vrule}height=6.44444pt,depth=-5.15558pt}\left(\mathrm{d}{\vartheta }_I+i\sin \vartheta \mathrm{d}{\phi }_I\right),$$

$$m^I{\overline{\omega }}_I=$$

$$\stackrel{\circ }{\gamma }{}_{IJ}=$$

$$\mathbf{m}\left[{}_{s}f\right]=$$

$$ð\left[{}_{s}f(\vartheta ,\phi )\right]=$$

$$\overline{ð}\left[{}_{s}f(\vartheta ,\phi )\right]=$$

$$\sqrt{2}\text{lower}0.4pt\text{vrule}height=6.44444pt,depth=-5.15558pt{m}^J{m}^{I_1}\cdots m^{I_n}{\stackrel{\circ }{\mathfrak{D}}}_J\left[W_{(I_1\cdots I_n)}\right]=$$

$$\gamma {}_{IJ}={}_{-2}\gamma {\overline{\omega }}_I{\overline{\omega }}_J+2{}_0\gamma {\omega }_{(I}{\overline{\omega }}_{J)}+{}_{+2}\gamma {\omega }_I{\omega }_J,$$

$${\gamma }^{IJ}={}_0\tilde{\gamma }\left(-{}_{-2}\gamma m^I{m}^J+2{}_0\gamma m^{(I}{\overline{m}}^{J)}-{}_{+2}\gamma {\overline{m}}^I{\overline{m}}^J\right),$$

$${}_0\tilde{\gamma }={\left({}_0{\gamma }^2-{}_{-2}\gamma {}_{+2}\gamma \right)}^{-1}=\det ({\stackrel{\circ }{\gamma }}_{IJ})/\det \left({\gamma }_{IJ}\right).$$

$${\overline{\omega }}^I={\gamma }^{IJ}{\omega }_J=$$

$${\gamma }_{IJ}=\left[\begin{array}{cc}{\Gamma }_1& {\Gamma }_2\\ {\Gamma }_2& {\Gamma }_3{\sin }^2\vartheta \end{array}\right],$$

$${}_{-2}\gamma =$$

$$C^K{}_{IJ}=C^K{}_{(IJ)}=\frac{1}{2}{\gamma }^{KL}\left({\stackrel{\circ }{\mathfrak{D}}}_I[{\gamma }_{LJ}]+{\stackrel{\circ }{\mathfrak{D}}}_J[{\gamma }_{IL}]-{\stackrel{\circ }{\mathfrak{D}}}_L[{\gamma }_{IJ}]\right).$$

$${\mathfrak{D}}_I[V_J]=$$

$${\mathfrak{D}}_I[T_{JK}]=$$

$$2\sqrt{2}\text{lower}0.4pt\text{vrule}height=6.44444pt,depth=-5.15558pt{C}_{KIJ}=$$

Full text

1. Introduction

Gluing techniques provide for a powerful method of geometric analysis (GA) which may be exploited to combine multiple, distinct, solutions to a (system of) PDE of interest through their gradual deformation over some open set $\Omega$ so as to furnish a new, composite solution that approximately coincides with the original solutions away from $\Omega$ [16, 25, 17].

In this work we focus on the vacuum Einstein constraint equations $\mathcal{C}^{\mathrm{ADM}}[g,\,K]=0$ [22, 5, 72, 44, 80]. For concreteness, recall that the constraints split into the scalar Hamiltonian constraint $\mathcal{H}^{\mathrm{ADM}}[g,\,K]=0$ and the vectorial momentum constraint $\mathcal{M}{}_{j}^{\mathrm{ADM}}[g,\,K]=0$. In general, these are to be satisfied on a Riemannian manifold $\Sigma$ by a spatial metric $g_{ij}$, together with extrinsic curvature $K_{ij}$. From the perspective of the (numerical) evolution problem, a triplet $(\Sigma,\,g_{ij},\,K_{ij})$ constitutes an initial data set.

Consider a moment-in-time (MIT) symmetry (where $K_{ij}=0$), such that the constraints reduce to the single, scalar-flat condition $\mathcal{R}[g]=0$. In this setting, by exploiting a GA based technique of scalar curvature deformation, Corvino [24] has shown the following: let $g$ be an arbitrary asymptotically flat metric on $\mathbb{R}^{3}$, satisfying the Hamiltonian constraint $\mathcal{R}[g]=0$, with positive ADM-mass $m_{0}$. Then there exists another asymptotically flat metric $\hat{g}$ satisfying the constraint which agrees with $g$ on a compact set $K_{1}$ and is identical to a Schwarzschild solution (with, in general, a different ADM-mass and a shifted centre of mass) outside another compact set $K_{2}$ with $K_{1}\Subset K_{2}$. Thus, the new metric may be regarded as a composite metric obtained by gluing the original metric to the Schwarzschild metric in the transition region $\Omega:=K_{2}\backslash K_{1}$. The truly novel feature of this exterior asymptotic gluing (EAG) construction is that the new composite initial data set exactly coincides with its respective constituents outside the “transition region” $\Omega$ (which, as an example, may be imagined to be a spherical shell of finite thickness). This possibility of local gluing, where the region over which two initial data sets are spliced together is of compact support, is entirely due to the underdeterminedness of the constraint equations. The assumption of asymptotic flatness for the interior metric seems to have been made for technical reasons in order to guarantee that the non-linear operator in question is surjective. In the present work we stick with this assumption and glue only metrics which are asymptotically flat. It would be interesting, however, to see whether Corvino’s approach would also allow us to glue arbitrary scalar-flat metrics.

By relaxing the MIT condition and applying a similar gluing strategy it has been shown that the Corvino-Schoen technique may also be used to glue to exact Kerr exteriors [26]. Thus, quite general interior gravitational configurations may be glued to exterior Schwarzschild or Kerr regions forming a composite solution with precise asymptotics111Sacrificing local control on solution character whilst engineering asymptopia without gluing has also been investigated [4]..

A striking variant of the above, where $\Omega$ is replaced by a conical region of infinite extent, is Carlotto-Schoen gluing [15] (see also [17]). In principle, effective screening is allowed for by manipulation of vacuum initial data alone. Furthermore, the Corvino-Schoen and Carlotto-Schoen gluings may be utilised so as to construct $N$-body initial data sets [18].

Related to the above is the method of connected sum or IMP gluing [53, 52, 54]. Here the conformal (Lichnerowicz-York) framework is adopted and consequently a determined elliptic system results. Given any two solutions of the constraints: $(\check{\Sigma}_{0},\,\check{g}_{ij},\,\check{K}_{ij})$ and $(\hat{\Sigma}_{0},\,\hat{g}_{ij},\,\hat{K}_{ij})$, say, a new solution may be produced by first removing small neighbourhoods $\check{N}$ and $\hat{N}$ about the points $\check{p}\in\check{\Sigma}_{0}$ and $\hat{p}\in\hat{\Sigma}_{0}$ respectively. Then, new data $(\Sigma_{0},\,g_{ij},\,K_{ij})$ is found by connecting $\partial\check{N}$ along an interpolating tube to $\partial\hat{N}$ with $\Sigma_{0}$ resulting in a connected sum manifold with the topology of $\check{\Sigma}_{0}\#\hat{\Sigma}_{0}$. By suitable interpolation of the pairs $(\check{g}_{ij},\,\check{K}_{ij})$ and $(\hat{g}_{ij},\,\hat{K}_{ij})$ a composite constraint satisfying solution $(\Sigma_{0},\,g_{ij},\,K_{ij})$ may be found.

Alternatively, identifying $\check{\Sigma}_{0}$ and $\hat{\Sigma}_{0}$ allows for a handle (wormhole) to be introduced to a given initial data set. On account of the determinedness this leads to a global deformation of the initial data set which is small away from the gluing site. By combining the results of [24, 19] it was shown in [20, 21] how the deformation may be localised. To date there do not appear to have been attempts made to prepare numerical initial data based on the IMP approach.

Aside from the ability to control asymptopia of initial data sets, engineering of exotic properties is interesting in its own right and, indeed, a robust scheme for fashioning numerical solutions could prove useful in endowing the associated space-time with a particular, desired phenomenology. For example, inspired by the result of [24], a potential path towards minimisation of so-called “spurious” gravitational radiation content was provided in [41]. We refer the reader there for further details.

While the above gluing results have intriguing properties, an unfortunate aspect is that the GA flavour of proof technique is quite technical in nature. It is not entirely clear how to proceed if direct numerical preparation of an initial data set based on such methods is desired. This is evidenced by the fact that there exists only a single attempt [41] based on formal perturbation theory to “embed” within the conformal framework a problem that seeks to mimic the setup of Corvino’s result [24].

Briefly, the idea in [41] was to work at an MIT symmetry, assume axisymmetry and fix internal Brill-Lindquist (BL) data. Then, over an annular $\Omega$ a conformally transformed Brill-Wave [12] ansatz on the form of the conformal factor $\psi$ is made. It was claimed that composite solutions exist to this problem when the exterior is a suitably chosen Schwarzschild initial data set. This approach requires a further, ad hoc treatment of the decay rates of $\psi$ as $\partial\Omega$ is approached. Following this programme, it appears that a numerical solution may be constructed [28] (see however, the modified, Newton-Krylov based approach of [74]).

An additional insight is offered in [28] as to how consistent selection of exterior data (or parameters on valid internal data) may be made by exploitation of an integrability condition. Such arguments are not required in the proof of [24]. It does not appear that numerical evolution has been performed based on the results of [28, 74]. Indeed, we are not aware of any numerical evolution of initial data sets which have been prepared based on gluing techniques.

More broadly, the technique of scalar curvature deformation may be of potential interest in studies involving geometric curvature flow. Such flows were introduced to general relativity in [37] and the idea built upon in [55] to rule out a class of counterexamples to the cosmic censorship hypothesis proposed by Penrose in [70] as encapsulated by an inequality relating black hole (ADM) mass and the area of its apparent horizon. The veracity of this inequality in a special case was first established rigorously in [51] by exploiting the inverse mean curvature flow for MIT data sets. The more general problem without this restriction remains open and numerical investigation utilising the weak formulation approach of [51] may help shed light on the matter. The related Ricci flow [47] has also been studied numerically [35, 78] where in [35] preliminary evidence for critical behaviour along the flow was presented. Another potentially novel scenario to consider may be whether scalar curvature deformation can be employed as a mechanism to control the appearance of such critical behaviour.

Our goal in this work is to provide some insight as to how the proof in [24] may be more directly adapted to a numerical technique itself without making use of the conformal programme for reformulation of the constraint equations as in [41, 28, 74]. In so doing, we shall numerically construct initial data as composite solutions with an MIT symmetry. To begin, we elaborate upon Corvino’s method at a formal level in §2. Scalar curvature deformation over $\Omega\Subset\Sigma$ and construction of a solution metric describing it, is effected iteratively, through solution of a sequence of linear sub-problems. The basic ingredient of this is described in §2.1. How the iteration is to proceed, together with an obstruction that occurs in the particular case of solving the constraints themselves and our proposed remedy is detailed in §2.2.

Having outlined the problem at the abstract level we next turn our attention to providing a robust description of geometric quantities required for the problem that is suitable for numerical work. To this end, a frame based approach is introduced in §3. In particular, $\Sigma$ is viewed as foliated by topological $2$-spheres. With a view towards efficient numerical implementation, the intrinsic geometry is described through the $\eth$-formalism, which is briefly recounted in §3.1. Details on how it may be adapted to treat topological $2$-spheres are provided in §3.2. For convenience, the relation between intrinsic and ambient quantities adapted to our discussion is touched upon in §3.3.

The success and versatility of pseudo-spectral methods [45] motivates our numerical approach in §4. In particular, function approximation of intrinsic quantities cast in the $\eth$-formalism is discussed in §4.1. Approximation of more general quantities over $\Sigma$ is described in §4.2. For the deformation problem at hand, we supplement the discussion with some complex analytic considerations that can assist in improving numerical solution quality in §4.3.

With particulars of the physical problem and numerical technique fixed we subsequently investigate prototype problems in §5. As an initial test, the case of scalar curvature deformation in the context of spherical symmetry is initially investigated in §5.1, and self-consistent convergence tests performed in §5.2. Following this, a relaxation to the class of axisymmetric problems is set up in §5.3 and explored in §5.4. In §5.5 all previously introduced material is brought together and we perform gluing of internal BBH data (for Brill-Lindquist and Misner initial data) to exterior Schwarzschild initial data. Numerical performance of the approach together with properties of the physical construction are investigated. Finally §6 concludes.

2. Corvino’s method

The argument for solving the Einstein constraints $\mathcal{C}^{\mathrm{ADM}}[g,\,K]=0$ presented in [24] by virtue of exterior asymptotic gluing (EAG) is quite technical in nature and consequently how one should proceed in order to fashion a numerical technique is somewhat opaque. Our goal here is to provide a sketch of the idea adapted to the aforementioned context (cf. the general discussions of [3, 17, 16, 25]). The physical setting is vacuum with vanishing cosmological constant at an MIT symmetry $(K_{ij}=0)$ and in what follows $\Sigma$ is to be understood as an initial Cauchy slice222Here particularised to $\mathrm{dim}(\Sigma)=3$.. Under these assumptions $\mathcal{C}^{\mathrm{ADM}}[g,\,K]=0$ reduces to the single, non-trivial, scalar-flat condition $\mathcal{R}[g]=0$ and $(\Sigma,\,g)$ is sought.

To explain EAG and fix the desired behaviour of $g$ recall that asymptotically flat (AF) data are characterised by the existence of a diffeomorphism between the “end” of $\Sigma$ and $\mathbb{R}^{3}$ with a ball $\mathbb{B}$ removed. Let $\delta_{\mathrm{Euc}}$ be the Euclidean metric. For AF $(\Sigma,\,g_{E})$, end coordinates $\{x^{i}\}_{i=1:3}$ may be introduced such that decay of $g_{E}$ (and derivatives thereof) to $\delta_{\mathrm{Euc}}$ is controlled by negative powers of $|x|$ (see [3, 17]).

The result of [24] concerns an equivalence class of Schwarzschild initial data where a representative in isotropic form is provided by:

[TABLE]

and the $(1+3)$-parameter tuple $(\mathrm{M}_{\mathrm{ADM}},\,C^{i})$ describes the ADM mass and centre of mass. We identify $\Sigma$ with its image in $\mathbb{R}^{3}$. Let $\mathbb{B}_{\rho}\subset\Sigma$ be the ball of radius $\rho>0$. Introduce the compactly contained domain $\Sigma\Supset\Omega_{\rho}:=\mathbb{B}_{2\rho}\setminus\mathbb{B}_{\rho}$ the closure of which is a spherical shell of thickness $\rho$ and serves as a “transition region”. A selection of sufficiently large $\rho$ allows one to smoothly combine any AF $(\mathbb{B}_{\rho},\,{g}_{E})$ satisfying the scalar-flat condition ${\mathcal{R}}[g_{E}]=0$ with $(\Sigma\setminus\mathbb{B}_{2\rho},\,{g}_{S})$ over $\Omega_{\rho}$ via judicious selection of $g_{S}$. This latter is accomplished through tuning of the parameters $(\mathrm{M}_{\mathrm{ADM}},\,C^{i})$ and iterative correction of a smooth, interpolating “background metric” $\overline{g}_{\Omega_{\rho}}$. To understand the procedure, introduce the smooth cut-off function $\chi$ equal to $1$ on $\mathbb{B}_{\rho}$ and [math] outside $\mathbb{B}_{2\rho}$, and on $\Omega_{\rho}$ set:

[TABLE]

Clearly, ${\mathcal{R}}[\overline{g}_{\Omega_{\rho}}]=0$ on $\Sigma\setminus\Omega_{\rho}$ whereas on $\Omega_{\rho}$ we have ${\mathcal{R}}[\overline{g}_{\Omega_{\rho}}]=\delta$ where $\delta$ is a compactly supported function. Furthermore, we shall assume that ${\mathcal{R}}[\overline{g}_{\Omega_{\rho}}]$ is non-constant on $\Omega_{\rho}$ to avoid a technical issue outlined in §7.

To proceed further, the problem is now viewed as a local (i.e., compactly supported) deformation of the scalar curvature. Consider the change $\delta\rightarrow\check{\delta}:={\mathcal{R}}[\overline{g}_{\Omega_{\rho}}]+\tilde{\delta}$ for $\tilde{\delta}$ sufficiently small. The idea is to seek a suitable correction $h$ to the background metric such that ${\mathcal{R}}[\overline{g}_{\Omega_{\rho}}+h]=\check{\delta}$ is satisfied. The approach of [24] is to linearise about the background metric:

[TABLE]

where the linear problem ${L}_{\overline{g}_{\Omega_{\rho}}}[h]=\tilde{\delta}$ is investigated so as to characterise properties of the underdetermined elliptic operator $L_{\overline{g}_{\Omega_{\rho}}}$. Unfortunately, ${L}_{\overline{g}_{\Omega_{\rho}}}$ fails to be injective and the question of surjectivity of ${L}_{\overline{g}_{\Omega_{\rho}}}$ is addressed with a demonstration of injectivity of the formal adjoint ${L}^{*}_{\overline{g}_{\Omega_{\rho}}}$ which is overdetermined. This latter is then utilised, working within weighted function spaces yielding a so-called “basic estimate” over $\Omega_{\rho}$ where a certain growth (decay) rate of functions near $\partial\Omega_{\rho}$ is permitted. Thus, boundary behaviour of functions is implicitly controlled by a weight-function $\omega$ allowing for a variational based solution to the above linear problem. The control on properties of the linear solution turns out to be sufficiently strong to also allow for a Picard iteration scheme to obtain the solution $h$ to the nonlinear local deformation problem.

In principle, the problem of finding $h$ with ${\mathcal{R}}[\overline{g}_{\Omega_{\rho}}+h]=0$ may also be pursued in this way. However, an obstruction exists in that an approach to the scalar-flat condition, as for instance when ${\overline{g}}_{\Omega_{\rho}}\rightarrow{g}_{S}$, induces an approximate, non-trivial kernel $\mathscr{K}[{\overline{g}_{\Omega_{\rho}}}]:=\ker({L}_{\overline{g}_{\Omega_{\rho}}}^{*})$ of the formal linear adjoint. This may be ameliorated through judicious selection of $(M_{\mathrm{ADM}},C^{i})$ so as to work in a space transverse to $\ker({L}_{\delta_{\mathrm{Euc}}}^{*})$ when solving the previously described variational problem at the linearised level. Unfortunately, while [24] demonstrates that such a selection exists a method for a priori specification of the parameters is not provided and hence we instead adopt a direct, numerical linear-algebraic strategy.

We now proceed to provide further details of the variational approach to solving the linearised problem in §2.1 with sufficient detail for our numerical scheme. The nonlinear deformation shall be addressed in §2.2 together with our method for approaching the issue of non-trivial kernel.

2.1. Linear corrections via weak-formulation

For the sake of exposition we shall assume geometric quantities to be defined with respect to $\Omega:=\Omega_{\rho}$. Background quantities will be denoted by an over-bar. Furthermore, we shall assume that $\overline{\mathcal{R}}$ is non-constant. Suppose that a sufficiently small, smooth, local deformation (i.e. of compact support) of the scalar curvature is made $\mathcal{R}[\overline{g}]\rightarrow\mathcal{R}[g]:=\mathcal{R}[\overline{g}]+\tilde{\delta}$ and a metric correction $h$ satisfying $g=\overline{g}+h$ is sought. The problem may be investigated perturbatively by noting that formally $\overline{g}_{ij}\rightarrow g_{ij}=\overline{g}_{ij}+\varepsilon h_{ij}$ induces a corresponding linear-order correction to the scalar curvature $\mathcal{R}[g+\varepsilon h]=\mathcal{R}[\overline{g}]+\varepsilon L_{\overline{g}}[h]$ where standard methods yield [84]:

[TABLE]

Thus, solution of $L_{\overline{g}}[h]=\tilde{\delta}$ is required and hence properties of the linear operator $L_{\overline{g}}$ are investigated in [24]. While it turns out that $L_{\overline{g}}$ is underdetermined elliptic one may instead work with the formal $\mathrm{L}^{2}$ adjoint $\mathrm{L}{}^{*}_{\overline{g}}[\cdot]_{ij}$ as identified from the inner product $\langle\mathrm{L}{}_{\overline{g}}[h],\,f\rangle_{\mathrm{L}^{2}(\Omega)}=\langle h,\,\mathrm{L}{}^{*}_{\overline{g}}[f]\rangle_{\mathrm{L}^{2}(\Omega)}$:

[TABLE]

which is injective [24] (see also §7). Introduce the weighted, Sobolev space functional $\mathcal{V}:H^{2}_{\omega}(\Omega)\rightarrow\mathbb{R}$ defined by:

[TABLE]

where $\tilde{\delta}\in\mathrm{L}^{2}_{\omega}(\Omega)$, $\omega$ is a weight function (to be defined) and $\mathrm{d}{\mu}_{\overline{g}}$ is the integration measure induced by $\overline{g}$. To find the unique $u$ satisfying (2.6), we introduce the test-function $\eta\in C^{\infty}_{c}(\Omega)$ and consider the variation:

[TABLE]

where the compact support of $\eta$ (or alternatively a presumed decay rate for $\omega$ towards $\partial\Omega$) enables us to drop all boundary terms. Equation (2.7) is the so-called weak-formulation [30, 67, 38] of the following strong-form problem [24]):

[TABLE]

where $C^{k,\,\alpha}_{\omega^{-1}}(\Omega)$ is a weighted Hölder space [24]. In light of this equivalence and later use of Eq.(2.7) to iteratively construct $h$ we shall refer to $u$ as a “potential function”. In order to ensure future enforcement of $u\in H{}^{2}_{\omega}(\Omega)$ in Eq.(2.7) at the numerical level — indeed allowing for controlled growth of $u$ towards $\partial\Omega$ — an explicit rewriting exploiting the decay properties of the weight term $\omega$ can be made and viewed as a solution ansatz.

Suppose $x$ is a boundary defining function in a neighbourhood of $\partial\Omega$, i.e., $x\geq 0$ with $x=0$ and $\mathrm{d}x\neq 0$ on $\partial\Omega$. Suppose $\omega\sim x^{2N}$ for $N$ sufficiently large. The condition $\partial{}^{k}_{x}[u]\in L^{2}_{\omega}(\Omega)$ where $k\in\{0,\,1,\,2\}$ leads to:

[TABLE]

where $\mathcal{N}$ is a function with quadratic decay in $x$ towards $\partial\Omega$ and $\tilde{u}$ shall be assumed to be bounded and smooth. We shall defer explicit specification of $\omega$ to §5.

2.2. Nonlinear local $\mathcal{R}$ deformation and gluing

The problem we would now like to solve is: On $\Sigma$ fix a choice of $\overline{g}$ and hence $\overline{\mathcal{R}}$ and $L{}_{\overline{g}}$ together with $\delta$ of compact support on $\Omega$. Assume that $\overline{\mathcal{R}}$ is non-constant. Find a symmetric $2$-tensor $h\in\mathcal{S}_{2}(\Sigma)$ with compact support on $\Omega$ such that $\mathcal{R}[\overline{g}+h]=\mathcal{R}[\overline{g}]+\delta$.

Theorem $1$ of [24] allows us to proceed as follows: Set ${}^{(0)}\delta:=\mathcal{R}[g]-\mathcal{R}[\overline{g}]$, solve Eq.(2.7) for ${}^{(0)}u$, and hence construct ${}^{(0)}h_{ij}=\omega L{}_{\overline{g}}^{*}[{}^{(0)}u]_{ij}$. This yields ${}^{(1)}g_{ij}=\overline{g}_{ij}+{}^{(0)}h_{ij}$ with ${}^{(0)}h_{ij}$ small in an appropriate Hölder space. Now, it would be natural to apply Newton’s method and linearise about the new metric ${}^{(1)}g_{ij}$. However, it turns out [24] that this would apparently result in a loss of differentiability. Instead333It would be interesting to see whether this analytical problem manifests itself also on the numerical level. However, we have not pursued this any further, yet. the proof technique leverages Picard iteration with the linearisation fixed at the background $\overline{g}_{ij}$ and the approximate solution is iteratively improved via:

[TABLE]

where ${}^{(m)}\tilde{h}_{ij}:=L_{\overline{g}}^{*}\left[{}^{(m)}u\right]$ and the solution (metric) is given by $g_{ij}=\overline{g}_{ij}+\lim_{N\rightarrow\infty}\omega\sum_{m=0}^{N}{}^{(m)}\tilde{h}_{ij}$.

An issue remains when $\overline{\mathcal{R}}$ is constant and correspondingly a non-trivial kernel $\mathscr{K}[\overline{g}]:=\ker(L_{\overline{g}}^{*})$ exists. This must be addressed if exterior asymptotic gluing (EAG) to a time-symmetric slice of Schwarzschild $g_{S}$ is to be achieved and a solution to the constraints found. Consider $\overline{g}_{\Omega_{\rho}}$ of Eq.(2.2). The dimension of $\mathscr{K}[\overline{g}_{\Omega_{\rho}}]$ is related to the underlying (approximate) symmetries of $\overline{g}_{\Omega_{\rho}}$ and due to the assumption of $\overline{g}_{\Omega_{\rho}}$ being asymptotically Euclidean must approach that of $\mathscr{K}[\delta_{\mathrm{Euc}}]$ in the asymptotic regime (described below) [16, 25]. It is also known that in linearisation of the full vacuum constraints (no longer at an MIT) a further contribution to the kernel of the corresponding formal linear adjoint arises, which is comprised of the generators of translation and rotation of $\mathbb{R}^{3}$ [16, 25]. Collectively, elements of the non-trivial kernel in this latter case are called Killing initial data (KID) due to the one-to-one correspondence with Killing vectors in the vacuum space-time obtained by evolving the initial data set [66].

Thus for EAG on Schwarzschild we clearly encounter an obstruction as $\rho$ is increased due to the fall-off properties of $\overline{g}_{\Omega_{\rho}}$ and approach to the non-trivial $\mathscr{K}[g_{S}]$. As a preliminary, notice that we can identify $\mathscr{K}_{0}:=\mathscr{K}[\delta_{\mathrm{Euc}}]=\mathrm{span}(1,\,x^{1},\,x^{2},\,x^{3})$. This can be seen by observing Eq.(2.5) implies:

[TABLE]

and hence $L_{\delta_{\mathrm{Euc}}}^{*}[\cdot]_{ij}$ annihilates affine functions of the form $f=a+b_{k}x^{k}$ where $a\in\mathbb{R}$ and $b_{k}\in\mathbb{R}^{3}$. That these are the only possible functions then follows from considering the finite-dimensional space of initial data for Eq.(7.8).

Returning to EAG on Schwarzschild, note that we are only approximately approaching a non-trivial kernel. To account for this [24] proceeds by introducing an approximating kernel $\mathscr{K}_{*}:=\zeta\mathscr{K}_{0}$ where $\zeta$ is a smooth, spherically symmetric bump function of compact support on $\Omega_{\rho}$. The idea is then to solve a projected nonlinear local deformation problem based on $\overline{g}_{\Omega_{\rho}}$ as above but working with functions in the $\mathrm{L}^{2}(\Omega_{\rho})$ orthogonal complement of $\mathscr{K}_{*}$. When carried out at sufficiently large $\rho$ this yields a glued solution $g$ with $\mathcal{R}[g]\in\mathscr{K}_{*}$ and a choice of $(M_{\mathrm{ADM}},\,C^{i})$ is shown to exist (though it is not demonstrated how to select these parameters a priori) such that a $g$ with $\mathcal{R}[g]=0$ may be found.

For our numerical approach a slightly different strategy shall be pursued. An alternative way to construct an appropriately projected problem that treats non-trivial $\mathscr{K}_{0}$ indirectly is provided by linear algebraic techniques. The idea here is to consider evaluation of the weak-formulation statement of Eq.(2.7) with a suitably chosen dense, approximating collection of test space and solution (trial) space functions. A singular value decomposition (SVD) of the ensuing linear system may then be inspected and any (approximate) kernel directly removed [79].

During numerical calculations involving EAG on Schwarzschild (to be performed in §5.5) symmetry conditions shall be imposed. A precise identification of the dimension of the non-trivial kernel in this context may be motivated as follows: Consider the affine functions $f$ as annihilated by $L{}^{*}_{\delta_{\mathrm{Euc}}}[\cdot]_{ij}$. According to [3] the parameters $a$ and $b_{k}$ entering $f$ as above primarily affect how $M_{\mathrm{ADM}}$ and centre of mass $C^{i}$ should be chosen in the composite (numerical) solution. Thus in a context with a high degree of symmetry the effective dimension of the kernel may be reduced.

3. Frame-formalism treatment of $\Sigma$-geometry

With the physical problem and geometric preliminaries outlined in §2 we now turn our attention to concretising the formulation for numerical work. Given a $\Sigma$ with underlying symmetries a chart selection exploiting this property allows for a description of geometric quantities that can lead to more efficient numerical schemes (see §4). In order to accomplish this in a robust fashion, such that issues of regularity do not arise from the choice of coordinatisation we adopt a frame based approach that leverages the so-called $\eth$-formalism. It shall be assumed that $\Sigma$ is endowed with metric $g_{ij}$ and that $(\Sigma,\,g_{ij})$ may be smoothly foliated by a one-parameter family of non-intersecting topological $2$-spheres $\mathbb{S}{}^{2}_{\rho}$ which are to be viewed as the level surfaces of a smooth function $\rho:\Sigma\rightarrow\mathbb{R}$. Denote the Levi-Civita connection associated with $g_{ij}$ by ${\nabla}$.

Following the standard prescription of ADM decomposition adapted to a spatial manifold, the normalised 1-form $n_{i}=N{\nabla}_{i}[\rho]$ provides a normal to $\mathbb{S}{}^{2}_{\rho}$. Recall that the ambient metric induces the metric $\gamma_{ij}$ on the submanifolds $\mathbb{S}{}^{2}_{\rho}$ via $g_{ij}=\gamma_{ij}+n_{i}n_{j}$ and gives rise to the projector $\mathcal{P}{}^{i}{}_{j}:=\gamma{}^{i}{}_{j}-n^{i}n_{j}$. Supplementation with $\mathcal{N}{}^{i}{}_{j}:=n^{i}n_{j}$ allows for decomposition of type $(q,\,r)$ tensor fields, collectively denoted $\mathcal{T}{}^{q}_{r}(\Sigma)$, into intrinsic and normal parts. Introduce a smooth vector field $\rho^{i}\in\mathfrak{X}(\Sigma)$ satisfying $\rho^{i}{\nabla}{}_{i}[\rho]=1$. Then $\rho^{i}=Nn^{i}+N^{i}$ where $N^{i}=\mathcal{P}{}^{i}{}_{j}\rho^{j}$ and consequently the ambient metric may be decomposed via:

[TABLE]

where capital Latin indices take values in $\{2,\,3\}$ and here $\delta$ is the Kronecker delta.

3.1. Intrinsic $2$-geometry and spin-weight

In order to further adapt the intrinsic $\mathbb{S}{}^{2}_{\rho}$ part of the geometry we take the view of [8, 7, 29]. Without going into too much detail we mention that the $\eth$-formalism is based on the fundamental relationships between the unit 2-sphere $\mathbb{S}^{2}$, its frame bundle and the group $\mathrm{SO}(3)$444Strictly speaking its simply connected cover $SU(2)$ is more fundamental because it allows us to also describe spinorial quantities but for the present purpose it is enough to consider the vectorial aspects related to the rotation group $\mathrm{SO}(3)$.. In what follows we regard the 2-sphere as the unit-sphere equipped with the usual Euclidean metric. The bundle of frames over $\mathbb{S}^{2}$ is diffeomorphic to the rotation group since every rotation matrix consists of three orthonormal vectors which form an oriented basis of $\mathbb{R}^{3}$. Interpreting the first vector as a point $\mathbf{e}$ on $\mathbb{S}^{2}$, the other two vectors yield an orthonormal basis in the tangent space $T_{\mathbf{e}}\mathbb{S}^{2}$. Keeping $\mathbf{e}$ fixed we see that all frames at $\mathbf{e}$ are related by a 2-dimensional rotation, i.e., an element of $\mathrm{SO}(2)$. It is easily seen that this correspondence between frames on the 2-sphere and a rotation matrix is bijective and that it allows us to regard the $2$-sphere as the factor space $\mathrm{SO}(3)/{\mathrm{SO}(2)}$. The projection map $\pi:\mathrm{SO}(3)\rightarrow\mathbb{S}^{2}$ is called the Hopf map.

Every tensor field defined at a point $\mathbf{e}\in\mathbb{S}^{2}$ can be decomposed into components with respect to a basis in $T_{\mathbf{e}}\mathbb{S}^{2}$ and we may regard these components as functions defined at a particular point on $\mathrm{SO}(3)$. Since they are components of a tensor field they change in a very characteristic way when we change the basis in $T_{\mathbf{e}}\mathbb{S}^{2}$. In this way we can describe every tensor field on the sphere by a set of functions with special behaviour under change of basis. By regarding this set as a whole we have eliminated the need for referring to a particular choice of basis on the 2-sphere. This is the main advantage in this formalism since it is well known that there are no globally well defined frames on $\mathbb{S}^{2}$ — a fact, which creates many problems for numerical simulations involving the 2-sphere.

Next, we introduce appropriate Euler angles for rotations and polar coordinates on the 2-sphere so that we can express these well defined global relationships in local coordinates. The Hopf map can be expressed in these coordinates as $\pi:(\theta,\,\psi,\,\phi)\mapsto(\vartheta,\,\varphi)=(\theta,\,\phi)$.

Consider the open subset $U\subset\mathbb{S}^{2}$ away from the poles ($\vartheta=0,\,\pi$) such that the Hopf map with respect to the given coordinates is well-defined. A smooth (real) orthonormal frame $\hat{\mathbf{e}}_{(I)}$ on $U$ may be introduced where the parentheses indicate distinct frame fields. Define the complex field \mathbf{m}:=\left(\hat{\mathbf{e}}_{(2)}+i\hat{\mathbf{e}}_{(3)}\right)/\mathchoice{{\hbox{\displaystyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\textstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\scriptstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=4.51111pt,depth=-3.6089pt}}}{{\hbox{\scriptscriptstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}. In terms of this complex linear combination we can express the action of $\mathrm{SO}(2)$ as multiplication with a phase $\mathbf{m}\mapsto\mathbf{m}^{\prime}=e^{i\psi}\mathbf{m}$ inducing rotation of the complex frame $(\mathbf{m},\,\overline{\mathbf{m}})$ together with its dual coframe ($\boldsymbol{\omega},\,\overline{\boldsymbol{\omega}})$ which leads to the notion of spin-weight [8]. Given a smooth tensor field $T\in\mathcal{T}{}^{q}_{r}(U)$ its (equivalent) spin-weighted representation is provided by:

[TABLE]

where $s=r_{1}-r_{2}-q_{1}+q_{2}$ is the spin-weight.

The unit-sphere metric $\overset{\circ}{\gamma}$ on $\mathbb{S}^{2}$ when expressed in these coordinates acquires the form:

[TABLE]

where the choice of coordinatisation555This selection is made here to align with the numerical scheme we utilise in §4. An equivalent construction may be performed in (for example) complex stereographic coordinates [43, 76, 42, 71]. entails that orthonormal frame vectors may be selected as $\hat{\mathbf{e}}_{(2)}=\partial_{\vartheta}$ and $\hat{\mathbf{e}}_{(3)}=\csc\vartheta\,\partial_{\varphi}$. The associated complex reference (co)frame becomes:

[TABLE]

subject to the complex orthonormality conditions666To avoid later confusion we shall keep $\mathbf{m}$ and $\boldsymbol{\omega}$ distinct as $\overset{\circ}{\gamma}$ shall be demoted from the status of a metric.:

[TABLE]

On account of Eqs. 3.3, 3.4 and 3.5 we thus have:

[TABLE]

Denote the Levi-Civita connection associated with $\overset{\circ}{\gamma}$ by $\overset{\circ}{\mathfrak{D}}$. We will now use this to define derivative operators which are adapted to the notion of spin-weight, mapping spin-weighted quantities to spin-weighted quantities. Suppose ${}_{s}f$ is a spin-weighted quantity in the sense of Eq.(3.2). We define the $\eth$ operators as components of the corresponding tensor field in the direction of $\mathbf{m}$ and $\overline{\mathbf{m}}$ as follows:

[TABLE]

where $\overset{\circ}{\Gamma}=m^{I}\overline{m}^{J}\overset{\circ}{\mathfrak{D}}_{J}\left[\omega_{I}\right]$ and for the choice of Eqs. 3.3 and 3.4 we have \overset{\circ}{\Gamma}=\cot\vartheta/\mathchoice{{\hbox{\displaystyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\textstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\scriptstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=4.51111pt,depth=-3.6089pt}}}{{\hbox{\scriptscriptstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}} together with:

[TABLE]

Explicit translation formulae for covariant derivatives may be arrived at by virtue of Eq.(3.7) (see appendix of [76], but note conventions differ):

[TABLE]

Finally, we note that if the tensor field $W$ is real then under complex conjugation ${}_{+s}W^{*}={}_{-s}W$ and furthermore the operator actions $\eth\leftrightarrow\overline{\eth}$ are interchanged.

3.2. Topological $2$-spheres

In order to relax our treatment to more general geometries the assumption of §3.1 shall be modified and instead we shall consider working with a manifold $\mathbb{S}^{2}_{\rho}$ which is diffeomorphic to $\mathbb{S}^{2}$ but is equipped with a different metric. The approach we follow is based on [43, 76] and hence we shall only briefly summarise the idea here.

Consider the manifold $(\mathbb{S}{}^{2}_{\rho},\,\gamma{}_{IJ})$ endowed with metric:

[TABLE]

where the coframe is that of Eq.(3.4) and the expression follows from consideration of the irreducible decomposition of a type $(0,\,2)$ tensor field [43, 71]. The $\overset{\circ}{\gamma}$ of Eq.(3.3) is now demoted to the status of an auxiliary field. Thus, to be explicit, while the conditions of Eq.(3.5) continue to hold, indicial manipulations of tensorial quantities are now to be performed with $\gamma{}_{IJ}$. The inverse metric is given by:

[TABLE]

where [43]:

[TABLE]

Note that:

[TABLE]

If coordinate components $(\vartheta,\,\varphi)$ are specified as:

[TABLE]

then comparison with Eq.(3.4) and Eq.(3.10) yields:

[TABLE]

and $\Gamma_{i}\in\mathbb{R}$ implies ${}_{-2}\gamma={}_{+2}\gamma^{*}$. In order to complete the ingredients for a manifestly regular treatment of quantities in $\mathcal{T}{}^{q}_{r}(\mathbb{S}^{2}_{\rho})$ we require description of the covariant derivative operator in this new setting. Denote the Levi-Civita connection associated with $(\mathbb{S}{}^{2}_{\rho},\,\gamma_{IJ})$ by ${\mathfrak{D}}$. Then ${\mathfrak{D}}$ and $\overset{\circ}{\mathfrak{D}}$ may be uniquely related by introduction of a $(1,\,2)$ tensor field $C{}^{K}{}_{IJ}$ [84]:

[TABLE]

The tensor field $C^{K}{}_{IJ}$ arises as the difference between Christoffel symbols associated with each connection and consequently in evaluating the action of ${\mathfrak{D}}$ on a given field the pattern of additional terms matches that of the usually required $\Gamma{}^{K}{}_{IJ}$. For example, let $V\in\mathcal{T}{}^{0}_{1}(\mathbb{S}^{2}_{\rho})$ and $T\in\mathcal{T}{}^{0}_{2}(\mathbb{S}^{2}_{\rho})$ then:

[TABLE]

Furthermore, $C{}^{K}{}_{IJ}=\gamma{}^{KL}C{}_{LIJ}$ may itself be described in terms of spin-weighted components of the metric $\gamma$ and $\eth$ derivatives thereof via projection exploiting Eq.(3.2) and Eq.(3.9):

[TABLE]

Thus translation formulae for construction of manifestly regular expressions (in the sense of coordinates) may also be derived in the present context for description of Eq.(3.17) and more general geometric quantities. In particular, see [43] for explicit calculations and expressions involving the scalar curvature $\mathcal{R}[\gamma]$.

3.3. $\Sigma$ decomposition

A frame formalism based description of geometric quantities exploiting the preferred selections made in §3.1 and §3.2 may now be constructed as follows. An element of the foliation of $\Sigma$ by $\mathbb{S}{}^{2}_{\rho}$ is fixed by selecting some $\rho_{0}$ wherein local coordinates $x^{I}=(\vartheta,\,\varphi)$ may be chosen. These coordinates may then be Lie dragged along the integral curves of $\rho^{i}$ to other leaves of the foliation [76, 6] resulting in $x^{i}=(\rho,\,\vartheta,\,\varphi)$. The preferred orthonormal complex (co)frame introduced in Eq.(3.4) is extended as $\omega_{i}:=\delta^{I}_{i}\omega_{I}$ and $m^{i}:=\delta^{i}_{I}m^{I}$ and further supplemented with $\hat{\mathbf{e}}_{(1)}=\partial_{\rho}$ which allows for spin-weighted decomposition of fields in $\mathcal{T}{}^{q}_{r}(\Sigma)$.

We briefly demonstrate how this pieces together in decomposition of the ambient metric of $(\Sigma,\,g_{ij})$. The normalisation condition on $n_{i}$ together with the fact that $N^{i}$ is an intrinsic vector leads to:

[TABLE]

which may be written by virtue of Eqs. 3.2, 3.4 and 3.5 as:

[TABLE]

where we have set ${}_{0}N:=N$. Similarly,

[TABLE]

where:

[TABLE]

We expand $g_{ij}$ (or indeed any covariant symmetric tensor field) with respect to the coframe as:

[TABLE]

and with Eq.(3.1) it is found that:

[TABLE]

A similar, though more laborious approach of projection may be used to find explicit decompositions for more general elements of $\mathcal{T}{}^{q}_{r}(\Sigma)$. In particular, Eq.(3.17) and Eq.(3.18) lead to a representation of the action of the ambient Levi-Civita connection $\nabla:\mathcal{T}{}^{q}_{r}(\Sigma)\rightarrow\mathcal{T}{}^{q}_{r+1}(\Sigma)$ in terms of the (complex) frame and thus spin-weighted components together with terms involving the extrinsic curvature $K_{IJ}:=\frac{1}{2}\pounds_{\mathbf{n}}[\gamma_{IJ}]$. Consequently manifestly regular, frame representations of the formal, linear adjoint $L_{\overline{g}}^{*}[\cdot]_{ij}$ appearing in Eq.(2.5) and indeed all related, required quantities for the scalar curvature deformation problem may be constructed [27].

4. Numerical method

In considering the numerical solution of the deformation problem described in §2 as adapted to the frame-formalism of §3 we exploit (pseudo)-spectral (PS) methods [81, 14, 49, 11] as they give rise to highly efficient techniques for solution approximation as the differentiability class of functions increases. In brief, the idea is that given a square-integrable function $f\in L^{2}(\Omega)$, global approximation of $f$ over $\Omega$ is made by truncating a representation of $f$ in terms of a suitably chosen complete orthonormal basis $(\Phi_{n}){}_{n=0}^{\infty}$ of $L^{2}(\Omega)$ as $\tilde{f}:=\sum_{n=0}^{N}f_{n}\Phi_{n}$. Numerical derivatives may thus be evaluated directly through their action on basis functions or embedded via recursion relations involving the expansion coefficients $(f_{n}){}_{n=0}^{N}$. The details of how the approximation is enforced and hence how the aforementioned coefficients are to be selected is controlled through a choice of test functions $(\Psi_{n}){}_{n=0}^{N}$ and the inner product associated with the natural function space for $f$.

4.1. Function approximation on $\mathbb{S}{}^{2}_{\rho}$

In order to numerically treat functions over $\Sigma$ we begin by considering (as in §3.2) the submanifold with metric $(\mathbb{S}^{2}_{\rho^{*}},\,\gamma_{IJ})$ of $\Sigma$ where $\rho^{*}$ has been fixed. We shall assume square integrability with respect to the measure induced by $\overset{\circ}{\gamma}$ (Eq.(3.3)) for sufficiently regular scalar fields ${}_{0}f:=f$ or more generally, upon projection via Eq.(3.2) and Eq.(3.4), spin-weighted components ${}_{s}f$ of tensor fields $f\in\mathcal{T}{}^{q}_{r}(\mathbb{S}^{2}_{\rho^{*}})$. Leveraging the well-known spin-weighted spherical harmonics (SWSH) $({}_{s}Y_{lm})_{l,m}$ allows us to write [42, 71]:

[TABLE]

which converges in the $L^{2}$ sense described in [8]. The band-limit $L$ appearing in Eq.(4.1) is fixed at some finite value to provide a truncated approximation ${}_{s}\tilde{f}$. On account of the SWSH orthonormality relation (note commensurate $s$) [42, 71]:

[TABLE]

an invertible map $\mathcal{F}:{}_{s}f\mapsto{}_{s}f_{lm}$ may be constructed allowing one to transform between nodal (sampled function) and modal (coefficient) descriptions of an approximated function. Due to {}_{0}Y_{00}(\vartheta,\,\varphi)=1/\mathchoice{{\hbox{\displaystyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\textstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\scriptstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=4.51111pt,depth=-3.6089pt}}}{{\hbox{\scriptscriptstyle\sqrt{2,}}\lower 0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}} [8], the relation of Eq.(4.2) together with $\mathcal{F}$ also may be viewed as supplying a general quadrature rule for functions of spin-weight [math].

During the course of our numerical work, the fast Fourier transformation (FFT) based algorithm of [8] is used to compute $\mathcal{F}$ and its inverse with an overall algorithmic complexity of $\mathcal{O}(L^{3})$. In this approach the function ${}_{s}f(\vartheta,\,\varphi)$ is sampled on a finite, product grid in $\vartheta$ and $\varphi$ of uniform spacing, and subsequently, data is periodically extended – the details of which are controlled by the value of $s$ and the choice of $L$. With a view towards later imposition of axisymmetry in §5.3 note that if the $\varphi$ dependence appearing in ${}_{s}f(\vartheta,\,\varphi)$ is trivial777Or equivalently only the $m=0$ mode need be considered in Eq.(4.1) and the angular dependence of the SWSH becomes trivial [42] motivating the definition ${}_{s}Y_{l}(\vartheta):={}_{s}Y_{l0}(\vartheta,\,\varphi)$. then algorithmic complexity may be further improved in accordance with the usual $\mathcal{O}(L\log L)$ scaling associated with a one-dimensional FFT [23]. While it is straightforward to modify the SWSH transformation algorithm of [8] such that the sampling in the $\varphi$ direction is reduced or varied adaptively, a few points888With the $\varphi$ sampling of [8] the sum over $m$ in Eq.(4.1) may be reduced to $|m|\leq\min(L_{\varphi},\,l)$ where $L_{\varphi}:=|s|$. must be sampled on account of various auxiliary quantities appearing in the implementation.

The transformation gives us the freedom to describe functions in two ways. This freedom is crucial for our scheme insofar as the $\eth$ and $\overline{\eth}$ operators introduced in §3.1 when evaluated on SWSH reduce to an algebraic action [8] (see also Eq.(4.15.122) of [71]):

[TABLE]

which in turn allows for numerical derivative calculation to be embedded in the modal representation of a numerically sampled spin-weighted function. Consequently, by making use of the SWSH in a truncated approximation ${}_{s}\tilde{f}(\vartheta,\,\varphi)$ based on Eq.(4.1) the adapted action of Eq.(4.3) shunts away any issues that may have arisen due to apparent singularities introduced by our choice of coordinatisation. Given two spin-weighted functions ${}_{s_{1}}f$ and ${}_{s_{2}}g$ where $s_{1}$ and $s_{2}$ may be distinct the SWSH transformation also allows for decomposition of the nodal, point-wise product ${}_{s_{1}}f{}_{s_{2}}g$ in terms of a linear combination of SWSH with spin-weight $s=s_{1}+s_{2}$ [8].

A final remark with respect to efficiency mentioned in the section introduction when working with spin-weighted functions is in order. Define the averaged coefficient $F_{l}:=\langle{}_{s}f_{lm}\rangle_{m}$ where ${}_{s}f_{lm}$ is as in Eq.(4.1). Suppose ${}_{s}f(\vartheta,\,\varphi)$ is smooth, then there exist $\mathcal{A},\,\mathcal{B}>0$ such that for sufficiently large $l$ we have $|F_{l}|\sim\mathcal{A}\exp(-\mathcal{B}l)$ [56, 11].

4.2. Function approximation on $\Sigma\simeq\mathbb{R}\times\mathbb{S}{}^{2}_{\rho}$

Since our goal is the numerical solution of the deformation problem outlined in §2 and formulated with respect to $\Omega\Subset\Sigma$ we shall now focus on the domain $\Omega:=[\rho_{\mathrm{min}},\,\rho_{\mathrm{max}}]\times\mathbb{S}^{2}_{\rho}$. In order to approximate fields in $\mathcal{T}{}^{q}_{r}(\Omega)$ the expansion of Eq.(4.1) is generalised by allowing $\rho^{*}$ to vary over the closed interval $[\rho_{\mathrm{min}},\,\rho_{\mathrm{max}}]$ such that the expansion coefficients acquire an additional univariate $\rho$ dependence for each fixed $l$ (and $m$ if axisymmetry is not assumed). As a preliminary we map the interval $[\rho_{\mathrm{min}},\,\rho_{\mathrm{max}}]$ to the standard interval $[-1,1]$ with coordinate $\nu$ via:

[TABLE]

If a function $f$ is continuous and either of bounded-variation or satisfies a Dini-Lipschitz condition on $[-1,\,1]$ then the Chebyshev series converges uniformly [63, 11]:

[TABLE]

Note that in Eq.(4.5) a factor of $1/2$ has been absorbed into $f_{0}$:

[TABLE]

where $\delta_{0n}=1$ if $n=0$ and is [math] otherwise. The rate of convergence of the truncated approximant $f_{N}$ is controlled by function differentiability. For $f\in C^{m+1}([-1,\,1])$ the bound $\left|f(\nu)-f_{N}(\nu)\right|=\mathcal{O}(N^{-m})$ for all $\nu\in[-1,\,1]$ holds [63]. Combining Eq.(4.1) and Eq.(4.5) allows for the spin-weighted representation of $X\in\mathcal{T}{}^{q}_{r}(\Omega)$ to be approximated as:

[TABLE]

where we have allowed for the possibility of adaptivity in $\varphi$ (indeed axisymmetry) when using [8]. As we shall actually restrict to axisymmetry in §5.3 it is convenient to further rewrite Eq.(4.7) as:

[TABLE]

where ${}_{s}Y_{l}(\vartheta)$ are real functions [8] and this expansion is to be understood as implicitly evaluated via Eq.(4.7).

4.3. Complex analytic tools

In principle we now have the ingredients required to turn Eq.(2.7) subject to the ansatz on $u$ offered by Eq.(2.9) into a numerical, linear-algebraic problem. However, imposing $u=\mathcal{N}\tilde{u}\omega^{\beta}$ requires various ratios of weight function terms (and their derivatives) to be computed. Additionally, a method is required for accurate determination of \mathcal{N}\delta/\mathchoice{{\hbox{\displaystyle\sqrt{\omega,}}\lower 0.4pt\hbox{\vrule height=4.30554pt,depth=-3.44446pt}}}{{\hbox{\textstyle\sqrt{\omega,}}\lower 0.4pt\hbox{\vrule height=4.30554pt,depth=-3.44446pt}}}{{\hbox{\scriptstyle\sqrt{\omega,}}\lower 0.4pt\hbox{\vrule height=3.01389pt,depth=-2.41113pt}}}{{\hbox{\scriptscriptstyle\sqrt{\omega,}}\lower 0.4pt\hbox{\vrule height=2.15277pt,depth=-1.72223pt}}} when only the numerical result of the product $\delta=f\omega$ is known. Recall that $\omega\rightarrow 0$ as $\partial\Omega$ is approached (see §2.1). Consequently division of two quantities that vanish towards $\partial\Omega$ in a manner that is known from analytical results to yield a quotient of well-defined (finite) value must be computed using only numerical data. A further issue occurs in that high-order derivatives (up to fourth order in Eq.(2.8) for example) are required which is known to be ill-posed999Ill-posed in the sense that small perturbations in the function to be differentiated may lead to large errors in the differentiated result [57, 64, 77]. when formulae are restricted to finite precision calculations with real arithmetic.

The concern of derivative accuracy for analytic functions may be mitigated by transformation to integrals in the complex plane through the use of the Cauchy representation formula (CRF) [61, 33, 9, 36]. This approach also potentially provides a solution to the division problem. To concretise the idea suppose that the real function $f$ possesses a complex analytic extension such that $f:U\rightarrow\mathbb{C}$ is holomorphic on an open set $U\subset\mathbb{C}$. Suppose the closed disc of radius $R$ satisfies $\overline{D}_{R}\subset U$. Recall the CRF allows for the value of $f$ (and complex derivatives thereof) to be calculated at a base point $z_{0}\in D_{R}$ through integration over a piecewise $C^{1}$ closed curve $\Gamma$ equipped with counter-clockwise orientation in $U\setminus\{z_{0}\}$ that can be continuously deformed in $U\setminus\{z_{0}\}$ to $\partial D_{R}$ [46]. In particular, if $\Gamma$ circumscribes a base point $z_{0}$ on the real line then $f$ may also be described implicitly at $z_{0}$ without recourse to direct sampling at the point. Immediately this provides a potential mechanism to avoid the numerically unstable division process.

Pursuing the problem of derivative conditioning further, Bornemann [9] investigates stability properties of computing the Taylor series coefficients $(\tilde{f}_{n}){}_{n=0}^{\infty}$ of $f$, which when evaluated by the CRF on an origin centered, circular contour $\Gamma=\Gamma_{C}$ of radius $r>0$ take the form:

[TABLE]

Since the integrand in Eq.(4.9) is periodic and analytic its approximation via the $m$-point trapezoidal rule:

[TABLE]

converges exponentially [83] (cf. the final remark of §4.1). This approach for calculating $\tilde{f}_{n}$ was advocated for by [61] together with the identification of $r^{n}\tilde{f}_{n}$ as got from Eq.(4.10) being readily evaluated via the FFT [62]. A delicate question of how to select101010We shall assume that for fixed $n$ the value of $m$ has been selected to satisfy the Nyquist condition so as to avoid spurious aliasing. the order dependent $r(n)\in(0,\,R)$ now arises. From the perspective of the CRF any choice is valid however numerical stability degrades in the limits $r\rightarrow 0$ and $r\rightarrow R$ [9]. While an early algorithm exists for determination of $r$ by a search procedure [33, 32] it has disadvantages due to the assumption that $(r^{n}\tilde{f}_{n})_{n=0}^{m-1}$ be approximately proportional to a geometric sequence (which may not be the case generally) and a requirement for judicious selection of starting value in the search [32].

The stability issues of evaluating Eq.(4.10) when working at finite precision arise from small, finite error in evaluation of $f$ amplifying to large error in evaluation of the sum. To analyse this [9] examines both the absolute and relative error associated with calculating the coefficients $\tilde{f}_{n}$. By considering a perturbation $\hat{f}$ of $f$ within a bound of the absolute error $\epsilon$ with respect to the $L^{\infty}$ norm over the contour $\Gamma_{C}$ it is found that the normalised coefficients $r^{n}\tilde{f}_{n}$ of Eq.(4.9) and their approximations $r^{n}\tilde{f}_{n}(m,r)$ by Eq.(4.10) remain within the same $\epsilon$ bound which follows from noting that the integral and sum are both rescaled mean values of $f$. Thus normalised Taylor coefficients are well conditioned with respect to absolute error. It is the relative error of coefficients that is shown to be crucial [9]; indeed the relative condition number $\kappa$ of the CRF111111For our purposes this may be computed via Eq.(4.9) with the definition of [50]. evaluated over $\Gamma_{C}$ for each coefficient is considered and through minimisation of $\kappa$ the existence of and a method for identification of an optimal $r_{*}(n)\in(0,\,R)$ is provided. Optimal here entails selection of $r_{*}(n)$ such that round-off error is minimised during numerical work.

The recent work of [85] extends this analysis to the case of Chebyshev expansion coefficients which may be considered to be embedded as the Taylor coefficients of a particular integral transformation there described. As a preliminary, define the Bernstein ellipse $\Gamma_{E}$ with foci at $\pm 1$ and major and minor semi-axis lengths summing to the “radius” parameter $r_{B}$:

[TABLE]

where it is assumed that $r_{B}\geq 1$. Set $u=r_{B}e^{i\theta}$ then for $z\in\Gamma_{E}(r_{B})$ and $|u|\geq 1$ the relation u(z)=z+\mathchoice{{\hbox{\displaystyle\sqrt{z^{2}-1,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\textstyle\sqrt{z^{2}-1,}}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{\scriptstyle\sqrt{z^{2}-1,}}\lower 0.4pt\hbox{\vrule height=4.51111pt,depth=-3.6089pt}}}{{\hbox{\scriptscriptstyle\sqrt{z^{2}-1,}}\lower 0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}} holds [63]. The Chebyshev polynomials of the first kind of degree $n$ may be defined by [63]:

[TABLE]

Introducing $w:=e^{i\theta}$ and making use of Eq.(4.12) yields:

[TABLE]

Which motivates extension of the domain of definition for $T_{n}$ to $\Gamma_{E}(r_{B})$ as provided by [63, 85]:

[TABLE]

Suppose that $f$ is analytic on the domain interior to $\Gamma_{E}(r_{B})$, i.e., $\mathrm{int}(\Gamma_{E}(r_{B}))$. Then the Chebyshev series (Eq.(4.5)) is convergent on $\mathrm{int}(\Gamma_{E}(r_{B}))$ [82] and Chebyshev coefficients $f_{n}$ can be given a complex analytic representation [85]:

[TABLE]

Periodicity and analyticity of the integrand again allow for efficient approximation of $f_{n}$ through the $m$-point trapezoidal rule (cf. Eq.(4.10)):

[TABLE]

In complete analogy to [9] the absolute and relative stability of the above approximation are then considered by [85]. It is shown that evaluation of the Chebyshev coefficients is absolutely stable. However, the relative error depends on the $r_{B}$ selected and is controlled by the relative condition number $\kappa(\Gamma_{E}(r_{B}),\,n)$.

The determination of $r_{*}(n)$ for general $f$ which optimises the relative stability is crucial for the computation of approximations to derivatives of (the truncated approximation of) $f$ in terms of the coefficients directly which involves the evaluation of the recursion relation [11]:

[TABLE]

initialised as $f_{n}^{(0)}:=f_{n}$ ($0\leq n\leq N$) and subject to the condition $f_{N-k+2}^{(k)}=f_{N-k+1}^{(k)}=0$. Unfortunately, the search for $r_{*}(n)$ requires extensive use of asymptotic approximations to infer the condition number $\kappa$ directly [9, 85].

In order to provide a practical, numerical method for the approximate determination of $\kappa$ we propose to instead approximate $\kappa$ on some $\Gamma_{E}(r_{B})$ via the $m$-point trapezoidal rule of Eq.(4.16). The convexity of $\log(\kappa)$ together with monotonicity of $r_{*}(n)$ [9, 85] allows for us to employ the downhill simplex minimisation algorithm [75]. We initialise the search at $n=0$ with $r_{B}=1$ constructing $r_{*}(0)$. For $n\geq 1$ the search is initialised with $r_{*}(n-1)$, which once complete yields $r_{*}(n)$. Thus an order-dependent sequence $(r_{*}(n))_{n=0}$ is iteratively generated.

If the replacement $r_{B}\rightarrow r_{*}(n)$ is made in the expression for ${f}_{n}$ of Eq.(4.16) then usage of the FFT to simultaneously compute the result for all orders $n$ is precluded – thus we propose a compromise: on account of the rapid convergence of Chebyshev series for smooth functions we consider instead the average $\langle r_{*}(n)\rangle_{n\leq n_{\sigma}}$ where:

[TABLE]

and $\sigma$ is a tolerance corresponding to a normalised coefficient magnitude. As only the scaled, absolute value of ${f}_{n}$ is required in Eq.(4.18) we may determine $n_{\sigma}$ by making use of Eq.(4.16) with $r_{B}=1$ and subsequently recalculate for an improved (relative) accuracy.

5. Prototype problems and EAG for interior BBH data

We are now in a position to numerically carry out scalar curvature deformation and provide composite, scalar-flat, initial data.

Explicit expressions for the cut-off functions $\chi$ appearing in §2, together with prototype weight functions $\omega$ are required. Based on the discussion in [68], define $\chi:[0,\,1]\rightarrow[0,\,1]$ by:

[TABLE]

where in order to avoid steep numerical gradients $f(x):=x^{N}$ $(N>0)$ shall be selected here121212It is also possible to select (for example) $f(x):=\exp(-1/x)$ however this may potentially degrade numerical properties of the scheme.. This serves (approximately) the role of a cut-off function. Let $\rho\in\Omega_{\rho}:=[\rho_{\mathrm{min}},\,\rho_{\mathrm{max}}]$ and define $\Delta\rho:=\rho_{\mathrm{max}}-\rho_{\mathrm{min}}$. For later convenience, we immediately (linearly) map so as to introduce $\chi_{L}:\Omega_{\rho}\to\mathbb{R}$ growing from [math] to $1$ over $[\rho_{\mathrm{min}},\,\rho_{\mathrm{min}}+\mathfrak{f}\Delta\rho]$ and, similarly, $\chi_{R}$ decaying from $1$ to [math] over $[\rho_{\mathrm{max}}-\mathfrak{f}\Delta\rho,\,\rho_{\mathrm{max}}]$ where $\mathfrak{f}>0$. Consequently, we may model a univariate, normalised, weight function $\hat{\omega}_{C}$ through:

[TABLE]

Thus, explicit selection of $N$, together with $\mathfrak{f}$ allows for implicit control on the behaviour of the potential $u$ in the vicinity of $\partial\Omega_{\rho}$ when Eq.(2.7) is solved numerically.

In order to close the details required to specify the ansatz of Eq.(2.9) we also introduce:

[TABLE]

where the prefactor choice is motivated through integration of the polynomial terms over $\rho\in\Omega_{\rho}$ so as to mitigate the dependence of the overall magnitude of $\mathcal{N}$ on the extent of $\Omega_{\rho}$. Unless otherwise stated, $\alpha=2$ will be selected in Eq.(5.3) throughout.

To demonstrate the numerical properties of our scheme we begin by considering the simpler setting of spherical symmetry in §5.1, which allows for self-consistent, convergence tests during numerical construction of the potential in §5.2 to be performed. The more physically interesting case of axisymmetry is described in §5.3 and a test problem investigated in §5.4. In §5.5 we demonstrate the gluing construction numerically in axisymmetry.

5.1. Spherical symmetry reduction

We now fix the region over which the deformation takes place as $\Omega:=[\rho_{\mathrm{min}},\,\rho_{\mathrm{max}}]\times\mathbb{S}^{2}$. Spherical symmetry is imposed via the metric ansatz:

[TABLE]

One finds that upon inserting this $\overline{g}_{ij}$ into Eq.(2.7) (i.e., the weak-formulation) together with the assumption that $u$ and $\eta$ have a univariate dependence on $\rho$ an effective, one-dimensional problem results due to angular dependence integrating out. This observation motivates formal expansion of test and trial space functions respectively through:

[TABLE]

In order to incorporate the solution ansatz of Eq.(2.9) the function families are taken to be:

[TABLE]

where $\mathcal{N}$ is defined in Eq.(5.3), $T_{n}$ is a Chebyshev polynomial and $\nu(\rho)$ is the grid mapping of Eq.(4.4).

Description of $L^{*}_{\overline{g}}[\cdot]_{ij}$ in the frame formalism (conventions of Eq.(3.23)) gives rise to coframe coefficient terms $\mathfrak{l}_{\rho\rho}$ and ${}_{s}\mathfrak{l}$ with integer $s$ satisfying $|s|\leq 2$. With the $\overline{g}_{ij}$ of Eq.(5.4) fixed during evaluation of $L^{*}_{\overline{g}}[\cdot]_{ij}$ only the coefficients131313Explicit expressions for which are provided in [27]. $\mathfrak{l}_{\rho\rho}$ and ${}_{0}\mathfrak{l}$ are non-zero. Schematically the weak formulation subject to Eq.(5.5) becomes:

[TABLE]

where $\mathfrak{m}[\cdot]$ and $\mathfrak{n}[\cdot]$ are linear functionals depending on $\overline{g}_{ij}$ and contain up to second order derivative operators in $\rho$. If the deformation is constructed based on a potential via $\delta[u]=\tilde{\delta}[u]\omega=L_{\overline{g}}\left[\omega L{}_{\overline{g}}^{*}[u]\right]$ then up to fourth order derivatives in $\rho$ are also required.

Once $A_{ij}$ and $\tilde{\delta}_{i}$ are assembled the solution coefficients ${}_{0}u_{j}$ may be extracted via standard, numerical, linear-algebraic techniques. Unfortunately the function family $(\Phi_{n})_{n=0}$ involves $\omega(\rho)^{\beta}$ with $\beta<0$ and hence some care is required in the assembly process itself so as to preserve numerical stability during the course of evaluation. Define the weighted operator:

[TABLE]

Substitution of Eq.(5.6) into $\mathfrak{m}[\Phi_{i}]\omega^{|\beta|}$ (or $\mathfrak{n}[\Phi_{i}]\omega^{|\beta|}$) appearing in $A_{ij}$ of Eq.(5.7) and expansion allows for a refactoring of expressions into products of $D{}^{n}_{\rho}[\zeta,\,\eta,\,\theta][\cdot]$ with manifestly regular functions involving the background metric coefficient terms and polynomials but excluding $\omega^{(n)}(\rho)$. Though involved, the manipulations are straightforward and provide for a mechanism to individually regularise terms containing the weight function.

During solution of the (local) nonlinear deformation problem as described in §2.2 the background metric of Eq.(5.4) is fixed and $g$ satisfying ${\mathcal{R}}[g]-{\mathcal{R}}[\overline{g}]=\delta$ for a given choice of $\delta$ is sought. The iterative scheme of §2.2 is implemented through construction of a sequence of solutions ${}^{(k)}_{\hphantom{(}0}u_{i}$ to Eq.(5.7) with $(\overline{F}(\rho),\,\overline{G}(\rho))$ fixed throughout. At each iterate the replacement $\tilde{\delta}\omega=\delta\rightarrow{}^{(k)}\delta$ is made, where ${}^{(k)}\delta$ is defined in accordance with Eq.(2.10). Corrections to the potential allow for updated metric functions to be formed through:

[TABLE]

where we have emphasised the background dependence of $\mathfrak{l}_{\rho\rho}$ and ${}_{0}\mathfrak{l}$.

During numerical construction of an update it is ${}^{(m)}\delta(\rho)={}^{(m)}\tilde{\delta}(\rho)\omega(\rho)$ that is known and hence a term of the form ${}^{(m)}\delta(\rho)\omega(\rho)^{\beta}$ with $\beta<0$ must be explicitly evaluated. This may potentially lead to numerical instability as $\rho\rightarrow\partial\Omega$ on account of the behaviour of $\omega$ in this limit. One method to alleviate this is provided in the tools of §4.3. Numerical calculation of the truncated family $\left(\partial_{\rho}^{m}\left[T(\nu(\rho))\right]\right)_{n=0}^{L_{\rho}-1}$ we continue to perform with real arithmetic based on recursion. The background metric coefficients $(\overline{F},\,\overline{G})$ however will be represented by sampling on a mapped Bernstein ellipse $\Gamma_{E}$ (see Eq.(4.11)) so as to provide a spectral representation (as in Eq.(4.16)) with radius parameter $r_{B}$ selected for each function according to the averaged, optimal radius $\langle r_{*}(n)\rangle_{n\leq n_{\sigma}}$. Derivatives of the background metric coefficients are to be prepared via the recursion relation of Eq.(4.17). Products of weight function terms appearing in $D{}^{n}_{\rho}[\zeta,\,\eta,\,\theta]$ together with polynomials shall be evaluated on $\Gamma_{E}$ with a radius parameter $r_{\circ}$ selected (uniformly for all basis function orders). During construction of terms such as the corrected metric coefficients $\left({}^{(i+1)}F(\rho(z)),\,{}^{(i+1)}G(\rho(z))\right)$ appearing in Eq.(5.9) or updated scalar curvature ${}^{[3]}{\mathcal{R}}\left[{}^{(i+1)}g\right]$ individual terms may initially be sampled on contours with distinct radii. In order to combine such terms an initial transformation to their respective modal representations is made, which allows for a subsequent, complex, nodal representation on a single, contour of commensurate radius (i.e., $r_{\circ}$) to be computed. When required, numerical quadrature is computed based on the real nodal representation of functions via the Clenshaw-Curtis rule [81] with the number of samples selected as $2L_{\rho}+2$.

5.2. Spherical symmetry: SCCT and local nonlinear deformation

We now perform self-consistent convergence tests (SCCT) on prototype problems. At the linear level, this entails selection of a background metric $\overline{g}_{ij}$, weight function $\omega$, and a “seed” potential function $u$ which allows for generation of a deformation $\delta$ analytically via Eq.(2.8). We now demonstrate that our numerical scheme is robust by showing that solution of the weak formulation yields $\tilde{u}$ which converges to $u$ as resolution is increased.

Introduce the background metric functions:

[TABLE]

the selection of which is motivated by both simplicity and construction of a prototype problem with non-constant background scalar curvature such that for $M\neq 0$ and $P\neq 0$ non-triviality of the kernel of $L_{\overline{g}}^{*}$ is avoided.

Define the seed potentials:

[TABLE]

Furthermore, we supplement the usual linear SCCT with direct specification of a target scalar curvature so as to provide prototypical scalar curvature deformation problems by introducing:

[TABLE]

where with Eq.(5.12) the target scalar curvature becomes:

[TABLE]

For convenience, remaining parameters are collected into the map:

[TABLE]

Results of numerical calculations involving a variety of numerical parameters with the complex analytic approach are provided in Fig.1. We find that while linear SCCT may be carried out with excellent accuracy the sequence of linear solutions entering the deformation problem is far more susceptible to instability. We ascribe this latter to the numerical division process involved in the calculation of ${}^{(m)}\tilde{\delta}/\omega^{|\beta|}$ where (small) local error in the vicinity of $\partial\Omega$ accumulates and is represented by spuriously populating high-order modes which in turn grow in scale with each iterate and gradually pollute low-order modes. Suppression of this is provided by filtering. While it is the case that either choice of $\beta=1$ or $\beta=-1/2$ in the ansatz on the potential $u$ appears to lead to convergence, unfortunately, as can be seen in Fig.1 (right) saturation in convergence still presents before numerical round-off.

As an alternative we pursue a hybrid scheme where terms involving $\partial_{\rho}^{n}[\omega]$ and $\partial_{\rho}^{n}[\mathcal{N}]$ that enter the factorisation of the integrand describing $A_{ij}$ in Eq.(5.7) are computed using arbitrary precision. All other quantities are calculated using standard, complex, floating-point arithmetic with the techniques of §4. An upshot of this approach is that for more generic weight functions such as:

[TABLE]

entering the deformation term $\delta$ no inconveniences due to complex analytic extensions or essential singularities arise.

We introduce further metric coefficient functions (cf. Eq.(5.10)):

[TABLE]

to which the metrics $\overline{g}_{A}$, $\overline{g}_{B}$ and $\overline{g}_{C}$ are associated. For convenience, set:

[TABLE]

Results of numerical calculations making use of the hybrid scheme are shown in Fig.2 for various test deformation problems.

As in Fig.1 we find that the results of calculations based on the hybrid approach presented in Fig.2 all lead to exponential convergence with a saturation in $\max_{j}\left|{}^{(I)}\delta(\rho_{j})\right|$ that is near numerical round-off. In contrast we find that no filtering is required and there does not appear to be much sensitivity with respect to how parameters of $\hat{\omega}_{C}$ are selected. Indeed, even with $\delta$ prepared such that $\omega_{B}$ is utilised we find convergence in the hybrid approach that does not degrade with increasing $L_{\rho}$. Due to these properties we henceforth shall only make use of this hybrid scheme and fix $\beta=-1/2$. We emphasise however that in the case of assembling $A_{ij}$ (see Eq.(5.7)) explicit refactoring of the integrand as described previously is required in order for numerical solutions to be found (linear or otherwise) based on both of the approaches investigated.

5.3. Axisymmetric deformation

Having investigated our numerical approach under the imposition of spherical symmetry in §5.1 we now turn our attention to scalar curvature deformation when the underlying metric is axisymmetric. It shall be assumed that this metric is of the form of Eq.(3.23) and that the coefficients appearing in Eq.(3.24) carry no $\varphi$ dependence. On account of the success of the mixed complex-analytic floating-point and arbitrary precision arithmetic hybrid approach a similar strategy shall be pursued here. As metric coefficients now carry a $\vartheta$ dependence that does not integrate out decompositions of fields shall be made by leveraging the SWSH functions and transformation algorithm described in §4.1.

Set $\Omega_{\rho}:=[\rho_{\mathrm{min}},\,\rho_{\mathrm{max}}]$, $\Omega_{\vartheta}:=[0,\,\pi]$ and $\Omega_{\varphi}:=\mathbb{S}^{1}$ then the full domain of interest is $\Omega_{\rho}\times\Omega_{\vartheta}\times\Omega_{\varphi}$; however, the $\varphi$ dependence is trivial and shall henceforth be suppressed. Taking the view that $\omega$ serves to impose boundary conditions by inducing decay towards $\partial\Omega_{\rho}$ on salient fields we shall continue to assume the dependence $\omega=\omega(\rho)$.

Thus, immediately we expand test and trial space functions respectively:

[TABLE]

with expansion functions of both spaces treated symmetrically:

[TABLE]

where $\Phi_{n}(\rho)$ is defined in Eq.(5.6) and ${}_{0}Y_{l}:={}_{0}Y_{l0}$ is an axisymmetric SWSH function as in §4.1. In the present context, the weak formulation of Eq.(2.7) becomes:

[TABLE]

where we have defined:

[TABLE]

and $\overline{g}(\rho,\,\vartheta)$ is the determinant of the background metric. To evaluate the internal contraction between $L_{\overline{g}}^{*}$ operators and implement a regularisation scheme analogous to that of §5.1 define the vector operator:

[TABLE]

and introduce:

[TABLE]

where the ($\omega$ independent) functions ${}_{s}\mathcal{C}_{qr}$ now complete specification of the contraction: they are the coefficients in front of all possible products of the derivatives of ${}_{0}u$. Each ${}_{s}\mathcal{C}_{qr}$ carries a spin-weight $s$ such that when combined with both $\boldsymbol{\mathcal{L}}[\cdot]$ the product has resultant spin-weight [math].

The aforementioned factoring serves an additional purpose beyond numerical regularisation in the assembly of $A_{mknl}$. As the expansions of Eq.(5.18) are truncated such that for $m$ and $n$ together there is a storage requirement of $L_{\rho}^{2}$ elements, each of which in turn requires $k$ and $l$ to be specified, the number of elements appearing in $A_{mknl}$ (ignoring symmetry) scales as $\mathcal{O}(L_{\rho}^{2}(L_{\vartheta}+1)^{2})$. Thus, if all elements are immediately constructed and sampled then naive intermediate calculations involving $A_{mknl}$ result in a storage requirement scaling as $\mathcal{O}\left(2L_{\rho}^{2}(2L_{\rho}+2)(L_{\vartheta}+1)^{2})(L_{\vartheta}+2)\right)\sim\mathcal{O}(L_{\rho}^{3}L_{\vartheta}^{3})$.

Embedding a quadrature evaluation at the intermediate stage is more efficient. To accomplish this we perform a further regrouping of individual terms in the integrand of $A_{mknl}$. On account of the tensor product basis utilised the action of the weighted operator $\tilde{\boldsymbol{\mathcal{L}}}[\cdot]:=\omega(\rho)^{1/2}\boldsymbol{\mathcal{L}}[\cdot]$ may be decoupled to $\rho$ and $\vartheta$ specific subspaces where with Eq.(5.6) and Eq.(5.19):

[TABLE]

and the individual components of $\mathcal{\boldsymbol{\mathcal{L}}}_{\rho}[\Phi_{n}]_{q}$ are of the form of $D^{m}_{\rho}\left[\frac{1}{2},\,0,\,0\right][\Phi_{n}]$ (see Eq.(5.8)). Set:

[TABLE]

then $A_{mknl}$ of Eq.(5.21) becomes:

[TABLE]

The inner quadrature ${}_{s}A_{mn}(\vartheta)$ is numerically evaluated with a Clenshaw-Curtis rule [81] whereupon expansion with the family $({}_{s}Y_{j}(\vartheta))_{j=0}^{L_{\vartheta}}$ allows for evaluation of the outer quadrature.

Linear SCCT requires evaluation of $\delta(\rho,\,\vartheta)=\tilde{\delta}(\rho,\,\vartheta)\omega(\rho)$ appearing in the integrand of $\tilde{\delta}_{nl}$ of Eq.(5.21) for a given choice of seed potential $u(\rho,\,\vartheta)=\mathcal{N}(\rho)\tilde{u}(\rho,\,\vartheta)\omega(\rho)^{-1/2}$. This we accomplish by expressing $\partial_{\rho}^{m}\left[\omega(\rho)L_{\overline{g}}^{*}[u(\rho,\,\vartheta)]\right]$ for $m=0,\,1,\,2$ via the (non-zero) spin-weighted terms $\partial_{\rho}^{m}\left[\omega(\rho)\mathfrak{l}_{\rho\rho}[u(\rho,\,\vartheta)]\right]$ and $\partial_{\rho}^{m}\left[\omega(\rho){}_{s}\mathfrak{l}[u(\rho,\,\vartheta)]\right]$ based on the decomposition technique described in §3.3. In accordance with Eq.(2.8), $\delta(\rho,\,\vartheta)$ is formed by application of the frame representation of $L_{\overline{g}}[\cdot]$. Finally, resultant terms are expanded and regrouped such that all $\omega$ containing terms are collected and represented solely via the weighted operators $D^{n}_{\rho}$ introduced in Eq.(5.8). This is possible due to the assumption of the univariate $\rho$ dependence of $\omega$.

To close this section, we provide an update rule for metric coefficient functions when represented in terms of the spin-weighted components $\left({}_{0}\overline{N},\,{}_{\pm 1}\overline{N},\,{}_{0}\overline{\gamma},\,{}_{\pm 2}\overline{\gamma}\right)$. On account of the underlying axisymmetry we may drop the distinction between $\pm|s|$. This is a consequence of the particular properties of the coordinate representation of the SWSH and the $\eth$ operators: in axisymmetry the representations of ${}_{\pm s}Y_{l0}$ are real and agree even though, abstractly, these quantities have different spin-weight and therefore lie in different spaces.

Given a sequence of potential function solutions ${}^{(k)}u_{nl}$ define the update functional:

[TABLE]

where $\mathfrak{l}$ is a general component of $L^{*}_{\overline{g}}[\cdot]_{ij}$ in the frame formalism (conventions of Eq.(3.23)). We may now write:

[TABLE]

In order to update ${}_{0}\overline{N}$ first compute:

[TABLE]

together with:

[TABLE]

Finally, based on Eq.(3.22) and Eq.(3.24) set:

[TABLE]

whereupon the positive root is taken.

5.4. Axisymmetric deformation: Test problem

On account of the restriction $\omega:=\omega(\rho)$ convergence properties in the axisymmetric case are largely controlled by the resolution selected in $\rho$. Essentially, for a sufficiently large, fixed $L_{\vartheta}$ behaviour as in §5.2 was observed. Hence, for the sake of expediency we will only provide an illustrative test problem here.

Introduce the background metric:

[TABLE]

where:

[TABLE]

We now represent $\overline{g}$ in terms of the spin-weighted components $\left({}_{0}\overline{N},\,{}_{-1}\overline{N},\,{}_{0}\overline{\gamma},\,{}_{-2}\overline{\gamma}\right)$. According to the decomposition of Eq.(3.23) and Eq.(3.24) we may immediately take ${}_{-1}\overline{N}=0$, and it follows that {}_{0}\overline{N}(\rho,\,\vartheta)=\rho\mathchoice{{\hbox{\displaystyle\sqrt{F(\vartheta),}}\lower 0.4pt\hbox{\vrule height=7.5pt,depth=-6.00003pt}}}{{\hbox{\textstyle\sqrt{F(\vartheta),}}\lower 0.4pt\hbox{\vrule height=7.5pt,depth=-6.00003pt}}}{{\hbox{\scriptstyle\sqrt{F(\vartheta),}}\lower 0.4pt\hbox{\vrule height=5.25pt,depth=-4.20003pt}}}{{\hbox{\scriptscriptstyle\sqrt{F(\vartheta),}}\lower 0.4pt\hbox{\vrule height=3.75pt,depth=-3.00002pt}}}/5. The intrinsic metric expression provided by Eq.(3.14) together with the maps of Eq.(3.15) yields:

[TABLE]

The target scalar curvature shall be defined by:

[TABLE]

where in this section we take $\delta$ as:

[TABLE]

Spin-weighted metric coefficients shall be sampled in $\Omega_{\rho}$ with mapped $\Gamma_{E}$ at fixed $r_{B}=1.6$ in order to numerically determine partial derivatives in $\rho$ based on the techniques discussed in §4.3 which are then sampled back to the real, mapped Chebyshev-Gauss-Lobatto grid (see [81] for a definition). Approximation in $\Omega_{\vartheta}$ is based on the axisymmetric SWSH algorithm discussed in §4.2.

A representative calculation for local scalar-curvature deformation is inspected in Fig.3 where geometric quantities are updated as described at the end of §5.3.

A further remark is in order: while we have selected $\overline{g}$ such that ${}_{-1}\overline{N}=0$ there is no reason a priori to inhibit non-zero updates to this quantity during the iterative construction of ${\mathcal{R}[g]}$. Indeed we find this is the case for the present example (see Fig.4 (left)). Furthermore, on account of $\left|{}^{(i)}\delta(\rho,\,\vartheta)\right|$ accumulating towards $\partial\Omega_{\vartheta}$ when $\rho\simeq\frac{1}{2}(\rho_{\mathrm{max}}-\rho_{\mathrm{min}})$ we inspect how modal representations in $\rho$ decay when averaging is performed over $\vartheta$ and vice versa for background and updated metric coefficients together with the spin-weighted contraction coefficients ${}_{s}\mathcal{C}_{mn}$ of Eq.(5.23). In order to compactly represent ${}_{s}\mathcal{C}_{mn}$ an additional average over $m$ and $n$ is taken over all coefficients of fixed $s$. Coefficient decay is displayed in Fig.4 (middle, right). It is clear that on average coefficients do not display any spurious growth – this was also verified by inspecting individual ${}_{s}\mathcal{C}_{mn}$.

5.5. Gluing: Internal binary black holes and external Schwarzschild

We finally turn our attention to a problem of physical interest, namely the gluing of binary black hole (Brill-Lindquist [13] and Misner [65]) data to an exterior asymptotic Schwarzschild end. As in [28] our approach shall be construction of initial data on a spherical shell $\Omega=\Omega_{\rho}\times\Omega_{\vartheta}\times\Omega_{\varphi}$. On the interior ball bounded by $\Omega$ for which $\rho<\rho_{\mathrm{min}}$ a vacuum constraint (at a MIT symmetry) satisfying, asymptotically Euclidean metric $g_{E}$ is prescribed whereas to the exterior of $\Omega$ where $\rho>\rho_{\mathrm{max}}$ Schwarzschild initial data are chosen. For $g_{\Omega}$ where $\rho\in\Omega_{\rho}$ we select a suitably truncated combination of these choices (see Eq.(5.41)). In contrast to [28] our numerical scheme does not follow the proposal of [41]. We rather attempt to follow the construction of Corvino’s proof [24] as closely as possible.

For convenience, recall the Euclidean metric with $\dim(\Sigma)=3$ in spherical coordinates:

[TABLE]

We can use conformal transformations so as to provide an interesting $g_{E}$ by rescaling $\delta_{\mathrm{Euc}}$ via the factor (function) $\psi$ as [6, 1]:

[TABLE]

Selection of initial data that can be interpreted as corresponding to a quantity of $\Xi$ black holes is provided by the Brill-Lindquist (BL) choice [13, 6]:

[TABLE]

where $r_{\xi}$ is the (coordinate) separation from the centre $C^{i}_{\xi}$ of the $\xi^{\mathrm{th}}$ black hole. In order to compare with [28] we work within the context of axisymmetry where symmetrically spaced, on-axis, equal mass, binary black hole data, i.e., $\Xi=2$, $m=m_{1}=m_{2}$ and (in Cartesian coordinates) $C_{1}^{i}=(0,\,0,\,d/2)=-C_{2}^{i}$ is chosen.

Free parameters appearing in $g_{E}$ are fixed as $m=2$ and $d=10$ so as to facilitate comparison with [28]. A further reason for this selection is to have a scenario where the two interior black hole horizons do not intersect and to avoid the formation of a tertiary outer horizon which is the case if the inequality $m/d\lesssim 0.64$ is satisfied (see also [13]).

External to $\Omega$ we follow [24] and select Schwarzschild initial data in isotropic form [6]:

[TABLE]

where $M_{\mathrm{ADM}}$ is the mass of the full $g$ on $\Sigma$ which is to satisfy ${\mathcal{R}}[g]=0$. The underlying axisymmetry together with invariance under $z\rightarrow-z$ for $g_{E}$ entails that $g_{S}$ need not be shifted from the origin and the only physical parameter to be adjusted for the gluing construction is $M_{\mathrm{ADM}}$.

On $\Omega$ put:

[TABLE]

where $\chi_{R}$ is the mapped cut-off function described at the start of §5. The $\omega$ entering the metric coefficient update formulae together with the weak formulation of Eq.(5.21) is selected as $\hat{\omega}_{C}$ of Eq.(5.2) with $\mathfrak{f}=6/5$ and in both $\chi_{R}$ and $\hat{\omega}_{C}$ polynomial decay with $N=4$ is chosen.

In the present context a further complication resulting from constant (i.e. zero) ${\mathcal{R}}$ exists. Specifically, if $g\rightarrow\delta_{\mathrm{Euc}}$ then as explained in §2.2 the kernel of the formal adjoint $L_{g}^{*}[\cdot]_{ij}$ becomes non-trivial and a modified, projected problem must instead be treated. In order to avoid extensive changes to our numerical scheme we propose to solve Eq.(5.20) via truncated singular-value decomposition (TSVD) [79]. For a more general choice of axisymmetric $g_{E}$ (no longer invariant under $z\rightarrow-z$) the centre of $g_{S}$ may also need adjustment during the gluing process thus we have two degrees of freedom if we view $g_{S}$ as a parameterised family of candidate solutions – accordingly all but the two smallest singular values $\sigma_{i}$ shall be retained.

In Fig.5 we display the results of a numerical calculation where the gluing constructed is implemented as previously described. In addition to determination of all updated geometric quantities we must further fix the “optimal” $M_{\mathrm{ADM}}$ in the sense that the resulting scalar curvature is minimised (and ideally [math]). This is accomplished by varying about $2m$.

In agreement with [28] we find that the mass parameter entering $g_{S}$ must satisfy $M_{\mathrm{ADM}}\geq 2m$ for the gluing construction to proceed. Qualitatively, similar behaviour is found when parameters selected for $\hat{\omega}_{C}$ and $\chi_{R}$ are modified.

In consideration of an MIT symmetry an alternative option for $g_{E}$ is possible. By making use of Eq.(5.38) and Eq.(5.39) with $\Xi=2$ in construction of $g_{E}$ we implicitly assumed a three-sheeted topology, i.e., black hole throats are disconnected and not isometric [6]. Instead, one may work with Misner data [65], which, based on the technique of spherical inversion images allows for a symmetric identification of the throats resulting in a “wormhole” within what is now a single, asymptotically flat, multiply connected manifold. For an observer external to a horizon the consequence of this topological manipulation is a modification to the interaction energy [39, 40] and hence we investigate this within the context of the gluing construction.

For concreteness, in cylindrical coordinates $(r,\varphi,z)$ the Euclidean metric takes the form:

[TABLE]

Misner data representing two equal-mass black holes aligned with $z$ and symmetrically situated about the origin is provided by [65, 6]:

[TABLE]

where $z_{n}:=\coth(n\mu)$ and $\mu$ is a free parameter which may be identified with the total mass:

[TABLE]

In fact any representative in this family of data is completely characterised by selection of $\mu$ on account of the proper length $L$ of a geodesic loop through the wormhole being [6]:

[TABLE]

To provide a direct comparison with the previous setup we solve Eq.(5.44) for $\mu$ when $\mu_{\mathrm{ADM}}=4$ using standard numerical techniques to find $\mu=1.14960525757536$. This in turn fixes $\psi_{M}$ of Eq.(5.43) which, upon mapping to spherical coordinates allows us to take $g_{E}=\psi_{M}^{4}\delta_{\mathrm{Euc}}$. The results of this numerical calculation are shown in Fig.6.

While it is the case that the new internal data reduces the required $\delta M$ and thus may be thought of as being more “efficient” we again find that the $M_{\mathrm{ADM}}$ parameterising the external Schwarzschild representative must be tuned to exceed the mass of $g_{E}$, that is, the metric on $\Omega$ tends to introduce additional energy to the gluing construction.

For both Brill-Lindquist and Misner data we numerically determined optimising masses (see Fig.5 and Fig.6 respectively) that allowed for gluing to exterior Schwarzschild to proceed at a variety of parameters. Recall that the kernel of $L_{\overline{g}_{\Omega}}^{*}[\cdot]_{ij}$ is only approximate and our TSVD procedure always discards the two smallest $\sigma_{i}$ associated with $A_{mknl}$ of Eq.(5.21) constructed based on $\overline{g}_{\Omega}$. It is thus important to inspect the full singular value spectrum of $A_{mknl}$ directly. Doing so (with values $\sigma_{i}$ ordered in descending magnitude) reveals distinct, discrete jumps in magnitude for the smallest two values however these are not particularly pronounced and as $\rho_{\mathrm{min}}$ or the extent of $[\rho_{\mathrm{min}},\,\rho_{\mathrm{max}}]$ is reduced a gradual decay is instead found. This is not unexpected for it is the case that ${\mathcal{R}}[\overline{g}_{\Omega}]\rightarrow 0$ only when $\rho\rightarrow\infty$, i.e., we are only working with an approximate kernel for $L_{\overline{g}_{\Omega}}^{*}[\cdot]_{ij}$. This may be responsible for the larger values of $\max_{\rho,\,\vartheta}\left|\mathcal{R}\left[{}^{(I)}g\right]\right|$ observed during use of smaller gluing regions. On account of this, a potential alternative approach to treat the kernel numerically may be to make use of the controlled filtering offered by a Tikhonov regularisation scheme [86, 69, 48], which we shall consider elsewhere.

6. Discussion and conclusion

In this work we have demonstrated a new numerical technique directly inspired by and based on the exterior asymptotic gluing (EAG) construction result of Corvino [24] that does not rely on a conformal Lichnerowicz-York decomposition of the constraints nor the Brill-wave ansatz approach of [41, 28, 74]. Our technique enabled fashioning of new solutions to the Einstein constraints in vacuum at a moment-in-time (MIT) symmetry based on a choice of internal Brill-Lindquist (BL) or Misner data glued over a transition region to a Schwarzschild exterior $g_{S}$. It appears that quite general asymptotically Euclidean, internal data may be glued in this sense. Unfortunately, for all calculations we performed, $M_{\mathrm{ADM}}$ of the interior set appeared as a lower bound in the sense that to construct composite initial data $(\Sigma,\,g)$ the parameter $M=M_{\mathrm{ADM}}+\delta M$ entering $g_{S}$ and enabling the gluing to proceed satisfied $\delta M\geq 0$. Thus a reduction of $M_{\mathrm{ADM}}$ based on BL internal data as claimed by [41] to be possible could not be found. This conclusion agrees with the general indications provided by the numerical results of [28, 74].

It would be of considerable interest to extend our numerical technique to incorporate a generalisation of Corvino’s result to EAG on Kerr as in [26]. For EAG on Kerr the proof technique remains similar albeit the MIT symmetry condition is relaxed. In particular, this means that the full constraint system must be considered inasmuch as the momentum constraint is no longer trivially satisfied due to the appearance of extrinsic curvature $K_{ij}$. From the point of view of numerical technique it should be feasible to employ a similar approach as in the EAG Schwarzschild scenario demonstrated here. However, clearly the system is considerably more involved. Potentially, while the technique of truncated singular value decomposition may still be feasible in treatment of the approximate kernel appearing in the adjoint linearisation of the full constraints a more geometric approach based on the Killing initial data interpretation (briefly described in §2.2) may be required.

An upshot of the increase in intricacy is a reduction in the rigidity of the possible composite $(\Sigma,\,g_{ij},\,K_{ij})$ forming initial data sets. For instance an analogous investigation to that made in this work could be based on internal Bowen-York initial data [10] and a similar question as to whether spurious gravitational wave content may be reduced could be asked. As we have not made use of conformal techniques (other than for the sake of convenience in specifying data to glue) this question is not obstructed by the results of [34, 59, 58] and may be worthwhile exploring further in this new setting.

Finally, composite data sets based on EAG would be of great interest to evolve numerically in order to better understand their dynamical properties. A potential path towards this end has been proposed in [28] where the property of an exact Schwarzschild exterior is exploited to allow for a hyperboloidal evolution scheme to proceed. We leave such investigations open to future work.

Acknowledgements

The authors are grateful for a University of Otago PhD scholarship to BD and a University of Otago Research Grant to JF.

7. Appendix

Suppose $\overline{g}\neq\delta_{\mathrm{Euc}}$ then the $L_{\overline{g}}^{*}[\cdot]_{ij}$ of Eq.(2.5) can be seen to have trivial kernel provided that $\overline{\mathcal{R}}$ is non-constant with the following formal calculation based on [31, 24]. Assume $f\in\ker(L_{\overline{g}}^{*})$ then by contraction of Eq.(2.5):

[TABLE]

Whereas taking the divergence yields:

[TABLE]

To rewrite this in terms of the scalar curvature we make use of the Bianchi identity [84]:

[TABLE]

which once contracted yields:

[TABLE]

and once further:

[TABLE]

Thus Eq.(7.2) and Eq.(7.5) show:

[TABLE]

That is, at points where $f$ does not vanish the gradient of $\overline{\mathcal{R}}$ must vanish. Consider now the behaviour of $f$ along an affinely parametrised geodesic $\gamma(s)$ with tangent vector $t^{i}$. Then directional derivatives $\overline{D}$ of $f$ along $\gamma$ are:

[TABLE]

where we made use of the geodesic equation [73]. With Eqs. 2.5 and 7.1 we find an ODE for the behaviour of $f$ along $\gamma$:

[TABLE]

Now when $f$ and $\overline{\nabla}{}_{i}[f]$ vanish at $x_{0}:=\gamma(0)$ then $f(\gamma(0))=0=f^{\prime}(\gamma(0))$ and hence by Eq.(7.8) $f(\gamma(s))=0$. It follows that $f$ must vanish in an entire neighbourhood of $x_{0}$. Due to the elliptic condition of Eq.(7.1) Aronszajn’s unique continuation theorem [2] implies that $f$ must vanish everywhere which would result in a trivial $\ker(L_{\overline{g}}^{*})$.

Suppose instead $f(x_{0})=0$ and $\overline{\nabla}_{i}[f](x_{0})\neq 0$. Then $x_{0}$ is a regular value and the zero-set of $f$ is an embedded submanifold $\mathscr{S}$ of co-dimension $1$, i.e., an embedded surface in $\Omega$ [60]. Finally, this implies that $\overline{\nabla}_{i}[\overline{\mathcal{R}}]=0$ on $\Omega\setminus\mathscr{S}$ and by continuity $\overline{\mathcal{R}}$ is constant on all $\Omega$.

Therefore when $\overline{\mathcal{R}}$ is not constant $\ker(L_{\overline{g}}^{*})$ is trivial and consequently $L_{\overline{g}}^{*}$ must be injective.