{'fraction_non_alphanumeric': 0.033328506612607396, 'fraction_numerical': 0.018896611642050393, 'mean_word_length': 4.303596569819532, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ineffectiveness of Pad\'e resummation techniques in post-Newtonian approximations

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 the gravitational wave energy flux do not accelerate the convergence of the standard Taylor approximants even in the test mass limit, and there is no reason why Pad\'e transformations should help in estimating parameters better in data analysis. Moreover, adding a pole to the flux seems unnecessary in the construction of these Pad\'e-approximated flux formulas. Pad\'e approximants may be useful in suggesting the form of fitting formulas. We compare a 15-orbit numerical waveform of the Caltech-Cornell group to the suggested Pad\'e waveforms of Damour et al. in the equal mass, nonspinning quasi-circular case. The comparison suggests that the Pad\'e waveforms do not agree better with the numerical waveform than the standard Taylor based waveforms. Based on this result, we design a simple EOB model by modifiying the ET EOB model of Buonanno et al., using the Taylor series 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 with the numerical waveform, much smaller than 0.04 radians the phase uncertainty in the numerical data itself. An EOB Hamiltonian can make use of a Pad\'e transformation in its construction, but this is the only place Pad\'e transformations seem useful.


I. INTRODUCTION
Even though general relativity was developed at the beginning of the twentieth century, no analytical solution is known for the two-body problem. Until recently, attempts to find a numerical solution failed because of the complexity of the mathematical equations and the instabilities inherent in the analytical formulations being used.
Studying the late dynamical evolution of these inspiraling compact binaries is important because they are among the most promising source of gravitational waves for the network of laser interferometric detectors such as LIGO and VIRGO. The detection of these gravitational waveforms (GW) is important for testing general relativity in the strong field limit. Moreover, these detectors can extract from the waves physical data about these sources such as the component masses and spins and the orbital eccentricity. For an unbiased extraction of these parameters, a large bank of accurate waveforms needs to be constructed. Numerical relativity alone cannot compute all the waveforms needed because of the computational cost. Instead, the waveforms are based on post-Newtonian (PN) approximations [14,15].
The post-Newtonian approximation is a slow-motion, weak-field approximation to general relativity. In order to produce a post-Newtonian waveform, the PN equations of motion of the binary are solved to yield explicit expressions for the accelerations of each body in terms of 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 gravitational waveform h and the gravitational wave flux F in terms of radiative multipole moments [26]. These radiative multipole moments are in turn related to the source multipole moments, which can be given in terms of the relative position and relative velocity of the binary [27]. Instead of comparing the post-Newtonian waveform with a numerical waveform along a certain direction with respect to the source, the comparison can be done in all directions by decomposing the waveform in terms of spherical harmonic modes. For an equal-mass non-spinning binary, the (2, 2) mode h 22 [28,29,30,31] is often used to compare numerical and post-Newtonian waveforms because it is the dominant mode. Its time derivativeḣ 22 is used to compute the gravitational wave flux. The resulting expressions for the orbital energy E, the gravitational energy flux F and the amplitude h 22 are given as Taylor series of the frequency-related parameter where M is the total mass of the binary and G = c = 1. The invariantly defined "velocity" another dimensionless parameter, is often used in writing these Taylor series.
Computing PN series to high order is difficult and time consuming. Since the various PN expressions are given as slowly convergent Taylor series, the Padé transformation [32,33] was suggested in Ref. [34] to accelerate the convergence of these series. The Padé transformation, P m n , consists of writing a Taylor series, T k , of order k as the ratio of two polynomials, one of order m in the numerator, and another of order n in the denominator, such  that m + n ≤ k. If well-behaved, this method accelerates the convergence of a Taylor series as the order of the Padé transformation, m + n, is increased. For example, in Table I we compare the convergence of the Taylor expansion of the exponential function, Exp n (v)(≡ e v ), at order n to its Padé approximant, Exp m+ǫ m (v) = P m+ǫ m [Exp n (v)], along the diagonal, where m = ⌊n/2⌋ and ǫ = 0 or 1. After twelve terms (n = 11), the last two partial sums of the Taylor expansion converge to 4 significant figures. However, the last two Padé approximants Exp 5 5 (v) and Exp 6 5 (v) converge to 6 significant figures. The error between the exact value of the exponential, 7.46331734, and the Padé approximant Exp 6 5 (v = 2.01) is 6 × 10 −8 , while the error between the eleventh order partial sum and the exact value is 10 −5 . Fig. 1 shows the convergence of the Taylor expansion of the exponential function and its Padé approximant.
The hope of accelerating the convergence of the post-Newtonian Taylor series of the energy and flux motivated the use of their Padé approximants to construct Padé approximant waveforms [34,35,36,37,38,39,40,41,42,43,44,45,46]. If these resummation techniques accelerate the convergence of the Taylor series in PN approximations, the range of validity of PN approximations suggested by Ref [47] could be extended. Moreover, the work of Refs. [48,49] in the test mass limit motivated the addition of a simple pole to the flux F of a binary system as the bodies approach the light ring orbit. By mathematical continuity, the existence of a pole in the equal mass case was anticipated [34].
More recently, waveforms are constructed by including these ideas in Effective One Body (EOB) models. The EOB approach [35, 36, 38, 39, 40, 41, 42, 43, 44, 50, 51,52,53,54,55,56,57] aims at providing an accurate analytical description of the motion and radiation of coalescing binary black holes. The approach consists of three separate ingredients: 1) a description of the conservative Hamiltonian part of the dynamicsĤ, 2) a formulation of the radiation reaction force F from the radiated flux F and 3) an expression of the GW waveform amplitude emitted by the coalescing binary system (i.e h 22 ).
The flux plays an important role in approximating the radiation 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. Since the equations of motion are known only to 3.5PN order, one has to rely on the assumed balance between energy loss in the system and radiated flux at infinity [63,64] to generate an approximate expression of the radiation reaction force at 3.5PN order beyond the leading term.
Ref. [65] computes the GW energy flux and GW frequency derivative from a highly accurate numerical simulation of an equal-mass, non-spinning black hole binary. By assuming energy balance, the (derivative of the) center-of-mass energy is estimated. These quantities are then compared to the numerical values using various Taylor, Padé and EOB models. The main goal of Ref. [65] is taking a set of well-established proposals in the literature for approximating waveforms and seeing how well they work in practice. Another goal of Ref. [65] is to examine some modifications of those proposals. The main goal of this paper, by contrast, is to show that a key ingredient in those proposals does not appear to be necessary.
In Ref. [66], Blanchet gave an argument that Padé and EOB resummations are unjustified because for two comparable-mass bodies there is no equivalent of the Schwarzschild light-ring orbit at the radius r = 3M . His argument is based on the PN coefficients of the binary's energy and their relation to predicting the innermost circular orbit. He finds that the radius of convergence of the PN series, which is related to the radius of the light-ring orbit, is around one (instead of 1/3 as for Schwarzschild). Blanchet concluded that Taylor series converge well for equal masses and that templates based on Padé/EOB are not justified because the dynamics of two bodies in General Relativity does not appear as a small "deformation" of the motion of a test particle in Schwarzschild. This paper arrives at similar conclusions but not by considering the innermost circular orbit, which is not precisely defined in the full nonlinear case. Instead, we compare Pade approximants of the flux and Pade/EOB waveforms to the numerical data of Refs. [67,68].
In this paper, we focus on testing two main techniques involved in building EOB models: the systematic use of Padé approximants, and the addition of a pole to the flux. The goal is to simplify these models by removing any unnecessary procedures in designing waveforms that provide good agreement with numerical waveforms.
Damour et al. [34,35] first suggested techniques for resumming the Taylor expansions of the energy and flux functions. Starting from the PN expansions of the energy E and the flux F , they proposed a new class of waveforms called P -approximants, based on three essential ingredients. The first step is the introduction of new energy-type (Eq. 4) and flux type (Eq. 16) functions, called e(v) and f (v) respectively. The second step is to Padé-approximate the Taylor expansion of these functions. The third step is to use these Padé transforms in the definition of the energy E (Eq. 6) and Padéapproximated flux (Eq. 20). The last step is to construct either the Padé-approximated waveform as in Sec. IV or the EOB waveform as in Sec. V. Schematically, the suggested procedure is summarized by the following map: Our notation is to denote by T m n (x) the Padé approximant of a k-th order Taylor series T k (x) with an m-th order polynomial in the numerator and an n-th order polynomial in the denominator such that m + n ≤ k, i.e. the Padé approximant of e k (x) is e m n (x). In Section II, we compare the 3PN Taylor series of the energy function to its possible Padé approximants using the intermediate energy function e(x), as suggested by Damour et al. [34]. We compute the last stable orbit frequency, defined as the frequency for which the energy reaches a minimum as a function of frequency, and also the poles of the energy in the complex plane corresponding to each possible Padé approximant. The large variation of last stable orbit frequency and poles does not  suggest good convergence of the Padé-approximated intermediate energy function e(x). The energy function E(x) is strongly dependent on the choice of the Padé approximant of e(x). Accordingly, the Padé waveform will also be strongly dependent on the choice of the Padé approximant.
In Section III, we present two possible methods for calculating the Padé approximant of the flux function. The first method simply takes the Padé approximant of the Taylor series treating the logarithmic contribution as constant. Following [34], the second method adds a pole to the Taylor series, factors out the logarithmic contribution to the series, and then computes the Padé approximant of the resulting Taylor series. We test the convergence of the Padé approximant for both methods versus their Taylor series. We find that the Padé approximants of the flux do not converge any faster than their Taylor counterpart.
A simple example that illustrates the problem is shown in Table II. There we compare the partial sums of the Taylor series for the flux with the corresponding Padé approximants in the test mass limit. 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. Even for a relatively small value of x, namely x = 0.04 (v = 0.2), the Taylor series is converging very slowly. After 12 terms, only about 4 or 5 significant digits seem reliable. Moreover, the Padé resummation shows very similar behavior; there is no improvement in the convergence. We will return to this example in Fig. 3.
In Section IV, we generate all the possible Padé wave-forms as suggested by Damour et al. [34] corresponding to 3 and 3.5 PN order. The waveform approximation requires the choice of a pole. We use the only physical pole, found from the 2PN Padé-approximated energy E 1 1 . We also use the last stable orbit from the 3PN energy Taylor series E 3 . The results are not very sensitive to this choice. We compare the Padé waveforms to a 15-orbit numerical waveform in the equal mass, nonspinning quasi-circular case [67]. The phase difference in these comparisons ranges between 0.05 and a few radians for well-defined Padé approximants (not having a pole in the frequency domain of interest) when the matching of the numerical and Padé waveforms is done at the gravitational wave frequency M ω = 0.1 [67]. None of the Padé waveforms agrees with the numerical waveform better than the Taylor series T4-3.5/3.0PN, which has an error of 0.02 radians. (We identify post-Newtonian approximants with three pieces of information: the label introduced by [35] for how the orbital phase is evolved; the PN order to which the orbital phase is computed; and the PN order at which the amplitude of the waveform is computed. See Ref. [67] for more details.) Our conclusion is that the Padé approximant might be helpful in suggesting fitting formulas but it does not provide a more rapidly convergent method. Note that the Padé transform also fails to accelerate the convergence of the T2, T3 and h 22 Taylor series (see Refs. [35,67] for the definition of these Taylor series).
In Section V, based on the results of the previous sections, we design a simple EOB model (closely related to the ET EOB model of Ref. [40]) using the Taylor series of the flux. We add one unknown 4PN term that we fit for by maximizing the agreement between the EOB model waveform and the numerical waveform. The model does not require adding a pole to the flux, nor an a priori knowledge of the last stable orbit from the energy function. This simple EOB model, with only one parameter to fit for, agrees with the numerical waveform to within 0.002 radians (3 × 10 −4 cycles). (This is six times smaller than the claimed numerical accuracy of [39], smaller by an even larger factor than the claimed numerical accuracy of [45], and 25 times smaller than the gravitational wave phase uncertainty of the numerical waveform. See Table III in Ref [67] for more details.) This model agrees with the numerical waveform better than any previously suggested Taylor, Padé or EOB waveform.

