Gravitational Waves from a Particle in Circular Orbits around a Rotating Black Hole to the 11th Post-Newtonian Order

We compute the energy flux of the gravitational waves radiated by a particle of mass $\m$ in circular orbits around a rotating black hole of mass $M$ up to the 11th post-Newtonian order (11PN), i.e. $v^{22}$ beyond the leading Newtonian approximation where $v$ is the orbital velocity of the particle. By comparing the PN results for the energy flux with high precision numerical results in black hole perturbation theory, we find the region of validity in the PN approximation becomes larger with increasing PN orders. If one requires the relative error of the energy flux in the PN approximation to be less than $10^{-5}$, the energy flux at 11PN (4PN) can be used for $v\lessapprox 0.33$ ($v\lessapprox 0.13$). The region of validity can be further extended to $v\lessapprox 0.4$ if one applies a resummation method to the energy flux at 11PN. We then compare the orbital phase during two-year inspiral from the PN results with the high precision numerical results. We find that for late (early) inspirals when $q\le 0.3$ ($q\le 0.9$), where $q$ is the dimensionless spin parameter of the black hole, the difference in the phase is less than 1 ($10^{-4}$) rads and hence these inspirals may be detected in the data analysis for space detectors such as eLISA/NGO by the PN templates. We also compute the energy flux radiated into the event horizon for a particle in circular orbits around a non-rotating black hole at 22.5PN, i.e. $v^{45}$ beyond the leading Newtonian approximation, which is comparable to the PN order derived in our previous work for the energy flux to infinity at 22PN.

in Sec. 4. Section 5 is devoted to a summary and discussions. Finally in the appendices we describe supplemental equations to practically compute the formulas in Sec. 2. Throughout this paper, we use geometrized units with c = G = 1.

Basic formulation
We solve the Teukolsky equation to calculate gravitational radiation from a particle orbiting around a Kerr black hole. In the Teukolsky formalism, the gravitational perturbation of the Kerr black hole is described by a master variable. If we consider the outgoing radiation to infinity, the master variable, the Newman-Penrose quantity Ψ 4 , is related to the amplitude of the gravitational wave at infinity by The Teukolsky equation can be solved by the decomposition of Ψ 4 as where a is the angular momentum of the black hole and the angular function −2 S aω ℓm (θ) is the spin-weighted spheroidal harmonic with spin weight s = −2, normalized as π 0 | −2 S aω ℓm (θ)| 2 sin θ dθ = 1.
The decomposition of Ψ 4 Eq. (2) leads to the separation of the Teukolsky equation into radial and angular parts, where λ is the eigenvalue of −2 S aω ℓm , ∆ = r 2 − 2M r + a 2 , K = (r 2 + a 2 )ω − ma and T ℓmω is the source term of the particle.

3/30
Using the two independent solutions Eq. (6), with the Green function method one can construct a solution of the radial Teukolsky equation that is purely outgoing at infinity and purely incoming at the horizon where the Wronskian W ℓmω is given as Then, the solution has the asymptotic form at the horizon as while the solution has the following asymptotic form at infinity as Using the formula of the source term T ℓmω [2,3], Z ∞,H ℓmω are expressed as  , and where A nn0 and other terms are given in Appendix A. When the particle follows bound geodesics of Kerr spacetime, there exist three fundamental frequencies for the orbits [34] and hence the frequency spectrum of T ℓmω becomes discrete. In the case of circular orbits,Z ∞,H ℓmω in Eq. (12) takes the form where Ω = v r /(r 0 (1 + qv 3 r )) is the angular frequency of the particle, r 0 is the orbital radius, q = a/M , and v r = M/r 0 is the orbital velocity. 4/30 The time-averaged gravitational wave luminosity at infinity is then given by [35] where · · · represents the time average, ω = mΩ, and (dE/dt) N is the Newtonian quadrupole formula defined by with v ≡ (M Ω) 1/3 . Similarly, the time-averaged gravitational wave luminosity at the horizon becomes [35] dE where α ℓmω = 256(2M r + ) 5 k(k 2 + 4ǫ 2 )(k 2 + 16ǫ 2 )ω 3 |C| 2 , Finally, the gravitational waveforms are given in terms ofZ ∞ ℓmω as In this paper, using Eqs. (15), (17) and (19) we compute the gravitational energy flux and waveforms in the post-Newtonian approximation, i.e., in the expansion with respect to v = (M Ω) 1/3 . For this purpose, it is necessary to compute the asymptotic amplitudes Z ∞,H ℓmω , which involve calculations of the angular Teukolsky function −2 S aω ℓm (θ) and the radial Teukolsky functions R in,up ℓmω (r). To this end, in Appendices B and C, we give a short review of the calculations of the series expansions of −2 S aω ℓm (θ) and R in,up

