# Relativistic orbits around spinning supermassive black holes. Secular evolution to 4.5 post-Newtonian order

Clifford M. Will, Matthew Maitra

► To cite this version:

Clifford M. Will, Matthew Maitra. Relativistic orbits around spinning supermassive black holes. Secular evolution to 4.5 post-Newtonian order. Physical Review D, 2017, 95 (6), pp.064003. 10.1103/Phys-RevD.95.064003. hal-01903081

HAL Id: hal-01903081

https://hal.science/hal-01903081v1

Submitted on 5 Oct 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Relativistic orbits around spinning supermassive black holes: Secular evolution to 4.5 post-Newtonian order

Clifford M. Will<sup>1,2,*</sup> and Matthew Maitra<sup>3,†</sup>

<sup>1</sup>Department of Physics, University of Florida, Gainesville, Florida 32611, USA

<sup>2</sup>GReCO, Institut d'Astrophysique de Paris, CNRS, Université Pierre et Marie Curie, 98 bis Bd. Arago, 75014 Paris, France

<sup>3</sup>Trinity College, University of Cambridge, Cambridge CB2 1TQ, United Kingdom

(Received 25 November 2016; published 3 March 2017)

We derive the secular evolution of the orbital elements of a stellar-mass object orbiting a spinning massive black hole. We use the post-Newtonian (PN) approximation in harmonic coordinates, with test-body equations of motion for the conservative dynamics that are valid through 3PN order, including spin-orbit, quadrupole and (spin)<sup>2</sup> effects, and with radiation-reaction contributions linear in the mass of the body that are valid through 4.5PN order, including the 4PN damping effects of spin-orbit coupling. The evolution equations for the osculating orbit elements are iterated to high PN orders using a two-time-scale approach and averaging over orbital time scales. We derive a criterion for terminating the orbit when its Carter constant drops below a critical value, whereupon the body plunges across the event horizon at the next closest approach. The results are valid for arbitrary eccentricities and arbitrary inclinations. We then analyze numerically the orbits of objects injected into high-eccentricity orbits via interactions within a surrounding star cluster, obtaining the number of orbits and the elapsed time between injection and plunge, and the residual orbital eccentricity at plunge as a function of inclination. We derive an analytic approximation for the time to plunge in terms of initial orbital variables. We show that, if the black hole is spinning rapidly, the flux of gravitational radiation during the final orbit before plunge may be suppressed by as much as 3 orders of magnitude if the orbit is retrograde on the equatorial plane compared to its prograde counterpart.

DOI: 10.1103/PhysRevD.95.064003

# I. INTRODUCTION AND SUMMARY

The relativistic motion of a star or small black hole in the field of a spinning massive black hole is a subject of much current interest, on two fronts. Such orbits may lead to gravitational-wave (GW) emission at a level that is potentially detectable by ground-based laser-interferometric detectors, such as LIGO-VIRGO, or by a future space-based detector, such as LISA [1]. In addition, relativistic effects on the orbits may play an important role in the evolution of dense nuclear clusters of stars and black holes orbiting a massive black hole at the center of a galaxy (for a recent review, see [2]). Indeed, these two problems are intertwined: interactions among the bodies in the cluster may inject bodies into orbits that pass very close to the black hole, leading to the emission of waves and a consequent inspiral into the hole.

At the one extreme—very close to the massive black hole—the analysis of gravitational-wave emission and inspiral is a notoriously difficult problem, particularly if one seeks sufficiently accurate predictions to be used effectively in the analysis of data from gravitational-wave detectors. A leading approach to this problem is

the “self-force” program, which attempts to go beyond simple black hole perturbation theory by incorporating the backreaction effects on the geometry of the small mass of the orbiting particle, while keeping the background Kerr geometry of the black hole exact (see [3] for a review). While great progress has been made, concrete results are currently limited to orbits around Schwarzschild black holes and equatorial orbits around Kerr black holes (see, e.g., [4]). Other approaches are hybrid in nature, combining information on the flux of radiation from black hole perturbation theory with assumptions of energy and angular momentum balance, separations of time scales, and adiabatic evolution of orbital parameters [5–8].

At the other extreme—very far from the black hole—relativistic effects are small, but can have important long-term effects on the evolution of the cluster. A leading example of this is the so-called “Schwarzschild barrier,” in which the relativistic advance of the pericenter of a given orbit can, over long time scales, reduce the probability of the star being injected into an orbit that becomes a gravitational-wave driven, “extreme mass-ratio” inspiral (EMRI) [9–14]. In this realm, post-Newtonian (PN) theory, the weak-field, slow-motion approximation to general relativity, is adequate for including the dominant relativistic effects. These can include the standard PN effects of the Schwarzschild geometry, notably the pericenter advance, as

<footnote>* cmw@phys.ufl.edu</footnote>

<footnote>† mam221@cam.ac.uk</footnote>
well as the effects of frame dragging and the quadrupole moment of the Kerr geometry, and post-Newtonian “cross-terms” [15] arising from the nonlinear interactions between the field of the black hole and the fields of the orbiting bodies.

Our interest lies between these extremes. We wish to use post-Newtonian theory to study a star or black hole that has just been injected into an orbit of such low angular momentum that it will pass close to the black hole, and thus will have its motion strongly influenced by gravitational radiation. Ignoring further perturbations by other stars, we will then follow its motion until its orbital parameters require the body to cross the event horizon at the next closest approach. To get a sense of what this means in terms of orbital parameters, consider a $10^6 M_\odot$ black hole, surrounded by a cluster of stellar-mass objects at a distance of a few milliparsecs. The Schwarzschild radius of such a black hole is 50 nanoparsecs. Thus, in order to get close to the black hole, the body must be injected into a very low angular momentum orbit, $L \sim [GMa(1-e^2)]^{1/2} \sim bv$ where *a* is the semimajor axis, of order a few milliparsecs; *e* is the eccentricity; *b* is the impact parameter, which has the scale of Schwarzschild radii; and $v \sim (GM/a)^{1/2}$ is the orbital velocity. This gives $(1-e^2) \sim (b/a)^2$. For $b \sim 100$ Schwarzschild radii, and $a \sim 1$ mpc, this implies $1-e \sim 10^{-5}$. Such extreme eccentricities are not the usual realm for applications of PN theory.

Still, initially the star is not in the final, highly relativistic inspiral regime, and thus we might expect PN theory to be a reasonable approach to the study of such orbits, particularly if it is carried out to high orders in the PN expansion. The PN approximation has proven to be “unreasonably effective” [16] in describing the quasicircular inspiral of comparable mass-compact bodies, such as black holes or neutron stars. In regimes where it might be expected to fail, it agrees surprisingly well with results from numerical relativity. Whether this effectiveness will carry over to the highly eccentric, extreme mass-ratio inspiral problem will depend in part on what questions are being addressed. If the question is “What is the precise, late-time gravitational waveform including phase evolutions?” then there is reason to be skeptical of any answer based on PN theory, absent some corroboration from self-force theory, numerical relativity or some other technique. However if the questions include

(i) How many orbits or how much time does it take to go from initial injection to plunge into the black hole?

(ii) What is the distribution of late-time orbital eccentricity as a function of inclination of the orbit relative to the black hole's equatorial plane?

(iii) How do the late-time energy flux and gravitational-wave frequency depend on inclination?

then PN theory might give quantitative estimates. Time to plunge, for example, is a useful input into schemes for

evolving the cluster, including the effects of the loss of stars to the inspiral process. Whether the late-time orbits are highly circularized or could have significant residual eccentricities, depending on inclination, could be important guides for research on such schemes as the self-force program. And the dependence of the gravitational-wave signal on inclination could impact estimates of the detectability of EMRI radiation by space detectors such as LISA. Questions such as these will be the main focus of this paper.

We first develop the necessary analytic formulas to describe the evolution of very eccentric orbits from injection to plunge, for a stellar-mass object of mass *m* in the gravitational field of a massive Kerr black hole of mass *M* and dimensionless spin parameter *χ*, derived to high orders in the PN approximation. In Sec. II we derive and display the basic two-body equations of motion to be used. They are valid through 3PN order for the conservative part of the orbital dynamics of a test particle, including frame-dragging, quadrupole-moment and (spin)$^2$ effects. They are valid through 4.5PN order in the dissipative part caused by gravitational radiation reaction, including the 4PN-order effects of spin-orbit coupling. The key results are Eqs. (2.3) and (2.4) (together with detailed formulas displayed in Appendix A), expressed in harmonic coordinates. The dissipative terms are given to first order in the reduced mass parameter $\eta = mM/(m+M)^2$, which is assumed to be very small. They are presented in a ready-to-use form that permits them to be integrated into an *N*-body code for evolving the stellar cluster.

Section III uses the formalism of osculating orbit elements and the “Lagrange planetary equations” for evolving the orbit elements in the presence of the relativistic, non-Keplerian perturbations (Sec. III A). We use a two-scale analysis to separate each orbit element into an orbit-averaged piece and a piece that varies on an orbital time scale, and derive equations both for the “secular” variations of the average elements on longer time scales, and for the periodic variations (Sec. III B and Appendix B). Because we are working to high orders in a PN expansion, we must carefully iterate the planetary equations to high orders. From the conservative sector of the equations of motion (Sec. III C), we obtain the usual advance of the pericenter of the orbit through 3PN order, and PN and higher-order variations in the angle of nodes. We also find, beginning at 3PN order, quadrupole and (spin)$^2$-induced variations in the eccentricity, semilatus rectum and inclination [Eqs. (3.22) and (3.23)]. We obtain 3PN accurate expressions for the conserved energy, angular momentum component along the black hole spin axis, and Carter constant, in terms of the average orbit elements [Eqs. (3.24)]. From the radiation-reaction sector (Sec. III D), we obtain equations of evolution for the eccentricity and semilatus rectum that include the conventional “Newtonian” or quadrupole formula contributions (Peters-Mathews terms) at 2.5PN order, along with 3.5PN
<table><caption>TABLE I. Critical semilatus rectum for capture.</caption><thead><tr><th rowspan="2">Inclination <i>i</i></th><th colspan="3"><i>p</i><sub>c</sub>/(GM/c<sup>2</sup>)</th></tr><tr><th><i>χ</i> = 0</th><th><i>χ</i> = 0.5</th><th><i>χ</i> = 1</th></tr></thead><tbody><tr><td>0</td><td>9.04</td><td>6.09</td><td>2.71</td></tr><tr><td>45°</td><td>9.04</td><td>6.78</td><td>4.05</td></tr><tr><td>90°</td><td>9.04</td><td>8.77</td><td>7.84</td></tr><tr><td>135°</td><td>9.04</td><td>11.04</td><td>12.90</td></tr><tr><td>180°</td><td>9.04</td><td>12.03</td><td>14.98</td></tr></tbody></table>

and 4.5PN terms, as well as 4PN spin-orbit terms [Eqs. (3.27)]. The latter terms change sign between prograde and retrograde orbits. We also find spin-orbit induced, radiation-reaction variations at 4PN order in the inclination, pericenter and nodal angle.

Section IV addresses an issue that arises with the osculating element formalism when the eccentricity tends to unity. This limit is not well behaved, in a manner similar to the well-known singular behavior of the eccentricity-pericenter pair of orbital elements in the limit of circular orbits. We define a new eccentricity related to the old by PN corrections, and display the conserved quantities and the secular evolution equations in terms of this new variable.

In Sec. V, we use these results to numerically evolve eccentric orbits from injection to plunge. We recast the equations for gravitational radiation damping of the orbit in terms of dimensionless variables $e \equiv GM/c^2p_i$, $x = p/p_i$, $e$, $i$, and $\theta$, where $p$, $e$ and $i$ are the semilatus rectum ($p_i$ is its initial value), eccentricity and inclination of the orbit, and $\theta$ is the orbital phase. These equations depend only on the spin parameter $\chi$, and thus the evolution in terms of $\theta$ depends only on $e$, $\chi$ and the initial conditions $e_i$ and $i_i$. The black hole mass does not appear. In Sec. VA we use our PN

expression for the Carter constant along with results of [17] to establish a critical value of $p$ such that for $p < p_c$, the orbit will have no inner turning point and will plunge across the event horizon [Eq. (5.5)]. This will mark the termination of our numerical integrations. Table I shows representative values of $p_c$ in units of $GM/c^2$ for various values of $\chi$ and $i$.

In Sec. V B we evolve the orbits with respect to phase, or number of orbits. Figure 1 illustrates the results and places them in an astrophysical context. It is a “phase space” diagram (in the terminology of those who study cluster dynamics) of energy vs angular momentum. Here, semi-major axis $a = p/(1 - e^2)$ is a proxy for energy and $1 - e^2$ is a proxy for angular momentum. The scale of the a-axis assumes a $10^6 M_\odot$ black hole; we also assume that $\chi = 1$. The dashed line (red in the online version) shows schematically the Schwarzschild barrier, below which stellar orbits effectively decouple from resonant and nonresonant relaxation effects caused by perturbations from the other stars that drive random-walk behavior in eccentricity. The precise location and slope of this barrier depend on the properties of the surrounding star cluster, and thus the curve shown is merely schematic. Below this barrier, a stellar orbit is still subject to perturbations by other stars, of course, but we have chosen to ignore such perturbations.

The three dotted curves show the critical values of $a$ vs $1 - e^2$ for plunge for equatorial prograde (lower red), polar (middle green), and equatorial retrograde (top blue) orbits, respectively, as established by our capture criterion of Sec. VA. These are curves of roughly constant $p$ or $a \propto (1 - e^2)^{-1}$. The three solid curves show the orbits in phase space of a body with $m = 50 M_\odot$, with initial values of $p_i = 100, 50$ and $20GM/c^2$, and with $e_i = 0.999$ for equatorial prograde (dashed red) and retrograde orbits

<figure><img src="image_2.png" /><figcaption>FIG. 1. Tracks of inspiral orbits with initial values of $p/(GM/c^2) = 100, 50$ and $20$, for a black hole of $M = 10^6 M_\odot$ and $\chi = 1$. Dotted curves show the critical values for capture by the black hole from equatorial prograde (lower red), polar (middle green), and equatorial retrograde (top blue) orbits. Right: Late-time deviation of equatorial prograde (dashed red) from equatorial retrograde (solid black) orbits, and their termination on the respective capture boundary.</figcaption></figure>
<table><caption>TABLE II. Residual eccentricity, $\iota = 180^{\circ}$.</caption><thead><tr><th>Semilatus rectum $p_i/(GM/c^2)$</th><th>Residual eccentricity $e$</th></tr></thead><tbody><tr><td>20</td><td>0.48</td></tr><tr><td>40</td><td>0.18</td></tr><tr><td>60</td><td>0.10</td></tr><tr><td>80</td><td>0.07</td></tr><tr><td>100</td><td>0.05</td></tr></tbody></table>

(solid black). Because $\eta = 5 \times 10^{-5}$, radiation reaction is a tiny effect initially. Thus for all but the late phase of the evolution, $p$ and $e$ are approximately constant. Thus $a$ for these orbits also varies inversely as $1 - e^2$. The difference between prograde and retrograde orbits is too small to be seen initially. The right-hand panel of Fig. 1 provides a blowup of the late evolution of the orbits, showing the termination of each orbit when the condition for plunge is reached. Retrograde orbits terminate at larger $a$ and larger $e$ than do prograde orbits. We also obtain the number of orbits to plunge and the residual eccentricity as functions of $p_i/(GM/c^2)$ and $\iota$. For initial values $e_i$ so close to unity, the results are independent of $e_i$. While prograde orbits are highly circularized by the time of plunge because of their longer inspiral, retrograde orbits may be left with significant eccentricity at the time of plunge. Table II shows residual values for $\iota = 180^{\circ}$ (equatorial retrograde) orbits for a range of initial $p_i$.

Using an expression for the gravitational-wave flux accurate to 2PN order beyond the quadrupole approximation including spin-orbit contributions, we estimate the flux from the final orbit before plunge as a function of inclination. Using a 2PN accurate expression for the orbital period $P$, we also estimate the gravitational-wave frequency at plunge. The results are shown in Fig. 2. We note that, for rapidly spinning black holes, there is a strong suppression of the energy flux for highly inclined and

retrograde orbits compared to their prograde counterparts. This is mostly due to the fact that inclined orbits plunge from greater distances than do equatorial prograde orbits, as can be seen in Fig. 1. For rapidly spinning black holes, the frequency can vary by as much as a factor of 3 between prograde and retrograde equatorial orbits, again because the former orbits survive to much closer to the spinning black hole than do the latter.

These effects of inclination on the flux and frequency may have an impact on estimates of detectability of gravitational radiation from EMRIs by a space-based detector such as LISA. Stars are likely to be injected into EMRI orbits with random inclinations. Thus an estimate that uses only the gravitational-wave flux calculated by black hole perturbation theory for an equatorial prograde orbit in Kerr before plunge, multiplied by a rate of injection of stars into EMRI orbits, may overestimate the total flux by failing to take the suppression induced by orbit inclination into account.

In Sec. V D we convert from $\theta$ to time using a 2PN accurate expression for orbital period $P$ and evolve the orbits with respect to time. Figure 3 shows typical results, here for $M = 10^6 M_{\odot}$, $\chi = 1$, and $e_i = 0.999$. Unlike the evolution with respect to the orbital phase, the time to plunge is sensitive to $M$ and $e_i$ for a given initial $p_i$, because $P \propto M^{-1/2}p_i^{3/2}(1 - e_i^2)^{-3/2}$. In fact, the dependence is sufficiently simple that we were able to find an analytic approximation for the time to plunge, that agrees very well with the numerical results over many orders of magnitude of times, given by

$$T_{\text{plunge}} = \frac{1}{74.3\eta} \left(\frac{GM}{c^3}\right) G'(e_i) e^{-3.96} \times (1 + 3\epsilon + 8\epsilon^{3/2}\chi \cos \iota)^4, \quad (1.1)$$

where $\epsilon = GM/c^2p_i$ and

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 2. Final-orbit gravitational-wave energy flux and frequency vs inclination for various Kerr spin parameters. Fluxes are normalized to the value for 0° inclination, while frequencies are normalized to the value for $\chi = 0$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 3. Evolution of semilatus rectum *p* and eccentricity *e* vs time, for three initial values of *p*, for a 10⁶ M☉ black hole; in all cases, *e*ᵢ = 0.999 and *χ* = 1.</figcaption></figure>

$$G'(e_i) = \frac{3.35}{\sqrt{1 - e_i^2}} - 5 + 8\sqrt{1 - e_i^2}. \quad (1.2)$$

Equation (1.1) can be incorporated into numerical codes that evolve stellar clusters around rotating black holes, providing a simple way to quickly determine the fate of a star that is injected into an orbit with an initial *p*, and *e*, and that suffers no further significant stellar perturbations.

The remainder of this paper provides the details underlying these results. Section VI makes concluding remarks.

## II. POST-NEWTONIAN EQUATIONS OF MOTION FOR A BODY ORBITING A SUPERMASSIVE ROTATING BLACK HOLE

We study the motion of a body with a very small mass *m* in the field of a supermassive rotating black hole of mass *M* and spin *S*. At zeroth order in the mass of the orbiting body, the motion is described by the standard geodesic equation for the Kerr geometry. These equations of motion are conservative, in that they admit conservation laws for energy *E*, angular momentum component *L*<sub>e</sub> in the direction of the spin axis *e*, and the so-called “Carter constant” *C*.

We also include the effects of gravitational radiation reaction, which at lowest order are linear in the mass of the small body. This introduces a new parameter into the problem, conventionally denoted by the “symmetric mass ratio” *η*, defined by

