Intermediate-mass-ratio black hole binaries II: Modeling Trajectories and Gravitational Waveforms
Hiroyuki Nakano, Yosef Zlochower, Carlos O. Lousto, Manuela Campanelli
IIntermediate-mass-ratio black hole binaries II:Modeling Trajectories and Gravitational Waveforms
Hiroyuki Nakano, Yosef Zlochower, Carlos O. Lousto, Manuela Campanelli
Center for Computational Relativity and Gravitation,and School of Mathematical Sciences, Rochester Institute of Technology,85 Lomb Memorial Drive, Rochester, New York 14623
We revisit the scenario of small-mass-ratio ( q ) black-hole binaries; performing new, more accurate,simulations of mass ratios 10:1 and 100:1 for initially nonspinning black holes. We propose fittingfunctions for the trajectories of the two black holes as a function of time and mass ratio (in the range1 / ≤ q ≤ /
10) that combine aspects of post-Newtonian trajectories at smaller orbital frequenciesand plunging geodesics at larger frequencies. We then use these trajectories to compute waveformsvia black hole perturbation theory. Using the advanced LIGO noise curve, we see a match of ∼ . (cid:96), m ) = (2 ,
2) mode between the numerical relativity and perturbative waveforms.Nonleading modes have similarly high matches. We thus prove the feasibility of efficiently generatinga bank of gravitational waveforms in the intermediate-mass-ratio regime using only a sparse set offull numerical simulations.
PACS numbers: 04.25.dg, 04.30.Db, 04.25.Nx, 04.70.Bw
I. INTRODUCTION
Numerical relativity (NR) has come a long way sincethe breakthroughs of 2005 [1–3] that allowed, for the firsttime, long-term evolution of black hole binaries (BHBs).Among NR’s significant achievements are its contribu-tions towards the modeling of astrophysical gravitationalwave sources that will be relevant for the first direct de-tection and parameter estimation by gravitational waveobservatories [4]. NR has also made contributions to themodeling of astrophysical sources, notably with the dis-covery of very large recoil velocities [5–8], and the appli-cation of the numerical techniques to combined systemsof black holes and neutron stars [9]. More mathematicalaspects of relativity have also recently been investigated,including the evolution of N-black holes [10], the explo-ration of the no-hair theorem [11, 12], and cosmic [13]and topological censorship [14], as well as BHBs in di-mensions higher than four [15].Among the remaining challenges are the exploration ofthe extremes of the BHB parameter space. The currentstate of the art simulations can simulate BHBs with massratios as small as q = 1 /
100 [16, 17] and highly spinningBHBs with intrinsic spins α = S H /M H up to (at least)0 .
97 [18]. Currently these runs are very costly and it ishard to foresee the possibility of completely covering theparameter space densely enough for match filtering thedata coming from advanced laser interferometric detec-tors by the time they become operational. It is thereforeimperative to develop interpolation techniques that allowfor astrophysical parameter estimations, at a reasonablelevel of accuracy, based on a sparse set of numerical sim-ulations.In this spirit, we designed a set of prototypical runsfor initially non spinning BHBs with small mass ratios q in the range 0 . ≤ q ≤ .
01. Since we expect that,for small enough mass ratios, this BHB system will be describable by perturbation theory, we compare the fullnumerical waveforms with those produced by perturba-tive evolutions via the Regge-Wheeler [19] and Zerilli [20]equations, supplemented by linear corrections in the spinof the large black hole [21]. The key ingredients pertur-bation theory needs is the relative trajectory of the smallblack hole with respect to the larger one and the back-ground mass and spin. In our previous tests we usedthe full numerical tracks and proved the perturbativewaveforms agree reasonably well with the full numeri-cal ones [21]. In this paper we develop a model with freefitting parameters for these trajectories based on post-Newtonian and geodesic input, and fit to full numericaltracks. For these fits, we use full numerical evolutions of q = 1 /
10 and q = 1 /
100 BHBs [16, 21] and perform new,more accurate simulations. We compare the waveformsfor the modeled tracks with the full numerical waveformsto confirm matching agreement within a fraction of a per-cent. Hence this method paves the way for further gen-eralizations and simulations to provide an approximate,yet accurate, bank for waveforms for second generationgravitational wave detectors.This paper is organized as follows. In Sec. II we re-view the numerical methodology to perform the smallmass ratio runs. In Sec. III we describe the runs andresults. In Sec. IV we describe the track modeling onpost-Newtonian expansions and fits to the full numericalresults. In Sec. V we give the results of using those tracksto generate perturbative waveforms and compares themwith those extracted from the full numerical evolutions.We discuss the consequences and future application ofthis techniques in Sec. VI. We also include an appendix Ato briefly discuss perturbative theory in numerical coor-dinates (1+log, trumpet coordinates). a r X i v : . [ g r- q c ] A ug II. NUMERICAL RELATIVITY TECHNIQUES
To compute the numerical initial data, we use thepuncture approach [22] along with the
TwoPunc-tures [23] thorn. In this approach the 3-metric on theinitial slice has the form γ ab = ( ψ BL + u ) δ ab , where ψ BL is the Brill-Lindquist conformal factor, δ ab is theEuclidean metric, and u is (at least) C on the punc-tures. The Brill-Lindquist conformal factor is given by ψ BL = 1+ (cid:80) ni =1 m pi / (2 | (cid:126)r − (cid:126)r i | ) , where n is the total num-ber of ‘punctures’, m pi is the mass parameter of puncture i ( m pi is not the horizon mass associated with puncture i ),and (cid:126)r i is the coordinate location of puncture i . We evolvethese black-hole-binary data-sets using the LazEv [24]implementation of the moving puncture approach [2, 3]with the conformal function W = √ χ = exp( − φ ) sug-gested by Ref. [25]. For the runs presented here, we usecentered, eighth-order finite differencing in space [10] anda fourth-order Runge-Kutta time integrator. (Note thatwe do not upwind the advection terms.)Our code uses the Cactus / EinsteinToolkit [26, 27]infrastructure. We use the
Carpet [28] mesh refinementdriver to provide a “moving boxes” style of mesh refine-ment. In this approach refined grids of fixed size are ar-ranged about the coordinate centers of both holes. The
Carpet code then moves these fine grids about the com-putational domain by following the trajectories of the twoblack holes.We use
AHFinderDirect [29] to locate apparenthorizons. We measure the magnitude of the horizon spinusing the Isolated Horizon algorithm detailed in Ref. [30].Note that once we have the horizon spin, we can calculatethe horizon mass via the Christodoulou formula m H = (cid:113) m + S / (4 m ) , (1)where m irr = (cid:112) A/ (16 π ) and A is the surface area of thehorizon. We measure radiated energy, linear momentum,and angular momentum, in terms of ψ , using the for-mulae provided in Refs. [31, 32]. However, rather thanusing the full ψ , we decompose it into (cid:96) and m modesand solve for the radiated linear momentum, droppingterms with (cid:96) ≥
5. The formulae in Refs. [31, 32] arevalid at r = ∞ . Typically, we would extract the radi-ated energy-momentum at finite radius and extrapolateto r = ∞ . However, for the smaller mass ratios examinedhere, noise in the waveform introduces spurious effectsthat make these extrapolations inaccurate. We there-fore use the average of these quantities extracted at radii r = 70, 80, 90, 100 and use the difference between thesequantities at different radii as a measure of the error.We extrapolate the waveform to r → ∞ using the per-turbative formula [21]lim r →∞ [ r ψ (cid:96)m ( r, t )]= (cid:20) r ψ (cid:96)m ( r, t ) − ( (cid:96) − (cid:96) + 2)2 (cid:90) t dt ψ (cid:96)m ( r, t ) (cid:21) r = r Obs + O ( R − ) , (2) where r Obs is the approximate areal radius of the sphere R Obs = const [Add a factor (1 / − M/r ) multiplying thesquare bracket to correct for a difference in normaliza-tion between the Psikadelia (numerical) and Kinnersleytetrads at large distances.] We have found that this for-mula gives reliable extrapolations for R Obs > ∼ M .Recent Cauchy-Characteristic extraction (CCE) stud-ies [33] (see also Ref. [34]) showed that this perturba-tive extrapolation formula can be more accurate thana linear extrapolation of the waveforms at finite radiusto infinite resolution (provided that the observer R is inthe far zone). Those studies compared the extrapolatedwaveforms with the gauge invariant waveform on I + ob-tained using a nonlinear characteristic evolution from afinite radius to I + along outgoing null slices. The au-thors of Ref. [33] measured the errors in rψ for an equal-mass BHB simulation when extracting at R = 50 M and R = 100 M , the corresponding extrapolation to ∞ , aswell as the error in rψ obtained by applying the per-turbative extrapolation formula (2) to the R = 100 M waveform. The errors in the perturbative extrapolationwhere the smallest over the entire waveform (both in am-plitude and phase).We note that multi-patch [35] and pseudospectral [36]techniques allow extraction radii very far from the source,leading to very small extrapolation errors. A. Gauge
We obtain accurate, convergent waveforms and horizonparameters by evolving this system in conjunction witha modified 1+log lapse and a modified Gamma-drivershift condition [2, 37], and an initial lapse α ( t = 0) =2 / (1 + ψ BL ). The lapse and shift are evolved with( ∂ t − β i ∂ i ) α = − αK, (3a) ∂ t β a = (3 / a − η ( x a , t ) β a , (3b)where different functional dependences for η ( x a , t ) havebeen proposed in Refs. [24, 38–42]. Here we use a modi-fication of the form proposed in Ref. [39], η ( x a , t ) = R (cid:112) ∂ i W ∂ j W ˜ γ ij (1 − W a ) b , (4)where we chose R = 1 .
31. The above gauge conditionis inspired by, but differs from Ref. [39] between the BHsand in the outer zones when a (cid:54) = 1 and b (cid:54) = 2. Oncethe conformal factor settles down to its asymptotic ψ = C/ √ r + O (1) form near the puncture, η will have theform η = ( R /C )(1 + b ( r/C ) a ) near the puncture and η = R r b − M/ ( aM ) b as r → ∞ . In practice we used a = 2 and b = 2, which reduces η by a factor of 4 atinfinity when compared to the original version of thisgauge proposed by Ref. [39]. We note that if we set b = 1then η will have a 1 /r falloff at r = ∞ as suggested byRef. [41]. Our tests indicate that the choices ( a = 2, TABLE I: Initial data parameters for the numerical simu-lations. Note that the q = 1 /
15 simulations are older andused a CFL factor twice as large. Here m p and m p are thetwo puncture mass parameters, the puncture were located at( x , ,
0) and ( x , ,
0) with momentum ± ( P r , P t ,
0) and zerospin. The measured horizon masses and total ADM mass arealso provided.Param q = 1 / q = 1 / q = 1 / m p m p x x -0.7531758055 -0.4438775230 -0.04743736368 P t P r -0.000168519 -0.000160518 -0.00001026521884 m H m H M ADM b = 1) and ( a = 1 , b = 1) lead to more noise in thewaveform than ( a = 2 , b = 2). III. SIMULATIONS AND RESULTS
The initial data parameters for the q = 1 / q = 1 / q = 1 /
100 simulations are given in Table I. Note thatthe q = 1 /
15 simulations are older and suffer from themass loss error discussed below.We used a base (coarsest) resolution of h = 4 M with11 levels of refinement for q = 1 /
10 simulations and 15levels of refinement for q = 1 /
100 for the low resolu-tion simulations. The outer boundaries were at 400 M .The higher resolution simulations were all based on thesegrids, but with correspondingly more gridpoints per level.The grid structure was chosen by studying the behaviorof the background potential for the propagation of per-turbations [16]. In the appendix A we show this poten-tial, both in the isotropic coordinates of the initial dataand final ’trumpet’ coordinates.We previously evolved a set of q = 1 / q = 1 /
15 [21],and q = 1 /
100 [16] BHBs using the standard choice ofCourant-Friedrichs-Lewy (CFL) factor dt/h ∼ .
5. Al-though we found the waveform at lower resolution appearto converge, an unphysical mass loss led to incorrect dy-namics at later times (see Fig. 1). This, in turn, led tooscillations in the errors as a function of grid resolution h . We found that reducing the CFL factor significantlyreduces these unphysical effects [67]. In Figs. 1 and 2we plot the horizon mass as a function of time for threeresolutions using both the old and new CFL factor forboth the q = 1 /
10 and q = 1 /
100 simulations. Note themuch better conservation of the mass and the correspond-ing reduction in the lifetime of binary. Post-Newtoniantrajectory evolutions (see Fig. 3) indicate that the masslosses as small as 1 part in 10 are dynamically impor-tant, and the error introduced by this mass loss is sig-nificantly reduced by the new integration. The effects of ho r i z on m a ss ( m H / M ) h (dt/h=0.5)h /1.2 (dt/h=0.5)h /1.2 (dt/h=0.5)h (dt/h=0.25)h /1.2 (dt/h=0.25)h /1.2 (dt/h=0.25) FIG. 1: The horizon mass conservation using CFL factors of0 . .
25 for the q = 1 /
10 simulations. A factor of 10better results are obtained using a CFL factor of 0 .
25. Theeffect on the trajectory is significant because the mass lossis dynamically important (it significantly delays the merger).The mass changes post merger (sharp spikes near t ∼ M and t ∼ M ) are not dynamically important. the time integrator on the numerical error, and in par-ticular the mass loss and constraint violation errors, willbe the subject of an upcoming paper by Ponce, Lousto,and Zlochower.For the q = 1 /
10 nonspinning BHB, we performed sim-ulations with central resolutions of h = M/ h = M/ . M/ . M/ .
64, and M/ .
48, with 11levels of refinement in each case (the coarsest grid reso-lutions where h , h / . h / . h / . h / .
58, re-spectively). We note that the waveform showed oscilla-tions in error as a function of resolution. We are thereforeusing the most widely spaced (in resolution) runs for theconvergence plot. In particular, we use the h = M/ M/ .
9, and M/ .
48 simulations.In Fig. 4 we show the convergence of the amplitudeand phase of ψ for these three resolutions. Convergenceappears to break just at the start of the final plunge.Figure 5 shows the highest resolution waveform and theRichardson extrapolation (assuming a second-order er-ror). The extrapolation error becomes large just near theplunge, at a frequency of ω = 0 . /M . We also plot thephase deference between the highest resolution run andthe Richardson extrapolated waveform. The vertical lineshows the point when the frequency is ω = 0 . /M . Thephase error, up to t = 850 M is within 0 .
04 radian, butincreases exponentially to 1 . ω = 0 . /M .Determining convergence for the q = 1 /
100 simula-tions proved challenging for two reasons. There was asmall random jump in the mass of the smaller BH be-tween resolutions that had a small effect on the trajec-tory. From Fig. 6 we can see a difference of one partin 10 in the small BH mass. The two medium reso-lution have lower masses, and therefore merge slightlylater than expected. The net effect is to confuse the or- −6 −5 −4 −3 −2 | d m / m | dt/h=0.5dt/h=0.25 FIG. 2: The horizon mass conservation using CFL factors of0 . .
25 for the q = 1 /
100 simulations. Here | δm | is theabsolute value of the difference between the expected horizonmass and the measured horizon mass for the smaller BH. Inthe case of q = 1 / q = 1 /
10, the horizonmass decreases (hence | δm | increases) with time when using aCFL factor of 0 .
5. A factor of 20 better results are obtainedusing a CFL factor of 0 .
25. The effect on the trajectory issignificant because the mass loss is dynamically important (itdelays the merger). The mass changes post merger (sharpspikes near t = 115 M and t = 172 M ) are not dynamicallyimportant. der when the mergers happens (we expected the higherresolutions runs to merge slightly later than the lowerresolutions runs). We see from Fig. 7 that the second-lowest resolution merges out of order with the other runs.Similarly, the third lowest resolution appears to mergelate, as well. Another source of confusion for our con-vergence study has to do with lower-order errors (notnecessarily associated with the mass) becoming impor-tant at higher resolutions. This leads to oscillations inthe error as a function of resolution that can cause themerger time to oscillate with resolution. Nevertheless,the actual deviations in the trajectories, for the higherresolutions, is small. The improvement over the origi-nal q = 1 /
100 simulation is substantial. The mass con-servation is a factor of 20 better. The old simulationshad an unphysical extra orbit inside the ISCO due tosignificant mass loss. Convergence measurements of thewaveform do not appear to suffer much from the abovementioned mass loss. However, due to low amplitude ofthe waveform compared to the grid noise (induced by therefinement boundaries), the differences in the waveformsat the various resolution were smaller than the grid-noisefluctuations. Figures 9 and 10 show the waveform andphase for the three highest resolutions. Note that thegridnoise apparent in Fig. 9 is due to reflections of theinitial data pulse at the refinement boundaries. For ourconfiguration, the imaginary part of the ( (cid:96) = 2 , m = 2)components of this pulse is small compared to the realpart. This leads to a significant reduction (about a factor P N o r b i t a l s epa r a t i on q=1/10 d M/M = 0q=1/10 d M/M =−1.8x10 −3 q=1/10 d M/M =−3.3x10 −4 FIG. 3: The 3.5PN orbital separations for the q = 1 /
10 bi-nary (the eccentricity is due to the fact that we use identicalparameters as the full numerical simulation, which are notquasi-circular PN parameters) for slightly different small BHmasses. A mass loss of 1 part in 10 is dynamically impor-tant, leading to a significantly delayed merger. Reducing themass change by a factor of 10 significantly reduces this error.TABLE II: Final remnant intrinsic spin, radiated energy, andrecoil velocity for the q = 1 / q = 1 /
15, and q = 1 / q = 1 / q = 1 / q = 1 / α . ± .
002 0 . ± .
006 0 . ± . E rad /M . ± . . ± . . ± . × − V kick ± − ± − . ± . − V pred
62 km s −
34 km s − . − of 10) in the noise in the imaginary part when comparedto the real part. Figure 8 shows the trajectories for thefour highest resolutions. Note that good agreement be-tween the curves when using the new smaller CFL factor.For reference, we report the final remnant propertiesfor the three configurations in Table II.The above gauge conditions (3) have a drawback nearmerger. The conformal function W goes to zero at thetwo punctures and there is a local maximum betweenthe two punctures. Due to this local maximum, η → η get progressivelycloser to the smaller BH. This, in turn, causes the coor-dinate radius of the horizon to shrink near merger. Asseen in Fig. 11, the horizon radius for the new q = 1 / r H /m H ) ofthe smaller BH is less than that of the larger BH. So therequired resolution (gridpsacing) for each hole does not
500 700 900 1100t/M10 −2 −1 pha s e_d i ff e r en c e s | F _1 − F _2|| F _2 − F _3|| F _2 − F _3| * 3.96993 (4th)| F _2 − F _3| * 2.45799(2nd)| F _2 − F _3| * 1.94462 (1st)500 700 900 1100t/M10 −8 −7 −6 −5 −4 a m p li t ude d i ff e r en c e s || y | −| y | ||| y | −| y | ||| y | −| y | | * 3.96993 (4th)|| y | −| y | | * 45799 (2nd)|| y | −| y | | * 1.94462 (1st) FIG. 4: Convergence of the ( (cid:96) = 2 , m = 2) mode of ψ forthe q = 1 /
10 waveform. The waveform starts out fourth-order accurate (due to our fourth-order waveform extractionalgorithm), but then reduces to second order during the lateinspiral, and dropping to first-order during the ringdown. Thevertical line indicates the time when the frequency is ω =0 . /M . The central resolutions for the three cases were h = M/ h = M/ .
9, and h = M/ .
48, respectively, andeach configuration used 11 levels of refinement. scale exactly with the BH’s mass. This lower effectiverelative size of the smaller BH may be partially respon-sible for the reduced conservation of the horizon massof the smaller BH when compared to the larger on (seeFig. 12).
IV. PERTURBATIVE TECHNIQUES
In Ref. [21], we extended the Regge-Wheeler-Zerilli(RWZ) formalism [19, 20] of black hole perturbation the-ory to include, perturbatively, a term linear in the spinof the larger black hole. This Spin-Regge-Wheeler-Zerilli(SRWZ) formalism has second-order perturbations withlinear dependence on the spin in the wave equations for
100 300 500 700 900t/M−2e−05−1e−0501e−052e−05 Re[ y ] e Re[ y ]
950 1050−0.0002−0.000100.00010.0002 F F (extrap)700 800 900 1000 11001e−021e−011e+001e+011e+02 | F − F (extrap)| FIG. 5: (TOP) A comparison of the highest resolution runwith the Richardson extrapolated waveform. At the plunge,the phase error in the waveform is sufficiently large that theextrapolation give erroneous results (note the unphysical os-cillation near t = 1050 M ). The extrapolation breaks down at ω = 0 . /M . (BOTTOM) The phase difference between theRichardson extrapolated waveform and the highest resolution. Schwarzschild perturbations.In the SRWZ formalism, the wave functions for theeven and odd parity perturbations are expressed as acombined wave function of the first and second pertur-bative orders,Ψ (cid:96)m ( t, r ) = Ψ (1) (cid:96)m ( t, r ) + Ψ (2) (cid:96)m ( t, r ) , Ψ (o) (cid:96)m = Ψ (o , (cid:96)m + 2 (cid:90) dt Ψ (o , Z , (cid:96)m , (5)where the superscript 1 and 2 denote the perturbativeorder, Ψ (cid:96)m denotes the even parity Zerilli function, andΨ (o , (cid:96)m and Ψ (o , Z , (cid:96)m denote the Regge-Wheeler (Cunning-ham et al.) and the Zerilli functions for the odd parityperturbations (see Ref. [44] for the relation between thewave functions). These functions satisfy the following ho r i z on m a ss h /1.32h /1.44h /1.58h /1.74 FIG. 6: The horizon mass for the q = 1 /
100 configurations.Jumps in the horizon mass of order 1 / are apparent be-tween resolutions. this has a small but noticeable dynamiceffect. pun c t u r e s epa r a t i on h /1.32h /1.44h /1.58h /1.74116.3 116.4 116.5 116.6 116.70.70.80.91 FIG. 7: The puncture separation as a function of time forthe q = 1 /
100 configurations. The second lowest resolutionrun appears to merge later than expected, consistent with itsreduced mass (apparent in Fig. 6) −5.5 −3.5 −1.5 0.5 2.5 4.5(x −x )/M−5.5−3.5−1.50.52.54.5 ( y − y ) / M h /1.32h /1.44h /1.58h /1.74 FIG. 8: The puncture trajectories for the q = 1 /
100 configu-ration.
50 100 150 200 250t/M−4e−05−2e−0502e−054e−05 I m [ y ( l = , m = ) ]
50 100 150 200 25001e−072e−073e−074e−075e−07 | y (h /1.44)− y (h /1.58)|| y (h /1.58)− y (h /1.74)| FIG. 9: The imaginary part (which is much less noisy than thereal part) of the ( (cid:96) = 2 , m = 2) mode of ψ for the q = 1 /
50 100 150 200t/M−40−30−20−10010 y ( l = , m = ) P ha s e F [h /1.3] F [h /1.44] F [h /1.58] F [h /1.73]205 206 207 208−27.2−27−26.8−26.6−26.4−26.2 FIG. 10: The phase of the ( (cid:96) = 2 , m = 2) mode of ψ for the q = 1 /
100 configurations. The differences in phases betweenresolutions are dominated by grid noise. ho r i z on c oo r d i na t e r ad i u s
100 r H (small BH)r H (big BH) FIG. 11: The horizon radius of both BHs for the new q =1 /
100 simulations. The horizon radius of the larger BH barelychanges, while the radius of the smaller BH decreases sharplyat merger. ho r i z on m a ss
100 m H (small BH)M H (big BH) FIG. 12: The horizon mass of both BHs for the new q = 1 / three equations. − ∂ ∂t Ψ (cid:96)m ( t, r ) + ∂ ∂r ∗ Ψ (cid:96)m ( t, r ) − V (even) (cid:96) ( r )Ψ (cid:96)m ( t, r ) + i m α ˆ P (even) (cid:96) Ψ (cid:96)m ( t, r )= S (even) (cid:96)m ( t, r ; r p ( t ) , φ p ( t )) , − ∂ ∂t Ψ (o , (cid:96)m ( t, r ) + ∂ ∂r ∗ Ψ (o , (cid:96)m ( t, r ) − V (odd) (cid:96) ( r )Ψ (o , (cid:96)m ( t, r )= S (odd , (cid:96)m ( t, r ; r p ( t ) , φ p ( t )) , − ∂ ∂t Ψ (o , Z , (cid:96)m ( t, r ) + ∂ ∂r ∗ Ψ (o , Z , (cid:96)m ( t, r ) − V (odd) (cid:96) ( r )Ψ (o , Z , (cid:96)m ( t, r )= S (odd , Z , (cid:96)m ( t, r ; r p ( t ) , φ p ( t )) , (6)where r ∗ = r + 2 M ln[ r/ (2 M ) −
1] is a characteris-tic coordinate, α is the non-dimensional spin parame-ter, V (even / odd) (cid:96) denotes the Zerilli and Regge-Wheelerpotential, respectively, S (cid:96)m is the source term, and( r p ( t ) , φ p ( t )) are the particle separation and orbital phaseas a function of time. The differential operator ˆ P (even) (cid:96) arises from the coupling between the BH’s spin and thefirst-order wave functions. We note that S (even) (cid:96)m and S (odd , Z , (cid:96)m include local and extended source terms arisingfrom the second perturbative order associated with prod-ucts of terms linear in α and the first-order perturbationfunctions. See Appendix A in Ref. [21] for more details.In order to evolve the perturbative equations (6), weneed the particle trajectory ( r p ( t ) , φ p ( t )) as a functionof time. Here we use a post-Newtonian inspired modelfor the BH trajectories, which will approximate the nu-merical trajectories from intermediate-mass-ratio BHBinspirals. In order to model the numerical trajectoriesfor the late inspiral phase of a BHB merger using a rel-atively straightforward fitting function, we start with anadiabatic quasicircular evolution, i.e., ignoring the eccen-tricity of the orbit and the component of radial velocitynot due to the radiation reaction, and then parametrizedeviations from a 3.5PN-TaylorT4 adiabatic inspiral. Wemodel the time derivative of the orbital frequency d Ω /dt and the orbital separation R as a function of the orbitalfrequency Ω based on an extension of standard post-Newtonian calculations, described below. A. The orbital frequency Ω evolution of the NRtrajectories First, we focus on the d Ω /dt evolution with respect tothe orbital frequency Ω. In post-Newtonian (PN) cal-culations for the quasicircular case, the evolution for Ωis obtained from the energy loss dE/dt and the relationbetween energy and frequency E (Ω) given by d Ω dt = dEdt (cid:18) dEd Ω (cid:19) − , (7) where dE/dt and dE/d Ω are obtained by appropriate PNexpansions. In the TaylorT1 waveform, we simply sub-stitute these quantities into Eq. (7). By expanding theTaylorT1 in the PN series again, we obtain the TaylorT4waveform, which has been shown to give good agreementwith numerical simulations [45]. Therefore, we develop afitting function based on the TaylorT4 evolution. TheseTaylor series waveforms are summarized in Ref. [46].We note here that for Schwarzschild geodesics, d Ω /dt diverges at the innermost stable circular orbit (ISCO), R Sch = 6 M in Schwarzschild coordinates. This is be-cause dE/d Ω becomes zero at the ISCO frequency M Ω =(1 / / ∼ . d Ω /dt .When fitting of the trajectories in the NR coordinatesfor various mass ratios, we use the following modifiedTaylorT4 evolution. d Ω dt = 965 Ω / M / η (cid:16) B (Ω / Ω ) β/ (cid:17) − (cid:32) (cid:18) − − η (cid:19) ( M Ω) / + 4 π M Ω+ (cid:18) η + 5918 η (cid:19) ( M Ω) / + (cid:18) − π − η π (cid:19) ( M Ω) / + (cid:18) π − γ − M Ω) − η + 45148 η π + 541896 η − η (cid:19) ( M Ω) + (cid:18) − π + 3586756048 η π + 914951512 η π (cid:19) ( M Ω) / + A (Ω / Ω ) α/ (cid:33) , (8)where d Ω /dt and Ω are obtained from the NR trajecto-ries, and M and η denote the total mass and symmetricmass ratio of the binary system, M = m + m ,η = m m M , (9)respectively. Here, A , α , B and β are the fitting param-eters, and the power in the parameters should be α > β > / is the frequency at R Sch = 3 M forcircular orbits, i.e. M Ω = (1 / / ∼ . q = 1 / q = 1 /
15 (from Ref. [21]), and q = 1 / q = 1 /
10 case, which is validup to M Ω = 0 . α and β , for this case are consistent with thecurrent PN formula, and we see that the fitting curve isclose to the TaylorT4 evolution for small M Ω.We show the fitting curve for the q = 1 /
15 case, whichis valid up to M Ω = 0 .
17, in Fig. 14. Again, the fittingparameters, α and β , for this case are consistent with thePN formula.For the q = 1 /
100 case, which is shown in Fig. 15, wefirst found by a direct fit of all parameters that α <
TABLE III: The summary of the fitting parameters for the“PN+” (a modification of the TaylorT4 evolution) in Eq. (8).The first line for the q = 1 /
100 case ( ∗ ) gives the fittingparameters obtained by allowing all 4 constants to vary, thesecond line ( † ) for the q = 1 /
100 case gives the parametersfrom a fitting of A and B only by assuming α = 8 and β = 15,and the third line ( ‡ ) gives the parameters extrapolated fromthe q = 1 /
10 and q = 1 /
15 results using Eq. (11).Mass ratio
A α B βq = 1 /
10 17 . . . . q = 1 /
15 26 . . . . q = 1 / ∗ . . . . q = 1 / † .
888 8 16 . q = 1 / ‡ .
985 9 . . . q = 1 /
10 case using Eq. (8)and orbital frequencies up to M Ω = 0 . R Sch = 2 M , and M Ω = 0 .
175 is the NR orbital frequencyaround R Sch = 3 M , assuming the NR coordinate system canbe approximated by a “trumpet” stationary 1 + log slice ofthe Schwarzschild spacetime. The inset shows the zoom-in ofthe initial part (around M Ω = 0 . Ω M d Ω / d t NR (q=1/10)FittingPN 0.0001 extend them to lower frequencies, we would expect thatthe slope increases, leading to a larger α . Below, we trytwo alternate fitting methods in order to try to improvethe fit. This first fit gives the best agreement with thenumerical trajectory so this is the one we use in Sec. Vto calculate the gravitational waveforms. However, thesecond fit reproduces the expected behavior at large sep-arations, so we expect that it gives a better waveform inthat regime.Note that, as seen in Fig. 15, the end point of theNR trajectory, which corresponds to R Sch = 2 M , has d Ω /dt <
0. This is a feature of Schwarzschild geodesics inSchwarzschild coordinates (e.g., a quasi-circular geodesicslightly inside the ISCO has d Ω /dt = − . q = 1 /
100 case the trajectories aremuch closer to geodesic motion than in the other cases.Due to the fact that the q = 1 /
100 simulations werevery short, obtaining accurate fitting parameters for the
FIG. 14: The fitting curve for the q = 1 /
15 case using Eq. (8)up to orbital frequencies of M Ω = 0 .
17. The solid (black),dashed (red) and dotted (blue) curves show the NR, fittingand TaylorT4 PN evolutions, respectively. In the NR evo-lution, the end point of the solid (black) curve correspondsto R Sch = 2 M , and M Ω = 0 .
17 is the NR orbital frequencyaround R Sch = 3 M , assuming the NR coordinate system canbe approximated by a “trumpet” stationary 1 + log slice ofthe Schwarzschild spacetime. Ω M d Ω / d t NR (q=1/15)FittingPN
FIG. 15: The fitting curve for the q = 1 /
100 case usingEq. (8) with orbital frequencies up to M Ω = 0 .
15. The solid(black) and dotted (blue) curves show the NR and TaylorT4PN evolutions, respectively. The Fitting ( ∗ , dashed (red)),(1, dot-dashed (green)) and (2, dot-dot-dashed (magenta))curves show the fitting curve with the ( ∗ ), ( † ) and ( ‡ ) cases inTable III. In the NR evolution, M Ω = 0 .
15 is the NR orbitalfrequency around R Sch = 3 M , assuming the NR coordinatesystem can be approximated by a “trumpet” stationary 1+logslice of the Schwarzschild spacetime. Ω M d Ω / d t NR (q=1/100)Fitting (*)Fitting (1)Fitting (2)PN trajectory is error prone. We therefore consider severaldifferent techniques. If we assume a simple fitting withrespect to the mass ratio,(
A, α, B, β ) = c η c , (10)for the fitting parameters in the case of the q = 1 /
10 and0 q = 1 /
15 cases, we obtain A = 0 . η − . ,α = 5 . η − . ,B = 5 . η − . ,β = 6 . η − . . (11)Then, inserting the symmetric mass ratio η = 10 /
121 for q = 1 /
100 gives fitting parameters for q = 1 / q = 1 / α = 8 and β = 15 and fit the remaining parameters using the NRtrajectory, we obtain the “Fitting (1)” parameters (sec-ond line labeled q = 1 /
100 in Table III). This fit smoothlyinterpolates between the PN trajectory at small frequen-cies and the start of the numerical trajectory at higherfrequencies. Since both the standard fitting procedure, and “Fitting(1)” give reasonable fits, but different pre-dictions, in order to accurately determine the constant inEqs. (11), we need both more accurate and longer evolu-tions NR simulations of q = 1 / B. The orbital radius R versus orbital frequency Ω in the NR coordinates In order to complete our description of the trajectorieswe need a fitting function that will allow us to expressthe separation R as a function of Ω. We use a func-tion R (Ω) inspired by the ADM-TT PN approach [49],because the isotropic coordinates (used in the ADM-TTPN approach) and the radial isotropic “trumpet” station-ary 1 + log slices of the Schwarzschild spacetime are verysimilar [50, 51].The fitting function based on the ADM-TT PN ap-proach is given by R = M ( M Ω) / (cid:18) (cid:18) − η (cid:19) ( M Ω) / + (cid:18) −
14 + 98 η + 19 η (cid:19) ( M Ω) / + (cid:18) − − η + 167192 η π − η + 281 η (cid:19) ( M Ω) (cid:19) / (1 + a (Ω / Ω ) a ) + C , (12)
TABLE IV: The fitting parameters for the “PN+” (a modi-fication of the ADM-TT PN calculation) in Eq. (12) for therelation between R and Ω in the NR coordinates.Mass ratio C a a q = 1 /
10 0 . . . q = 1 /
15 0 . . . q = 1 /
100 0 . . . where R and Ω are in the NR coordinates, and a , a and C are fitting parameters, and M Ω = (1 / / . Note that a should be greater than 2 to be consistent with the PNcalculation.Initially, we expected that C , which is the differencebetween the NR radial coordinate and “trumpet” coordi-nate, would be very small because this term is inconsis-tent with the ADM-TT PN formula. However, we foundthat this term is non-trivial for all three mass-ratio cases.The results of the fits are shown in Figs. 16-18 for q = 1 / q = 1 /
15, and q = 1 / a for all cases are con-sistent with the current PN calculation ( a > C ∼ . C contribution, there is a finite difference between the or- bital radius evaluated in the ADM-TT PN approach andthat by the fitting function in Eq. (12) even for smaller M Ω.As a result of the above fitting analysis, it is necessaryto consider a radial transformation to remove the offset C between the NR coordinates ( R NR = R for Eq. (12))and the “trumpet” coordinates, R NR → R Log = R NR − C . (13)After removing this offset, we then transform to the stan-dard Schwarzschild coordinates, where the explicit timeand radial coordinate transformations,( T Log , R
Log ) → ( T Sch , R
Sch ) , (14)are given in Ref. [52]. Here, we assume T NR = T Log . Wenote that after applying the C transformation discussedabove, we find that the frequency when R Sch = 3 M isgiven by M Ω = 0 . .
158 and 0 .
152 for the q = 1 / q = 1 /
15 and q = 1 /
100 cases, respectively.Note that this procedure of fitting NR trajectories di-rectly in the NR coordinates and then transforming toSchwarzschild coordinates, appears simpler than the al-ternative of evolving the perturbations directly in the“trumpet” coordinate system. We nevertheless studythis alternative in the Appendix A, where we provide thesourceless part of the perturbative equation and study itscharacteristic speeds in “trumpet” coordinates.1
FIG. 16: The fitting of the orbital radius R vs. orbital fre-quency Ωin the NR coordinates for the q = 1 /
10 case. Thesolid (black), dashed (red) and dotted (blue) curves show theNR, fitting and ADM-TT PN evolutions, respectively. Thefitting by using Eq. (12) is valid up to M Ω = 0 . R Sch = 2 M , and M Ω = 0 .
175 is the NR orbital frequencyaround R Sch = 3 M , assuming the NR coordinate system canbe approximated by a “trumpet” stationary 1 + log slice ofthe Schwarzschild spacetime. Ω R / M NR (q=1/10)FittingPN
FIG. 17: The fitting of the orbital radius R vs. orbital fre-quency Ω in the NR coordinates for the q = 1 /
15 case. Thesolid (black), dashed (red) and dotted (blue) curves show theNR, fitting and ADM-TT PN evolutions, respectively. Thefitting by using Eq. (12) is valid up to M Ω = 0 .
17. In the NRevolution, the end point of the solid (black) curve correspondsto R Sch = 2 M , and M Ω = 0 .
17 is the NR orbital frequencyaround R Sch = 3 M , assuming the NR coordinate system canbe approximated by a “trumpet” stationary 1 + log slice ofthe Schwarzschild spacetime. Ω R / M NR (q=1/15)FittingPN
V. RESULTS
We obtain an approximation to Ω( T NR ) by integrat-ing Eq. (8). We then obtain R ( T NR ) using Eq. (12) andthe approximate Ω( T NR ). The initial orbital phase isnot fixed, and we set it here such that the matchingbetween the perturbative and numerical waveforms aremaximized. FIG. 18: The fitting of the orbital radius R vs. orbital fre-quency Ω in the NR coordinates for the q = 1 /
100 case. Thesolid (black), dashed (red) and dotted (blue) curves show theNR, fitting and ADM-TT PN evolutions, respectively. Thefitting by using Eq. (12) is valid up to M Ω = 0 .
15. In the NRevolution, the end point of the solid (black) curve correspondsto R Sch = 2 M , and M Ω = 0 .
15 is the NR orbital frequencyaround R Sch = 3 M , assuming the NR coordinate system canbe approximated by a “trumpet” stationary 1 + log slice ofthe Schwarzschild spacetime. Ω R / M NR (q=1/100)FittingPN
Note that part of the difference between the NR wave-form and the perturbative waveforms is due to eccentric-ity in the numerical trajectories. This affects the timestep of the orbital evolution through dt = (cid:18) d Ω dt (cid:19) − d Ω . (15)This may help explain why in the q = 1 /
10 case, weobserve a slower time evolution in the fitting trajectorythan the NR evolution.In order to generate the final plunge trajectory andregularize the behavior of the source near the horizon inthe SRWZ formalism [21], we attach the model for thetrajectories to a plunging geodesic in a continuous way.To do this, we choose a matching radius R M ∼ M inSchwarzschild coordinates (the light ring) and then de-termine the energy and angular momentum of the tra-jectory at that point. We then complete the trajectoryfrom R ∼ R M to R ∼ M using the geodesic with thesame energy and angular momentum.When we set the matching radius to R M = 3 M , wefind that the energy E and angular momentum L in theSchwarzschild coordinates are given by E = 13 u t ,L = 9 M u φ , (16)where u µ = dx µ /dτ is the four velocity, and we focus onequatorial orbits. To evaluate u t in the above equation,we use g µν u µ u ν = − , (17)2where g µν denotes the Schwarzschild metric, and at R M = 3 M this becomes u t = (cid:20) −
3( ˙ R Sch ( T Sch )) − M ( ˙Φ Sch ( T Sch )) (cid:21) − / . (18)Here, ˙ R Sch = u r /u t = dR Sch /dT
Sch and ˙Φ
Sch = u φ /u t = d Φ Sch /dT
Sch are the velocity components obtained fromthe numerical data after the coordinate transformationfrom the “trumpet” coordinates to the Schwarzschild co-ordinates. We also use Eq. (17) to calculate u t in thesource term of the SRWZ formalism.In the SRWZ formalism, the gravitational waveformsat a sufficiently distant location R Obs is given by R Obs M ( h + − i h × ) = (cid:88) (cid:96)m (cid:112) ( (cid:96) − (cid:96) ( (cid:96) + 1)( (cid:96) + 2)2 M × (cid:16) Ψ (even) (cid:96)m − i Ψ (odd) (cid:96)m (cid:17) − Y (cid:96)m . (19)We use the normalized waveform, i.e., the expression inthe right hand side of the above equation to calculate thewave amplitude.On the other hand, to convert the NR ψ , which is ex-trapolated to infinity via Eq. (2), to the waveform h , weuse the “pyGWAnalysis” code in EinsteinToolkit [27](see Ref. [53] for details). Then, we calculate the matchbetween the NR and SRWZ waveforms in the advancedLIGO (Zero Det, High Power) noise curve [54]. In Ta-ble V, we summarize the match for the three mass ratios.In the following subsections, we treat each mass ratio caseseparately.
TABLE V: The match between the NR and SRWZ waveformsfor various total masses using the advanced LIGO noise curve.For the smaller masses, the entire waveform is in the advancedLIGO band, while for the larger masses, only the final ring-down part of the waveform is in band. For these larger masses,we integrate over frequencies f ≥ (cid:96) = 2 , m = 2) mode to calculate the match.Total mass ( M (cid:12) ) q = 1 / q = 1 / q = 1 / A. The q = 1 / case For the q = 1 /
100 case, the fitting formulas in Eqs. (8)and (12) are valid for orbital frequencies up to M Ω = 0 . M Ω = 0 .
15 by attaching these trajectories to aplunging geodesic at R M (as explained above).Ideally, we would like to use as much of the fitting tra-jectory as possible. However, when attempting to usethe fitting trajectory to calculate the 4-velocity at smallradii, we can not evaluate u t for R Sch ≤ . M fromEq. (17) because the approximation that the backgroundmetric is Schwarzschild in “trumpet” coordinates breaksdown and the vector (1 , ˙ R Sch , , ˙Φ Sch ) becomes spacelikein the background metric (but not in the numerical met-ric). This breakdown can be thought to arise from theinteraction of the singular part of the conformal factordue to the smaller BH (when it is close to horizon of thelarger one) with the singular part of the conformal fac-tor due to the larger BH. This, in turn, changes ˜Γ i , andtherefore the gauge, when the two BHs approach eachother. Also, since we see some delay in the phase evolu-tion for a small matching radius near the horizon, we set R M ∼ M for the start of the geodesic approximation.At this matching radius, E and L for the Schwarzschildgeodesic are evaluated from the orbital separation andthree velocity using Eq. (16). We find E = 0 . ,L/M = 3 . . (20)In the SRWZ formalism, we set the non-dimensional spinparameter α = 0 .
033 (see Table II). Although this valueis obtained in the NR simulation, we can also obtain itfrom an empirical formula [55].We next calculate the match in the advanced LIGOnoise curve. For the integration in the frequency domain,we consider gravitational wave frequencies M Ω ≥ . M = 484 M (cid:12) , where we have the lowerfrequency cut off (corresponding to the start of the nu-merical waveform) at f low = 10 . M = 0 . , (21)for the ( (cid:96) = 2 , m = 2) mode. The amplitude and phaseevolutions for the NR and SRWZ waveforms are shownin Fig. 19. To obtain this figure, we translated the wave-forms so that the maximum in the amplitude occurs at t = 0 and used the freedom in the phase to set it to zeroat t = − M .The ( (cid:96) = 2 , m = 1) mode of h is the leading contri-bution from the odd parity in the RWZ calculation. Theamplitude and phase evolutions for this mode are shownin Fig. 20. After adding a time translation so that themaximum in the amplitude occurs at t = 0 and a phasetranslation so that the phase φ = 0 is zero at t = − M ,we still see a small difference between the phases. For atotal mass of M = 242 M (cid:12) (this mass is chosen such thatthe lowest frequency is just in the advanced LIGO band f low = 10 . M = 0 . . (22)3 FIG. 19: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 2 , m = 2) mode of h for the q = 1 /
100 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -50 0 50t/M00.0050.010.0150.02 NRSRWZ -50 0 50t/M-30-20-100 NRSRWZ For the ( (cid:96) = 3 , m = 3) mode of h , the amplitude andphase evolutions are shown in Fig. 21. For M = 726 M (cid:12) (mass chosen so that f low = 10 . M = 0 . . (23) B. The q = 1 / case For the q = 1 /
10 case, the fitting formulae is valid upto larger orbital frequencies than for the q = 1 /
100 case.Nevertheless, u t can not be evaluated for R Sch ≤ . M from Eq. (17). Therefore, we need to start the geodesicapproximation at radii larger than this, and again we setthe matching radius R M ∼ M . We find that the energyand angular momentum for the Schwarzschild geodesicare given by E = 1 . ,L/M = 3 . , (24) FIG. 20: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 2 , m = 1) mode of h for the q = 1 /
100 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -50 0 50t/M00.0010.0020.0030.0040.0050.006 NRSRWZ -50 0 50t/M-30-25-20-15-10-50 NRSRWZ at the matching radius. Not that the energy is E > q = 1 / α = 0 .
26, consistent with the measuredspin from the numerical simulation (see Table II).The amplitude and phase evolutions of the ( (cid:96) = 2 , m =2) mode of h are shown in Fig. 22. While the phases ofthe NR and SRWZ waveforms are in good agreement,there are relatively large differences in the amplitudesboth for the inspiral and merger phases. We will discussthe amplitude differences in more detail in Sec. V D.We consider the match calculation for M Ω ≥ . M = 242 M (cid:12) , i.e., the lower frequency cutoff at f low = 10 . M = 0 . , (25)for the ( (cid:96) = 2 , m = 2) mode in the advanced LIGO noisecurve.4 FIG. 21: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 3 , m = 3) mode of h for the q = 1 /
100 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -50 0 50t/M00.0010.0020.0030.0040.0050.006 NRSRWZ -50 0 50t/M-60-50-40-30-20-100 NRSRWZ The amplitude and phase evolutions for the ( (cid:96) =2 , m = 1) mode of h are shown in Fig. 23., For M =121 M (cid:12) ( f low = 10 . M = 0 . . (26)The amplitude and phase evolutions for the ( (cid:96) =3 , m = 3) mode of h are shown in Fig. 24. For M =363 M (cid:12) ( f low = 10 . M = 0 . . (27)We have a similar mismatch in the ( (cid:96) = 2 , m = 1) and( (cid:96) = 3 , m = 3) modes.We note that although the ( (cid:96) = 2 , m = 1) and ( (cid:96) =3 , m = 3) modes have a slightly large mismatch (1 −M ), these modes are the sub-dominant contributions tothe gravitational wave. The maximum amplitudes are ∼ .
12, 0 .
035 and 0 .
04 for the ( (cid:96) = 2 , m = 2), ( (cid:96) = 2 , m =1) and ( (cid:96) = 3 , m = 3) modes, respectively. Therefore,we find that from a rough estimation the contributionof the mismatch for each mode are almost same in thegravitational wave.
FIG. 22: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 2 , m = 2) mode of h for the q = 1 /
10 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -800 -600 -400 -200 0t/M00.020.040.060.080.10.120.14 NRSRWZ -800 -600 -400 -200 0t/M-100-80-60-40-200 NRSRWZ Also, as discussed in Ref. [21], unlike the case of q =1 / C. The q = 1 / case For the q = 1 /
15 case, the fitting formula is valid upto larger orbital frequencies than for the q = 1 /
100 case.Here too, we can not calculate u t for R Sch ≤ . M fromEq. (17), and again we use the geodesic approximationand a matching radius R M ∼ M . We find that theenergy and angular momentum are given by E = 0 . ,L/M = 3 . . (28)In the SRWZ formalism, we set the non-dimensional spinparameter α = 0 .
19, consistent with the numerical results5
FIG. 23: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 2 , m = 1) mode of h for the q = 1 /
10 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -800 -600 -400 -200 0t/M00.010.020.030.040.05 NRSRWZ -800 -600 -400 -200 0t/M-60-50-40-30-20-100 NRSRWZ (see Table II).The amplitude and phase evolutions of the ( (cid:96) = 2 , m =2) mode of h are shown in Fig. 25. The match for M Ω ≥ .
09 with a total mass of M = 290 M (cid:12) ( f low =10 . M = 0 . , (29)using the advanced LIGO noise curve.The amplitude and phase evolutions for the ( (cid:96) =2 , m = 1) mode of h are shown in Fig. 26. For M =145 M (cid:12) ( f low = 10 . M = 0 . . (30)The amplitude and phase evolutions for the ( (cid:96) =3 , m = 3) mode of h are shown in Fig. 27. For M =435 M (cid:12) ( f low = 10 . M = 0 . . (31)As was seen in the q = 1 /
10 case, we have a somewhatlarge mismatch for the sub-dominant ( (cid:96) = 2 , m = 1) and( (cid:96) = 3 , m = 3) modes.
FIG. 24: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 3 , m = 3) mode of h for the q = 1 /
10 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -800 -600 -400 -200 0t/M00.010.020.030.040.050.06 NRSRWZ -800 -600 -400 -200 0t/M-150-100-500 NRSRWZ D. Amplitude differences
Here, we focus on the amplitude of the ( (cid:96) = 2 , m = 2)mode of h obtained using NR and the SRWZ formalism.In Fig. 28, we show the amplitude difference between theNR and SRWZ waveforms for the q = 1 /
10 case. Thereis a large ( ∼ δA = | A (SRWZ)22 − A (NR)22 | , (32)we find that the amplitude differences for the q = 1 / q = 1 /
15 cases around the maximum amplitude obey δA ( q =1 / ∼ . × δA ( q =1 / ∼ . . × δA ( q =1 / , (33)where 1 .
41 is the ratio of the symmetric mass ratios η .6 FIG. 25: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 2 , m = 2) mode of h for the q = 1 /
15 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -600 -500 -400 -300 -200 -100 0 100t/M00.020.040.060.080.10.12 NRSRWZ -600 -500 -400 -300 -200 -100 0 100t/M-80-60-40-200 NRSRWZ Figure 29 shows the amplitude differences (note that t = 0 denotes the time of the maximum amplitude) for q = 1 /
10 and q = 1 /
15. We see that the amplitudedifference before merger is proportional to the ratio ofthe symmetric mass ratio. This behavior is different fromour previous analysis in Ref. [21], where we directly usedthe NR trajectories in the perturbative calculation of thewaveforms, and found that the amplitude difference inthe inspiral phase scales like η . Since the fitting formulais based on the quasicircular evolution, any eccentricityof the orbit may create amplitude differences during theinspiral phase.Figure 30 shows the amplitude differences both for the q = 1 /
10 and q = 1 /
100 waveforms. Near the maximumamplitude, we find that the amplitude differences scalelike δA ( q =1 / ∼ . × δA ( q =1 / ∼ . . × δA ( q =1 / , (34)where 8 .
43 is the ratio of the symmetric mass ratio.Note that this rescaling does not capture the behavior
FIG. 26: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 2 , m = 1) mode of h for the q = 1 /
15 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -600 -500 -400 -300 -200 -100 0 100t/M00.010.020.030.04 NRSRWZ -600 -500 -400 -300 -200 -100 0 100t/M-50-40-30-20-100 NRSRWZ of the amplitude differences during the inspiral phase.Although it is difficult to obtain any exact relation be-tween the amplitude differences during the inspiral (pos-sibly due to effects of eccentricity), we do observe non-linear effects in the amplitude difference between the NRand SRWZ waveforms around the maximum amplitude,which we infer from the fact that the amplitude differencescales nonlinearly with the mass ratio.Figure 31 shows the amplitude differences both for the q = 1 /
15 and q = 1 /
100 cases. Near the maximum am-plitude, we find that the amplitude differences scale like δA ( q =1 / ∼ . × δA ( q =1 / ∼ . . × δA ( q =1 / , (35)where 5 .
98 is the ratio of the symmetric mass ratios.Again, we see the similar behavior to q = 1 /
10 and q = 1 /
100 comparison.We note that the C transformation in Eq. (13) alsoaffects the amplitude difference. As seen in Fig. 32, ifwe directly use the orbital radius R NR to calculate thewaveforms, we obtain an even larger amplitude than that7 FIG. 27: The amplitude (TOP) and phase (BOTTOM) ofthe ( (cid:96) = 3 , m = 3) mode of h for the q = 1 /
15 case. Thewaveforms were translated in time such that the maximum inthe amplitude occurs at t = 0 and the phases were adjustedby a constant offset such that φ = 0 at t = − M . Thesolid (black) and dashed (red) curves show the NR and SRWZwaveforms, respectively. -600 -500 -400 -300 -200 -100 0 100t/M00.010.020.030.04 NRSRWZ -600 -500 -400 -300 -200 -100 0 100t/M-140-120-100-80-60-40-200 NRSRWZ FIG. 28: The amplitude difference of the ( (cid:96) = 2 , m = 2)mode of h for the q = 1 /
10 case. Here, the difference isnormalized by the NR amplitude. We see ∼
10% differencein the inspiral and merger phases. -800 -600 -400 -200 0t/M00.050.1 obtained using the corrected orbital radius R Log . Thisbehavior is expected based on the fact that while theorbital radius is larger, the orbital frequency remains the
FIG. 29: The amplitude difference δA of the ( (cid:96) = 2 , m = 2)mode of h for the q = 1 /
10 and 1 /
15 cases. -500 -400 -300 -200 -100 0 100t/M00.0050.010.015 (q=1/10)(q=1/15)(q=1/15) * 1.41(q=1/15) * 2.3
FIG. 30: The amplitude difference δA of the ( (cid:96) = 2 , m = 2)mode of h for the q = 1 /
10 and 1 /
100 cases. -80 -60 -40 -20 0 20 40 60t/M00.0050.010.015 (q=1/10)(q=1/100)(q=1/100) * 8.43(q=1/100) * 21.5
FIG. 31: The amplitude difference δA of the ( (cid:96) = 2 , m = 2)mode of h for the q = 1 /
15 and 1 /
100 cases. -80 -60 -40 -20 0 20 40 60t/M00.0010.0020.0030.0040.0050.0060.007 (q=1/15)(q=1/100)(q=1/100) * 5.98(q=1/100) * 9.5 same. Although the orbital radius R NR gives a similaramplitude to the NR waveform in the initial inspiral part,we have a much larger amplitude in the merger phase.8 FIG. 32: The amplitude of the ( (cid:96) = 2 , m = 2) mode of h forthe q = 1 /
10 case. We have set the maximum amplitudes at t = 0 by the time translation. The solid (black), dashed (red)and dotted (blue) curves show the waveforms from the NR,SRWZ and SRWZ without the C transformation in Eq. (12),respectively. -800 -600 -400 -200 0t/M00.050.10.150.2 NRSRWZSRWZ (without C) VI. DISCUSSION
In the first paper of this series [21] we established thatperturbation theory of black holes in the extreme-mass-ratio expansion q = m /m (cid:28)
1, as in Eq. (6), canfaithfully represent an approximation to the full numer-ical simulations, if we provide the full numerical relativetrajectory of the holes as a function of time. In thissecond paper we proceed to model these trajectories us-ing full numerical simulations and a parametrized exten-sion of 3.5 PN quasi-adiabatic evolutions, Eqs. (8) and(12). This form of the expansion allows us to incorpo-rate the PN trajectory at large distances and smoothlyfit to the full numerical ones up to separations as smallas R Sch = 3 M , at which point we can describe the sub-sequent trajectory by a Schwarzschild plunging geodesic.In this approximation, we set the values of the energy E and angular momentum L of the geodesic to those of thepreviously fitted trajectory evaluated at R Sch = 3 M .Once we verified the procedure for each of the testruns, i.e. for q = 1 / , /
15 and 1 / q , asgiven in Eqs. (11). The other piece of information neededby the perturbative approach, in addition to the relativetrajectory, are the mass and spin of the final black holeformed after merger to provide the background metric.The spin and mass parameters for the SRWZ equation(6) can be estimated from empirical formulae [55, 56].Note that this formula correctly predicted the remnantblack hole’s mass and spin [16] for the three cases of massratio studied here.The current study can be considered a successful proofof principle, and needs to be supplemented by a thor-ough set of runs in the intermediate mass ratio regime0 . < ∼ q < ∼ . R > ∼ M and go through the merger of the holes until the forma-tion of a common horizon. The fitting formulae Eqs. (8)and (12) with Eqs. (11) should provide a good approx-imation within the fitting interval, but extrapolation toboth the comparable masses or extreme mass ratio couldquickly lose accuracy. The fitted tracks from our workcan be compared to alternative approaches to the prob-lem such as the EOB [57] and EOBNR [58, 59] andthe final waveforms to direct phenomenological ones asgiven in Refs. [60, 61]. Our approach provides a comple-mentary model for the intermediate mass- ratio regimeto the above approaches, which mostly apply to thecomparable-mass regime.The extension of our approach to initially highly-spinning black holes can be done essentially through thesame lines depicted here. Numerical simulations of highlyspinning holes with small mass ratio are feasible [62], andthe trajectories would already include the effect of thespins of both holes. The perturbative waveform calcula-tion would use the Teukolsky formalism [63] and its sec-ond order extension [31]. Fitting formulae based on theextension of spinning post-Newtonian trajectories withfree parameter dependence on the mass ratio and effec-tive spins should be used together with the final matchingto plunging geodesics on an spinning background. Thefinal remnant formulae for spinning binaries have alreadybeen proposed and fitted to available simulations [55, 56]. Acknowledgments
We thank Marcelo Ponce for providing the methodto reduce the horizon mass loss required to make ac-curate simulations in the small mass ratio regime. Wegratefully acknowledge the NSF for financial supportfrom Grants No. PHY-0722315, No. PHY-0653303, No.PHY-0714388, No. PHY-0722703, No. DMS-0820923,No. PHY-0929114, PHY-0969855, PHY-0903782, andNo. CDI-1028087; and NASA for financial support fromNASA Grants No. 07-ATFP07-0158 and No. HST-AR-11763. Computational resources were provided bythe Ranger cluster at TACC (Teragrid allocation TG-PHY060027N) and by NewHorizons at RIT.
Appendix A: Perturbative Equations in 1+log“trumpet” coordinates
In this work we used standard Schwarzschild coordi-nates to describe the background geometry and trajec-tory of the particle. To obtain the latter, we transformedthe particle’s location from the NR coordinates (assum-ing they were represented by 1 + log trumpet coordi-nates of a Schwarzschild spacetime) into the standardSchwarzschild coordinates. It is also useful to see theform of the perturbative equations in these 1 + log trum-pet coordinates. This will help to assess potential regionswhere this approach might have limitations. In partic-9ular, since the trumpet coordinates are obtained in thelate-time, stationary limit, wherever there still are strongdynamics, the approach could be inaccurate.The generalized Regge-Wheeler-Zerilli equations for anarbitrary spherically symmetric background were foundin Ref. [64]. Here, for the sake of simplicity, we useonly the generalized Regge-Wheeler equation given inRef. [65]. Formally, this equation can be written as¨Φ RW = c ˙Φ (cid:48) RW + c Φ (cid:48)(cid:48) RW + c ˙Φ RW + c Φ (cid:48) RW + c Φ RW , (A1)where the dot and dash denote the time and radial deriva-tives, and the coefficients c i ( i = 1 , , , ,
5) are deter-mined by the background metric.For the background metric, we use the notation of [52],which has the form ds = − ( α − β ) dt + 2 βf dtdr + 1 f dr + R ( dθ + sin θdφ ) . (A2)where α , β , f and R are functions of the radial coordinate r only. The coefficients c i are then given by c = 2 f ( r ) β ( r ) ,c = ( α ( r ) − β ( r ) ) f ( r ) ,c = − ( β ( r ) α (cid:48) ( r ) − α ( r ) β (cid:48) ( r )) f ( r ) α ( r ) ,c = f ( r ) α ( r ) ( − α ( r ) β ( r ) β (cid:48) ( r ) + α ( r ) α (cid:48) ( r )+ β ( r ) α (cid:48) ( r )) + f ( r ) f (cid:48) ( r ) ( α ( r ) − β ( r ) ) ,c = − α ( r ) V ( r ) , (A3)where V denotes the Regge-Wheeler potential, V ( r ) = 1 R ( r ) (cid:18) (cid:96) ( (cid:96) + 1) − MR ( r ) (cid:19) . (A4)When using the standard Schwarzschild coordinates, theabove equation becomes the standard Regge-Wheelerequation. We note that for the Zerilli equation, we mayset [64] V ( r ) = 2 R ( r ) ( λR ( r ) + 3 M ) (cid:2) λ ( λ + 1) R ( r ) +3 λ M R ( r ) + 9 λM R ( r ) + 9 M (cid:3) , (A5) where λ = ( (cid:96) − (cid:96) + 2) / α ( r ) = f ( r ) R (cid:48) ( r ) ,α ( r ) − β ( r ) = 1 − MR ( r ) . (A6)From the stationary 1 + log slice condition, we have Ce α = (cid:18) α − MR (cid:19) R , (A7)where the constant C is determined by requiring that α (cid:48) be regular at the “critical” point (where the denominatorvanishes). CM = 1128 (3 + √ exp(3 − √ . (A8)In isotropic coordinates we obtain α ( r ) = r R (cid:48) ( r ) R ( r ) . (A9)Since the above equations are transcendental, we haveto find numerical solutions. However, it is also useful tostudy their asymptotic behavior for large r . At large r we have, R ( r ) = r + M + 14 M r − C er + O (1 /r ) ,α ( r ) = 1 − Mr + 12 M r − M r + (cid:18)
18 + 12 CM e (cid:19) M r + O (1 /r ) , (A10)and the other two functions β and f can be found fromEqs. (A6).Inserting these functions into Eqs. (A3), we find thatthe coefficients c i are given by0 c = 2 √ C er + O (1 /r ) ,c = 1 − Mr + 172 M r − M r + (cid:18) CM e (cid:19) M r + O (1 /r ) ,c = − √ C er + O (1 /r ) ,c = 2 Mr − M r + 392 M r + (cid:18) − CM e (cid:19) M r + O (1 /r ) ,c = (cid:96) ( (cid:96) + 1) r − M (1 + (cid:96) ( (cid:96) + 1)) r + 12 M (84 + 37 (cid:96) ( (cid:96) + 1)) r − M (303 + 79 (cid:96) ( (cid:96) + 1)) r + (cid:18)
378 + 105916 (cid:96) ( (cid:96) + 1) + 94 CM e (cid:96) ( (cid:96) + 1) (cid:19) M r + O (1 /r ) . (A11)If we set C = 0, the above coefficients are the sameas those derived in the standard isotropic coordinates( r Sch = r Iso (1+ M/ (2 r Iso )) ) for the Schwarzschild space-time.An important piece of information in setting up thegrid refinements, particularly for the smaller black holeis the perturbative potential [16]. Using the numericalmethod given in Ref. [52], we can obtain the potentialterm c . The “Regge-Wheeler” potentials for the stan-dard Schwarzschild ( r Sch ), isotropic ( r Iso ), and the trum-pet radial coordinates ( r Log ) (associated to the NR co-ordinates) are shown in Fig. 33, while the radial deriva-tives of these potentials are shown in Fig. 34. In thetrumpet coordinates, the radial coordinate covers up to r Log = 0, and α and R around r Log = 0 have the followingforms [52]. α ( r Log ) ∼ r . ,R ( r Log ) ∼ R + (cid:16) r Log . (cid:17) . , (A12)where R (cid:39) . M . From the above equation, we seethat the sphere corresponding to r Log = 0 has a finitecircumference.We observe that the differences between the “trumpet”and isotropic coordinates (outside the horizon), as deter-mined by differences in the potentials, occurs in a regionbetween the horizon and the maximum of the potential.This lends credence to the postulate that the main differ-ences between the numerical coordinates and the “trum-pet” coordinates lie in this region, as well. Indeed, as thetwo BHs approach each other, the conformal function W near the larger BH must change (since W = 0 at the lo-cation of the smaller BH, and varies smoothly outside thetwo punctures). Thus, as the small BH falls through themaximum of the potential, the coordinate system mustbecome distorted in this important region. Therefore,we may expect that the assumption that the trajectoryin the background spacetime is well described by assum-ing the coordinate radius is the trumpet radial coordinate breaks down. This may explain the need to match tra-jectories to plunging geodesics with slightly larger E and L values than one might expect (i.e. E > c + = f ( r ) ( α ( r ) − β ( r )) ,c − = − f ( r ) ( α ( r ) + β ( r )) . (A13)Figure 35 shows the radial dependence in the trum-pet coordinates. Note that inside the horizon, ( r Log =0 . M ), both speeds remain negative and that the out-going characteristic speed vanishes precisely at the hori-zon. Both modes reach their asymptotic speeds ( ± FIG. 33: The solid (black), dashed (blue), and dotted (red)curves show the standard Regge-Wheeler (TOP) and Zerilli(BOTTOM) potentials in the standard the Schwarzschild co-ordinates, isotropic coordinates, and “trumpet” coordinates(which are the stationary 1 + log slices of the Schwarzschildspacetime in isotropic coordinates), respectively, for the (cid:96) = 2mode. Note that the location of the horizon is different in eachsystem. r Sch < M , r Iso < M/ r Log < . M are insidethe horizon. Note that the potential is finite at the horizon,and is always positive. Sch /M, r
Iso /M, r
Log /M-0.2-0.100.10.20.3 - M c SchwarzschildIsotropicTrumpet0 2 4 6 8 10r
Sch /M, r
Iso /M, r
Log /M-0.2-0.100.10.20.3 - M c SchwarzschildIsotropicTrumpet [1] F. Pretorius, Phys. Rev. Lett. , 121101 (2005), gr-qc/0507014.[2] M. Campanelli, C. O. Lousto, P. Marronetti, andY. Zlochower, Phys. Rev. Lett. , 111101 (2006), gr-qc/0511048.[3] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, andJ. van Meter, Phys. Rev. Lett. , 111102 (2006), gr-qc/0511103.[4] B. Aylott et al., Class. Quant. Grav. , 165008 (2009),0901.4399.[5] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-ritt, Astrophys. J. , L5 (2007), gr-qc/0701164.[6] J. A. Gonz´alez, M. D. Hannam, U. Sperhake, B. Brug- mann, and S. Husa, Phys. Rev. Lett. , 231101 (2007),gr-qc/0702052.[7] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-ritt, Phys. Rev. Lett. , 231102 (2007), gr-qc/0702133.[8] C. O. Lousto and Y. Zlochower (2011), 1108.2009.[9] Y. Sekiguchi and M. Shibata, Astrophys. J. , 6(2011), 1009.5303.[10] C. O. Lousto and Y. Zlochower, Phys. Rev. D77 , 024034(2008), 0711.1165.[11] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.Rev.
D79 , 084012 (2009), 0811.3006.[12] R. Owen, Phys. Rev.
D81 , 124042 (2010), 1004.3768.[13] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys. FIG. 34: The solid (black), dashed (blue), and dotted (red)curves show the radial derivative of the standard Regge-Wheeler (TOP) and Zerilli (BOTTOM) potentials in the stan-dard Schwarzschild coordinates, isotropic coordinates, and“trumpet” coordinates, respectively, for the (cid:96) = 2 mode. r Sch < M , r Iso < M/
2, and r Log < . M are insidethe horizon. The vertical dashed lines denotes the location ofthe horizon in the “trumpet” coordinates. Sch /M, r
Iso /M, r
Log /M-0.500.5 - M d c / d r SchwarzschildIsotropicTrumpet0 2 4 6 8 10r
Sch /M, r
Iso /M, r
Log /M-0.500.5 - M d c / d r SchwarzschildIsotropicTrumpet
Rev.
D74 , 041501(R) (2006), gr-qc/0604012.[14] M. Ponce, C. Lousto, and Y. Zlochower, Class. Quant.Grav. , 145027 (2011), 1008.2761.[15] M. Shibata and H. Yoshino, Phys. Rev. D81 , 104035(2010), 1004.4970.[16] C. O. Lousto and Y. Zlochower, Phys. Rev. Lett. ,041101 (2011), 1009.0292.[17] U. Sperhake, V. Cardoso, C. D. Ott, E. Schnetter, andH. Witek (2011), 1105.5391.[18] G. Lovelace, M. A. Scheel, and B. Szilagyi, Phys. Rev.
D83 , 024010 (2011), 1010.2777.[19] T. Regge and J. A. Wheeler, Phys. Rev. , 1063(1957).[20] F. J. Zerilli, Phys. Rev. D2 , 2141 (1970).[21] C. O. Lousto, H. Nakano, Y. Zlochower, and M. Cam-panelli, Phys. Rev. D82 , 104057 (2010), 1008.4360.[22] S. Brandt and B. Br¨ugmann, Phys. Rev. Lett. , 3606(1997), gr-qc/9703066. FIG. 35: The ingoing and outgoing characteristic speeds.We see that c + = 0 but c − is finite at the horizon (denotedby the vertical dashed line, r Log = 0 . M ), and that bothspeeds are negative inside the horizon. The thin curves showthe characteristic speeds for the perturbation in the isotropiccoordinates. Note that in this latter case, both speeds arezero on the horizon ( r Iso = M/ Log /M (r
Iso /M)-0.500.5 c + c - [23] M. Ansorg, B. Br¨ugmann, and W. Tichy, Phys. Rev. D70 , 064011 (2004), gr-qc/0404056.[24] Y. Zlochower, J. G. Baker, M. Campanelli, and C. O.Lousto, Phys. Rev.
D72 , 024021 (2005), gr-qc/0505055.[25] P. Marronetti, W. Tichy, B. Br¨ugmann, J. Gonzalez, andU. Sperhake, Phys. Rev.
D77 , 064010 (2008), 0709.2160.[26] Cactus Computational Toolkit home page: .[27] Einstein Toolkit home page: http://einsteintoolkit.org .[28] E. Schnetter, S. H. Hawley, and I. Hawke, Class. Quan-tum Grav. , 1465 (2004), gr-qc/0310042.[29] J. Thornburg, Class. Quant. Grav. , 743 (2004), gr-qc/0306056.[30] O. Dreyer, B. Krishnan, D. Shoemaker, and E. Schnetter,Phys. Rev. D67 , 024018 (2003), gr-qc/0206008.[31] M. Campanelli and C. O. Lousto, Phys. Rev.
D59 ,124022 (1999), gr-qc/9811019.[32] C. O. Lousto and Y. Zlochower, Phys. Rev.
D76 ,041502(R) (2007), gr-qc/0703061.[33] M. Babiuc, B. Szilagyi, J. Winicour, and Y. Zlochower(2010), 1011.4223.[34] C. Reisswig, N. T. Bishop, D. Pollney, and B. Szilagyi,Phys. Rev. Lett. , 221101 (2009), 0907.2637.[35] D. Pollney, C. Reisswig, E. Schnetter, N. Dorband, andP. Diener, Phys. Rev.
D83 , 044045 (2011), 0910.3803.[36] B. Szilagyi, L. Lindblom, and M. A. Scheel, Phys. Rev.
D80 , 124010 (2009), 0909.3557.[37] M. Alcubierre, B. Br¨ugmann, P. Diener, M. Koppitz,D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev.
D67 ,084023 (2003), gr-qc/0206072.[38] M. Alcubierre et al. (2004), gr-qc/0411137.[39] D. M¨uller and B. Br¨ugmann, Class. Quant. Grav. ,114008 (2010), 0912.3125. [40] D. M¨uller, J. Grigsby, and B. Br¨ugmann, Phys. Rev. D82 , 064004 (2010), 1003.4681.[41] E. Schnetter, Class. Quant. Grav. , 167001 (2010),1003.0859.[42] D. Alic, L. Rezzolla, I. Hinder, and P. Mosta, Class.Quant. Grav. , 245023 (2010), 1008.2212.[43] J. A. Gonz´alez, U. Sperhake, B. Brugmann, M. Hannam,and S. Husa, Phys. Rev. Lett. , 091101 (2007), gr-qc/0610154.[44] C. O. Lousto, Class. Quant. Grav. , S543 (2005), gr-qc/0503001.[45] M. Boyle et al., Phys. Rev. D76 , 124038 (2007),0710.0158.[46] P. Ajith, M. Boyle, D. A. Brown, S. Fairhurst, M. Han-nam, I. Hinder, S. Husa, B. Krishnan, R. A. Mercer,F. Ohme, et al. (2007), 0709.0093.[47] A. Ori and K. S. Thorne, Phys. Rev.
D62 , 124022 (2000),gr-qc/0003032.[48] M. Kesden, Phys. Rev.
D83 , 104011 (2011), 1101.3749.[49] G. Sch¨afer,
Post-Newtonian methods: Analytic results onthe binary problem (Springer, New York, 2010), chap. 6,0910.2857.[50] B. J. Kelly, W. Tichy, M. Campanelli, and B. F. Whiting,Phys. Rev.
D76 , 024008 (2007), 0704.0628.[51] W. Tichy, B. Br¨ugmann, M. Campanelli, and P. Diener,Phys. Rev.
D67 , 064008 (2003), gr-qc/0207011.[52] B. Br¨ugmann, Gen. Rel. Grav. , 2131 (2009),0904.4418.[53] C. Reisswig and D. Pollney (2010), 1006.1632. [54] Advanced ligo anticipated sensitivity curves , URL https://dcc.ligo.org/cgi-bin/DocDB/ShowDocument?docid=2974 .[55] C. O. Lousto, M. Campanelli, Y. Zlochower, andH. Nakano, Class. Quant. Grav. , 114006 (2010),0904.3541.[56] L. Rezzolla, Class. Quant. Grav. , 094023 (2009),0812.2325.[57] T. Damour and A. Nagar, Phys. Rev. D79 , 081503(2009), 0902.0136.[58] Y. Pan et al., Phys. Rev.
D81 , 084041 (2010), 0912.3466.[59] N. Yunes, A. Buonanno, S. A. Hughes, Y. Pan, E. Ba-rausse, et al., Phys. Rev.
D83 , 044044 (2011), 1009.6013.[60] L. Santamaria, F. Ohme, P. Ajith, B. Bruegmann,N. Dorband, et al., Phys. Rev.
D82 , 064016 (2010),1005.3306.[61] R. Sturani, S. Fischetti, L. Cadonati, G. Guidi, J. Healy,et al. (2010), 1012.5172.[62] C. O. Lousto and Y. Zlochower, Phys. Rev.
D79 , 064018(2009), 0805.0159.[63] S. A. Teukolsky, Astrophys. J. , 635 (1973).[64] O. Sarbach and M. Tiglio, Phys. Rev.
D64 , 084016(2001), gr-qc/0104061.[65] E. Pazos et al., Class. Quant. Grav. , S341 (2007),gr-qc/0612149.[66] A. Zenginoglu, Class. Quant. Grav.27