II. ENERGY FUNCTION
Damour et al. [34] introduced a new energy-type function e(x), where x is the PN frequency related parameter. This assumed more "basic" energy function e(x) is constructed out of the total relativistic energy E tot (x) of the binary system. Explicitly where m 1 , m 2 are the masses of the bodies. The total relativistic energy function E tot is related to the post Newtonian energy function E(x) through where M is the total mass (M = m 1 + m 2 ). Solving for E(x) in terms of e(x) using Eqs. (4) and (5), we get [34] E(x) = 1 + 2ν where the symmetric mass ratio is ν = m 1 m 2 /M 2 . The orbital energy function E(x) is known as a Taylor series E k up to 3PN order as a function of x and ν [15] Using the above equations, we compute the Taylor series expansion, e k (x), of e(x) up to 3PN order: In the test mass limit (ν → 0), the exact function e(x) coincides with the Padé approximant P 1 1 (x) of its Taylor expansion in Eq. (8) This quantity has a pole at x pole = 1/3. The orbital energy is then and it derivative is The last stable orbit occurs where so in the limit ν → 0 the last stable orbit is at exactly x lso = 1/6. On the grounds of mathematical continuity between the test mass limit ν → 0 and the finite mass ratio case, Damour et al. [34] argued that the exact function e(x) should be meromorphically extendable in at least part of the complex plane and should have a simple pole on the real axis. They suggested that Padé approximants would be excellent tools for giving accurate representations of functions having such poles.
Once we know the Taylor series of the new energy function, e k (x), we compute its Padé approximant e m n (x), with m + n ≤ k. The Padé-approximated energy, E m n (x) is obtained by replacing e(x) in Eq. (6) with e m n (x). In the equal mass case (ν = 1/4), we can define several Padé approximants of e k (x). The most interesting Padé approximants have a maximal sum of their indices, since they should be closest to the unknown exact function if the Padé resummation is converging. In Fig. 2, we show a plot of the PN energy function E 3PN (x) and its Padé approximants Although the Padé approximants of the energy are of maximal order, they differ significantly. Good convergence of the Padé approximants requires good agreement between approximants of the same order n + m, if there is no pole in the region of interest (0 < x 0.4). For example, there is no a priori reason why one should prefer either E 2 1 or E 1 2 . Although both have the same order and are equally close to the diagonal, the difference between these functions is quite large.
In Table III we compute the locations of the poles and the last stable orbits for all of these Padé approximants. The ill-convergence of the Padé transform is again seen by looking at the variation of the last stable orbit positions. In Table III for example, x lso of E 2 1 differs by about 8% from x lso of E 1 2 . Moreover, for finite ν, the poles are all complex or not in the interval [0, 1] except for the case x pole = 52/109, corresponding to the Padé-approximated energy E 1 1 . There is no reason why this should be the "exact" pole that should be used in the formalism, since none of the third-order Padé approximants of the 3PN energy has a physical pole.
In summary, using Padé approximants for the energy function in the equal mass case does not seem to provide any benefit. The differences between the various Padé approximants of the energy are large. The quantities x pole and x lso do not show any regular behavior that could be a sign of a physical pole that could be found by using the Padé transform.