$$\eta \equiv \frac{Mm}{(M+m)^2} \sim \frac{m}{M} \quad \text{when } m \ll M. \quad (2.1)$$

Formally this is inconsistent with our test-body conservative equations of motion, which assume *m* = 0.

However, in post-Newtonian theory, the equations of motion for nonspinning bodies are known to high PN orders for arbitrary masses, and in those equations the coefficients are polynomials in powers of *η* < 1/4. The problems that we wish to address with this work all involve stars or small black holes orbiting supermassive black holes so that *η* ≪ 1. Thus corrections due to finite masses in the conservative terms will be subdominant. But the gravitational radiation-reaction terms are proportional to *η*. They are also dissipative; i.e. they violate the conservation of quantities associated with the orbit. We will only be interested in the long-term dissipative effects of gravitational radiation at lowest order in *η*, and therefore we are justified in ignoring finite mass contributions to the conservative equations of motion.

Throughout, we will work in harmonic coordinates, since those are the most suitable for discussing gravitational radiation and radiation reaction, and are the basis for post-Newtonian theory [18]. Appendix A provides details for transforming the Kerr geometry from the conventional Boyer-Lindquist coordinates to harmonic coordinates, expanding the metric to 3PN order, and obtaining the test-body equations of motion and the conserved quantities *E*, *L*<sub>e</sub> and *C* to 3PN order. We make use of the fact that the black hole spin *S* and the Kerr parameter *a* can be written as

$$S \equiv \frac{GM^2}{c} \chi e, \quad a \equiv \frac{S}{M} = \frac{GM}{c} \chi, \quad (2.2)$$

where *χ* is the dimensionless Kerr spin parameter bounded by 0 ≤ *χ* ≤ 1 and *e* is a unit vector parallel to the spin axis. The geodesic equation to 3PN order is then given by
$$ \frac{dv}{dt} = -\frac{GMn}{r^2} - \frac{GM}{c^2 r^2} \left[ \left( v^2 - 4 \frac{GM}{r} \right) \mathbf{n} - 4 \dot{r} \mathbf{v} \right] + \frac{GM^2}{c^3 r^3} \chi \left[ 6\mathbf{e} \cdot (\mathbf{n} \times \mathbf{v}) \mathbf{n} + 6 \dot{r} (\mathbf{n} \times \mathbf{e}) - 4 (\mathbf{v} \times \mathbf{e}) \right] \\ - \frac{G^2 M^2}{c^4 r^3} \left\{ \left[ \left( \frac{9GM}{r} - 2 \dot{r}^2 \right) \mathbf{n} + 2 \dot{r} \mathbf{v} \right] - \frac{3}{2} \frac{GM}{r} \chi^2 \left[ 5 \mathbf{n} (\mathbf{e} \cdot \mathbf{n})^2 - 2 \mathbf{e} (\mathbf{e} \cdot \mathbf{n}) - \mathbf{n} \right] \right\} \\ - \frac{G^2 M^2}{c^5 r^3} \chi \left\{ \frac{GM}{r} \left( 20 \mathbf{e} \cdot (\mathbf{n} \times \mathbf{v}) \mathbf{n} + 16 \dot{r} (\mathbf{n} \times \mathbf{e}) - 12 (\mathbf{v} \times \mathbf{e}) \right) + 6 \dot{r} \mathbf{e} \cdot (\mathbf{n} \times \mathbf{v}) \mathbf{v} \right\} \\ + \frac{G^3 M^3}{c^6 r^4} \left\{ \left[ \left( \frac{16GM}{r} - \dot{r}^2 \right) \mathbf{n} + 4 \dot{r} \mathbf{v} \right] + \chi^2 \left[ \frac{3}{2} \left( 5 \mathbf{n} (\mathbf{e} \cdot \mathbf{n})^2 - 2 \mathbf{e} (\mathbf{e} \cdot \mathbf{n}) - \mathbf{n} \right) \left( v^2 - \frac{4GM}{r} \right) \right. \right. \\ \left. \left. - 6 \mathbf{v} \left( 5 \dot{r} (\mathbf{e} \cdot \mathbf{n})^2 - 2 (\mathbf{v} \cdot \mathbf{e}) (\mathbf{e} \cdot \mathbf{n}) - \dot{r} \right) + \frac{2GM}{r} (\mathbf{n} - 6 \mathbf{n} (\mathbf{e} \cdot \mathbf{n})^2 + (\mathbf{e} \cdot \mathbf{n}) \mathbf{e}) \right] \right\}. \quad (2.3) $$

“Monopole” or pure mass PN terms occur at the usual 1PN ($c^{-2}$), 2PN ($c^{-4}$) and 3PN ($c^{-6}$) orders. Spin-orbit terms linear in $\chi$ occur at 1.5PN ($c^{-3}$) and 2.5PN ($c^{-5}$) orders, shifted as usual by a half PN order relative to their counterparts for normal slowly rotating bodies. The 2PN term quadratic in $\chi$ is the leading-order contribution of the black hole’s Newtonian quadrupole moment, while the 3PN terms quadratic in $\chi$ are a combination of (spin)$^2$ terms and of cross-terms between the quadrupole and monopole terms.

We now supplement these equations with the contributions from the gravitational radiation reaction. These include 2.5PN terms ($c^{-5}$), 3.5PN terms ($c^{-7}$) [19], a 4PN spin-orbit contribution ($c^{-8}$) [20] and 4.5PN terms ($c^{-9}$) [21]. They are given, in harmonic coordinates, by

$$ \frac{dv}{dt_{RR}} = \frac{8}{5} \eta \frac{G^2 M^2}{c^5 r^3} \left[ \left( 3v^2 + \frac{17}{3} \frac{GM}{r} \right) \dot{\mathbf{r}}\mathbf{n} - \left( v^2 + 3 \frac{GM}{r} \right) \mathbf{v} \right] - \frac{8}{5} \eta \frac{G^2 M^2}{c^7 r^3} \left[ \left( \frac{183}{28} v^4 + \frac{519}{42} v^2 \frac{GM}{r} - \frac{285}{4} v^2 \dot{r}^2 + \frac{147}{4} \dot{r}^2 \frac{GM}{r} + 70 \dot{r}^4 + \frac{989}{14} \frac{G^2 M^2}{r^2} \right) \dot{\mathbf{r}}\mathbf{n} - \left( \frac{313}{28} v^4 - \frac{205}{42} v^2 \frac{GM}{r} - \frac{339}{4} v^2 \dot{r}^2 + \frac{205}{12} \dot{r}^2 \frac{GM}{r} + 75 \dot{r}^4 + \frac{1325}{42} \frac{G^2 M^2}{r^2} \right) \mathbf{v} \right] \\ - \frac{1}{5} \eta \frac{G^3 M^3}{c^8 r^4} \chi \left\{ \mathbf{e} \cdot (\mathbf{n} \times \mathbf{v}) (A_{4\text{SO}} \dot{\mathbf{r}}\mathbf{n} + B_{4\text{SO}} \mathbf{v}) - \frac{2}{3} (\mathbf{v} \times \mathbf{e}) \dot{r} C_{4\text{SO}} + \frac{1}{2} (\mathbf{n} \times \mathbf{e}) D_{4\text{SO}} \right\} + \frac{8}{5} \eta \frac{G^2 M^2}{c^9 r^3} [E_{4.5} \dot{\mathbf{r}}\mathbf{n} - F_{4.5} \mathbf{v}]. \quad (2.4) $$

The expressions for $A_{4\text{SO}}$, $B_{4\text{SO}}$, $C_{4\text{SO}}$, $D_{4\text{SO}}$, $E_{4.5}$ and $F_{4.5}$ are given in Appendix A. For simplicity we do not include the nonlocal 4PN “tail” term.

In studying the evolution of the system under radiation reaction, we cannot simply ignore the conservative PN test-body terms displayed in Eq. (2.3), and treat the problem as if it were Newtonian gravity plus radiation reaction. That is because, in solving for the secular evolution of the orbital elements, there will be cross-term interactions between, say, the 1PN conservative perturbations or the 1.5PN spin-orbit terms and the 2.5PN radiation-reaction terms, that could produce dissipative contributions at 3.5PN order or 4PN spin-orbit order, respectively. These interactions must be taken into account to get the complete higher-order evolution.

### III. LAGRANGE PLANETARY EQUATIONS FOR THE OSCULATING ORBIT ELEMENTS

#### A. Basic equations

We begin with a brief review of standard orbital perturbation theory, used to compute deviations from

Keplerian two-body motion induced by perturbing forces, described by the equation of motion

$$ \frac{d^2 \mathbf{x}}{dt^2} = - \frac{GM}{r^3} \mathbf{x} + \delta \mathbf{a}, \quad (3.1) $$

where $\delta \mathbf{a}$ is a perturbing acceleration. For a general orbit described by $\mathbf{x}$ and $\mathbf{v} = dx/dt$, we define the “osculating” Keplerian orbit using a set of orbit elements: the semilatus rectum $p$, eccentricity $e$, inclination $i$, nodal angle $\Omega$ and pericenter angle $\omega$, defined by the following set of equations:

$$ \mathbf{x} \equiv \mathbf{r} \mathbf{n}, $$

$$ r \equiv p/(1 + e \cos f), $$

$$ \mathbf{n} \equiv (\cos \Omega \cos \phi - \cos i \sin \Omega \sin \phi) \mathbf{e}_X $$

$$ + (\sin \Omega \cos \phi + \cos i \cos \Omega \sin \phi) \mathbf{e}_Y + \sin i \sin \phi \mathbf{e}_Z, $$

$$ \lambda \equiv \partial \mathbf{n} / \partial \phi, \quad \hat{\mathbf{h}} = \mathbf{n} \times \lambda, $$

$$ \mathbf{h} \equiv \mathbf{x} \times \mathbf{v} \equiv \sqrt{GM} \hat{\mathbf{h}}, \quad (3.2) $$
where $f \equiv \phi - \omega$ is the true anomaly, $\phi$ is the orbital phase measured from the ascending node, and $e_A$ are chosen reference basis vectors. From the given definitions, we see that $\mathbf{v} = \dot{\mathbf{r}}\mathbf{n} + (h/r)\lambda$ and $\dot{\mathbf{r}} = (he/p)\sin f$. We also define the alternative orbit elements:

$$ \alpha = e \cos \omega, \quad \beta = e \sin \omega, \qquad (3.3) $$

so that now

$$ r = p/(1 + \alpha \cos \phi + \beta \sin \phi). \qquad (3.4) $$

One then defines the radial $\mathcal{R}$, cross-track $\mathcal{S}$ and out-of-plane $\mathcal{W}$ components of the perturbing acceleration $\delta\mathbf{a}$, defined respectively by $\mathcal{R} \equiv \mathbf{n} \cdot \delta\mathbf{a}$, $\mathcal{S} \equiv \lambda \cdot \delta\mathbf{a}$ and $\mathcal{W} \equiv \hat{\mathbf{h}} \cdot \delta\mathbf{a}$, and writes down the “Lagrange planetary equations” for the evolution of the orbit elements,

$$ \frac{dp}{dt} = 2\sqrt{\frac{p}{GM}}r\mathcal{S}, $$

$$ \frac{d\alpha}{dt} = \sqrt{\frac{p}{GM}} \left[ \mathcal{R} \sin \phi + \mathcal{S}(\alpha + \cos \phi) \left( 1 + \frac{r}{p} \right) + \mathcal{W} \frac{r}{p} \beta \cot \iota \sin \phi \right], $$

$$ \frac{d\beta}{dt} = \sqrt{\frac{p}{GM}} \left[ -\mathcal{R} \cos \phi + \mathcal{S}(\beta + \sin \phi) \left( 1 + \frac{r}{p} \right) - \mathcal{W} \frac{r}{p} \alpha \cot \iota \sin \phi \right], $$

$$ \frac{dt}{dt} = \sqrt{\frac{p}{GM}} \mathcal{W} \left(\frac{r}{p}\right) \cos \phi, $$

$$ \frac{d\Omega}{dt} = \sqrt{\frac{p}{GM}} \mathcal{W} \left(\frac{r}{p}\right) \frac{\sin \phi}{\sin \iota}, $$

$$ \frac{d\phi}{dt} = \frac{h}{r^2} - \cos \iota \frac{d\Omega}{dt}, \qquad (3.5) $$

(see [18] for further discussion.) The final equation is not really an orbit element equation, but serves to close the system; in some formulations the sixth equation determines the “time of pericenter passage.” These equations are *exact*: they are simply a reformulation of Eq. (3.1) in terms of new variables. However, they are particularly useful when the perturbations are small, so that solutions can be obtained by a process of iteration. At lowest order (no perturbations), the elements $p, \alpha, \beta, \Omega$ and $\iota$ are constant, and $d\phi/dt = \sqrt{GMp}/r^2$; those solutions can be plugged into the right-hand side, the equations integrated to find corrections, and so on. The corrections to the orbit elements tend to be of two classes: periodic corrections, which vary on an orbital time scale, and secular corrections, which vary on a longer time scale, depending upon the nature of the perturbations. There are many approaches to separating

secular from periodic effects (see [22], for example), but we will adopt a specific approach, known as “multiple-scale analysis,” that is used in many problems in physics [23]. One advantage of this method is that it can be carried out systematically to higher orders in perturbation theory.

We have formulated the Lagrange planetary equations using $\alpha$ and $\beta$ instead of $e$ and $\omega$ because of the well-known problem that the latter formulation is singular in the limit $e \to 0$ (the pericenter angle $\omega$ is formally undefined in this limit). At 1PN order, for example, the osculating orbit that corresponds to a true circular orbit ($r$ = constant) is one in which $e \approx 3GM/c^2p$; the pericenter advances at the same rate as the orbital motion, $\omega \propto \phi$; and the true anomaly is perpetually fixed at the apocenter, $f = \pi$ [24]. The fact that $\omega$ no longer experiences a small precession illustrates the singular nature of this limit. By contrast, the $\alpha - \beta$ formulation is completely regular in the limit $e = (\alpha^2 + \beta^2)^{1/2} = 0$. An alternative approach involves defining two eccentricities, one associated with the radial motion and one with the angular motion [25,26]. However, we will stick with the $\alpha, \beta$ orbit elements.

We now use the sixth planetary equation to express the remaining five equations as differential equations in $\phi$, in the general form

$$ \frac{dX_{\alpha}(\phi)}{d\phi} = \epsilon Q_{\alpha}(X_{\beta}(\phi), \phi), \qquad (3.6) $$

where $\alpha, \beta$ label the orbit element; $\epsilon$ is a small parameter that characterizes the perturbation; and

$$ Q_{\alpha}(X_{\beta}(\phi), \phi) = \frac{(r^2/h)Q_{\alpha}^{(t)}}{1 - (r^2/h) \cos t Q_{\Omega}^{(t)}}, \qquad (3.7) $$

where the $Q_{\alpha}^{(t)}$ denote the right-hand sides of Eqs. (3.5). Note that, because we are working to higher PN orders, we must, at least in principle, include the additional factor involving $Q_{\Omega}^{(t)}$.

## B. Two-time-scale analysis of the planetary equations

We anticipate that the solutions for the $X_\alpha$ will have pieces that vary on a “short” orbital time scale, corresponding to the periodic functions of $\phi$, but may also have pieces that vary on a long time scale, of order $1/\epsilon$ times the short time scale. In a two-time-scale analysis [7,23,24,27], one treats these two times formally as independent variables, and solves the ordinary differential equations as if they were partial differential equations for the two variables. We define the long-time-scale variable

$$ \theta \equiv \epsilon\phi, \qquad (3.8) $$

and write the derivative with respect to $\phi$ as
$$ \frac{d}{d\phi} = \epsilon \frac{\partial}{\partial\theta} + \frac{\partial}{\partial\phi}. \qquad (3.9) $$

We make an ansatz for the solution for $X_\alpha(\theta, \phi)$:

$$ X_{\alpha}(\theta, \phi) = \tilde{X}_{\alpha}(\theta) + \epsilon Y_{\alpha}(\tilde{X}_{\beta}(\theta), \phi). \qquad (3.10) $$

On short time scales, this would represent the standard solution: an initial value of the variable $\tilde{X}_\alpha$, plus a solution that depends on the initial conditions. But on longer time scales, these “initial” values are allowed to vary.

This split of $X_\alpha$ is not unique. However, if the explicit $\phi$ dependence in $Q_\alpha$ is periodic, a natural choice to define the split is to assume that

$$ \tilde{X}_{\alpha}(\theta) = \langle X_{\alpha}(\theta, \phi) \rangle, \quad \langle Y_{\alpha}(\tilde{X}_{\beta}(\theta), \phi) \rangle = 0, \qquad (3.11) $$

where the “average” $\langle...\rangle$ is defined by

$$ \langle A \rangle = \frac{1}{2\pi} \int_{0}^{2\pi} A(\theta, \phi) d\phi, \qquad (3.12) $$

holding $\theta$ fixed. For any function $A(\theta, \phi)$ we define the “average-free” part as

$$ \mathcal{A}\mathcal{F}(A) = A(\theta, \phi) - \langle A \rangle. \qquad (3.13) $$

We now substitute Eq. (3.8) into (3.6), and divide by the parameter $\epsilon$, to obtain

$$ \frac{d\tilde{X}_{\alpha}(\theta)}{d\theta} + \frac{\partial Y_{\alpha}}{\partial\phi} + \epsilon \frac{\partial Y_{\alpha}}{\partial\tilde{X}_{\gamma}} \frac{d\tilde{X}_{\gamma}}{d\theta} = Q_{\alpha}(\tilde{X}_{\beta} + \epsilon Y_{\beta}, \phi), \qquad (3.14) $$

where we sum over the range of $\gamma$. Taking the average and average-free parts of this equation, we obtain

$$ \frac{d\tilde{X}_\alpha}{d\theta} = \langle Q_\alpha(\tilde{X}_\beta + \epsilon Y_\beta, \phi) \rangle, \qquad (3.15\text{a}) $$

$$ \frac{\partial Y_{\alpha}}{\partial \phi} = \mathcal{A}\mathcal{F}(Q_{\alpha}(\tilde{X}_{\beta} + \epsilon Y_{\beta}, \phi)) - \epsilon \frac{\partial Y_{\alpha}}{\partial \tilde{X}_{\gamma}} \frac{d\tilde{X}_{\gamma}}{d\theta}. \qquad (3.15\text{b}) $$

These equations can then be iterated in a straightforward way. At zeroth order, Eq. (3.15a) yields $d\tilde{X}_\alpha/d\theta = \langle Q_\alpha^0 \rangle$ where $Q_\alpha^0 = Q_\alpha(\tilde{X}_\beta, \phi)$, which is the conventional result

whereby one averages the perturbation holding the orbit elements fixed. We write the expansion $Y_\alpha = Y_\alpha^{(0)} + \epsilon Y_\alpha^{(1)} + \epsilon^2 Y_\alpha^{(2)} + ...$. We then integrate Eq. (3.15b) holding $\theta$ fixed to obtain $Y_\alpha^0$. Appendix B provides details on integrating this equation with boundary conditions chosen to ensure that the answer is average free. The iteration continues until one obtains all contributions to $d\tilde{X}_\alpha/d\theta$ compatible with the order in $\epsilon$ to which $Q_\alpha$ is known. The final solution including periodic terms is given by Eq. (3.10), with the secular evolution of the $\tilde{X}_\alpha$ given by solutions of Eqs. (3.15a). From these solutions one can reconstruct the instantaneous orbit using Eqs. (3.2).

If one is interested only in the secular evolutions to a chosen order, then one can Taylor expand $Q_\alpha$ in Eqs. (3.15) in powers of $\epsilon$ and carry out the iterations explicitly. Through order $\epsilon^2$ beyond the lowest-order term, the result is

