{'fraction_non_alphanumeric': 0.031005295795048207, 'fraction_numerical': 0.014144751731317611, 'mean_word_length': 4.49589552238806, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improvements to the construction of binary black hole initial data

Construction of binary black hole initial data is a prerequisite for numerical evolutions of binary black holes. This paper reports improvements to the binary black hole initial data solver in the Spectral Einstein Code, to allow robust construction of initial data for mass-ratio above 10:1, and for dimensionless black hole spins above 0.9, while improving efficiency for lower mass-ratios and spins. We implement a more flexible domain decomposition, adaptive mesh refinement and an updated method for choosing free parameters. We also introduce a new method to control and eliminate residual linear momentum in initial data for precessing systems, and demonstrate that it eliminates gravitational mode mixing during the evolution. Finally, the new code is applied to construct initial data for hyperbolic scattering and for binaries with very small separation.


Introduction
Almost a century ago the existence of gravitational waves was first predicted [1]. Gravitational radiation offers an exciting new observational window [2,3] and the enticing possibility of multimessenger astronomy. With the second generation of gravitational wave detectors poised to come online [4,5,6], it is more important than ever to model the likely sources of gravitational waves. Some of the most promising are binary black holes, with predicted detection rates of 0.4 − 1000 per year for Advanced LIGO [7]. To detect such systems, matched filtering techniques must be used in order to extract the signal from the noise [8]. This requires accurate models of binary black hole inspiral, merger and ringdown. A vast amount of work has been done in this direction in full numerical relativity which is necessary to describe the very dynamic plunge and merger regimes (see e.g. [9,10,11,12,13] for overviews of the field). While many groups now successfully simulate binary black hole systems [14,15,16,17,18], much of the vast 7-dimensional parameter space consisting of the mass ratio q and the dimensionless spins χ A,B remains unexplored. Most of the attention has been focused on binaries close to equal mass (q 8) and modest spin (χ A,B 0.8) (although see [19,20,21,22,23,24]) For stellar mass black hole binaries, one can expect mass ratios 15 and arbitrary spin magnitudes and orientations, which leads to precession of the spins and the orbital plane. Precessing, high mass-ratio binaries have interesting dynamics, causing large modulations of the gravitational waveform. One can expect even higher mass ratios (q 30) for neutron star-black hole (NSBH) binaries (see [25] a BH-Wolf-Rayet system with BH mass 30M ). At high mass ratios, BBH systems can be used as proxies for NSBH systems(e.g. [26]). One would thus like to simulate high-mass ratio BBH systems.
Intermediate mass black holes (IMBH) with masses m = 10 2 − 10 4 M have been hypothesised to exist to complete the BBH mass hierarchy (e.g., the review [27]). Searches for IMBH have been performed and several candidates have been identified (see e.g. [28,29] for recent observations). Higher mass ratio (10 q 100) systems may serve as models for binaries containing an IMBH and a stellar mass black hole or neutron star. Advanced era gravitational wave detectors might be able to observe gravitational waves from such systems, with a detection rate of up to 10 events per year for stellar-mass -IMBH binaries [7]. It is thus important to explore these systems in numerical relativity.
The first step to numerically evolving a binary black hole spacetime is the construction of appropriate data on the initial hypersurface [30]. This involves the solution of the elliptic constraint equations with free data that corresponds to a binary in quasi-equilibrium, ideally allowing for arbitrary masses, spins and velocities of the two black holes. The Spectral Einstein Code (SpEC) [31] includes a BBH initial data solver [32] based on the extended conformal thin sandwich equations [33,34], incorporating quasi-equilibrium black hole boundary conditions [35,36,37]. This solver has been used to construct BBH for a wide range of configurations [38]. Construction of BBH with increasing mass-ratio, increasing spin magnitudes and the desire to construct initial data for highly spinning BBH with arbitrary spin axes have necessitated a variety of improvements to the initial data code compared to its original presentation [32,36,39,40,41]. This paper summarizes these improvements and extends the original code even further, in anticipation of future desire to study even more generic BBH systems.
Specifically, here, we present: (i) Flexible domain-decomposition to allow a wider range of mass-ratios, spins and separations.
(ii) Adaptive mesh-refinement to enhance computational efficiency and to ensure robust numerical convergence for mass-ratios q 5 and dimensionless spins 0.9. (iii) Improved updating formulae for iterative determination of the free parameters. These formulas allow one to achieve very high spins and mass ratios, for example an equal-mass binary with aligned spins of 0.9999, and a q = 50 single-spinning binary with spin of 0.95 on the large black hole. (iv) Building on previous work [42,43], we control of the ADM linear momentum to avoid drifts of the center of mass in BBH evolutions. This eliminates gravitational mode mixing due to the motion of the centre of mass with respect to a fixed extraction sphere.
(v) Control of the center of mass.
This paper is organized as follows. In Sec. 2 we describe in detail the numerical enhancements and additions to the code. In Sec. 3 we present the results of initial data construction for several challenging configurations as well as an exploratory evolution of a new data set that demonstrates that the control of linear momentum in initial data leads to the elimination of gravitational wave mode mixing. Finally we summarize the results in Sec. 4 and introduce the construction of initial data for closely separated binaries and binaries on hyperbolic orbits as applications of the techniques developed in this paper.

Numerical techniques
The main task of constructing initial data is twofold: first, to solve the elliptic constraint equations on the initial hypersurface; and then, to ensure that the solution represents the astrophysical situation of interest (in our case, a black-hole binary in quasi-equilibrium). In SpEC, the former is achieved by using a pseudo-spectral multidomain method; see [32]. The number of subdomains is kept fixed, but the resolution of each subdomain is dynamically adjusted to obtain low truncation error. To enforce quasi-equilibrium conditions, SpEC employs the extended conformal thin sandwich (XCTS) formalism [34]. Before solving the conformal thin-sandwich equations, various free parameters must be chosen -for example, the sizes of the excision regions, and certain other parameters that affect the location, spin or motion of the black holes. The free parameters differ from the physical parameters one desires to control, such as the masses and spins of the black holes, or the linear momentum P ADM of the initial data hypersurface. Therefore, iterative root-finding is needed, as described in Buchman et al [41]. To minimize the computational cost associated with many iterations of high resolution solves, we adopt a hybrid approach. The resolution of the domain and the free parameters are adjusted simulateneously based on the current estimated truncation error and the differences between the desired and obtained physical quantities.
In the remainder of this section, we describe in detail the improvements to the initial data code. Figure 1 indicates the geometry of the domain-decomposition employed here. There are two inner spherical shells (thick black circles labeled A and B), which are surrounded by a set of cylinders (light blue) that are aligned with the axis connecting the two black holes.

Domain decomposition
Along the axis of the cylinders there are three subdomains with rectangular crosssection (indicated in green). One of these is located between the two excision spheres, and is a truncated square pyramid. The other two are rectangular blocks. In earlier work [32] the two inner spherical shells were restricted to have the same outer radius, and all cylinders were restricted to have the same inner radius. This restriction results in a comparatively larger shell around the smaller black hole (B). For very unequal mass systems, m B m A , in particular, it may be preferable to have a smaller outer radius of shell B, roughly comparable with the sphere of influence of black hole B. This would maximize the agreement of the geometry of the domain decomposition with the structure of the solution. Therefore, here, we allow unequal radii of the two inner shells, as indicated in Fig. 1. This has the largest impact when we consider small separations in initial data (for example, for studying remnant properties) where the old domain decomposition requires a larger separation between the two black holes than the new domain-decomposition in order for the solver to converge.
The new domain-decomposition uses several parameters from which the placement and dimension of each subdomain follow unambiguously. We begin by specifying the inner and outer spherical shells: • The centres of the inner spherical shells, c A and c B , and of the outer spherical shell, c C . Note that c C is not required to lie on the line connecting c A and c B . • The inner and outer radii of the inner spherical shells and the outer spherical shell, The remaining parameters α, f cyl , f block , and f C determine the relative sizes of the cylinders and rectangular blocks: • The rectangular blocks and cylinders end on planes orthogonal to the axis connecting the centers of the excision spheres. The location of these planes is determined by the parameter α, through the requirement that these planes intersect the inner spherical shells A and B in circles of radius R A,B sin α. The opening angle of these circles as viewed from the center of the spheres is chosen to have the same value for all four planes. • The inner radii of the cylinders are determined by the parameter f cyl via Note that f cyl < 1 is required for the cylinders 1 and 3 to cover all volume outside the spheres A and B. • The size of the blocks orthogonal to the line connecting the two spheres is determined by the parameter f block , The multiplier f block must satisfy f block > f cyl to ensure that the blocks cover the entire open region within cylinders 0, 2, and 4.
• The multipler f C , which measures how much larger the outer size of the cylinders is compared to the inner edge of the outer spherical shell: To ensure complete overlap between the cylinders and the sphere C, f C > 1+C ⊥ /r C , with C ⊥ being the distance from point C to the axis of the cylinders.
The value of f block will determine the relative size of the face of the blocks to the inner spheres: If f block > 1, then the edge of the block will be entirely outside the inner spherical shell. Conversely, if f block < 1/ √ 2, then the face of the rectangular block is completely contained within the inner spherical shell. These considerations will impact which subdomain (sphere of cylinder) will provide boundary data for the blocks.
Our standard values for the grid-internal geometry coefficients are α = 45 • , f cyl = 0.95, f block = 1.05, and f C = 1.1. We have found these choices to be robust for a wide variety of component masses, spins and separations.

Adaptive mesh refinement
An important factor in efficiently generating high-accuracy initial data is the choice of resolution in each of the subdomain used in our domain decomposition (see Fig. 1). Typically, we want our representation of the solution to have about the same accuracy in all subdomains. Unfortunately, we do not know a priori what resolution is needed in a given subdomain to reach a target accuracy. Furthermore, the optimal resolution varies significantly with the physical parameters of the binary. The old initial data solver [32,36] used hard-coded resolutions, tuned to equal-mass low spin BBH. For unequal mass systems, rapidly spinning black holes, and/or widely separated binaries the old resolutions are less efficient and can even prevent convergence of the elliptic solver when a high accuracy is requested.
To generate initial data, we generally go through multiple intermediate solves, progressively improving the accuracy of the solution while converging towards the desired binary parameters. So instead of predetermining the resolution which will be used in each subdomain at each level of refinement, we can use the preceeding numerical solution to predict the optimal resolution in each subdomain to reach a target accuracy. This significantly improves the efficiency of the initial data solver, with computing times decreased by about an order of magnitude for challenging configurations. And it also allows us to push the binary parameters to more extreme values.
Our multi-domain spectral solver represents the solution in each subdomain as a tensor-product of basis-functions. Depending on the topology of the subdomain, the basis functions are Chebyshev polynomials, and/or Fourier series, and/or spherical harmonics (see [32] for details).
Following Szilágyi [44], for a given subdomain and a given basisfunction, we define the power P i in the i-th mode by the root-mean-square value of all the coefficients of the i-th mode across all spectral coefficients of the other basisfunctions. For instance, in a spherical shell with spectral expansion of the form u(r, θ, φ) = , where N θφ = (L + 1) 2 represents the number of angular coefficients ‡. ‡ For spherical harmonic basis-functions, the top two modes are filtered [32] and are therefore not included in the data P i .
For the expected spectral convergence, P i should decay exponentially as a function of i [45,32], i.e. log 10 P i when plotted vs. i should be a straight line. The slope f of this line represents the decrease in the magnitude of the spectral coefficients when going from mode i to mode i + 1. We estimate f using Eq. (53) of Szilágyi [44].The current truncation error of the spectral expansion is approximated as the highest retained coefficient [45].
Given the current estimate of the error as and the estimate of the convergence rate as f , we can reach a target accuracy * by adding modes to the spectral expansion (recall f < 0 and a higher accuracy means a lower ).The answer is rounded up so that ∆N > 1 if the current accuracy is worse than the target accuracy, and we set ∆N = 0 if < * , i.e. the resolution is not allowed to decrease. For the configuration q3 from Table 1 the resolution was allowed to decrease without noticable impact on the convergence behaviour, cf. Figure 7.
The outer spherical shell needs comparatively small angular resolution ∼ 10, and sometimes AMR yields the same resolution at neighboring E DT . Because the ADMquantities are exclusively evaluated in the outer spherical shell (cf. Sec. 2.4 below), this would result in apparent non-convergence of ADM linear and angular momentum. Therefore, we increase the angular resolution of the outer sphere by one extra gridpoint in the θ direction and the corresponding two extra grid-points in the φ direction, whenever AMR triggers an adjustment to the domain decomposition.

Iterative determination of free parameters
When constructing initial data, we wish to achive desired masses M * A , M * B and desired black hole spin vectors χ * A and χ * B . The free data, however, is instead given by the radii and angular frequencies of the apparent horizons r A,B and Ω H A,B , which we write as Therefore, one needs to determine values of the free parameters that result in the desired physical parameters. Thus we must solve the system of equations The standard approach to the problem would be to use Newton's method; however, evaluating the Jacobian J n is too expensive numerically as every evaluation of the function F requires an elliptic solve. We instead use the following approach: make an initial guess u 0 based on the Kerr expressions for both black holes, and perform an elliptic solve for F 0 . We then construct an analytic Jacobian J 0 by using Eqs. (7,8) to evaluate the partial derivatives, and update the initial guess by After this we update the Jacobian using Broyden's method [46]: where ∆u k = u k − u k−1 . This corresponds to the "secant" approximation for a function of one variable. Finally we set The major advantage of this approach lies in the use of numerical information in the update of the Jacobian. This is important in the regime where the simple analytic Jacobian becomes inadequate. Broyden's method is applied to the intrinsic physical properties of each black hole, i.e. the eight parameters listed in (5). We also control more general properties of the binary, such as the total linear momentum and the position of its centre of mass. As discussed in Sec. 2.5 this is done with explicit updating formulae that are applied simulateneously at every step of Broyden's method.
We are now faced with two intertwined iterations: AMR to tune grid-sizes to a desired truncation error; and root-finding to adjust free parameters to achieve the desired physical masses, spins, etc. When the physical parameters are still far away from the desired values, very stringent AMR resolution would waste computing time, so we aim to tighten the AMR resolution while simultaneously decreasing root-finding errors. We do so by using an overall truncation error target E DT for AMR. We start with a large value for E DT , corresponding to a small grid-size. As root-finding residuals decrease, we will decrease E DT . We proceed as follows: At iteration k = 0, 1, 2, . . ., we compute two measures of progress in root finding: First, the residual R k which quantifies how close the physical parameters are to their desired values. R k is simply the rms error in the physical parameters: Second, the improvement I k that indicates how quickly root-finding converges, defined as where We monitor 2 conditions: where E DT is the desired truncation error, and R and I are tunable parameters. The first condition assures that the resolution is increased if the root-finding convergence becomes "flat" (e.g., due to the inability to measure the masses accurately enough at the current resolution). The second condition ensures AMR resolution is sufficiently high to ensure the physical parameters can be computed more accurately than the current R k , with R being a safety factor. If either condition is satisfied and we have already reached our termination truncation error then the initial data construction is completed. Otherwise, we divide E DT by a factor of 10 and continue with the next itertion. For all cases we have encountered, the choices R = 10 2 and I = 1.5 have proven to be robust.

Calculation of asymptotic quantities
Accurate knowledge of the total energy, linear momentum and angular momentum of the constructed initial data sets aid their characterization. Even more important, accurate control of the total linear momentum is essential to avoid a drift of the center of mass of the binary during long evolutions, cf. Fig. 6. We define the linear and angular momenta on a slice Σ intersecting spatial infinity on the surface S ∞ using the Arnowitt-Deser-Misner (ADM) prescription. Our initial data satify the asymptotic gauge conditions [47] ∂γ needed to remove ambiguities in the definition of the ADM angular momentum, as well as the boundary condition g µν = η µν on S ∞ . The old code [32,48] directly evaluated the resulting surface integrals at infinity [49,47], using extrapolation in powers of 1/r to infinite radius [48]. P ADM is then found to be a combination of 1/r 2 terms of K ij , and J ADM a combination of 1/r 3 terms. The old technique, therefore, is very sensitive to small errors in K ij in the outermost sphere of our computational domain (the outer boundary is typically located at r out ∼ 10 10 M ) and particularly to the presence of constraint violating modes in that sphere. Typically, this leads to large errors in P ADM at low resolution, and large errors in J ADM even at our highest resolution. Higher accuracy can be obtained by assuming that the constraints are satisfied on our computational domain, and utilizing Gauss' law to recast the surface integrals on S ∞ as the sum of a surface integral on a sphere S 0 located at a smaller radius and a volume integral. Utilizing Ψ(S ∞ ) = 1, we write Here the normal dS j to S 0 points into the interior of V 0 (e.g. along +r if it is a coordinate sphere) and the factor Ψ 10 was inserted to eliminate terms with spatial derivatives of Ψ from Eq. (21). Using the momentum constraint in the absence of sources, the volume integral can be simplified to Here, whereΓ i jk are the connections derived from the conformal metricγ ij . Note that for conformal flatness and maximal slicing, G i = 0 and the volume integral disappears (see [50]).
In practice, for conformally curved initial data, The outer spherical shell extends to outer radius ∼ 10 10 M . Therefore, in the numerical evaluation of the volume integral in Eq. (19), the volume element associated with the outermost grid-point becomes very large and introduces numerical noise. To avoid this, we roll off the integrand G i beyond a certain radius R c , i.e. we replace G i byG i given bỹ We choose R c = 1000 max(w A , w B ), where w A,B are the widths of the Gaussians that enforce exponential falloff to conformal flatness (cf. Eqs. 45 and 46 of Lovelace et al [40]). The ADM angular momentum is also rewritten using Gauss' law as with cyclical permutations of (x, y, z) yielding the other components. For maximal slicing and conformal flatness in V 0 , Eq. (23) simplifies to Because Eq. (23) relies on the cancellation of large volume terms, it can be sensitive to errors in K ij . Accordingly, we use Eq. (24) using a surface S 0 at sufficiently large radius such that in V 0 the metric is conformally flat and K = 0.
To illustrate the importance of the transformations applied to the ADM integrals, we consider the convergence test for configuration q50. We evaluate P ADM using Eq. (15) and Eq. (19), and we evaluate J ADM using Eq. (16) and Eq. (23). Figure 2 shows the results.
The calculation of P ADM is improved by about one order of magnitude when utilizing Gauss' law, whereas J ADM improves by several orders of magnitude. We point out  Figure 2. Accuracy of the calculation of P ADM and J ADM for two different methods of evaluation. We evaluate P ADM and J ADM when truncation error E DT = 10 −n−3 is reached, and plot differences to the next lower resolution n − 1. Data shown for case q50 in Table 1.
that, in order to achieve any convergence for the old J ADM calculation, we had to manually increase the radial resolution in the outer sphere by 1 whenever the domain decomposition is adjusted. We also compute a new diagnostic, the centre-of-mass C CoM of the initial data sets using the formalism developed in Ref. [51]. In conformal flatness, the expressions from [51] reduce to where n is the outward-pointing unit normal, n = r/r. Equation (25) is numerically evaluated by expanding the conformal factor Ψ in a power-series in 1/r. We read off the (angle-dependent) coefficient of the 1/r 2 term, and expand this coefficient in spherical harmonics. Each individual spherical harmonic term can be integrated against n analytically, so that the integral (25) collapses to a linear combination of sphericalharmonic coefficients.

Control of linear momentum and centre of mass
The quasi-equilibrium conformal thin-sandwich formalism to construct binary black hole initial data was developed in a series of papers [36,37,39,40,41]. In this formalism, one chooses two excised regions (usually taken to be coordinate spheres) with centres c A,B , and solves the extended conformal thin sandwich equations [33,34] in the exterior. Boundary conditions on the excised regions ensure that they are apparent horizons, and control the spin of each black hole. The locations and the sizes of the excised regions correlate with the position and masses of the two black holes. Orbital rotation is induced by the requirement that certain time-derivatives vanish in a frame rotating with orbital velocity Ω 0 about the orign. One finally incorporates a radial expansion factorȧ 0 , which allows fine control of the orbital eccentricity [39,52,53,41]. By a suitable choice of the conformal quantities, the quasi-equilibrium approach can generate initial data with black hole spins of order 0.9998 [40]. One shortcoming of the formalism presented in [41] lies in a lack of control of the center of mass of the binary, and only incomplete control of the ADM linear momentum P ADM . The past implementations use the location of the black holes to partially control P ADM . Consider a small displacement δc applied to the centres of both excision regions. Through the orbital rotation Ω 0 about the origin, the displacement δc induces a change in velocity of the black holes of Ω 0 ×δc, with a corresponding change in P ADM . Therefore, δc could be used to cancel the components of P ADM orthogonal to Ω 0 ; however, the cross-product in Ω 0 × δc prevented any correction parallel to Ω 0 . For head-on collisions with Ω 0 = 0, no control of P ADM is possible at all. For the non-precessing simulations presented in [41], the component of P ADM parallel to Ω 0 vanishes by symmetry, and no problems arose. However, for generic precessing binaries, there will be a non-zero linear momentum orthogonal to the orbital plane, which results in a drift of the center of mass for very long simulations (see [54] for an extreme example).
Here, we propose a different means to control the full P ADM , while simultaneously allowing us to control the center of mass as well. We fix the relative separation of the centres of the excision spheres, where the separation vector D is user-specified. We use the choice of c A to control the center-of-mass C CoM of the binary. Once a first initial data set is computed (with, in general, C CoM = 0), we can update With the black-hole centres now used to control the centre of mass, we need a different means to control P ADM . We add in the outer boundary condition on the shift (Eq. (38c) of [40]) a constant velocity v 0 :  Table 1) . Shown are the magnitude of C CoM (black solid lines) and P ADM (red dashed lines) as a function of root-finding iteration k.
Here B represents the outer boundary, a sphere with radius R = 10 10 M . The velocity v 0 will effect the overall motion of the binary, and will be reflected in a corresponding change in P ADM by E ADM v 0 , where E ADM is the ADM-energy of the binary. During iterative root-finding of the free parameters, we adjust v 0 to achieve P ADM = 0 §. To motivate the updating formula for v 0 , consider a perturbation of v 0 by δv 0 , and a perturbation of c A by δc. If we allow the masses to vary, then a Newtonian-inspired formula is To summarize, relative to earlier initial-data sets, we modify the outer boundary condition for the shift by the term v 0 , cf. Eq. (28), and use updating formulae (27) and (29) to adjust c A and v 0 . Section 2.4 describes how we compute P ADM and C CoM .
We demonstrate the efficiency of the updating formulas Eqs. (27,29) in Fig. 3 that shows the magnitude of C CoM and P ADM as a function of root-finding iteration for a q = 10 precessing binary (case q10 in Table 1). The convergence is evidently very fast, with the final values of ∼ 10 −6 and ∼ 10 −8 respectively. This means that even for an inspiral lasting 10 6 M, the drift of the centre of mass due to residual linear momentum in initial data will be only ∼ 0.01M . § Using the obtained vector β i as the shift-vector in an evolution results in a translating outer boundary; this effect is eliminated by evolving with a shift vector of β i − v i 0 .