III. FLUX FUNCTION
The general form of the PN flux at order N is: where the normalized flux F is a Taylor expasion in v with logarithmic terms where the post-Newtonian coefficients A i and B i are functions of the mass ratio parameter ν. They are given in the test mass limit in Ref. [69] and in the equal mass quasi-circular case in Ref. [15]. The flux series has a logarithmic contribution starting at 3PN. Padé approximants, however, are well defined only for pure polynomials. Two possible methods are therefore used to compute the Padé approximant of the flux. The first method simply treats the logarithmic terms as constants and resums the series as a pure polynomial such that the Padé- The second method, suggested by Ref. [34], defines a new flux function f by adding a pole, factoring the logarithmic terms from the series, and finally computing the Padé approximant of the pure polynomial. Since we would like to check the convergence of the Padéapproximated flux versus its Taylor series, we sketch the definitions of the various functions involved. According to Ref. [34], two ideas are needed for a good representation of the analytic structure of the flux. First, since in the test mass limit F is thought to have a simple pole at the light ring [49], one might expect it by continuity to have a pole in the comparable mass case. This motivates the introduction of the following factored flux function, f (v; ν):  TABLE III: Values of the poles and last stable orbit (lso) of the energy for the case ν = 1/4. The poles x pole and last stable orbit frequency of the function E m n (x) depend significantly on which Padé approximant is constructed from the Taylor series e k (x). The only physical pole is x pole = 52/109, which is at a larger value than the pole in the test mass limit. The position of the last stable orbit also varies significantly.
where v pole is the pole of the Padé-approximanted energy function used.
Second, the logarithmic term that appears in the flux function needs to be factored out so we can use the standard Padé transformation. After factoring the logarithmic terms out, the flux function f becomes 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éapproximated energy. Then the Taylor series of the flux with a pole is defined as The Padé approximant of the intermediate flux function where v lso (e m n ; ν) denotes the last stable orbit velocity for the Padé approximant P m n e(x) . Finally, the corresponding Padé approximant of the flux F (v) is given by where v pole (e m n ; ν) denotes the pole velocity defined by e m n (x).