$$ \frac{d\tilde{X}_\alpha}{d\phi} = \epsilon \langle Q_\alpha^{(0)} \rangle + \epsilon^2 \left[ \langle Q_{\alpha,\beta}^{(0)} \int_0^\phi Q_\beta^{(0)} d\phi' \rangle + \langle Q_{\alpha,\beta}^{(0)} \rangle \langle \phi Q_\beta^{(0)} \rangle - \langle \phi Q_{\alpha,\beta}^{(0)} \rangle \langle Q_\beta^{(0)} \rangle - \pi \langle Q_{\alpha,\beta}^{(0)} \rangle \langle Q_\beta^{(0)} \rangle \right] + O(\epsilon^3), \qquad (3.16) $$

where the subscript $\beta$ denotes $\partial/\partial\tilde{X}_\beta$, and where we have converted from $\theta$ back to $\phi$. Since we are working to 3PN order in the conservative dynamics and to 4.5PN order in the radiation-reaction sector, we will need to include the $O(\epsilon^3)$ contributions; Appendix B provides details.

### C. Secular conservative dynamics

We first carry out this procedure for the conservative part of the equations of motion. At zeroth order in $\epsilon$, we obtain

$$ \frac{d\tilde{p}}{d\theta} = \frac{d\tilde{\tau}}{d\theta} = \frac{d\tilde{\Omega}}{d\theta} = 0, \qquad \frac{d\tilde{\alpha}}{d\theta} = -\frac{3GM}{c^2\tilde{p}}\tilde{\beta}, \qquad \frac{d\tilde{\beta}}{d\theta} = \frac{3GM}{c^2\tilde{p}}\tilde{\alpha}, \qquad (3.17) $$

from which we obtain $d\tilde{e}/d\theta = 0$ and $d\tilde{\omega}/d\theta = 3GM/c^2\tilde{p}$, the standard 1PN results for Schwarzschild. For the average-free part, we obtain, including the 1.5PN-order spin-orbit contribution,

$$ \begin{aligned} Y_p^0 &= -8 \frac{GM}{c^2} (\tilde{\alpha} \cos \phi + \tilde{\beta} \sin \phi) + 4 \tilde{p} \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \cos \tilde{\tau} (\tilde{\alpha} \cos \phi + \tilde{\beta} \sin \phi), \\ Y_\alpha^0 &= -\frac{GM}{2c^2 \tilde{p}} [5 \tilde{\beta} \sin 2\phi + 5 \tilde{\alpha} \cos 2\phi + 2(3 + 7\tilde{\alpha}^2 - \tilde{\beta}^2) \cos \phi + 16 \tilde{\alpha} \tilde{\beta} \sin \phi] \\ &\quad + \frac{1}{2} \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \cos \tilde{\tau} [\tilde{\beta}^2 \cos 3\phi - \tilde{\alpha} \tilde{\beta} \sin 3\phi - 2 \tilde{\beta} \sin 2\phi - (4 - 4\tilde{\alpha}^2 + 9\tilde{\beta}^2) \cos \phi + 11 \tilde{\alpha} \tilde{\beta} \sin \phi]. \end{aligned} $$
$$
\begin{align}
Y_{\beta}^{0} &= -\frac{GM}{2c^2 \tilde{p}} [5\tilde{\alpha} \sin 2\phi - 5\tilde{\beta} \cos 2\phi + 2(3 + 7\tilde{\beta}^2 - \tilde{\alpha}^2) \sin \phi + 16\tilde{\alpha} \tilde{\beta} \cos \phi] \nonumber \\
&\quad + \frac{1}{2} \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \cos \tilde{\iota} [\tilde{\alpha}^2 \sin 3\phi - \tilde{\alpha} \tilde{\beta} \cos 3\phi + 2\tilde{\alpha} \sin 2\phi - (4 - 4\tilde{\beta}^2 + 7\tilde{\alpha}^2) \sin \phi + 13\tilde{\alpha} \tilde{\beta} \cos \phi], \nonumber \\
Y_{\tilde{\iota}}^{0} &= -\frac{1}{2} \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \sin \tilde{\iota} [\tilde{\alpha} \cos 3\phi + \tilde{\beta} \sin 3\phi + 2 \cos 2\phi + 3\tilde{\alpha} \cos \phi + \tilde{\beta} \sin \phi], \nonumber \\
Y_{\tilde{\Omega}}^{0} &= -\left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi [\sin 2\phi (1 + \tilde{\alpha} \cos \phi + \tilde{\beta} \sin \phi) - 2\tilde{\alpha} \sin \phi + 2\tilde{\beta} \cos \phi]. \tag{3.18}
\end{align}
$$

The solution for the orbital separation as a function of $\phi$ is
given by

$$
r^{-1} = \frac{1 + (\tilde{\alpha} + Y_{\alpha}^{0}) \cos \phi + (\tilde{\beta} + Y_{\beta}^{0}) \sin \phi}{\tilde{p} + Y_{p}^{0}} . \quad (3.19)
$$

In the limit $\tilde{\alpha} = \tilde{\beta} = 0$, this yields

$$
r^{-1} = \frac{1 - 3(GM/c^2 \tilde{p}) - 2(GM/c^2 \tilde{p})^{3/2} \chi \cos \tilde{\iota}}{\tilde{p}} \quad (3.20)
$$

corresponding to a circular orbit, without any of the
singular behavior associated with $\tilde{\omega}$. On the other hand,
in the limit $\tilde{e} \to 1$, we obtain

$$
r^{-1} = \frac{1 + \cos f}{\tilde{p}} \left[ 1 + \frac{11}{2} \frac{GM}{c^2 \tilde{p}} - \left( \frac{GM}{c^2 \tilde{p}} \right)^{3/2} \chi \cos \tilde{\iota} (6 - 2 \cos f - \sin f \sin 2\phi) - \frac{15}{2} \frac{GM}{c^2 \tilde{p}^2} \right] \quad (3.21)
$$

where $f = \phi - \tilde{\omega}$. Because of the presence of the final 1PN
term, $r^{-1}$ does not vanish in the obvious manner at $f = \pi$,
as might be expected for an unbound orbit, but rather at
$f = \pi + (15GM/c^2\tilde{p})^{1/2}$. We will address this and other
unusual behavior of the solutions as $\tilde{e} \to 1$ in Sec. IV.

Carrying out the iterations to an order consistent with the
3PN order of the equations of motion, we obtain finally

$$
\frac{d\tilde{p}}{d\theta} = -6\tilde{p} \left(\frac{GM}{c^2 \tilde{p}}\right)^3 \tilde{\alpha} \tilde{\beta} \chi^2 \sin^2 \tilde{\iota}, \qquad (3.22a)
$$

$$
\begin{equation}
\begin{split}
\frac{d\tilde{\alpha}}{d\theta} = & -\frac{3GM}{c^2\tilde{p}}\tilde{\beta} + 6\left(\frac{GM}{c^2\tilde{p}}\right)^{3/2}\tilde{\beta}\chi\cos\tilde{\iota} + \frac{3}{4}\left(\frac{GM}{c^2\tilde{p}}\right)^2\tilde{\beta}[10 - \tilde{\alpha}^2 - \tilde{\beta}^2 - \chi^2(5\cos^2\tilde{\iota} - 1)] \\
& - 3\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2}\tilde{\beta}(8 - 3\tilde{\alpha}^2 - 3\tilde{\beta}^2)\chi\cos\tilde{\iota} \\
& - \frac{3}{2}\left(\frac{GM}{c^2\tilde{p}}\right)^3\tilde{\beta}\left[(29 + 34\tilde{\alpha}^2 + 34\tilde{\beta}^2) + \frac{1}{4}\chi^2(7 + 2\tilde{\alpha}^2 - 18\tilde{\beta}^2 - [23 - 62\tilde{\alpha}^2 - 90\tilde{\beta}^2]\cos^2\tilde{\iota})\right],
\end{split}
\tag{3.22b}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
\frac{d\tilde{\beta}}{d\theta} = & \frac{3GM}{c^2\tilde{p}}\tilde{\alpha} - 6\left(\frac{GM}{c^2\tilde{p}}\right)^{3/2}\tilde{\alpha}\chi\cos\tilde{\iota} - \frac{3}{4}\left(\frac{GM}{c^2\tilde{p}}\right)^2\tilde{\alpha}[10 - \tilde{\alpha}^2 - \tilde{\beta}^2 - \chi^2(5\cos^2\tilde{\iota} - 1)] \\
& + 3\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2}\tilde{\alpha}(8 - 3\tilde{\alpha}^2 - 3\tilde{\beta}^2)\chi\cos\tilde{\iota} \\
& + \frac{3}{2}\left(\frac{GM}{c^2\tilde{p}}\right)^3\tilde{\alpha}\left[(29 + 34\tilde{\alpha}^2 + 34\tilde{\beta}^2) - \frac{1}{4}\chi^2(3 + 6\tilde{\alpha}^2 + 26\tilde{\beta}^2 + [13 - 70\tilde{\alpha}^2 - 98\tilde{\beta}^2]\cos^2\tilde{\iota})\right].
\end{split}
\tag{3.22c}
\end{equation}
$$

$$
\frac{d\tilde{\iota}}{d\theta} = -3 \left(\frac{GM}{c^2 \tilde{p}}\right)^3 \tilde{\alpha} \tilde{\beta} \chi^2 \sin \tilde{\iota} \cos \tilde{\iota}, \qquad (3.22\text{d})
$$

$$
\frac{d\tilde{\Omega}}{d\theta} = 2 \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi - \frac{3}{2} \left(\frac{GM}{c^2 \tilde{p}}\right)^2 \chi^2 \cos \tilde{\iota} - 3 \left(\frac{GM}{c^2 \tilde{p}}\right)^{5/2} \chi (4 - \tilde{\alpha}^2 - \tilde{\beta}^2) + \frac{3}{2} \left(\frac{GM}{c^2 \tilde{p}}\right)^3 (8 - 7 \tilde{\alpha}^2 - 9 \tilde{\beta}^2) \chi^2 \cos \tilde{\iota}. \quad (3.22\text{e})
$$

We also obtain, but do not display, the average-free, periodic contributions $Y_\alpha$ through 3PN order. From the equations for $\alpha$
and $\beta$ we obtain
$$ \frac{d\tilde{\omega}}{d\theta} = 3\frac{GM}{c^2\tilde{p}} - 6\left(\frac{GM}{c^2\tilde{p}}\right)^{3/2} \chi \cos \tilde{t} - \frac{3}{4}\left(\frac{GM}{c^2\tilde{p}}\right)^2 [10 - \tilde{e}^2 + \chi^2(1 - 5\cos^2\tilde{t})] + 3\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2} (8 - 3\tilde{e}^2)\chi \cos \tilde{t} \\ + \frac{3}{8}\left(\frac{GM}{c^2\tilde{p}}\right)^3 \{4(29 + 34\tilde{e}^2) - \chi^2[4(4 - 17\tilde{e}^2) - 2(9 - 40\tilde{e}^2)\sin^2\tilde{t} + (4\tilde{e}^2 + 5(1 - 2\tilde{e}^2)\sin^2\tilde{t})\cos(2\tilde{\omega})]\}, \quad (3.23a) $$

$$ \frac{d\tilde{e}}{d\theta} = -\frac{3}{4} \left( \frac{GM}{c^2 \tilde{p}} \right)^3 \tilde{e} (5 + 4\tilde{e}^2) \chi^2 \sin^2 \tilde{t} \sin \tilde{\omega} \cos \tilde{\omega}. \quad (3.23b) $$

Substituting the definitions (3.2) of the osculating orbits into the expressions (A8), (A9), and (A11) for $E$, $L_e$ and $C$, we obtain $E = -GM(1 - \tilde{e}^2)/2\tilde{p} + \dots$, $L_e = (GM\tilde{p})^{1/2} \cos \tilde{t} + \dots$, and $C = (GM\tilde{p}) + \dots$, where the PN corrections depend in general on the phase $\phi$. However, now substituting the solutions for each orbit element $X_\alpha = \tilde{X}_\alpha + Y_\alpha(\tilde{X}_\beta, \phi)$, with the $Y_\alpha$ computed through 3PN order, we find that all $\phi$ dependence cancels, and the conserved quantities are given purely in terms of the $\tilde{X}_\alpha$:

$$ c^{-2}E = -\frac{1}{2} \frac{GM(1 - \tilde{e}^2)}{c^2 \tilde{p}} + \frac{1}{8} \left(\frac{GM}{c^2 \tilde{p}}\right)^2 (19 + 38\tilde{e}^2 + 3\tilde{e}^4) - \frac{1}{16} \left(\frac{GM}{c^2 \tilde{p}}\right)^3 \{(197 - 205\tilde{e}^2 - 98\tilde{e}^4 - 5\tilde{e}^6) \\ - \chi^2[(4 + 6\tilde{e}^2)(1 - 3\cos^2\tilde{t}) - 9\tilde{e}^2\sin^2\tilde{t}\cos 2\tilde{\omega}]\} - \frac{1}{4} \left(\frac{GM}{c^2 \tilde{p}}\right)^{7/2} [2(20 + 70\tilde{e}^2 + 33\tilde{e}^4) - \tilde{e}^2(9 + 4\tilde{e}^2)\cos(2\tilde{\omega})] \chi \cos \tilde{t} \\ + \frac{1}{128} \left(\frac{GM}{c^2 \tilde{p}}\right)^4 \{9923 + 12756\tilde{e}^2 + 11478\tilde{e}^4 + 952\tilde{e}^6 + 35\tilde{e}^8\} + \chi^2[8(220 + 347\tilde{e}^2 + 254\tilde{e}^4) \\ - 4(644 + 1099\tilde{e}^2 + 695\tilde{e}^4)\sin^2\tilde{t} - 8\tilde{e}^2(94 + 80\tilde{e}^2 - 235\sin^2\tilde{t} - 177\tilde{e}^2\sin^2\tilde{t})\cos^2\tilde{\omega}]\}, \quad (3.24a) $$

$$ L_e = \sqrt{m\tilde{p}} \cos \tilde{t} \left\{ 1 + \frac{GM}{2c^2\tilde{p}}(7 + \tilde{e}^2) - \frac{1}{8}\left(\frac{GM}{c^2\tilde{p}}\right)^2(37 - 18\tilde{e}^2 - 3\tilde{e}^4) + \frac{1}{16}\left(\frac{GM}{c^2\tilde{p}}\right)^3[411 + 481\tilde{e}^2 + 86\tilde{e}^4 + 5\tilde{e}^6 \\ + \chi^2[(4(46 + 61\tilde{e}^2) - 5(40 + 63\tilde{e}^2)\sin^2\tilde{t} - 2\tilde{e}^2(32 - 111\sin^2\tilde{t})\cos^2\tilde{\omega}]]\right\} \\ - \sqrt{m\tilde{p}}\chi \left\{ \left(\frac{GM}{c^2\tilde{p}}\right)^{3/2}(1 + \cos^2\tilde{t}) + \frac{1}{4}\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2}[4(5 + 12\tilde{e}^2) - (30 + 41\tilde{e}^2)\sin^2\tilde{t} - 2\tilde{e}^2(4 - 5\sin^2\tilde{t})\cos^2\tilde{\omega}] \right\}, \quad (3.24b) $$

$$ C = m\tilde{p} \left\{ 1 + \frac{GM}{c^2\tilde{p}}(7 + \tilde{e}^2) - 4\left(\frac{GM}{c^2\tilde{p}}\right)^{3/2} \chi \cos \tilde{t} + \frac{1}{2}\left(\frac{GM}{c^2\tilde{p}}\right)^2 [6 + 16\tilde{e}^2 + 2\tilde{e}^4 + \chi^2 \sin^2\tilde{t}(1 - 2\tilde{e}^2 \cos^2\tilde{\omega})] \\ - 2\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2} \chi \cos \tilde{t}(12 + 13\tilde{e}^2 - 2\tilde{e}^2 \cos^2\tilde{\omega}) + \frac{1}{8}\left(\frac{GM}{c^2\tilde{p}}\right)^3 [152 + 570\tilde{e}^2 + 125\tilde{e}^4 + 8\tilde{e}^6 \\ + \chi^2[4(54 + 61\tilde{e}^2) - (280 + 347\tilde{e}^2)\sin^2\tilde{t} - 2\tilde{e}^2(32 - 49\sin^2\tilde{t} + 4\tilde{e}^2 \sin^2\tilde{t})\cos^2\tilde{\omega}]] \right\}. \quad (3.24c) $$

Thus $E$, $L_e$ and $C$ are constant over an orbital time scale. However, a closer inspection of $E$ and $C$ indicates a potential problem. Because $\tilde{p}$, $\tilde{e}$ and $\tilde{t}$ are constant through 2.5PN order in the conservative part of the dynamics [see Eqs. (3.22a), (3.22d) and (3.23b)], the PN corrections in $E$, $L_e$ and $C$ are constant through 3PN order, except for two: the 2PN terms in $E$ and $C$, proportional to $\chi^2$, that depend on $\tilde{\omega}$. The quadrupole moment of the black hole is responsible for these $\tilde{\omega}$-dependent terms. The pericenter $\tilde{\omega}$ varies at 1PN order at a rate $d\tilde{\omega}/d\theta = 3GM/c^2\tilde{p}$,

suggesting that $E$ and $C$ might not really be constant through 3PN order. However, the variations induced by the advance of $\tilde{\omega}$ in those 2PN terms in $E$ and $C$ are exactly canceled by the 3PN variations in $\tilde{p}$ and $\tilde{e}$ [Eqs. (3.22a) and (3.23b)] acting on the Newtonian terms $-GM(1 - \tilde{e}^2)/c^2\tilde{p}$ and $m\tilde{p}$, respectively. In $L_e$, the 3PN variations in $\tilde{p}$ and $\tilde{t}$ acting on the Newtonian term $(GM\tilde{p})^{1/2} \cos \tilde{t}$ automatically cancel each other. Thus $E$, $L_e$ and $C$ are also conserved on a pericenter precession time scale, as expected. But to verify this, it was essential to include the
relevant quadrupole-monopole cross-terms in the equations of motion, and to include the higher-order corrections in solving the Lagrange planetary equations. This is an example of how the “quadrupole conundrum,” discussed in [15], is resolved.

We also need the orbital period, given by

$$P \equiv 2\pi \left\langle \frac{dt}{d\phi} \right\rangle = 2\pi \left\langle \frac{(r^2/h)}{1 - (r^2/h) \cos \iota Q_{\Omega}^{(\iota)}} \right\rangle. \quad (3.25)$$

Substituting the full solutions (3.10) for the orbital elements, we obtain, to 2PN order,

