Improvements to the construction of binary black hole initial data
Abstract
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 massratio above 10:1, and for dimensionless black hole spins above 0.9, while improving efficiency for lower massratios 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.
pacs:
04.25.Dm, 04.30.Db, 04.70.Bw1 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 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 7dimensional parameter space consisting of the mass ratio and the dimensionless spins remains unexplored. Most of the attention has been focused on binaries close to equal mass () and modest spin () (although see [19, 20, 21, 22, 23, 24]) For stellar mass black hole binaries, one can expect mass ratios and arbitrary spin magnitudes and orientations, which leads to precession of the spins and the orbital plane. Precessing, high massratio binaries have interesting dynamics, causing large modulations of the gravitational waveform. One can expect even higher mass ratios () for neutron starblack hole (NSBH) binaries (see [25] a BHWolfRayet system with BH mass ). At high mass ratios, BBH systems can be used as proxies for NSBH systems(e.g. [26]). One would thus like to simulate highmass ratio BBH systems.
Intermediate mass black holes (IMBH) with masses 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 () 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 stellarmass  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 quasiequilibrium, 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 quasiequilibrium 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 massratio, 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 domaindecomposition to allow a wider range of
massratios, spins and separations.
(ii) Adaptive meshrefinement to enhance computational efficiency and
to ensure robust numerical convergence
for massratios and dimensionless spins .
(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 equalmass binary with aligned spins
of , and a singlespinning binary with
spin of 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.
2 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 blackhole binary in quasiequilibrium). In SpEC, the former is achieved by using a pseudospectral 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 quasiequilibrium conditions, SpEC employs the extended conformal thin sandwich (XCTS) formalism [34]. Before solving the conformal thinsandwich 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 of the initial data hypersurface. Therefore, iterative rootfinding 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.
2.1 Domain decomposition
Figure 1 indicates the geometry of the domaindecomposition 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.
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, , 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 domaindecomposition in order for the solver to converge.
The new domaindecomposition 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, and , and of the outer spherical shell, . Note that is not required to lie on the line connecting and .

The inner and outer radii of the inner spherical shells and the outer spherical shell, , , , and , , .
The remaining parameters and 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 . 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 via
(1) Note that 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 ,
(2) The multiplier must satisfy to ensure that the blocks cover the entire open region within cylinders 0, 2, and 4.