A. Flux for the test mass case
The exact gravitational wave luminosity F is not known analytically in the test particle limit. It has been computed numerically by Poisson [68]. The post-Newtonian expansion of the flux is known in the test mass limit to 5.5PN order [69]. This allows us to test the rate of convergence of the Taylor series of the normalized flux F n (Eq. 14) and its Padé-approximant F m n constructed treating the logarithmic term as a constant (Eq. 14). We also test the convergence of the flux function F n (Eq. 18) and its Padé approximant F m n (Eq. 20). These convergence tests use the known values v pole = 1/ √ 3 and v lso = 1/ √ 6 for the test mass limit as discussed in Sec. II.
In Fig. 3, we test the convergence of the various flux functions at the velocity value 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. We use the Padé approximant along the diagonal P m+ǫ m where ǫ = 0 or 1. The rates of convergence of the Taylor expansion and its Padé approximant are nearly equal for the two methods, whether or not we include a pole. As the PN order increases, the Taylor series and its Padé approximant alternate in which provides 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 3 PN order the Padé approximant of the flux F m n (Eq. 20) fits the numerical data the best, while at 3.5 and 5PN order the Taylor series of the flux (Eq. 14) is the best. At 5.5PN the Padé approximant of the flux (Eq. 20) gives the best agreement. The results are similar for other values of v. No method has the best convergence rate.
According to Padé theory, the convergence of the Padé approximant is best along the diagonal, but it is equally good along the off-diagonal terms if no pole exists in the region of interest (i.e. no zeroes appear in the denominator of the Padé approximant.) For this reason, we show the error between all the possible maximal Padéapproximated fluxes F   This suggests that the Padé approximation should only be used to suggest a fitting formula for the numerical data, since there is no internal self-consistency in the agreement. The offdiagonal approximants do not show any regular pattern of convergence to the numerical data nor are they better than the Taylor series. , Fn and F m+ǫ m are given in Eqs. 14, 15, 18 and 20 respectively. The Padé approximants do not converge faster than their Taylor series counterparts. The Padé and Taylor series alternate at providing the best agreement with the "exact" data as the PN order increases. Contrast the behavior here with Fig. 1.