$$P = 2\pi \left( \frac{\tilde{p}^3}{GM(1-\tilde{e}^2)^3} \right)^{1/2} \left[ 1 + \frac{3GM}{2c^2\tilde{p}} \left( \frac{4+9\tilde{e}^2+2\tilde{e}^4}{1-\tilde{e}^2} \right) \right. \\ \left. + 6 \left( \frac{GM}{c^2\tilde{p}} \right)^{3/2} \chi \cos\tilde{t} - \frac{3}{16} \left( \frac{GM}{c^2\tilde{p}} \right)^2 \frac{1}{(1-\tilde{e}^2)^2} \right. \\ \left. \times \{56-1214\tilde{e}^2-941\tilde{e}^4-151\tilde{e}^6-40(1-\tilde{e}^2)^{7/2} \right. \\ \left. + \chi^2(1-\tilde{e}^2)(24-4\tilde{e}^2-(32+7\tilde{e}^2-18\tilde{e}^2\cos^2\tilde{\omega})\sin^2\tilde{t}) \right]. \quad (3.26)$$

Higher-order contributions to P will not be needed.

## D. Secular evolution with radiation reaction

We now include the radiation-reaction terms in the equations of motion. At zeroth order in our two-scale analysis, we have the usual orbital average of the $Q_{\alpha}^0$, holding the orbit elements fixed on the right-hand side. But since we are working to 4.5PN order, we must again include higher-order corrections in the two-scale analysis. Because we work only to linear order in the reduced-mass parameter $\eta$, we need to include only cross-terms between the conservative and radiation-reaction sectors. For example, periodic terms induced by 1PN, 1.5PN spin-orbit and 2PN conservative perturbations in the orbit elements occurring in the 2.5PN radiation-reaction expressions could generate secular contributions at 3.5PN, 4PN spin-orbit and 4.5PN orders, while periodic terms induced by 2.5PN radiation-reaction perturbations in the orbit elements occurring in the 1PN, 1.5PN spin-orbit, and 2PN conservative expressions could also generate secular contributions at the same order. These cross-term effects are why it is essential to include conservative terms at high PN orders when analyzing radiation-reaction effects at high PN orders.

We carry out the iterations of the two-scale equations to an order consistent with the 4.5PN order of the radiation-reaction terms in the equations of motion. After transforming from $\tilde{\alpha}$ and $\tilde{\beta}$ to $\tilde{e}$ and $\tilde{\omega}$, we obtain

$$\frac{d\tilde{p}}{d\theta_{RR}} = -\frac{8}{5}\eta\tilde{p}\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2}(8+7\tilde{e}^2) + \frac{1}{210}\eta\tilde{p}\left(\frac{GM}{c^2\tilde{p}}\right)^{7/2}(22072-6064\tilde{e}^2-1483\tilde{e}^4) \\ +\frac{2}{15}\eta\tilde{p}\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\cos\tilde{t}(968+2280\tilde{e}^2+297\tilde{e}^4) - \frac{1}{11340}\eta\tilde{p}\left(\frac{GM}{c^2\tilde{p}}\right)^{9/2}(8272600+777972\tilde{e}^2-947991\tilde{e}^4-4743\tilde{e}^6), \quad (3.27a)$$

$$\frac{d\tilde{e}}{d\theta_{RR}} = -\frac{1}{15}\eta\tilde{e}\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2}(304+121\tilde{e}^2) + \frac{1}{840}\eta\tilde{e}\left(\frac{GM}{c^2\tilde{p}}\right)^{7/2}(144392-34768\tilde{e}^2-2251\tilde{e}^4) \\ +\frac{1}{30}\eta\tilde{e}\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\cos\tilde{t}(9400+10548\tilde{e}^2+789\tilde{e}^4) \\ -\frac{1}{30240}\eta\tilde{e}\left(\frac{GM}{c^2\tilde{p}}\right)^{9/2}(43837360+4062372\tilde{e}^2-1866490\tilde{e}^4-250065\tilde{e}^6), \quad (3.27b)$$

$$\frac{d\tilde{t}}{d\theta_{RR}} = -\frac{1}{30}\eta\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\sin\tilde{t}[88-1248\tilde{e}^2-441\tilde{e}^2+24\tilde{e}^2(20+13\tilde{e}^2)\cos^2\tilde{\omega}], \quad (3.27c)$$

$$\frac{d\tilde{\omega}}{d\theta_{RR}} = \frac{4}{5}\eta\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\cos\tilde{t}\tilde{e}^2(20+13\tilde{e}^2)\sin\tilde{\omega}\cos\tilde{\omega}, \quad (3.27d)$$

$$\frac{d\tilde{\Omega}}{d\theta_{RR}} = -\frac{4}{5}\eta\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\tilde{e}^2(20+13\tilde{e}^2)\sin\tilde{\omega}\cos\tilde{\omega}. \quad (3.27e)$$

Notice that the 4PN spin-orbit terms act to reduce the orbital decay for prograde orbits ($\tilde{t} < \pi/2$) and act to enhance the decay for retrograde orbits ($\tilde{t} > \pi/2$). This phenomenon has also been seen in numerical simulations of the late stage of inspiral of spinning black holes.
Substituting Eqs. (3.27) into Eqs. (3.24), keeping contributions through 2PN order beyond the leading terms, and dropping $\chi^2$ terms, since they are not present in the radiation-reaction terms, we obtain the rates of change of $E$, $L_e$ and $C$:

$$ \frac{dE}{d\theta} = -\frac{32}{5}\eta \left(\frac{GM}{c^2\tilde{p}}\right)^{7/2} \left\{ \left(1 + \frac{73}{24}\tilde{e}^2 + \frac{37}{96}\tilde{e}^4\right) - \left(\frac{GM}{c^2\tilde{p}}\right) \left(\frac{95216 + 73224\tilde{e}^2 - 18074\tilde{e}^4 - 2393\tilde{e}^6}{5376}\right) \right. \\ \left. - \left(\frac{GM}{c^2\tilde{p}}\right)^{3/2} \chi \cos\tilde{\iota} \left(\frac{1936 + 12024\tilde{e}^2 + 6582\tilde{e}^4 + 195\tilde{e}^6}{192}\right) \right. \\ \left. + \left(\frac{GM}{c^2\tilde{p}}\right)^2 \left(\frac{121274560 + 71538080\tilde{e}^2 + 191064\tilde{e}^4 - 317664\tilde{e}^6 - 427473\tilde{e}^8}{580608}\right) \right\}, \quad (3.28a) $$

$$ \frac{dLe}{d\theta} = -\frac{32}{5}\eta\sqrt{GM\tilde{p}}\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2} \left\{ \left(1 + \frac{7}{8}\tilde{e}^2\right)\cos\tilde{\iota} - 5\left(\frac{GM}{c^2\tilde{p}}\right)\left(\frac{6296 - 1000\tilde{e}^2 - 739\tilde{e}^4}{2688}\right)\cos\tilde{\iota} \right. \\ - \frac{1}{192}\left(\frac{GM}{c^2\tilde{p}}\right)^{3/2} \chi[2(584 + 1944\tilde{e}^2 + 297\tilde{e}^4) - 3(488 + 1744\tilde{e}^2 + 293\tilde{e}^4)\sin^2\tilde{\iota} + 12\tilde{e}^2(20 + 13\tilde{e}^2)\cos 2\tilde{\omega}\sin^2\tilde{\iota}] \\ \left. + \left(\frac{GM}{c^2\tilde{p}}\right)^2\left(\frac{14458192 - 819888\tilde{e}^2 + 38886\tilde{e}^4 + 147537\tilde{e}^6}{145152}\right)\cos\tilde{\iota} \right\}, \quad (3.28b) $$

$$ \frac{dC}{d\theta} = -\frac{64}{5}\eta(GM\tilde{p})\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2} \left\{ \left(1 + \frac{7}{8}\tilde{e}^2\right) - \left(\frac{GM}{c^2\tilde{p}}\right) \left(\frac{22072 - 14576\tilde{e}^2 - 4871\tilde{e}^4}{2688}\right) \right. \\ \left. - \left(\frac{GM}{c^2\tilde{p}}\right)^{3/2} \chi \cos\tilde{\iota} \left(\frac{776 + 2112\tilde{e}^2 + 297\tilde{e}^4}{96}\right) + \left(\frac{GM}{c^2\tilde{p}}\right)^2 \left(\frac{7837144 - 985668\tilde{e}^2 + 1212441\tilde{e}^4 + 294930\tilde{e}^6}{145152}\right) \right\}. \quad (3.28c) $$

## E. Comparison with other work

We compare these results with other approaches in situations where there is overlap. Mora and Will [27] found the secular evolutions of the orbit elements in the conservative sector to 3PN order in the spinless ($\chi = 0$), arbitrary mass case, using the $\alpha - \beta$ formulation and the same two-scale approach. In that limit, our Eqs. (3.22) and (3.23) reduce to $d\tilde{p}/d\theta = d\tilde{e}/d\theta = d\tilde{\iota}/d\theta = d\tilde{\Omega}/d\theta = 0$, and

$$ \frac{d\tilde{\omega}}{d\theta} = 3u - \frac{3}{4}u^2(10 - \tilde{e}^2) + \frac{3}{2}u^3(29 + 34\tilde{e}^2), \quad (3.29) $$

where $u = GM/c^2\tilde{p}$. On the other hand, Mora and Will obtained [Eq. (2.29) of [27]]

$$ \frac{d\tilde{\omega}}{d\theta} = 3\tilde{u} - \frac{3}{4}\tilde{u}^2(10 - \tilde{e}^2) + \frac{3}{2}\tilde{u}^3(29 - 30\tilde{e}^2), \quad (3.30) $$

where $\tilde{u}$ was defined to be $\langle GM/c^2\tilde{p} \rangle$. These expressions appear to differ. However, $u$ and $\tilde{u}$ are not quite the same. Using our solutions for $Y_p$ to the appropriate order, it is straightforward to show that $\tilde{u} = \langle GM/c^2(\tilde{p} + Y_p) \rangle = u + 32\tilde{e}^2u^3 + O(u^5)$, and thus that the two expressions are in fact equivalent. This illustrates the caution that, when

working to high orders in perturbation theory, the precise definition of “average” is important. Mora and Will also obtained the radiation-reaction driven evolution of the elements to 3.5PN order in the spinless black hole case. Equations (3.27) agree completely with Eq. (2.28) of [27] to the relevant order; in this case the transformation between $\tilde{u}$ of Mora-Will and our $u = GM/c^2\tilde{p}$ is a higher-order effect.

Our results also agree with calculations of the rate of evolution of the orbital frequency $\Omega_{\text{circ}}$ of a quasicircular inspiral orbit using high PN-order calculations of the energy flux combined with the assumption of energy balance [28,29]. Through 2PN order beyond the quadrupole approximation (corresponding to our 4.5PN order in the equations of motion), including the leading spin-orbit contributions, the result in the small-$\eta$ limit is [Eq. (3) of [29]]

$$ \frac{d\Omega_{\text{circ}}}{dt} = \frac{96}{5}\eta\left(\frac{c^3}{GM}\right)^2 x^{11/3} \left[1 - \frac{743}{336}x^{2/3} - \frac{113}{12}x\chi\cos\tilde{\iota} + \frac{34103}{18144}x^{4/3}\right], \quad (3.31) $$

where $x \equiv GM\Omega_{\text{circ}}/c^3$. In [29], $\Omega_{\text{circ}}$ was defined to be
$$ \begin{aligned} \Omega_{\text{circ}} &= |v|/r \\ &= \frac{(GM)^{1/2}}{p^{3/2}} (1 + \alpha \cos \phi + \beta \sin \phi) \\ & \quad \times (1 + \alpha^2 + \beta^2 + 2\alpha \cos \phi + 2\beta \sin \phi)^{1/2}. \end{aligned} \quad (3.32) $$

Substituting our solutions for the orbit elements $p$, $\alpha$ and $\beta$ including periodic terms, and taking the $\tilde{e} = 0$ limit, we obtain

$$ x = \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \times \left[1 - 6 \frac{GM}{c^2 \tilde{p}} - 4 \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \cos \tilde{t} + 39 \left(\frac{GM}{c^2 \tilde{p}}\right)^2\right]. \quad (3.33) $$

For $\tilde{e} = 0$, the orbital period (3.26) becomes

$$ P = 2\pi \left(\frac{\tilde{p}^3}{GM}\right)^{1/2} \left[1 + 6 \frac{GM}{c^2 \tilde{p}} + 6 \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \cos \tilde{t} - 3 \left(\frac{GM}{c^2 \tilde{p}}\right)^2\right], \quad (3.34) $$

after dropping the 2PN contribution involving $\chi^2$. Then calculating $d\Omega_{\text{circ}}/dt = (2\pi/P)(c^3/GM)dx/d\theta$, inserting $(d\tilde{p}/d\theta)_{RR}$ from Eq. (3.27) in the $\tilde{e} = 0$ limit, and then using Eq. (3.33) to convert the result from $\tilde{p}$ back to $x$, we obtain exactly Eq. (3.31).

Gergely et al. [30] derived the 4PN spin-orbit effects on the radiation-reaction evolution of the orbit elements by calculating fluxes of energy and angular momentum at infinity including the contributions from the current multipole moments of a rotating body (not necessarily a black hole), and including 1.5PN spin-orbit effects in the conservative equations of motion, while Sago and Fujita [31] calculated the fluxes in black hole perturbation theory and then expanded in a PN sequence. Both papers used a specific definition of the orbit elements $a_0$, $e_0$ and $i_0$, given by [see for example Eqs. (2.15) and (5.1) of [30]]

$$ \begin{aligned} r_{\max(\min)} &\equiv a_0(1 \pm e_0), \\ \cos(i_0) &= \frac{L_e}{L}, \end{aligned} \quad (3.35) $$

and calculated $\langle dE/dt \rangle$, $\langle dL_e/dt \rangle$, $\langle dL/dt \rangle$ in terms of those orbit elements. Using our solutions (3.18) through 1.5PN order for the orbital elements (including the periodic contributions), and using the approximation $L = \sqrt{C}$, it is straightforward to show that

$$ \begin{aligned} a_0 &= \frac{\tilde{p}}{1 - \tilde{e}^2} + \left(\frac{GM}{c^2}\right) \frac{3 + 13\tilde{e}^2 - \tilde{e}^4}{(1 - \tilde{e}^2)^2} \\ & \quad + 2\tilde{p} \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \cos \tilde{t} + \dots, \\ e_0 &= \tilde{e} + \tilde{e} \frac{GM}{2c^2 \tilde{p}} (17 - 2\tilde{e}^2) \\ & \quad - 2\tilde{e} (1 - \tilde{e}^2) \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \cos \tilde{t} + \dots, \\ \cos i_0 &= \cos \tilde{t} - \left(\frac{GM}{c^2 \tilde{p}}\right)^{3/2} \chi \sin^2 \tilde{t} + \dots \end{aligned} \quad (3.36) $$

Converting from $d/dt$ to $d/d\theta$ using the corrections to the orbital period in Eq. (3.26) through 1.5PN order, it is straightforward to show that the Newtonian and the 1.5PN spin-orbit contributions to the fluxes of $E$, $L_e$ and $L$ displayed in Eqs. (5.6), (5.7) and (5.8) of [30] and in Eqs. (32)–(34) of [31] are completely equivalent to the corresponding contributions in our Eqs. (3.28). On the other hand, the 1PN and 2PN contributions in Eqs. (32)–(34) of [31] are not equivalent. This is undoubtedly another example of how different kinds of averages can affect higher-order PN terms. Furthermore, at leading PN order, $\langle dX/dt \rangle$ is the same as $\langle dX/d\phi \rangle/\langle dt/d\phi \rangle$. This also holds for the leading contribution of a specific type of effect, such as the spin-orbit contribution. But because of cross-term effects, it does not necessarily hold for 1PN and 2PN corrections to the leading Newtonian flux.

We have not attempted a direct comparison with other related results [32–36], largely because they use rather different orbital parametrizations than those used here.

## IV. HIGHLY ECCENTRIC ORBITS

In this section, we prepare the groundwork for the numerical analysis of the evolution of highly eccentric orbits around a massive black hole, to be addressed in Sec. V. We have already noted the unexpected behavior of $r^{-1}$ in the $\tilde{e} \to 1$ limit in Eq. (3.21). Examination of the expression for $E$ [Eq. (3.24a)] reveals that, for fixed semilatus rectum $\tilde{p}$, the Newtonian energy vanishes as the eccentricity $\tilde{e}$ tends toward unity, but the PN corrections remain finite. Similarly, while the Newtonian part of the period $P$ [Eq. (3.26)] tends to infinity in this limit as $(1 - \tilde{e}^2)^{-3/2}$ as expected for the transition from a bound to an unbound orbit, the PN and 2PN corrections to the period blow up even faster. It suggests that the true boundary between bound and unbound orbits in PN theory corresponds to something other than $\tilde{e} = 1$. This behavior is reminiscent of the singular behavior of the $\tilde{e} - \tilde{\omega}$ parametrization of osculating orbits in the limit of circular orbits in PN theory, which forces either the $\alpha - \beta$ parametrization of the orbit elements or a parametrization in terms of two
eccentricities [25,26]. But to our knowledge, this behavior at $\tilde{e} = 1$ does not seem to have received much attention.

In an attempt to cure this problem, we shall define a post-Newtonian corrected eccentricity $e$, subject to the following constraints: (i) the limit $e \to 0$ still corresponds to $\tilde{e} \to 0$; (ii) the energy vanishes as $(1 - e^2)$ as $e \to 1$, holding $\tilde{p}$ fixed. We define the transformation to the new eccentricity by

$$ \tilde{e} \equiv e \left[ 1 + \frac{GM}{c^2 \tilde{p}} f_1(e) + \left( \frac{GM}{c^2 \tilde{p}} \right)^{3/2} f_{1.5}(e) \chi \cos \tilde{\tau} + \left( \frac{GM}{c^2 \tilde{p}} \right)^2 f_2(e) + \left( \frac{GM}{c^2 \tilde{p}} \right)^{5/2} f_{2.5}(e) \chi \cos \tilde{\tau} + \left( \frac{GM}{c^2 \tilde{p}} \right)^3 f_3(e) + \dots \right], \quad (4.1) $$

where $f_\alpha(e)$ are to be polynomials in $e$, constrained to be finite as $e \to 0$. We substitute this into Eq. (3.24a) and require that $E$ acquire the form

$$ E \to -\frac{1}{2} \frac{GM(1-e^2)}{\tilde{p}} \left[ 1 + \frac{GM}{c^2 \tilde{p}} g_1(e) + \left( \frac{GM}{c^2 \tilde{p}} \right)^{3/2} g_{1.5}(e) \chi \cos \tilde{\tau} + \left( \frac{GM}{c^2 \tilde{p}} \right)^2 g_2(e) + \left( \frac{GM}{c^2 \tilde{p}} \right)^{5/2} g_{2.5}(e) \chi \cos \tilde{\tau} + \left( \frac{GM}{c^2 \tilde{p}} \right)^3 g_3(e) + \dots \right]. \quad (4.2) $$

The solutions are

$$ f_1(e) = \frac{8g_1(e)(1 - e^2) - 19 - 38e^2 - 3e^4}{8e^2}, \quad f_{1.5}(e) = \frac{g_{1.5}(e)(1 - e^2)}{e^2}, $$

$$ f_2(e) = \frac{1}{e^4} \left\{ (8g_1(e) - 19)^2 - 4e^2 \left[ (8g_1(e) - 19)(4g_1(e) - 19) + 32g_2(e) + 394 - 8\chi^2(1 - 3\cos^2\tilde{\tau}) \right] + O(e^4) \right\}, $$

$$ f_{2.5}(e) = -\frac{1}{8e^4} \left\{ (8g_1(e) - 19)g_{1.5}(e) - e^2[g_{1.5}(e)(16g_1(e) - 57) + 80 + 8g_{2.5}(e)] + O(e^4) \right\}, $$