The multipler , which measures how much larger the outer size of the cylinders is compared to the inner edge of the outer spherical shell:
(3) To ensure complete overlap between the cylinders and the sphere C, , with being the distance from point to the axis of the cylinders.
The value of will determine the relative size of the face of the blocks to the inner spheres: If , then the edge of the block will be entirely outside the inner spherical shell. Conversely, if , 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 gridinternal geometry coefficients are , , , and . We have found these choices to be robust for a wide variety of component masses, spins and separations.
2.2 Adaptive mesh refinement
An important factor in efficiently generating highaccuracy 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 hardcoded resolutions, tuned to equalmass 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 multidomain spectral solver represents the solution in each subdomain as a tensorproduct of basisfunctions. 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 in
the th mode by the rootmeansquare value of all the coefficients
of the th mode across all spectral coefficients of the other
basisfunctions. For instance, in a spherical shell with spectral
expansion of the form
,
the radial power would be
,
where represents the number of angular
coefficients
For the expected spectral convergence, should decay exponentially as a function of [45, 32], i.e. when plotted vs. should be a straight line. The slope of this line represents the decrease in the magnitude of the spectral coefficients when going from mode to mode . We estimate 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 , we can reach a target accuracy by adding
(4) 
modes to the spectral expansion (recall and a higher accuracy means a lower ).The answer is rounded up so that if the current accuracy is worse than the target accuracy, and we set 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 , and sometimes AMR yields the same resolution at neighboring . Because the ADMquantities are exclusively evaluated in the outer spherical shell (cf. Sec. 2.4 below), this would result in apparent nonconvergence 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 gridpoints in the direction, whenever AMR triggers an adjustment to the domain decomposition.
2.3 Iterative determination of free parameters
When constructing initial data, we wish to achive desired masses , and desired black hole spin vectors and . The free data, however, is instead given by the radii and angular frequencies of the apparent horizons and , which we write as
(5) 
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
(6) 
The standard approach to the problem would be to use Newton’s method; however, evaluating the Jacobian is too expensive numerically as every evaluation of the function requires an elliptic solve. We instead use the following approach: make an initial guess based on the Kerr expressions for both black holes,
(7)  
(8) 
and perform an elliptic solve for . We then construct an analytic Jacobian 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]:
(9) 
where . This corresponds to the “secant” approximation for a function of one variable. Finally we set
(10) 
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 gridsizes to a desired truncation error; and rootfinding 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 rootfinding errors. We do so by using an overall truncation error target for AMR. We start with a large value for , corresponding to a small gridsize. As rootfinding residuals decrease, we will decrease . We proceed as follows: At iteration , we compute two measures of progress in root finding: First, the residual which quantifies how close the physical parameters are to their desired values. is simply the rms error in the physical parameters:
(11) 
Second, the improvement that indicates how quickly rootfinding converges, defined as
(12) 
where .
We monitor 2 conditions:
where is the desired truncation error, and and are tunable parameters. The first condition assures that the resolution is increased if the rootfinding 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 , with 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 by a factor of 10 and continue with the next itertion. For all cases we have encountered, the choices and have proven to be robust.
2.4 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 using the ArnowittDeserMisner (ADM) prescription. Our initial data satify the asymptotic gauge conditions [47]
(13)  
(14) 
needed to remove ambiguities in the definition of the ADM angular momentum, as well as the boundary condition on . The old code [32, 48] directly evaluated the resulting surface integrals at infinity [49, 47],
(15)  
(16) 
using extrapolation in powers of to infinite radius [48]. is then found to be a combination of terms of , and a combination of terms. The old technique, therefore, is very sensitive to small errors in in the outermost sphere of our computational domain (the outer boundary is typically located at ) and particularly to the presence of constraint violating modes in that sphere. Typically, this leads to large errors in at low resolution, and large errors in 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 as the sum of a surface integral on a sphere located at a smaller radius and a volume integral. Utilizing , we write
(17)  
Here the normal to points into the interior of (e.g. along if it is a coordinate sphere) and the factor was inserted to eliminate terms with spatial derivatives of from Eq. (21). Using the momentum constraint in the absence of sources,
(18) 
the volume integral can be simplified to
(19) 
Here,
(20)  
(21) 
where are the connections derived from the conformal metric . Note that for conformal flatness and maximal slicing, and the volume integral disappears (see [50]).
In practice, for conformally curved initial data, The outer spherical shell extends to outer radius . Therefore, in the numerical evaluation of the volume integral in Eq. (19), the volume element associated with the outermost gridpoint becomes very large and introduces numerical noise. To avoid this, we roll off the integrand beyond a certain radius , i.e. we replace by given by
(22) 
We choose , where 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
(23) 
with cyclical permutations of yielding the other components. For maximal slicing and conformal flatness in , Eq. (23) simplifies to
(24) 
Because Eq. (23) relies on the cancellation of large volume terms, it can be sensitive to errors in . Accordingly, we use Eq. (24) using a surface at sufficiently large radius such that in the metric is conformally flat and .
To illustrate the importance of the transformations applied to the ADM integrals, we consider the convergence test for configuration q50. We evaluate using Eq. (15) and Eq. (19), and we evaluate using Eq. (16) and Eq. (23). Figure 2 shows the results.
The calculation of is improved by about one order of magnitude when utilizing Gauss’ law, whereas improves by several orders of magnitude. We point out that, in order to achieve any convergence for the old calculation, we had to manually increase the radial resolution in the outer sphere by whenever the domain decomposition is adjusted.
We also compute a new diagnostic, the centreofmass of the initial data sets using the formalism developed in Ref. [51]. In conformal flatness, the expressions from [51] reduce to
(25) 
where is the outwardpointing unit normal, . Equation (25) is numerically evaluated by expanding the conformal factor in a powerseries in . We read off the (angledependent) coefficient of the term, and expand this coefficient in spherical harmonics. Each individual spherical harmonic term can be integrated against analytically, so that the integral (25) collapses to a linear combination of sphericalharmonic coefficients.
2.5 Control of linear momentum and centre of mass
The quasiequilibrium conformal thinsandwich 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 , 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 timederivatives vanish in a frame rotating with orbital velocity about the orign. One finally incorporates a radial expansion factor , which allows fine control of the orbital eccentricity [39, 52, 53, 41]. By a suitable choice of the conformal quantities, the quasiequilibrium approach can generate initial data with black hole spins of order [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 . The past implementations use the location of the black holes to partially control . Consider a small displacement applied to the centres of both excision regions. Through the orbital rotation about the origin, the displacement induces a change in velocity of the black holes of , with a corresponding change in . Therefore, could be used to cancel the components of orthogonal to ; however, the crossproduct in prevented any correction parallel to . For headon collisions with , no control of is possible at all. For the nonprecessing simulations presented in [41], the component of parallel to vanishes by symmetry, and no problems arose. However, for generic precessing binaries, there will be a nonzero 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 , while simultaneously allowing us to control the center of mass as well. We fix the relative separation of the centres of the excision spheres,
(26) 
where the separation vector is userspecified. We use the choice of to control the centerofmass of the binary. Once a first initial data set is computed (with, in general, ), we can update
(27) 
With the blackhole centres now used to control the centre of mass, we need a different means to control . We add in the outer boundary condition on the shift (Eq. (38c) of [40]) a constant velocity :
(28) 
Here represents the outer boundary, a sphere with radius
. The velocity will effect the overall
motion of the binary, and will be reflected in a corresponding change
in by , where
is the ADMenergy of the binary. During iterative
rootfinding of the free parameters, we adjust to
achieve
To motivate the updating formula for , consider a perturbation of by , and a perturbation of by . If we allow the masses to vary, then a Newtonianinspired formula is
(29) 
To summarize, relative to earlier initialdata sets, we modify the outer boundary condition for the shift by the term , cf. Eq. (28), and use updating formulae (27) and (29) to adjust and . Section 2.4 describes how we compute and .
We demonstrate the efficiency of the updating formulas Eqs. (27,29) in Fig. 3 that shows the magnitude of and as a function of rootfinding iteration for a precessing binary (case q10 in Table 1). The convergence is evidently very fast, with the final values of and respectively. This means that even for an inspiral lasting M, the drift of the centre of mass due to residual linear momentum in initial data will be only .
3 Numerical results
3.1 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 . 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 by 10, until a truncation error of 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, , closely tracks the truncation error target .
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 gridpoints . The top panel demonstrates the exponential decrease in the norms of the Hamiltonian and momentum constraints, which implies that this data set is constraintsatisfying. The constraints are given explicitly by:
(30)  
(31) 
where is the covariant derivative associated with the spatial metric. The norm is simply the normalized pointwise norm over all collocation points:
(32) 
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:
(33) 
Once again, the convergence is essentially exponential. The bottom panel of Fig. 5, finally examines the convergence of the ADM quantities and the centerofmass 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 .
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 . Figure
6 shows the components of
as a function of rootfinding iteration for both
the new and the old code.
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, 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 points, an improvement of a factor of
over the old code.
3.2 Rootfinding
It is also important to examine the performance of the updated rootfinding procedure based on Broyden’s method. Figure 8 shows the rootfinding 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 rootfinding procedure described in Sec. 2.3 results in errors of order . Finally, for case Spin0.9999, the results are drastically different: the old code has errors in the masses of order several and spins of order . Since we are attempting to construct a binary with dimensionless spins of it becomes clear that the old code is inadequate for this purpose. On the other hand, the new rootfinding procedure successfully reduces the errors in physical quantities to the level of . 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 2550% 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.
3.3 Exploratory evolution
We have emphasized above the importance of controlling . We now evolve initial data for case q3 constructed with the old and the new initialdata code, and compare the two evolutions in detail.
We being by considering the convergence of constraints and quasilocal quantities during the evolution. The top panel of Figure 9 shows the 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 constraintviolating 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 quasilocal 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 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 offcenter from the centerofmass 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 loworder 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
. 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: for the old,
for the new code.
The difference in eccentricity arises because the new term
in the outer boundary condition Eq. (28) does
slightly modify the content of the initial data. In this particular
case, , so that it is not unreasonable to expect the orbital
eccentricity to change by a comparable magnitude. The initial orbital
frequency ,
initial radial velocity , and initial separation
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 subdominant modes. Figure 12 shows the modes of the spinweighted sphericalharmonic (SWSH) decompositions of the waveforms measured from the old initial data and 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 to a frame in which the center of mass is moving as in the original initial data. The initial position of 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 of the old initial data as measured by Eq. (19). This transformation is applied entirely at future null infinity by the method 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 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.
4 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, highspin binaries. A more flexible domain decomposition allows for stable solution for highmass ratio and high spin binaries. An enhanced rootfinding 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 . Both systems have the same impact parameter , and differ only in the black hole spins: In one case the black holes are nonspinning, in the other both holes have dimensionless spins initially in the direction. Figure 13 shows the trajectories of the two black holes. In the presence of spin, the spinorbit 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 , , (oriented in random directions) and initial orbital frequency of , and initial coordinate separation M. We note that initial data for binaries near ISCO at high massratio is challenging and further work remains to be done to make it robust for regime.
Acknowledgements
We thank Geoffrey Lovelace, Larry Kidder and Mark Scheel for helpful
discussions. Calculations were performed with the SpECcode [31]. We gratefully acknowledge support
from NSERC of Canada, from the Canada Research Chairs Program, and
from the Canadian Institute for Advanced Research. FF gratefully
acknowledges support from the Vincent and Beatrice Tremaine
Postdoctoral fellowship at CITA. Support for this work was provided by
NASA through Einstein Postdoctoral Fellowship grant numbered
PF4150122. We further gratefully acknowledge support from the
Sherman Fairchild Foundation; from NSF Grants PHY1306125 and
AST1333129 at Cornell; and from NSF Grants No. PHY1440083 and
AST1333520 at Caltech. Calculations were performed at the Gravity
cluster and the GPC supercomputer at the SciNet HPC
Consortium [71]; SciNet is funded by: the Canada Foundation
for Innovation (CFI) under the auspices of Compute Canada; the
Government of Ontario; Ontario Research Fund (ORF) – Research
Excellence; and the University of Toronto. Further calculations were
performed on the Briarée cluster from Sherbrooke University,
managed by Calcul Québec and Compute Canada. The operation of this
supercomputer is funded by the Canada Foundation for Innovation (CFI),
Ministère de l’Économie, de l’Innovation et des Exportations
du Québec (MEIE), RMGA and the Fonds de recherche du Québec 
Nature et Technologies (FRQNT).
References
Footnotes
 : Class. Quantum Grav.
 For spherical harmonic basisfunctions, the top two modes are filtered [32] and are therefore not included in the data .
 Using the obtained vector as the shiftvector in an evolution results in a translating outer boundary; this effect is eliminated by evolving with a shift vector of .
 Both codes compute in the same way (via \crefeq:PADM), but differ in the way it is controlled.
 We note that the case q50 could not be constructed with the old code.
 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].