B. Flux for the equal mass case
For binaries of comparable mass on a quasi-circular orbit, the flux is known only to 3.5PN order [15]. In Ref. [65] for a quasi-circular non-spinning binary, the numerical flux was computed by integrating the spinweighted spherical harmonic components of the Weyl scalar Ψ 4 . The numerical flux data we use in this paper was provided by Harald P. Pfeiffer and Michael Boyle. The estimated error in measuring the flux data was about 0.2%. The velocity range for the simulation was from v ∼ 0.26 (x ∼ 0.06) to v ∼ 0.4 (x ∼ 0.16).
In the equal mass case, we cannot do an accurate convergence test early in the evolution as in Fig. 3 for two reasons. The first reason is the "junk radiation" (noise early in the evolution from imprecise initial data) during the first few orbits. The second reason is the inability to accurately define the numerical flux as a function of the orbital frequency of the binary. The numerical normalized flux is computed as a function of ω 22 /2, where ω 22 is the wave frequency of theḣ 22 mode. Instead, in Table IV we compare the convergence of the four flux functions F n , F converged to one additional significant figure over F n , however F m+ǫ m cannot reliably be considered more accurate than F n because it converges to a slightly different value. The Padé approximants do not seem to converge to a larger number of significant figures than the Taylor flux function F n .
In Fig. 5, we plot the numerical normalized flux F N R , the 3.5PN flux F 3.5 and the maximal Padé-approximated flux functions F 3 4 , F 4 3 , F 5 2 , F 6 1 and F 7 0 (≡ F 7 ). Although a good fit to the numerical flux during the last 15-orbit inspiral before merger, there is no guarantee that this is true at low velocities.