Initial data construction
We test the improvements described in the previous sections on several cases of interest, whose parameters are summarized in Table 1. The parameters were chosen to demonstrate the range of initial data sets that can be constructed with the new code and to provide some overlap with regions of parameter space which could be achieved previously.
We first illustrate the performance of the AMR outlined in Sec. 2.2 with the case q3, a configuration we will compare with the old BBH solver below. To demonstrate AMR in isolation, we fix initial data parameters, and start with target truncation error E DT = 10 −3 . We solve the constraint equations, estimate spectral truncation errors and update numerical resolution via Eq. (4). Whenever we reach the desired truncation error, we tighten the AMR error tolerances by dividing E DT by 10, until a truncation error of 10 −9 is reached. Figure 4 illustrates the behaviour of the AMR algorithm during this test. The top panel shows the total number of collocation points in the domain, which grows with each AMR iteration. The bottom panel demonstrates that the largest truncation error across all subdomains, max , closely tracks the truncation error target E DT . Figure 5 shows a convergence test of the AMR sequence shown in Fig. 4. Plotted are various quantities as a function of the effective number of grid-points N 1/3 . The top panel demonstrates the exponential decrease in the L 2 norms of the Hamiltonian and momentum constraints, which implies that this data set is constraint-satisfying. The constraints are given explicitly by: where D is the covariant derivative associated with the spatial metric. The L 2 norm is simply the normalized pointwise norm over all collocation points: The convergence of the masses and spins is shown in the middle panel. Here we plot the norms of the differences between the quantity at a given iteration and its value at the highest resolution: Once again, the convergence is essentially exponential. The bottom panel of Fig. 5, finally examines the convergence of the ADM quantities and the center-of-mass computation. Though convergence is not as clean as for the constraints, the bottom panel of Figure 5 shows that all the asymptotic quantities can be determined to better than 10 −6 . To conclude our detailed examination of the initial data set q3, we contrast the new code described here with the old code [32,41]. One of the most important upgrades lies in the control of P ADM . Figure 6 shows the components of P ADM as a function of root-finding iteration k for both the new and the old code. Both codes successfully drive P x ADM and P y ADM to zero as expected. But only the new code also drives P z ADM to zero, whereas the old code yields P z ADM 0.00138 . As we shall see in the next section, this produces non-trivial differences in the evolution.
Turning to the more challenging cases listed in Table 1, we have performed similar tests to case q3, with the free parameters fixed to their values at the end of root finding and only the resolution changing from iteration to iteration. As an example, Both codes compute P ADM in the same way (via ??), but differ in the way it is controlled.  Convergence of physical quantities with increasing resolution for case q3. Top: the norm of Hamiltonian and momentum constraints, middle: masses and spins, bottom: ADM quantities and center of mass. "∆" indicates the difference between the value of the respective quantity at the current resolution and the highest resolutions. The free parameters are fixed to their values at the end of root finding. Figure 7 shows a subset of the convergence data. This figure demonstrates that the exponential convergence shown previously for case q3 extends to all cases. In particular, the constraints are exponentially convergent. All four cases complete with a maximal resolution of less than 2.5 × 10 5 points, an improvement of a factor of 2 − 3 over the old code. ¶ Owing to the more challenging configurations, however, the constraints are 1-2 orders of magnitude larger. The physical parameters are also exponentially converegent with resolution, as illustrated in the lower panel of Fig. 7. We use ∆M B since it is frequently harder in a high mass ratio to resolve the smaller black hole, so this provides a conservative convergence test.