11PN results for the time-averaged energy flux
In this paper, we derive the 11PN formula for the energy flux in the case of a test particle in a circular orbit around the equatorial plane of a Kerr black hole. Since the expressions are very long, we exhibit the 7.5PN expression for the energy flux at infinity in Sec. 3.1 and the next new 7PN terms in the energy flux into the horizon in Sec. 3.2. We also compute the energy flux into the event horizon for a particle in a circular orbit around a Schwarzschild black hole up to 22.5PN beyond the Newtonian approximation. The complete expressions for the energy flux will be publicly available online [ where κ = 1 − q 2 , γ is the Euler constant, ζ(n) is the zeta function, (20) are the new terms derived by the post-Newtonian approximation in this paper. 1 Among these terms, the O(v 9 ) − O(v 11 ) terms agree with the analytic expressions in Ref. [36], which determined the post-Newtonian coefficients of the energy flux up to 20PN by fitting with very accurate, one part in 10 600 , numerical calculation of the energy flux. For the 6PN and higher PN order energy flux at infinity, in Ref. [36], some of the post-Newtonian coefficients are not extracted as analytic values but as numerical values. This is not only because it is difficult to numerically extract analytic coefficients for combinations of transcendental numbers such as π, Euler's constant, and logarithms of prime numbers, but also because numerical fitting of post-Newtonian coefficients is done by presenting these coefficients as a polynomial in q although irrational functions in q such as polygamma functions and logarithms appear from 6.5PN onward as shown in Eq. (20). Further, by performing a small q expansion of our 11PN expression, we also find that our 11PN energy flux to infinity is consistent with the one in Ref. [36] up to 11PN. 2 From Eq. (20), we find the coefficient in q (ln v) 0 at 6PN is given by The above analytic value is consistent with the numerical value of the coefficient in q (ln v) 0 at 6PN energy flux to infinity in Ref. [36], 83.16039023577041 . . ..

Horizon flux
The next new 7PN terms for the energy flux into the horizon are given by dE dt where dE dt (9) H is defined through and N PN is the PN order, i.e. N PN = 11 when the PN order is 11PN. Notice that the energy flux down into the horizon starts at O(v 5 ) (O(v 8 )), i.e. 2.5PN (4PN), beyond the quadrupole formula when q = 0 (q = 0) [32,37]. Again, performing a small q expansion of our 11PN expression, we find our 11PN energy flux to the horizon is consistent with the one in Ref. [36] up to 11PN. For the case of a particle in a circular orbit around a Schwarzschild black hole, we also derive the energy flux down the event horizon at 22.5PN, which is consistent with the one in Ref. [36] up to 22.5PN.

Comparisons between 11PN results and numerical results
To investigate the accuracy of the energy flux in the post-Newtonian approximation, we compare PN results with numerical results, based on a method in Refs. [38,39]. With this 9/30 numerical method, one can investigate gravitational waves with an accuracy of about 14 significant figures in double precision calculations. Hence one can use the numerical results to estimate the accuracy in the PN results by a comparison. For the comparison in this section, we need to compute the energy flux using Eqs. (15) and (17). To numerically compute the energy flux, we set the maximum value of ℓ to 15, which gives the relative error in the energy flux better than 10 −5 for the comparisons in this section. For the energy flux at 11PN, we need to compute ℓ up to 13 (5) for the energy flux to infinity (the horizon).
In Sec. 4.1, comparisons for the the energy flux are done for several values of the spin of the Kerr black hole. In Secs. 4.2 and 4.3, the same comparisons are done using resummation techniques, the factorized resummation introduced in Ref. [40] and the exponential resummation in Ref. [41], for the post-Newtonian approximation to the energy flux. We will see how resummation methods improve the performance in the post-Newtonian approximation for the energy flux. Finally, in Sec 4.4, we compare the total cycle of orbits during a two-year inspiral for representative binaries in the eLISA frequency band. . From these figures, one will find that the relative error becomes smaller with increasing PN order for v ≤ 0.3, except for accidental agreements for certain values of the velocity. However, the relative error around ISCO does not necessarily become smaller with increasing PN order when q > 0.3. The relative error for 11PN is smaller than 10 −5 when v 0.33, irrespective of the values of q investigated in the paper. Figure 3 shows the relative error in the energy flux down the horizon from numerical results and PN approximations as a function of the orbital velocity up to ISCO in the case of the Schwarzschild black hole. The agreement between the numerical energy flux and post-Newtonian energy flux becomes better when the PN order is higher even around ISCO. The relative error in the 22.5PN energy flux into the horizon around ISCO is about 10 −5 , which is comparable to the one for the 22PN energy flux to infinity in Ref. [27].

Energy flux: Factorized resummation to PN approximation
In this section, we compare the total energy flux from numerical results with PN results using a factorized resummation. The factorized resummation was introduced to improve the convergence in the PN energy flux to infinity for a test particle moving in Schwarzschild spacetime [40,42,43] and Kerr spacetime [31]. The factorized resummation was then extended to the PN energy flux down the horizon of the Schwarzschild black hole in Ref. [44] and the Kerr black hole in Ref. [45].

Factorization of the energy flux at infinity.
In the factorized resummation of the energy flux at infinity, we decompose the multipolar gravitational waveforms into five factors as where ǫ p denotes the parity of the multipolar waveforms, h Absolute values of the relative error in the total energy flux from numerical results and PN approximations as a function of the orbital velocity, v = (M/r 0 ) 1/2 [1 + q (M/r 0 ) 3/2 ] −1/3 , up to ISCO for q = 0.1, 0.3, 0.5, 0.7, 0.9 and 0.998. The relative error for 11PN is smaller than 10 −5 when v 0.33. formalism, T ℓm resums the leading logarithms of the tail effects,δ ℓm is the supplemental phase factor and ρ ℓm is the ℓth root of the amplitude of the waveforms, which takes care of a term linear in ℓ at 1PN in the waveforms and means that a better convergence in the factorized waveforms might be expected (for more details see, e.g., Refs. [40,42,43]).
The first factor h (N,ǫp) ℓm takes the form  where φ is the orbital phase and n and c ℓ+ǫp (ν) are functions of the symmetric mass ratio ν ≡ µ M/(M + µ) 2 , defined by 12/30 The second factorŜ whereẼ andL z are the specific energy and the angular momentum of the particle, given bỹ The third factor T ℓm in Eq. (23) is defined by where r 0s = 2M/ √ e is introduced to reproduce the test-particle limit waveforms [31].
The fourth and fifth factors in Eq. (23),δ ℓm and ρ ℓm , can be derived by comparing the multipolar waveforms Eq. (23) with those obtained from Eq. (19). For the comparison, it is useful to express waveforms Eq. (19) in terms of the −2 spin-weighted spherical harmonics where −2 P lm (θ) is defined as where we used the orthogonality condition of the −2 spin-weighted spherical harmonics, andX is the complex conjugate of X. Note that in Eq. (32) the mixing ofZ ∞ ℓmω happens among the same m and different ℓ modes [31]. Note that the infinite summation over ℓ ′ in Eq. (32) can be truncated at a certain ℓ ′ for a given post-Newtonian order sinceZ ∞ ℓ ′ mω = O(v ℓ ′ +2+ǫp ) (see, e.g., Ref. [30]).

Factorization of the energy flux down the horizon.
For the factorized resummation of the energy flux down the horizon, the modal energy flux is decomposed as [44,45] where the factor (1 − 2v 3 r + /a) is motivated by the factor k = ω − ma/(2M r + ) = m (Ω − a/(2M r + )) in Eq. (18), which is responsible for the sign of the modal energy flux to the horizon. The second factor η N,H ℓm represents the leading term in the modal energy flux into the horizon and takes the form where n (H,0) ℓm = − 5 32 and c ℓ+ǫp (ν) where (z) n = Γ(z + n)/Γ(z).
The definition for the third factorŜ  Figures 4 and 5 show the relative error in the total energy flux from numerical results and PN approximations as a function of the orbital velocity up to ISCO using the factorized resummation to PN approximations [40]. From these figures, one will find that the relative error becomes smaller as PN order becomes higher for v ≤ 0.3, except for accidental agreements for certain values of the velocity. However, the relative error around ISCO does not necessarily become smaller for higher PN orders when q > 0.3. The relative error for 11PN is smaller than 10 −5 when v 0.4, irrespective of the values of q investigated in the paper. We note that the region of the velocity, v 0.4, is larger than the one using the Taylor expanded PN energy flux, v 0.33.

Energy flux: Exponential resummation to PN approximation
In this section, we compare the total energy flux from numerical results with PN results using the exponential resummation [41,46]. In the exponential resummation, the modal energy fluxes to infinity, η ∞ ℓm , and the horizon, η H ℓm , are decomposed as where η N,∞ ℓm and η N,H ℓm are the leading terms for η ∞ ℓm and η H ℓm respectively, the denominator (1 − 3v 2 r + 2qv 3 r ) is the square of the denominator ofŜ The factorsη ∞ ℓm andη H ℓm in the exponential resummation, Eq. (40), can be derived by comparing with the Taylor expanded modal energy fluxes η ∞ ℓm and η H ℓm . Figures 6 and 7 show the relative error in the total energy flux from numerical results and PN approximations as a function of the orbital velocity up to ISCO using exponential 15/30   resummation to PN approximations [41]. From these figures, one will find that the relative error becomes smaller as the PN order becomes higher for v ≤ 0.3, except for accidental agreements for certain values of the velocity. However, the relative error around ISCO does not necessarily become smaller at higher PN orders when q > 0.3. The relative error for 11PN is smaller than 10 −5 when v 0.4, irrespective of values of q investigated in the paper. Again, we note that the region of the velocity, v 0.4, is larger than the one using the Taylor expanded PN energy flux, v 0.33. 16

Phase difference during the two-year inspiral
We compare the orbital phase from PN results with numerical results during two-year inspirals to estimate the applicability of the PN results in the data analysis. For the comparison, we choose two representative systems of EMRIs in the eLISA frequency band, System-I and System-II, following Refs. [26,27,47,48]. System-I is an early inspiral of an EMRI with masses (M, µ) = (     For the calculation of the orbital phase, we define the phase as Ψ ℓm (t) = m t 0 Ω(t ′ )dt ′ , where Ω(t) = M 1/2 /r(t) 3/2 /(1 + qM 3/2 /r(t) 3/2 ) is the angular frequency of the particle and r(t) is the orbital radius as a function of time. The orbital radius r(t) is derived as r(t) = t (dr/dt ′ ) dt ′ = t (∂r/∂Ẽ) ( Fig. 2 but using exponential resummation to the energy flux in the post-Newtonian approximation. The relative error for 11PN is less than 10 −5 when v 0.4, whose region is larger than v 0.33 for the Taylor expanded energy flux in Fig. 2. time to perform the numerical integration is less than a second if we use the cubic spline interpolation. Figures 8 and 9 show absolute values of the difference in the orbital phase for the dominant ℓ = m = 2 mode between the PN and the numerical results during two-year inspirals for several values of the spin of the black hole. As for the PN approximations, we show results using the factorized resummation in Sec. 4.2, which are better than those using the Taylor expanded PN energy flux and comparable to those using the exponential resummation. The dephases between the 11PN results and numerical results after the two-year inspiral are less than 10 −4 rad for System-I. However, the dephases after the two-year inspiral become larger than a radian for System-II when q > 0.3. Thus, one has to derive higher PN order results 19/30 for the energy flux to achieve a dephase of less than a radian for System-II when q > 0.3, which represents a stronger-field situation than the one for System-I.

Summary
We have investigated gravitational waves from a particle moving in circular orbits in Kerr spacetime using the post-Newtonian approximation and computed the energy flux up to 11PN. We have also computed the energy flux down the event horizon for a particle in circular orbits around a Schwarzschild black hole at 22.5PN beyond the Newtonian approximation to fill the gap in the PN order between the energy flux at infinity, currently known at 22PN, and the event horizon, previously known at 6.5PN beyond Newtonian approximation.
To investigate how higher PN order expressions improve the applicability to data analysis of eLISA/NGO, comparisons between PN results and high-precision numerical results in black hole perturbation theory have been done. We first compared PN energy flux to numerical energy flux and found that the region of validity in the PN energy flux becomes larger as the PN order becomes higher. If the relative error of the energy flux in the PN approximation should be less than 10 −5 , the energy flux at 11PN satisfies this requirement for 20/30 v 0.33, which clearly shows an improvement from v 0.13 in an earlier work at 4PN [24]. The region of validity in the 11PN energy flux can become larger, v 0.4, if one uses resummation techniques such as factorized resummation [40] and exponential resummation [41]. 3 The region of validity might become further larger if one takes account of the structure of homogeneous solutions of the Teukolsky equation more carefully [46].
Finally, we compared the orbital phase during the two-year inspiral using the factorized resummed PN flux and the high-precision numerical flux. We found that the dephase is less than 1 (10 −4 ) rad for late (early) inspirals when q ≤ 0.3 (q ≤ 0.9). This implies that the 11PN factorized resummed flux may be used to detect early inspirals in the data analysis of eLISA/NGO. To detect gravitational waves from late inspirals when q > 0.3, however, it is necessary to obtain higher PN order expressions than 11PN. From numerical calculations in black hole perturbation theory, it is estimated that we may need to compute at least up to ℓ = 30, i.e. 28PN, to obtain the relative error of 10 −5 in the energy flux at ISCO for q = 0.9 [27]. If it is not possible to perform such a high PN order calculation, it may be necessary to use other approaches that compute unknown PN coefficients by numerical fitting [36,47,48].
(B16) 24/30 According to the theory of three-term recurrence relations [54], A (1) is a minimal solution and A (1) requires that the eigenvalue s E ℓm (ξ) satisfies a certain transcendental equation, which is expressed in terms of continued fractions.
In order to obtain the equation that determines the eigenvalue s E ℓm (ξ), it is convenient to introduce the following quantities: . (B17) Using the three-term recurrence relation Eq. (B11) we can express R n as an infinite continued fraction, and L n as a finite continued fraction, The expression for R n is valid if this infinite continued fraction converges. Noting the properties of the three-term recurrence relations (see p. 35 in Ref. [54]), it can be proved that the continued fraction Eq. (B18) converges if the eigenvalue s E ℓm (ξ) is finite. We obtain the equation to determine the eigenvalue s E ℓm (ξ) dividing Eq. (B11) by the expansion coefficients s A (n) ℓm (ξ) where R n+1 and L n−1 are defined by the continued fractions Eqs. (B18) and (B19), which are convergent for finite values of s E ℓm (ξ). There are many roots in Eq. (B20) for given n, m, s and ξ. These roots are associated with the same m, s, and ξ but with different ℓ. In order to find the root for a given ℓ, it is useful to choose n = n ℓ ≡ ℓ − (α + β)/2 in Eq. (B20) since in the limit ξ → 0 all the terms in Eq. (B20) become O(ξ 2 ). This means that the choice naturally gives the leading term of a series expansion of the eigenvalue s E ℓm (ξ) in terms of ξ as ℓ(ℓ + 1). When | ξ | is not large, we can obtain the analytic expression of s E ℓm (ξ) in a series of ξ as where For the numerical calculation to determine s E ℓm (ξ), one can use the analytic expression of s E ℓm (ξ) above as an initial value to find the root in Eq. (B20).
2 Γ(n + α + 1)Γ(n + β + 1) (2n + α + β + 1)Γ(n + 1)Γ(n + α + β + 1) ℓm (ξ) is made by the requirement that s S aω ℓm (x) reduces to the spin-weighted spherical harmonics in the limit ξ → 0. 26/30 C. Homogeneous solutions of the radial Teukolsky equation In this paper, we use a formalism developed by Mano, Suzuki, and Takasugi (MST) to compute the homogeneous solutions of the radial Teukolsky equation [28,29]. In the formalism, analytic expressions of homogeneous solutions are given using two kinds of series expansions in terms of hypergeometric functions and Coulomb wave functions, which are, respectively, convergent at the horizon and infinity. One can obtain analytic expressions of the asymptotic amplitudes of the homogeneous solutions by analytic matching of the two kinds of solutions in the overlapping region of convergence. Compared to numerical integration methods to solve the Teukolsky equation, the formalism is quite powerful for very accurate numerical calculations of gravitational waves [36,38,39,55]. The formalism is also very powerful for the performance of post-Newtonian expansions of gravitational waves at higher orders since the series expansion of homogeneous solutions is closely related to the low-frequency expansion. Applying the formalism to the post-Newtonian approximation in black hole perturbation theory, the energy flux going down the horizon was calculated up to 6.5PN for a particle in a circular and equatorial orbit around a Kerr black hole [32] and the 2.5PN energy flux to infinity was computed for a particle in a slightly eccentric and inclined orbit around the Kerr black hole [56,57]. More recently, we applied the formalism to obtain the 5.5PN waveforms for a particle in a circular orbit around a Schwarzschild black hole [30] and the 4PN waveforms for a particle in a circular and equatorial orbit around the Kerr black hole [31], which, respectively, confirmed the 5.5PN energy flux in Ref. [23] and the 4PN energy flux in Ref. [24]. In Refs. [26,27], we extended the formalism to obtain very high PN expressions for the energy flux to infinity for the particle in a circular orbit around the Schwarzschild black hole. For more details of the formalism, we refer the reader to a recent review, Ref. [2].
In the MST formalism, one can expand a homogeneous solution of the radial Teukolsky equation in a series of Coulomb wave functions as whereẑ = ω(r − r − ), τ = (ǫ − m q)/κ, (a) n = Γ(a + n)/Γ(a) and F N (η, z) is a Coulomb wave function defined by where Φ(α, β; z) is the confluent hypergeometric function, regular at z = 0 (see Sec. 13 in Ref. [58]). Note that the so-called renormalized angular momentum ν is introduced in the homogeneous solution in a series of Coulomb wave functions, Eq. (C1). ν is a generalization of ℓ, which has a property such that ν → ℓ as ǫ → 0, and is determined through conditions that the series of Coulomb wave functions, Eq. (C1), converges and actually represents a homogeneous solution of the radial Teukolsky equation.
The series of Coulomb wave functions, Eq. (C1), converges and represents a homogeneous solution of the radial Teukolsky equation if ν satisfies the following equation: where R ν n and L ν n are defined in terms of infinite continued fractions, which can be derived from the three-term recurrence relation, Eq. (C3). Observe that one can obtain two kinds of expansion coefficients, a ν n , from two kinds of the continued fractions, R ν n and L ν n . If ν is chosen to satisfy Eq. (C5) for a certain n, the two kinds of the expansion coefficients coincide and the series of Coulomb wave functions, Eq. (C1), converges for r > r + .