IV. PADÉ WAVEFORMS
The construction of the post-Newtonian waveforms requires solving the post-Newtonian equations describing the motion of the binary and the generation of gravitational waves. Substituting the orbital evolution predicted by the equations of motion into the expressions for the waveform would not generate waveforms accurate enough for matched filtering in detecting gravitational waves [70]. To compute the waveform at 3PN order, it is necessary to solve the equations of motion at 5.5PN order, because the radiation reaction contributes to the equations of motion starting at 2.5PN order. However, for a non-spinning binary of equal mass and on a circular orbit, accurate waveforms at 3PN order can be constructed under two further assumptions. The first assumption is that the binary follows a slow adiabatic inspiral. The second assumption is that of energy balance between the orbital binding energy and the energy emitted by the gravitational waves, where the energy balance equation is defined as The procedure of constructing the standard Padé waveforms [34] is similar to one used to construct the TaylorT1 waveforms in Ref. [34,67]. The main difference is the use of Padé approximants of the energy and flux to compute the orbital phase, as described in Secs. II and III, instead of their Taylor expansions. The orbital phase used in the Padé waveforms is obtained by numerically integrating The fraction on the right side of Eq. (22) is retained as a ratio of the Padé approximants of the post-Newtonian expansions, and is not expanded further before numerical integration. The waveform is produced by substituting the orbital phase into the spherical harmonic mode h 22 of the post-Newtonian waveform, which is known up to 3PN order [28,29,30,31].
Given the expressions for the Padé-approximated energy and flux in Sections II and III, and the Taylor series of the waveform amplitude [28,29,30,31], there is still a set of choices that must be made in order to produce a Padé-approximated waveform that can be compared to our numerical waveform. These include 1. The Padé approximant of the orbital energy, E k l .
2. The flux function and its Padé approximant F n m .
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 waveform amplitude are kept.

A. Procedure
We consider numerical gravitational waves extracted with the Newman-Penrose scalar Ψ 4 , using the same procedure as in [71]. To minimize gauge effects, we compare its (2, 2) component extrapolated to infinite extraction radius according to Ref. [67]. The extracted waveform is split into real phase φ and real amplitude A, defined by Ref. [67] as The gravitational-wave frequency is given by The spherical harmonic component (2,2) of Ψ 4 is then compared to the numerically twice-differentiated post-Newtonian expression of h 22 , A 22 , as in Ref. [67]. Following [67,72,73], the matching procedure needed to set the arbitrary time offset t 0 and the arbitrary phase offset φ 0 is done by demanding that the PN and NR gravitational wave phase and gravitational wave frequency agree at some fiducial frequency ω M .

B. Results
In this section, we compare the numerical waveform to the Padé waveforms corresponding to the 3.5 PN order of energy and flux using the 3PN Taylor series of the post-Newtonian amplitude A 22 . The energy and flux functions used are those suggested by Ref. [34]. We do not generate all possible waveforms using different Padé approximants of the energy or the flux at low PN orders, since all these resummed series showed no improvement in the convergence rate.
As introduced in Sec. II, we use the Padéapproximated energy E 3 0 , E 2 1 , E 1 2 and E 0 3 corresponding to the PN Taylor series of the energy, and the Padéapproximated energy, E 1 1 , corresponding to its 2PN Taylor expansion. For the flux, the diagonal Padé approximant F 3 3 is used in addition to all possible Padé approximants of flux at 3.5PN order, F m 7−m , where 0 ≤ m ≤ 7 as described in Sec. III.
The Padé-approximated flux has two parameters, v lso and v pole as discussed in Sec. III. The value v pole = 52/109 is used. We also tested varying the pole location, but found that we could not improve the agreement significantly.
From Table III, any value of the velocity of the last stable orbit could be used. We use the 3PN value v lso = 0.254 and also use v lso = 0.199. The latter is used when the Padé approximant E 1 1 is employed in the construction of the waveform. In the remaining cases, we use v lso = 0.254 since it is quite close to the estimates from other Padé approximants of the energy. The effect of changing the value of v lso is not significant compared to changing the order of the Padé approximant for the energy or the flux.
To do the comparison, we match the Padéapproximated and numerical waveforms at the wave frequency M ω = 0.1. Then we measure the maximum phase difference between the numerical waveform and each of these Padé waveforms during the inspiral when the numerical wave frequency is between M ω = 0.035 and M ω = 0.1 (as in the upper panel of Fig. 7). Our results are summarized in Fig. 6, which shows the phase differences for each of the Padé approximants of energy E k l and flux F m 7−m . On the same Figure, we include phase differences for the waveforms generated using the Padéapproximated flux F 3 3 under the m = −1 entry. When E 1 1 is used, the phase error ranges between 2 and 5 radians as m increases from −1 to 7. Using all the possible Padé approximants of the 3PN energy, the estimated phase difference ranges from 0.05 to 2.5 radians. Using the Taylor series with a pole (m = 7) resulted in a large phase difference ranging between 1 and 1.5 radians. The diagonal Padé term F 3 4 of the flux generates similar phase differences, ranging from 0.06 to 0.2 radians as the Padé order of the energy changes.
The Padé-approximated waveforms do not fit the numerical data better than the waveforms using the Taylor expansion of the flux. Although the Padé-waveforms along the diagonal have a phase difference less than 0.25 radians, none of these waveforms fits the numerical waveforms better than TaylorT4 at 3.5PN order as shown in Ref. [67]. Moreover, the dependence of the phase difference on the Padé order suggests that there is no reason why it should help in estimating the parameters better in data analysis. This is as expected from the poor convergence of the the Padé approximant of the flux discussed in Sec. III.
The Padé resummation techniques were also tested on the Taylor series for the amplitude, and they showed no improvement in the convergence of the series. In adition, none of the tests that were performed on the Padé resummed Taylor series of the T2 and T3 waveforms showed a faster convergence rate. In fact, there is no improvement in convergence for any Taylor series in the PN approximation that we have investigated.