$$ \begin{aligned} f_3(e) = \frac{1}{1024e^6} & \left\{ (8g_1(e) - 19)^3 - 2e^2 \left[ (8g_1(e) - 19)^2(12g_1(e) - 19) \right. \right. \\ & \left. \left. + 4(8g_1(e) - 19)(16g_2(e) + 197 - 4\chi^2(1 - 3\cos^2\tilde{\tau})) + 256g_{1.5}(e)^2\chi^2\cos^2\tilde{\tau} \right] \right. \\ & \left. + e^4 \left[ 4(8g_1(e) - 19)^2(6g_1(e) - 7) + 8(16g_1(e) - 57)(16g_2(e) + 197 - 4\chi^2(1 - 3\cos^2\tilde{\tau})) \right. \right. \\ & \left. \left. + (8g_1(e) - 19)(755 + \chi^2(32 - 120\sin^2\tilde{\tau} + 144\cos^2\tilde{\omega}\sin^2\tilde{\tau})) \right] \right. \\ & \left. + 1024g_3(e) - 79384 + 6528\chi^2 - 20608\chi^2\cos^2\tilde{\tau} + 1024g_{1.5}(e)^2\chi^2\cos^2\tilde{\tau} + O(e^6) \right\}. \end{aligned} \quad (4.3) $$

A necessary and sufficient condition that the $f_\alpha(e)$ be regular at $e = 0$ is the set of constraints

$$ \begin{aligned} g_1(e) &= \frac{19}{8}, & g_{1.5}(e) &= 0, & g_2(e) &= -\frac{197}{16} + \frac{1}{4}\chi^2(1 - 3\cos^2\tilde{\tau}), & g_{2.5}(e) &= -10, \\ g_3(e) &= \frac{9923}{128} - \frac{1}{8}\chi^2(51 - 161\cos^2\tilde{\tau}), & & & & \end{aligned} \quad (4.4) $$

resulting in the transformation

$$ \tilde{e} = e \left\{ 1 - \frac{3GM}{8c^2\tilde{p}}(19 + e^2) + \frac{1}{128}\left(\frac{GM}{c^2\tilde{p}}\right)^2 \left[5351 + 698e^2 + 23e^4 + 8\chi^2(20 - \sin^2\tilde{\tau}(39 - 18\cos^2\tilde{\omega}))\right] \right. \\ \left. + \frac{1}{4}\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2} \chi \cos\tilde{\tau} \left[7(10 + 27e^2) - 2(9 + 4e^2)\cos^2\tilde{\omega}\right] - \frac{1}{1024}\left(\frac{GM}{c^2\tilde{p}}\right)^3 \left[343065 + 107609e^2 + 4243e^4 + 91e^6 \right. \right. \\ \left. \left. + \chi^2(32(1571 + 571e^2) - 8(10791 + 3293e^2)\sin^2\tilde{\tau} - 128(47 + 40e^2)\cos^2\tilde{\omega} + 16(2137 + 843e^2)\cos^2\tilde{\omega}\sin^2\tilde{\tau})\right] \right\}. \quad (4.5) $$

With this transformation, $E$ now takes the form
$$c^{-2}E = -\frac{1}{2} \frac{GM(1-e^2)}{c^2 \tilde{p}} \left[ 1 - \frac{19}{4} \frac{GM}{c^2 \tilde{p}} + \frac{1}{8} \left( \frac{GM}{c^2 \tilde{p}} \right)^2 (197 - 4\chi^2(1-3\cos^2\tilde{t})) \right. \\ \left. + 20 \left( \frac{GM}{c^2 \tilde{p}} \right)^{5/2} \chi \cos\tilde{t} - \frac{1}{64} \left( \frac{GM}{c^2 \tilde{p}} \right)^3 (9923 + 16\chi^2(110 - 161\sin^2\tilde{t})) \right]. \quad (4.6)$$

As expected, in the limit $e \to 0$, E agrees with Eq. (3.24a) in the limit $\tilde{e} \to 0$.

In addition, there is no longer any $\tilde{\omega}$ dependence in the 2PN terms in E, thus curing the “quadrupole conundrum” discussed in Sec. III C. This can be understood as follows: in the conservative dynamics, $\tilde{p}$ and $\tilde{e}$ are constant through 2.5PN order, but vary at 3PN order [see Eqs. (3.22a) and (3.23b)]. However, $\tilde{\omega}$ appears in the 2PN term in the transformation between $\tilde{e}$ and $e$ in Eq. (4.5). Thus through 3PN order

$$\frac{d\tilde{e}}{d\theta} = \frac{de}{d\theta} - \frac{9}{4} e \left(\frac{GM}{c^2 \tilde{p}}\right)^2 \chi^2 \sin^2\tilde{t} \sin\tilde{\omega} \cos\tilde{\omega} \frac{d\tilde{\omega}}{d\theta}. \quad (4.7)$$

Substituting Eq. (3.23b) for $d\tilde{e}/d\theta$ along with $d\tilde{\omega}/d\theta = 3GM/c^2\tilde{p}$, we obtain

$$\frac{de}{d\theta} = 3e \left(\frac{GM}{c^2 \tilde{p}}\right)^3 \chi^2 \sin^2\tilde{t} \sin\tilde{\omega} \cos\tilde{\omega}. \quad (4.8)$$

Combining this with Eq. (3.22a) for $d\tilde{p}/d\theta$, which is unaffected to this order by the change in $e$, we find that

$$\frac{d\tilde{a}}{d\theta} \equiv \frac{d}{d\theta} \left( \frac{\tilde{p}}{(1 - e^2)} \right) = O\left( \frac{GM}{c^2 \tilde{p}} \right)^{7/2}. \quad (4.9)$$

In other words, 3PN variations in $\tilde{a}$ are no longer needed to compensate for variations in E induced by the advance of the pericenter, since that orbit element no longer appears in E.

The transformation of $e$ simultaneously cures the bad behavior of $P$, which now takes the form

$$P = 2\pi \left( \frac{\tilde{p}^3}{GM(1-e^2)^3} \right)^{1/2} \left\{ 1 + \frac{3}{8} \frac{GM}{c^2 \tilde{p}} (16-5e^2) + 6 \left( \frac{GM}{c^2 \tilde{p}} \right)^{3/2} \chi \cos \tilde{t} \right. \\ \left. - \frac{3}{128} \left( \frac{GM}{c^2 \tilde{p}} \right)^2 [448 - 88e^2 + 35e^4 - 320(1-e^2)^{3/2} - 64\chi^2(1-4\cos^2\tilde{t})] \right\}. \quad (4.10)$$

It also cures the behavior of the orbit shown in Eq. (3.21); in terms of $e$, $r^{-1}$ now factorizes into

$$r^{-1} = \frac{1 + \cos f}{\tilde{p}} \left[ 1 - 2 \frac{GM}{c^2 \tilde{p}} - \left( \frac{GM}{c^2 \tilde{p}} \right)^{3/2} \chi \cos \tilde{t} (6 - 2 \cos f - \sin f \sin 2\phi) + \dots \right]. \quad (4.11)$$

In terms of the new eccentricity, $C$ and $L_e$ become

$$L_e = \sqrt{m\tilde{p}} \left\{ \cos \tilde{t} \left[ 1 + \frac{GM}{2c^2\tilde{p}}(7+e^2) - \frac{1}{8}\left(\frac{GM}{c^2\tilde{p}}\right)^2(37+39e^2) + \frac{1}{16}\left(\frac{GM}{c^2\tilde{p}}\right)^3(411+1043e^2+18e^4 \right. \right. \\ \left. + 2\chi^2[92+132e^2-(100+177e^2)\sin^2\tilde{t}-32e^2\cos^2\tilde{\omega}+120e^2\cos^2\tilde{\omega}\sin^2\tilde{t}]) \right] \\ \left. - \chi \left(\frac{GM}{c^2\tilde{p}}\right)^{3/2} \left[ (1+\cos^2\tilde{t}) + \frac{GM}{4c^2\tilde{p}}(20+48e^2-(30+41e^2)\sin^2\tilde{t}-8e^2\cos^2\tilde{\omega}+10e^2\cos^2\tilde{\omega}\sin^2\tilde{t}) \right] \right\}, \quad (4.12a)$$

$$C = m\tilde{p} \left\{ 1 + \frac{GM}{c^2\tilde{p}}(7+e^2) - 4\left(\frac{GM}{c^2\tilde{p}}\right)^{3/2}\chi\cos\tilde{t} + \frac{1}{4}\left(\frac{GM}{c^2\tilde{p}}\right)^2[12-25e^2+e^4+2\chi^2\sin^2\tilde{t}(1-2e^2\cos^2\tilde{\omega})] \right. \\ \left. - 2\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2}\chi\cos\tilde{t}(12+13e^2-2e^2\cos^2\tilde{\omega}) + \frac{1}{8}\left(\frac{GM}{c^2\tilde{p}}\right)^3[152+733e^2-21e^4 \right. \\ \left. + 2\chi^2(108+132e^2-(140+193e^2)\sin^2\tilde{t}-32e^2\cos^2\tilde{\omega}+(115-e^2)e^2\cos^2\tilde{\omega}\sin^2\tilde{t})] \right\}. \quad (4.12b)$$

The final step in transforming to the new eccentricity is to determine the evolution induced by radiation reaction. We invert Eq. (4.5) to obtain $e$ as a function of $\tilde{e}$, $\tilde{p}$, $\tilde{t}$ and $\tilde{\omega}$, find $de/d\theta$ by inserting the radiation-reaction equations (3.27)
for the $d\tilde{X}_a/d\theta$, and then convert back to the new $e$. In the remaining $d\tilde{X}_a/d\theta_{RR}$ expressions, we convert from $\tilde{e}$ to $e$. The results are

$$ \frac{d\tilde{p}}{d\theta_{RR}} = -\frac{8}{5}\eta\tilde{p}\left(\frac{GM}{c^2\tilde{p}}\right)^{5/2}(8+7e^2) + \frac{1}{210}\eta\tilde{p}\left(\frac{GM}{c^2\tilde{p}}\right)^{7/2}(22072+27452e^2+281e^4) \\ +\frac{2}{15}\eta\tilde{p}\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\cos\tilde{\iota}(968+2280e^2+297e^4) \\ -\frac{1}{810}\eta\tilde{p}\left(\frac{GM}{c^2\tilde{p}}\right)^{9/2}(590900+941316e^2-100860e^4-4383e^6), \qquad (4.13a) $$

$$ \frac{de}{d\theta_{RR}} = -\frac{1}{15}\eta e \left(\frac{GM}{c^2\tilde{p}}\right)^{5/2}(304+121e^2) + \frac{1}{840}\eta e \left(\frac{GM}{c^2\tilde{p}}\right)^{7/2}(221000+120086e^2+1277e^4) \\ +\frac{1}{30}\eta e \left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\cos\tilde{\iota}(9400+10548e^2+789e^4) \\ -\frac{1}{15120}\eta e \left(\frac{GM}{c^2\tilde{p}}\right)^{9/2}(39598064+26131872e^2-1139399e^4-150795e^6), \qquad (4.13b) $$

$$ \frac{d\tilde{\iota}}{d\theta_{RR}} = -\frac{1}{30}\eta\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\sin\tilde{\iota}[88 - 1008e^2 - 285e^2 + 12e^2(20 + 13e^2)\cos(2\tilde{\omega})], \qquad (4.13c) $$

$$ \frac{d\tilde{\omega}}{d\theta_{RR}} = \frac{2}{5}\eta\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi\cos\tilde{\iota}e^2(20 + 13e^2)\sin(2\tilde{\omega}), \qquad (4.13d) $$

$$ \frac{d\tilde{\Omega}}{d\theta_{RR}} = -\frac{2}{5}\eta\left(\frac{GM}{c^2\tilde{p}}\right)^4\chi e^2(20 + 13e^2)\sin(2\tilde{\omega}). \qquad (4.13e) $$

V. NUMERICAL EVOLUTIONS

We now wish to apply the analytic results obtained in Secs. II–IV to the long-term evolution of highly eccentric orbits around a massive spinning black hole. For a $10^6 M_\odot$ black hole, these would be orbits with semimajor axis $a \sim$ milliparsecs, and $1 - e \sim 10^{-5}$. Since $GM/c^2$ for the black hole is of order 50 nanoparsecs, this corresponds to an initial semilatus rectum $p \sim 2a(1 - e) \sim 20(GM/c^2)$. Since our PN expansion is in powers of $GM/c^2\tilde{p}$, this is a regime where we might hope to obtain reasonable results for the orbital evolution.

We first average the secular evolution equations (4.13) over $\tilde{\omega}$; since the pericenter advances on a time scale that is short compared to the radiation-reaction time scale, this is a reasonable approximation. Since the dependence on $\tilde{\omega}$ occurs only in the highest-order terms in Eqs. (4.13), we do not need to carry out an additional two-scale analysis, and can average over $\tilde{\omega}$ holding the other elements fixed. [Notice from Eqs. (3.22) and (4.8) that the 3PN conservative variations in $\tilde{p}$, $e$, and $\tilde{\iota}$ also average to zero over a pericenter time scale.] We then scale $\tilde{p}$ by $GM/c^2$, defining the dimensionless variables $e \equiv GM/c^2\tilde{p}_i$ and $x \equiv \tilde{p}/\tilde{p}_i$, where $\tilde{p}_i$ is the initial value of $\tilde{p}$. In terms of these variables, the evolution equations become (henceforth, we drop the tildes on all our orbital variables)

$$ \frac{dx}{d\theta_{RR}} = -\frac{8}{5}\eta e^{5/2}x^{-3/2}(8+7e^2) + \frac{1}{210}\eta e^{7/2}x^{-5/2}(22072+27452e^2+281e^4) \\ +\frac{2}{15}\eta e^4x^3\chi\cos\iota(968+2280e^2+297e^4) - \frac{1}{810}\eta e^{9/2}x^{-7/2}(590900+941316e^2-100860e^4-4383e^6), \qquad (5.1a) $$

$$ \frac{de}{d\theta_{RR}} = -\frac{1}{15}\eta e e^{5/2} x^{-5/2}(304 + 121e^2) + \frac{1}{840}\eta e e^{7/2} x^{-7/2}(221000 + 120086e^2 + 1277e^4) \\ +\frac{1}{30}\eta e e^4 x^{-4}\chi \cos\iota(9400 + 10548e^2 + 789e^4) - \frac{1}{15120}\eta e e^{9/2} x^{-9/2}(39598064 + 26131872e^2 \\ -1139399e^4 - 150795e^6), \qquad (5.1b) $$
$$ \frac{dt}{d\theta_{RR}} = -\frac{1}{30}\eta\epsilon^4 x^{-4}\chi \sin\iota(88 - 1008e^2 - 285e^2), \quad (5.1c) $$

$$ \frac{d\tilde{\omega}}{d\theta_{RR}} = 0, \qquad (5.1d) $$

$$ \frac{d\tilde{\Omega}}{d\theta_{RR}} = 0. \qquad (5.1e) $$

Notice that, when the initial $p$ is scaled in units of $GM/c^2$, and with $x_i = 1$ the orbital evolution as a function of phase $\theta$ depends only on $\epsilon$, $e_i$ and $\iota_i$ and is independent of the mass of the black hole. We shall see that the only place where the black hole mass plays a role is in setting the conversion from the orbital phase to time.

In addition, the effect of the spin-orbit radiation reaction on the inclination $\iota$ is so small that the inclination has been found to be constant in all our numerical solutions. Henceforth we shall drop Eq. (5.1c), and treat the inclination as strictly constant.

### A. Capture by the black hole

As the body inspirals toward the black hole, its angular momentum decreases to a point whereupon, even in the absence of additional radiation reaction, it finds itself in an orbit with no inner turning point, and it ultimately crosses the event horizon. In most previous work this critical angular momentum has been chosen to be $4GM/c$, corresponding to the angular momentum of a particle with zero energy at the innermost unstable peak of the effective potential in the Schwarzschild geometry. For rotating black holes, this is not a good approximation, because of the strong effects of frame dragging, among other things. In Ref. [17], we derived a general condition for capture by a rotating black hole, for a body with zero energy (or with relativistic energy $\tilde{E} = 1$), by finding the critical value of the Carter constant such that the effective radial potential for a geodesic simultaneously vanishes and has a vanishing radial derivative (with a suitable sign for the second radial derivative), corresponding to a turning point right at the unstable peak of the potential. Defining $L \equiv \sqrt{C}$, and $\cos \iota \equiv L_c/L$, the result was an eighth-order polynomial equation for the critical value $L_c$ as a function of $\iota$ and the Kerr parameter $\chi$. The numerical results are plotted in Fig. 1 of [17]. We then found an analytic approximation to the curves given by

$$ L_c = \frac{2GM}{c} \left[ 1 + \sqrt{1 - \chi \cos \iota - \frac{1}{8} \chi^2 \sin^2 \iota F(\chi, \cos \iota)} \right], \quad (5.2) $$

which has the correct behavior both when $\chi = 0$ (Schwarzschild) and when $\iota = 0$ or $\iota = \pi$ (equatorial orbits) for arbitrary $\chi$. The function $F$ is a power series in $\chi$:

$$ F(\chi, \cos \iota) = 1 + \frac{1}{2}\chi \cos \iota + \frac{1}{64}\chi^2(7 + 13\cos^2 \iota) + \frac{1}{128}\chi^3 \cos \iota (23 + 5\cos^2 \iota) + O(\chi^4). \quad (5.3) $$

This series solution agrees with the numerical solutions of the critical condition to better than 0.5 percent for $0 \le \chi \le 0.9$; and to better than 5 percent for $0.9 \le \chi \le 0.99$.

In our numerical evolutions, we will impose the condition $\sqrt{C} = L_c$, where $C$ is given by Eq. (4.12b) and $L_c$ is given by Eq. (5.2). However, in the capture condition (5.2) the inclination was defined by $\cos \iota \equiv L_c/L$; from Eqs. (4.12a) and (4.12b), we see that

$$ \cos \iota = \cos \tilde{\iota} - \left(\frac{GM}{c^2p}\right)^{3/2} \chi \sin^2 \tilde{\iota} + O\left(\left(\frac{GM}{c^2p}\right)^2\right). \quad (5.4) $$

The difference between $\iota$ and $\tilde{\iota}$ is a maximum for polar orbits, but then we expect $GM/c^2p$ to be less than 1/16, so we expect the difference between the two measures of inclination to be only a few percent. We will ignore that difference.

Using the expression for $C$ to 2PN order and converting to our variables $x$ and $\epsilon$, we obtain the capture condition

$$ x^{1/2} \left[ 1 + \frac{\epsilon}{2x}(7 + e^2) - 2 \left(\frac{\epsilon}{x}\right)^{3/2} \chi \cos \iota - \frac{1}{8} \left(\frac{\epsilon}{x}\right)^2 (37 + 39e^2 - 2\chi^2(1 - e^2)\sin^2 \iota) \right] = \epsilon^{1/2} \left[ 1 + \sqrt{1 - \chi \cos \iota - \frac{1}{8} \chi^2 \sin^2 \iota F(\chi, \cos \iota)} \right]. \quad (5.5) $$

Solutions for $p_c/(GM/c^2) = x_c/\epsilon$ as a function of $\iota$ and $\chi$ are shown in Fig. 4. Since the orbits will have circularized to varying degrees by the time of plunge, we have chosen

<figure><img src="image_6.png" /><figcaption>FIG. 4. Critical values of $p$ for capture by a rotating black hole as a function of the Kerr parameter $\chi$, for various orbital inclinations. Orbital energy $E$ is assumed to be zero.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 5. Evolution of semilatus rectum *p* and eccentricity *e* vs number of orbits, for three initial values of *p*; in all cases, *e*<sub>i</sub> = 0.999 and *χ* = 1.</figcaption></figure>

*e* = 0 in that figure; for *e* = 0.5, for example, the curves are virtually the same.

We shall use the capture condition (5.5) to terminate
the numerical evolutions of *x* and *e*. An important caveat is
that this condition strictly applies only for orbits with
*E* = 0; however, since *E* ~ *GM*/*p*, we might expect
the error in our capture condition to be *O*(*GM*/*c*<sup>2</sup>*p*),
smaller for retrograde orbits, larger for prograde orbits.
Improvement of the capture criterion for *E* ≠ 0 is a topic for
further work.

