{'fraction_non_alphanumeric': 0.0364930376482723, 'fraction_numerical': 0.011057246003094379, 'mean_word_length': 4.385623819575603, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical construction of initial data for Einstein's equations with static extension to space-like infinity

We describe a numerical method to construct Cauchy data extending to space-like infinity based on Corvino's (2000) gluing method. Adopting the setting of Giulini and Holzegel (2005), we restrict ourselves here to vacuum axisymmetric spacetimes and glue a Schwarzschildean end to Brill-Lindquist data describing two non-rotating black holes. Our numerical implementation is based on pseudo-spectral methods, and we carry out extensive convergence tests to check the validity of our numerical results. We also investigate the dependence of the total ADM mass on the details of the gluing construction.

We describe a numerical method to construct Cauchy data extending to space-like infinity based on Corvino's (2000) gluing method. Adopting the setting of Giulini and Holzegel (2005), we restrict ourselves here to vacuum axisymmetric spacetimes and glue a Schwarzschildean end to Brill-Lindquist data describing two non-rotating black holes. Our numerical implementation is based on pseudo-spectral methods, and we carry out extensive convergence tests to check the validity of our numerical results. We also investigate the dependence of the total ADM mass on the details of the gluing construction.

I. INTRODUCTION
Many situations of astrophysical interest can be described to good approximation as isolated systems: an asymptotically flat spacetime containing a compact self-gravitating source such as a collapsing star, a black hole binary, etc. A fundamental problem in the numerical solution of the Einstein equations for such systems is the treatment of the far field. Access to the asymptotic region known as conformal infinity [1] is important for several reasons. Firstly, gravitational radiation is only defined in an unambiguous way at future null infinity. Including this region in the computational domain enables extraction of the gravitational radiation emitted by the source in a straightforward way. This is important for the modelling of astrophysical sources of gravitational radiation. Secondly, many open problems in mathematical relativity such as black hole stability and cosmic censorship are statements about the global structure of spacetime. If numerical studies are to shed light on these questions then access to conformal infinity is indispensable.
The standard approach to numerical relativity is based on the Cauchy formulation of Einstein's equations. The t = const slices are truncated at a finite distance from the source, where boundary conditions are imposed. These must ensure that the resulting initial-boundary value problem is well posed, they must be compatible with the constraints that hold on the individual t = const slices, and ideally they should be absorbing, i.e. the artificial boundary should be transparent to gravitational radiation. Despite much progress in this direction (see [2] for a review article), this approach is necessarily limited because exact absorbing boundary conditions cannot be defined at a finite distance in the full nonlinear theory of general relativity so that linearisation about a given background spacetime is typically assumed. And imperfect boundary conditions can easily destroy relevant features of the solutions such as late-time power-law tails caused by the backscattering of gravitational radiation.
An alternative to evolution on truncated Cauchy slices is evolution on hyperboloidal slices extending to future null infinity I + . (Examples of hyperboloidal slices are the slices Σ 1 and Σ 2 in Fig. 1.) In this approach a conformal transformation is applied to the spacetime metric, combined with a compactifying coordinate transformation that maps infinity to a finite location. The conformal boundary of the slices becomes a pure outflow boundary so that no boundary conditions are required there. Hyperboloidal evolution was first advocated in general relativity by Friedrich in the context of his regular conformal field equations [3], a symmetric hyperbolic formulation of the (suitably augmented) Einstein equations that is completely regular up to the conformal boundary. For reviews of the theoretical development as well as numerical implementations based on this system, see e.g. [1,4,5]. An alternate method is based on a straightforward ADM [6] split of the conformally transformed Einstein equations on hyperboloidal surfaces of constant mean curvature [7]. The resulting equations are formally singular at I + but can nevertheless be evaluated there in terms of regular conformal data. Based on this system, stable numerical evolutions of a gravitationally perturbed Schwarzschild black hole in axisymmetry were arXiv:1411.7878v1 [gr-qc] 28 Nov 2014 achieved [8]; later matter fields were also included [9,10]. Further proposals for hyperboloidal evolution systems that, as far as we know, have not been implemented numerically yet can be found in [11,12].
The hyperboloidal surfaces are only partial (in our case, future) Cauchy surfaces. The problem remains how to evolve entire spacetimes from Cauchy data extending to space-like infinity. The main difficulty here is that part of the Cauchy data-namely some of the components of the Weyl tensor-are singular at space-like infinity if the ADM mass is not zero [13]. In [14] Friedrich proposed a way to render these Cauchy data regular while guaranteeing the regularity of the conformal field equations at space-like infinity. The basic ingredient of this approach is the blowing up of space-like infinity i 0 to a cylinder I = [−1, 1] × S 2 that serves as a link of finite length (along the time direction) between past I − and future I + null infinity. The 2-spheres I ± = I ∩ I ± where the cylinder meets future and past null infinity are called critical sets. The equations that propagate the data from I − to I + along the cylinder acquire an extremely simple form in Friedrich's representation that makes them ideal for numerical implementation, see [15][16][17][18][19] for some recent numerical work. On the cylinder all the spatial derivatives drop out. Therefore, the cylinder is a total characteristic of the system and hence no boundary conditions are required there. However, this intrinsic system of propagation equations degenerates at the critical sets I ± and develops logarithmic singularities there that are expected to travel along null infinity and spoil its smoothness. In Friedrich's approach this generic singular behaviour is successfully reproduced. Its appearance has been made explicit and related to the structure of the initial data. In other words, there is a possibility that by choosing appropriately the initial data the occurrence of non-smooth features in the solutions at null infinity can be avoided. A possible solution proposed already in [14] is to prescribe initial data that respect a set of regularity conditions involving the Cotton tensor. However it turned out [20] that these conditions are not sufficient to prevent the occurrence of the logarithmic singularities in higher order expansions of the solutions of the intrinsic system of propagation equations. In [20] Valiente Kroon proposed a new regularity condition in the form of the following conjecture: Conjecture. If an initial data set which is time symmetric and conformally flat in a neighbourhood of infinity yields a development with a smooth null infinity, then the initial data are in fact Schwarzschildean in that neighbourhood.
Recently, the results in [21,22] have pointed in favour of the conjecture, but there is still work to be done in order to fully prove it. What has been shown is that initial data that are exactly Schwarzschildean in a neighbourhood of infinity guarantee the smoothness of the solution at the critical sets. It remains to be proved that the smoothness of null infinity follows from the smoothness at the critical sets. If true, the conjecture unveils the special role that static data play in the smooth development of Cauchy data extending to space-like infinity.
One might object that initial data that are static in a neighbourhood of space-like infinity are overly restrictive. However, a powerful result by Corvino [23] suggests that this is not the case. He showed that any given asymptotically flat and conformally flat initial data can be truncated and glued along an annulus to a Schwarzschild metric in the exterior, provided the radius of the gluing annulus is sufficiently large and the mass of the exterior Schwarzschild metric is chosen appropriately. There are otherwise no additional restrictions on the metric in the interior, in particular non-static spacetimes including gravitational radiation are allowed. The method has been generalised to stationary rotating ends described by the Kerr metric, and a cosmological constant has been included [24][25][26][27].
Corvino's result can be used for the evolution problem as follows (see also [28]). Since his initial data are Schwarzschild in a neighbourhood of space-like infinity i 0 on the initial Cauchy slice Σ 0 (see Fig. 1), the future development of these initial data will also be Schwarzschild in a neighbourhood of i 0 (the shaded region in Fig. 1). By placing an artificial timelike boundary in this region, the data on Σ 0 can be evolved to the future for some time using standard Cauchy evolution with exact boundary conditions taken from the known Schwarzschild solution. From this evolution, data on a hypersurface Σ 1 are obtained, e.g. a hypersurface of constant mean curvature. Outside the artificial boundary, the solution on Σ 1 is known analytically (Schwarzschild), so we obtain data on a complete hyperboloidal surface. These can then be taken as initial data for a hyperboloidal evolution code. For the problem studied in the present paper (vacuum axisymmetric spacetimes), the code developed in [8] can in principle be used.
The use of our initial data for the evolution problem. The data on the Cauchy surface Σ0 are glued to a Schwarzschild end (shaded region) and evolved to the surface Σ1 using Cauchy evolution with a finite boundary (dashed line) in the Schwarzschild region. Given the known Schwarzschild solution, the data on Σ1 can be completed to a hyperboloidal surface reaching I + . This then serves as initial data for a hyperboloidal evolution code.
The present paper deals with the first step of this proposal, namely the construction of initial data based on Corvino's gluing method. It should be stressed that the proof of Corvino's theorem is not explicit, i.e. it does not provide us with a prescription for how to actually construct the glued initial data. One of the aims of this paper is to compute such data numerically, at least in a simple setting. We assume here that spacetime is vacuum and axisymmetric. Corvino's method under these assumptions was first studied analytically by Giulini and Holzegel [29]. An important achievement of this paper was to turn Corvino's idea into an explicit PDE problem that can, in principle at least, be solved to obtain the glued data. The 3-metric at a moment of time symmetry in a vacuum axisymmetric spacetime can be written in the form of a Brill wave [30]. This comprises both the Schwarzschild solution in isotropic coordinates and, by superposition, Brill-Lindquist data [31] for an axisymmetric configuration of two non-rotating black holes (not in equilibrium). Giulini and Holzegel took the metric in the interior to be Brill-Lindquist and glued it to a Schwarzschild metric in the exterior using a general Brill wave metric on the gluing annulus. They were mainly interested in the question whether the ADM mass (i.e., the mass of the exterior Schwarzschild solution) can be smaller than the sum of the two Brill-Lindquist black hole masses, as they expected that this would reduce the (generally unwanted) gravitational radiation introduced in the gluing region. They claimed that this can be done at least to first order in the inverse gluing radius. Using numerical methods we are able to study the solution also for smaller gluing radii. This paper is organised as follows. In Sec. II we describe the details of the gluing construction and derive the equations to be solved. A novel ingredient is an integrability condition that fixes the relation between the masses of the Brill-Lindquist black holes and the exterior Schwarzschild solution (Sec. II D). Sec. III is devoted to the numerical implementation. We describe the pseudo-spectral method we use (Sec. III A) and test the code with an artificial exact solution (Sec. III B) before turning to the actual gluing problem in Sec. III C. Detailed convergence tests are carried out. Finally, we investigate how the total ADM mass depends on the details of the gluing procedure (Sec. III D). We conclude with a discussion of our results and an outlook on future work in Sec. IV.

II. THE GLUING CONSTRUCTION
Following the line of thought in [29], we set up here the mathematical framework on which our numerical study of the gluing construction in the subsequent section will be based. We will also derive an integrability condition that unveils the dependence of the ADM mass on the details of the gluing construction. Fig. 2 encapsulates the basic features of the construction proposed in [29]: the interior spacetime consists of Brill-Lindquist data, the exterior spacetime extending to space-like infinity is Schwarzschild, and the transition between the two data sets takes place along a gluing annulus which is equipped with a Brill wave metric. The gluing annulus extends from r int to r ext . Graphical two-dimensional representation of the gluing construction we are going to consider in the following. A Schwarzschildean end will be glued to Brill-Lindquist data representing two non-rotating black holes in the interior along a transition region equipped with a Brill wave metric.

A. Basic ingredients
More specifically, in the interior (r ≤ r int ) we consider axisymmetric vacuum Brill-Lindquist data [31] describing two black holes at a moment of time symmetry, where δ = dr 2 + r 2 (dθ 2 + sin 2 θ dφ 2 ) denotes the three-dimensional Minkowski line element in spherical polar coordinates, and m k and c k , with k = 1, 2, are the masses and coordinate centres of the two black holes, respectively. In order to simplify our formulation, we will assume in the following that the two black holes are of equal mass, i.e. m 1 = m 2 = m, and that they lie symmetrically to the origin on the z-axis, i.e. c 1 = − c 2 = c = (0, 0, d 2 ), see Fig. 3. With these choices the line element (II.1) reduces to Notice that the above line element is written in conformally flat form, a feature that will play a key role in the subsequent development of the gluing construction. It can be readily confirmed that the ADM mass of the Brill-Lindquist data (II.2) is equal to 2m.
In the present work we will consider only Brill-Lindquist data where the horizons of the two black holes do not intersect. Also all cases where a third outer horizon [31], enclosing both black holes, forms-which appears when the black holes are very close to each other-would not be considered here. As shown in [31] both the above requirements are satisfied when the mass-to-distance ratio satisfies the inequality m/d 0.64. In this setting, the radius of the event horizon of each of the black holes is given by the formula [31] 3. The Brill-Lindquist data. We consider two black holes of equal mass positioned symmetrically to the origin on the z-axis. The two black holes are at a moment of time symmetry, i.e. they do not carry any spin and they are momentarily at rest relative to each other.
Therefore in the following, in order to keep the gluing annulus away from any possible horizons of the Brill-Lindquist data, the gluing radius r int will be chosen in such a way that the inequality r int > d/2 + r hor is always satisfied.
We intend to glue a Schwarzschildean end to the Brill-Lindquist data (II.2) residing in the interior of our construction. Thus, in the exterior (r ≥ r ext ) of the gluing annulus we consider the usual spherically symmetric Schwarzschild data, which when expressed in isotropic coordinates can be written in the following conformally flat form, By construction, the mass M is identical with the ADM mass of the entire glued initial data. The above two data sets (II.2) and (II.3) will be glued together using a Brill wave. This choice follows naturally from the axisymmetric nature of the Brill-Lindquist data considered in the interior of the construction. Brill waves [30] are the most general axisymmetric vaccum spacetimes with hypersurfaceorthogonal Killing vector. In spherical coordinates, the spatial metric at a moment of time symmetry is given by the Weyl-type line element The function q(r, θ) will be the unknown of our construction. It must satisfy the boundary conditions The latter condition follows from the fact that q is an even function of θ. To justify the former, one has first to write the metric (II.4) in Cartesian coordinates and inspect the behaviour of its metric coefficients on the z-axis; then the vanishing of q along the z-axis follows as a necessary regularity condition that guarantees the absence of any conical singularities on it [32]. The conformal factor ψ(r, θ) introduced above must be positive definite everywhere and must satisfy the asymptotic conditions lim r→∞ ψ = 1 at space-like infinity.
In summary, we want to construct a spacetime that is Brill-Lindquist (II.2) in the interior r ≤ r int , is of general Brill wave form (II.4) on the intermediate gluing annulus r int ≤ r ≤ r ext , and is Schwarzschild (II.3) in the exterior r ≥ r ext . In addition, all the transitions between the different regions must be smooth.

B. The recipe
The novelty of Giulini's and Holzegel's construction lies in the way they incorporated Corvino's original idea [23] solely into the definition of the conformal factor ψ, i.e. the metric on the entire threedimensional time-symmetric slice is given by the Brill wave metric (II.4) with Here β(r, θ) is the so-called gluing function, which apart of being smooth has the following properties: and all its r-derivatives must vanish at r = r int and r = r ext . An example of such a function-for simplicity we assume here θ-independence-is given in terms of the hyperbolic tangent, Let us see now how the gluing construction described in Sec. II A can be realised by the choice (II.6) of the conformal factor. Notice that the first and second term in (II.6) are of Brill-Lindquist and Schwarzschildean character, respectively. In the interior r ≤ r int the gluing function equals unity, β = 1, therefore the second term in (II.6) vanishes. Thus, the conformal factor ψ consists now only of its Brill-Lindquist part; inserting it into the Brill wave metric (II.4) and enforcing q to vanish in the interior region, the Brill wave coincides exactly with the Brill-Lindquist data (II.2). In a similar manner in the exterior r ≥ r ext only the Schwarzschildean part of ψ survives, as β = 0 there. Again inserting the resulting conformal factor in (II.4) and setting q = 0 also in the exterior region, the Brill wave (II.4) coincides with the Schwarzschildean data (II.3). In the intermediate region r int ≤ r ≤ r ext the conformal factor ψ, and consequently the function q in (II.4), have a more complicated form.
The function q(r, θ) in the gluing region will be determined by Einstein's equations. In addition to the boundary conditions (II.5) on the z-axis, smoothness requires that q and all its radial derivatives vanish at the boundaries of the gluing annulus: q = 0 and ∂ n q ∂r n = 0 at r = r int , r ext , (II.8) for all n ∈ N. The boundary conditions that q must satisfy are summarised in Fig. 4.

C. Mathematical fomulation
Having set up our gluing scheme in the previous sections, we now move on to Einstein's equations. On the initial slice these reduce to the momentum and Hamiltonian constraints. The former is identically satisfied as our data are time-symmetric, so we are left only with the Hamiltonian constraint, which in the time-symmetric case reduces to the vanishing of the Ricci scalar of the Brill wave metric (II.4), i.e.
Expanding the Ricci scalar in the above expression, the Hamiltonian constaint results in an inhomogeneous Poisson equation of the form The boundary conditions. The thick lines indicate the loci where boundary conditions on the function q must be imposed. On the z-axis the conditions (II.5) related to the axisymmetry of our construction must be satisfied while on the boundaries of the gluing annulus the conditions (II.8) must be implemented.
According to our construction in Sec. II B, the right-hand side of the above elliptic equation is given by the form of ψ in (II.6) and (II.7). Since this is fixed a priori, we will consider the right-hand side of (II.9) as an inhomogeneity and denote it by f . It should be noted that (II.9) reduces to a homogeneous Poisson equation outside the gluing annulus as the constancy of β enforces f to vanish there. Summarising, our goal in the following will be to numerically solve the second-order linear PDE (II.9) for q(r, θ) subject to the boundary conditions (II.5) and (II.8).

D. Integrability condition
At first sight it might seem that the choice of the two mass parameters m and M in the conformal factor (II.6) is unconstrained. If this were true then nothing would prevent us from gluing a Minkowskian end to the Brill-Lindquist data in the interior! This would obviously violate the positive mass theorem [33]. In fact Einstein's equations constrain the choice of the masses. One way to see this is by employing the machinery developed by Brill [30] in order to prove that the ADM mass of time-symmetric, axisymmetric, vacuum gravitational waves is positive definite. It turns out that in our setting this result can be used as a condition to determine the relation between the masses involved in our construction.
Following Brill's arguments in [30], we repeat here his original derivation adjusted to the details of our construction. Our starting point is the Poisson equation (II.9) expressed in cylindrical coordinates (ρ, φ, z): which when expressed in terms of the three-dimensional flat Laplace operator ∇ 2 in cylindrical coordinates takes the form Integrating over the interior V of a large sphere Σ of radius R centred at the origin, one gets The integration of the last term reads where in the last step we used the first of the boundary conditions (II.5) and expressed the remaining term in spherical coordinates. In the rest of the proof, the first two integrals in (II.10) will also be expressed in spherical coordinates. Inserting the result of the above integration into (II.10) and reexpressing the first and third term through the divergence theorem, one arrives at In the limit R → ∞ the last term of the above expression vanishes because q = 0 for r > r ext . In addition, according to (II.6), the conformal factor in the limit R → ∞ behaves like 1 + M 2 R ; thus the first term of the expession above reads Taking into account the last two results, (II.11) in the limit R → ∞ reduces to which is obviously positive definite. It is interesting that this expression for the ADM mass only depends on the conformal factor. Recall that the ADM mass M of our construction appears in the definition of the conformal factor (II.6) and consequently is also present into the integrand above. Based on this observation one can use (II.12) as an integrability condition for the ADM mass, namely the integral on the right-hand side of (II.12) for a specific choice of M must return the same value for the ADM mass.
In addition to M , the mass m of each of the black holes also enters the definition (II.6); thus (II.12) can be also viewed as a condition constraining the choice of these two masses. The integrability condition (II.12) will be used in the following way. After choosing the free parameters entering (II.6) and evaluating analytically the form of the expression inside the square brackets of (II.12), the integral in (II.12) will be computed numerically using the integrate sub-package of the Python SciPy library. The value of the integral computed in this way will be denoted by M I in contrast to the parameter M chosen originally.
To exemplify the use of the condition (II.12), we will use as a test case the scenario that the distance d between the black holes in the interior is taken to be zero. As we are gluing together two identical Schwarzschildean data sets in this case, one expects that the only allowed combination of the masses is M = 2 m. Indeed, the integrability condition (II.12) confirms this expectation. In Fig. 5 we have plotted the difference ∆M = M I − M between the integral (II.12) and the originally chosen value M of the ADM mass against the mass m of the individual black holes. We chose M = 4 so one expects that the integrability condition would fix the mass of the black holes to be m = 2. This is exactly the case, as depicted in Fig. 5.
We will return to the relation between M and m for general separations d of the two black holes in Sec. III D.

III. NUMERICAL IMPLEMENTATION OF THE GLUING CONSTRUCTION
In this section our numerical implementation of the gluing construction described in the previous section and some first numerical results are presented.

A. Setting up the numerical scheme
We choose to solve the Poisson equation (II.9) numerically using pseudo-spectral methods. Accordingly, the unknown function q(r, θ) is approximated by a truncated series of suitable specific polynomials. We choose to expand the r-dependence of q in Chebyshev polynomials T k and the θ-dependence in Fourier-sine series for reasons (in addition to the ones presented in [34]) that will soon become apparent.
Our two-dimensional physical domain is given by (r, θ) ∈ [r int , r ext ] × [0, π]. While the range of the angular coordinate θ is in accordance with the expansion in Fourier-sine series, the range of the radial coordinate r is not, as the Chebyshev polynomials are defined on the interval [−1, 1]. In order to map the original r-domain to [−1, 1], we use the mapping where K and L is the number of collocation points along the radial and θ-direction, respectively. Let us now turn to the boundary conditions (II.5) and (II.8). In fact this is by far the most involved part of our numerical implementation. In order to satisfy (II.8) and the first of (II.5) we make the following ansatz: whereq is an arbitrary function of its arguments and B(x) is a function of "bump" character on the gluing annulus, i.e B(x) and all its x-derivatives vanish on the boundaries of the gluing annulus. An example of a "bump" function with the above properties looks like where b 3 , b 4 are constants. The last boundary condition, i.e. the second of (II.5), is satisfied if one expands the newly introduced functionq(x, θ) in the way described in the first paragraph of this section, namelyq where K, L are as above and the constants a kl are the expansion coefficients of our series. So, by assuming that q in (III.1) is a multiple of a "bump" function B, the vanishing of q and all its xderivatives at x = ±1 is guaranteed. Similarly, by assuming that q is a multiple of sin θ, one ensures that q vanishes on the z-axis. Finally, the expansion ofq as a Fourier-sine series setsq to zero on the z-axis; consequently ∂ θ q also vanishes there as ∂ θ q = B(cos θq + sin θ ∂ θq ). The code has been written from scratch in Python.

B. Testing the code with an exact solution
Before we start using our code to study numerically the Poisson equation (II.9), we will carry out-as one should always do-some numerical tests to check the performance of our code. For this a family of exact solutions will be used. The exact solutions will be computed in the following way. First, we choose a q and compute analytically the outcome of the left-hand side of (II.9), then we equate the resulting expession with the inhomogeneity f . Now, having at hand the expression for f , one can solve numerically (II.9) for q and compare the outcome with the exact expression of q chosen originally. This procedure will give us hints about the accuracy and the convergence properties of the code.
As exact solutions we will use the following family of functions, where κ is a non-negative integer and B is the "bump" function (III.2). The main reason for choosing the above family of solutions is that it allows us to control, through the choice of κ, the differentiability, and consequently the smoothness, at x = 0. Obviously, if κ is zero or a multiple of three, then the functionq corresponding to (III.4) is a polynomial and thus infinitely differentiable C ∞ . For any other value of κ, (III.4) is finitely differentiable C j . In the following, we will assume that κ takes the values κ = 0, 7, 13, 19, 61 and as a consequence the solution (III.4) will be C ∞ , C 2 , C 4 , C 6 , C 20 at x = 0, respectively. Our goal of doing all this is not only to show that the numerical solutions converge to the exact ones, but also to observe the expected relation, see e.g. [35], between the convergence of the numerical solutions and the smoothness of the exact solution, i.e. the smoother the solution (III.4), the faster the convergence of the code. Our findings are presented in Fig. 6. Both graphs therein depict the log 10 of the absolute value of the maximum error (in other worlds the L ∞ norm) between the numerical and the corresponding exact solution for different numbers of grid points N , where here and in the following we have chosen K = L =: N . Fig. 6(a) illustrates the case of smooth functions (κ = 0). Here one observes the typical "step" behaviour of the convergence plots corresponding to polynomial functions [35]; this is because the Chebyshev polynomials form a complete basis for the polynomials, so that (III.4) is represented exactly for N > 20 (the error settles down to numerical roundoff ∼ 10 −15 ). On the other hand, Fig. 6(b) shows the case of finitely differentiable functions. It can be easily seen that in all the cases considered the numerical solutions converge to the exact ones, but with different speed. A detailed inspection of the individual plots shows that, as expected, the speed of convergence is faster the smoother our solution [35].
C. Numerical realisation of the gluing construction

Results
The results of the previous section constitute strong evidence that our code can reproduce successfully the exact solutions (III.4), and its convergence behaviour is as expected. Thus we are confident enough to proceed further in the numerical study of the gluing construction and look for general solutions of (II.9).
In order to do so, one has first to choose appropriately the free parameters entering the definition of the conformal factor (II.6) and then to compute the inhomogeneity f by evaluating the right-hand side of (II.9). Recall that according to its definition, the conformal factor depends on the mass m of the individual Brill-Lindquist black holes, their mutual distance d, the mass M of the exterior Schwarzschild region, the location of the gluing annulus r int , r ext , and the form of the gluing function (II.7). Except for a couple of conditions that constrain their choice, the above parameters can be freely chosen. The first condition follows from the fact that the gluing annulus has to be placed away from any horizons of the Brill-Lindquist data; for this the inequality r int > d/2 + r hor must always be satisfied-see Sec. II A for the details. The second condition constrains the relation of the masses m and M , as discussed in Sec. II D. Fig. 7 shows several numerical solutions of the system (II.9), (II.5), (II.8) for the following choice of the free parameters: m = 2, d = 10, r ext = 2 r int , and the ADM mass M has been chosen in accordance with the values of Fig. 11. Starting from Fig. 7(a), the distance of the gluing annulus from the origin has been gradually increased from r int = 10 to r int = 1000. As expected, the further away one places the gluing annulus, the smaller the numerically computed values of q become. This behaviour follows naturally from the fact that the Brill-Lindquist data (II.2) look more and more like Schwarzschild data the further away one goes from the origin; consequently, the Brill wave-essentially the function q-does not have "to do a lot of work" to glue the two sets of data together. Similar behaviour is observed when the distance of the annulus from the origin is kept fixed but its width is gradually increased. Now, the magnitude of q gradually decreases as it has "more and more space" to perform the gluing between the two data sets. The results of Fig. 7 are the first evidence that the gluing constructions proposed in [23,29] can be realised numerically. Whereas the analysis of [23,29] applies only to the case when the gluing annulus is placed at large distances, our numerical findings here demonstrate that these results can be extended to smaller gluing radii.
At this point, it is worth checking again what happens in the case that the distance between the two black holes is taken to be d = 0 so that there is only a single black hole of mass 2m in the centre. One would expect that as long as the condition M = 2m is satisfied, the function q must vanish; for in this setting the Brill-Lindquist data (II.2) are already in Schwarzschild form. It turns out that our code correctly reproduces the trivial solution for arbitrary position of the gluing annulus.

Convergence analysis
Let us turn now to the convergence analysis of our numerical solutions. In contrast to Sec. III B, here we do not have an exact solution to compare our numerical findings with. Thus, we have to follow a different approach to check the convergence of our numerical solutions. The usual way to proceed in such a situation is to study the decay of the expansion coefficients a kl in (III.3), see [36]. The expansion coefficients a kl must gradually decay to zero for increasingly large indices in order for the series expansion (III.3) to converge. Once q has been computed numerically, the expansion coefficients can be readily evaluated by inverting (III.3). Fig. 8 depicts the results of our convergence analysis for the numerical solution of Fig. 7(b). Because of the two-dimensional nature of the series expansion (III.3), one has to choose along which direction to study a kl . We chose here to study the convergence behaviour of the diagonal expansion coefficients a N N as they provide a good indication of the overall decay of a kl . The fall-off behaviour of |a N N | is depicted in Fig. 8(a) on a logarithmic scale; the observed approximately linear behaviour suggests an exponential decay. To make this statement more quantitative, one has to study the ratio − log 10 (|a N N |)/N in the limit N → ∞. Therefore, following [36], if the limit is a non-negative number then the expansion coefficients converge to zero exponentially. In Fig. 8(b) one clearly sees a tendency of the ratio − log 10 (|a N N |)/N to asymptote to a small positive number, which is a strong indication of exponential decay.
We will conclude the present section by presenting another indication that the numerical solutions produced in Sec. III C 1 converge exponentially. In Fig. 9 we compare numerical solutions of different resolutions to the one with the highest resolution for the solution of Fig. 7(b). Specifically, the numerical values of q for each resolution are interpolated onto the same grid and compared with the solution of highest resolution there (here an 80 × 80 grid). Finally, the L 2 -norm of the absolute value of the error for each resolution has been plotted on a logarithmic scale, see Fig. 9. The curve falls off in an approximately linear fashion.

D. Behaviour of the ADM mass
We will now investigate the dependence of the ADM mass on the details of the gluing construction. Namely, we examine if it is possible to choose the free parameters entering the definition of the conformal factor (II.6) in such a way that the ADM mass M can take values different from the sum of the two Brill-Lindquist black holes, i.e. M = 2 m. The case 0 < M < 2 m corresponds to a reduction of the ADM mass, while the case M > 2 m to an increase. In other words, we explore the possibility of gluing together the spacetimes (II.2) and (II.3) under the assumption that their asymptotic behaviour at space-like infinity (when considered separately) is different.
As already mentioned in Sec. II D, the integrability condition (II.12) can be used to study the dependence of the ADM mass on the details of the gluing construction. Depending on the choice of the free parameters, the right-hand side of (II.12), i.e. M I , can take on values that do not necessarily agree with M . In this case the integrability condition would be violated, ∆M = M I − M = 0. Here, we will only be interested in the case that M I = M holds, corresponding to a true physical solution. There are two equivalent ways to check if the integrability condition allows for a reduction (increase) of the ADM mass. In the first approach, one fixes m and looks for a value of M < 2m (M > 2m) such that M I = M ; in the second approach, one fixes M and looks for a value of m > M/2 (m < M/2) such that again M I = M . Although both approaches are completely equivalent, it turns out that the numerical study of the latter is more convenient for our purposes here. Using the second of the aforementioned approaches one can easily determine, without having to compute the exact value of M , if the integrability condition allows for an increase or decrease of the ADM mass. Therefore, in Fig. 10(a) we study the dependence of the difference ∆M = M I − M on the mass m of the individual black holes for different locations of the gluing annulus. If the violation ∆M of the integrability condition has different signs for two different values of the mass m, then according to the intermediate value theorem ∆M must vanish somewhere in between these two values of m. In Fig. 10(a), the free parameters were chosen to be M = 4, d = 10, r ext = 2 r int , and the gluing annulus has been progressively placed closer and closer to the origin-but not too close as we have to keep away from any horizons of the Brill-Lindquist data. According to the discussion above, these choices entail that the ADM mass is reduced (increased) when the integrability condition is satisfied for m > 2 (m < 2). It is clearly visible that the closer one places the gluing annulus to the origin, the smaller is the mass m for which the integrability condition is satisfied. The two upper curves of Fig. 10(a), i.e. r int = 20 and 25, lead to an increase of the ADM mass as the condition ∆M = 0 is evidently satisfied for values m < 2. Based on the data of Fig. 10(a) one cannot decide with certainty if the ADM mass is altered for the remaining r int = 50 case. To clarify this point, we have plotted in Fig. 10(b) the difference ∆M as a function of the gluing radius r int for fixed m = 2. Based on the results of Fig. 10(a), one can safely conclude that ∆M increases with m; therefore, if ∆M is positive for m = 2 then an appropriate decrease of m will cause ∆M to vanish-a setting that leads to an increase of the ADM mass. Fig. 10(b) provides strong evidence that the ADM mass is increased for any position of the gluing annulus (no matter how far out). For gluing radii larger than r int = 6000 the violation ∆M becomes of the same order of magnitude of the numerical error, i.e. 10 −11 , which indicates that it is not possible to draw any decisive conclusions about the behaviour of ∆M there. By increasing the precision one could in principle continue further out. From Fig. 10(b) it seems that the presence of the Brill wave causes the increase of the ADM mass; one expects that ∆M asymptotes to zero from positive values as the gluing annulus is progressively placed further out. In the limiting case that the gluing is performed at infinity, where the two spacetimes become indistinguishable, the Brill wave becomes trivial and ∆M vanishes. gluing annulus ∈ [6,12] gluing annulus ∈ [10,20] FIG. 11. Increase of the ADM mass as a function of the gluing radius for two choices of the gluing function: the θ-independent function β(r) (II.7) and the θ-dependent functionβ(r, θ) (III.5). For rint = 6 the increase amounts to 0.02131 ± 2 × 10 −7 and 0.0875 ± 0.0004, respectively. By just quintupling the gluing radius (from 6 to 30), these values decrease by two orders of magnitude.
Let us look a little more closely into the details of the increase of the ADM mass and try to determine it quantitatively. As already indicated in Fig. 10, the increase is larger for smaller gluing radii r int . It is only limited by our requirement that the gluing annulus lies outside the black hole horizons, r int > d/2 + r hor (see Sec. II A). (For our choice of the free parameters, i.e. m = 2, d = 10, r ext = 2 r int , the mass-to-distance ratio is m/d = 0.2 and for this the event horizon of each one of the Brill-Lindquist black holes reads r hor = 10/11, see Sec. II A for the details. Therefore, the gluing radius must always satisfy the condition r int > 65/11.) In Fig. 11 the actual increase of the ADM mass, M I − 2m, for different locations of the gluing annulus is presented. In addition to the gluing function β(r) (II.7) used before, we have also included a θ-dependence here in a form similar to the one considered in [29], β(r, θ) = β(r) + B(r) sin 2 θ It is apparent that the inclusion of the θ-dependence into the gluing function results in a further increase of the ADM mass. Notice that in both cases, the amount of increase, M I −2m, reduces extremely fast to zero with increasing gluing radius: increasing the gluing radius from r int = 6 to 30 results in a decrease of M I − 2m of two orders of magnitude. It was mentioned above that the increase of the ADM mass can be attributed to the presence of the Brill wave responsible for the gluing. To further clarify this point, we will consider the integrability condition (II. 12 Fig. 11. Thus, it seems that indeed the Brill wave is responsible for the increase of the ADM mass. peaks inside the gluing annulus, indicating that the presence of the Brill wave is mainly responsible for the increase of the ADM mass. Here the θ-independent gluing function (II.7) has been used; for the θ-dependent gluing function (III.5) the increase of the ADM mass is slightly larger, apart from that the curves look similar.
We will conclude by presenting a setup that can lead to a reduction of the ADM mass. The key point here is to find a way to reduce the contribution of the Brill wave to the ADM mass. There are three possible ways to do that. One can either increase the gluing radius (i.e. place the gluing annulus further and further out), widen the gluing annulus, or decrease the distance between the black holes. According to the findings of Fig. 10(b), the first method alone cannot lead to a reduction of the ADM mass. It turns out that an appropriate combination of the latter two does the job. Specifically, for the choice m = 2, there exists a small range of values d ∈ [3.2, 5] of the distance d between the black holes for which reduction is possible when the gluing annulus is wide enough and placed very close, but still outside, the horizons of the Brill-Lindquist data. The lower limit d ≈ 3.2 follows from the condition m/d 0.64, see Sec. II A, which prevents the appearance of a third outer horizon enclosing both black holes. Let us clarify the discussion above by giving an example of such an (admittedly quite exceptional) setup. We choose d = 4, m = 2, and r ext = 100-for this choice of the free parameters the mass-to-distance ratio m/d = 1/2 respects the aforementioned condition, the horizon of each black hole reads r hor = 4/5, and the gluing radius must always satisfy r int > 2.8. For this setup, we plot in Fig. 13 the change-increase or reduction-of the ADM mass for different positions of the gluing radius r int . Reduction is possible for a very limited range of the gluing radius, i.e. roughly 2.8 < r int < 3.15. Thus, to achieve reduction of the ADM mass, one has to place the black holes close to each other and arrange that the gluing annulus starts very close to their horizons and is wide enough. Our findings here indicate that reduction of the ADM mass, although possible, cannot be considered as a generic feature of the gluing construction.

IV. DISCUSSION
The purpose of this paper was to demonstrate for the first time how Corvino's gluing construction [23] can be implemented numerically in order to compute nontrivial Cauchy data that are Schwarzschild in a neighbourhood of space-like infinity.
Our numerical implementation is based on the analytical work by Giulini and Holzegel [29], who applied Corvino's method to axisymmetric vacuum spacetimes. In their setting, spacetime is Brill-Lindquist (II.2) out to some radius, is described by a general Brill wave (II.4) along an intermediate gluing region, and is Schwarzschild (II.3) outside this region. Einstein's equations determine the equation to be solved numerically, namely the second-order linear PDE (II.9) subject to the boundary conditions (II.5) and (II.8). In order to obtain physically meaningful solutions, one has to constrain the choice of the two mass parameters m and M appearing in the definition of the conformal factor. It turns out that Einstein's equations imply an integrability condition (II.12) that can be used for this purpose. In addition, we make sure that the gluing region lies outside of any black hole horizons.
To solve numerically the elliptic equation describing the gluing construction, we chose to use pseudospectral methods. An extensive convergence analysis, both for an artificial exact solution (Sec. III B) and for the actual gluing problem (Sec. III C 2), demonstrates the accuracy and convergence of our numerical solutions. Our results confirm the behaviour that one would intuitively expect: the numerically computed values of q decrease with the distance of the gluing annulus from the origin and with its width, see Fig. 7.
Giulini and Holzegel [29] wondered whether it is possible to choose the gluing parameters in such a way that the ADM mass M is smaller than 2m, the sum of the two Brill-Lindquist black hole masses. By reducing the ADM mass, one might hope to reduce the amount of gravitational radiation in the gluing region. Our findings in Sec. III D suggest that the presence of the Brill wave in the gluing region generically tends to increase the ADM mass. We have only been able to reduce the ADM mass in a rather special setup where the gluing region extends from close to the black hole horizons to a large distance. It should be stressed that there is a lot of freedom in the choice of the gluing function β. We tried both a θ-independent gluing function (II.7) and one with a simple θ-dependence (III.5). It could be that there exist gluing functions that lead to a more generic reduction of the ADM mass, even though we think this is unlikely.
We remark that there are other proposals for constructing Cauchy data extending to space-like infinity that are not based on Corvino's gluing method. For example, Avila [37] considered initial data that are only asymptotically static up to a given order at space-like infinity. It would be interesting to implement this approach numerically as well. Evolving such data to future null infinity is likely to be more complicated than in our approach, where spacetime is known a priori in a whole neighbourhood of space-like infinity.
Our ultimate goal is to compute an entire spacetime from the Cauchy data constructed using the methods described in this paper. As a first step, we will evolve our data to a first hyperboloidal surface reaching future null infinity; this can then be used as initial data for a hyperboloidal evolution code based on either the regular conformal field equations or the alternate approaches described in Sec. I.

V. ACKNOWLEDGMENTS
We are grateful to Carla Cederbaum, Helmut Friedrich, Domenico Giulini, Gustav Holzegel and Martín Reiris for helpful discussions. O.R.'s research is supported by a Heisenberg Fellowship and grant RI 2246/2 from the German Research Foundation (DFG).