Root-finding
It is also important to examine the performance of the updated root-finding procedure based on Broyden's method. Figure 8 shows the root-finding results for cases q3, q10   and Spin0.9999 done with the old and new versions of the code. Note that during rootfinding the resolution of the subdomains is also allowed to change to achieve the desired truncation error. For low mass ratio both codes show similar rates of convergence and final errors. The situation changes for case q10, where the old code has trouble achieving low errors in masses and spins, while the new root-finding procedure described in Sec. 2.3 results in errors of order 10 −6 . Finally, for case Spin0.9999, the results are drastically different: the old code has errors in the masses of order several ×10 −3 and spins of order 10 −4 . Since we are attempting to construct a binary with dimensionless spins of 0.9999 it becomes clear that the old code is inadequate for this purpose. On the other hand, the new root-finding procedure successfully reduces the errors in physical quantities to the level of 10 −6 . Thus, the new algorithm allows us to achieve the desired values of the physical quantities which is especially important as we push to higher spin magnitudes. On average, the new code is about 25-50% as fast as the old one. For example, for the case q10, the old code took 12.4 hours to complete, whereas the new took 6 hours on 12 cores of a Westmere node of the Briarée compute cluster. Therefore, the new code is indeed more efficient than the old while achieving the same or better accuracy.