B. Evolution in terms of orbital phase

We have integrated Eqs. (5.1a) and (5.1b) numerically, for initial values of *p* ranging from 10 to 100 *GM/c*², for 1 − *e*ᵢ ranging from 10⁻² to 10⁻⁵, and for τᵢ ranging from 0 (equatorial prograde) to π (equatorial retrograde). We choose η = 5 × 10⁻⁵, corresponding to a 50 *M*☉ object orbiting a 10⁶ *M*☉ black hole. We will frequently choose χ = 1, in order to maximize the effects due to the spin of the black hole.

Figure 5 shows the evolution of *p*/(*GM/c*²) and *e* as a
function of the number of orbits, for initial values
*p*ᵢ/(*GM/c*²) 100, 50 and 20, and for *e*ᵢ = 0.999. We find,
not surprisingly, that, for initial values of *e* very near unity,
the evolution with phase is independent of *e*ᵢ. Retrograde
orbits are seen to decay faster than prograde orbits, a
phenomenon that is clear from the spin-orbit terms in
Eqs. (5.1a) and (5.1b), and that has also been observed in
numerical relativity simulations of inspiraling spinning
black holes. In some mysterious way, this may be related
to the cosmic censorship hypothesis; since any residual
angular momentum that the particle carries into the black

hole increases the hole’s angular momentum if the inspiral
orbit is prograde, the particle must inspiral for longer to
ensure that enough angular momentum is radiated away
to preclude creating a black hole with $\chi > 1$ (a naked
singularity). No such issue is present for retrograde orbits.
Because of the shorter inspiral time, retrograde orbits are
also seen to plunge from orbits with larger eccentricity than
prograde orbits. In all cases, the larger the initial $p$, the
longer the inspiral, and hence the smaller the final eccen-
tricity. The residual values of $e$ are plotted against $p_i$ for
various inclinations in the left panel of Fig. 6. The right
panel of Fig. 6 shows the number of orbits to plunge vs the
initial $p$, for $e_i = 0.999$. The curves for $i = 135^\circ$ and $180^\circ$
drop precipitously to zero at values of $p_i$ that are already at
the critical value for immediate plunge.

C. Inclined orbits and the late-time flux
of gravitational waves

Here we study the effect of orbital inclination on the flux
and frequency of gravitational waves near the end point of
the orbital evolution. It is not our intent to provide accurate
predictions for the flux itself—in the highly relativistic
regime near the onset of plunge, this may be pushing the
PN approximation beyond its regime of validity. However,
our PN results may give useful insights into the dependence
of the gravitational waves on the inclination of the orbit.
To investigate this we estimate the flux as

$$
\mathcal{F} = \frac{2\pi}{P} \left( \frac{dE}{dp} \frac{dp}{d\theta_{RR}} + \frac{dE}{de} \frac{de}{d\theta_{RR}} \right), \quad (5.6)
$$

where *E* is given by Eq. (4.6) multiplied by η*M*, *P* is given
by (4.10), and *dp/dθ*<sub>RR</sub> and *de/dθ*<sub>RR</sub> are given by (4.13a)
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6. Left panel: Orbital eccentricity at time of plunge vs initial $p$. Right panel: Number of orbits vs initial $p$. In both cases $e_i = 0.999$ and $\chi = 1$.</figcaption></figure>

and (4.13b). Since the radiation-reaction expressions do not include $\chi^2$ terms, we set $\chi^2 = 0$ in the expressions for E and P. The result is

$$ \mathcal{F} = C^{-1} \frac{32\eta c^5}{5 G} \left(\frac{GM}{c^2 p}\right)^5 (1 - e^2)^{3/2} \left\{ \left(1 + \frac{73}{24}e^2 + \frac{37}{96}e^4\right) - \frac{GM}{c^2 p} \left(\frac{95216 + 306240e^2 + 53242e^4 + 715e^6}{5376}\right) \right. \\ \left. - \left(\frac{GM}{c^2 p}\right)^{3/2} \chi \cos \iota \left(\frac{1936 + 12024e^2 + 6582e^4 + 195e^6}{192}\right) \right. \\ \left. + \left(\frac{GM}{c^2 p}\right)^2 \left(\frac{121274560 + 421538216e^2 + 84768510e^4 - 1355193e^6 - 659322e^8}{580608}\right) \right\}, \quad (5.7) $$

where

$$ C = 1 + \frac{3GM}{8c^2p}(16 - 5e^2) + 6\left(\frac{GM}{c^2p}\right)^{3/2}\chi\cos\iota - \frac{3}{128}\left(\frac{GM}{c^2p}\right)^2[448 - 88e^2 + 35e^4 - 320(1 - e^2)^{3/2}]. \quad (5.8) $$

We estimate the gravitational-wave frequency as twice the orbital frequency, or $\omega_{GW} = 4\pi/P$. We then evaluate the flux and frequency at the point of plunge for various orbital evolutions. Both quantities depend almost entirely on the value of $p$ at plunge and on the orbital inclination and $\chi$. Because the eccentricity is generally less than 0.5 in all cases, the dependence on $e$ is very weak, and therefore largely independent of the initial values of $p$ and $e$. The top panel of Fig. 2 shows the flux for $p_i = 100GM/c^2$ and $e_i = 0.999$, as a function of $\chi$, for a range of inclinations, all normalized to the flux for $\iota = 0$, i.e. for equatorial prograde orbits. For the GW frequency, we normalize to the frequency for $\chi = 0$, so that the true GW frequency is given by

$$ \omega_{GW} = 0.057 \left( \frac{10^6 M_{\odot}}{M} \right) g(\chi, \iota) \text{ Hz,} \quad (5.9) $$

where $g(\chi, \iota)$ is the function plotted in the right panel of Fig. 2.

The strong suppression of the energy flux for highly inclined and retrograde orbits is mostly due to the fact that inclined orbits plunge from greater distances than do equatorial prograde orbits. For example, for $\chi = 0.6$, the factor of 140 suppression of the flux between 0° and 180° inclination is accounted for within a factor of 2 by the factor of 2.3 between the values of $p$ at plunge for the two types of orbit (see Fig. 4), raised to the fifth power. We have not plotted the corresponding curve for $\chi = 1$; while the PN expansion should be a reasonable approximation for the flux for retrograde orbits, with $p$ between 8 and $15GM/c^2$, it is not likely to be reasonable for the prograde equatorial orbit, with $p \sim 3GM/c^2$.

### D. Evolution in terms of time

We now address the evolution of our orbits with respect to time. We assume that, for each orbit element $\tilde{X}_\alpha$, we can approximate
<figure><img src="image_10.png" /><figcaption>FIG. 7. Left: Time to plunge in years vs initial $p$ for a $10^6 M_\odot$ black hole; $e_i = 0.999$ and $\chi = 1$. Right: Comparison between calculated plunge time and analytic fit of Eq. (1.1); for $p$ from 20 to $200GM/c^2$, and $\chi = 1$.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

$$ \frac{d\tilde{X}_\alpha}{dt} = \frac{2\pi}{P} \frac{d\tilde{X}_\alpha}{d\theta}. \qquad (5.10) $$

$$ \frac{dx}{dt} = -24\eta \left(\frac{c^3}{GM}\right) \left(\frac{GM}{c^2 p_i}\right)^4 \frac{(1 - e_i^2 x^{85/36})^{3/2}}{x^3}. \quad (5.13) $$

This is not strictly the same as having solved the original Lagrange planetary equations (3.5) in terms of time directly, since time-averaged orbit elements are not the same as angle-averaged elements, beyond first order in perturbation theory. But since the differences occur at higher PN orders, Eq. (5.10) will serve our present purposes. Figure 3 shows the evolution of $p$ and $e$ vs time for various initial values of $p$, for a $10^6 M_\odot$ black hole, and the left panel of Fig. 7 shows the time to plunge as a function of $p_i$ for various inclinations. Because the initial periods for these highly eccentric orbits are so long for a given $p$, a relatively small number of orbits consumes a great deal of time, so that nothing dramatic happens until very late (compare Figs. 3 and 5).

Because of this fact, it is possible to derive a simple analytic formula that encapsulates the time to plunge to reasonable accuracy, and then to refine it empirically to obtain a formula that also accounts for the differences due to spin-orbit effects. We first note from Eqs. (4.13a) and (4.13b) at leading quadrupole order, and for $e \approx 1$, that

$$ \frac{dp}{de} \approx \frac{72}{85} \frac{p}{e}, \qquad (5.11) $$

so that $e \approx e_i(p/p_i)^{85/72}$, as long as $e \sim 1$. Then, again to quadrupole order,

$$ \frac{dp}{dt} \approx \frac{2\pi}{P} \frac{dp}{d\theta} = -24\eta c \left(\frac{GM}{c^2 p}\right)^3 (1 - e^2)^{3/2}, \quad (5.12) $$

where we have used the Newtonian approximation for $P$. Substituting $x = p/p_i$ and $e = e_i x^{85/72}$, we obtain

The time to plunge is then given by inverting Eq. (5.13) and integrating over $x$ from $x_{\text{plunge}}$ to unity. Since $x_{\text{plunge}} \ll 1$ and since the integrand varies as $\sim x^3$ for small $x$, we can safely integrate over $x$ from 0 to 1, to obtain

$$ T_{\text{plunge}} = \frac{1}{96\eta} \left(\frac{GM}{c^3}\right) e^{-4} G(e_i), \quad (5.14) $$

where $e = GM/c^2 p_i$, and

$$ \begin{aligned} G(e_i) &= 4 \int_0^1 \frac{x^3 dx}{(1 - e_i^2 x^{85/36})^{3/2}} \\ &= {}_2F_1\left(\frac{3}{2}, \frac{144}{85}; \frac{229}{85}; e_i^2\right) \\ &= -\frac{11328}{935} \sqrt{\pi} \frac{\Gamma(59/85)}{\Gamma(33/70)} {}_2F_1\left(\frac{3}{2}, \frac{144}{85}; \frac{3}{2}; Z\right) \\ &\quad + \frac{288}{85} Z^{-1/2} {}_2F_1\left(1, \frac{203}{170}; \frac{1}{2}; Z\right), \end{aligned} \quad (5.15) $$

where ${}_2F_1(a, b; c; z)$ is the hypergeometric function and $Z = 1 - e_i^2$. Thus for a chosen $e$ and $e_i$, $T_{\text{plunge}}$ scales linearly with the mass of the black hole. The final representation of $G(e_i)$ is the most useful for values of $e_i$ near unity; it can be expressed as a series expansion of the form

$$ G(e_i) \approx \frac{288/85}{\sqrt{1 - e_i^2}} - 5.925739 + 8.091903 \sqrt{1 - e_i^2} + \dots \qquad (5.16) $$
Including only the first term in Eq. (5.16), we obtain $T_{\text{plunge}} = (3/85\eta)(GM/c^3)e^{-4}(1 - e_i^2)^{-1/2}$, which is in agreement with the results of Peters [37] [see also Eq. (26) of [38]]. We will adjust the coefficients slightly from the nominal values provided by Eq. (5.16) to get a decent fit to the actual dependence of $T_{\text{plunge}}$ on $e_i$. Including PN corrections to Eq. (5.14) particularly to account for spin-orbit effects, and adjusting the power of $\epsilon$ slightly from 4 to 3.96, we obtain the fit given in Eq. (1.1). As shown in the right panel of Fig. 7, this gives a surprisingly good fit to the time to plunge, to better than a percent for $p_i$ ranging from 80 to many hundred $GM/c^2$, to a few percent for $p_i$ between 40 and $80GM/c^2$ and to better than 20 percent down to $p_i = 20GM/c^2$. For example, for $p_i = 100$ and $e_i = 0.999$, the difference in plunge time between $\tau = 0^\circ$ and $\tau = 180^\circ$ is about seven percent, while Eq. (1.1) agrees with the calculated time for each inclination to 0.3 percent. The agreement as a function of $e_i$ for a given $p_i$ is at the level of fractions of a percent over values of $e_i$ ranging from 0.99 to 0.99999 that imply a range of over 4 orders of magnitude in orbital period. We have not attempted to tweak the formula to improve the fit for smaller values of $p_i$, since the orbits in this regime are already sufficiently relativistic that the PN approximation is becoming less accurate anyway.

Amaro-Seoane et al. [39] (hereafter ASF) have also endeavored to analyze the effects of black hole spin on the evolution of orbits. They employed a hybrid approach that combined 2PN results for the conservative orbital motion, black hole perturbation theory results for

gravitational radiation, and an adiabatic assumption that corrected orbits in response to the changes in $E$, $L_e$ and $C$ [5,6,40]. In particular, ASF reported results for the time to plunge (Table 1 of [39]) for a range of black hole masses and spins, for a range of initial semimajor axes and eccentricities, and for inclinations of $0^\circ$ and $57.3^\circ$ (one radian). As they pointed out, all the orbits discussed in that table actually correspond to the single choice $p_i = 8GM/c^2$, corresponding to the critical value for plunge for a Schwarzschild black hole, but above the critical value for prograde orbits around a Kerr black hole (see Fig. 4, where our value for the critical $p$ for Schwarzschild is closer to $9GM/c^2$). With this common choice of initial condition for $p$, then for all cases with the same values of $e_i$, $\tau_i$ and $\chi$, the time to plunge should scale with black hole mass [see Eq. (5.14)]. Table III shows a selection of results comparing our computation of time to plunge with those from Table 1 of ASF. The first thing to notice is that the scaling with mass is not apparent in the results of ASF. In addition, our times to plunge are consistently shorter than those of ASF and the dependence on the angle of inclination is stronger in our case. The starred entries in Table III are cases where the initial orbits are already below our plunge criterion of Eq. (5.5), and thus the plunge time is zero. We must acknowledge that, because $p$ in our calculations of $T_{\text{plunge}}$ is evolving from $9GM/c^2$ to around $3GM/c^2$, the PN approximation is being pushed up to or beyond its limit of validity, so we do not wish to claim too much accuracy for our values of $T_{\text{plunge}}$ in

<table><caption>TABLE III. Time to plunge: a comparison.</caption><thead><tr><th rowspan="2">MBH mass (M<sub>☉</sub>)</th><th rowspan="2">MBH spin χ</th><th rowspan="2">Semimajor axis a<sub>i</sub> (pc)</th><th rowspan="2">Eccentricity 1 − e<sub>i</sub></th><th rowspan="2">Inclination τ (radians)</th><th colspan="2">Time to plunge (yr)</th></tr><tr><th>This work</th><th>Ref. [39]</th></tr></thead><tbody><tr><td>5.0 × 10<sup>4</sup></td><td>0.30</td><td>9.57 × 10<sup>−3</sup></td><td>10<sup>−6</sup></td><td>0.0</td><td>93</td><td>392</td></tr><tr><td>5.0 × 10<sup>4</sup></td><td>0.30</td><td>9.57 × 10<sup>−3</sup></td><td>10<sup>−6</sup></td><td>0.7</td><td>***</td><td>392</td></tr><tr><td>5.0 × 10<sup>5</sup></td><td>0.30</td><td>9.57 × 10<sup>−2</sup></td><td>10<sup>−6</sup></td><td>0.0</td><td>934</td><td>4.31 × 10<sup>4</sup></td></tr><tr><td>5.0 × 10<sup>5</sup></td><td>0.30</td><td>9.57 × 10<sup>−2</sup></td><td>10<sup>−6</sup></td><td>1.0</td><td>***</td><td>4.31 × 10<sup>4</sup></td></tr><tr><td>1.0 × 10<sup>7</sup></td><td>0.30</td><td>1.91</td><td>10<sup>−6</sup></td><td>0.0</td><td>1.87 × 10<sup>4</sup></td><td>4.02 × 10<sup>6</sup></td></tr><tr><td>1.0 × 10<sup>7</sup></td><td>0.30</td><td>1.91</td><td>10<sup>−6</sup></td><td>1.0</td><td>***</td><td>3.78 × 10<sup>6</sup></td></tr><tr><td>1.0 × 10<sup>6</sup></td><td>0.70</td><td>1.91 × 10<sup>−2</sup></td><td>10<sup>−5</sup></td><td>0.0</td><td>960</td><td>1.35 × 10<sup>4</sup></td></tr><tr><td>1.0 × 10<sup>6</sup></td><td>0.70</td><td>1.91 × 10<sup>−2</sup></td><td>10<sup>−5</sup></td><td>1.0</td><td>644</td><td>1.20 × 10<sup>4</sup></td></tr><tr><td>5.0 × 10<sup>6</sup></td><td>0.70</td><td>9.57 × 10<sup>−2</sup></td><td>10<sup>−5</sup></td><td>0.0</td><td>4.80 × 10<sup>3</sup></td><td>1.55 × 10<sup>5</sup></td></tr><tr><td>5.0 × 10<sup>6</sup></td><td>0.70</td><td>9.57 × 10<sup>−2</sup></td><td>10<sup>−5</sup></td><td>1.0</td><td>3.22 × 10<sup>3</sup></td><td>1.36 × 10<sup>5</sup></td></tr><tr><td>5.0 × 10<sup>7</sup></td><td>0.70</td><td>9.57 × 10<sup>−1</sup></td><td>10<sup>−5</sup></td><td>0.0</td><td>4.80 × 10<sup>4</sup></td><td>1.55 × 10<sup>7</sup></td></tr><tr><td>5.0 × 10<sup>7</sup></td><td>0.70</td><td>9.57 × 10<sup>−1</sup></td><td>10<sup>−5</sup></td><td>1.0</td><td>3.22 × 10<sup>4</sup></td><td>1.35 × 10<sup>7</sup></td></tr><tr><td>1.0 × 10<sup>6</sup></td><td>0.99</td><td>1.91 × 10<sup>−2</sup></td><td>10<sup>−5</sup></td><td>0.0</td><td>1.66 × 10<sup>3</sup></td><td>1.43 × 10<sup>4</sup></td></tr><tr><td>1.0 × 10<sup>6</sup></td><td>0.99</td><td>1.91 × 10<sup>−2</sup></td><td>10<sup>−5</sup></td><td>1.0</td><td>782</td><td>1.17 × 10<sup>4</sup></td></tr><tr><td>1.0 × 10<sup>7</sup></td><td>0.99</td><td>1.91 × 10<sup>−1</sup></td><td>10<sup>−5</sup></td><td>0.0</td><td>1.66 × 10<sup>4</sup></td><td>1.44 × 10<sup>6</sup></td></tr><tr><td>1.0 × 10<sup>7</sup></td><td>0.99</td><td>1.91 × 10<sup>−1</sup></td><td>10<sup>−5</sup></td><td>1.0</td><td>7.82 × 10<sup>3</sup></td><td>1.18 × 10<sup>6</sup></td></tr><tr><td>5.0 × 10<sup>7</sup></td><td>0.99</td><td>9.57 × 10<sup>−1</sup></td><td>10<sup>−5</sup></td><td>0.0</td><td>8.32 × 10<sup>4</sup></td><td>1.65 × 10<sup>7</sup></td></tr><tr><td>5.0 × 10<sup>7</sup></td><td>0.99</td><td>9.57 × 10<sup>−1</sup></td><td>10<sup>−5</sup></td><td>1.0</td><td>3.91 × 10<sup>4</sup></td><td>1.35 × 10<sup>7</sup></td></tr></tbody></table>
Table III. Nevertheless, the scaling with black hole mass and the dependence on inclination should be robust. Reconciling the differences between these two approaches is a subject for future research.