V. SIMPLE EOB MODEL
We have described the failure of the Padé resummation techniques to accelerate the convergence of any PN Taylor series, the absence of any signature of a pole in the flux in the equal mass case, and the erratic pattern of agreement between the Padé waveforms and the numerical waveform. It seems one might as well simply use the 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 good agreement with the numerical waveform by using a simple EOB model. The only parameter we introduce and fit for is an unknown 4PN contribution to the flux.

A. EOB waveforms
The EOB formalism [53] is a non-perturbative analytic approach that handles the relative dynamics of two relativistic bodies. This approach of resumming the PN theory is expected to extend the validity of the PN results into the strong-field limit. The procedure for generating an EOB waveform follows closely the steps in Sec. IV. Instead of using the energy balance equation, we compute the orbital phase by numerically integrating Hamilton's equations. The EOB waveform is generated by substituting the orbital phase into the waveform amplitude A 22 at 3PN order. The two fundamental ingredients that allow computing the orbital phase are the real HamiltonianĤ and the the radiation reaction F φ .

B. Hamilton's equations
In terms of the canonical position variables r and φ and their conjugate canonical momenta p r and p φ , where r is the relative separation and φ is the orbital phase, the real dynamical Hamiltonian is defined as [54]: where and where the radial potential A function is defined as the series The Taylor series of the A function is replaced by its Padé approximant A 1 3 . Here the Padé approximant is not used to 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, the EOB Hamiltonian is non-physical for the last few orbits; the orbital frequency stays nearly constant for several orbits before merger. For the B function, the Taylor expansion suffices: Then Hamilton's equations of motion are given in the quasi-circular case by where F φ is the radiation reaction in the φ direction representing the nonconservative part of the dynamics. In Eq. 32, ∂ φĤ = 0 sinceĤ is independent of φ. The radiation reaction is deduced from the post-Newtonian flux as in Refs. [41,58,59] In this equation, we have introduced an unknown 4PN flux term, F 8 , the only parameter that we fit for in this EOB model.

C. Initial conditions
To integrate Hamilton's equations, we need appropriate initial conditions for a quasi-circular orbit. Refs. [41,42,52] indicate how to define some "post-adiabatic" initial conditions. However, these initial conditions do not generate an orbit with as low an eccentricity as the numerical simulation, roughly 5×10 −5 . At a given radius r, starting from the post-adiabatic initial conditions of p r and p φ , we therefore reduce the eccentricity iteratively in two steps. The first step includes evolving Hamilton's equations in the conservative regime (F = 0) and iteratively changing the value of p φ until the eccentricity measured from the evolution of the orbital separation is of the order 10 −9 . The second step is based on evolving the nonconservative Hamilton's equations with the 4PN flux and iteratively changing the p r momentum until the eccentricity is again of the order 10 −5 . This circularization procedure is repeated as we iterate F 8 to maximize the agreement between the waveforms.