Exploratory evolution
We have emphasized above the importance of controlling P ADM . We now evolve initial data for case q3 constructed with the old and the new initial-data code, and compare the two evolutions in detail.
We being by considering the convergence of constraints and quasi-local quantities during the evolution. The top panel of Figure 9 shows the L 2 norm of the normalized constraint violations during the evolution (see Eq.(71) of [55]). It is obvious that both codes show similar convergence properties, as expected. Further, the initial spike of constraint violations due to junk radiation is virtually indistinguishable, which indicates that the new method of constructing initial data does not introduce additional constraint-violating modes. The middle and bottom panels of Fig. 9 show the evolution of the Christodoulou mass and the spin magnitude of the large black hole. The differences between the evolutions of the old and new initial data sets are consistent with truncation error. Thus we conclude that the quasi-local quantities are the same in both data sets.
Turning attention to the trajectories of the black holes, we find a stark difference in the evolutions. Figure 10 shows the motion of the large black hole in inertial coordinates for both runs. The uncontrolled residual linear momentum P z ADM in the old initial data causes the centre of mass of the binary to drift linearly during the evolution, as shown in the right panel of Figure 10. Such a drift may have multiple undesirable consequences. Most immediately, it causes the gravitational wave extraction spheres to be off-center from the center-of-mass of the binary, which will cause mixing of the spherical harmonic modes of the gravitational radiation, an effect discussed in more detail below. Moreover, SpEC's constraint preserving outer boundary conditions [56,57,58] are designed to work best for low-order spherical harmonic modes of the outgoing radiation. If the binary is offset relative to the outer boundary (for instance due to a drift of the center of mass), higher order spherical harmonic components will become more important, possibly leading to an additional runaway acceleration of the center of mass [59].
To examine the dynamics of the binary, we study the orbital frequency vector Ω ≡ r ×ṙ/|r| 2 . The left panel of Figure 11 shows the projection of Ω onto the unit sphere, making it apparent that the precession and nutation dynamics are very similar until very close to merger. The right panel shows a plot of Ω ≡ |Ω| from which several features are apparent. The evolution of Ω is qualitatively the same in both cases, consistent with expectation that removing a coordinate motion of the centre of mass does not change the binary dynamics. Likewise the initial pulse of junk radiation (inset A) appears quite similar. However, small oscillations in Ω are more pronounced in the new code (inset C). This is reflected in the measured values of the eccentricities: e = 10 −4 for the old, e = 2.5×10 −4 for the new code. The difference in eccentricity arises because the new term v i o in the outer boundary condition Eq. (28) does slightly modify the content of the initial data. In this particular case, |v 0 | ∼ 10 −3 , so that it is not unreasonable to expect the orbital eccentricity to change by a comparable magnitude. The initial orbital frequency Ω 0 , initial radial velocityȧ 0 , and initial separation D 0 listed in Table 1 were tuned to result in essentially vanishing eccentricity in the old initial data [53]. The new initial data constructed from the identical initial data parameters must therefore have a slightly larger eccentricity. If we had tuned to vanishing eccentricity with the new initial data, then the old initial data would exhibit the larger eccentricity. The evolutions of the old and new initial data also result in a different time to merger, cf. panel B of Fig. 11. This difference could be caused either by the slightly different inspiral dynamics like eccentricity, or could simply be due to truncation error of our low resolution evolution.
Finally, we examine the waveforms for the two runs. Most strikingly, the movement of the coordinate centre of mass shown in Figure 10 is also reflected in the sphericalharmonic decomposition of the waveform. This is most easily seen in the sub-dominant modes. Figure 12 shows the ( , m) = (3, 1) modes of the spin-weighted sphericalharmonic (SWSH) decompositions of the waveforms h old measured from the old initial data and h new measured from the new initial data. Since gravitational waves in SpEC are extracted on a coordinate sphere centered on the origin, a drifting source mixes the modes of the SWSH decomposition. As seen in the lower panel of the figure, this mixing introduces very large effects. To verify that these effects are primarily due to the motion of the center of mass, we have also transformed h new to a frame in which the center of mass is moving as in the original initial data. The initial position of h new is transformed to agree with the center of mass of the old initial data as measured by Eq. (25), and its velocity is transformed to agree with P ADM /M ADM of the old initial data as measured by Eq. (19). This transformation is applied entirely at future null infinity by the method The h 3,1 waveform modes, as measured in evolutions of the original and the new initial data, and extrapolated to future null infinity, I + . The upper panel shows the waveform h new from the new initial data, measured in a frame that is centered at the origin of the simulation coordinates. The lower panel shows the same data with a transformation applied on I + as described in the text, as well as the waveform h old from the original initial data measured in its simulation coordinatesin which the black holes are moving as shown in Fig. 10. Essentially, the center of mass is stationary at the origin in the upper panel, and is moving in the lower panel.
described in [60], and is a special case of a BMS transformation [61,62]. It will thus be seen in any waveforms, whether extrapolated [63] (as seen here) or extracted by Cauchy characteristic methods [64,65,66,67]. As shown in the lower panel of Fig. 12, the transformation reproduces the features seen in h old very well.
Mode decompositions like this one are used very frequently for analyzing numerical models, and for constructing analytical models. If they are unmodeled and uncontrolled, effects like those seen in the lower panel will simply appear to be errors in the waveform. This could negatively impact uncertainty estimates of numerical simulations, error estimates for analytical waveforms, or calibration of waveform models to numerical results. These effects will also be present in any calculation that uses the waveforms to compute physical quantities such as the flux of linear and angular momentum. By removing extraneous displacements and boosts, this new initial data code simplifies such analyses. + + The drift described here is a linear motion due to residual linear momentum in initial data. Controlling this drift will not help for other types of motion present in very long simulations; see [59].

