Calculation of radiation reaction effect on orbital parameters in Kerr spacetime

We calculate the secular changes of the orbital parameters of a point particle orbiting a Kerr black hole, due to the gravitational radiation reaction. For this purpose, we use the post-Newtonian (PN) approximation in the first order black hole perturbation theory, with the expansion with respect to the orbital eccentricity. In this work, the calculation is done up to the fourth post-Newtonian (4PN) order and to the sixth order of the eccentricity, including the effect of the absorption of gravitational waves by the black hole. We confirm that, in the Kerr case, the effect of the absorption appears at the 2.5PN order beyond the leading order in the secular change of the particle's energy and may induce a superradiance, as known previously for circular orbits. In addition, we find that the superradiance may be suppressed when the orbital plane inclines with respect to the equatorial plane of the central black hole. We also investigate the accuracy of the 4PN formulae by comparing to numerical results. If we require that the relative errors in the 4PN formulae are less than $10^{-5}$, the parameter region to satisfy the condition will be $p\gtrsim 50$ for $e=0.1$, $p\gtrsim 80$ for $e=0.4$, and $p\gtrsim 120$ for $e=0.7$ almost irrespective of the inclination angle nor the spin of the black hole, where $p$ and $e$ are the semi-latus rectum and the eccentricity of the orbit. The region can further be extended using an exponential resummation method to $p\gtrsim 40$ for $e=0.1$, $p\gtrsim 60$ for $e=0.4$, and $p\gtrsim 100$ for $e=0.7$. Although we still need the higher order calculations of the PN approximation and the expansion with respect to the orbital eccentricity to apply for data analysis of gravitational waves, the results in this paper would be an important improvement from the previous work at the 2.5PN order, especially for large $p$ region.