D. Best Fit of F8
To find the best fit for F 8 , we iteratively solve for the minimum in the phase difference between the numerical Phase and amplitude differences between the EOB waveform and the numerical waveform. After fitting for the best value of F8, the phase difference is less than 0.002. The early noise is due to junk radiation at the early stage of the numerical simulation. r * is the tortoise coordinate defined in [67].
and EOB waveforms. The waveforms are matched as in Sec. IV at the wave frequency mω = 0.1, and the phase difference is defined as the maximal phase difference during the inspiral phase up to the wave frequency mω = 0.1. We find a best fit value F 8 = −333.75 corresponding to the initial conditions r = 17, φ = 0, p r = −0.0008, p φ = 4.53235. A change of 1% in F 8 changes the maximal phase difference from less than 0.002 radians to about 0.01 radians. Note that without adding the fitting parameter F 8 , the phase difference is about 1.7 radians during the 15-orbit inspiral.

E. Results
In the upper panel of Fig. 7, we plot the phase difference between the numerical waveform and the EOB waveform computed using the 3PN Taylor series of the amplitude A 22 . The phase difference is less than 0.002 radians after maximizing the agreement between the waveforms in the region where mω ≤ 0.1. The early noise is due to junk radiation at the early stage of the numerical simulation as described in Sec II C of Ref. [67]. The phase uncertainty in the simulation was estimated to be 0.05 radians; See Table III in [67].
In the lower panel of Fig. 7, we plot the relative difference between the amplitude of the numerical waveform and the EOB waveform. The EOB waveform amplitude does not fit the numerical waveform amplitude as well as the wave phase does. This is expected because the waveform amplitude is known to 3PN order only, and no free parameter in the amplitude was fitted for. The agreement between the amplitude of this EOB model and the numerical waveform is similar to the agreement between the amplitude of TaylorT4 3.5/3.0 and the numerical waveform in Fig. 21 in [67].
This EOB model is a modification of the ET EOB model of Ref. [74]. It fits the numerical phase very well without using the Padé resummation techniques nor a pole in the flux.
Even though we have found very good agreement between the waveforms, these results only suggest that the EOB model is a very good fitting model. Moreover, having fit a particular waveform, there is no guarantee the model will have predictive power for a more general case.

VI. CONCLUSIONS
Convergence tests show that none of the Taylor series in the PN approximation, such as the energy or the flux, could be replaced by a Padé approximant that converges faster. Other attempts we tried to accelerate the convergence of these series also failed, as for example using the Levin method to accelerate convergence [33]. As a result, more reliable waveforms could not be constructed using a Padé resummation scheme. Moreover, the Padé waveforms also do not fit numerical simulation data better than the Taylor waveforms. Thus, they do not seem to be better than the Taylor waveforms in building templates for waveforms. This conclusion is independent of the Padé approximants used to test the convergence. Taking for example the sub-diagonal Padé approximant does not show any improvement in the convergence rate. In addition, this conclusion is independent of the numerical data we used. We can simply take the highest PN order of the Taylor series or the Padé approximant and use it as the "exact" value of the function to test the convergence at low frequency.
Based on the dependence of the flux on the velocity in the equal mass case, we do not find it helpful to add a pole to the flux. Therefore, we recommend using Taylor series instead of the Padé approximant to generate waveforms both in the time and frequency domains. The simple EOB model used in this paper agrees with the numerical data very well; the phase difference during the inspiral is much less then the estimated phase uncertainty in the numerical data. This model does not use Padé approximants or poles except in one place to enforce a last stable orbit. Since Padé approximation does not accelerate the convergence of any PN Taylor series, there is no reason why it should estimate parameters better in data analysis of waveforms. and Nicolás Yunes. We thank Jihad Touma for helpful discussions about Padé approximants, Harald P. Pfeiffer and Michael Boyle for providing the numerical data of the flux in the equal mass case, and Eric Poisson for providing the numerical data of the flux in the test mass limit. This work was supported in part by grants from the Sherman Fairchild Foundation to Cornell; by NSF grants PHY-0652952, DMS-0553677, PHY-0652929, and NASA grant NNG05GG51G at Cornell.