References
 A. Einstein. Über gravitationswellen. In Sitzber. Preuss. Akad. Wiss., page 154, Jan 1918.
 B.S. Sathyaprakash and B.F. Schutz. Physics, Astrophysics and Cosmology with Gravitational Waves. Living Rev.Rel., 12:2, 2009.
 G. Hobbs, A. Archibald, Z. Arzoumanian, D. Backer, M. Bailes, et al. The international pulsar timing array project: using pulsars as a gravitational wave detector. Class.Quant.Grav., 27:084013, 2010.
 Gregory M. Harry. Advanced LIGO: The next generation of gravitational wave detectors. Class.Quant.Grav., 27:084006, 2010.
 J. Aasi et al. Advanced LIGO. Class.Quant.Grav., 32:074001, 2015.
 The Virgo Collaboration. Advanced Virgo Baseline Design, 2010. VIRâ027Aâ09.
 J. Abadie et al. Predictions for the Rates of Compact Binary Coalescences Observable by Groundbased Gravitationalwave Detectors. Class. Quant. Grav., 27:173001, 2010.
 Lee S. Finn. Detection, measurement, and gravitational radiation. Phys. Rev. D, 46(12):5236, December 1992.
 L. Lehner and F. Pretorius. Numerical Relativity and Astrophysics. Ann.Rev.of Astron. & Astroph., 52:661–694, August 2014.
 Harald P. Pfeiffer. Numerical simulations of compact object binaries. Class. Quant. Grav., 29:124004, 2012.
 Ian Hinder. The Current Status of Binary Black Hole Simulations in Numerical Relativity. Class. Quant. Grav., 27:114004, 2010.
 Sean T. McWilliams. The status of blackhole binary merger simulations with numerical relativity. Class. Quantum Grav., 28:134001, 2011.
 Mark Hannam. Status of blackholebinary simulations for gravitational wave detection. Class. Quant. Grav., 26:114001, 2009.
 Larne Pekowsky, Richard O’Shaughnessy, Jim Healy, and Deirdre Shoemaker. Comparing gravitational waves from nonprecessing and precessing black hole binaries in the corotating frame. Phys. Rev. D, 88:024040, 2013.
 B. Aylott, J. G. Baker, W. D. Boggs, M. Boyle, P. R. Brady, D. A. Brown, B. Brügmann, L. T. Buchman, A. Buonanno, L. Cadonati, J. Camp, M. Campanelli, J. Centrella, S. Chatterji, N. Christensen, T. Chu, P. Diener, N. Dorband, Z. B. Etienne, J. Faber, S. Fairhurst, B. Farr, S. Fischetti, G. Guidi, L. M. Goggin, M. Hannam, F. Herrmann, I. Hinder, S. Husa, V. Kalogera, D. Keppel, L. E. Kidder, B. J. Kelly, B. Krishnan, P. Laguna, C. O. Lousto, I. Mandel, P. Marronetti, R. Matzner, S. T. McWilliams, K. D. Matthews, R. A. Mercer, S. R. P. Mohapatra, A. H. Mroué, H. Nakano, E. Ochsner, Y. Pan, L. Pekowsky, H. P. Pfeiffer, D. Pollney, F. Pretorius, V. Raymond, C. Reisswig, L. Rezzolla, O. Rinne, C. Robinson, C. Röver, L. Santamaría, B. Sathyaprakash, M. A. Scheel, E. Schnetter, J. Seiler, S. L. Shapiro, D. Shoemaker, U. Sperhake, A. Stroeer, R. Sturani, W. Tichy, Y. T. Liu, M. van der Sluys, J. R. van Meter, R. Vaulin, A. Vecchio, J. Veitch, A. Viceré, J. T. Whelan, and Y. Zlochower. Testing gravitationalwave searches with numerical relativity waveforms: Results from the first Numerical INJection Analysis (NINJA) project. Class. Quantum Grav., 26(16):165008, 2009.
 P. Ajith, Michael Boyle, Duncan A. Brown, Bernd Brugmann, Luisa T. Buchman, et al. The NINJA2 catalog of hybrid postNewtonian/numericalrelativity waveforms for nonprecessing blackhole binaries. Class. Quantum Grav., 29(12):124001, 2012.
 P. Ajith, Michael Boyle, Duncan A. Brown, Bernd Brugmann, Luisa T. Buchman, et al. Addendum to ’The NINJA2 catalog of hybrid postNewtonian/numericalrelativity waveforms for nonprecessing blackhole binaries’. Class. Quantum Grav., 30(19):199401, 2013.
 Ian Hinder et al. Erroranalysis and comparison to analytical models of numerical waveforms produced by the NRAR Collaboration. Classical and Quantum Gravity, 31(2):025012, 2014.
 Carlos O. Lousto and Yosef Zlochower. Orbital evolution of extrememassratio blackhole binaries with numerical relativity. Phys. Rev. Lett., 106:041101, January 2011.
 Ulrich Sperhake, Vitor Cardoso, Christian D. Ott, Erik Schnetter, and Helvi Witek. Extreme black hole simulations: collisions of unequal mass black holes and the point particle limit. 2011.
 Mark A. Scheel, Matthew Giesler, Daniel A. Hemberger, Geoffrey Lovelace, Kevin Kuper, Michael Boyle, and Béla Szilágyi. 2014. submitted to Class. Quantum Grav.
 Vijay Varma, Parameswaran Ajith, Sascha Husa, Juan Calderon Bustillo, Mark Hannam, and Michael Puerrer. Gravitationalwave observations of binary black holes: Effect of nonquadrupole modes. 2014. arXiv:1409.2349.
 Juan CalderÃ³n Bustillo, Alejandro BohÃ©, Sascha Husa, Alicia M. Sintes, Mark Hannam, et al. Comparison of subdominant gravitational wave harmonics between postNewtonian and numerical relativity calculations and construction of multimode hybrids. 2015.
 Deirdre Shoemaker, Karan Jani, Lionel London, and Larne Pekowsky. Connecting Numerical Relativity and Data Analysis of Gravitational Wave Detectors. Astrophys.Space Sci.Proc., 40:245–258, 2015.
 J. M. Silverman and A. V. Filippenko. On IC 10 X1, the Most Massive Known StellarMass Black Hole. Astrophys. J. Lett., 678:L17–L20, 2008.
 F. Foucart, L. Buchman, M. D. Duez, M. Grudich, L. E. Kidder, I. MacDonald, A. Mroue, H. P. Pfeiffer, M. A. Scheel, and B. Szilagyi. First direct comparison of nondisrupting neutron starblack hole and binary black hole merger simulations. Phys. Rev. D, 88(6):064017, September 2013.
 M. Coleman Miller and E. J. M. Colbert. IntermediateMass Black Holes. Int. J. Mod. Phys. D, 13:1–64, 2004.
 M. Mezcua, T. P. Roberts, A. P. Lobanov, and A. D. Sutton. The powerful jet of an offnuclear intermediatemass black hole in the spiral galaxy ngc 2276. Monthly Notices of the Royal Astronomical Society, 448(2):1893–1899, 2015.
 Dheeraj R. Pasham, Tod E. Strohmayer, and Richard F. Mushotzky. A 400solarmass black hole in the galaxy m82. Nature, 513(7516):74–76, Sep 2014. Letter.
 Greg Cook. Initial data for numerical relativity. Living Rev. Rel., 3, November 2000. 5.
 http://www.blackholes.org/SpEC.html.
 H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A. Teukolsky. A multidomain spectral method for solving elliptic equations. Comput. Phys. Commun., 152:253–273, 2003.
 James W. York. Conformal “thinsandwich” data for the initialvalue problem of general relativity. Phys. Rev. Lett., 82(7):1350–1353, Feb 1999.
 Harald P. Pfeiffer and James W. York. Extrinsic curvature and the Einstein constraints. Phys. Rev. D, 67(4):044022, Feb 2003.
 Gregory B. Cook. Corotating and irrotational binary black holes in quasicircular orbits. Phys. Rev. D, 65(8):084003, Mar 2002.
 Gregory B. Cook and Harald P. Pfeiffer. Excision boundary conditions for blackhole initial data. Phys. Rev. D, 70(10):104016, Nov 2004.
 Matthew Caudill, Greg B. Cook, Jason D. Grigsby, and Harald P. Pfeiffer. Circular orbits and spin in blackhole initial data. Phys. Rev. D, 74(6):064011, 2006.
 Abdul H. Mroue, Mark A. Scheel, Bela Szilagyi, Harald P. Pfeiffer, Michael Boyle, Daniel A. Hemberger, Lawrence E. Kidder, Geoffrey Lovelace, Sergei Ossokine, Nicholas W. Taylor, Anil Zenginoglu, Luisa T. Buchman, Tony Chu, Evan Foley, Matthew Giesler, Robert Owen, and Saul A. Teukolsky. A catalog of 174 binary black hole simulations for gravitational wave astronomy. Phys. Rev. Lett., 111:241104, 2013.
 Harald P. Pfeiffer, Duncan A. Brown, Lawrence E. Kidder, Lee Lindblom, Geoffrey Lovelace, and Mark A. Scheel. Reducing orbital eccentricity in binary black hole simulations. Class. Quantum Grav., 24(12):S59–S81, 2007.
 Geoffrey Lovelace, Robert Owen, Harald P. Pfeiffer, and Tony Chu. Binaryblackhole initial data with nearlyextremal spins. Phys. Rev. D, 78:084017, 2008.
 Luisa T. Buchman, Harald P. Pfeiffer, Mark A. Scheel, and Béla Szilágyi. Simulations of unequal mass binary black holes with spectral methods. Phys. Rev. D, 86:084033, 2012.
 Francois Foucart, Matthew D. Duez, Lawrence E. Kidder, and Saul A. Teukolsky. Black holeneutron star mergers: effects of the orientation of the black hole spin. Phys. Rev. D, 83:024005, 2011.
 Katherine Henriksson, FranÃÂ§ois Foucart, Lawrence E. Kidder, and Saul A. Teukolsky. Initial data for highcompactness black holeneutron star binaries. 2014.
 Béla Szilágyi. Key Elements of Robustness in Binary Black Hole Evolutions using Spectral Methods. Int.J.Mod.Phys., D23(7):1430014, 2014.
 John P. Boyd. Chebyshev and Fourier Spectral Methods. Dover Publications, second edition, 1999.
 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes: The Art of Scientific Computing (3rd Ed.). Cambridge University Press, New York, 2007.
 J. W. York, Jr. Kinematics and dynamics of general relativity. In L. L. Smarr, editor, Sources of Gravitational Radiation, pages 83–126, 1979.
 Harald P. Pfeiffer, Gregory B. Cook, and Saul A. Teukolsky. Comparing initialdata sets for binary black holes. Phys. Rev. D, 66:024047, 2002.
 R. Arnowitt, S. Deser, and Charles W. Misner. The dynamics of general relativity. In L. Witten, editor, Gravitation: An Introduction to Current Research, pages 227–265. Wiley, New York, 1962.
 Francois Foucart, Lawrence E. Kidder, Harald P. Pfeiffer, and Saul A. Teukolsky. Initial data for black holeneutron star binaries: a flexible, highaccuracy spectral method. Phys. Rev. D, 77:124051, 2008.
 D. Baskaran, S.R. Lau, and A.N. Petrov. Centerofmass integral in canonical general relativity. Annals Phys., 307:90–131, 2003.
 Michael Boyle, Duncan A. Brown, Lawrence E. Kidder, Abdul H. Mroué, Harald P. Pfeiffer, Mark A. Scheel, Gregory B. Cook, and Saul A. Teukolsky. Highaccuracy comparison of numerical relativity simulations with postNewtonian expansions. Phys. Rev. D, 76(12):124038, 2007.
 Alessandra Buonanno, Lawrence E. Kidder, Abdul H. Mroué, Harald P. Pfeiffer, and Andrea Taracchini. Reducing orbital eccentricity of precessing blackhole binaries. Phys. Rev. D, 83:104034, 2011.
 Serguei Ossokine, Lawrence E. Kidder, and Harald P. Pfeiffer. Precessiontracking coordinates for simulations of compactobjectbinaries. Phys. Rev. D, 88:084031, 2013.
 Lee Lindblom, Mark A. Scheel, Lawrence E. Kidder, Robert Owen, and Oliver Rinne. A new generalized harmonic evolution system. Class. Quantum Grav., 23:S447–S462, 2006.
 Oliver Rinne, Luisa T. Buchman, Mark A. Scheel, and Harald P. Pfeiffer. Implementation of higherorder absorbing boundary conditions for the Einstein equations. Class. Quantum Grav., 26:075009, 2009.
 Oliver Rinne, Luisa T. Buchman, Mark A. Scheel, and Harald P. Pfeiffer. Implementation of absorbing boundary conditions for the Einstein equations. In Proceedings of Spanish Relativity Meeting, volume 1122, pages 384–387. AIP Conference Proceedings, 2008.
 L. T. Buchman and O. C. A. Sarbach. Improved outer boundary conditions for Einstein’s field equations. Class. Quantum Grav., 24:S307–S326, 2007.
 Bela Szilagyi, Jonathan Blackman, Alessandra Buonanno, Andrea Taracchini, Harald P. Pfeiffer, et al. Numerical relativity reaching into postNewtonian territory: a compactobject binary simulation spanning 350 gravitationalwave cycles. 2015.
 Michael Boyle. Transformations of asymptotic gravitationalwave data. Not yet published, 2015.
 R. K. Sachs. Gravitational waves in general relativity. VIII. waves in asymptotically flat spacetime. Proc. R. Soc. Lond. A, 270(1340):103–126, October 1962.
 Ezra Newman and Roger Penrose. An approach to gravitational radiation by a method of spin coefficients. J. Math. Phys., 3(3):566–578, 1962.
 Michael Boyle and Abdul H. Mroué. Extrapolating gravitationalwave data from numerical simulations. Phys. Rev. D, 80(12):124045–14, December 2009.
 Nigel T. Bishop, Roberto Gómez, Luis Lehner, and Jeffrey Winicour. Cauchycharacteristic extraction in numerical relativity. Phys. Rev. D, 54(10):6153–6165, Nov 1996.
 Maria Babiuc, Béla Szilágyi, Ian Hawke, and Yosef Zlochower. Gravitational wave extraction based on Cauchycharacteristic extraction and characteristic evolution. Class. Quantum Grav., 22(23):5089–5107, 2005.
 Nicholas W. Taylor, Michael Boyle, Christian Reisswig, Mark A. Scheel, Tony Chu, Lawrence E. Kidder, and Béla Szilágyi. Comparing gravitational waveform extrapolation to Cauchycharacteristic extraction in binary black hole simulations. Phys. Rev. D, 88:124010, Dec 2013.
 C. J. Handmer, B. Szilágyi, and J. Winicour. Gauge invariant spectral characteristic extraction.
 Thibault Damour, Federico Guercilena, Ian Hinder, Seth Hopper, Alessandro Nagar, et al. StrongField Scattering of Two Black Holes: Numerics Versus Analytics. Phys.Rev., D89(8):081503, 2014.
 U. Sperhake, V. Cardoso, F. Pretorius, E. Berti, T. Hinderer, and N. Yunes. Cross section, final spin and zoomwhirl behavior in highenergy black hole collisions. Phys. Rev. Lett., 2009. in press.
 James Healy, Frank Herrmann, Ian Hinder, Deirdre M. Shoemaker, Pablo Laguna, and Richard A. Matzner. Superkicks in hyperbolic encounters of binary black holes. Phys. Rev. Lett., 102:041101, 2009.
 Chris Loken, Daniel Gruner, Leslie Groer, Richard Peltier, Neil Bunn, Michael Craig, Teresa Henriques, Jillian Dempsey, ChingHsing Yu, Joseph Chen, L Jonathan Dursi, Jason Chong, Scott Northrup, Jaime Pinto, Neil Knecht, and Ramses Van Zon. SciNet: Lessons Learned from Building a Powerefficient Top20 System and Data Centre. J. Phys.: Conf. Ser., 256:012026, 2010.