Introduction
The gravitational two-body problem is a fundamental issue in general relativity. This also attracts great interest in gravitational wave physics because binary inspirals are promising sources of gravitational waves which are expected to be detected directly by ongoing gravitational wave observatories in the world. Understanding the dynamics of binary system is required to predict the emitted gravitational waveforms accurately for efficient searches of the signal in observed data.
One of major approaches for this purpose is the gravitational self-force (GSF) picture in the black hole perturbation theory. In this picture, a binary is regarded as a point mass orbiting a black hole and the dynamics can be described by the equation of motion of the mass including the effect of the interaction with the self-field, that is, the GSF. After the formal expression of the GSF was presented by Mino, Sasaki and Tanaka [1] and Quinn and Wald [2], a lot of efforts have been devoted to develop practical formulations and methods to calculate the GSF (for example, refer to [3] for the formulation of GSF, [4,5] for the recent progress in practical calculations of GSF).
Although a lot of progress has been made, however, it is still challenging to calculate the GSF directly for general orbits, especially in Kerr spacetime. Practical calculations of the GSF with high accuracy will require a huge amount of time and computer resources mainly because of the regularization problem induced by the point mass limit. Therefore it is important to develop a way to reduce the cost of computing the GSF. The two-timescale expansion method [6] gives a hint for it: assuming that a point mass does not encounter any transient resonances (e.g. shown in [7]), the orbital phase, which is the most important information to predict the waveform, can be expressed in the expansion with respect to the mass ratio, η, as where Φ (0) and Φ (1) are quantities of order unity. The leading term, Φ (0) , can be calculated from the knowledge up to the time-averaged dissipative piece of the first order GSF, corresponding to the secular growth. The calculation of this secular contribution can be simplified significantly by using the radiative field defined as half the retarded solution minus half the advanced solution for the equation of the gravitational perturbation [8][9][10], i.e. the adiabatic approximation method, because the radiative field is the homogeneous solution free from the divergence induced by the point mass limit. This method allows us to calculate the leading term accurately without spending huge computational resources. On the other hand, the calculation of Φ (1) requires the rest of the first order GSF (the oscillatory part of the dissipative GSF and the conservative GSF) and the time-averaged dissipative piece of the second order GSF. There is no simplification in calculating these post-1 adiabatic pieces at present. Since Φ (1) is subleading, however, the requirement of the accuracy is not so high compared to that of the leading term. This fact suggests that it is possible to reduce the computational cost by using a suitable method with an appropriate error tolerance to calculate each piece of the GSF (for example, a hybrid approach is proposed in [11]). In this work, we focus on the time-averaged dissipative part of the first order GSF, which has the dominant contribution to the evolution of inspirals, and present the analytic post-Newtonian (PN) formulae. So far, several works in this direction had been done for two restricted classes of orbits: circular orbits and equatorial orbits. (See [12] and references therein for early works in 1990's). Recently, thanks to the progress of computer technology, the very higher order post-Newtonian calculations can be possible for circular equatorial orbits: the 22PN calculation of the energy flux is demonstrated in Schwarzschild case [13], and the 11PN calculation in Kerr case [14]. There is also the calculation of the secular GSF effects for slightly eccentric and slightly inclined (non-equatorial) orbits [10], and later it had been extended to orbits with arbitrary inclination [15], where the PN formulae of the secular GSF effects are presented in the expansion with respect to the orbital eccentricity. However, the calculation in [15] has been done only up to the 2.5PN order with the second order correction of the eccentricity. Also the absorption to the black hole is ignored there. The main purpose of this work is to update the results in [15] up to the 4PN order and the sixth order correction of the eccentricity, including the effect of the absorption to the black hole. 2/32 This paper is organized as follows. In Sec.2, we give a brief review of the geodesic motion of a point particle in Kerr spacetime, the gravitational perturbations induced by the particle, and the adiabatic approximation method of calculating the secular effect of the GSF. In Sec.3.1, we present the PN formulae of the secular changes of the the energy, azimuthal angular momentum, and Carter parameter of the particle due to the gravitational radiation reaction in the expansion with respect to the orbital eccentricity. In Sec.3.2, we investigate the accuracy of our PN formulae by comparing to numerical results given by the method in [16][17][18], which can give each modal flux at the accuracy about 14 significant figures. In Sec.3.3, we implement a resummation method to the PN formulae given in this work in order to improve the accuracy. In Sec.3.4, we discuss the convergence of the analytic formulae as the PN expansion and the expansion with respect to the eccentricity. Finally, we summarize the paper in Sec.4. For the readability of the main text, we present the PN formulae for the orbital parameters, the fundamental frequencies, the orbital motion in Appendices A and B, which are used in calculating the secular changes of the orbital parameters. We also present the PN formulae for the secular changes of an alternative set of the orbital parameters in Appendix C. Throughout this paper we use metric signature (− + ++) and "geometrized" units with c = G = 1.

Review of formulation: Adiabatic radiation reaction
The orbital evolution of a point particle due to the time-averaged dissipative part of the GSF is often described in terms of the secular changes of the orbital parameters. In order to calculate the changes, we need the information on the first order gravitational perturbations induced by the particle when it moves along the background geodesics. In this section, we review the geodesic dynamics of a point particle in Kerr spacetime, the gravitational perturbations induced by the particle, and the adiabatic evolution of the orbital parameters.

Geodesic motion
The Kerr metric in the Boyer-Lindquist coordinates, (t, r, θ, ϕ), is given by where Σ = r 2 + a 2 cos 2 θ, ∆ = r 2 − 2M r + a 2 , M and aM are the mass and angular momentum of the black hole, respectively. There are two Killing vectors related to the stationarity and axisymmetry of Kerr spacetime, which are expressed as ξ µ (t) = (1, 0, 0, 0) and ξ µ (ϕ) = (0, 0, 0, 1). In addition, it is known that Kerr spacetime possesses a Killing tensor, K µν = 2Σl (µ n ν) + r 2 g µν , where l µ and n µ are the Kinnersley's null vectors given by For the geodesic motion of a particle in Kerr geometry, there are three constants of motion related to the symmetries: where u α is the four velocity of the particle.Ê andL correspond to the specific energy and azimuthal angular momentum of the particle respectively.Q is called as the Carter constant, which corresponds to the square of the specific total angular momentum in Schwarzschild case. These specific variables are measured in units of the particle's mass, µ. One can recover the expressions in the standard units as There is another definition of the Carter constant, C ≡ Q − (aE − L) 2 , which vanishes for equatorial orbits. In this paper, we use C as one of the orbital parameters, instead of Q. By using these constants of motion, the geodesic equations can be expressed in the following form as where we introduced a new parameter λ through the relation dλ = dτ /Σ, and some functions as V tr (r) := r 2 + a 2 ∆ P (r), V tθ (θ) := −a(aÊ sin 2 θ −L), A generic geodesic orbit in Kerr spacetime can be characterized by three parameters, {E, L, C} 1 . In the case of a bound orbit, we can use an alternative set of parameters, {r p , r a , θ min }, instead of {E, L, C}, where r p and r a are the values of r at the periapsis and apoapsis and θ min is the minimal value of θ, respectively. Using this set of parameters, we can describe the range in which the motion takes place as r p ≤ r ≤ r a and θ min ≤ θ ≤ π − θ min . There is another useful choice of parameters used in [19], {p, e, ι}, defined by By analogy to the parametrization used in celestial mechanics, p, e, ι are referred as semi-latus rectum, orbital eccentricity, orbital inclination angle, respectively. For later convenience, we also introduce Y = cos ι and v = 1/p. Since v corresponds to the magnitude of the orbital velocity, it can be used as the post-Newtonian parameter. For example, we call the O(v 8 )-correction from the leading order as the fourth order post-Newtonian (4PN) correction.
It is worth noting that, by introducing λ, the radial and longitudinal equations of motion in Eq. (6), are completely decoupled. For an bound orbit, therefore, the radial and longitudinal motions are periodic with the periods, {Λ r , Λ θ }, defined by This means that these motions can be expressed in terms of Fourier series as where Ω r and Ω θ are the radial and longitudinal frequencies given by and we choose the initial values so that r(λ = 0) = r p and θ(λ = 0) = π/2. 2 Since the temporal and azimuthal equations of motion in Eq. (7) are divided into the rand θ-dependent parts, the solutions can be divided into three parts: the linear term with respect to λ, the oscillatory part with period of Λ r , and the oscillatory part with period of Λ θ . They can be expressed as where the index A runs over {r, θ}, and with · · · λ ≡ lim T →∞ (2T ) −1 T −T dλ · · · , representing the time average along the geodesic. We choose the initial conditions as t(λ = 0) = ϕ(λ = 0) = 0. Ω ϕ corresponds to the frequency of the orbital rotation.
In Appendices A and B, we present the PN formulae of the orbital parameters, {E, L, C}, the fundamental frequencies, {Ω t , Ω r , Ω θ , Ω φ }, and the Fourier coefficients of the motions in Eqs. (15), (16), (18) and (19). 2 If the ratio of the radial and longitudinal frequencies is irrational, we can adjust the origin of λ approximately so that the radial and longitudinal oscillations reach the minima simultaneously at λ = 0 [8]. On the other hand, it is not the case if the ratio is rational, i.e. the resonance case. This implies that the secular evolution of a resonant orbit cannot be described only by the PN formulae derived in this work [20].

Secular evolution of orbital parameters
The gravitational perturbations in Kerr spacetime can be described by the Weyl scalar, Ψ 4 , which satisfies the Teukolsky equation [21]. To solve the Teukolsky equation, the method of separation of variables is often used, in which Ψ 4 is decomposed in the form as where S Λ (θ) is the spin-2 spheroidal harmonics and Λ represents a set of indices in the Fourier-harmonic expansion, {ℓ, m, ω}. The separated equation for the radial function is given by where T Λ is the source term constructed from the energy-momentum tensor of the point particle, andλ is the eigenvalue determined by the equation for S Λ (To find the basic formulae for the Teukolsky formalism used in this paper, refer to the section 2 in [22] for example). The amplitudes of the partial waves at the horizon and at infinity are defined by the asymptotic forms of the solution of the radial equation as with r + ≡ M + √ M 2 − a 2 and k = ω − ma/(2M r + ). Since the spectrum with respect to ω gets discrete in the case of a bound orbit, Z H,∞ Λ take the form whereΛ denotes the set of indices, {ℓ, m, n r , n θ }, and ω mnrnθ ≡ Ω −1 t (mΩ ϕ + n r Ω r + n θ Ω θ ) .
With these amplitudes, the secular changes of the orbital parameters, {E, L, C}, can be expressed by where and C S is the Starobinsky constant given by [23] |C S | 2 = (λ + 2) 2 + 4aωm − 4a 2 ω 2 λ 2 + 36aωm − 36a 2 ω 2 6/32 It should be noted that, in these formulae, the averaged rates of change are expressed with respect to the Boyer-Lindquist time, which can be related to those with respect to λ [24] as for a function of time, I(t). Also it should be noted that each formula in Eqs. (26)- (28) can be divided into the infinity part and the horizon part: the former consists of the terms including the amplitudes of the partial waves at the infinity,Z ∞ Λ , the latter consists of the terms including the amplitudes at the horizon,Z H Λ . As for the energy and azimuthal angular momentum, the infinity parts are balanced with the corresponding fluxes radiated to infinity and the horizon parts with the absorption of the gravitational waves into the central black hole [23,25].
The practical calculation ofZ H,∞ Λ involves solving the geodesic equations, calculating two independent homogeneous solutions of Eq. (22) and the spin-2 spheroidal harmonics, and the Fourier transformation of functions consisting of them. In this work, we followed the same procedure proposed in [15] to perform these calculations analytically.
In performing the summation in Eqs. (26)-(28) practically, we need to truncate the summation to finite ranges ofΛ = {ℓ, m, n r , n θ }. To obtain the accuracy of the 4PN and O(e 6 ), it is necessary to sum ℓ in the range 2 ≤ ℓ ≤ 6 (2 ≤ ℓ ≤ 3), n r in the range −3 ≤ n r ≤ 3 (−2 ≤ n r ≤ 3) and n θ in the range −8 ≤ n θ ≤ 12 (−4 ≤ n θ ≤ 6) for the infinity (horizon) part. The other modes out of these ranges are the higher PN corrections than the 4PN order or the higher order corrections than O(e 6 ).

PN formulae of the secular changes of orbital parameters
In this work, we derived the analytic 4PN order formulae of Eqs.(26)- (28) in the expansion with respect to the orbital eccentricity, e, up to O(e 6 ) (we simply call them as the 4PN O(e 6 ) formulae). Since the full expressions of the 4PN O(e 6 ) formulae are too lengthy to show in the text, we show the infinity parts up to the 3PN order and the horizon parts up to the 3.5PN order (while we keep the expansions with respect to e up to O(e 6 )). The complete expressions of the 4PN O(e 6 ) formulae will be publicly available online [26].
The infinity parts of Eqs.(26)- (28) where the leading contributions are given by The horizon parts of Eqs.(26)-(28) are given by 16 + 120 e 2 + 90 e 4 + 5 e 6 8 + 9 From the leading order expressions in Eq. (35), one will find that the Carter parameter, C, does not change due to the radiation of the gravitational waves when Y = 1 (equatorial orbits) because (dC/dt) N = 0.
In the Schwarzschild case, the Carter parameter corresponds to the square of the equatorial angular momentum (the normal component to the rotational axis of the central black hole). Then there is expected to exist the duality between L 2 and C due to the spherical symmetry. In fact, from Eqs. (33) and (34), (and also from (37) and (38)), one can find that dL 2 /dt t vanishes in Y = 0 (polar orbits) while dC/dt t for Y = 0 coincides with dL 2 /dt t for Y = 1. This can be also realized by seeing that the secular change of the total angular momentum, d(L 2 + C)/dt t , is independent of Y . Then, it might be possible to understand that dL/dt ∞ t becomes 1.5PN from the leading order when q = 0 and Y = 0 (polar orbits) due to the spin-orbit coupling.
From the expressions of the horizon parts shown in Eqs. (36)- (38), we find that the absorption of the gravitational waves to the central black hole contributes at O(v 5 ) from the leading order in Eq. (35) t can be positive for q > 0, which means that the particle can gain the energy through a superradiance phenomenon. These observations are consistent with the results for circular, equatorial orbits shown in Refs. [27,28].
We also find that the superradiance terms in Eq. (36) vanish for Y = 0, and that dE/dt H t has only the 4PN and higher order corrections. The superradiance terms may come from the coupling between the black hole spin and the orbital angular momentum, like ∝ L · S ∝ q cos ι. Hence, when the orbital inclination increases (Y gets small correspondingly), the superradiance is suppressed [29].

Comparison to numerical results
To investigate the accuracy of the 4PN O(e 6 ) formulae derived in this work, we compare them to the corresponding numerical results given by the method established in Ref. [16][17][18], which enables one to compute the modal fluxes with the relative error of ∼ 10 −14 in double precision computations. In the practical computations, as well as in deriving the analytic expressions, we need to truncate the summation to finite ranges ofΛ = {ℓ, m, n r , n θ } in Eqs. (26)- (28). In order to save the computation time in the numerical calculation, we sum ℓ up to 7. We can check that the error due to neglecting terms for ℓ ≥ 8 is smaller than the relative error in the 4PN O(e 6 ) formulae from the corresponding numerical results up to ℓ = 7. We also truncate n r and n θ to achieve the relative error of ∼ 10 −7 in numerical results up to ℓ = 7. For the parameters investigated in the comparison, the relative error of ∼ 10 −7 achieved by truncating n r and n θ is again smaller than the relative error in the 4PN O(e 6 ) formulae from the numerical results up to ℓ = 7. Thus, we can regard numerical results as benchmarks to investigate the accuracy in our analytic formulae.

11/32
Here we define the relative error in the analytic formula of dE/dt t by where dE/dt Ana t denotes the analytic formula in order to distinguish it from the corresponding numerical result, dE/dt Num t . We also define the relative errors in the analytic formulae of dL/dt t and dC/dt t in a similar manner and denote them as ∆ L and ∆ C respectively. Fig. 1 shows several plots of ∆ E for the 4PN O(e 6 ) formula as a function of p for several sets of (e, ι) with q = 0.9. In the plots, we also show the relative errors in the 2.5PN O(e 2 ) and 3PN O(e 4 ) formulae for reference. From the plots for e = 0.1 (three on the top), one can find that ∆ E falls off faster than p −4 for p 10 (Similarly, the relative errors in the 2.5PN O(e 2 ) and 3PN O(e 4 ) formulae fall off faster than p −5/2 and p −3 ). Noting v = 1/p, this would be a good indication that our PN formula has been derived correctly up to required order. ∆ E is expected to contain not only higher order corrections than the 4PN order, but also the higher order corrections of eccentricity than O(e 6 ) in the lower PN terms, which will become dominant when p and e get larger. In fact, seeing the plots for e = 0.7 in Fig. 1, one can find that the relative error strays out of the expected power law line for large p. This behavior is clearer in the plots of the relative error in the 2.5PN O(e 2 ) formula. From Eqs. (32) and (36), we know that the relative error in the 2.5PN O(e 2 ) formula contains the O(e 4 ) correction in the O(v 0 ) term. The effect of this correction appears as large-p plateaus in the plots (also see Fig. 6). This may motivate us to perform the higher order expansion with respect to the orbital eccentricity in the PN formulae or to derive the PN formulae without performing the expansion with respect to the orbital eccentricity [30][31][32][33]. In addition, it might be noted that the behavior of the relative error does not strongly depend on the inclination angle ι for fixed q and e.
In Fig. 2, we show the relative errors in the 4PN O(e 6 ) formulae for the secular changes of the three orbital parameters, {E, L, C}, for several sets of (q, e) and ι = 50 • . As in the case of ∆ E shown in Fig. 1, the relative errors, ∆ L and ∆ C , fall off faster than p −4 when p 10, except for the large p region (p 100) in the case of e = 0.7. Thus, the 4PN O(e 6 ) formulae for the secular changes of the orbital parameters are expected to be valid up to O(v 8 ). From Fig. 2, one might think that it is enough to investigate only ∆ E to discuss the accuracy of our formulae since there are not large differences in the relative errors, ∆ E , ∆ L and ∆ C . Fig. 3 shows contour plots for ∆ E as a function of p and e for several sets of (ι, q). From these plots, one may be able to comprehend the accuracy of our PN formulae more easily than using Figs. 1 and 2. One will find that the relative error becomes smaller (larger) for larger (smaller) p and smaller (larger) e. Moreover, it might be noticed that the relative error does not strongly depend on the inclination angle ι for fixed q as expected from Fig. 1. If one requires ∆ E < 10 −5 as an error tolerance, one can use the contour line with the label 10 −5 to find the region of validity in the figure. For example, one will find that ∆ E < 10 −5 for p 50 and e = 0.1, p 80 and e = 0.4, and p 120 and e = 0.7 when q = 0.9. 12/32 We truncated the plots at p = 6 because the relative errors get too large (nearly or more than unity) in p < 6 to be meaningful. One finds that the relative error becomes smaller with increasing orders of the PN approximation and the expansion with respect to the eccentricity. The relative error in our 4PN O(e 6 ) formula falls off faster than p −4 when the eccentricity is small, e.g. e 0.4. Since v = 1/p, this would imply that our 4PN formula is correctly representing the secular change up to the 4PN order. Note, however, that the relative error in the 4PN O(e 6 ) formula for e = 0.7 falls off slower than p −4 when the semi-latus rectum becomes larger,e.g. p > 100. This might be because of the higher order corrections of e than O(e 6 ), which will contain the lower PN terms than the 4PN order. We also note that changing the inclination angle, ι, does not change the dependence on p of the relative error for fixed q and e so much. This might be checked more easily in contour plots in Fig. 3, which show the relative error as a function of p and e for fixed q and ι.

Implementation of an exponential resummation method
In order to improve the accuracy in the analytic PN formulae, one may apply some resummation methods such as Padé approximation [34], the factorized resummation [35][36][37] and the exponential resummation [38]. Since the exponential resummation may be the simplest 13/32 Edt dLzdt dCdt Fig. 2 The relative errors in the analytic PN formulae for the secular changes of the three orbital parameters, {E, L, C}, as functions of the semi-latus rectum p for ι = 50 • , q = 0.9, 0.5, 0.1 and −0.9 (from top to bottom) and e = 0.1, 0.4 and 0.7 (from left to right). We truncated the plots at p = max{6, p s (e, ι)}, where p s (e, ι) is the value of p at the "separatrix" (the boundary between stable and unstable orbits), because the relative errors get too large in p < 6 to be meaningful and the orbit is not stable for p < p s (e, ι). As pointed out in Fig. 1, the relative errors in our 4PN O(e 6 ) formulae fall off faster than p −4 when the eccentricity is small, e.g. e 0.4, while the fall-off gets slower when p is larger for e = 0.7. There are not large differences in the behaviors of ∆ E , ∆ L and ∆ C . This suggests that it might be enough to focus only on dE/dt t to investigate the accuracy and convergence of our 4PN formulae. one to implement among them, we here choose to implement the exponential resummation. We apply it to our 4PN formulae and check how the accuracy is improved.
To introduce the exponential resummation, we make use of the following identity The relative error in the 4PN O(e 6 ) formula for the secular change of the particle's energy, ∆ E , as a function of the semi-latus rectum p and the eccentricity e for q = 0.9, 0.5, 0.1 and −0.9 (from top to bottom) and ι = 20 • , 50 • and 80 • (from left to right). We truncated the plots at p = max{6, p s (e, ι)} because the relative errors get too large in p < 6 to be meaningful and the orbit is not stable for p < p s (e, ι). From the figures, it is easily found that the relative error becomes smaller (larger) for larger (smaller) p and smaller (larger) e with fixed q and ι. If one requires the relative error to be less than 10 −5 , the region in the semi-latus rectum p and the eccentricity e will be p 50 and e = 0.1, p 80 and e = 0.4, and p 120 and e = 0.7 when q = 0.9. It might be noticed that the relative error does not strongly depend on the inclination angle ι for fixed q as pointed out in Fig. 1.
where I = {E, L, C}. The exponential resummation can be obtained by replacing the exponent in (40) to the expansion with respect to v, 15/32 where we do not perform the expansion with respect to e. Since our PN formulae for dI/dt t are given at the 4PN order, we truncate F I n after O(v 8 ). Finally, the exponential resummed form is expressed as dI dt Fig. 4 shows the relative errors in the exponential resummed forms of the secular changes of E, L and C, estimated by using Eq. (39). We also show the relative errors in the Taylortype formulae in the same graphs for comparison. One will find that the relative errors in the exponential resummed forms are less than those in the Taylor-type formulae in most cases, except for dC/dt t in the case of q = 0.9, (e, ι) = (0.1, 50 • ). Using the exponential resummation when q = 0.9 and ι = 50 • , the region to satisfy ∆ E < 10 −5 is extended to p 40 from p 50 for e = 0.1, p 60 from p 80 for e = 0.4, and p 100 from p 120 for e = 0.7. This might motivate us to use the resummation method to improve the accuracy of Taylor-type formulae even in the case of general orbits.

Convergence with respect to v and e of the analytic formulae
Apart from comparisons to the numerical results, we may also discuss the convergence property in our PN formulae with respect to v and e by investigating the contribution of each order of v and e in the formulae although this is a rough estimation.
Since ∆ n shows the relative importance of the O(v n ) term in the PN formulae, it can be used to investigate the convergence property with respect to v: it is expected that |∆ n+1 | < |∆ n | for moderately large n if the PN formula converges. In Fig. 5, we plot the relative contribution of each order, ∆ n , as a function of p for several sets of (e, ι) and q = 0.9. From this figure, one may find that ∆ n does not strongly depend on the inclination angle, ι, as shown in Sec. 3.2, while it strongly depends on e. The convergence gets worse when the orbital eccentricity becomes larger. This tendency is particularly evident in the small-p region. Fixing the value of p, the orbit with larger e passes by closer to the central black hole and will be affected by the stronger gravitational field. Hence the PN convergence is expected to be worse when the eccentricity becomes larger.
Next, in order to investigate the convergence of the expansion with respect to the orbital eccentricity in the PN formula, we introduce A n as where the term A 0 coincides with the energy flux for circular orbits and A n = 0 when n is odd. One may ask whether the condition, |A 2n+2 e 2n+2 | < |A 2n e 2n |, is satisfied for moderately large integer n if the series converges. From Fig. 6, it is found that the condition is satisfied 16/32  . We truncated the plots at p = 6 because the relative errors in the PN formulae get too large in p < 6 to be meaningful. Using the exponential resummation, the accuracy is improved in most cases. For example, the region to satisfy ∆ E < 10 −5 is improved from p 50 to p 40 for e = 0.1, p 80 to p 60 for e = 0.4, and p 120 to p 100 for e = 0.7. This would suggest us to try to apply resummation methods to the PN formulae even in the case of general orbits.
in most cases. As expected, the convergence becomes slower when the eccentricity is larger. Especially, the convergence gets worse when p 10 in e = 0.7 case. The calculation of the higher PN corrections will be necessary to improve the bad convergence for small p. We also note that A n does not strongly depend on ι for a fixed q as in Sec. 3.2.

Summary
We have derived the secular changes of the orbital parameters, the energy, azimuthal angular momentum, and Carter parameter of a point particle orbiting a Kerr black hole, by using the post-Newtonian approximation in the first order black hole perturbation theory. We have extended the previous work [15], which derived formulae up to the 2.5PN order with the second order correction with respect to the eccentricity, to the 4PN order with the sixth order correction with respect to the eccentricity. We have also included the contribution due to the black hole absorption, which has not been included in [15]. As shown in the case of equatorial, circular orbits [27,28], we have found that the secular changes of the 17/32 when q = 0.9. We truncated the plots at p = 6 because the relative contributions get too large in p < 6 to be meaningful. It is expected that |∆ n+1 | < |∆ n | for moderately large n if the PN formula converges. As shown in Sec. 3.2, ∆ n does not strongly depend on ι for a fixed e although ∆ n strongly depends on e. In fact, the convergence seems worse the orbital eccentricity becomes larger. This tendency is clear for small p, e.g. p 10.
three orbital parameters due to the absorption, appear at the 2.5PN (4PN) from the leading order in the Kerr (Schwarzschild) case, and that the 2.5PN and 3.5PN contributions of the absorption to the secular change of the particle's energy can be positive for q > 0, which implies that a superradiance can be realized in the Kerr case. We have also found that the superradiant contributions in the secular change of the energy get smaller when the inclination angle becomes larger and they vanishes for polar (Y = 0) orbits. This means that the superradiant scattering may be suppressed for inclined orbits [29].
To investigate the accuracy in our 4PN formulae, we have compared the formulae to highprecision numerical results [18] in Sec. 3.2. We have found that the accuracy gets worse when the orbital velocity and the orbital eccentricity become larger, as expected. If the relative error in the 4PN O(e 6 ) formula for the secular change of the energy is required to be less than 10 −5 , the parameter region to satisfy it might be p 50 for e = 0.1, p 80 for e = 0.4, and p 120 for e = 0.7 when q = 0.9. The region does not strongly depend on the orbital inclination angle. From Fig. 1, one can clearly find the improvement of the accuracy in our 18/32 . We plot the absolute value of A n e n as a function of the semi-latus rectum p for ι = 20 • , 50 • and 80 • (from left to right) when q = 0.9. We truncated the plots at p = 6 because the relative contributions get too large in p < 6 to be meaningful. It is expected that |A 2n+2 e 2n+2 | < |A 2n e 2n | for moderately large n if the series with respect to e converges. This condition is satisfied in most cases shown in this figure. The convergence becomes slower when the eccentricity is larger. Especially, the convergence for p 10 is quite bad in e = 0.7 case. We also note that A n does not strongly depend on ι for a fixed q as mentioned in Sec. 3.2.
PN formulae from the previous work at the 2.5PN order and the second order correction in the orbital eccentricity [15] whose relative error is larger than 10 −2 for p 100 when e 0.4 since, in this region, the error due to the truncation of the expansion with respect to the orbital eccentricity is larger than the one of the PN expansion. One may improve the accuracy of our PN formulae by using resummation methods. In this paper, we have applied the exponential resummation [38] to our 4PN formulae and confirmed that the resummation method improves the accuracy in most cases investigated here. For example, we found that the region in which the relative errors are less than 10 −5 can be extended from p 50 to p 40 for e = 0.1, p 80 to p 60 for e = 0.4, and p 120 to p 100 for e = 0.7.
We also investigate the convergence properties of the PN expansion and the expansion with respect to the orbital eccentricity, respectively. Both convergences get worse when the semi-latus rectum is smaller; in other words, the gravitational field becomes stronger. This tendency gets clearer in the case of large eccentricity, in which the particle passes by closer to the central black hole.
In order to improve the accuracy and convergence of the 4PN O(e 6 ) formulae near the central black hole and to obtain the physical information of the source in the strong-field region, it is necessary to derive the higher order corrections of the PN expansion and the expansion with respect to the eccentricity. It may be possible to avoid the expansion with respect to the eccentricity and to derive the PN formulae applicable to arbitrary eccentricity. So far the PN formulae of the rate of the energy loss without performing the expansion with respect to the eccentricity had been derived for equatorial orbits in [30][31][32][33]. The extension of these results to the case of inclined orbits is challenging: we can obtain analytic expressions for general bound geodesic orbits in Kerr spacetime without performing the expansion with respect to the eccentricity nor the inclination by using results in Ref. [39], while we need to reformulate the source term of the Teukolsky equation and the derivation of the partial waves constructed form the source term. We would like to leave it to the future work. computations were performed at the cluster "Baltasar-Sete-Sóis" in CENTRA/IST. Some analytic calculations were carried out on HA8000/RS440 at Yukawa Institute for Theoretical Physics in Kyoto University.

A. PN formulae for the orbital parameters and fundamental frequencies
In this section, we present the PN formulae of the orbital parameters, {Ê,L,Ĉ}, and the fundamental frequencies, {Ω t , Ω r , Ω θ , Ω φ }. Here we show the formulae up to the 3PN O(e 6 ) order to save space although it is possible to calculate them to the higher order. The higher order results will be publicly available online [26].

B. Fourier coefficients of bound orbits
Here we show the PN formulae of the Fourier coefficients in Eqs. (15), (16), (18) and (19) up to the 3PN O(e 6 ) order. The 4PN O(e 6 ) results obtained in this work will be available online [26].
In this work, we follow the same procedure as in [15] to derive the amplitudes of the partial waves, Z H,∞ Λ in (23). In the formal expression, the dependence of ϕ (θ) appears in the form of the combination as X ≡ sin θe iϕ (θ) , which can be expressed in the Fourier series as Therefore we show the Fourier coefficients of X instead of ϕ (θ) .

Fig. C1
The relative errors in the analytic formulae for the secular change of the orbital eccentricity due to the gravitational waves to infinity. We plot ∆ e , defined in a similar manner to Eq. (39), as a function of the semi-latus rectum p for q = 0.9, e = 0.1, 0.4 and 0.7 (from left to right) and ι = 50 • . We truncated the plots at p = 6 because the relative errors get too large in p < 6 to be meaningful. The relative error in the previous 2.5PN O(e 2 ) formula given in [15] strays off the p −3 line earlier than the 2.5PN O(e 2 ) formula in this paper. This trend is clearer for larger e. The relative errors in the 3PN O(e 4 ) and 4PN O(e 4 ) formulae fall off faster than p −3 and p −4 for small e cases as expected, while this is not the case for e = 0.7 because of the higher order correction of e than O(e 4 ).
However, the O(1/e) contribution turns out to vanish due to a cancellation in taking the combination, and hence de/dt t is O(e), which corresponds to the well-known fact that circular orbits remain circular [40,41], i.e. de/dt t = 0 when e = 0. This cancellation reduces the reliable order in de/dt t by O(e 2 ), compared to the order of dI/dt t for I = {E, L, C}.
Since we calculate dI/dt t up to O(e 6 ) in this paper, we can obtain de/dt t correctly up to O(e 4 ) from the leading order. dv/dt ∞ t and dY /dt ∞ t in Eqs. (C2) and (C4) are consistent up to the 2.5PN O(e 2 ) order with the previous results in Ref. [15], while we find inconsistency in the O(e 2 ) terms of the formula for de/dt ∞ t in [15]. This may be explained by the reduction in the reliable order mentioned above: the calculations of dI/dt t in [15] are done up to O(e 2 ), and therefore the resultant formula of de/dt t is reliable only at the leading order. We can also confirm it numerically. In Fig. C1 we show the relative errors in the two analytic formulae by comparing to numerical results [18] in a similar manner to Eq. (39). It can be found that the relative error in the previous 2.5PN O(e 2 ) formula strays out of the expected power law, p −3 , earlier than that in our 2.5PN O(e 2 ) formula. This trend is clearer for larger eccentricity. We can also confirm the validity of our formula by seeing that the relative errors in our 4PN O(e 4 ) formula falls off faster than p −4 . (The leading PN order of the difference in the O(e 2 ) terms between the previous and our formulae is 39 38 e 2 v 2 . If our formula contains any error in the e 2 v 2 term, the relative error will not fall off faster than p −4 for large p.) A similar reduction in the PN order occurs in the calculation of dY /dt ∞ t : although each term in the linear combination of Eq. (C1) for J = Y is O(v 8 ), the terms at the first two orders, O(v 8 ) and O(v 10 ), vanish due to a cancellation in taking the combination. As a result, the leading order of dY /dt t is O(v 11 ) and hence the reliable order relative to the leading term is reduced to O(v 5 ) (2.5PN order) when we have dI/dt t for I = {E, L, C} up to O(v 8 ) (4PN order).  Fig. C2 The relative errors in the 4PN formulae for the secular changes of the orbital parameters, {v, e, Y }. We plot the relative errors, ∆ J defined in a similar manner to Eq. (39), as functions of the semi-latus rectum p for q = 0.9, e = 0.1, 0.4 and 0.7 (from left to right) and ι = 50 • . We truncated the plots at p = 6 because the relative errors get too large in p < 6 to be meaningful. ∆ v and ∆ e fall off faster than p −4 , while ∆ Y approximately fall off as p −3 , slower than O(p −4 ). This confirms that the relative order of the PN correction of the analytic formula for dY /dt t is reduced from 4PN to 2.5PN because of the cancellation of the low PN terms.
In Fig. C2, we show the relative errors in the analytic PN formulae for the secular changes of the orbital parameters, {v, e, Y }, derived from the 4PN O(e 6 ) formulae of dI/dt t for I = {E, L, C}. Similarly in Fig. 1, the relative errors in the analytic formulae for dv/dt t and de/dt t as functions of the semi-latus rectum p fall off faster than p −4 when the eccentricity is small. Observe, however, that the relative error in the analytic formula for dY /dt t falls off faster than p −5/2 , but slower than p −4 .
From the leading order expressions in Eq. (C5), one will find the well known fact that equatorial orbits stay in the equatorial plane [10,15,42], i.e. dY /dt t = 0 when Y = 1. In the Schwarzschild case (q = 0), the secular changes of v and e do not depend on Y in addition to dY /dt t = 0. This implies that the orbital plane can be fixed on the equatorial plane (θ = π/2) due to the spherical symmetry of Schwarzschild spacetime.
One will also find that the radiation reaction reduces the orbital eccentricity and increases the orbital velocity since (de/dt) N ≤ 0 and (dv/dt) N ≥ 0 [40,41], while the radiation reaction increases (decreases) the inclination angle since (dY /dt) N ≤ 0 ((dY /dt) N ≥ 0) when q ≥ 0 (q ≤ 0) [10,15,42]. Moreover, the secular change of the inclination angle is smaller than those of the other orbital parameters since d ln e/d ln v t = O(v 0 ) and d ln Y /d ln v t = O(v 3 ).