## VI. DISCUSSION AND CONCLUSIONS

We have used post-Newtonian theory to analyze the motion of stars orbiting rotating black holes in orbits that are initially very eccentric and with arbitrary inclinations relative to the black hole's equatorial plane. We incorporated conservative terms from the test body equations of motion in the Kerr geometry, including all spin effects through 3PN order, and radiation-reaction contributions to the equations of motion through 4.5PN order, including spin-orbit contributions, and carefully solved the evolution equations for the osculating orbit elements to the corresponding PN order using a two-time-scale approach. The orbits were terminated when the Carter constant associated with the orbit dropped below a critical value. We found that retrograde orbits terminate farther from the black hole and with larger residual orbital eccentricities than do prograde orbits, and consequently that the flux of gravitational radiation from the final stage of retrograde orbits can be substantially suppressed relative to the flux from the comparable prograde orbit. We also provided a number of results in forms that are “ready to use” in numerical simulations of star clusters orbiting a spinning central black hole, including two-body equations of motion and an analytic formula for the time to plunge that take into account the effect of the spin of the black hole.

A key question is how well the PN approximation can be trusted in the relativistic regimes we consider here. The inverse of the relativistic expansion parameter $GM/c^2p$ ranges from values as large as the hundreds when the star is injected into a high-eccentricity orbit from the surrounding cluster to as small as 3 when a prograde orbit terminates at the plunge point. Consequently, some caution is in order in interpreting these results. While the PN approximation has proven to be “unreasonably effective” in describing the dynamics of comparable-mass binary systems [16], its effectiveness is not so clear in the extreme mass-ratio limit, particularly in the absence of alternative approaches such as numerical relativity or self-force theory against which detailed comparisons could be made. It would be interesting to see if various resummations of the PN sequence, along the lines of the “effective one-body” (EOB) approach [41,42] might be useful for the eccentric orbits being considered here (see e.g. [43]).

## ACKNOWLEDGMENTS

This work was supported in part by the National Science Foundation, Grants No. PHY 13-06069 and

No. PHY 16-00188. We are grateful for the hospitality of the Institut d'Astrophysique de Paris where much of this work was carried out. We thank Tal Alexander, Anna Heffernan, David Merritt and Bernard Whiting for useful discussions. One of us (C. M. W.) thanks Alexandre le Tiec for permission to incorporate unpublished joint work into Appendix A.

## APPENDIX A: MOTION OF A PARTICLE AROUND A SPINNING BLACK HOLE IN HARMONIC COORDINATES

In this Appendix, we provide some of the details underlying the equations of motion used in this paper.

### 1. Test-body equations to 3PN order

The conservative equations of motion used in this paper are those of a test body orbiting a Kerr black hole, to 3PN order and in harmonic coordinates. We begin with the Kerr line element in standard Boyer-Lindquist coordinates [44], given by

$$ds^2 = -\left(1 - \frac{2mr_{\text{bl}}}{\rho^2}\right)c^2 dt_{\text{bl}}^2 + \frac{\rho^2}{\Delta} dr_{\text{bl}}^2 + \rho^2 d\theta_{\text{bl}}^2 \\ + \left(r_{\text{bl}}^2 + \frac{a^2}{c^2} + \frac{2mr_{\text{bl}}a^2\sin^2\theta_{\text{bl}}}{c^2\rho^2}\right)\sin^2\theta_{\text{bl}}d\phi_{\text{bl}}^2 \\ - \frac{4mar_{\text{bl}}}{\rho^2}\sin^2\theta_{\text{bl}}dt_{\text{bl}}d\phi_{\text{bl}}, \quad (\text{A1})$$

where $m \equiv GM/c^2$ is the geometrized mass and $a = GM\chi/c$ is the Kerr parameter, related to the spin angular momentum $S$ by $a \equiv S/M$, with $|\chi| \le 1$ required for the metric to describe a black hole; $\rho^2 = r_{\text{bl}}^2 + (a/c)^2\cos^2\theta_{\text{bl}}$, and $\Delta = r_{\text{bl}}^2 - 2mr_{\text{bl}} + (a/c)^2$. The transformation from Boyer-Lindquist coordinates to Cartesian-like harmonic coordinates satisfying $\partial_{\alpha}(\sqrt{-gg^{\alpha\beta}\partial_{\beta}x_{\text{H}}^{(y)}}) = 0$ is given by [45–47]

$$\begin{aligned} t_{\text{H}} &= t_{\text{bl}}, \\ x_{\text{H}} + iy_{\text{H}} &= (r_{\text{bl}} - m + ia/c)e^{i\psi} \sin \theta_{\text{bl}}, \\ z_{\text{H}} &= (r_{\text{bl}} - m) \cos \theta_{\text{bl}}, \end{aligned} \quad (\text{A2})$$

where

$$\psi = \phi_{\text{bl}} + \frac{a}{c(r_{+} - r_{-})} \ln \left| \frac{r_{\text{bl}} - r_{+}}{r_{\text{bl}} + r_{-}} \right| \quad (\text{A3})$$

and $r_{\pm} = m \pm \sqrt{m^2 - (a/c)^2}$. Inverting these transformations iteratively to the requisite PN order and inserting into Eq. (A1), we obtain the metric in harmonic coordinates to 3PN order:
$$
\begin{align}
g_{00} &= -\left(\frac{1 - GM/c^2 r}{1 + GM/c^2 r}\right) - \left(\frac{GM}{c^2 r}\right)^3 \chi^2[3(\mathbf{n} \cdot \mathbf{e})^2 - 1] \nonumber \\
&\quad + 2\left(\frac{GM}{c^2 r}\right)^4 \chi^2[4(\mathbf{n} \cdot \mathbf{e})^2 - 1] + O(\epsilon^5), \\
g_{0j} &= 2\frac{GM^2}{c^4 r^2} \chi(\mathbf{n} \times \mathbf{e})^j \left(1 - \frac{GM}{c^2 r}\right) + O(\epsilon^4), \\
g_{ij} &= \left(1 + \frac{GM}{c^2 r}\right)^2 (\delta_{ij} - n_i n_j) + \left(\frac{1 + GM/c^2 r}{1 - GM/c^2 r}\right) n_i n_j \nonumber \\
&\quad - \left(\frac{GM}{c^2 r}\right)^3 \chi^2[3(\mathbf{n} \cdot \mathbf{e})^2 - 1]\delta_{ij} \nonumber \\
&\quad + 2\left(\frac{GM}{c^2 r}\right)^3 \chi(\mathbf{n} \times \mathbf{e})^{(i}n^{j)} + O(\epsilon^4), \tag{A4}
\end{align}
$$

where **e** is a unit vector in the direction of the black hole spin. The first term in $g_{00}$ is the usual Schwarzschild piece in harmonic coordinates; the second term is the contribution of the body's quadrupole moment, and the third term is a 3PN cross term between the monopole and quadrupole terms combined with a (spin)$^2$ term. That the terms in the time-space component $g_{0j}$ are of order $c^{-4}$ and $c^{-6}$ rather than $c^{-3}$ and $c^{-5}$ is a reflection of the well-known fact that for black holes, effects linear in spin are effectively 1/2 PN order higher than they would be for normal slowly rotating bodies. In $g_{ij}$ the first two terms are again the standard Schwarzschild contribution, while the third is the quadrupole term. The final term may be dropped: it is a pure gauge

term, as can be checked by making the spatial coordinate transformation

$$x^j = \bar{x}^j + \frac{1}{3} \left( \frac{GM}{c^2 \bar{r}} \right)^3 \chi(\bar{x} \times \mathbf{e})^j, \quad (A5)$$

which maintains harmonic gauge, eliminates that term, and leaves the rest of the metric unchanged, to 3PN order.

The equations of motion for a test body can then be obtained from the geodesic equation, expressed using coordinate time as the parameter

$$\frac{d^2 x^i}{dt^2} + \left(\Gamma^j_{\beta\gamma} - \Gamma^0_{\beta\gamma} \frac{v^j}{c}\right) v^\beta v^\gamma = 0, \quad (A6)$$

where $v^\alpha \equiv (c, v^j)$, and $\Gamma^\alpha_{\beta\gamma}$ are the Christoffel symbols computed from the metric. The result through 3PN order is Eq. (2.3). Since the Kerr metric is stationary and axisymmetric, it admits the two conserved quantities of energy (per unit mass-energy) $\tilde{E}$ and the **e**-component of angular momentum per unit mass, $L_e$, given by

$$\tilde{E} \equiv -u_0/c, \quad L_e \equiv u_\phi, \quad (A7)$$

where $u^\alpha \equiv dx^\alpha/d\tau$ is the four-velocity of the test particle. It is straightforward to express these to 3PN order in harmonic coordinates using Eq. (A4). Defining $\tilde{E} \equiv 1 + E/c^2$, we obtain

$$
\begin{align}
E &= \frac{1}{2} v^2 - \frac{GM}{r} + \frac{1}{8c^2} \left( 3v^4 + 12v^2 \frac{GM}{r} + 4 \frac{G^2 M^2}{r^2} \right) + \frac{1}{c^4} \left( \frac{5}{16} v^6 + \frac{21}{8} v^4 \frac{GM}{r} + \frac{7}{4} v^2 \frac{G^2 M^2}{r^2} + \frac{1}{2} v^2 \frac{G^2 M^2}{r^2} - \frac{1}{2} \frac{G^3 M^3}{r^3} \right) \nonumber \\
  &\quad + \frac{1}{2c^4} \left( \frac{GM}{r} \right)^3 \chi^2[3(\mathbf{n} \cdot \mathbf{e})^2 - 1] - \frac{2}{c^5} \left( \frac{GM}{r} \right)^2 v^2 \chi \frac{L_z}{r} \nonumber \\
  &\quad + \frac{1}{c^6} \left( \frac{35}{128} v^8 + \frac{55}{16} v^6 \frac{GM}{r} + \frac{135}{16} v^4 \frac{G^2 M^2}{r^2} + \frac{3}{4} v^2 v^2 \frac{G^2 M^2}{r^2} + \frac{5}{4} v^2 \frac{G^3 M^3}{r^3} + \frac{3}{2} v^2 \frac{G^3 M^3}{r^3} + \frac{3}{8} \frac{G^4 M^4}{r^4} \right) \nonumber \\
  &\quad - \frac{1}{4c^6} \left( \frac{GM}{r} \right)^3 \chi^2 \left( 3v^2[3(\mathbf{n} \cdot \mathbf{e})^2 - 1] + 2 \frac{GM}{r} [5(\mathbf{n} \cdot \mathbf{e})^2 - 1] \right), \tag{A8}
\end{align}
$$

$$
\begin{align}
L_e = L_z & \left\{ 1 + \frac{1}{2c^2} \left( v^2 + 6 \frac{GM}{r} \right) + \frac{1}{8c^4} \left( 3v^4 + 28v^2 \frac{GM}{r} + 28 \frac{G^2 M^2}{r^2} \right) - \frac{2}{c^5} \left( \frac{GM}{r} \right)^2 \frac{L_z}{r} \chi \right. \nonumber \\
& \left. + \frac{1}{c^6} \left( \frac{5}{16} v^6 + \frac{33}{8} v^4 \frac{GM}{r} + \frac{45}{4} v^2 \frac{G^2 M^2}{r^2} + \frac{1}{2} v^2 \frac{G^2 M^2}{r^2} + \frac{5}{2} \frac{G^3 M^3}{r^3} \right) - \frac{3}{2c^6} \left( \frac{GM}{r} \right)^3 \chi^2[3(\mathbf{n} \cdot \mathbf{e})^2 - 1] \right\} \nonumber \\
& + \frac{2}{c^3} \frac{G^2 M^2}{r} \chi[(\mathbf{n} \cdot \mathbf{e})^2 - 1] \left( 1 + \frac{v^2}{2c^2} \right), \tag{A9}
\end{align}
$$

where $L_z \equiv (\mathbf{x} \times \mathbf{v}) \cdot \mathbf{e}$. The geometry admits an additional constant of the motion known as the Carter constant [48], given by $C_0 \equiv K_{\mu\nu} u^\mu u^\nu$, where $K_{\mu\nu}$ is the Killing tensor of the metric. However, we will find it more useful to employ a related “Carter” constant, given by

$$C \equiv C_0 - (a\tilde{E})^2 + 2a\tilde{E}L_e = \rho_{\text{bl}}^4(u^{\theta_{\text{bl}}})^2 + \sin^{-2}\theta_{\text{bl}}L_e^2 + a^2\cos^2\theta_{\text{bl}}(1 - \tilde{E}^2). \quad (A10)$$
This definition has the advantage that, in the Schwarzschild limit ($a \to 0$), $\sqrt{C}$ is equal to the conserved total angular momentum, and in the Newtonian limit, it is equal to $L \equiv |\mathbf{x} \times \mathbf{v}|$. Thus this version of C will have a more useful post-Newtonian expansion. Converting to harmonic coordinates using Eqs. (A2), inserting the expressions for $\tilde{E} = 1 + E/c^2$ and $L_e$, and using the fact that, in our harmonic coordinates, $(d\theta/dt)^2 = r^{-4}(L^2 - \sin^{-2}\theta L_z^2)$, we obtain, to 3PN order,

$$C = L^2 \left\{ \begin{aligned} & 1 + \frac{1}{c^2} \left( v^2 + 6 \frac{GM}{r} \right) + \frac{1}{c^4} \left( v^4 + 10v^2 \frac{GM}{r} + 16 \frac{G^2 M^2}{r^2} \right) \\ & + \frac{1}{c^4} \left( \frac{GM}{r} \right)^2 \chi^2 [2(\mathbf{n} \cdot \mathbf{e})^2 - 1] - \frac{4}{c^5} \left( \frac{GM}{r} \right)^2 \frac{L_z}{r} \chi \\ & + \frac{1}{c^6} \left( v^6 + 14v^4 \frac{GM}{r} + 47v^2 \frac{G^2 M^2}{r^2} + \dot{r}^2 \frac{G^2 M^2}{r^2} + 26 \frac{G^3 M^3}{r^3} \right) \\ & + \frac{1}{c^6} \left( \frac{GM}{r} \right)^2 \chi^2 \left( v^2 [2(\mathbf{n} \cdot \mathbf{e})^2 - 1] + 3 \frac{GM}{r} [(\mathbf{n} \cdot \mathbf{e})^2 - 1] \right) \end{aligned} \right\} \\ & - \frac{4}{c^3} \frac{G^2 M^2}{r} L_z \chi \left\{ 1 + \frac{1}{c^2} \left( v^2 + 3 \frac{GM}{r} \right) \right\} + \frac{1}{c^4} \left( \frac{GM}{r} \right)^2 L_z^2 \chi^2 \left\{ 1 + \frac{1}{c^2} \left( v^2 + 6 \frac{GM}{r} \right) \right\} \\ & - \frac{G^2 M^2}{c^4} \chi^2 (\mathbf{n} \cdot \mathbf{e})^2 \left\{ v^2 - 2 \frac{GM}{r} + \frac{1}{c^2} \left( v^4 + 2v^2 \frac{GM}{r} + 2 \frac{G^2 M^2}{r^2} \right) \right\} \\ & - \frac{2}{c^4} L \dot{r} \frac{G^2 M^2}{r} \chi^2 (\mathbf{\lambda} \cdot \mathbf{e})(\mathbf{n} \cdot \mathbf{e}) \left\{ 1 + \frac{1}{c^2} \left( v^2 + 6 \frac{GM}{r} \right) \right\} - \frac{4}{c^6} \frac{G^4 M^4}{r^2} \chi^2 [(\mathbf{n} \cdot \mathbf{e})^2 - 1], \end{aligned} \right. \quad (A11)$$

where $\lambda \equiv \hat{h} \times \mathbf{n}$, $\hat{h} = L/L$ and we use the fact that $r^2 \sin\theta\dot{\theta} = -h(\lambda \cdot \mathbf{e})$. In these expressions the terms of odd orders in $1/c$ arise from spin-orbit coupling, whose effects are always shifted by half a PN order relative to monopole or quadrupole contributions.

## 2. Radiation-reaction terms at 4PN spin orbit and 4.5 PN order

Here we display the coefficients in the 4PN spin orbit and 4.5PN contributions to radiation reaction in Eq. (2.4). The spin-orbit terms were derived from first principles in [20]; the coefficients are given by

$$ \begin{aligned} A_{4\text{SO}} &= 120v^2 + 280\dot{r}^2 + 453\frac{GM}{r}, \\ B_{4\text{SO}} &= 87v^2 - 675\dot{r}^2 - \frac{901}{3}\frac{GM}{r}, \\ C_{4\text{SO}} &= 48v^2 + 15\dot{r}^2 + 364\frac{GM}{r}, \\ D_{4\text{SO}} &= 31v^4 - 260v^2\dot{r}^2 + 245\dot{r}^4 - \frac{689}{3}v^2\frac{GM}{r} + 537\dot{r}^2\frac{GM}{r} + \frac{4}{3}\left(\frac{GM}{r}\right)^2. \end{aligned} \quad (A12) $$

Because we are dealing with a single spinning body at rest, these coefficients are independent of the “spin-supplementary condition” adopted.

The 4.5PN contributions were derived by Gopakumar et al. [21], extending a method developed by Iyer and Will [49,50]. The idea is to write down the most general expression for 4.5PN radiation-reaction terms for a binary system of point masses (in the center-of-mass frame), and to fix the coefficients in that expression by requiring that the equations yield the proper expressions for energy and angular momentum flux, accurate to 2PN order beyond the quadrupole approximation. The procedure does not fix the coefficients uniquely, but the remaining free parameters can be shown to be related to the freedom to make coordinate or gauge changes at the relevant 4.5PN order. At 4.5PN order, there are 12 degrees of freedom, parametrized by $\psi_1...\psi_9, \chi_6, \chi_8$ and $\chi_9$, in the notation of Ref. [21]. The resulting functions $E_{4.5}$ and $F_{4.5}$ in Eq. (2.4) are given by
$$ 
\begin{aligned}
E_{4.5} = & \left(\frac{779}{168} + 3\psi_2 - 3\chi_6\right) v^6 - \left(\frac{295}{84} + 5\psi_2 - 5\chi_6 - 5\psi_4 + 5\chi_8\right) v^4 \dot{r}^2 - 9\psi_7 \dot{r}^6 \\
& + \left(\frac{145}{6} - 7\psi_4 + 7\chi_8 + 7\psi_7\right) v^2 \dot{r}^4 + \left(\frac{6793}{84} - 2\psi_1 - 3\psi_2 + 3\chi_6 + 3\psi_6 - 3\chi_9\right) v^4 \frac{GM}{r} \\
& - \left(\frac{218401}{504} + 4\psi_2 + 5\psi_4 + 6\psi_6 - 5\psi_8 - 2\chi_6 - 5\chi_8 - 6\chi_9\right) v^2 \dot{r}^2 \frac{GM}{r} \\
& + \left(\frac{54161}{126} - 2\psi_4 - 7\psi_7 - 8\psi_8\right) \dot{r}^4 \frac{GM}{r} - \left(83 + 2\psi_3 + 3\psi_6 - 3\chi_9 - 3\psi_9\right) v^2 \frac{G^2 M^2}{r^2} \\
& + \left(\frac{83407}{252} - 2\psi_6 - 5\psi_8 - 7\psi_9\right) \dot{r}^2 \frac{G^2 M^2}{r^2} + \left(\frac{41297}{108} - 2\psi_5 - 3\psi_9\right) \frac{G^3 M^3}{r^3},
\end{aligned}
$$

