Ineffectiveness of Padé resummation techniques in post-Newtonian approximations
aa r X i v : . [ g r- q c ] M a y Ineffectiveness of Pad´e resummation techniques in post-Newtonian approximations
Abdul H. Mrou´e, Lawrence E. Kidder, and Saul A. Teukolsky Center for Radiophysics and Space Research, Cornell University, Ithaca, New York, 14853 (Dated: September 15, 2018)We test the resummation techniques used in developing Pad´e and Effective One Body (EOB)waveforms for gravitational wave detection. Convergence tests show that Pad´e approximants of thegravitational wave energy flux do not accelerate the convergence of the standard Taylor approxi-mants even in the test mass limit, and there is no reason why Pad´e transformations should help inestimating parameters better in data analysis. Moreover, adding a pole to the flux seems unnec-essary in the construction of these Pad´e-approximated flux formulas. Pad´e approximants may beuseful in suggesting the form of fitting formulas. We compare a 15-orbit numerical waveform of theCaltech-Cornell group to the suggested Pad´e waveforms of Damour et al. in the equal mass, non-spinning quasi-circular case. The comparison suggests that the Pad´e waveforms do not agree betterwith the numerical waveform than the standard Taylor based waveforms. Based on this result, wedesign a simple EOB model by modifiying the ET EOB model of Buonanno et al., using the Taylorseries of the flux with an unknown parameter at the fourth post-Newtonian order that we fit for.This simple EOB model generates a waveform having a phase difference of only 0.002 radians withthe numerical waveform, much smaller than 0.04 radians the phase uncertainty in the numericaldata itself. An EOB Hamiltonian can make use of a Pad´e transformation in its construction, butthis is the only place Pad´e transformations seem useful.
PACS numbers: 04.25.dg, 04.25.Nx, 04.30.Db
I. INTRODUCTION
Even though general relativity was developed at thebeginning of the twentieth century, no analytical solutionis known for the two-body problem. Until recently, at-tempts to find a numerical solution failed because of thecomplexity of the mathematical equations and the insta-bilities inherent in the analytical formulations being used.In the past few years, breakthroughs in numerical rela-tivity [1, 2, 3, 4] allowed a system of two inspiraling blackholes to be evolved through merger and the ringdown ofthe remnant black hole [5, 6, 7, 8, 9, 10, 11, 12, 13].Studying the late dynamical evolution of these inspi-raling compact binaries is important because they areamong the most promising source of gravitational wavesfor the network of laser interferometric detectors suchas LIGO and VIRGO. The detection of these gravita-tional waveforms (GW) is important for testing generalrelativity in the strong field limit. Moreover, these de-tectors can extract from the waves physical data aboutthese sources such as the component masses and spinsand the orbital eccentricity. For an unbiased extractionof these parameters, a large bank of accurate waveformsneeds to be constructed. Numerical relativity alone can-not compute all the waveforms needed because of thecomputational cost. Instead, the waveforms are basedon post-Newtonian (PN) approximations [14, 15].The post-Newtonian approximation is a slow-motion,weak-field approximation to general relativity. In orderto produce a post-Newtonian waveform, the PN equa-tions of motion of the binary are solved to yield explicitexpressions for the accelerations of each body in termsof the binary’s orbital frequency Ω [14, 16, 17, 18, 19,20, 21, 22, 23, 24, 25]. Then solving the post-Newtonian wave generation problem yields expressions for the gravi-tational waveform h and the gravitational wave flux F interms of radiative multipole moments [26]. These radia-tive multipole moments are in turn related to the sourcemultipole moments, which can be given in terms of therelative position and relative velocity of the binary [27].Instead of comparing the post-Newtonian waveform witha numerical waveform along a certain direction with re-spect to the source, the comparison can be done in all di-rections by decomposing the waveform in terms of spher-ical harmonic modes. For an equal-mass non-spinningbinary, the (2 ,
2) mode h [28, 29, 30, 31] is often usedto compare numerical and post-Newtonian waveforms be-cause it is the dominant mode. Its time derivative ˙ h isused to compute the gravitational wave flux. The result-ing expressions for the orbital energy E , the gravitationalenergy flux F and the amplitude h are given as Taylorseries of the frequency-related parameter x = ( M Ω) / , (1)where M is the total mass of the binary and G = c = 1.The invariantly defined “velocity” v = x / , (2)another dimensionless parameter, is often used in writingthese Taylor series.Computing PN series to high order is difficult and timeconsuming. Since the various PN expressions are givenas slowly convergent Taylor series, the Pad´e transforma-tion [32, 33] was suggested in Ref. [34] to accelerate theconvergence of these series. The Pad´e transformation, P mn , consists of writing a Taylor series, T k , of order k asthe ratio of two polynomials, one of order m in the nu-merator, and another of order n in the denominator, such n Exp n ( v ) P m + ǫm [Exp n ( v )]0 1.0000000 1.00000001 3.0099999 3.00999992 5.0300499 -401.00003 6.3834834 9.13136364 .0635838 .06014925 .3369841 n = P nk =0 v k /k ! of the exponential function Exp( v ) and its Pad´eapproximant Exp m + ǫm at v = 2 . m = ⌊ n/ ⌋ . The Pad´e ap-proximant converges to six significant figures while the Taylorseries converges to four significant figures at v = 2 .
01. Theerror between the exact value of the exponential, 7 . ( v = 2 .
01) is 6 × − , whilethe error between the Taylor approximant Exp ( v = 2 . − . that m + n ≤ k . If well-behaved, this method acceleratesthe convergence of a Taylor series as the order of the Pad´etransformation, m + n , is increased. For example, in Ta-ble I we compare the convergence of the Taylor expansionof the exponential function, Exp n ( v )( ≡ e v ), at order n to its Pad´e approximant, Exp m + ǫm ( v ) = P m + ǫm [Exp n ( v )],along the diagonal, where m = ⌊ n/ ⌋ and ǫ = 0 or 1.After twelve terms ( n = 11), the last two partial sumsof the Taylor expansion converge to 4 significant figures.However, the last two Pad´e approximants Exp ( v ) andExp ( v ) converge to 6 significant figures. The error be-tween the exact value of the exponential, 7 . ( v = 2 .
01) is 6 × − , whilethe error between the eleventh order partial sum and theexact value is 10 − . Fig. 1 shows the convergence of theTaylor expansion of the exponential function and its Pad´eapproximant.The hope of accelerating the convergence of the post-Newtonian Taylor series of the energy and flux motivatedthe use of their Pad´e approximants to construct Pad´eapproximant waveforms [34, 35, 36, 37, 38, 39, 40, 41,42, 43, 44, 45, 46]. If these resummation techniques ac-celerate the convergence of the Taylor series in PN ap-proximations, the range of validity of PN approximationssuggested by Ref [47] could be extended. Moreover, thework of Refs. [48, 49] in the test mass limit motivatedthe addition of a simple pole to the flux F of a binarysystem as the bodies approach the light ring orbit. Bymathematical continuity, the existence of a pole in theequal mass case was anticipated [34].More recently, waveforms are constructed by includingthese ideas in Effective One Body (EOB) models. TheEOB approach [35, 36, 38, 39, 40, 41, 42, 43, 44, 50, n -6-4-202 L og | E xp - E xp ( . ) | Exp n Exp mm+ ε FIG. 1: Convergence of the Taylor expansion, Exp n = P nk =0 v k /k ! of the exponential function Exp( v ) and its Pad´eapproximant Exp m + ǫm at v = 2 . m = ⌊ n/ ⌋ . The Pad´eapproximant converges faster than the Taylor series.
51, 52, 53, 54, 55, 56, 57] aims at providing an accurateanalytical description of the motion and radiation of coa-lescing binary black holes. The approach consists of threeseparate ingredients: 1) a description of the conservativeHamiltonian part of the dynamics ˆ H , 2) a formulationof the radiation reaction force F from the radiated flux F and 3) an expression of the GW waveform amplitudeemitted by the coalescing binary system (i.e h ).The flux plays an important role in approximating theradiation reaction force F in the EOB models [41, 58, 59].The leading-order radiation reaction force F [60, 61, 62]enters the equations of motion at 2.5PN order. Sincethe equations of motion are known only to 3.5PN order,one has to rely on the assumed balance between energyloss in the system and radiated flux at infinity [63, 64]to generate an approximate expression of the radiationreaction force at 3.5PN order beyond the leading term.Ref. [65] computes the GW energy flux and GW fre-quency derivative from a highly accurate numerical sim-ulation of an equal-mass, non-spinning black hole bi-nary. By assuming energy balance, the (derivative ofthe) center-of-mass energy is estimated. These quan-tities are then compared to the numerical values usingvarious Taylor, Pad´e and EOB models. The main goalof Ref. [65] is taking a set of well-established proposalsin the literature for approximating waveforms and seeinghow well they work in practice. Another goal of Ref. [65]is to examine some modifications of those proposals. Themain goal of this paper, by contrast, is to show that akey ingredient in those proposals does not appear to benecessary.In Ref. [66], Blanchet gave an argument that Pad´eand EOB resummations are unjustified because for twocomparable-mass bodies there is no equivalent of theSchwarzschild light-ring orbit at the radius r = 3 M . Hisargument is based on the PN coefficients of the binary’senergy and their relation to predicting the innermost cir-cular orbit. He finds that the radius of convergence of thePN series, which is related to the radius of the light-ringorbit, is around one (instead of 1 / E and the flux F , they proposed a new class ofwaveforms called P -approximants, based on three essen-tial ingredients. The first step is the introduction ofnew energy-type (Eq. 4) and flux type (Eq. 16) func-tions, called e ( v ) and f ( v ) respectively. The second stepis to Pad´e-approximate the Taylor expansion of thesefunctions. The third step is to use these Pad´e trans-forms in the definition of the energy E (Eq. 6) and Pad´e-approximated flux (Eq. 20). The last step is to constructeither the Pad´e-approximated waveform as in Sec. IV orthe EOB waveform as in Sec. V. Schematically, the sug-gested procedure is summarized by the following map: (cid:20) E n , F n (cid:21) → (cid:20) e n , f n (cid:21) → (cid:20) e mn , f mn (cid:21) → (cid:20) E ( e mn ) , F ( f mn ) (cid:21) → h . (3)Our notation is to denote by T mn ( x ) the Pad´e approx-imant of a k -th order Taylor series T k ( x ) with an m -thorder polynomial in the numerator and an n -th orderpolynomial in the denominator such that m + n ≤ k , i.e.the Pad´e approximant of e k ( x ) is e mn ( x ).In Section II, we compare the 3PN Taylor series of theenergy function to its possible Pad´e approximants usingthe intermediate energy function e ( x ), as suggested byDamour et al. [34]. We compute the last stable orbitfrequency, defined as the frequency for which the energyreaches a minimum as a function of frequency, and alsothe poles of the energy in the complex plane correspond-ing to each possible Pad´e approximant. The large vari-ation of last stable orbit frequency and poles does not PN order F n F m + ǫm F n F m + ǫm v = 0 . x = 0 . F n , F m + ǫm , F n and F m + ǫm are given in Eqs. 14, 15, 18 and 20 respectively. Even in thetest mass limit, the Pad´e approximant of the flux fails to con-verge faster that its 5 . v = 0 .
2. After 12 terms, only about 4 or 5 signifi-cant digits seem reliable for the Taylor expansions and theirPad´e approximants. The lack of improvement in the conver-gence of the Pad´e approximants should be contrasted withthe example in Table I. suggest good convergence of the Pad´e-approximated in-termediate energy function e ( x ). The energy function E ( x ) is strongly dependent on the choice of the Pad´eapproximant of e ( x ). Accordingly, the Pad´e waveformwill also be strongly dependent on the choice of the Pad´eapproximant.In Section III, we present two possible methods for cal-culating the Pad´e approximant of the flux function. Thefirst method simply takes the Pad´e approximant of theTaylor series treating the logarithmic contribution as con-stant. Following [34], the second method adds a pole tothe Taylor series, factors out the logarithmic contributionto the series, and then computes the Pad´e approximantof the resulting Taylor series. We test the convergenceof the Pad´e approximant for both methods versus theirTaylor series. We find that the Pad´e approximants of theflux do not converge any faster than their Taylor coun-terpart.A simple example that illustrates the problem is shownin Table II. There we compare the partial sums of theTaylor series for the flux with the corresponding Pad´eapproximants in the test mass limit. The four flux func-tions F n , F m + ǫm , F n and F m + ǫm are given in Eqs. 14, 15, 18and 20 respectively. Even for a relatively small value of x ,namely x = 0 .
04 ( v = 0 . E . We also use the last stable orbit from the 3PN en-ergy Taylor series E . The results are not very sensi-tive to this choice. We compare the Pad´e waveforms toa 15-orbit numerical waveform in the equal mass, non-spinning quasi-circular case [67]. The phase differencein these comparisons ranges between 0.05 and a few ra-dians for well-defined Pad´e approximants (not having apole in the frequency domain of interest) when the match-ing of the numerical and Pad´e waveforms is done at thegravitational wave frequency M ω = 0 . h Taylor series (see Refs. [35, 67] for thedefinition of these Taylor series).In Section V, based on the results of the previous sec-tions, we design a simple EOB model (closely related tothe ET EOB model of Ref. [40]) using the Taylor series ofthe flux. We add one unknown 4PN term that we fit forby maximizing the agreement between the EOB modelwaveform and the numerical waveform. The model doesnot require adding a pole to the flux, nor an a prioriknowledge of the last stable orbit from the energy func-tion. This simple EOB model, with only one parameterto fit for, agrees with the numerical waveform to within0.002 radians (3 × − cycles). (This is six times smallerthan the claimed numerical accuracy of [39], smaller byan even larger factor than the claimed numerical accu-racy of [45], and 25 times smaller than the gravitationalwave phase uncertainty of the numerical waveform. SeeTable III in Ref [67] for more details.) This model agreeswith the numerical waveform better than any previouslysuggested Taylor, Pad´e or EOB waveform. II. ENERGY FUNCTION
Damour et al. [34] introduced a new energy-type func-tion e ( x ), where x is the PN frequency related parameter.This assumed more “basic” energy function e ( x ) is con-structed out of the total relativistic energy E tot ( x ) of thebinary system. Explicitly e ( x ) ≡ (cid:18) E − m − m m m (cid:19) − , (4) where m , m are the masses of the bodies. The totalrelativistic energy function E tot is related to the postNewtonian energy function E ( x ) through E tot ( x ) = M [1 + E ( x )] , (5)where M is the total mass ( M = m + m ). Solving for E ( x ) in terms of e ( x ) using Eqs. (4) and (5), we get[34] E ( x ) = n ν hp e ( x ) − io / − , (6)where the symmetric mass ratio is ν = m m /M . Theorbital energy function E ( x ) is known as a Taylor series E k up to 3PN order as a function of x and ν [15] E ( x ) = − ν x (cid:26) −
112 (9 + ν ) x − (cid:0) − ν + 13 ν (cid:1) x + h − (cid:18) − π (cid:19) ν − ν − ν i x (cid:27) . (7)Using the above equations, we compute the Taylor seriesexpansion, e k ( x ), of e ( x ) up to 3PN order: e ( x ) = − x (cid:26) − (1 + 13 ν ) x − (3 − ν ) x − h (cid:0) − π (cid:1) ν + 10336 ν − ν i x (cid:27) . (8)In the test mass limit ( ν → e ( x )coincides with the Pad´e approximant P ( x ) of its Taylorexpansion in Eq. (8) e ( x ; ν →
0) = − x − x − x . (9)This quantity has a pole at x pole = 1 /
3. The orbitalenergy is then E ( x ; ν →
0) = ν r − x − x − x − ! , (10)and it derivative is dE ( x ; ν → dx = − ν − x − x ) / . (11)The last stable orbit occurs where dEdx = 0 , (12)so in the limit ν → x lso = 1 /
6. On the grounds of mathematical continu-ity between the test mass limit ν → E n e r gy E E E E E E PSfrag replacements x FIG. 2: Post Newtonian Energy at 3PN and its Pad´e approx-imants for the case ν = 1 /
4. The plot includes the high valueof x lso = 0 .
36, the numerical data available is at x = 0 . E , E , E and E vary significantly althoughthey all correspond to the 3PN Taylor series of the energyfunction. E is very different from the other functions, whichsuggests a poorly convergent Pad´e approximant. mass ratio case, Damour et al. [34] argued that the ex-act function e ( x ) should be meromorphically extendablein at least part of the complex plane and should have asimple pole on the real axis. They suggested that Pad´eapproximants would be excellent tools for giving accuraterepresentations of functions having such poles.Once we know the Taylor series of the new energy func-tion, e k ( x ), we compute its Pad´e approximant e mn ( x ),with m + n ≤ k . The Pad´e-approximated energy, E mn ( x )is obtained by replacing e ( x ) in Eq. (6) with e mn ( x ). Inthe equal mass case ( ν = 1 / e k ( x ). The most interesting Pad´eapproximants have a maximal sum of their indices, sincethey should be closest to the unknown exact function ifthe Pad´e resummation is converging. In Fig. 2, we showa plot of the PN energy function E ( x ) and its Pad´eapproximants E , E , E , E and E as a function of x .Although the Pad´e approximants of the energy are ofmaximal order, they differ significantly. Good conver-gence of the Pad´e approximants requires good agreementbetween approximants of the same order n + m , if thereis no pole in the region of interest (0 < x . . E or E . Although both have the same order andare equally close to the diagonal, the difference betweenthese functions is quite large.In Table III we compute the locations of the poles andthe last stable orbits for all of these Pad´e approximants.The ill-convergence of the Pad´e transform is again seen bylooking at the variation of the last stable orbit positions.In Table III for example, x lso of E differs by about 8% from x lso of E . Moreover, for finite ν , the poles are allcomplex or not in the interval [0 ,
1] except for the case x pole = 52 / E . There is no reason why this should be the“exact” pole that should be used in the formalism, sincenone of the third-order Pad´e approximants of the 3PNenergy has a physical pole.In summary, using Pad´e approximants for the energyfunction in the equal mass case does not seem to pro-vide any benefit. The differences between the variousPad´e approximants of the energy are large. The quanti-ties x pole and x lso do not show any regular behavior thatcould be a sign of a physical pole that could be found byusing the Pad´e transform. III. FLUX FUNCTION
The general form of the PN flux at order N is: F ( v ) = 325 ν v × F N , (13)where the normalized flux F is a Taylor expasion in v with logarithmic terms F N ( v ) = N X k =0 A k v k + N X k =6 B k v k ! log v , (14)where the post-Newtonian coefficients A i and B i arefunctions of the mass ratio parameter ν . They are givenin the test mass limit in Ref. [69] and in the equal massquasi-circular case in Ref. [15]. The flux series has alogarithmic contribution starting at 3PN. Pad´e approx-imants, however, are well defined only for pure polyno-mials. Two possible methods are therefore used to com-pute the Pad´e approximant of the flux. The first methodsimply treats the logarithmic terms as constants and re-sums the series as a pure polynomial such that the Pad´e-approximated flux F mn is F mn ( v ) = P mn (cid:2) F N ( v ) (cid:3) . (15)The second method, suggested by Ref. [34], defines anew flux function f by adding a pole, factoring the log-arithmic terms from the series, and finally computingthe Pad´e approximant of the pure polynomial. Sincewe would like to check the convergence of the Pad´e-approximated flux versus its Taylor series, we sketch thedefinitions of the various functions involved. Accordingto Ref. [34], two ideas are needed for a good representa-tion of the analytic structure of the flux. First, since inthe test mass limit F is thought to have a simple pole atthe light ring [49], one might expect it by continuity tohave a pole in the comparable mass case. This motivatesthe introduction of the following factored flux function, f ( v ; ν ): f ( v ; ν ) ≡ (cid:18) − vv pole ( ν ) (cid:19) F ( v ; ν ) , (16) Energy x pole x lso E − . E /
109 = 0 .
477 0 . E − . E − .
41 0 . E . ± . i . E . ± . i , − .
696 0 . ν = 1 /
4. The poles x pole and laststable orbit frequency of the function E mn ( x ) depend signif-icantly on which Pad´e approximant is constructed from theTaylor series e k ( x ). The only physical pole is x pole = 52 / where v pole is the pole of the Pad´e-approximanted energyfunction used.Second, the logarithmic term that appears in the fluxfunction needs to be factored out so we can use the stan-dard Pad´e transformation. After factoring the logarith-mic terms out, the flux function f becomes f n ( v ; ν ) = " vv lso N X k =6 ℓ k v k ! × N X k =0 f k v k ! , (17)where the coefficients l k and f k are given in Ref. [34], and v lso is the velocity of the last stable orbit of the Pad´e-approximated energy. Then the Taylor series of the fluxwith a pole is defined as F n ( v ; ν ) ≡ f n ( v ; ν )1 − v/v pole ( ν ) . (18)The Pad´e approximant of the intermediate flux function f ( v ) is defined as f mn ( v ) ≡ " vv lso ( e mn ; ν ) N X k =6 ℓ k v k ! × P mn " N X k =0 f k v k , (19)where v lso ( e mn ; ν ) denotes the last stable orbit velocityfor the Pad´e approximant P mn (cid:2) e ( x ) (cid:3) . Finally, the cor-responding Pad´e approximant of the flux F ( v ) is givenby F mn ( v ; ν ) ≡ f mn ( v ; ν )1 − v/v pole ( e mn ; ν ) , (20)where v pole ( e mn ; ν ) denotes the pole velocity defined by e mn ( x ). A. Flux for the test mass case
The exact gravitational wave luminosity F is notknown analytically in the test particle limit. It hasbeen computed numerically by Poisson [68]. The post-Newtonian expansion of the flux is known in the test masslimit to 5.5PN order [69]. This allows us to test the rate ofconvergence of the Taylor series of the normalized flux F n (Eq. 14) and its Pad´e-approximant F mn constructed treat-ing the logarithmic term as a constant (Eq. 14). We alsotest the convergence of the flux function F n (Eq. 18) andits Pad´e approximant F mn (Eq. 20). These convergencetests use the known values v pole = 1 / √ v lso = 1 / √ v = 0 .
2. The four fluxfunctions F n , F m + ǫm , F n and F m + ǫm are given in Eqs. 14,15, 18 and 20 respectively. We use the Pad´e approximantalong the diagonal P m + ǫm where ǫ = 0 or 1. The ratesof convergence of the Taylor expansion and its Pad´e ap-proximant are nearly equal for the two methods, whetheror not we include a pole. As the PN order increases, theTaylor series and its Pad´e approximant alternate in whichprovides a better fit to the numerical data for the flux.For example, at 2PN order the Taylor flux with a pole(Eq. 18) fits the numerical data the best. At 2.5 and 3PN order the Pad´e approximant of the flux F mn (Eq. 20)fits the numerical data the best, while at 3.5 and 5PN or-der the Taylor series of the flux (Eq. 14) is the best. At5.5PN the Pad´e approximant of the flux (Eq. 20) givesthe best agreement. The results are similar for othervalues of v . No method has the best convergence rate.According to Pad´e theory, the convergence of the Pad´eapproximant is best along the diagonal, but it is equallygood along the off-diagonal terms if no pole exists in theregion of interest (i.e. no zeroes appear in the denom-inator of the Pad´e approximant.) For this reason, weshow the error between all the possible maximal Pad´e-approximated fluxes F − nn (Eq. 14) and the numericalflux for three values of v (= 0 .
2, 0 .
25, 0 .
35) ( x = 0 . .
06, 0 .
12) in Fig 4. The 5.5PN Taylor series, denoted by F , fits the exact numerical data better than the Pad´eapproximants F , F , F , F . In the other cases, thePad´e approximants provide a better agreement (i.e. F , F , F and F ) for the three values of v . This suggeststhat the Pad´e approximation should only be used to sug-gest a fitting formula for the numerical data, since thereis no internal self-consistency in the agreement. The off-diagonal approximants do not show any regular patternof convergence to the numerical data nor are they betterthan the Taylor series. PN order -6-5-4-3-2-1 L og | F - F N R | F n F mm+ ε F n F mm+ ε FIG. 3: Convergence of the flux approximations in the testmass limit for v = 0 .
2. The four flux functions F n , F m + ǫm , F n and F m + ǫm are given in Eqs. 14, 15, 18 and 20 respectively. ThePad´e approximants do not converge faster than their Taylorseries counterparts. The Pad´e and Taylor series alternate atproviding the best agreement with the “exact” data as thePN order increases. Contrast the behavior here with Fig. 1. B. Flux for the equal mass case
For binaries of comparable mass on a quasi-circularorbit, the flux is known only to 3.5PN order [15]. InRef. [65] for a quasi-circular non-spinning binary, thenumerical flux was computed by integrating the spin-weighted spherical harmonic components of the Weylscalar Ψ . The numerical flux data we use in this paperwas provided by Harald P. Pfeiffer and Michael Boyle.The estimated error in measuring the flux data was about0.2%. The velocity range for the simulation was from v ∼ .
26 ( x ∼ .
06) to v ∼ . x ∼ . ω /
2, where ω is the wave frequency of the ˙ h mode. Instead, inTable IV we compare the convergence of the four fluxfunctions F n , F m + ǫm , F n and F m + ǫm (defined in Eqs. 14,15, 18 and 20 respectively as a function of PN order)for v = 0 . x = 0 . v pole = 0 .
69 ( x pole = 52 / v lso = 0 .
50 ( x lso = 0 . E . Theconvergence does not depend on these values althoughthe flux values listed in Table IV do depend somewhat L og | F n11 - n - F N R | v = 0.20v = 0.25v = 0.35 PSfrag replacements n v FIG. 4: Error between maximal Pad´e approximants of the flux F (Eq. 15) and the numerical flux in the test mass limit at v =0 . , . , .
35. The 5.5PN Taylor series, denoted by F , fitsthe exact numerical data better than the Pad´e approximants F , F , F , F . In the other cases, the Pad´e approximantsprovide a better agreement (i.e. F , F , F , F and F , F and F ). on the values of v pole and v lso . We choose a medium ve-locity ( v = 0 .
2) to make the rate of convergence clear.At 3.5PN order, all four flux functions agree to 2 signif-icant figures. However after seven terms, F n convergedto 3 significant figures, F m + ǫm converged to 4 significantfigures, while F n and F m + ǫm converged to 2 significantfigures. The flux function F m + ǫm converged to one addi-tional significant figure over F n , however F m + ǫm cannotreliably be considered more accurate than F n because itconverges to a slightly different value. The Pad´e approx-imants do not seem to converge to a larger number ofsignificant figures than the Taylor flux function F n .In Fig. 5, we plot the numerical normalized flux F NR ,the 3.5PN flux F . and the maximal Pad´e-approximatedflux functions F , F , F , F and F ( ≡ F ). Although F . diverges from the numerical flux early at v ∼ . F , F and F .The quantity F has a pole and fails to capture the nu-merical flux behavior completely. The quantity F is bydefinition the Taylor flux with a pole, F . This functionshows that adding a pole to the Taylor expansion of theflux F . degrades the fit with the numerical flux. More-over, the numerical flux does not suggest the existence ofa pole at a large velocity ( v ∼ . v ∼ .
4. Adding a pole does not seem a usefulidea in this case at least. On the other hand, F and F are a better fit to the numerical data during most of thevelocity range of the 15-orbit data. The flux function F is especially a good fit to the numerical normalized fluxat high velocities. However, even though F and F are PN order F n F m + ǫm F n F m + ǫm v = 0 . x = 0 . v pole = 0 .
69 ( x pole = 52 / v lso = 0 .
50 ( x lso = 0 . F n , F m + ǫm , F n and F m + ǫm are given in Eqs. 14, 15, 18 and 20 respectively.At 3.5PN order, all four flux functions agree to 2 significantfigures. After seven terms, F n converges to 3 significant fig-ures, F m + ǫm converges to 4 significant figures, while F m and F m + ǫm converge to 2 significant figures. a good fit to the numerical flux during the last 15-orbitinspiral before merger, there is no guarantee that this istrue at low velocities. IV. PAD´E WAVEFORMS
The construction of the post-Newtonian waveforms re-quires solving the post-Newtonian equations describingthe motion of the binary and the generation of gravita-tional waves. Substituting the orbital evolution predictedby the equations of motion into the expressions for thewaveform would not generate waveforms accurate enoughfor matched filtering in detecting gravitational waves [70].To compute the waveform at 3PN order, it is necessary tosolve the equations of motion at 5.5PN order, because theradiation reaction contributes to the equations of motionstarting at 2.5PN order. However, for a non-spinning bi-nary of equal mass and on a circular orbit, accurate wave-forms at 3PN order can be constructed under two furtherassumptions. The first assumption is that the binary fol-lows a slow adiabatic inspiral. The second assumption isthat of energy balance between the orbital binding en-ergy and the energy emitted by the gravitational waves,where the energy balance equation is defined as dEdt = − F. (21)The procedure of constructing the standard Pad´ewaveforms [34] is similar to one used to construct theTaylorT1 waveforms in Ref. [34, 67]. The main differenceis the use of Pad´e approximants of the energy and flux tocompute the orbital phase, as described in Secs. II andIII, instead of their Taylor expansions. The orbital phaseused in the Pad´e waveforms is obtained by numericallyintegrating d Ω dt = 325 ν v F mn dE kl /d Ω . (22) N o r m a li ze d F l ux F NR F F F F F F PSfrag replacements v FIG. 5: Normalized flux for an equal mass nonspinning bi-nary. We plot the numerical flux F NR , the 3.5PN flux F . and the maximal Pad´e-approximated flux functions F , F , F , F and F ( ≡ F ). The early noise is caused by the junkradiation. The fraction on the right side of Eq. (22) is retained asa ratio of the Pad´e approximants of the post-Newtonianexpansions, and is not expanded further before numericalintegration. The waveform is produced by substitutingthe orbital phase into the spherical harmonic mode h of the post-Newtonian waveform, which is known up to3PN order [28, 29, 30, 31].Given the expressions for the Pad´e-approximated en-ergy and flux in Sections II and III, and the Taylor seriesof the waveform amplitude [28, 29, 30, 31], there is stilla set of choices that must be made in order to producea Pad´e-approximated waveform that can be compared toour numerical waveform. These include1. The Pad´e approximant of the orbital energy, E kl .2. The flux function and its Pad´e approximant F nm .3. The velocity of the pole and the last stable orbit, v pole and v lso .4. The PN order through which terms in the waveformamplitude are kept. A. Procedure
We consider numerical gravitational waves extractedwith the Newman-Penrose scalar Ψ , using the same pro-cedure as in [71]. To minimize gauge effects, we compareits (2 ,
2) component extrapolated to infinite extractionradius according to Ref. [67]. The extracted waveform issplit into real phase φ and real amplitude A , defined by m | δ φ | F :E F :E F :E F :E F :E FIG. 6: Phase difference between the 3 and 3.5 PN Pad´eapproximated and numerical waveforms matched at the wavefrequency Mω = 0 .
1. We use the Pad´e-approximated flux F m − m (Eq. 20) and energy E kl . We include in the figure thewaveforms using the Pad´e-approximated flux F using m = −
1. There is no entry for m = 4 since the Pad´e-approximatedflux F has a pole in the frequency range of the simulation. Ref. [67] as Ψ ( r, t ) = A ( r, t ) e − iφ ( r,t ) . (23)The gravitational-wave frequency is given by ω = dφdt . (24)The spherical harmonic component (2,2) of Ψ is thencompared to the numerically twice-differentiated post-Newtonian expression of h , A , as in Ref.[67]. Follow-ing [67, 72, 73], the matching procedure needed to setthe arbitrary time offset t and the arbitrary phase offset φ is done by demanding that the PN and NR gravita-tional wave phase and gravitational wave frequency agreeat some fiducial frequency ω M . B. Results
In this section, we compare the numerical waveformto the Pad´e waveforms corresponding to the 3.5 PN or-der of energy and flux using the 3PN Taylor series ofthe post-Newtonian amplitude A . The energy and fluxfunctions used are those suggested by Ref. [34]. We donot generate all possible waveforms using different Pad´eapproximants of the energy or the flux at low PN orders,since all these resummed series showed no improvementin the convergence rate.As introduced in Sec. II, we use the Pad´e-approximated energy E , E , E and E corresponding to the PN Taylor series of the energy, and the Pad´e-approximated energy, E , corresponding to its 2PN Tay-lor expansion. For the flux, the diagonal Pad´e approxi-mant F is used in addition to all possible Pad´e approx-imants of flux at 3.5PN order, F m − m , where 0 ≤ m ≤ v lso and v pole as discussed in Sec. III. The value v pole =52 /
109 is used. We also tested varying the pole loca-tion, but found that we could not improve the agreementsignificantly.From Table III, any value of the velocity of the laststable orbit could be used. We use the 3PN value v lso =0 .
254 and also use v lso = 0 . E is employed in the constructionof the waveform. In the remaining cases, we use v lso =0 .
254 since it is quite close to the estimates from otherPad´e approximants of the energy. The effect of changingthe value of v lso is not significant compared to changingthe order of the Pad´e approximant for the energy or theflux.To do the comparison, we match the Pad´e-approximated and numerical waveforms at the wave fre-quency M ω = 0 .
1. Then we measure the maximumphase difference between the numerical waveform andeach of these Pad´e waveforms during the inspiral whenthe numerical wave frequency is between
M ω = 0 . M ω = 0 . E kl and flux F m − m . On the same Figure, we include phasedifferences for the waveforms generated using the Pad´e-approximated flux F under the m = − E is used, the phase error ranges between 2and 5 radians as m increases from − .
05 to 2 . m = 7) resulted in alarge phase difference ranging between 1 and 1.5 radians.The diagonal Pad´e term F of the flux generates similarphase differences, ranging from 0 .
06 to 0 . V. SIMPLE EOB MODEL
We have described the failure of the Pad´e resumma-tion techniques to accelerate the convergence of any PNTaylor series, the absence of any signature of a pole inthe flux in the equal mass case, and the erratic patternof agreement between the Pad´e waveforms and the nu-merical waveform. It seems one might as well simply usethe Taylor series at all steps of computing waveforms.Also it does not seem that the parameters v pole and v lso are useful. In this section, we show how to get goodagreement with the numerical waveform by using a sim-ple EOB model. The only parameter we introduce andfit for is an unknown 4PN contribution to the flux. A. EOB waveforms
The EOB formalism [53] is a non-perturbative ana-lytic approach that handles the relative dynamics of tworelativistic bodies. This approach of resumming the PNtheory is expected to extend the validity of the PN resultsinto the strong-field limit. The procedure for generatingan EOB waveform follows closely the steps in Sec. IV. In-stead of using the energy balance equation, we computethe orbital phase by numerically integrating Hamilton’sequations. The EOB waveform is generated by substitut-ing the orbital phase into the waveform amplitude A at3PN order. The two fundamental ingredients that allowcomputing the orbital phase are the real Hamiltonian ˆ H and the the radiation reaction F φ . B. Hamilton’s equations
In terms of the canonical position variables r and φ andtheir conjugate canonical momenta p r and p φ , where r isthe relative separation and φ is the orbital phase, the realdynamical Hamiltonian is defined as [54]:ˆ H = 1 ν p ν ( H EOB − , (25)where H EOB = s A (cid:16) p φ r + p r B + 2 ν (4 − ν ) p r r (cid:17) , (26)and where the radial potential A function is defined asthe series A = 1 − r + 2 νr + (cid:16) − π (cid:17) νr . (27) The Taylor series of the A function is replaced by its Pad´eapproximant A . Here the Pad´e approximant is not usedto accelerate the convergence of the Taylor expansion of A . Instead, it leads to the existence of a last stable orbit(see Ref. [39] and references therein). Otherwise, theEOB Hamiltonian is non-physical for the last few orbits;the orbital frequency stays nearly constant for severalorbits before merger. For the B function, the Taylorexpansion suffices: B = 1 A h − νr + 2(3 ν − νr i . (28)Then Hamilton’s equations of motion are given in thequasi-circular case by ∂ t r = ∂ p r ˆ H , (29) ∂ t φ = ∂ p φ ˆ H , (30) ∂ t p r = − ∂ r ˆ H , (31) ∂ t p φ = −F φ , (32)where F φ is the radiation reaction in the φ direction rep-resenting the nonconservative part of the dynamics. InEq. 32, ∂ φ ˆ H = 0 since ˆ H is independent of φ . The radi-ation reaction is deduced from the post-Newtonian fluxas in Refs. [41, 58, 59] F φ = F + F v νv . (33)In this equation, we have introduced an unknown 4PNflux term, F , the only parameter that we fit for in thisEOB model. C. Initial conditions
To integrate Hamilton’s equations, we need appropri-ate initial conditions for a quasi-circular orbit. Refs. [41,42, 52] indicate how to define some “post-adiabatic” ini-tial conditions. However, these initial conditions do notgenerate an orbit with as low an eccentricity as the nu-merical simulation, roughly 5 × − . At a given radius r ,starting from the post-adiabatic initial conditions of p r and p φ , we therefore reduce the eccentricity iterativelyin two steps. The first step includes evolving Hamilton’sequations in the conservative regime ( F = 0) and it-eratively changing the value of p φ until the eccentricitymeasured from the evolution of the orbital separation isof the order 10 − . The second step is based on evolvingthe nonconservative Hamilton’s equations with the 4PNflux and iteratively changing the p r momentum until theeccentricity is again of the order 10 − . This circulariza-tion procedure is repeated as we iterate F to maximizethe agreement between the waveforms. D. Best Fit of F To find the best fit for F , we iteratively solve for theminimum in the phase difference between the numerical1 (t-r*)/m | A E O B - A N R | / A N R δ φ (r a d i a n s ) FIG. 7: Phase and amplitude differences between the EOBwaveform and the numerical waveform. After fitting for thebest value of F , the phase difference is less than 0 . r ∗ is the tortoise coordinate definedin [67]. and EOB waveforms. The waveforms are matched as inSec. IV at the wave frequency mω = 0 .
1, and the phasedifference is defined as the maximal phase difference dur-ing the inspiral phase up to the wave frequency mω = 0 . F = − .
75 correspondingto the initial conditions r = 17, φ = 0, p r = − . p φ = 4 . F changes the max-imal phase difference from less than 0.002 radians toabout 0.01 radians. Note that without adding the fittingparameter F , the phase difference is about 1.7 radiansduring the 15-orbit inspiral. E. Results
In the upper panel of Fig. 7, we plot the phase dif-ference between the numerical waveform and the EOBwaveform computed using the 3PN Taylor series of theamplitude A . The phase difference is less than 0 .
002 ra-dians after maximizing the agreement between the wave-forms in the region where mω ≤ .
1. The early noise isdue to junk radiation at the early stage of the numeri-cal simulation as described in Sec II C of Ref. [67]. Thephase uncertainty in the simulation was estimated to be0 .
05 radians; See Table III in [67].In the lower panel of Fig. 7, we plot the relative differ-ence between the amplitude of the numerical waveformand the EOB waveform. The EOB waveform amplitudedoes not fit the numerical waveform amplitude as wellas the wave phase does. This is expected because thewaveform amplitude is known to 3PN order only, andno free parameter in the amplitude was fitted for. The agreement between the amplitude of this EOB model andthe numerical waveform is similar to the agreement be-tween the amplitude of TaylorT4 3.5/3.0 and the numer-ical waveform in Fig. 21 in [67].This EOB model is a modification of the ET EOBmodel of Ref. [74]. It fits the numerical phase very wellwithout using the Pad´e resummation techniques nor apole in the flux.Even though we have found very good agreement be-tween the waveforms, these results only suggest that theEOB model is a very good fitting model. Moreover, hav-ing fit a particular waveform, there is no guarantee themodel will have predictive power for a more general case.
VI. CONCLUSIONS
Convergence tests show that none of the Taylor seriesin the PN approximation, such as the energy or the flux,could be replaced by a Pad´e approximant that convergesfaster. Other attempts we tried to accelerate the conver-gence of these series also failed, as for example using theLevin method to accelerate convergence [33]. As a result,more reliable waveforms could not be constructed usinga Pad´e resummation scheme. Moreover, the Pad´e wave-forms also do not fit numerical simulation data betterthan the Taylor waveforms. Thus, they do not seem to bebetter than the Taylor waveforms in building templatesfor waveforms. This conclusion is independent of thePad´e approximants used to test the convergence. Takingfor example the sub-diagonal Pad´e approximant does notshow any improvement in the convergence rate. In addi-tion, this conclusion is independent of the numerical datawe used. We can simply take the highest PN order of theTaylor series or the Pad´e approximant and use it as the“exact” value of the function to test the convergence atlow frequency.Based on the dependence of the flux on the velocityin the equal mass case, we do not find it helpful to adda pole to the flux. Therefore, we recommend using Tay-lor series instead of the Pad´e approximant to generatewaveforms both in the time and frequency domains. Thesimple EOB model used in this paper agrees with thenumerical data very well; the phase difference during theinspiral is much less then the estimated phase uncertaintyin the numerical data. This model does not use Pad´e ap-proximants or poles except in one place to enforce a laststable orbit. Since Pad´e approximation does not acceler-ate the convergence of any PN Taylor series, there is noreason why it should estimate parameters better in dataanalysis of waveforms.
Acknowledgments
It is a pleasure to acknowledge useful discussions withEmanuele Berti, Michael Boyle, Alessandra Buonanno,Lee Lindblom, Harald P. Pfeiffer, Yi Pan, Mark A. Scheel2and Nicol´as Yunes. We thank Jihad Touma for helpfuldiscussions about Pad´e approximants, Harald P. Pfeifferand Michael Boyle for providing the numerical data ofthe flux in the equal mass case, and Eric Poisson forproviding the numerical data of the flux in the test mass limit. This work was supported in part by grants fromthe Sherman Fairchild Foundation to Cornell; by NSFgrants PHY-0652952, DMS-0553677, PHY-0652929, andNASA grant NNG05GG51G at Cornell. [1] F. Pretorius, Phys. Rev. Lett. , 121101 (2005).[2] F. Pretorius, Class. Quant. Grav. , S529 (2006).[3] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo-chower, Phys. Rev. Lett. , 111101 (2006).[4] J. G. Baker, J. Centrella, D. I. Choi, M. Koppitz, andJ. van Meter, Phys. Rev. Lett. , 111102 (2006).[5] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.Rev. D , 061501(R) (2006), gr-qc/0601091.[6] F. Herrmann, I. Hinder, D. Shoemaker, and P. Laguna,Class. Quantum Grav. , S33 (2007), gr-qc/0601026.[7] P. Diener, F. Herrmann, D. Pollney, E. Schnetter, E. Sei-del, R. Takahashi, J. Thornburg, and J. Ventrella, Phys.Rev. Lett. , 121101 (2006), gr-qc/0512108.[8] M. A. Scheel, H. P. Pfeiffer, L. Lindblom, L. E. Kidder,O. Rinne, and S. A. Teukolsky, Phys. Rev. D , 104006(2006), gr-qc/0607056.[9] U. Sperhake, Phys. Rev. D , 104015 (2007), gr-qc/0606079.[10] B. Br¨ugmann, J. A. Gonzalez, M. Hannam, S. Husa,U. Sperhake, and W. Tichy, Phys. Rev. D , 024027(2008), gr-qc/0610128.[11] P. Marronetti, W. Tichy, B. Br¨ugmann, J. Gonzalez,M. Hannam, S. Husa, and U. Sperhake, Class. QuantumGrav. , S43 (2007), gr-qc/0701123.[12] Z. B. Etienne, J. A. Faber, Y. T. Liu, S. L. Shapiro, andT. W. Baumgarte, Phys. Rev. D , 101503(R) (2007),arXiv:0707.2083 [gr-qc].[13] B. Szil´agyi, D. Pollney, L. Rezzolla, J. Thornburg, andJ. Winicour, Class. Quantum Grav. , S275 (2007), gr-qc/0612150.[14] L. Blanchet, T. Damour, and G. Esposito-Far`ese, Phys.Rev. D , 124007 (2004).[15] L. Blanchet, Living Rev. Relativity , 4 (2006).[16] P. Jaranowski and G. Sch¨afer, Phys. Rev. D , 7274(1998).[17] P. Jaranowski and G. Sch¨afer, Phys. Rev. D , 124003(1999).[18] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Rev.D , 021501(R) (2000), , 029903(E) (2000).[19] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Rev.D , 044021 (2001), , 029901(E) (2002).[20] L. Blanchet and G. Faye, Phys. Lett. A , 58 (2000).[21] L. Blanchet and G. Faye, Phys. Rev. D , 062005 (2001).[22] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Lett.B , 147 (2001).[23] Y. Itoh, T. Futamase, and H. Asada, Phys. Rev. D ,064038 (2001).[24] Y. Itoh and T. Futamase, Phys. Rev. D , 121501(R)(2003).[25] Y. Itoh, Phys. Rev. D (2004), 064018.[26] K. S. Thorne, Rev. Mod. Phys. , 299 (1980).[27] L. Blanchet, Class. Quantum Grav. , 1971 (1998).[28] L. E. Kidder, Phys. Rev. D , 044016 (2008),arXiv:0710.0614. [29] K.G. Arun, L. Blanchet, B.R. Iyer and M.S.S.Qusailah, Class. Quantum Grav. , 3771 (2004),arXiv:0404085[gr-qc].[30] K.G. Arun, L. Blanchet, B.R. Iyer and M.S.S. Qusailah,Class. Quantum Grav. , 3115 (2005).[31] L. Blanchet, G. Faye, B.R. Iyer and S. Sinha (2008),arXiv:0802.1249[gr-qc].[32] C. M. Bender and S. A. Orszag, Advanced MathematicalMethods for Scientists and Engineers: Asymptotic Meth-ods and Perturbation Theory, 2nd ed. (Springer, NewYork, 1999).[33] William H. Press, Saul A. Teukolsky, William T. Vetter-ling and Brian P. Flannery,
Numerical Recipes: The Artof Scientific Computing, 3rd ed. (Cambridge UniversityPress, New York, 2007).[34] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Phys.Rev. D , 885 (1998).[35] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Phys.Rev. D , 044023 (2001), , 029902(E) (2005).[36] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Phys.Rev. D , 027502 (2002), , 029901(E) (2005).[37] T. Damour, P. Jaranowski, and G. Sch¨afer, Phys. Rev.D , 084011 (2000), arXiv:gr-qc/0005034.[38] T. Damour and A. Nagar , Phys. Rev. D , 064028(2007).[39] T. Damour and A. Nagar , Phys. Rev. D , 024043(2008), gr-qc/0803.0915.[40] A. Buonanno, Y.-B. Chen, and M. Vallisneri, Phys.Rev. D , 024016 (2003), , 029903(E) (2006), gr-qc/0205122.[41] A. Buonanno and T. Damour, Phys. Rev. D , 064015(2000).[42] A. Buonanno, Y. Chen, and T. Damour, Phys. Rev. D , 104005 (2006).[43] A. Buonanno, Y. Pan, J. G. Baker, J. Centrella, B. J.Kelly, S. T. McWilliams, and J. R. van Meter, Phys.Rev. D , 104049 (2007), arXiv:0706.3732v2 [gr-qc].[44] A. Buonanno and T. Damour, Phys. Rev. D , 084006(1999).[45] T. Damour, A. Nagar, M. Hannam, S. Husa and B. Brug-mann (2008), arXiv:0803.3162[gr-qc].[46] E.K. Porter and B.S. Sathyaprakash, Phys. Rev. D ,024017 (2005).[47] N. Yunes and E. Berti (2008), arXiv:0803.1853[gr-qc].[48] E. Poisson, Phys. Rev. D , 1497 (1993).[49] C. Cutler, L. S. Finn, E. Poisson, and G. J. Sussman,Phys. Rev. D , 1511 (1993).[50] A. Buonanno, Y. Chen, Y. Pan, and M. Vallisneri, Phys.Rev. D , 104003 (2004).[51] T. Damour, Phys. Rev. D , 124013 (2001).[52] T. Damour, B. R. Iyer, P. Jaranowski, and B. S.Sathyaprakash, Phys. Rev. D , 064028 (2003).[53] T. Damour (2008), arXiv:gr-qc/0802.4047v1.[54] T. Damour and A. Nagar, Phys. Rev. D , 044003 (2007), arXiv:0704.3550 [gr-qc].[55] T. Damour, P. Jaranowski and G. Sch¨afer (2008),arXiv:0803.0915 [gr-qc].[56] T. Damour and A. Gopakumar, Phys. Rev. D , 124006(2006).[57] T. Damour, A. Nagar, E. N. Dorband, D. Pollney, andL. Rezzolla (2007), arXiv:0712.3003 [gr-qc].[58] B. R. Iyer and C. M. Will, Phys. Rev. Lett. , 113(1993).[59] B. R. Iyer and C. M. Will, Phys. Rev. D , 6882 (1995).[60] T. Damour, N. Deruelle, Phys. Rev. Lett. , 81(1981).[61] T. Damour, in Gravitational Radiation , edited byN. Deruelle and T. Piran (North-Holland, Amsterdam,1983), pp. 59–144.[62] G. Sch¨afer, Gen. Relativ. Gravit. , 255 (1986).[63] L. Blanchet, T. Damour, B.R. Iyer, C.M. Will and A.G.Wiseman, Phys. Rev. Lett. , 3515 (1995).[64] L. Blanchet, Phys. Rev. D , 714 (1997).[65] M. Boyle, A. Buonanno, L.E. Kidder, A.H. Mroue,Y. Pan, H.P. Pfeiffer and M.A. Scheel (2008),arXiv:0804.4184 [gr-qc].[66] L. Blanchet (2002), arXiv:0207037[gr-qc]. [67] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P.Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky,Phys. Rev. D , 124038 (2007), arXiv:0710.0158 [gr-qc].[68] E. Poisson, Phys. Rev. D , 5719 (1995).[69] T. Tanaka, H. Tagoshi, and M. Sasaki, Prog. Theor.Phys. , 1087 (1996).[70] C. Cutler, T. A. Apostolatos, L. Bildsten, L. S. Finn,E. E. Flanagan, D. Kennefick, D. M. Markovic, A. Ori,E. Poisson, G. J. Sussman and K. S. Thorne, Phys. Rev.Lett. , 2984 (1993).[71] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom,G. Lovelace, and M. A. Scheel, Class. Quantum Grav. , S59 (2007).[72] J. G. Baker, J. R. van Meter, S. T. McWilliams, J. Cen-trella, and B. J. Kelly, Phys. Rev. Lett. , 181101(2007).[73] M. Hannam, S. Husa, U. Sperhake, B. Br¨ugmann,and J. A. Gonzalez, Phys. Rev. D (2008),arXiv:0706.1305v2 [gr-qc].[74] A. Buonanno, Y. Chen, Y. Pan, and M. Vallisneri, Phys.Rev. D70