Improved methods for simulating nearly extremal binary black holes
Mark A. Scheel, Matthew Giesler, Daniel A. Hemberger, Geoffrey Lovelace, Kevin Kuper, Michael Boyle, Bela Szilagyi, Lawrence E. Kidder
IImproved methods for simulating nearly extremal binary black holes
Mark A. Scheel, Matthew Giesler,
1, 2
Daniel A. Hemberger, Geoffrey Lovelace,
2, 1
Kevin Kuper, Michael Boyle, B´ela Szil´agyi, and Lawrence E. Kidder Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125, USA Gravitational Wave Physics and Astronomy Center, California State University Fullerton, Fullerton, California 92834, USA Center for Radiophysics and Space Research, Cornell University, Ithaca, New York 14853, USA (Dated: September 3, 2018)Astrophysical black holes could be nearly extremal (that is, rotating nearly as fast as possi-ble); therefore, nearly extremal black holes could be among the binaries that current and futuregravitational-wave observatories will detect. Predicting the gravitational waves emitted by mergingblack holes requires numerical-relativity simulations, but these simulations are especially challengingwhen one or both holes have mass m and spin S exceeding the Bowen-York limit of S/m = 0 . S/m = 0 .
99. We also use these methods tosimulate a non-precessing binary black hole coalescence, where both black holes have
S/m = 0 . I. INTRODUCTION
Second-generation interferometers such as AdvancedLIGO, Virgo, and KAGRA [1–4] will soon begin search-ing for gravitational waves. To increase the number ofgravitational-wave detections and to maximize what wecan learn about the waves’ sources, we require accuratetheoretical models of the sources and the emitted gravi-tational radiation.The inspiral, merger, and ringdown of binary blackholes (BBHs) are among the most promising astrophys-ical sources of gravitational waves. As the black holesorbit, they lose energy to gravitational radiation, inspi-raling until they collide and merge to form a final blackhole (the “remnant”) that eventually settles to a station-ary Kerr state.A BBH is characterized by 7 intrinsic parameters: thespin angular momenta (cid:126)S of each hole and the mass ratio q . The spin magnitude of a black hole is often charac-terized by the dimensionless quantity χ ≡ S/m , where S = | (cid:126)S | , m is the black hole mass, and we use geometrizedunits where c = G = 1. A black hole with the maxi-mum possible dimensionless spin χ = 1 is called extremal.There is considerable uncertainty in the expected massratios and spins of astrophysical BBHs that are likely tobe detected by gravitational-wave interferometers; how-ever, there is evidence that nearly extremal black holesexist in nature. For instance, recent measurements ofstellar-mass black holes (such as Cygnus X-1 [5–7], GRS1915+105 [8], and GX 339-4 [9]) and supermassive blackholes (such as Swift J0501.9-3239 [10]) suggest that therecould be a population of black holes with spins of χ ∼ χ = 1 remainsalmost completely unexplored. Numerical simulationsof nearly extremal, merging black holes are especiallychallenging. One reason for this is that initial data fora BBH must satisfy the Einstein constraint equations,but the simplest method for constructing constraint-satisfying initial data, the Bowen-York method [27–29],cannot yield initial data with nearly extremal black holes.This is because the Bowen-York construction assumesthat the initial spatial geometry is conformally flat (i.e.,that the initial spatial metric is proportional to the met-ric of flat space). Conformally flat spacetimes cannotrepresent black holes that i) are in equilibrium, and ii)possess linear [30] or angular [31, 32] momentum; there-fore, conformally-flat spinning black holes are out of equi-librium and will quickly relax to an equilibrium state.Specifically, Bowen-York puncture initial data can pro-duce BBHs with initial spins as large as χ = 0 . χ = 0 . a r X i v : . [ g r- q c ] D ec data through inspiral, merger, and ringdown is especiallychallenging for two reasons, as discussed in Refs. [36–38]. First, the portion of the spacetime near the horizonsrequires very high resolution (and thus high computa-tional cost), since metric gradients are much larger thanfor lower spins. Second, for evolution methods that ex-cise the singularities inside each black hole and evolveonly the exterior region, constructing and maintaining asuitable computational domain that keeps each excisionboundary just inside the corresponding apparent hori-zon becomes more and more challenging as the spin ap-proaches extremality.In this paper, we use the phrase “nearly extremal”to refer to χ > .
93, i.e., to a black hole with a spinabove the Bowen-York limit. Note that a black hole with χ = 0 .
93 is significantly less extremal than a black holewith χ = 0 . χ . For instance, if the rotational energy of a Kerrblack hole with a fixed mass is denoted E rot ( χ ), then E rot (0 . /E rot (1) is only 59%, while E rot (0 . /E rot (1)is 92 .
5% (c.f., Fig. 1 of Ref. [36]). Furthermore, the to-tal energy that a BBH emits in gravitational waves alsoscales nonlinearly with χ . For example, for equal massesand equal spins aligned with the orbital angular momen-tum, a BBH with χ = 1 radiates 10% more energy thana BBH with χ = 0 .
93, whereas a BBH with χ = 0 . χ = 0(Eq. (9) of Ref. [41]). Nonlinear scaling with χ is alsoseen for binaries consisting of a black hole and a neutronstar: the amount of neutron-star matter remaining out-side the black hole just after tidal disruption increasesvery rapidly with black-hole spin (Fig. 10 of Ref. [38]).Several groups have constructed and evolved Bowen-York puncture initial data with spins near (but below)the Bowen-York limit [42–45]. Recently, Ruchlin etal. [46] constructed and evolved puncture initial data fora head-on collision of two black holes with equal massand spins of magnitude χ = 0 .
99. Only four previ-ously published simulations [24, 36, 37] out of hundredspublished to date contain the quasi-circular coalescenceof BBHs with χ > .
93. These four simulations wereevolved using the Spectral Einstein Code (SpEC) [47]from “superposed-Kerr-Schild” excision initial data [48]and have equal masses and equal spins aligned or anti-aligned with the orbital angular momentum.In this paper, we present technical improvements thathave enabled us to simulate BBHs with black-hole spinsup to χ = 0 .
994 ( i.e., E rot (0 . /E rot (1) = 87 . II. TECHNIQUES
We carry out numerical simulations with the SpectralEinstein Code (SpEC) [47]. We construct [51] quasi-equilibrium [48, 52] constraint-satisfying [53] initial databased on a weighted superposition of two boosted, spin-ning Kerr-Schild black holes [48]. We use an iterativemethod to produce initial data with low eccentricity [54–56].We use a generalized harmonic formulation [57–60] ofEinstein’s equations and damped harmonic gauge [61–63]to evolve the initial data. The adaptively-refined [36, 64]grid extends from pure-outflow excision boundaries justinside the apparent horizons [63, 65–67] to an artificialouter boundary, where we enforce constraint-preservingboundary conditions [60, 68, 69]. The grid has only oneexcision boundary after the holes merge [65, 66]. We usea pseudospectral fast-flow algorithm [70] to find appar-ent horizons, and we compute spins on these apparenthorizons using the approximate Killing vector formalismof Cook, Whiting, and Owen [71, 72].In the remainder of this section, we describe new tech-niques that allow simulations of binaries with large black-hole spins. Large spins are difficult for two reasons. First,the metric gradients near the black-hole horizons becomelarger with larger spin, making increased numerical res-olution necessary in this region. Second, black-hole exci-sion is more difficult: in SpEC, we remove the physicalsingularity inside each black hole by placing an artifi-cial excision boundary just inside each apparent horizonand evolve only the region exterior to all excision bound-aries. We find that the maximum required coordinatedistance between an excision boundary and the corre-sponding horizon becomes smaller with larger spin, sothat our algorithm for dynamically adjusting the exci-sion boundaries to track the size, shape, and motion ofthe horizons must be more accurate. We consider bothof these difficulties below.Not all of the improvements discussed here were nec-essary for all of the simulations described in Sec. III. Forexample, the simulation discussed in Sec. III A succeededwithout some of the grid and control system improve-ments; however, these improvements became necessarywhen simulating even larger black-hole spins (Sec. III B)or allowing generic spin directions and unequal masses(Sec. III C).
A. Grid improvements
Meeting the need for high resolution near the hori-zons is accomplished via spectral adaptive mesh refine-ment [64]. This includes both p -type refinement (chang-ing the number of collocation points in a given spectralsubdomain) and h -type refinement (adding, removing, orchanging the distribution of subdomains). The simula-tions described here used the algorithm detailed in [64],with adjustments to default parameters so as to allow forhigher resolution. In particular: We increased the num-ber of radial collocation points in a spherical subdomainthat forces h -refinement from 20 points to 40, we dis-abled angular h -refinement in the spherical subdomainsthat touch the excision boundary so as to retain a singlespherical subdomain at this boundary, and we increasedthe allowed number of spherical-harmonic coefficients inspherical shells from L = 40 to L = 80. Note that most ofthese changes (such as allowing up to L = 80) were nec-essary only for a small portion of the simulation when thehorizons are highly distorted, such as during the initial“junk radiation” transients (when spurious gravitationalradiation is emitted as the BBH relaxes to equilibrium)and near the moment of merger.We also reduced the initial distance between the appar-ent horizons and the excision boundaries. To understandthis change, note that when solving elliptic equations forinitial data, the excision boundaries are made to coin-cide exactly with the apparent horizons via boundaryconditions imposed on those surfaces. But for the evolu-tion, the excision boundaries must be slightly inside thehorizons, so that the horizons are fully contained in thecomputational domain and therefore can be determinedby the apparent horizon finder. To accomplish this, afterthe initial data have been determined, these data are ex-trapolated slightly inside the horizons to a new excisionboundary, before the evolution begins. For large spins,this extrapolation occurs in the region where metric gra-dients are growing rapidly as r decreases, so placing theexcision boundary at a larger r reduces those gradients. To carry out some of the simulations shown here, wemoved the initial excision boundary radius from 94% to98% of the initial horizon radius. B. Control system improvements
Several of the improvements necessary for handlinghigh spins involve control systems used to adjust map-pings between coordinate systems. These control systemsand the mappings are described in detail in [66]. Here webriefly summarize important points, and we discuss keydifferences from [66].
1. Summary of size and shape control systems
In SpEC, we remove the physical singularity insideeach black hole by placing an artificial excision boundaryjust inside each apparent horizon, evolving only the re-gion exterior to all excision boundaries. We use multiplecoordinate systems to handle excision of black holes thatare moving and changing shape [63, 65, 66, 73–76]. Wecall “inertial coordinates” ¯ x i those asymptotically iner-tial coordinates in which the black holes orbit each other,have a distorted and dynamical shape, and approach eachother as energy is lost to gravitational radiation. We ap-ply spectral methods in a different coordinate system,“grid coordinates” x i , in which the excision boundariesare spherical and time-independent. We connect thesetwo coordinate systems with an analytic mapping func-tion M : x i → ¯ x i that depends on a set of time-dependentparameters λ ( t ). These parameters λ ( t ) are adjusted au-tomatically by feedback control systems so that, as theapparent horizons of the black holes move and changeshape (in the inertial frame), the excision boundaries aremapped to inertial-coordinate surfaces that follow thismotion and remain just inside the apparent horizons.The control of all parameters λ ( t ) is accomplished inthe same way, using a general control system we have de-veloped, as described in [66]. The part of the algorithmthat distinguishes one λ ( t ) from another is the specifica-tion of the control error Q ( t ), which is different for eachcontrol parameter. For example, the λ Scaling ( t ) that rep-resents the distance between the excision boundaries hasa different Q ( t ) than the matrix λ Rotation ( t ) that repre-sents the rotation of the inertial coordinates with respectto the grid coordinates. If there exists a desired value of λ ( t ), call it λ target , which depends on observables A (suchas the positions or shapes of the apparent horizons) butdoes not depend on λ itself, then we define Q ( t ) = λ target ( A ) − λ ( t ) . (1)For more general situations in which λ target depends on λ itself, we generalize the above definition: we require that λ takes on its desired value when Q = 0, and we requirethat ∂Q∂λ = − O ( Q ) . (2)Given Q ( t ), our algorithm adjusts the corresponding λ ( t )so that Q ( t ) is driven towards zero; this driving occurs ona timescale τ d that is determined dynamically and thatis different for each control system.The full map from grid to inertial coordinates is ¯ x i = M x i , where M = M Translation ◦ M
Rotation ◦ M
Scaling ◦M Skew ◦ M
CutX ◦ M
Shape . (3)Each of these maps is described in detail in Sec. 4 of [66]. Shape control.
Here we are concerned only with thelast map, M Shape , which is defined as: x i (cid:55)→ x i (cid:32) − (cid:88) H f H ( r H , θ H , φ H ) r H (cid:88) (cid:96)m Y (cid:96)m ( θ H , φ H ) λ H(cid:96)m ( t ) (cid:33) . (4)The index H goes over each of the two excised regions A and B , and the map is applied to the grid-frame co-ordinates. The polar coordinates ( r H , θ H , φ H ) centeredabout excised region H are defined in the usual way,the quantities Y (cid:96)m ( θ H , φ H ) are spherical harmonics, and λ H(cid:96)m ( t ) are expansion coefficients that parameterize themap near excision region H ; these λ H(cid:96)m ( t ) are the co-efficients that we adjust using a control system. Thefunction f H ( r H , θ H , φ H ) is chosen to be unity near ex-cision region H and zero near the other excision region,so that the distortion maps for the two black holes aredecoupled; see Eq. 72 and Fig. 4 of [66] for a precise def-inition of f H ( r H , θ H , φ H ). In the following, the controlsystems for each excised region H , while independent,are identical in operation, so we will omit the H labelsfor simplicity.We control λ (cid:96)m ( t ) so that each excision boundary isdriven to the same shape as the corresponding apparenthorizon; this results in conditions on λ (cid:96)m ( t ) for (cid:96) > λ ( t ) unconstrained [66]. Size control.
The size of the excision boundary, asencoded in the remaining coefficient λ ( t ), must satisfytwo conditions. Horizon tracking.
The first is that the excision bound-ary remains inside the apparent horizon. To satisfy thiscondition, we first write the shape of each apparent hori-zon as an expansion in spherical harmonics, parameter-ized in terms of polar coordinates about the center of thecorresponding excision boundary,ˆ r AH (ˆ θ, ˆ φ ) = (cid:88) (cid:96)m ˆ S (cid:96)m Y (cid:96)m (ˆ θ, ˆ φ ) , (5)where the intermediate frame ˆ x i is connected to the gridframe by the map M Distortion = M CutX ◦ M
Shape . (6) By construction, M Distortion leaves invariant the centersof the excision boundaries, and the angles with respectto these centers. Then we choose Q = ˙ˆ S (∆ r − − ˙ λ (7)where ∆ r = 1 − (cid:104) ˆ r EB (cid:105)(cid:104) ˆ r AH (cid:105) (8)is the relative difference between the average radius of theapparent horizon (in the intermediate frame) and the av-erage radius of the excision boundary. The angle brack-ets in Eq. (8) represent averaging over angles. Choosing Q ( t ) according to Eq. (7) drives d/dt (∆ r ) towards zero,so that the excision boundary remains a fixed (relative)distance inside the apparent horizon. Characteristic speed tracking.
The second conditionthat must be satisfied by λ ( t ) involves characteristicspeeds of the evolved Einstein equations: well-posednessof our system of equations requires that all of the char-acteristic speeds be non-negative, i.e. characteristic fieldsmust flow into the black hole. The minimum character-istic speed at each excision boundary is given by v = − α − ¯ β i ¯ n i − ¯ n i ∂ ¯ x i ∂t , (9)where α is the lapse, ¯ β i is the shift, and ¯ n i is the normalto the excision boundary pointing out of the computa-tional domain , i.e., toward the center of the hole. It ispossible to write [66] v = v + ˆ n i x i r Y ˙ λ , (10)where v collects all terms that are independent of ˙ λ .Therefore, a control system that controls ˙ λ and drives v to some target value v T can be constructed by definingthe control error Q = (min( v ) − v T ) / (cid:104)− Ξ (cid:105) , (11)where Ξ = ˆ n i x i r Y , (12)and the minimum is over the excision boundary. Notethat Ξ < n i and x i /r point in opposite direc-tions. Switching between horizon and characteristic speedtracking.
Note that Eqs. (7) and (11) specify two differ-ent control systems that control the same degree of free-dom, λ : the first control system, which we call “hori-zon tracking”, adjusts ˙ λ to control ∆ r , and the other,which we call “characteristic speed tracking”, adjusts ˙ λ to control v . Both ∆ r and v must remain nonnegative fora successful evolution, but we cannot use both Eqs. (7)and (11) simultaneously. Furthermore, changes in ˙ λ af-fect ∆ r and v in the opposite direction: if ˙ λ increases,∆ r increases, but the characteristic speed v decreases.In practice (Sec. II B 3), we now alternate between thetwo control systems, Eqs. (7) and (11). We monitor both v and ∆ r as functions of time and predict whether eitherof these quantities is likely to become negative in theimmediate future; if so, we estimate the timescale τ v or τ ∆ r on which this will occur. If τ ∆ r is small enough, i.e.∆ r is in imminent danger of becoming negative, we useEq. (7) to control ∆ r . If τ v is small enough that v is indanger of becoming negative, we use Eq. (11) to control v . The details of how we make these decisions have beenimproved since the description in Sec. 5.3 of [66], so wedescribe the improved algorithm below.
2. Improvements in gain scheduling
We now describe improvements in the control systemsthat were necessary for our new high-spin simulations tosucceed.
Comoving characteristic speed as a control sys-tem diagnostic.
We define a new quantity v c which wecall the comoving characteristic speed : v c = − α − ¯ β i ¯ n i − ¯ n i ∂ ¯ x i ∂ ˆ t + ˆ n i x i r (cid:34) Y ˙ˆ S (∆ r −
1) + (cid:88) (cid:96)> Y (cid:96)m ( θ, φ ) ˙ λ (cid:96)m ( t ) (cid:35) . (13)The comoving characteristic speed v c is what the charac-teristic speed v would be if Q ( t ) in Eq. (7) were exactlyzero, i.e. if ∆ r were constant in time. In other words,if we turn on horizon tracking, the control system drives v toward v c . This tells us (for instance) that if we find v c <
0, we should not use horizon tracking, since horizontracking would drive v to a negative value. The instan-taneous value of v c is independent of ˙ λ and roughlyindependent of λ ; the only dependence on λ comesfrom the smooth spatial variation of the metric functions.Hence, v c is a useful quantity for separating the effectsof the control system for ˙ λ from the effects of othercontrol systems.One way we use v c is in determining whether our con-trol system for ˙ λ will fail. During a simulation, v c isusually positive, but it routinely becomes negative forshort periods of time, particularly when the shapes ofthe horizons are changing rapidly, for example near t = 0when the black holes are ringing down from initial “junkradiation” transients. However, if v c becomes negativeand remains so indefinitely, our control system for ˙ λ must eventually fail. This is because for v > v c < r must be decreasing, so if we keep v > Control error damping timescale improve-ments.
For many of the high-spin SpEC simulations that failed before we made the improvements describedin this paper, we observed that v c < λ ( t ) parameters other than λ ; in par-ticular, the shape parameters λ (cid:96)m for (cid:96) >
0. In otherwords, the shape and position of the excision boundarydiffered from the shape and position of the horizon by asufficient amount that it was not possible to make both v > r > λ .This particular problem was fixed by changing the al-gorithm for setting the tolerance on the control error Q ( t ), for all Q ( t ) except Q . Associated with each ofour control systems is a timescale parameter τ d which isadjusted dynamically. The control error Q ( t ) is dampedlike e − t/τ d , under the assumption that τ d is smaller thanall other timescales in the problem. Therefore decreasing τ d results in smaller values of Q ( t ). The previous methodof adjusting τ d is described in Sec. 3.3 of [66]: at regulartime intervals t i , the timescale is changed according to τ i +1 d = βτ id , (14)where β = . , if ˙ Q/Q > − / τ d and | Q | or | ˙ Qτ d | > Q Max t . , if | Q | < Q Min t and | ˙ Qτ d | < Q Min t , otherwise . (15)Here Q Min t and Q Max t are thresholds for the control error Q , set to constant values Q Max t = 2 × − m A /m B + m B /m A (16) Q Min t = 14 Q Max t . (17)In the new algorithm, we make three changes. The firstis that Q Max t is no longer a constant: instead, it is chosento be Q Max t = a ∆ r min , where a is a constant (typicallychosen to be 0.05( m A + m B ) for those Q values withdimensions of length, and 0.05 for those Q values thatare dimensionless) and ∆ r min is the minimum relativedistance between the excision boundary and the apparenthorizon: ∆ r min = min ˆ θ, ˆ φ (cid:32) − ˆ r EB (ˆ θ, ˆ φ )ˆ r AH (ˆ θ, ˆ φ ) (cid:33) . (18)The second change is that we define an estimate of thetime that the horizon will cross the excision surface τ ∆ r cross = − ∆ r min (cid:18) ddt ∆ r min (cid:19) − , (19)and if τ ∆ r cross > τ id > τ ∆ r cross , then we set τ i +1 d = τ ∆ r cross instead of using Eq. (14).Both of the above changes force each Q ( t ) to be closerto zero when the excision boundary approaches the hori-zon. A third, minor, change we make in the algorithm tMax = constQ tMax = const ∆ r min FIG. 1. The control error Q ( t ) associated with one particularcontrol system for two different S ++0 . simulations which differonly in the treatment of Q Max t for that control system. Theblack dashed curve shows the case in which Q Max t is a con-stant, given by Eq. (16), and the red solid curve shows the casein which Q Max t is chosen to be 0 . r min , with ∆ r min givenby Eq. (18). The former simulation crashes at t ∼ M . affects only the behavior of Q ( t ) at early times: the ini-tial values of each τ d were decreased so that each Q ( t )is smaller at earlier times; these initial values are speci-fied separately from the tolerances Q Max t that determinewhen τ d is modified. The effect of the first change, set-ting Q Max t proportional to ∆ r min , is illustrated in Fig. 1for one particular control system. In Fig. 1 and the re-mainder of the paper, M ≡ m A + m B is the sum of theChristodoulou masses of the two black holes at the time t relax when the initial “junk radiation” transients havedecayed away.
3. Size control: switching between Eqs. (7) and (11).
At every time step, the control system for λ is gov-erned by a Q given by either Eq. (7) or Eq. (11), with anassociated damping timescale τ d and (if using characteris-tic speed control) a target speed v T . At regular intervals(typically every time step), the algorithm has an oppor-tunity to change from using Eq. (7) to using Eq. (11) The Q ( t ) illustrated here is the one for the control system thatcomputes a smooth approximation ˆ r appx AH ( t ) to the average hori-zon radius; this approximate value is used to compute ˙ˆ S and∆ r in Eq. (7), in order to reduce the number of calls to the com-putationally expensive horizon finder (see section 7 and Eq. (108)of [66] for details). or vice versa, and to choose a new value of τ d and (ifusing characteristic speed control) v T . Here we describehow we make these choices. A previous version of thisalgorithm was described in [66], but many improvementshave been made since then.Because the goal of the λ control system is to keepboth v and ∆ r min positive, we regularly monitor v and∆ r min as functions of time. We predict whether eitherof these quantities is likely to become negative in theimmediate future, and if so, we estimate the timescale τ v or τ ∆ r min on which this will occur, using the methoddescribed in Appendix C of [66]. Because the sign of v c is important to the success of horizon tracking, we alsomonitor v c as a function of time, and if it is positive anddecreasing, we predict the timescale τ v c on which it willbecome negative. If v , v c , or ∆ r min are increasing insteadof decreasing, we define the corresponding timescale τ v , τ v c , or τ ∆ r min to be infinite.We begin by determining whether v is in imminentdanger of becoming negative, so that some immediateaction must be taken to prevent this from occurring. Weregard v to be in danger if τ v < τ d and τ v < τ ∆ r min .Furthermore, if characteristic speed control is in effect,we additionally require τ v < σ τ d and v < σ v T to deem v in danger; here σ (cid:46) σ ∼ . Thefirst requirement, τ v < σ τ d , prevents the algorithm fromswitching back and forth between characteristic speedcontrol and horizon tracking on each time step. The sec-ond requirement, v < σ v T , prevents the control sys-tem from rapidly decreasing the characteristic speed toachieve a target v T that is less than v .If v is deemed to be in danger, the action taken toprevent v from becoming negative depends on the currentstate of the control system. If characteristic speed controlis in effect, then it remains in effect, and τ d is set equal to τ v in order to drive v towards v T more quickly. If horizontracking is in effect, and if v c < v c is decreasing, thencharacteristic speed control goes into effect, with v T = σ v , and τ d is left unchanged. The constant σ , typically1.01, prevents the control system from switching backand forth on each timestep. Finally, if horizon trackingis in effect, and if v c > v c is nondecreasing, thenhorizon tracking remains in effect and we reduce τ d by afactor of σ < v toward v c , which is in no danger of becoming negative.If v is deemed not to be in danger, then we checkwhether ∆ r min is in danger of soon becoming negative.We regard ∆ r min to be in danger if τ ∆ r min < τ v and if τ ∆ r min < σ τ d , where σ is a constant typically chosento be 20. Furthermore, if horizon tracking is in effect,we additionally require τ ∆ r min < σ τ d to deem τ ∆ r min indanger, where σ < Labels for control system constants like σ i and η are consistentwith the notation in Ref. [66]. C h a r ac t e r i s ti c s p ee d v (allow horizon tracking when v c <0)v (prevent horizon tracking when v c <0)v c FIG. 2. Characteristic speed v and comoving characteristicspeed v c for two different S ++0 . simulations that differ onlyin the algorithm for treating the situation in which ∆ r min is deemed in danger while characteristic speed control is ineffect and v c <
0. The dashed red curve shows v for a sim-ulation in which horizon tracking becomes active in this sit-uation; the code crashes early, at only t ∼ . M . The solidred curve shows v for a simulation in which for this situationcharacteristic speed control remains in effect, but the targetcharacteristic speed is reduced as described in the text. Thequantity v c is the same for both simulations. condition prevents the control system from switching onevery time step.If ∆ r min is in danger, the action again depends on thestate of the control system and other variables. If horizontracking is in effect, then it remains in effect, and τ d is setequal to τ ∆ r min in order to drive ∆ r min to a constant morequickly. If characteristic speed control is in effect and if v c >
0, then horizon tracking goes into effect, and τ d isset equal to τ ∆ r min . We require v c > v towards v c , and we wish to maintain v >
0; if horizon trackingbecomes active even if v c <
0, the simulation often fails,as shown in Fig. 2. To solve the problem illustrated byFig. 2, when the code finds that ∆ r min is in danger whilecharacteristic speed control is in effect and if v c <
0, thenthe code allows characteristic speed control to remain ineffect, but it sets the new τ d to min( τ d , τ ∆ r min ), and itreduces v T to ηv , where η < v T will reduce v but will increase ∆ r min . If v c < v T will occuras needed. As mentioned above, if v c < otherthan the one for λ to attempt to make v c positive, asdiscussed in Sec. II B 2.If neither v nor ∆ r min are in imminent danger ofbecoming negative, then the system attempts to findan equilibrium using horizon tracking. If characteristicspeed control is in effect, and if v c >
0, ˙ v c ≥
0, and either v > v T or v c > v , then horizon tracking goes into effect,using the current τ d . However, if both v and v c are de- creasing, horizon tracking does not go into effect unless v is decreasing faster than v c and τ v c > σ τ d , where σ is a constant we usually set to 5. The purpose of thesevarious conditions on v , v c , and their derivatives and pre-dicted zero-crossing times is to prevent horizon trackingfrom going into effect when it is likely that a switch backto characteristic speed control will soon be necessary. Forexample, if ˙ v c < v c is decreasing faster than v , thenwe anticipate that v c will soon become negative, in whichcase horizon tracking is inappropriate because it woulddrive v towards zero.The behavior of the control system depends on var-ious constants σ i (1 < i <
7) and η described above;these constants govern decisions made by the algorithm.These constants have restricted values (e.g. η should notbe greater than unity), but they were chosen without anyfine tuning. Changing their values slightly will change de-tails such as the exact value of τ d at a particular timestep,but we expect that small changes in parameters will notchange whether a simulation succeeds or fails, and willchange physical results only at the level of truncationerror (because the control system changes the grid coor-dinates).Occasionally when horizon tracking is in effect, we findthat the value of ∆ r min is excessively large or small. If itis excessively small, then τ d becomes small, and we areforced to reduce the timestep in the evolution equationsto keep the control system stable, resulting in a largecomputational expense. If it is too large, then the exci-sion boundary lies deep inside the horizon, and excessivecomputational resources are needed to resolve the largegradients. Therefore, we allow a drift term to sometimesbe added to Eq. (7), as discussed in [66]. III. SIMULATIONS
We present three new simulations, summarized in Ta-ble I. We will refer to quantities defined in Table Ithroughout the remainder of this paper. The techniquesdescribed in Sec. II were essential to the successful com-pletion of these simulations.
A. Equal-mass, aligned spins χ = 0 . The first simulation we present, and refer to as S ++0 . , isan equal-mass case in which each black hole has a spin of χ = 0 .
99 aligned with the orbital angular momentum. At t = t relax the simulation has M ω orb = 0 . M isthe sum of the relaxed Christodoulou masses. The binarythen evolves through 25 orbits, merger and ringdown.This simulation took 83 days on 48 cores for the highestresolution.To assess numerical convergence, we perform severalsimulations that are identical except for the numericalresolution, which we label by an integer N . Larger N corresponds to finer resolution, but the absolute scale Name Catalog ID t relax q r m rA m rB Mω r orb χ rA θ rA /π φ rA /π χ rB θ rB /π φ rB /π e N M f χ f S ++0 . SXS:BBH:0177 320.0 1.0 0.5 0.5 0.0154 0.989 0.00 – 0.989 0.00 – 12.6 25.4 0.888 0.949 S ++0 . SXS:BBH:0178 640.0 1.0 0.5 0.5 0.0157 0.994 0.00 – 0.994 0.00 – 8.6 25.4 0.887 0.950 S . . SXS:BBH:0179 380.0 1.5 0.6 0.4 0.0148 0.991 0.00 0.73 0.200 0.24 0.23 322.4 23.8 0.922 0.897TABLE I. Summary of physical simulation parameters. Data are publicly available online [77] indexed by their Catalog ID.Quantities with an r superscript are reported at time t = t relax , the time after the initial “junk radiation” transients havesettled down: q is the mass ratio, m H is the Christodoulou mass of an individual black hole (where H represents black hole A or B ), Mω orb is the orbital frequency, χ H is the dimensionless spin, θ H is the angle between (cid:126)ω orb and (cid:126)χ H , and φ H is theangle between the separation vector and the component of (cid:126)χ H in the orbital plane. The remaining quantities are eccentricity e , number of orbits N from t = 0 to merger, final Christodoulou mass M f , and final spin magnitude χ f . -5 0 5x/M-505 y / M FIG. 3. The trajectories of the centers of the individual ap-parent horizons for the highest resolution of S ++0 . . of N is different for different physically distinct simula-tions. The value of N enters the simulation through thetolerance in adaptive mesh refinement (AMR): the AMRtruncation error tolerance is chosen to be proportionalto e − N . For each value of N , we compute the complexphase φ of the (cid:96) = 2 , m = 2 component of Ψ . Wethen take the difference ∆ φ between φ computed usingotherwise-identical simulations using different values of N .Figure 4 shows these differences for S ++0 . . No align-ment of the waveforms in time or phase has been per-formed. Note the rapid convergence: ∆ φ between N = 3and N = 4 (labeled “4-3”) is significantly smaller than∆ φ between the two lower resolutions. Also note thatthe difference “3-2” is nearly the same as “4-2”, indicat-ing that this difference effectively measures the numericaltruncation error in the N = 2 simulation. Similarly, thedifference “4-3” represents the numerical truncation er-ror in the N = 3 simulation. Furthermore, one wouldexpect that the truncation error in the N = 4 simulation t/M − − ∆ φ FIG. 4. Convergence test for S ++0 . . Shown are gravitational-wave phase differences between Ψ computed using differentvalues of the numerical resolution parameter N . Several dif-ferences are shown, and labeled by the values of N that arecompared, e.g. “3-2” means N = 3 versus N = 2. Waveformsare extracted at a finite radius r = 465 M , and no alignmentof waveforms was performed. is smaller than the “4-3” curve by another order of mag-nitude (although it would be necessary to run an N = 5simulation to actually measure this).In the S ++0 . initial data, the spin of each black holeis 0 .
99. When the system is evolved, the spins decreasevery slightly for the first ∼ M as initial transients prop-agate away from the horizons, as shown in the upper insetof Fig. 5. Then the spins level off and become roughlyconstant, but with a small negative slope. All values ofresolution N agree quite well, and the higher two resolu-tions are indistinguishable in Fig. 5. The spins decreasemore rapidly just before merger ( t ∼ M ). The com-mon horizon first appears with a spin greater than thefinal value, and then relaxes as the remnant black holesettles down, as shown in the lower inset of Fig. 5. Thefinal spin is χ f = 0 . χ χ indiv χ merg FIG. 5. Spin magnitude as a function of time for S ++0 . . Atearly times, the spin of one of the apparent horizons is shownat resolutions N = 2 (black solid), N = 3 (red dotted) and N = 4 (blue dashed). A closeup of early times is shown in theupper inset. At late times, the spin of the merged apparenthorizon is shown as a function of time for the same resolutions,and a closeup of late times is shown in the lower inset. the difference between the two highest resolution simula-tions.The radiated energy fraction E rad is the relative changein energy of the binary from t = −∞ to t = ∞ and canbe computed from E rad ≡ − E ∞ E −∞ = 1 − M f M . (20)The final Christodoulou mass M f is the energy of thesystem at t = ∞ , because the remnant is in equilibriumat the end of the simulation; the total Christodouloumass M at t = t relax is the energy of the system at t = −∞ , because the individual black-hole masses changeby less than one part in 10 between t = −∞ and t = t relax (see, e.g. Eq. 14 in Ref. [78]). We find that E rad = 11 . χ f = 0 . E rad = 11 . χ = 0 .
97 is expected to overestimatethe final spin (see Fig. 6 in Ref. [41]) and underestimate t/M − − ∆ φ FIG. 6. Convergence test for S ++0 . . Labels are the sameas for Fig. 4. For N (cid:54) = 5, the simulations were started at t branch = 1414 M , using the N = 5 solution as initial data. the final radiated energy (see Fig. 8 in Ref. [41]), and thisis what we find with S ++0 . . B. Equal-mass, aligned spins χ = 0 . We repeated the equal-mass aligned-spin simulationabove, but with a larger spin. We refer to this case as S ++0 . . The initial data were chosen with χ = 0 .
995 foreach black hole, but the spins drop to χ = 0 . t = 10 M of evolution time, a much smallertimescale than the relaxation time t relax (this rapid ini-tial decrease in spin can also be seen for S ++0 . in theupper inset of Fig. 5). The simulation S ++0 . representsthe largest spin ever simulated for a black-hole binary. Ithas M ω orb = 0 . t = t relax , and then proceedsthrough 25 orbits, merger, and ringdown. The high-est resolution completed in approximately 71 days on 48cores. Note that this simulation, S ++0 . , was computa-tionally cheaper than the lower-spin simulation, S ++0 . ,and achieved a smaller overall phase error (see Figs. 4and 6). This is due to code optimization that was donebetween the time that the S ++0 . and S ++0 . simulationswere carried out; for the same version of SpEC, thereis actually a steep increase in computational cost as afunction of spin.Obtaining convergence was more difficult for this sim-ulation than for S ++0 . . The reason is that it is difficultto fully resolve the initial transients, sometimes called“junk radiation”, that result from imperfect initial data.If these transients are unresolved, then the small changesin masses, spins, and trajectories caused by these tran-0sients are effectively random, and therefore otherwise-identical simulations run with different values of resolu-tion N will differ by random small amounts that will notconverge with increasing N . So to investigate conver-gence, we remove the initial transients in the followingway. We first carry out a simulation with one value of N , call it N base . In the case of S ++0 . , N base representsthe highest resolution. Then we choose some fiducial time t = t branch > t relax at which we decide that the transientshave decayed away. We then carry out simulations with N (cid:54) = N base starting at t = t branch , using the N = N base solution as initial data. This procedure removes the ef-fects of the transients from our convergence tests.However, this procedure alone was insufficient toachieve convergence. When convergence is rapid enoughin a particular subdomain so that adding a single gridpoint results in a large decrease in truncation error, itis possible for two different AMR truncation error toler-ances, e.g. e N and e N − , to result in the same numberof grid points for that subdomain. This makes the trun-cation error in that subdomain identical for two differentvalues of N , which spoils convergence tests for simula-tions with those values of N . To remedy this problemin such cases, we increase the spacing in truncation errortolerance as a function of level N : the truncation er-ror tolerance is proportional to 10 N instead of e N . This,combined with the procedure to remove the effect of tran-sients, results in good convergence, as shown in Fig. 6.The spin of the remnant black hole is χ f = 0 . E rad = 11 . χ f = 0 . E rad = 11 . C. Unequal-mass, precessing
The final simulation we present is an unequal-mass casewith q = 1 .
5, in which the larger black hole has a spinof χ = 0 .
99 aligned with the orbital angular momentum,while the smaller black hole has a spin magnitude of χ =0 . S . . ,using a notation similar to that introduced earlier. Thesimulation has M ω orb = 0 . t = t relax , and thenproceeds through 23 orbits, merger, and ringdown. Thissimulation took approximately 26 days on 48 cores forthe highest resolution using the same optimized versionof SpEC as the S ++0 . case described in Sec. III B.We found that for S . . we needed to remove the effectof unresolved initial transients and increase the spacingin AMR truncation error tolerance to obtain acceptableconvergence results. To do this we followed the sameprocedure as for S ++0 . , described in Sec. III B. Figure 7shows good convergence of the gravitational-wave phasedifference when using this procedure. t/M − − ∆ φ FIG. 7. Convergence test for S . . . Labels are the same asfor Figure 4. For N (cid:54) = 4, the simulations were started at t branch = 1362 M , using the N = 4 solution as initial data.FIG. 8. Coordinate trajectories (green and purple lines) of theblack holes and coordinate shapes of the individual and com-mon apparent horizons (surfaces) at the moment of merger,for S . . . The horizons are colored according to their vortic-ity [79]. Figure 8 shows the trajectories of the centers of the ap-parent horizons for this simulation, as well as the individ-ual apparent horizons and the common apparent horizonat the moment when the common horizon first appears.Trajectories and horizon shapes are shown in the asymp-totically inertial coordinate system used in the simula-tion. Because the spin of the smaller hole (cid:126)χ B is notaligned with the orbital angular momentum, the systemprecesses, so the trajectories do not lie in a plane.Figure 9 shows the precession of the spin and orbitalfrequency vectors in S . . . The spin (cid:126)χ A and orbital fre-quency (cid:126)ω orb initially point along the z -axis. Because themisaligned spin (cid:126)χ B is on the smaller black hole and ismuch smaller in magnitude than (cid:126)χ A , it has a minimaleffect on the orbital dynamics, so (cid:126)χ A and (cid:126)ω orb remainnear the z -axis throughout the simulation. Therefore,we consider the precession to be mild. As angular mo-mentum is carried away by gravitational radiation, theopening angles of the precession cones change. The an-1 ˆ ! orb ˆ A ˆ B FIG. 9. Precession of the spins and orbital frequency forthe highest-resolution simulation N = 4 of S . . . The unitvector spins, ˆ χ A and ˆ χ B , and orbital frequency ˆ ω orb trace theprecession on the unit sphere. The precession curves for N =2 and N = 3 converge to the N = 4 curves shown here, andthe curves for N = 3 and N = 4 are nearly indistinguishable. gles of (cid:126)ω orb and (cid:126)χ A with respect to the z -axis increasefrom 0 ◦ at t = 0 to 6 ◦ and 12 ◦ , respectively, at the timeof merger. In contrast, the angle of (cid:126)χ B with respect tothe z -axis decreases from 45 ◦ to 12 ◦ . The spins (cid:126)χ A and (cid:126)χ B complete 2.1 and 2.5 precession cycles, respectively,and (cid:126)ω orb completes 2.4 precession cycles.The spin of the remnant black hole is χ f = 0 . E rad = 7 . χ f = 0 . E rad = 7 . IV. RESULTSA. Spin evolution during inspiral
During the inspiral, the tidal field of each black hole af-fects its companion, and this interaction slowly changesthe black-hole masses and spins as a function of time.For aligned spins, Alvi [78] has derived perturbative ex-pressions for the time rate of change of the mass andspin of a black hole in a binary. Chatziioannou, Poisson, To evaluate the quantities S || and ∆ || in Ref. [25], we used the z -component of (cid:126)S and (cid:126) ∆ at t relax , which should be strictly validonly for non-precessing binaries. Also, the formula for χ f inRef. [25] requires evaluating certain quantities at the innermoststable circular orbit (ISCO) of a Kerr black hole with a spinof χ f , so that χ f is not given in closed form; for simplicity weevaluate the ISCO quantities using the measured χ f from thesimulation. t/M -8 -6 -4 | d S / d t | / M Low resMedium resHigh resAlvi (0PN)AlviCPY (1PN) × -7 × -7 FIG. 10. Magnitude of dS/dt of one of the black holes from S ++0 . . Shown are three numerical resolutions, Alvi’s expres-sion as written (Eq. (11) of [78]), Alvi’s expression truncatedto leading order, and the CPY expression[80, 81] to 1PN or-der. The inset zooms closer to the high-resolution numericalcurve. and Yunes (hereafter CPY) [80], have recently extendedthese expressions to higher order in perturbation theory.Although CPY’s expressions in Ref. [80] are computedto 1.5PN beyond leading order (i.e. terms in dS/dt pro-portional to v and terms in dM/dt proportional to v ,where v = M/r is the PN expansion parameter), their1.5PN terms are incorrect and will be corrected soon [81];so here we will truncate CPY’s expressions to 1PN order.In our simulations we track the apparent horizons as afunction of time, and at frequent time intervals we mea-sure both the surface area and the spin of the horizons.The spin computation is carried out using the approx-imate Killing vector formalism of Cook, Whiting, andOwen [71, 72]. The mass of the black hole is then com-puted using Christodoulou’s formula. We compare ournumerical results to the analytic results of Alvi and CPY.To compare a black-hole mass or spin from a numer-ical simulation to that of a perturbative expression, thetwo quantities must be compared at the same event alongthe black hole trajectory. Although waveform quantitiesat future null infinity computed by numerical simulationsare routinely compared with waveforms computed by PNexpansions, it is not straightforward to compare near-zone quantities like black-hole masses and spins becauseof gauge ambiguities. Here we make two comparisons.The first compares quantities at the same numerical andperturbative t coordinate. The second assumes that theorbital angular velocity ω orb = dφ/dt of the black holein the numerical simulation can be equated with that ofthe perturbative expression. Note that in both the nu-2 x = (M ω orb ) | d S / dx | / M t/M Low resMedium resHigh resAlvi (0PN)AlviCPY (1PN)
FIG. 11. Magnitude of dS/dx , where x ≡ ( Mω orb ) / , ofone of the black holes from S ++0 . . The top horizontal axisshows t/M of the highest-resolution numerical simulation, forcomparison with values of x . merical and perturbative cases, the t coordinate becomesthe Minkowski t at infinity, and the φ coordinate is pe-riodic. Because of the approximate helical Killing vector d/dt + ω orb d/dφ , ω orb is approximately an angular veloc-ity at infinity. Therefore, one might hope that equatingthe perturbative and numerical ω orb yields better agree-ment than, e.g., equating the radial coordinate r of thesimulation with that of perturbation theory.Figure 10 compares the magnitude of dS/dt of oneof the black holes for S ++0 . with the expressions ofboth CPY and Alvi. We include numerical results forthree resolutions in Fig. 10 because the magnitude of dS/dt is extremely small and difficult to resolve. Indeed,the lowest resolution fails to resolve dS/dt until around t = 6000 M , when dS/dt grows to about 10 − M , andthe medium resolution fails to resolve dS/dt only slightlyearlier. Note that Alvi’s expression includes some 1.5PNterms, but ignores 1PN effects such as magnetic-typetidal perturbations and the difference between the globalPN time coordinate and the local time coordinate of aframe moving along with one of the black holes. There-fore, we plot both Alvi’s expression in its entirety, andAlvi’s expression truncated to lowest (0PN) order. TheCPY expression includes 0PN and 1PN terms. The CPYand Alvi expressions agree to 0PN order.Figure 10 shows overall excellent agreement betweenthe PN and numerical simulation results. All the pertur-bative curves agree within our our numerical error upto t ∼ M , but for t > M none of the per-turbative approximations agree with the numerical re-sult within numerical error. This disagreement at late t/M -10 -8 -6 -4 | d m / d t | Low resMedium resHigh resAlvi (0PN)AlviCPY (1PN) × -8 × -8 × -8 FIG. 12. Magnitude of dm/dt of one of the black holes from S ++0 . . Shown are three numerical resolutions, Alvi’s expres-sion as written (Eq. (11) of [78]), Alvi’s expression truncatedto leading order, and the CPY expression[80, 81] truncatedto 1PN order. The inset zooms closer to the high-resolutionnumerical curve. times is not surprising since all the perturbative expres-sions should lose accuracy shortly before merger. Wecan eliminate the time coordinate, a possible source ofgauge dependence, by instead plotting dS/dx versus x ,where x ≡ ( M ω orb ) / . This is shown in Fig. 11. Toobtain dS/dx from dS/dt and to obtain x from t , it isnecessary to have some function x ( t ). For the numericalcurves, this function is obtained from the numerical timecoordinate and the numerical orbital frequency. For theperturbative curves, this function is the PN expressionfor x ( t ) derived from Eq. (4.14) of Ref. [82]. Thus, allthe numerical curves in Figs. 10 and 11 are independentof any perturbative assumptions, and all the perturbativecurves in Figs. 10 and 11 are independent of the numer-ical data, except that the perturbative and numerical t coordinates are both represented by the same horizon-tal axis of Fig. 10, and the perturbative and numerical ω orb are both represented by the same horizontal axis ofFig. 11.In Fig. 11, the perturbative and numerical expressionsagree early in the inspiral, but not at late times; this isexpected because perturbative expressions become inac-curate for large x . Alvi’s full expression appears to agreewith the numerical simulations slightly better than theothers for small x , but that expression diverges from thenumerical result at larger x earlier than the others. Notethat Fig. 11 emphasizes late times because the frequencyincreases very rapidly with time.Figures 12 and 13 are similar to Figs. 10 and 11 exceptthat they show the change in Christodoulou mass insteadof the change in dimensionful spin. As was the case for3 x = (M ω orb ) | d m / dx | / M Low resMedium resHigh resAlvi (0PN)AlviCPY (1PN)
FIG. 13. Magnitude of dm/dx , where x ≡ ( Mω orb ) / , of oneof the black holes from S ++0 . . the spin comparisons, both the Alvi and CPY formulasagree well with each other and with the numerical resultearly in the inspiral, but do not agree at late times. Notethat since the derivative of the mass is smaller (by afactor of v in PN) than the derivative of the spin, dm/dt is more difficult to resolve numerically than dS/dt , asseen by the larger numerical errors in Figs. 12 and 13compared with the numerical errors in Figs. 10 and 11. B. Orbital hangup
During a BBH inspiral, the orbital frequency ω orb secularly evolves along with the black-hole masses andspins. For equal-mass binaries with equal spins aligned(or antialigned) with the orbital angular momentum, thenumber of orbits until merger increases as a function of S · L . Damour [83] observed this effect, today commonlycalled “orbital hangup”, in an effective-one-body modelof the holes’ motion; the effect is a consequence of post-Newtonian spin-orbit coupling [84]. Campanelli, Lousto,and Zlochower [85] first demonstrated orbital hangup innumerical simulations of merging BBHs.Instead of examining the number of orbits from the tra-jectories, we infer the number of orbits from the dominant (cid:96) = m = 2 mode of the emitted gravitational waves . Wedo this because it is easier to define a gauge-invariant Specifically, we extrapolate the gravitational waves measured ona series of concentric shells to r → ∞ , as discussed in detail inSec. IV C. χ = χ = χ = χ = χ = χ = χ = χ = M ω M ω FIG. 14. The evolution of the derivative of the gravitational-wave frequency ˙ ω = dω /dt , for simulations S ++0 . and S ++0 . and (for comparison) simulations S ++0 . [41], S ++0 . [41], S ++0 . [41], S ++0 . [41], S ++0 . [37], and S ++0 . [24]. time of merger from the waveforms than from the trajec-tories; specifically, we define the time of merger as thetime when the waveform amplitude is at a maximum.Let h ( t ) be the − Y spin-weighted spherical har-monic mode of the gravitational wave strain h ( t ), and let ω be the frequency of h ( t ). Figure 14 shows the timeevolution of dω /dt for simulations S ++0 . and S ++0 . . Forcomparison, we also show results for other simulationswith equal masses and equal spins aligned with the or-bital angular momentum [24, 37, 41]. Note that dω /dt is positive and steadily increasing: the frequency doesnot slow down or momentarily remain constant, as a lit-eral interpretation of the term “orbital hangup” mightsuggest.Figure 15 shows the gravitational-wave cycles accu-mulated between an initial gravitational-wave frequencyof M ω = 0 .
036 (i.e., an initial orbital frequency of
M ω orb = 0 . h peaks). Simulations S ++0 . and S ++0 . reveal that the or-bital hangup depends approximately linearly on the ini-tial spin χ , even at spins that are nearly extremal; how-ever, most of our simulations only agree with the linearfit to O (0 . C. Comparison with analytic approximants
We compare the gravitational waveforms from oursimulations to several analytic waveform approximants.The numerical waveforms were computed by perform-ing Regge-Wheeler-Zerilli extraction [86, 87] at a se-quence of radii between 100 M and 465 M , and thenextrapolating to I + using the open-source GWFrames N u m b e r o f c y c l e s χ -2-1012 r e s i du a l FIG. 15. The number of gravitational-wave cycles as a func-tion of the initial spin χ , measured after the initial relax-ation, for simulations S ++0 . and S ++0 . and (for compari-son) simulations S ++0 . [41], S ++0 . [41], S ++0 . [41], S ++0 . [41], S ++0 . [50], S ++0 . [37], and S ++0 . [24]. Upper panel:
The num-ber of gravitational-wave cycles of h accumulated between agravitational-wave frequency Mω = 0 .
036 and merger (i.e.,the time when the amplitude of h peaks). The dashed lineis a linear fit to the data. Lower panel:
Fractional difference(“residual”) between our results and the linear fit, with un-certainties for simulations except S ++0 . (which we ran at onlyone resolution) estimated as differences between medium andhigh numerical resolutions. software package [88–90]. The TaylorT1, TaylorT4,and TaylorT5 approximants were constructed usingthe PostNewtonian module in
GWFrames . The EOBapproximants were constructed using
SEOBNRv2 [49]from the LIGO Algorithm Library, with the func-tion
SimIMRSpinAlignedEOBWaveform modified to re-turn h ( t ). Physical parameters for the approximantswere taken from the highest resolution from each sim-ulation at the relaxation time. Because SEOBNRv2 is strictly valid only for non-precessing systems, andtherefore accepts only scalar values of the spins as in-put, it is not obvious what to input for the case of S . . . We pass the z -component of the spins into SimIMRSpinAlignedEOBWaveform . If instead we pass the To our knowledge, the
PostNewtonian module includes all termscurrently found in the literature. Non-spin terms are given upto 4.0 PN order for the binding energy [13, 91]; 3.5 PN withincomplete 4.0 PN information for the flux [13]; and 3.5 PN forthe waveform modes [92–94]. The spin-orbit terms are given to4.0 PN in the binding energy [95]; 3.5 PN with incomplete 4.0PN terms in flux [96]; and 2.0 PN in the waveform modes [89].Terms quadratic in spin are given to 2.0 PN order in the bindingenergy and flux [97, 98], and waveform modes [89, 97, 99]. − − − − − − t/M − − − − ∆ φ TaylorT5TaylorT4TaylorT1SEOBNRv2NR low resNR medium res
FIG. 16. Phase differences ∆ φ of h as a function of retardedtime before merger for S ++0 . . Shown are differences betweenthe highest numerical resolution and several analytic approx-imants. Differences between the highest numerical resolutionand other numerical resolutions are shown for comparison.The waveforms are aligned in the time interval delimited bythe black triangles. − − − − − t/M − − − − ∆ φ TaylorT5TaylorT4TaylorT1SEOBNRv2NR low resNR medium res
FIG. 17. Phase differences ∆ φ of h between numericaland approximant data for S ++0 . . Labels are the same as forFig. 16. spin magnitudes, we see larger disagreements between theEOB and numerical waveforms for S . . , likely due to achange in the strength of spin-orbit coupling. We will seebelow that non-precessing EOB agrees remarkably wellwith S . . despite the mild precession of this simulation.In Figs. 16, 17, and 18, we show for S ++0 . , S ++0 . , and S . . (respectively) the phase difference ∆ φ of h be-tween the highest numerical resolution and the PN andEOB approximants. We also include ∆ φ between thehighest numerical resolution and other numerical reso-lutions for comparison. To compute ∆ φ , we first aligneach waveform with the highest resolution numerical-relativity (NR) waveform using the procedure prescribed5 − − − − − t/M − − − − ∆ φ TaylorT5TaylorT4TaylorT1SEOBNRv2NR low resNR medium res
FIG. 18. Phase differences ∆ φ of h between numerical andapproximant data for S . . . Labels are the same as for Fig. 16.Note that the Taylor models include precession but SEOBNRv2 does not. However, the precession of S . . is mild so the nu-merical waveform still agrees reasonably well with SEOBNRv2 . in Ref. [100]: we find the time offset δt and phase off-set δφ that minimize Φ( δt, δφ ), a measure of the phasedifference in h , given byΦ( δt, δφ ) ≡ (cid:90) t t [ φ a ( t ) − φ b ( t + δt ) + δφ ] dt, (21)where δφ can be computed analytically from δtδφ ( δt ) = 1 t − t (cid:90) t t [ φ a ( t ) − φ b ( t + δt )] dt. (22)The alignment interval t ∈ [ t , t ] is the same for all com-parisons with a particular simulation. The lower bound t is chosen such that the junk radiation has left the com-putational domain for all numerical resolutions, specifi-cally t = max[ t + 3( t relax − t )], where t is the timeat the beginning of the waveform. The upper bound t is chosen such that the gravitational-wave frequencychanges by at least 10% during the interval [ t , t ], assuggested in Ref. [101].We have also computed ∆ φ with a few other alignmentmethods, including the three-dimensional minimizationof complex h differences in Ajith et al. 2008 (Eq. 4.9in Ref. [102]) and the four-dimensional minimization overtime and frame-rotation degrees of freedom in Boyle 2013(Eq. 22 in Ref. [88]). We have found that our results arequalitatively independent of alignment method.The TaylorT family of PN approximants shows thelargest discrepancy with our highest numerical resolu-tion. Phase errors between PN and NR waveforms growto several radians before the merger in every case. Thesmallest phase errors outside the alignment interval oc-cur for S . . , which is likely a consequence of the smallerblack hole having a moderate spin. We find the best agreement with TaylorT1, in contrast to PN compar-isons for other nearly extremal systems [37], which foundthe best agreement with TaylorT4 for spins aligned withthe orbital angular momentum; note that the PN wave-forms considered in Ref. [37] include fewer higher-orderPN terms than we do here. This is further evidence thatagreement with a particular PN approximant in the Tay-lorT family depends sensitively on the PN order. Agree-ment with a particular PN approximant also depends onthe parameters of the simulation (e.g., Ref. [37]).The EOB approximant performs significantly betterthan the PN approximants for S ++0 . and S ++0 . , whichis impressive considering that the parameters of thesewaveforms are outside the range in which SEOBNRv2 wascalibrated to NR. Only about 5 radians of phase error isaccumulated in S ++0 . and S ++0 . . Phase error increasesto a little over 10 radians in S . . , but this case is pre-cessing, and SEOBNRv2 is only valid for non-precessingsystems. However, the precession is mild (cf. Figs. 8and 9), which could account for the relatively good agree-ment.The analytic approximants show much larger ∆ φ atearly times for S . . (see Fig. 18) than for S ++0 . and S ++0 . . We conjecture that this is due to the relativelylarge eccentricity of S . . (see Table I), whereas the PNand EOB models used here are non-eccentric. Note thatwe use precessing PN models for comparing to S . . .The phase errors between numerical waveforms com-puted at different resolutions are convergent. Because ofthe rapid convergence, the difference between the twohighest numerical resolutions represents the numericalerror in the second highest resolution; to determine thenumerical error of the highest resolution waveform, wewould need to perform a simulation at an even higherresolution. As a conservative estimate of the numeri-cal error of the highest resolution waveform, we use thedifference between the two highest-resolution waveformsas an upper bound. The upper bound of the numericalphase error of the highest resolution simulation, com-puted in this way, is thus about 0.2 radians for S ++0 . and S ++0 . and about 1 radian for S . . .In Figs. 16, 17, and 18, the larger numerical phaseerrors in the lower resolutions of S ++0 . and S . . areexpected, because these simulations use a larger spac-ing in AMR truncation error tolerance as described inSec. III B. The larger spacing increases relative phaseerrors between successive numerical resolutions. Never-theless, our comparisons show that numerical errors aremuch smaller than the errors in the PN and EOB wave-forms for systems with nearly extremal black holes, in-dicating that these numerical waveforms will be usefulfor calibrating and extending the regime of validity forapproximate waveforms. Note that
SEOBNRv2 was calibrated by minimizing unfaithfulnessrather than phase error; it is possible to have relatively largephase errors even when the unfaithfulness is small [49]. V. CONCLUSION
We have presented improved methods for simulatingthe binary evolution of nearly extremal black holes, i.e.,black holes with spins above the Bowen-York limit of χ = 0 .
93. These techniques enable robust simulationsin the portion of BBH parameter space where the blackholes have very large spins. Because nearly extremalblack holes might exist in astrophysical binaries, thesesimulations will be important for helping to maximizewhat we can learn from gravitational-wave experiments.We have applied our new methods to carry out thefirst unequal-mass, mildly-precessing BBH simulationcontaining a nearly extremal black hole, and to extendaligned-spin BBH simulations to spin magnitudes thatbegin to approach the Novikov-Thorne limit of χ = 0 . χ in anextremely nonlinear fashion, we find that the numberof orbits starting from a chosen orbital frequency (i.e.,the orbital hangup) scales approximately linearly with χ . Finally, after demonstrating numerical convergence,we have found that our numerical waveforms agree with SEOBNRv2 much better than with TaylorT PN approxi-mants, even though the parameters for these simulationsare outside the range in which
SEOBNRv2 was calibrated.However, even the
SEOBNRv2 waveforms disagree with our numerical waveforms by more than our numerical trun-cation error. This indicates that these simulations aresufficiently accurate to validate and further improve an-alytical waveform approximants for future gravitational-wave observations. How significant these improvementswill be for Advanced LIGO is the subject of future work.
ACKNOWLEDGMENTS
We are grateful to Eric Poisson, Nicolas Yunes, andKaterina Chatziioannou for detailed discussions aboutperturbative expressions for tidal torquing and aboutthe problems inherent in comparing numerical and post-Newtonian expressions for near-field quantities. Wethank Alessandra Buonanno and Sebastiano Bernuzzi forhelpful discussions. Simulations used in this work werecomputed with SpEC [47]. This work was supported inpart by the Sherman Fairchild Foundation; NSF grantsPHY-1440083 and AST-1333520 at Caltech, NSF grantsPHY-1306125 and AST-1333129 at Cornell, and NSFgrant PHY-1307489 at California State University Fuller-ton; a 2013–2014 California State University FullertonJunior Faculty Research Grant. Computations were per-formed on the Zwicky cluster at Caltech, which is sup-ported by the Sherman Fairchild Foundation and by NSFaward PHY-0960291; on the NSF XSEDE network undergrant TG-PHY990007N; on the Orca cluster supportedby NSF award NSF-1429873 and by California State Uni-versity Fullerton; and on the GPC supercomputer at theSciNet HPC Consortium [103]. SciNet is funded by: theCanada Foundation for Innovation under the auspices ofCompute Canada; the Government of Ontario; OntarioResearch Fund–Research Excellence; and the Universityof Toronto. [1] G. M. Harry (LIGO Scientific Collaboration), Class.Quant. Grav. , 084006 (2010).[2] The Virgo Collaboration, “Advanced Virgo Baseline De-sign,” (2009), [VIR-0027A-09].[3] The Virgo Collaboration, “Advanced Virgo TechnicalDesign Report,” (2012), [VIR-0128A-12].[4] K. Somiya and the KAGRA Collaboration, Class. Quan-tum Grav. , 124007 (2012).[5] L. Gou, J. E. McClintock, M. J. Reid, J. A. Orosz, J. F.Steiner, R. Narayan, J. Xiang, R. A. Remillard, K. A.Arnaud, and S. W. Davis, Astrophys. J. , 85 (2011),arXiv:1106.3690 [astro-ph.HE].[6] A. Fabian, D. Wilkins, J. Miller, R. Reis, C. Reynolds, et al. , MNRAS , 217 (2012), arXiv:1204.5854 [astro-ph.HE].[7] L. Gou, J. E. McClintock, R. A. Remillard, J. F. Steiner,M. J. Reid, et al. , Astrophys.J. , 29 (2014).[8] J. E. McClintock, R. Shafee, R. Narayan, R. A. Remil-lard, S. W. Davis, and L.-X. Li, Astrophys. J. , 518(2006).[9] J. Miller, C. Reynolds, A. Fabian, G. Miniutti, and L. Gallo, Astrophys.J. , 900 (2009), arXiv:0902.2840[astro-ph.HE].[10] D. Walton, E. Nardini, A. Fabian, L. Gallo, andR. Reis, MNRAS , 2901 (2013), arXiv:1210.4593[astro-ph.HE].[11] J. E. McClintock, R. Narayan, and J. F. Steiner,Space Sci.Rev. , 295 (2014), arXiv:1303.1583 [astro-ph.HE].[12] C. S. Reynolds, (2013), arXiv:1302.3260 [astro-ph.HE].[13] L. Blanchet, Living Rev.Rel. , 4 (2006).[14] F. Pretorius, Phys. Rev. Lett. , 121101 (2005),arXiv:gr-qc/0507014 [gr-qc].[15] M. Campanelli, C. Lousto, P. Marronetti, and Y. Zlo-chower, Phys. Rev. Lett. , 111101 (2006), arXiv:gr-qc/0511048 [gr-qc].[16] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz,and J. van Meter, Phys. Rev. Lett. , 111102 (2006),arXiv:gr-qc/0511103 [gr-qc].[17] J. Centrella, J. G. Baker, B. J. Kelly, and J. R. vanMeter, Rev. Mod. Phys. , 3069 (2010).[18] H. P. Pfeiffer, Class. Quant. Grav. , 124004 (2012), arXiv:1203.5166 [gr-qc].[19] M. Hannam, (2013), arXiv:1312.3641 [gr-qc].[20] A. Le Tiec, Int. J. Mod. Phys. D , 1430022 (2014),arXiv:1408.5505 [gr-qc].[21] P. Ajith, M. Boyle, D. A. Brown, B. Brugmann, L. T.Buchman, et al. , Class. Quantum Grav. , 124001(2012).[22] I. Hinder et al. (The NRAR Collaboration), Clas-sical and Quantum Gravity , 025012 (2014),arXiv:1307.5307 [gr-qc].[23] L. Pekowsky, R. O’Shaughnessy, J. Healy, andD. Shoemaker, Phys. Rev. D , 024040 (2013),arXiv:1304.3176 [gr-qc].[24] A. H. Mroue, M. A. Scheel, B. Szilagyi, H. P. Pfeiffer,M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace,S. Ossokine, N. W. Taylor, A. Zenginoglu, L. T. Buch-man, T. Chu, E. Foley, M. Giesler, R. Owen, andS. A. Teukolsky, Phys. Rev. Lett. , 241104 (2013),arXiv:1304.6077 [gr-qc].[25] J. Healy, C. O. Lousto, and Y. Zlochower, Phys. Rev.D , 104052 (2014), arXiv:1406.7295 [gr-qc].[26] J. Clark, L. Cadonati, J. Healy, I. S. Heng, J. Logue, et al. , (2014), arXiv:1406.5426 [gr-qc].[27] J. M. Bowen, Gen. Relativ. Gravit. , 227 (1979).[28] J. M. Bowen and J. W. York, Jr., Phys. Rev. D , 2047(1980).[29] S. Brandt and B. Br¨ugmann, Phys. Rev. Lett. , 3606(1997).[30] J. W. York, Jr., in Essays in General Relativity , editedby F. J. Tipler (Academic, New York, 1980) pp. 39–58.[31] A. Garat and R. H. Price, Phys. Rev. D , 124011(2000).[32] J. A. Valiente Kroon, Phys. Rev. Lett. , 041101(2004).[33] G. B. Cook and J. W. York, Jr., Phys. Rev. D , 1077(1990).[34] S. Dain, C. O. Lousto, and R. Takahashi, Phys. Rev.D , 104038 (2002).[35] M. Hannam, S. Husa, and N. O. Murchadha, Phys.Rev. D , 124007 (2009).[36] G. Lovelace, M. A. Scheel, and B. Szil´agyi, Phys. Rev.D , 024010 (2011), arXiv:1010.2777 [gr-qc].[37] G. Lovelace, M. Boyle, M. A. Scheel, andB. Szil´agyi, Class. Quant. Grav. , 045003 (2012),arXiv:arXiv:1110.2229 [gr-qc].[38] G. Lovelace, M. D. Duez, F. Foucart, L. E. Kidder, H. P.Pfeiffer, M. A. Scheel, and B. Szil´agyi, Class. QuantumGrav. , 135004 (2013), arXiv:1302.6297 [gr-qc].[39] I. D. Novikov and K. S. Thorne, “Black holes,” (Gordonand Breach, New York, 1973) p. 343.[40] K. S. Thorne, Astrophys. J. , 507 (1974).[41] D. A. Hemberger, G. Lovelace, T. J. Loredo, L. E. Kid-der, M. A. Scheel, B. Szil´agyi, N. W. Taylor, andS. A. Teukolsky, Phys. Rev. D , 064014 (2013),arXiv:1305.5991 [gr-qc].[42] L. Rezzolla et al. , Astrophys. J. , 1422 (2008),arXiv:0708.3999 [gr-qc].[43] M. Hannam, S. Husa, F. Ohme, D. M¨uller, andB. Br¨ugmann, Phys. Rev. D , 124008 (2010),arXiv:1007.4789.[44] P. Marronetti, W. Tichy, B. Br¨ugmann, J. Gonz´alez,and U. Sperhake, Phys. Rev. D , 064010 (2008).[45] S. Dain, C. O. Lousto, and Y. Zlochower, Phys. Rev.D , 024039 (2008), arXiv:0803.0351v2 [gr-qc]. [46] I. Ruchlin, J. Healy, C. O. Lousto, and Y. Zlochower,(2014), arXiv:1410.8607 [gr-qc].[47] .[48] G. Lovelace, R. Owen, H. P. Pfeiffer, and T. Chu, Phys.Rev. D , 084017 (2008).[49] A. Taracchini, A. Buonanno, Y. Pan, T. Hinderer,M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace,A. H. Mroue, H. P. Pfeiffer, M. A. Scheel, B. Szilagyi,and A. Zenginoglu, Phys.Rev. D89 , 061502 (2014),arXiv:1311.2544 [gr-qc].[50] G. Lovelace, M. A. Scheel, R. Owen, M. Giesler,R. Katebi, B. Szil´agyi, T. Chu, N. Demos, D. A.Hemberger, L. E. Kidder, H. P. Pfeiffer, and N. Af-shari, (2014), submitted to Class. Quantum Grav.,arXiv:1411.7297 [gr-qc].[51] H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A.Teukolsky, Comput. Phys. Commun. , 253 (2003).[52] M. Caudill, G. B. Cook, J. D. Grigsby, and H. P. Pfeif-fer, Phys. Rev. D , 064011 (2006), gr-qc/0605053.[53] J. W. York, Phys. Rev. Lett. , 1350 (1999).[54] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom,G. Lovelace, and M. A. Scheel, Class. Quantum Grav. , S59 (2007), gr-qc/0702106.[55] A. Buonanno, L. E. Kidder, A. H. Mrou´e, H. P. Pfeiffer,and A. Taracchini, Phys. Rev. D , 104034 (2011),arXiv:1012.1549 [gr-qc].[56] A. H. Mrou´e and H. P. Pfeiffer, (2012), arXiv:1210.2958[gr-qc].[57] H. Friedrich, Commun. Math. Phys. , 525 (1985).[58] D. Garfinkle, Phys. Rev. D , 044029 (2002).[59] F. Pretorius, Class. Quantum Grav. , 425 (2005).[60] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, andO. Rinne, Class. Quantum Grav. , S447 (2006).[61] L. Lindblom and B. Szil´agyi, Phys. Rev. D , 084019(2009), arXiv:0904.4873.[62] M. W. Choptuik and F. Pretorius, Phys. Rev. Lett. ,111101 (2010), arXiv:0908.1780 [gr-qc].[63] B. Szil´agyi, L. Lindblom, and M. A. Scheel, Phys. Rev.D , 124010 (2009), arXiv:0909.3557 [gr-qc].[64] B. Szil´agyi, Int.J.Mod.Phys. D23 , 1430014 (2014),arXiv:1405.3693 [gr-qc].[65] M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D.Matthews and H. P. Pfeiffer, Phys. Rev. D , 024003(2009), arXiv:gr-qc/0810.1767.[66] D. A. Hemberger, M. A. Scheel, L. E. Kidder,B. Szil´agyi, G. Lovelace, N. W. Taylor, and S. A.Teukolsky, Class. Quantum Grav. , 115001 (2013),arXiv:1211.6079 [gr-qc].[67] S. Ossokine, L. E. Kidder, and H. P. Pfeiffer,arXiv:1304.3067 (2013), arXiv:1304.3067 [gr-qc].[68] O. Rinne, Class. Quantum Grav. , 6275 (2006).[69] O. Rinne, L. Lindblom, and M. A. Scheel, Class. Quan-tum Grav. , 4053 (2007).[70] C. Gundlach, Phys. Rev. D , 863 (1998).[71] G. B. Cook and B. F. Whiting, Phys. Rev. D ,041501(R) (2007).[72] R. Owen, Topics in Numerical Relativity: The peri-odic standing-wave approximation, the stability of con-straints in free evolution, and the spin of dynamicalblack holes , Ph.D. thesis, California Institute of Tech-nology (2007).[73] M. A. Scheel, H. P. Pfeiffer, L. Lindblom, L. E. Kidder,O. Rinne, and S. A. Teukolsky, Phys. Rev. D , 104006(2006). [74] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mrou´e, H. P.Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky,Phys. Rev. D , 124038 (2007), arXiv:0710.0158 [gr-qc].[75] L. T. Buchman, H. P. Pfeiffer, M. A. Scheel,and B. Szil´agyi, Phys. Rev. D , 084033 (2012),arXiv:1206.3015 [gr-qc].[76] T. Chu, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D , 124051 (2009), arXiv:0909.1313 [gr-qc].[77] .[78] K. Alvi, Phys. Rev. D , 104020 (2001).[79] R. Owen, J. Brink, Y. Chen, J. D. Kaplan, G. Lovelace,K. D. Matthews, D. A. Nichols, M. A. Scheel, F. Zhang,A. Zimmerman, and K. S. Thorne, Phys. Rev. Lett. , 151101 (2011).[80] K. Chatziioannou, E. Poisson, and N. Yunes, Phys.Rev. D , 044022 (2013), arXiv:1211.1686 [gr-qc].[81] K. Chatziioannou, E. Poisson, and N. Yunes, In prepa-ration.[82] L. E. Kidder, Phys. Rev. D , 821 (1995).[83] T. Damour, Phys. Rev. D , 124013 (2001), arXiv:gr-qc/0103018 [gr-qc].[84] L. E. Kidder, Phys. Rev. D52 , 821 (1995), arXiv:gr-qc/9506022.[85] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.Rev. D , 041501(R) (2006), gr-qc/0604012.[86] O. Sarbach and M. Tiglio, Phys. Rev. D , 084016(2001).[87] O. Rinne, L. T. Buchman, M. A. Scheel, and H. P.Pfeiffer, Class. Quantum Grav. , 075009 (2009).[88] M. Boyle, Phys. Rev. D , 104006 (2013).[89] M. Boyle, L. E. Kidder, S. Ossokine, and H. P. Pfeiffer,(2014), arXiv:1409.4431, arXiv:1409.4431.[90] S. Ossokine, M. Boyle, L. E. Kidder, A. H. Mrou´e, H. P.Pfeiffer, M. A. Scheel, and B. Szil´agyi, (2014), in prepa-ration.[91] D. Bini and T. Damour, “Analytical determination of the two-body gravitational interaction potentialat the 4th post-Newtonian approximation,” (2013),arXiv:1305.4884 [gr-qc].[92] L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, Class.Quant. Grav. , 165003 (2008), arXiv:0802.1249 [gr-qc].[93] G. Faye, S. Marsat, L. Blanchet, and B. R. Iyer, Class.Quant. Grav. , 175004 (2012).[94] G. Faye, L. Blanchet, and B. R. Iyer, “Non-linearmultipole interactions and gravitational-wave octupolemodes for inspiralling compact binaries to third-and-a-half post-Newtonian order,” (2014), arXiv:1409.3546[gr-qc].[95] A. Boh´e, S. Marsat, G. Faye, and L. Blanchet, Class.Quant. Grav. , 075017 (2013).[96] S. Marsat, L. Blanchet, A. Boh´e, and G. Faye, “Grav-itational waves from spinning compact object binaries:New post-Newtonian results,” (2013), arXiv:1312.5375[gr-qc].[97] C. M. Will and A. G. Wiseman, Phys. Rev. D , 4813(1996).[98] K. Arun, A. Buonanno, G. Faye, and E. Ochsner, Phys.Rev. D , 104023 (2009), arXiv:0810.5336 [gr-qc].[99] A. Buonanno, G. Faye, and T. Hinderer, Phys. Rev. D , 044009 (2013).[100] M. Boyle, A. Buonanno, L. E. Kidder, A. H. Mrou´e,Y. Pan, et al. , Phys. Rev. D , 104020 (2008),arXiv:0804.4184 [gr-qc].[101] I. MacDonald, S. Nissanke, and H. P. Pfeiffer, Class.Quantum Grav. , 134002 (2011), arXiv:1102.5128 [gr-qc].[102] P. Ajith, S. Babak, Y. Chen, M. Hewitson, B. Krishnan, et al. , Phys. Rev. D , 104017 (2008), arXiv:0710.2335[gr-qc].[103] C. Loken, D. Gruner, L. Groer, R. Peltier, N. Bunn,M. Craig, T. Henriques, J. Dempsey, C.-H. Yu, J. Chen,L. J. Dursi, J. Chong, S. Northrup, J. Pinto, N. Knecht,and R. V. Zon, J. Phys.: Conf. Ser.256