$$ 
\begin{aligned}
F_{4.5} = & \left(\frac{417}{28} - \psi_1\right) v^6 - \left(\frac{380}{3} - 3\psi_1 + 3\chi_6\right) v^4 \dot{r}^2 - \left(\frac{485}{6} - 7\chi_8\right) \dot{r}^6 + \left(\frac{34445}{168} + 5\chi_6 - 5\chi_8\right) v^2 \dot{r}^4 \\
& + \left(\frac{1859}{56} + \psi_1 - \psi_3\right) v^4 \frac{GM}{r} - \left(\frac{16687}{42} - 4\psi_1 - 4\psi_3 - 3\chi_6 + 3\chi_9\right) v^2 \dot{r}^2 \frac{GM}{r} \\
& + \left(\frac{99499}{252} + 2\chi_6 + 5\chi_8 + 6\chi_9\right) \dot{r}^4 \frac{GM}{r} - \left(\frac{2967}{28} - \psi_3 + \psi_5\right) v^2 \frac{G^2 M^2}{r^2} \\
& + \left(\frac{3166}{21} + 2\psi_3 + 5\psi_5 + 3\chi_9\right) \dot{r}^2 \frac{G^2 M^2}{r^2} + \left(\frac{395929}{2268} + \psi_5\right) \frac{G^3 M^3}{r^3}.
\end{aligned}
\quad (A13) $$

It turns out, not surprisingly, that the arbitrary parameters completely drop out of the orbit-averaged equations for the orbital elements.

## APPENDIX B: TWO-SCALE ANALYSIS OF THE LAGRANGE PLANETARY EQUATIONS

Here we provide some of the details of the two-scale analysis that we used to find the equations for the secular evolutions of the orbit elements. The Lagrange planetary equations for the orbit elements $X_\alpha(t)$ take the general form

$$ \frac{dX_\alpha(t)}{dt} = \epsilon Q_\alpha(X_\beta(t), t), \quad (B1) $$

where $\alpha$ labels the orbit element, $\epsilon$ is a small parameter that characterizes the perturbation, and the $Q_\alpha(X_\beta(t), t)$ are functions of the orbit elements $X_\beta$ as well as explicit functions of $t$ (in practice, an angular variable such as the phase $\phi$ or the true anomaly $f$ is often a proxy for $t$). Those functions are assumed to be periodic with period $P$ or $2\pi$. To zeroth order in $\epsilon$, $dX_\alpha/dt = 0$, and thus the $X_\alpha$ are constants.

In general there are six orbital elements, including one related to the “time of pericenter passage,” but for many problems for which the actual time of events is not important, it is useful to eliminate time from the problem and to work in terms of an angular variable expressing the orbital phase. The possibilities include the phase $\phi$

measured from the ascending node, the true anomaly $f$ measured from the pericenter of the orbit, or the “eccentric anomaly” $u$, related to $f$ by the transformations

$$ 
\begin{aligned}
\cos u &= \frac{\cos f + e}{1 + e \cos f}, \\
\sin u &= \frac{\sqrt{1 - e^2} \sin f}{1 + e \cos f},
\end{aligned}
\quad (B2) $$

with $r = p/(1 + e \cos f) = a(1 - e \cos u)$, where $a = p(1 - e^2)$. The transformations between these variables and time $t$ are given by

$$ 
\begin{aligned}
\frac{d\phi}{dt} &= \frac{h}{r^2} - \cos \iota \frac{d\Omega}{dt}, \\
\frac{df}{dt} &= \frac{h}{r^2} - \frac{d\omega}{dt} - \cos \iota \frac{d\Omega}{dt}, \\
\frac{du}{dt} &= \frac{\sqrt{1 - e^2}}{1 + e \cos f} \frac{df}{dt}.
\end{aligned}
\quad (B3) $$

Another possibility is the “mean anomaly,” $\ell \equiv 2\pi t/P$, where $P$ is the orbital period. After one of these transformations, we arrive at five orbit element equations in the generic form

$$ \frac{dX_\alpha(\phi)}{d\phi} = \epsilon Q_\alpha(X_\beta(\phi), \phi), \quad (B4) $$
where $Q_\alpha(X_\beta(\phi), \phi) = (dt/d\phi)Q_\alpha(X_\beta(t), t)$, where we use $\phi$ to represent generically the phase variable $\{\phi, f, u, \ell\}$ being used.

As outlined in Sec. III B, we define the long-time-scale variable $\theta \equiv \epsilon\phi$, write the derivative with respect to $\phi$ as $d/d\phi \equiv \epsilon\partial/\partial\theta + \partial/\partial\phi$ and define

$$X_{\alpha}(\theta, \phi) \equiv \tilde{X}_{\alpha}(\theta) + \epsilon Y_{\alpha}(\tilde{X}_{\beta}(\theta), \phi), \quad (B5)$$

where $\tilde{X}_\alpha(\theta)$ is the average of $X_\alpha$ over $\phi$, and $Y_\alpha$ is the average-free part, where the average and average-free parts are defined by

$$\langle A \rangle \equiv \frac{1}{2\pi} \int_0^{2\pi} A(\theta, \phi) d\phi, \quad \mathcal{A}\mathcal{F}(A) \equiv A(\theta, \phi) - \langle A \rangle, \quad (B6)$$

where the integrals are carried out holding $\theta$ fixed.

Substituting our definition of $X_\alpha$ into Eq. (B4), and taking the average and average-free parts, we obtain the two main equations of the procedure

$$\frac{d\tilde{X}_\alpha}{d\theta} = \langle Q_\alpha(\tilde{X}_\beta + \epsilon Y_\beta, \phi) \rangle, \quad (B7a)$$

$$\frac{\partial Y_{\alpha}}{\partial \phi} = \mathcal{A}\mathcal{F}(Q_{\alpha}(\tilde{X}_{\beta} + \epsilon Y_{\beta}, \phi)) - \epsilon \frac{\partial Y_{\alpha}}{\partial \tilde{X}_{\gamma}} \frac{d\tilde{X}_{\gamma}}{d\theta}. \quad (B7b)$$

Note that, by virtue of our assumption that $\theta$ and $\phi$ are independent, $\partial Y_\alpha/\partial \tilde{X}_\gamma$ is automatically average free.

In integrating Eq. (B7b), it is useful to consider the general equation

$$\frac{\partial A}{\partial \phi} = \mathcal{A}\mathcal{F}(B), \quad (B8)$$

for some functions A and B, where B is periodic with period $2\pi$. The solution is

$$A = \int_{0}^{\phi} \mathcal{A}\mathcal{F}(B) d\phi' + C, \quad (B9)$$

where C is a constant of integration, fixed by the condition that $\langle A \rangle = 0$. Note that for any function D,

$$\begin{aligned} \left\langle \int_{0}^{\phi} D d\phi' \right\rangle &= \frac{1}{2\pi} \int_{0}^{2\pi} d\phi \int_{0}^{\phi} D d\phi' \\ &= 2\pi \langle D \rangle - \langle \phi D \rangle, \end{aligned} \quad (B10)$$

after switching the order of integration. It is then simple to show that $C = \langle \phi B \rangle - \pi \langle B \rangle$, so that the general, average-free solution of Eq. (B9) is

$$A = \int_{0}^{\phi} B d\phi' - (\phi + \pi) \langle B \rangle + \langle \phi B \rangle. \quad (B11)$$

Note that $\langle \partial A/\partial \phi \rangle = 0$, and that $A(0) = A(2n\pi) = \langle (\phi - \pi)B \rangle$, for any integer $n$; thus A is also periodic. Averages of expressions involving integrals satisfy the useful property

$$\left\langle A \int_{0}^{\phi} B d\phi' \right\rangle = - \left\langle B \int_{0}^{\phi} A d\phi' \right\rangle + 2\pi \langle A \rangle \langle B \rangle. \quad (B12)$$

We now iterate Eqs. (B7) in powers of $\epsilon$. We first expand

$$Q_{\alpha}(\tilde{X}_{\beta} + \epsilon Y_{\beta}, \phi) \equiv \sum_{m=0}^{\infty} \frac{\epsilon^m}{m!} Q_{\alpha, \beta \dots \gamma}^{(0)} Y_{\beta} \dots Y_{\gamma}, \quad (B13)$$

where we sum over repeated indices $\beta, \gamma$, etc., and where

$$Q_{\alpha}^{(0)} \equiv Q_{\alpha}(\tilde{X}_{\beta}, \phi), \quad (B14)$$

$$Q_{\alpha, \beta \dots \gamma}^{(0)} \equiv \frac{\partial^m Q_{\alpha}^{(0)}}{\partial \tilde{X}_{\beta} \dots \partial \tilde{X}_{\gamma}}. \quad (B15)$$

We also expand

$$Y_{\alpha} \equiv Y_{\alpha}^{(0)} + \epsilon Y_{\alpha}^{(1)} + \epsilon^2 Y_{\alpha}^{(2)} + \dots \quad (B16)$$

Then, to order $\epsilon^0$, we have

$$\frac{d\tilde{X}_\alpha}{d\theta} = \langle Q_\alpha^{(0)} \rangle, \quad (B17)$$

and

$$\frac{dY_{\alpha}^{(0)}}{d\phi} = \mathcal{A}\mathcal{F}(Q_{\alpha}^{(0)}), \quad (B18)$$

with the solution

$$Y_{\alpha}^{(0)} = \int_{0}^{\phi} Q_{\alpha}^{(0)} d\phi' - (\phi + \pi) \langle Q_{\alpha}^{(0)} \rangle + \langle \phi Q_{\alpha}^{(0)} \rangle. \quad (B19)$$

At first order in $\epsilon$, we find, from Eq. (B7b)

$$\frac{dY_{\alpha}^{(1)}}{d\phi} = \mathcal{A}\mathcal{F}(Q_{\alpha, \beta}^{(0)} Y_{\beta}^{(0)}) - Y_{\alpha, \beta}^{(0)} \frac{d\tilde{X}_{\beta}}{d\theta}. \quad (B20)$$

Substituting for $Y_\alpha^{(0)}$ from Eq. (B19) and for $d\tilde{X}_\beta/d\theta$ from Eq. (B17), and making use of (B11), we obtain
$$Y_{\alpha}^{(1)} = \int_{0}^{\phi} Q_{\alpha,\beta}^{(0)} d\phi' \int_{0}^{\phi'} Q_{\beta}^{(0)} d\phi'' + [\langle\phi Q_{\beta}^{(0)}\rangle - (\phi + \pi)\langle Q_{\beta}^{(0)}\rangle] \int_{0}^{\phi} Q_{\alpha,\beta}^{(0)} d\phi' + (\phi + \pi) \langle Q_{\beta}^{(0)} \int_{0}^{\phi} Q_{\alpha,\beta}^{(0)} d\phi' \rangle \\ + \langle\phi Q_{\alpha,\beta}^{(0)} \int_{0}^{\phi} Q_{\beta}^{(0)} d\phi' \rangle + \frac{1}{2}(\phi^2 - 4\pi^2/3)\langle Q_{\beta}^{(0)} \rangle \langle Q_{\alpha,\beta}^{(0)} \rangle - (\phi + \pi)\langle\phi Q_{\beta}^{(0)}\rangle \langle Q_{\alpha,\beta}^{(0)}\rangle \\ - \frac{1}{2}\langle Q_{\beta}^{(0)} \rangle \langle\phi(\phi + 2\pi)Q_{\alpha,\beta}^{(0)}\rangle + \langle\phi Q_{\beta}^{(0)} \rangle \langle\phi Q_{\alpha,\beta}^{(0)}\rangle. \qquad (B21)$$

It is straightforward to show explicitly that $\langle Y_{\alpha}^{(1)} \rangle = 0$.

To obtain $d\tilde{X}_{\alpha}/d\theta$ to order $\epsilon^2$, we expand Eq. (B7a) to second order in $\epsilon$ and substitute Eqs. (B19) and (B21) for $Y_{\alpha}^{(0)}$ and $Y_{\alpha}^{(1)}$, to obtain an expression purely in terms of derivatives and averages of $Q_{\alpha}^{(0)}$. Converting back to the unscaled $\phi = \theta/\epsilon$, we obtain

$$\frac{d\tilde{X}_{\alpha}}{d\phi} = \epsilon\langle Q_{\alpha}^{(0)} \rangle + \epsilon^2 \left[ \langle Q_{\alpha,\beta}^{(0)} \int_{0}^{\phi} Q_{\beta}^{(0)} d\phi' \rangle + \langle Q_{\alpha,\beta}^{(0)} \rangle \langle \phi Q_{\beta}^{(0)} \rangle \right. \\ \left. - \langle \phi Q_{\alpha,\beta}^{(0)} \rangle \langle Q_{\beta}^{(0)} \rangle - \pi \langle Q_{\alpha,\beta}^{(0)} \rangle \langle Q_{\beta}^{(0)} \rangle \right] \\ + \epsilon^3 \left[ \langle Q_{\alpha,\beta}^{(0)} Y_{\beta}^{(1)} \rangle + \frac{1}{2} \langle Q_{\alpha,\beta}^{(0)} Y_{\beta}^{(0)} Y_{\gamma}^{(0)} \rangle \right] + O(\epsilon^4), \qquad (B22)$$

where $Y_{\beta}^{(0)}$ and $Y_{\beta}^{(1)}$ are given by Eqs. (B19) and (B21). The first term in Eq. (B22) is the standard lowest-order

result in which “constant” values of the orbit elements are inserted into $Q_{\alpha}$ and the result is averaged over one period. The higher-order terms result from the effect of periodic terms in the orbit elements on the behavior of the $Q_{\alpha}$.

For the conservative parts of the equations of motion used in this paper, the first term is of 1PN order and higher, the second is of 2PN order and higher and the third is of 3PN order and higher; thus the terms shown are sufficient to calculate the conservative evolution of the elements through 3PN order. For the radiation-reaction evolution, the first term is of 2.5PN order and higher. Since we are only concerned with effects that are linear in the mass ratio $\eta$, the higher-order terms will contribute only cross-terms between radiation-reaction and conservative effects. Thus the second term will be of 3.5PN order and higher, while the final term will be of 4.5PN order and higher. Thus the terms shown will be sufficient to determine the radiation-reaction evolution through 4.5PN order.

[1] S. Babak, J. R. Gair, and R. H. Cole, in *Equations of Motion in Relativistic Gravity*, edited by D. Pützfeld, C. Lämmerzahl, and B. Schutz, Fundamental Theories of Physics Vol. 179 (Springer, Switzerland, 2015), pp. 783–812.

[2] T. Alexander, Annu. Rev. Astron. Astrophys. 55, 1 (2017).

[3] E. Poisson, A. Pound, and I. Vega, Living Rev. Relativ. 14, 7 (2011).

[4] M. van de Meent, Phys. Rev. D 94, 044034 (2016).

[5] K. Glampedakis, S. A. Hughes, and D. Kennefick, Phys. Rev. D 66, 064005 (2002).

[6] J. R. Gair and K. Glampedakis, Phys. Rev. D 73, 064037 (2006).

[7] T. Hinderer and É. É. Flanagan, Phys. Rev. D 78, 064028 (2008).

[8] S. A. Hughes, arXiv:1601.02042.

[9] C. Hopman and T. Alexander, Astrophys. J. 645, 1152 (2006).

[10] D. Merritt, T. Alexander, S. Mikkola, and C. M. Will, Phys. Rev. D 84, 044024 (2011).

[11] A. S. Hamers, S. F. Portegies Zwart, and D. Merritt, Mon. Not. R. Astron. Soc. 443, 355 (2014).

[12] B. Bar-Or and T. Alexander, Classical Quantum Gravity 31, 244003 (2014).

[13] P. Brem, P. Amaro-Seoane, and C. F. Sopuerta, Mon. Not. R. Astron. Soc. 437, 1259 (2014).

[14] B. Bar-Or and T. Alexander, Astrophys. J. 820, 129 (2016).

[15] C. M. Will, Phys. Rev. D 89, 044043 (2014).

[16] C. M. Will, Proc. Natl. Acad. Sci. U.S.A. 108, 5938 (2011).

[17] C. M. Will, Classical Quantum Gravity 29, 217001 (2012).

[18] E. Poisson and C. M. Will, *Gravity: Newtonian, Post-Newtonian, Relativistic* (Cambridge University Press, Cambridge, England, 2014).

[19] M. E. Pati and C. M. Will, Phys. Rev. D 65, 104008 (2002).

[20] C. M. Will, Phys. Rev. D 71, 084027 (2005).

[21] A. Gopakumar, B. R. Iyer, and S. Iyer, Phys. Rev. D 55, 6030 (1997).

[22] D. Brouwer and G. M. Clemence, *Methods of Celestial Mechanics* (Academic Press, New York, 1961).

[23] C. M. Bender and S. A. Orszag, *Advanced Mathematical Methods for Scientists and Engineers* (McGraw-Hill, New York, 1978).

[24] C. W. Lincoln and C. M. Will, Phys. Rev. D 42, 1123 (1990).
[25] T. Damour and N. Deruelle, Ann. Inst. Henri Poincaré **43**, 107 (1985).

[26] T. Damour and N. Deruelle, Ann. Inst. Henri Poincaré **44**, 263 (1986).

[27] T. Mora and C. M. Will, Phys. Rev. D **69**, 104021 (2004).

[28] L. E. Kidder, Phys. Rev. D **52**, 821 (1995).

[29] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Phys. Rev. Lett. **74**, 3515 (1995).

[30] L. Á. Gergely, Z. I. Perjés, and M. Vasúth, Phys. Rev. D **57**, 876 (1998).

[31] N. Sago and R. Fujita, Prog. Theor. Exp. Phys. **2015**, 073E03 (2015).

[32] R. Rieth and G. Schäfer, Classical Quantum Gravity **14**, 2357 (1997).

[33] K. G. Arun, L. Blanchet, B. R. Iyer, and M. S. S. Qusailah, Phys. Rev. D **77**, 064035 (2008).

[34] K. G. Arun, L. Blanchet, B. R. Iyer, and S. Sinha, Phys. Rev. D **80**, 124018 (2009).

[35] E. Forseth, C. R. Evans, and S. Hopper, Phys. Rev. D **93**, 064058 (2016).

[36] C. R. Galley and I. Z. Rothstein, arXiv:1609.08268.

[37] P. C. Peters, Phys. Rev. **136**, B1224 (1964).

[38] R. M. O'Leary, B. Kocsis, and A. Loeb, Mon. Not. R. Astron. Soc. **395**, 2127 (2009).

[39] P. Amaro-Seoane, C. F. Sopuerta, and M. D. Freitag, Mon. Not. R. Astron. Soc. **429**, 3155 (2013).

[40] C. F. Sopuerta and N. Yunes, Phys. Rev. D **84**, 124060 (2011).

[41] A. Buonanno and T. Damour, Phys. Rev. D **59**, 084006 (1999).

[42] A. Buonanno and T. Damour, Phys. Rev. D **62**, 064015 (2000).

[43] D. Bini and T. Damour, Phys. Rev. D **86**, 124012 (2012).

[44] R. H. Boyer and R. W. Lindquist, J. Math. Phys. (N.Y.) **8**, 265 (1967).

[45] G. B. Cook and M. A. Scheel, Phys. Rev. D **56**, 4775 (1997).

[46] G. B. Cook, Living Rev. Relativ. **3**, 5 (2000).

[47] S. Hergt and G. Schäfer, Phys. Rev. D **77**, 104001 (2008).

[48] B. Carter, Phys. Rev. **174**, 1559 (1968).

[49] B. R. Iyer and C. M. Will, Phys. Rev. Lett. **70**, 113 (1993).

[50] B. R. Iyer and C. M. Will, Phys. Rev. D **52**, 6882 (1995).