Discussion
Numerical evolution of binary black hole spacetimes requires accurate initial data. In this work we have improved the initial data techniques in SpEC to allow access to a much wider parameter space of generically precessing high mass ratio, high-spin binaries. A more flexible domain decomposition allows for stable solution for high-mass ratio and high spin binaries. An enhanced root-finding algorithm is used to achieve desired physical parameters for the binary. This becomes important when a naive analytic Jacobian is not appropriate, which is precisely the case for high mass ratios and spins, see Figure 8. Adaptive mesh refinement drastically improves efficiency and robustness of the code, displaying exponential convergence of constraints, c.f. Figure 7. Finally, a new method to control the linear momentum is used to eliminate a linear drift of the centre of mass during evolution. This in turn nullifies spurious gravitational mode mixing, which is of paramount importance for construction of hybrid waveforms or calibration of phenomenological models as demonstrated by Figure 12.
An interesting application of the improved initial data code is the construction of initial data for hyperbolic encounters. Such systems have been studied in the past (e.g. [68,69,70]) and provide a laboratory for exploring strong field physics in a different regime than the binary inspiral. Using the new code, we have successfully constructed initial data for hyperbolic encounters for a selection of mass ratios and spins, which was not possible before in SpEC. As a simple example, we evolve two systems of two equal mass black holes that are initially separated by 60M and have a velocity of ≈ 0.14c. Both systems have the same impact parameter b N R = 15M , and differ only in the black hole spins: In one case the black holes are non-spinning, in the other both holes have dimensionless spins χ = 0.5 initially in the x direction. Figure 13 shows the trajectories of the two black holes. In the presence of spin, the spin-orbit interactions cause the plane of scattering to change and also change the deflection angle of the hyperbolic encounter. Exploration of other parameters is left to future investigations.
Another application is the construction of initial data for binaries with very small initial separation, corresponding to only a few orbits before merger. This is useful if one is interested in the properties of the merger remnant, e.g. for calibrating analytical waveform models but evolving a long inspiral is too computationally expensive. As an example, we construct initial data for a system with q = 21, χ 1 = 0.66, χ 2 = 0.41 (oriented in random directions) and initial orbital frequency of M Ω = 0.032, and initial coordinate separation D 0 = 8.82 M. We note that initial data for binaries near ISCO at high mass-ratio is challenging and further work remains to be done to make it robust for q > 10 regime.