Post-Newtonian expansions of extreme mass ratio inspirals of spinning bodies into Schwarzschild black holes

Abstract

Space-based gravitational-wave detectors such as LISA are expected to detect inspirals of stellar-mass compact objects into massive black holes. Modeling such inspirals requires fully relativistic computations to achieve sufficient accuracy at leading order. However, subleading corrections such as the effects of the spin of the inspiraling compact object may potentially be treated in weak-field expansions such as the post-Newtonian (PN) approach.

In this work, we calculate the PN expansion of eccentric orbits of spinning bodies around Schwarzschild black holes. Then we use the Teukolsky equation to compute the energy and angular momentum fluxes from these orbits up to the 5PN order. Some of these PN orders are exact in eccentricity, while others are expanded up to the tenth power in eccentricity. Then we use the fluxes to construct a hybrid inspiral model, where the leading part of the fluxes is calculated numerically in the fully relativistic regime, while the part linear in the small spin is analytically approximated using the PN series. We calculate LISA-relevant inspirals and respective waveforms with this model and a fully relativistic model. Through the calculation of mismatch between the waveforms from both models we conclude that the PN approximation of the linear-in-spin part of the fluxes is sufficient for lower eccentricities.

I Introduction and summary

I.1 Extreme mass ratio inspirals

Forthcoming space-based gravitational wave (GW) detectors such as LISA, TianQin, or Taiji [1, 2, 3, 4] will be able to detect signals from various sources, including extreme mass ratio inspirals (EMRIs) [5]. These systems consist of a stellar mass black hole or neutron star in orbit around a massive black hole with the mass ratio $\epsilon=\mu/M$ of the small (secondary) mass $\mu$ and large (primary) mass $M$ between $10^{-7}$ and $10^{-4}$. Because of gravitational radiation reaction, the orbit of the small body decays and it inspirals into the primary body while radiating GWs to infinity. Because the secondary body completes many orbits in the strong gravitational field of the primary body, the detection of GWs from such systems will give a unique insight into the strong-field regime around massive black holes, which will also allow us to test general relativity to high precision [6, 7]. Furthermore, the study of EMRI populations will provide new insights in cosmology and astrophysics [5, 8].

To achieve the aforementioned goals, the parameters of the detected systems must be estimated with high precision. Because signals from EMRIs and other astrophysical sources will overlap, detection and parameter estimation will be done through matched filtering, which is based on the comparison of the received signal with many waveform templates. For this purpose, the waveforms must be generated for a wide range of parameters with phase accurate to fractions of radian [9]. There are several methods for modeling binary systems, and the choice of the most suitable method depends on the parameters of the system, such as the mass ratio and compactness.

I.2 Black hole perturbation theory

In particular, for the modeling of EMRIs, black hole perturbation theory (BHPT) is often employed, where the spacetime is expanded in the mass ratio around a background spacetime of the primary [10]. Then, the system can be effectively described as a point particle moving in the background spacetime while inducing a perturbation of this spacetime. This perturbation acts on the particle with the so-called self-force, which can be expanded in the powers of the mass ratio. Because the mass ratio and, therefore, the perturbation is small, the inspiral timescale is much slower than the orbital timescale. Thus, to efficiently solve the problem, two-timescale expansion is often used, where the system is described using a set of orbital parameters $\mathcal{J}_{i}$, which evolve slowly ($\dot{\mathcal{J}}_{i}=\order{\epsilon}$) and a set of phases $\psi^{i}$ that evolve quickly ($\dot{\psi}^{i}=\order{1}$) [11]. The phases are directly related to the phase of the GW. When we consider an inspiral that sweeps through some finite range of frequencies such as a GW detector band, we can use the separation of scales to expand the phase elapsed during this process as

$$\Phi_{\text{GW}}=\epsilon^{-1}\Phi_{0}(t)+\Phi_{1}(t)+\order{\epsilon}\,.$$ (1)

The first term, which is called the adiabatic term, is of the order of $\epsilon^{-1}$ radians, while the second, postadiabatic term, is in the order of radians and cannot be neglected to achieve subradian accuracy. The adiabatic term consists of the contribution from the time-averaged dissipative (time-antisymmetric) part of the first-order in the mass ratio self-force, while the postadiabatic term consists of a number of contributions from different physical effects [10]. Namely, the postadiabatic term requires the inclusion of the oscillating part of the dissipative and conservative (time-symmetric) first-order self-force, time-averaged dissipative part of the second-order self-force, the force caused by the spin-curvature coupling of the secondary, and corrections to the dissipative self-force caused by the secondary spin.

To find all the contributions up to the postadiabatic term, one in principle has to find the metric perturbation up to the second order in the mass ratio, regularize it near the particle, and calculate the self-force from the regular part. However, thanks to flux-balance laws, the averaged dissipative part of the self-force can often be obtained from the asymptotic GW fluxes to infinity and through the horizon of the primary black hole [12, 13, 14, 15]. The first-order flux-balance laws for non-spinning secondaries were obtained by Mino [12] and Sago et al. [13]. For spinning secondaries, the flux-balance law was proven only for the evolution of energy and azimuthal angular momentum [14]. Second-order flux-balance laws for the energy and azimuthal angular momentum of non-spinning secondaries on quasi-circular orbits in Schwarzschild space-time were derived by Miller and Pound [16]; these derivations are expected to straightforwardly generalize to generic orbits. The currently open question is a concrete formulation of some sort of flux-balance law for the so-called Carter constant evolve at second order in the mass ratio and under secondary spin corrections to the motion (see Ref. [17] for some recent effort in this direction). A less obvious quantity that did not have a flux-balance law to date was the aligned component of the secondary spin $s_{\parallel}$; this question is resolved by us in Section III.

The first-order adiabatic fluxes must be calculated with high accuracy since the error will be enhanced by $\epsilon^{-1}$ compared to the postadiabatic terms. As a general rule, the error of the adiabatic term must be $\order{\epsilon}$ smaller than the error of the postadiabatic term. This opens up the possibility of using various approximations for the calculation of the postadiabatic effects.

I.3 Post-Newtonian expansions

As mentioned above, there are other techniques to model binary systems with different mass ratios and separations. One of such techniques is the post-Newtonian (PN) theory [18], which is valid for systems with large separations and low relative velocities. This method relies on expanding the Einstein equations in the inverse square of the speed of light in vacuum $1/c^{2}$, thus assuming that quantities such as the dimensionless speed squared $v^{2}/c^{2}$ or the dimensionless gravitational potential $GM/(rc^{2})$ are small. Currently, the state-of-the-art results for comparable-mass spinning binaries on eccentric orbits are expansions of the energy and angular momentum fluxes to 3PN order [19] beyond the Newtonian quadrupole formula [20] and for nonspinning objects on circular orbits to 4.5PN [21].

The regime in which both PN theory and BHPT are valid offers the possibility of cross-validating the results of both theories. In particular, when the results of BHPT are analytically expanded in a PN parameter (see a review of older results in [22, 23]), direct comparisons can be made with pure PN computations truncated at the first order in the mass ratio. Furthermore, the BHPT computations can typically be expanded to higher PN orders than the existing PN computations at finite mass ratio. Finally, careful considerations of the symmetries of the mass ratio expansion of the PN series reveal that the BHPT results can often have a “strategic” importance for obtaining unknown pieces of the equations of motion of binaries at any mass ratio [24, 25].

Such results can also be utilized to calibrate effective-one-body models, which is an approach to binary modeling that takes input from numerical relativity, PN theory, and BHPT [26, 27, 28]. In particular, the dynamics of spinning test particles in black hole space-times proved to also be useful in the development of effective-one-body models (see e.g. [29, 30, 31]).

The PN expansion of BHPT results was first used by Poisson [32], where the energy fluxes to infinity from circular orbits in the Schwarzschild spacetime were expanded to 1.5PN orders beyond the Newtonian order. These results were then extended to higher PN orders, to fluxes through the horizon and to the Kerr spacetime [33, 34, 35, 36, 37, 38]. The latest results are infinity fluxes and horizon fluxes from circular orbits in the Schwarzschild spacetime to 22PN [39] and in the Kerr spacetime to 11PN order [40].

The effects of the spin of the secondary body were first incorporated into the fluxes from circular orbits around a Kerr black hole by Tanaka et al. [41] to 2.5PN order, and later by Nagar et al. [30] and Akcay et al. [14] for circular orbits in the Schwarzschild spacetime to the 5.5PN order.

The formalism was also extended to generic orbits of non-spinning bodies in Kerr spacetime, where one needs the evolution of three constants of motion, namely the energy, angular momentum, and the Carter constant [13, 42, 43]. The latest results, i.e. 5PN fluxes with expansions in eccentricity to $e^{10}$ were used by Isoyama et al. [44] to generate generic adiabatic inspirals and waveforms.

Another direction in which this technique was utilized was to calculate PN expansions of energy and angular-momentum fluxes from highly eccentric orbits in the Schwarzschild spacetime. In Ref. [45] the authors identified singular factors in the form $(1-e^{2})^{-k}$ yielding convergent series in eccentricity after the factorization of such terms. In addition, they used highly accurate numerical calculations to find the coefficients of the series in eccentricity and the PN series to 7PN order. Later, an analytical form for the leading and subleading logarithmic terms was found by Munna and Evans [46, 47]. Finally, the energy and angular momentum fluxes to infinity and through the horizon were found up to the 19PN order [48, 49, 50].

These expansions not only provided validation of the results of the PN theory, but were used by Burke et al. [51] in the calculation of adiabatic inspirals, where the authors also tested the possibility of using the waveforms derived from such inspirals for accurate parameter estimation. It was found that the 9PN adiabatic fluxes from eccentric orbits in Schwarzschild spacetime introduce bias on the system parameters and, therefore, cannot be used instead of the fully relativistic fluxes. However, in the same work, a hybrid model with fully relativistic adiabatic (first-order in the mass ratio) fluxes and 3PN expansion of postadiabatic fluxes (second-order in the mass ratio) was used, which was proven to be sufficient for accurate parameter estimation in some cases.

The secondary spin corrections to the fluxes are of the order of the mass ratio and, consequently, contribute at postadiabatic order, the same as the PN-expanded pieces used in Burke et al. [51]. Therefore, it may be possible to approximate them using PN expansion and avoid computationally expensive numerical calculations of the fully relativistic fluxes such as was done in Refs. [52, 53, 14, 54, 55, 56].

I.4 Summary of results

• In this work, we PN-expanded the analytical expressions for eccentric, precessing trajectories of spinning bodies in Schwarzschild spacetime that were recently found by Witzany and Piovano [57]. The expanded trajectories and other relevant quantities can be found in a Supplemental notebook [58].

• Then, we employed the Teukolsky equation to find the energy and angular momentum fluxes from these orbits as a closed-form series in the PN parameter and eccentricity. We linearized the fluxes in the secondary spin and found the linear-in-spin correction up to 5PN and at least 10th power in eccentricity. We were able to fully factorize and resum the fluxes as a finite series in eccentricity up to 2.5 PN with partial resummation results also at higher orders. We demonstrated that the resulting eccentricity series converges even up to $e\to 1$ in Figures 1 and 2. The resulting spin corrections to fluxes are in equations (55),(56) and Appendix A as well as the Supplemental notebook.

• We tested the convergence of the PN series by analytically integrating the phase evolution of quasi-circular inspirals with the result in equation (LABEL:eq:Deltaphi). Using this general result, we computed the phase contributions of LISA-band inspirals of $100M_{\odot}$ spinning black holes into massive black holes of mass $10^{6}M_{\odot}$ in Table 1. This demonstrated that the 5PN expansion is not sufficiently accurate for the nonspinning part of the flux, but it is sufficient for the spin correction in the case of quasi-circular inspirals.

• Hence, we then used these flux corrections in a hybrid inspiral model, where the nonspinning part was calculated numerically in a fully relativistic regime and the linear-in-spin part is expressed analytically as a PN series. To be able to do so, we also derived that the aligned component of the secondary spin $s_{\parallel}$ will stay conserved during generic EMRIs.

• Using the hybrid model, we computed fiducial LISA-band eccentric inspirals using the same binary masses as for the quasi-circular case. Additionally, we used the FEW package [59, 60, 61, 62] to generate relativistic waveforms corresponding to the inspirals. As a test of the formalism, we computed LISA mismatches of the hybrid-model waveforms with waveforms corresponding to fully relativistic inspirals. The mismatches presented in Figure 8 imply that the hybrid model should be adequate for the detection of the vast majority of LISA EMRIs. Additionally, it should not introduce significant biases for parameter estimates of less eccentric events.

I.5 Organization of paper

This paper is organized as follows. Section II reviews the motion of spinning bodies in Schwarzschild spacetime and introduces PN and eccentricity expansions of these trajectories. This is followed by Section III where the self-torque acting on the spin vector is presented, which is then used to derive the adiabatic evolution of the parallel component of the secondary spin. Next, Section IV examines the GW fluxes from orbits described in the previous Section. First, the Teukolsky formalism is presented which is then used to calculate the PN expansions of the fluxes. Next, Section V presents the hybrid model for the adiabatic inspirals that includes the PN expansion of the fluxes and the calculation of inspirals using this model. Finally, Section VI provides a discussion of the importance of the results and outlooks.

I.6 Notation

Geometrized units, where the gravitational constant and the speed of light in vacuum are set to unity ($G=c=1$), are used throughout this paper. Spacetime indices are denoted with Greek letters, while tetrad indices are denoted with Latin letters. The signature of the metric is $(-,+,+,+)$, while the Riemann tensor is defined as $a_{\nu;\kappa\lambda}-a_{\nu;\lambda\kappa}\equiv R^{\mu}{}_{\nu\kappa\lambda}a_{\mu}$, where the semicolon denotes the covariant derivative and $a_{\mu}$ is a covector. The sign of the Levi-Civita pseudotensor is defined as $\sqrt{-g}\epsilon_{tr\theta\phi}=-\epsilon^{tr\theta\phi}/\sqrt{-g}=1$.

II PN expansion of eccentric equatorial trajectories of spinning bodies

II.1 Spinning-particle trajectory

Let us briefly summarize the properties of the closed analytical solution of the bound motion of spinning particles near Schwarzschild black holes as presented in Ref. [57]. We consider the motion in Schwarzschild spacetime given as

$$\differential s^{2}=-f(r)\differential t^{2}+\frac{1}{f(r)}\differential r^{2}+r^{2}(\differential\theta^{2}+\sin^{2}\theta\differential\phi^{2})\,,$$ (2)

where $f(r)=1-2M/r$. The motion of the spinning particle is described by Mathisson-Papapetrou-Dixon equations under the Tulczyjew-Dixon or covariant spin supplementary condition $s^{\mu\nu}P_{\nu}=0$, where $P_{\mu}$ is the particle momentum and $s^{\mu\nu}$ is the spin tensor per unit particle mass. The solution is valid up to $\mathcal{O}(s)$ corrections to the orbital motion, and to leading order in the spin sector. In this truncation one has $P_{\mu}=\mu u_{\mu}+\mathcal{O}(s^{2})$, where $\mu$ is the particle mass and $u_{\mu}$ is the four-velocity. One then equivalently parametrizes the spin by the spin vector and the spin tensor

$$s^{\mu}=-\frac{1}{2}\epsilon^{\mu\nu\kappa\lambda}s_{\nu\kappa}u_{\lambda}\,,$$ (3)

$$s^{\mu\nu}=\epsilon^{\mu\nu\kappa\lambda}u_{\kappa}s_{\lambda}\,,$$ (4)

where $\epsilon^{\mu\nu\kappa\lambda}$ is the Levi-Civita pseudotensor.

Note that the definition of $s^{\lambda}$ depends on the orientation of the basis. This is further complicated by the fact that raising or lowering indices formally changes the orientation. Here we fix the convention by

$$\epsilon_{tr\theta\phi}=-\frac{1}{(-g)}\epsilon^{tr\theta\phi}=\frac{1}{\sqrt{-g}}\,,$$ (5)

which yields a right-handed basis for upper-index quantities under the assumption of a conventional transformation from $r,\theta,\phi$ to Cartesian coordinates. However, this also means that our formulas have a relative minus sign in front of any spin correction as compared to Ref. [57].

The 3 rotational symmetries of the Schwarzschild space-time around the $x,y,z$ axes generate a conserved total angular momentum vector of the generically inclined spinning particle. Upon rotation of the coordinate equator $\theta=\pi/2$ into the plane perpendicular to this vector, the generic motion becomes near-equatorial in the resulting frame, $\theta=\pi/2+\delta\theta+\mathcal{O}(s^{2})$. As a result, the Mathisson-Papapetrou-Dixon equations fully separate and the solutions are parametrized by the three constants of motion

$$E=-u_{\mu}\xi_{(t)}^{\mu}+\frac{1}{2}\xi^{(t)}_{\mu;\nu}s^{\mu\nu}\,,$$ (6)

$$J_{z}=u_{\mu}\xi_{(\phi)}^{\mu}-\frac{1}{2}\xi^{(\phi)}_{\mu;\nu}s^{\mu\nu}\,,$$ (7)

$$s_{\parallel}=\frac{s^{\mu}l_{\mu}}{\sqrt{l_{\alpha}l^{\alpha}}}=\frac{Y_{\mu\nu}s^{\mu}u^{\nu}}{2\sqrt{K_{\alpha\beta}u^{\alpha}u^{\beta}}}\,,$$ (8)

$$l^{\mu}\frac{\partial}{\partial x^{\mu}}\equiv\frac{r\dot{\theta}}{\sin\theta}\frac{\partial}{\partial\phi}-r\sin\theta\dot{\phi}\frac{\partial}{\partial\theta}\,,$$ (9)

where $E$ has the meaning of total orbital and spin-orbital energy per unit mass, $J_{z}$ the orbital and spin-orbital angular momentum, and $s_{\parallel}$ is the component of spin aligned with the orbital angular momentum. Furthermore, $Y_{\mu\nu}=-Y_{\nu\mu},Y_{\mu\nu;\kappa}=-Y_{\mu\kappa;\nu}$ is the Killing-Yano tensor of the Schwarzschild spacetime, which means that $s_{\parallel}$ is related to the more general Rüdiger constants in Kerr spacetime and the separation constant for spinning particles in Kerr found by separation of the Hamilton-Jacobi equation [63, 64]. It should also be noted that in the aligned frame the magnitude of the total angular momentum is by construction equal to the single component $J_{z}$.

The solution is parametrized by Carter-Mino time $\differential\lambda=\differential\tau/r^{2}$, where $\tau$ is proper time. The radial solution is then expressed in the form

$$r(\lambda)=\frac{r_{3}(r_{1}-r_{2})\mathrm{sn}^{2}\left(\frac{K(k)}{\pi}q^{r},k\right)-r_{2}(r_{1}-r_{3})}{(r_{1}-r_{2})\mathrm{sn}^{2}\left(\frac{K(k)}{\pi}q^{r},k\right)-(r_{1}-r_{3})}\,,$$ (10)

$$k^{2}=\frac{(r_{1}-r_{2})r_{3}}{(r_{1}-r_{3})r_{2}}\,,$$ (11)

$$q^{r}\equiv\Upsilon^{r}\lambda+q^{r}_{0}\,,$$ (12)

where $K(k)$ is the complete elliptic integral of the first kind, $\mathrm{sn}()$ is the Jacobi sn function, and $q^{r}_{0}$ is an integration constant determined by initial conditions. $\Upsilon^{r}$ is given in Ref. [57] and represents the fundamental frequency of motion with respect to Mino time. The radii $r_{1}>r_{2}>r_{3}$ are then the non-zero roots of the polynomial $R(r)$ appearing in the radial equation of motion

$$\frac{\differential r}{\differential\lambda}=\pm\sqrt{R(r)}\,,$$ (13)

$$\begin{split}&R(r)\equiv r\Big{[}r^{3}(E^{2}-1)-rJ_{z}(J_{z}-2s_{\parallel}E)\\ &\phantom{R(r)=}+2M(r^{2}+J_{z}(J_{z}-3s_{\parallel}E))\Big{]}\,\\ &\phantom{R(r)}=(1-E^{2})r(r_{1}-r)(r-r_{2})(r-r_{3})\,.\end{split}$$ (14)

The roots $r_{1},r_{2}$ are the physical turning points of the bound motion and are conventionally parameterized by the orbital parameters dimensionless semilatus rectum $p$ and eccentricity $e$ defined through the relation

$$r_{1}=\frac{Mp}{1-e}\;,\qquad r_{2}=\frac{Mp}{1+e}\,.$$ (15)

One can then express $E(p,e,s_{\parallel}),J_{z}(p,e,s_{\parallel}),r_{3}(p,e,s_{\parallel})$ in closed form by examining equation (14).

The $t,\phi$ degrees of freedom are then given as

$$t(\lambda)=q^{t}+\Delta t(q^{r})\,,\;\phi(\lambda)=q^{\phi},$$ (16)

$$q^{t}\equiv\Upsilon^{t}\lambda+q^{t}_{0}\,,\;q^{\phi}\equiv\Upsilon^{\phi}\lambda+q^{\phi}_{0},\,$$ (17)

$$\Delta t(q^{r})=\tilde{T}_{r}\left(\mathrm{am}\left(\frac{q^{r}}{\pi}K(k),k\right)\right)-\frac{\tilde{T}_{r}\left(\pi\right)}{2\pi}q^{r}\,,$$ (18)

where $\mathrm{am}()$ is the Jacobi amplitude, $\tilde{T}_{r}$ is given in eq. (49) of Ref. [57], $\Upsilon^{t},\Upsilon^{\phi}$ are the $t,\phi$ Mino frequencies, and $q^{t}_{0},q^{\phi}_{0}$ are integration constants. It is interesting to note that unlike in Kerr, in Schwarzschild space-time the Carter-Mino time $\lambda$ is simply proportional to $\phi$, which means that the $\phi$ counterpart of $\Delta t(q^{r})$ vanishes in equation (16).

Finally, the spin degree of freedom and $\delta\theta$ depend on a precession angle $\psi$ with the evolution

$$\psi(\lambda)=q^{\psi}+\tilde{\Psi}_{r}\left(\mathrm{am}\left(\frac{q^{r}}{\pi}K(k),k\right)\right)-\frac{\tilde{\Psi}_{r}\left(\pi\right)}{2\pi}q^{r},$$ (19)

$$q^{\psi}\equiv\Upsilon^{\psi}\lambda+q^{\psi}_{0},\,$$ (20)

where $\Upsilon^{\psi},q^{\psi}_{0}$ again have analogous meanings as above and $\tilde{\Psi}_{r}$ is a known function.

The deviation from the equatorial plane is then given as

$$\delta\theta=-\frac{\sqrt{(s^{2}-s_{\parallel}^{2})(J_{z}^{2}+r^{2})}\sin\psi}{J_{z}r}\,,$$ (21)

and the spin vector can be expressed as

$$s^{t}=-\sqrt{\frac{s^{2}-s_{\parallel}^{2}}{f(J_{z}^{2}+r^{2})}}\left(\frac{EJ_{z}\cos\psi}{\sqrt{f}}+\dot{r}r\sin\psi\right)\,,$$ (22a)

$$s^{r}=-\sqrt{s^{2}-s_{\parallel}^{2}}\left(\frac{J_{z}\dot{r}\cos\psi}{r}+\frac{Er\sin\psi}{\sqrt{J_{z}^{2}+r^{2}}}\right)\,,$$ (22b)

$$s^{\phi}=-\frac{\sqrt{(s^{2}-s_{\parallel}^{2})(J_{z}^{2}+r^{2})}\cos\psi}{r^{2}}\,,$$ (22c)

$$s^{\theta}=\frac{s_{\parallel}}{r}\,,$$ (22d)

where $\dot{r}$ is the proper-time derivative of $r$ expressed as $r^{2}\sqrt{R(r)}$. Note that even though the spin is parametrized by $s_{\parallel}$, the orientation of the spin vector is generic, and we are thus dealing with absolutely generic bound orbits of spinning test particles in Schwarzschild spacetime in this paper.

II.2 PN expansion of the trajectory

The constants of motion, orbital frequencies, and trajectory $(t,r,\phi)$ as a function of the phase $q^{r}$ can be expanded in a formal PN parameter. In this work we use the parameter

$$v=\sqrt{\frac{1}{p}}\;.$$ (23)

Other choices include e.g. the gauge independent parameter $x=(M\Omega_{\phi})^{2/3}$, however, the parameter $v$ is convenient when the orbit is parametrized with $p$ and $e$ and one can reexpress the final result in different PN parameters. Every order in $v$ corresponds to one half of the PN order, that is, the expansion to 7 orders in $v$ NLO (next to the leading order) corresponds to $3.5$PN orders NLO.

Since the expansions of the geodesic quantities were calculated before and are available in the literature [43], here we present only the PN expansion of the linear-in-spin correction of any given quantity. We write the expansion as $E=E_{\rm(g)}+s_{\parallel}\delta E/M$ and $J_{z}=J_{z\rm(g)}+s_{\parallel}\delta J_{z}/M$, where $E_{\rm(g)},J_{z\rm(g)}$ are the geodesic expressions at fixed orbital parameters $e,v$. Then it is straightforward to expand the linear-in-spin part of energy $\delta E$ and angular momentum $\delta J_{z}$ from Eqs. (32)-(33) in Ref. [57]. The results read

$$\delta E=-\frac{(1-e^{2})^{2}v^{5}}{2}\sum_{k=0}^{\infty}\binom{-3/2}{k}(-(3+e^{2})v^{2})^{k}\,,$$ (24)

$$\begin{split}&\frac{\delta J_{z}}{M}=1-2v^{2}-\frac{27+34e^{2}+3e^{4}}{8}v^{4}\\ &\phantom{\delta J_{z}=}-\frac{81+131e^{2}+39e^{4}+5e^{6}}{8}v^{6}+\order{v^{8}}.\end{split}$$ (25)

$\delta E$ was expanded to the order $v^{14}$ which corresponds to 9 orders NLO, while $\delta J_{z}$ was expanded to $v^{11}$ (here we show only the expansion to $v^{6}$ for simplicity; the full expressions can be found in the Supplemental Material [58]).

Since the parameter of the elliptic integrals $K$, $E$, and $\Pi$ in Eqs. (S22)-(S24) in Ref. [57] is $k^{2}\sim\order{ev^{2}}$, we can expand the expression in $k^{2}$. We then write the linear-in-spin parts of the orbital frequencies in Carter-Mino time as $\Upsilon=\Upsilon_{\rm(g)}+s_{\parallel}\delta\Upsilon/M$ and obtain

$$\begin{split}\frac{\delta\Upsilon^{t}}{M}&=\frac{e^{2}}{(1-e^{2})^{3/2}}v^{-1}+\quantity(\frac{9+8e^{2}}{2(1-e^{2})^{3/2}}-6)v\\ &\phantom{=}-\quantity(48+\frac{9e^{2}}{2}-\frac{330+255e^{2}-88e^{4}+9e^{6}}{(1-e^{2})^{3/2}})v^{3}\\ \phantom{\delta\Upsilon^{t}}&\phantom{=}+\order{v^{5}}\,,\end{split}$$ (26)

$$\begin{split}\delta\Upsilon^{r}&=\frac{1}{2}\left(3-e^{2}\right)v^{2}+\frac{1}{4}\left(-3e^{4}+4e^{2}+33\right)v^{4}\\ &\phantom{=}+\frac{1}{16}\left(-15e^{6}-23e^{4}+357e^{2}+693\right)v^{6}\\ &\phantom{=}+\order{v^{8}}\,,\end{split}$$ (27)

$$\begin{split}\delta\Upsilon^{\phi}&=(e^{2}+3)\bigg{(}-\frac{1}{2}v^{2}-\frac{1}{4}\left(3e^{2}+5\right)v^{4}\\ &\phantom{=}-\frac{1}{16}\left(15e^{4}+50e^{2}+63\right)v^{6}+\order{v^{8}}\bigg{)}\,.\end{split}$$ (28)

From the Mino time frequencies, we can calculate the coordinate time frequencies

$$\Omega_{i}=\frac{\Upsilon^{i}}{\Upsilon^{t}}\,,$$ (29)

and their linear-in-spin parts

$$\delta\Omega_{i}=\frac{\delta\Upsilon^{i}\Upsilon^{t}_{\rm(g)}-\Upsilon^{i}_{\rm(g)}\delta\Upsilon^{t}}{(\Upsilon^{t}_{\rm(g)})^{2}}\;.$$ (30)

where $i=r,\phi$. A nice special formula is that $\delta\Omega_{\phi}=-3v^{6}/2$ to all orders in $v$ for $e=0$.

We expanded the expressions to 9 orders in $v$ beyond the leading order while keeping the eccentricity dependence exact. The results were then expanded to $e^{10}$ for later calculations. Similarly to $\delta J_{z}$, here we show only the expansion to 4 orders NLO. The full expressions can be found in the Supplemental Material [58].

By expanding the Jacobi elliptic function ${\rm sn}(u,k)$ in $k^{2}$ in Eq. (10), we found the radial coordinate parametrized with $q^{r}$ as a series in $v$ and $e$ and extracted the linear-in-spin part $\delta r(q^{r})$. Again, due to the length of the expression, it is not included here but can be found in the Supplemental Material [58].

Next, we focused on $\Delta t(q^{r})$ which is the oscillating part of $t=(\Upsilon^{t}/\Upsilon^{r})q^{r}+\Delta t(q^{r})$. Note that the oscillating part of $\phi=(\Upsilon^{\phi}/\Upsilon^{r})q^{r}$ is zero. Since the expression for $t(q^{r})$ is too long, its PN expansion is computationally expensive. Instead, we expanded the equation

$$\derivative{t}{q^{r}}=\derivative{t}{\lambda}\derivative{\lambda}{q^{r}}=\derivative{t}{\lambda}\Upsilon_{r}^{-1}$$ (31)

in $v$ and $e$, where

$$\derivative{t}{\lambda}=\frac{r^{2}E}{f}-s_{\parallel}\frac{J_{z}f^{\prime}r^{2}}{2f}\;.$$ (32)

(C.f. Eq. (29) in Ref. [57].) In this way, we obtain a Fourier series of $\cos(nq^{r})$ that is trivial to integrate to obtain $t(q^{r})$. Then we extracted the linear-in-spin part $\delta\Delta t(q^{r})$, which is available in the Supplemental Material.

III Adiabatic evolution of the constants of motion

The metric perturbations sourced by the spinning secondary will lead to a self-torque and a self-force, which will drive its motion away from the motion of the spinning test body in the Schwarzschild metric. In this section, we derive the equations governing the leading-order secular evolution of the spinning-secondary orbit under this perturbation.

The evolution of the secondary under self-force and self-torque can be cast in the form of MPD equations in the effective regularized metric $\hat{g}_{\mu\nu}=g_{\mu\nu}+h_{\mu\nu}$ [65, 66, 14, 55] (we drop the conventional “R” index on the regularized metric perturbation $h_{\mu\nu}$ here for notational simplicity). As a result, under the assumption that the spin tensor $s^{\mu\nu}$ is unchanged in the perturbed metric, using Eq. (3) we obtain different definitions of the spin vector $s^{\mu}$ with respect to the Schwarzschild metric, and $\hat{s}^{\mu}$ with respect to the effective metric. The two definitions are related as follows

$$s^{\mu}=\hat{s}^{\mu}-\epsilon\quantity(\frac{1}{2}h_{\alpha}{}^{\alpha}\hat{s}^{\mu}+\frac{1}{2}h_{\alpha\beta}u^{\alpha}u^{\beta}\hat{s}^{\mu}-h^{\mu}{}_{\nu}\hat{s}^{\nu})\,.$$ (33)

Because the spin tensor $\hat{s}^{\mu}$ is parallel-transported in the effective metric, the spin vector on the Schwarzschild metric $s^{\mu}$ then experiences the self-torque

$$\frac{{\rm D}s^{\mu}}{\differential\tau}=-\frac{1}{2}h_{\alpha\beta;\rho}N^{\mu\alpha\beta\rho}\,,$$ (34)

$$\begin{split}&N^{\mu\alpha\beta\rho}=2g^{\mu\alpha}u^{[\beta}s^{\rho]}-g^{\mu\rho}u^{\alpha}s^{\beta}\\ &\phantom{N^{\mu\alpha\beta\rho}=}+g^{\alpha\beta}u^{\rho}s^{\mu}+u^{\alpha}u^{\beta}u^{\rho}s^{\mu}\,.\end{split}$$ (35)

The Tulczyjew-Dixon SSC in the effective metric $u^{\mu}\hat{s}_{\mu}=0$ is conserved up to $\mathcal{O}(s^{2})$ due to the general properties of the MPD equations in any metric. From equation (33) it can be seen that $\hat{s}^{\mu}$ is always $\mathcal{O}(\epsilon s)$ close to $s^{\mu}$ without secularly growing terms. As a result, $u^{\mu}s_{\mu}=\mathcal{O}(\epsilon s,s^{2})$ will also hold during evolution. Similarly, the spin magnitude with respect to the effective metric $\hat{s}^{\mu}\hat{s}_{\mu}$ is conserved up to higher-order terms. The background spin magnitude will then also be conserved up to $\mathcal{O}(s\epsilon,s^{2})$ terms at all times.

The energy and angular momentum of the spinning secondary are generated by the Killing symmetries of the background, so it is not surprising that their evolution averaged over the orbital time-scale balances the corresponding gravitational-wave fluxes [14, 55]

$$\left\langle\frac{\mathrm{d}E}{\mathrm{d}\tau}\right\rangle=-\mathcal{F}^{E}\,,\;\,\left\langle\frac{\mathrm{d}J_{z}}{\mathrm{d}\tau}\right\rangle=-\mathcal{F}^{J_{z}}\,.$$ (36)

However, the system has an additional degree of freedom in the form of the direction of the spin vector $s^{\mu}$. Specifically, it is conceivable that the self-torque drives the spin vector towards a more aligned, counter-aligned, or orthogonal configuration with respect to the angular momentum of the orbit. In other words, we need to derive the evolution of the constant $s_{\parallel}$.

Using the definition of $s_{\parallel}$ from eq. (8) we express the time derivative as

$$\dot{s}_{\parallel}=\frac{\dot{l}_{\nu}s^{\nu}}{\sqrt{l_{\alpha}l^{\alpha}}}+\frac{l_{\nu}\dot{s}^{\nu}}{\sqrt{l_{\alpha}l^{\alpha}}}-\frac{l_{\nu}s^{\nu}}{(l_{\alpha}l^{\alpha})^{3/2}}l_{\sigma}\dot{l}^{\sigma}\,.$$ (37)

We now average the relation above over orbital timescales to obtain the secular contribution to the evolution. We also discard terms of order $\mathcal{O}(s^{2},\epsilon^{2})$ but keep terms of order $\mathcal{O}(s,\epsilon s,\epsilon)$ as is consistent with the order of the scheme.

After averaging, only the parallel part of spin remains since it can be seen from equation (22) that all the other components of $s^{\mu}$ are fully oscillating on the orbital timescale. Additionally, the third term can be written using $s_{\parallel}$

$$\langle\dot{s}_{\parallel}\rangle=\left\langle\frac{\dot{l}_{\nu}s_{\parallel}^{\nu}}{\sqrt{l_{\alpha}l^{\alpha}}}+\frac{l_{\nu}\dot{s}^{\nu}}{\sqrt{l_{\alpha}l^{\alpha}}}-\frac{s_{\parallel}}{l_{\alpha}l^{\alpha}}l_{\sigma}\dot{l}^{\sigma}\right\rangle\,.$$ (38)

The parallel part of the angular-momentum vector is expressed as $s_{\parallel}^{\mu}=s_{\parallel}l^{\mu}/\sqrt{l_{\alpha}l^{\alpha}}$, which yields

$$\langle\dot{s}_{\parallel}\rangle=\left\langle\frac{s_{\parallel}\dot{l}_{\nu}l^{\nu}}{l_{\alpha}l^{\alpha}}+\frac{l_{\nu}\dot{s}^{\nu}}{\sqrt{l_{\alpha}l^{\alpha}}}-\frac{s_{\parallel}}{l_{\alpha}l^{\alpha}}l_{\sigma}\dot{l}^{\sigma}\right\rangle\,.$$ (39)

The first and third terms above subtract. Finally, we substitute the self-torque $\dot{s}^{\mu}$ from Eq. (34) to obtain

$$\begin{split}&\langle\dot{s}_{\parallel}\rangle=\left\langle\frac{l_{\nu}\dot{s}^{\nu}}{\sqrt{l_{\alpha}l^{\alpha}}}\right\rangle\\ &=-\frac{1}{2}\left\langle\frac{-l^{\mu}\frac{{\rm D}}{\differential\tau}b_{\mu}+l^{\mu}h_{\mu\beta;\rho}u^{\beta}s^{\rho}-l^{\rho}h_{\alpha\beta;\rho}u^{\alpha}s^{\beta}}{\sqrt{l_{\alpha}l^{\alpha}}}\right\rangle,\end{split}$$ (40)

$$b^{\mu}\equiv\quantity(g^{\mu\alpha}s^{\beta}-s^{\mu}\quantity(g^{\alpha\beta}+u^{\alpha}u^{\beta}))h_{\alpha\beta}\,.$$ (41)

Because $\dot{l}^{\mu}=\order{\epsilon}+\order{s}$ and $b^{\mu}=\order{\epsilon s}$, the first term together with the denominator can be written as a total derivative that does not contribute to the average. Additionally, since $s_{\parallel}^{\mu}$ and $l^{\mu}$ are colinear, the second and third terms cancel under the average and we obtain $\langle\dot{s}_{\parallel}\rangle=0+\mathcal{O}(\epsilon^{2},s^{2})$.

In conclusion, the leading-order adiabatic evolution of the spinning secondary orbit will be only due to the decay of $E$ and $J_{z}$ as given by Eq. (36), and $s$ and $s_{\parallel}$ can be treated simply as constants for the purposes of 1PA inspirals.

This derivation holds for generic orbits in Kerr spacetime and extends the same result for circular orbits in Schwarzschild spacetime in [55]. This is because $l_{\mu}\equiv Y_{\mu\nu}u^{\nu}$ is parallel-transported also along Kerr geodesics and thus all the derivation steps above apply without any change also for the motion near spinning primary black holes.

IV Gravitational-wave fluxes

IV.1 Teukolsky formalism

For the calculation of the PN expansion of GW fluxes in the framework of black hole perturbation theory, we use a similar approach to the one we used in Refs. [54, 67, 56] where we solved the Teukolsky equation in the frequency domain. Because the radial motion is periodic, the strain at infinity $h=h_{+}-ih_{\crossproduct}$ can be written as a sum over $l$, $m$ multipoles and harmonic modes $n$ and $j$ as

$$h=\frac{2}{r}\sum_{lmnj}\frac{C^{+}_{lmnj}}{\omega_{mnj}^{2}}{}_{-2}Y_{lm}(\theta)e^{-i\omega_{mnj}(t-r^{\ast})+im\phi}\;.$$ (42)

where we sum over $2\leq l<+\infty$, $-l\leq m\leq l$, $-\infty<n<+\infty$, and $-1<j<1$, $C^{+}_{lmnj}$ are the Teukolsky amplitudes at infinity, $\omega_{mnj}=m\Omega_{\phi}+n\Omega_{r}+j\Omega_{\psi}$ is the frequency of given mode, ${}_{-2}Y_{lm}(\theta)$ is spin-weighted spherical harmonic, $(t,r,\theta,\phi)$ are the coordinates of the observer and $r^{\ast}=r+2M\log(r/(2M)-1)$ is the tortoise coordinate.

The orbit-averaged energy and angular momentum fluxes to infinity can be expressed as sums over the $l$, $m$, $n$, $j$ modes in the form

$$\mathcal{F}^{E} =\sum_{lmnj}\frac{\absolutevalue{C^{+}_{lmnj}}^{2}}{4\pi\omega_{mnj}^{2}}\;,$$ (43a)

$$\mathcal{F}^{J_{z}} =\sum_{lmnj}\frac{m\absolutevalue{C^{+}_{lmnj}}^{2}}{4\pi\omega_{mnj}^{3}}\;.$$ (43b)

Because the amplitudes for $j=\pm 1$ are proportional to $s_{\perp}=\sqrt{s^{2}-s_{\parallel}^{2}}$ and for $j=0$ are independent of $s_{\perp}$, to linear order in spin the fluxes are independent of $s_{\perp}$ and we can sum only over $l$, $m$, and $n$ with $j=0$ [41, 56]. Therefore, we will write $C^{+}_{lmn}\equiv C^{+}_{lmn0}$ and $\omega_{mn}\equiv\omega_{mn0}$. Furthermore, as discussed later in this Section, the horizon fluxes are of higher PN order and we do not consider them here.

The asymptotic amplitudes can be found from the integral over the radial phase $q^{r}$

$$C^{+}_{lmn}=\frac{1}{W_{lmn}\Upsilon^{t}}\int_{0}^{2\pi}I^{+}_{lmn}(q^{r})e^{i\psi_{mn}(q^{r})}\differential q^{r}\,,$$ (44)

where

$$I^{+}_{lmn} =r^{2}\sum_{ab}\quantity(B^{0}_{ab}F^{ab}_{lmn}+B^{r}_{ab}\partialderivative{F^{ab}_{lmn}}{r})\,,$$ (45)

$$\psi_{mn}(q^{r}) =\omega_{mn}\Delta t(q^{r})-m\Delta\phi(q^{r})+nq^{r}\,,$$ (46)

with the sum over Kinnersley tetrad legs $ab=nn,n\overline{m},\overline{m}\overline{m}$. Note that we have rearranged the expression for $I^{+}_{lmn}$ in Eq. (52) from [56] and introduced quantities

$$F^{ab}_{lmn} =\sum_{i=0}^{I_{ab}}(-1)^{i}f^{(i)}_{ab}\derivative[i]{R^{-}_{lmn}}{r}\,,$$ (47)

$$B^{0}_{ab} =A^{\rm m}_{ab}+A^{\rm d}_{ab}+i\quantity(\omega_{mn}B^{t}_{ab}-mB^{\phi}_{ab})\,.$$ (48)

The functions $F^{ab}_{lmn}$ depend on the spin-weighted spherical harmonics ${}_{2}Y_{lm}$ through the functions $f^{(i)}_{ab}$ defined in Eqs. (B4) in [56], and on the solution of homogeneous radial Teukolsky equation $R^{-}_{lmn}\equiv R^{-}_{lm\omega_{mn}}$ satisfying purely outgoing boundary condition at the horizon (sometimes called the “in” solution). The quantities $A^{\rm m,d}_{ab}$ and $B^{\mu}_{ab}$ are calculated from the trajectory, the four-velocity, and the spin tensor and are defined in Eqs. (49) of [56].

Similarly to Ref. [67], we expand the expression for $C^{+}_{lmn}$ in the secondary spin $s_{\parallel}$. However, here the amplitudes and fluxes are expanded with fixed orbital parameters $p$ and $e$ as opposed to fixed orbital frequencies $\Omega_{r,\phi}$. The linear-in-spin part of the amplitude can be written as

$$\delta C^{+}_{lmn}=-\quantity(\frac{\delta\Upsilon^{t}}{\Upsilon^{t}_{\rm(g)}}+\frac{\partial_{\omega}W_{lmn}\delta\omega_{mn}}{W^{\rm(g)}_{lmn}})\frac{C^{{\rm(g)}+}_{lmn}}{W^{\rm(g)}_{lmn}\Upsilon^{t}_{\rm(g)}}+\frac{1}{\Upsilon^{t}_{\rm(g)}}\int_{0}^{2\pi}\quantity(\delta I^{+}_{lmn}+iI^{{\rm(g)}+}_{lmn}\delta\psi_{mn})e^{i\psi_{mn}^{\rm(g)}(q^{r})}\differential q^{r}\,,$$ (49)

where

$$\delta I^{+}_{lmn} =\frac{r_{\rm(g)}^{2}}{W^{\rm(g)}_{lmn}}\quantity(\quantity(\frac{2\delta r}{r_{\rm(g)}}A^{\rm m}_{{\rm(g)}ab}+\delta B^{0}_{ab})F^{{\rm(g)}ab}_{lmn}+\quantity(B^{r}_{ab}+A^{\rm m}_{{\rm(g)}ab}\delta r)\partialderivative{F^{{\rm(g)}ab}_{lmn}}{r}+A^{\rm(m)}_{{\rm(g)}ab}\partialderivative{F^{{\rm(g)}ab}_{lmn}}{\omega_{mn}}\delta\omega_{mn})\,,$$ (50)

$$\delta\psi_{mn} =\delta\omega_{mn}\Delta t_{\rm(g)}+\omega^{\rm(g)}_{mn}\delta\Delta t$$ (51)

with

$$\delta B^{0}_{ab} =\delta A^{\rm m}_{ab}+A^{\rm d}_{ab}+i\quantity(\omega B^{t}_{ab}-mB^{\phi}_{ab})\,,$$ (52)

$$\delta A^{\rm m}_{ab} =u_{\mu}u_{\nu}\partial_{r}(\lambda^{\mu}_{a}\lambda^{\nu}_{b})\delta r+2u_{\mu}\delta u_{\nu}\lambda^{(\mu}_{a}\lambda^{\nu)}_{b}\,,$$ (53)

where $\lambda^{\mu}_{a}$ are the legs of the Kinnersley tetrad.

Then, the fluxes can be separated into the geodesic and spin part as $\mathcal{F}^{E,J_{z}}=\mathcal{F}^{E,J_{z}}_{\rm(g)}+s_{\parallel}\delta\mathcal{F}^{E,J_{z}}/M$ to obtain

$$\delta\mathcal{F}^{E} =\sum_{lmn}\frac{\omega_{mn}\Re{C^{{\rm(g)}+}_{lmn}\overline{\delta C^{+}_{lmn}}}-\absolutevalue{C^{{\rm(g)}+}_{lmn}}^{2}\delta\omega_{mn}}{2\pi\omega_{mn}^{3}}\,,$$ (54a)

$$\delta\mathcal{F}^{J_{z}} =\sum_{lmn}m\frac{\omega_{mn}2\Re{C^{{\rm(g)}+}_{lmn}\overline{\delta C^{+}_{lmn}}}-3\absolutevalue{C^{{\rm(g)}+}_{lmn}}^{2}\delta\omega_{mn}}{4\pi\omega_{mn}^{4}}\,.$$ (54b)

IV.2 PN expansion of the fluxes

The geodesic amplitude $C^{{\rm(g)}+}_{lmn}$ and the linear-in-spin part $\delta C^{+}_{lmn}$ can be calculated as a PN series and series in $e$ by substituting the expansions of the trajectory from Section II.2. However, we also need a weak-field and low-velocities expansion of the radial function $R^{-}_{lmn}$ and Wronskian $W_{lmn}$. This has been done in [33, 35] where these quantities were expanded in $z=\omega r=\order{v}$ and $\epsilon=2M\omega=\order{v^{3}}$ (see [23] for review). Therefore, after substituting these variables, we obtain the expansion of the function $R^{-}_{lmn}(r(q^{r}))$ in $v$ and $e$.

After the expansion in $v$ and $e$, the integrals for the geodesic part $C^{{\rm(g)}+}_{lmn}$ in Eq (44) and for the linear-in-spin part $\delta C^{+}_{lmn}$ (49) consist of a finite Fourier series in $q^{r}$; therefore, they are trivial to integrate. In this way, we obtain the amplitudes with their linear-in-spin parts from which we calculate the fluxes and their linear-in-spin parts.

In the PN approximation and after expansion in eccentricity, the sums over $l$, $m$, and $n$ in the geodesic fluxes (43) and their linear-in-spin parts (54) are finite, since higher terms contribute only to higher order in $v$ and $e$. Unlike the geodesic parts of the $l$, $m$ multipoles of the fluxes, which start at $(l-2)$PN order for even $l+m$ and at $(l-1)$PN order for odd $l+m$, the linear-in-spin parts start at $(l-1/2)$PN order for both even and odd $l+m$. Since the linear-in-spin parts of the fluxes start at 1.5PN order, which corresponds to the spin-orbit coupling, we need to expand them to 3.5PN order NLO to obtain a 5PN expansion. Therefore, the fluxes and their linear-in-spin parts are summed over $2\leq l\leq 5$ and $-l\leq m\leq l$ to obtain the geodesic fluxes in the 3.5PN order and the linear-in-spin parts to 5PN order. Because the $n$-modes of the fluxes $\mathcal{F}_{lmn}$ behave as $\order{e^{2n}}$ and thanks to the symmetry $\mathcal{F}_{l,m,n}=\mathcal{F}_{l,-m,-n}$, we sum over $n$ in the range $-m<n\leq 5$ to obtain expansion to $e^{10}$.

Therefore, when the geodesic fluxes are completed to the 5PN order from, e.g., [48], we obtain the full 5PN energy and angular momentum fluxes from a spinning body orbiting a Schwarzschild black hole up to linear order in spin. Note that during the calculation of the linear-in-spin part, nonzero terms appear in the 1PN position, which cancel out and the series then start at 1.5PN term. Therefore, because of the subtraction of the leading term, the trajectory must be expanded to one order higher than is the order of the final series. Furthermore, because the horizon fluxes for nonspinning secondary in Schwarzschild spacetime start at 4PN order, the linear-in-spin contribution to the horizon fluxes starts at 5.5PN order. Thus, we do not need to consider them here.

As discussed in [45], each PN term contains a factor with a certain power of $1-e^{2}$. When the fluxes are expressed using the parameter $v$, this factor reads $(1-e^{2})^{3/2}$ for all orders and can be factored out [43]. The resulting linear-in-spin parts of the energy and angular-momentum fluxes have the form

$$\begin{split}&\delta\mathcal{F}^{E}=\mathcal{F}^{E}_{\rm N}\quantity(1-e^{2})^{3/2}\Bigg{[}\delta f_{3}v^{3}+\delta f_{5}v^{5}+\delta f_{6}v^{6}+\delta f_{7}v^{7}+\delta f_{8}v^{8}+\left(\delta f_{9}+\delta f_{9}^{\log v}\left(\gamma-\frac{35\pi^{2}}{107}+\log v\right)\right)v^{9}\\ &\phantom{\delta\mathcal{F}^{E}=}+\delta f_{10}v^{10}+\order{v^{11}}\Bigg{]}\,,\end{split}$$ (55)

$$\begin{split}&\delta\mathcal{F}^{J_{z}}=\mathcal{F}^{J_{z}}_{\rm N}\quantity(1-e^{2})^{3/2}\Bigg{[}\delta g_{3}v^{3}+\delta g_{5}v^{5}+\delta g_{6}v^{6}+\delta g_{7}v^{7}+\delta g_{8}v^{8}+\quantity(\delta g_{9}+\delta g_{9}^{\log v}\quantity(\gamma-\frac{35\pi^{2}}{107}+\log v))v^{9}\\ &\phantom{\delta\mathcal{F}^{J_{z}}=}+\delta g_{10}v^{10}+\order{v^{11}}\Bigg{]}\,,\end{split}$$ (56)

where

$$\mathcal{F}^{E}_{\rm N}=\frac{32}{5}\quantity(\frac{\mu}{M})^{2}v^{10}\,,\quad\mathcal{F}^{J_{z}}_{\rm N}=\frac{32}{5}\frac{\mu^{2}}{M}v^{7}\,,$$ (57)

are the Newtonian fluxes, and $\gamma$ is the Euler–Mascheroni constant. The $\delta f_{i}(e),\delta g_{i}(e)$ are functions of eccentricity similar to the enhancement functions of Peters & Mathews and can be found in Appendix A.

Similarly to the geodesic part, we were able to resum the leading term $\delta f_{3}$ ($\delta g_{3}$), the 1PN and 2PN contributions $\delta f_{5}$, $\delta f_{7}$ ($\delta g_{5}$, $\delta g_{7}$) and the logarithmic term $\delta f_{9}^{\log v}$ ($\delta g_{9}^{\log v}$) and write them in closed form.

After expansion to $e^{12}$ and factorization of $(1-e^{2})^{3/2}$, in the functions $\delta f_{3}$, $\delta g_{3}$, $\delta f_{9}^{\log(v)}$, and $\delta g_{9}^{\log(v)}$, some of the last terms vanished. In particular, the series $\delta f_{3}$ ended at $e^{6}$, $\delta g_{3}$ at $e^{4}$, $\delta f_{9}^{\log(v)}$ at $e^{10}$, and $\delta g_{9}^{\log(v)}$ at $e^{8}$. Furthermore, after subtracting terms proportional to $\sqrt{1-e^{2}}$ from $\delta f_{5}$, $\delta g_{5}$, $\delta f_{7}$ and $\delta g_{7}$, in the remaining series, similarly, some terms vanished. These series ended at $e^{8}$, $e^{6}$, $e^{10}$, and $e^{8}$, respectively. Therefore, we did not verify some of the resummations to all orders in eccentricity, but we assume that they are true from the similar behavior of the geodesic part [46, 47].

Figure 1: Coefficients in the PN expansion and eccentricity expansion of the linear part of the energy flux $\delta\mathcal{F}^{E}$ from Eq. (55).

Figure 2: Coefficients in the PN expansion and eccentricity expansion of the linear part of the angular momentum flux $\delta\mathcal{F}^{J_{z}}$ from Eq. (56).

In Figures 1 and 2 we plot the coefficients of the PN series of the linear-in-spin parts of the energy and angular momentum flux. Each line shows the coefficients of the series in $e$ for a given PN order. The coefficients seem to decrease with eccentricity for all PN orders, which suggests that the truncation of the eccentricity series does not cause a large error. However, it may be improved by fitting and subtracting unknown terms proportional to $\sqrt{1-e^{2}}$ or $\log(1+\sqrt{1-e^{2}})$, which we know to appear in geodesic fluxes [46, 47].

Fluxes can also be expressed using the gauge-invariant quantity $x=(M\Omega_{\phi})^{2/3}$. Then, they can be linearized as

$$\mathcal{F}(x,e,s_{\parallel})=\mathcal{F}_{\rm(g)}(x,e)+s_{\parallel}\evaluated{\delta\mathcal{F}}_{x,e}(x,e)\,,$$ (58)

where the linear-in-spin part with fixed $x$ and $e$ can be obtained from the linear-in-spin part with fixed $v$ and $e$ as

$$\evaluated{\delta\mathcal{F}}_{x,e}(x,e)=\evaluated{\delta\mathcal{F}}_{v,e}(v_{\rm(g)}(x,e),e)+\partialderivative{\mathcal{F}_{\rm g}(v,e)}{v}\delta v\,,$$ (59)

where $v_{\rm(g)}(x,e)$ and $\delta v(x,e)$ can be found from $\Omega_{\phi}$ as inverse series of the series $x=(M(\Omega_{\phi}^{\rm(g)}(v(x,e),e)+s_{\parallel}\delta\Omega_{\phi}(v(x,e),e)))^{2/3}/M$. The results for $\delta v$ and $\evaluated{\delta\mathcal{F}}_{x,e}(x,e)$ are given in the Supplemental material.

After this transformation, the energy flux for zero eccentricity agrees with the results of Nagar et al. [30], where the PN expansion of energy fluxes from spinning bodies on circular orbits of a Schwarzschild black hole was derived.

To verify our results, we compare them with the results of Henry and Khalil [19], where the energy and angular momentum fluxes from eccentric spinning binaries were calculated using the PN theory. Their results to 3PN and $e^{8}$ are given as functions of $x$ and the time eccentricity $e_{t}$, which is used in the quasi-Keplerian parametrization described in Eqs. (2.26) in [19]. Therefore, we had to transform their fluxes to functions of $x$ and $e$ using a relation between $e_{t}$ and $e$ derived in Appendix B. After the transformation, the linear parts of the energy and angular momentum fluxes derived in [19] agree with our results up to the 3PN order, $e^{8}$ and the first order in the mass ratio.

Figure 3: Relative difference between the PN expansion of the linear-in-spin part of the energy (top) and angular momentum (bottom) flux $\delta\mathcal{F}_{\text{PN}}$ and the fully relativistic value of $\delta\mathcal{F}_{\text{num}}$ for different eccentricities. The dashed lines show dependence $\order{p^{-4}}$ which should be the order of the error.

To further validate our results, we compare the PN series with fully relativistic numerically calculated linear parts of the fluxes calculated in [67]. We calculate the relative errors

$$\absolutevalue{1-\frac{\delta\mathcal{F}_{\rm PN}}{\delta\mathcal{F}_{\rm num}}}$$ (60)

for $\delta\mathcal{F}^{E}$ and $\delta\mathcal{F}^{J_{z}}$ and plot them in Figure 3 as functions of $p$ for different values of the eccentricity. These plots verify that the relative difference decreases with increasing $p$. For comparison, we also plot the behavior $p^{-4}$ since it is the behavior of the first neglected PN term (because the fluxes are expanded to 3.5PN orders NLO). The relative differences seem to decrease with higher power of $p$ which is probably caused by the smaller magnitude of the 4PN NLO term compared to the 4.5PN NLO term. For higher $p$, the relative difference is dominated by the interpolation error of the numerical fluxes.

V Flux-driven inspirals

Once we obtained the energy and angular momentum fluxes, we can calculate the inspiral, i.e. the evolution of the orbital parameters. As discussed in Section III, the fluxes of energy and angular momentum are sufficient to calculate the evolution of $p$ and $e$ since $s_{\parallel}$ is conserved.

V.1 Analytical integration of quasi-circular inspirals

To obtain a first understanding of the convergence of the PN expansion, it is useful to examine the dynamics analytically. This is achieved by using the PN expansion of Schwarzschild geodesic fluxes as obtained in Refs [48, 50] and implemented in the PostNewtonianSelfForce Mathematica package [68] along with the spin fluxes derived here. While it is in principle possible to analytically integrate the dynamics at generic eccentricity, the symbolic computations become prohibitively expensive. For this reason, we restricted ourselves to quasicircular inspirals for the analytical convergence exploration (thus essentially restricting ourselves to the earlier flux formulas of Refs. [35, 38, 30]).

In that case, we can evolve the inspiral only in terms of the PN expansion parameter $v$. Furthermore, we can reparametrize the evolution with the azimuthal phase $\phi$

$$\begin{split}&\frac{\mathrm{d}v}{\mathrm{d}\phi}\Big{|}_{e=0}=\frac{1}{2}\sqrt{\frac{M}{p^{3}}}\left(\frac{\mathrm{d}E}{\mathrm{d}p}\right)^{-1}\!\mathcal{F}^{E}\,\Omega^{\phi}\\ &\phantom{\frac{\mathrm{d}v}{\mathrm{d}\phi}\Big{|}_{e=0}}=\frac{\mathrm{d}v}{\mathrm{d}\phi}\Bigg{|}_{\mathrm{(g)}}+\frac{s_{\parallel}}{M}\frac{\mathrm{d}v}{\mathrm{d}\phi}\Bigg{|}_{\mathrm{(s)}}\,,\end{split}$$ (61)

where the relations for the fluxes, frequencies, and $E(p)$ are evaluated at $e=0$ and receive $\mathcal{O}(s_{\parallel}/M)$ corrections as described above. For the quasi-circular inspiral, one could equivalently use the $\mathcal{F}^{J_{z}}$ flux and a $J_{z}(p)$ relation due to the identity $E=J_{z}\Omega^{\phi}$ at zero eccentricity.

Then, we get the equation for the evolution of azimuthal phase as a function of $v$ by $\phi^{\prime}(v)=(\mathrm{d}v/\mathrm{d}\phi)^{-1}$ as

$$\begin{split}&\phi^{\prime}(v)=\frac{5M}{32\mu v^{6}}\Bigg{[}1+\frac{743}{336}v^{2}-4\pi v^{3}+\frac{3}{2}\frac{s_{\parallel}}{M}v^{3}+\frac{3058673}{1016064}v^{4}-\frac{7729\pi}{672}v^{5}+\frac{743}{224}\frac{s_{\parallel}}{M}v^{5}+\Phi_{6}v^{6}-6\pi\frac{s_{\parallel}}{M}v^{6}\\ &\phantom{\phi^{\prime}(v)=}-\frac{15419335\pi}{1016064}v^{7}+\frac{3058673}{677376}\frac{s_{\parallel}}{M}v^{7}+\Phi_{8}v^{8}-\frac{7729\pi}{448}\frac{s_{\parallel}}{M}v^{8}+\Phi_{9}v^{9}+\Phi_{s9}\frac{s_{\parallel}}{M}v^{9}+\Phi_{10}v^{10}-\frac{15419335\pi}{677376}\frac{s_{\parallel}}{M}v^{10}\Bigg{]},\end{split}$$ (62)

where the $\Phi$ coefficients are defined as

$$\Phi_{6}\equiv-\frac{10817850546611}{93884313600}+\frac{1712}{105}(\gamma+\log v)+\frac{32}{3}\pi^{2}+\frac{3424}{105}\log 2\,,$$ (63)

$$\Phi_{8}\equiv-\frac{2500489942240134443}{3690780136243200}+\frac{9203}{210}(\gamma+\log v)+\frac{9049}{252}\pi^{2}+\frac{50551}{882}\log 2+\frac{47385}{1568}\log 3\,,$$ (64)

$$\Phi_{9}\equiv\pi\left(\frac{90036665674763}{187768627200}-\frac{6848}{105}(\gamma+\log v)-\frac{64}{3}\pi^{2}-\frac{13696}{105}\log 2\right),$$ (65)

$$\Phi_{9s}\equiv-\frac{10270192050611}{62589542400}+\frac{856}{35}(\gamma+\log v)+16\pi^{2}+\frac{1712}{35}\log 2\,,$$ (66)

$$\begin{split}&\Phi_{10}\equiv-\frac{1417220168422461061151}{505226791983513600}+\frac{6470582647}{110020680}(\gamma+\log v)+\frac{578223115}{12192768}\pi^{2}+\frac{53992839431}{220041360}\log 2\\ &\phantom{\Phi_{10}=}-\frac{5512455}{87808}\log 3\,.\end{split}$$ (67)

This relation can then be easily integrated term by term to obtain the change in $\phi$ between two referential values of $v$. We can take the end of the inspiral to be at the innermost stable circular orbit, which is at

$$v_{\rm ISCO}=\sqrt{\frac{1}{6}}+\frac{1}{18}\frac{s_{\parallel}}{M}\,.$$ (68)

Furthermore, we want to parameterize the initial condition by a referential initial frequency where the signal enters the band, $\Omega_{\phi}=\Omega^{\phi}_{\rm i}$. The perturbative inversion of the relation $\Omega_{\phi}(v)$ yields $v=\Omega_{\phi}^{1/3}+s_{\parallel}\Omega_{\phi}^{4/3}/(2M)$. As a result, we get the inspiral phase as

$$\begin{split}&\Delta\phi_{\rm ISCO}=\frac{M}{32\mu\omega_{\mathrm{i}}^{5/3}}\Bigg{[}1+\frac{3715}{1008}\omega_{\mathrm{i}}^{1/3}-10\pi\omega_{\mathrm{i}}+\frac{15293365}{1016064}\omega_{\mathrm{i}}^{4/3}+\Delta_{5/3}\omega_{\mathrm{i}}^{5/3}+\Delta_{2}\omega_{\mathrm{i}}^{2}+\frac{77096675\pi}{2032128}\omega_{\mathrm{i}}^{7/3}\\ &\phantom{\phi_{\rm ISCO}-\phi_{\rm I}=}+\Delta_{8/3}\omega_{\mathrm{i}}^{8/3}+\Delta_{3}\omega_{\mathrm{i}}^{3}+\Delta_{10/3}\omega_{\mathrm{i}}^{10/3}\Bigg{]}+\frac{s_{\parallel}}{M}\frac{5}{128\omega_{\mathrm{i}}^{2/3}}\Bigg{[}1+\Delta_{s2/3}\omega_{\mathrm{i}}^{2/3}+32\pi\omega_{\mathrm{i}}-\frac{15293365}{1016064}\omega_{\mathrm{i}}^{4/3}\\ &\phantom{\phi_{\rm ISCO}-\phi_{\rm I}=}+\frac{7729\pi}{168}\omega_{\mathrm{i}}^{5/3}+\Delta_{s2}\omega_{\mathrm{i}}^{2}+\frac{3083867\pi}{63504}\omega_{\mathrm{i}}^{7/3}\Bigg{]}+\mathcal{O}(\omega_{\mathrm{i}}^{6/3})\,,\end{split}$$ (69)

where we have defined $\omega_{\mathrm{i}}\equiv M\Omega^{\phi}_{\mathrm{i}}$ as the initial frequency in units of $M$. The $\Delta$ coefficients then read

$$\begin{split}&\Delta_{5/3}\equiv-\frac{223791298249051766200631}{163693480602658406400\sqrt{6}}+\frac{377580814447}{3960744480\sqrt{6}}\gamma+\frac{383368458940043\pi}{5407736463360}-\frac{428\pi\gamma}{189}+\frac{28366605835\pi^{2}}{438939648\sqrt{6}}\\ &\phantom{\Delta_{5/3}=}-\frac{20\pi^{3}}{27}+\frac{136763321753}{990186120\sqrt{6}}\log 2+\frac{102239\pi}{4032}\log 2-\frac{5198370032377}{126743823360\sqrt{6}}\log 3+\frac{51643\pi}{1728}\log 3\\ &\phantom{\Delta_{5/3}=}+\frac{38645}{672}\pi\log(\omega_{\mathrm{i}}^{1/3})\,,\end{split}$$ (70)

$$\Delta_{2}\equiv\frac{12348611926451}{18776862720}-\frac{1712}{21}\left(\gamma+\log(\omega_{\mathrm{i}}^{1/3})\right)-\frac{160\pi^{2}}{3}-\frac{3424}{21}\log 2\,,$$ (71)

$$\Delta_{8/3}\equiv\frac{2554404624135128353}{2214468081745920}-\frac{9203}{126}\left(\gamma+\log(\omega_{\mathrm{i}}^{1/3})\right)-\frac{45245\pi^{2}}{756}-\frac{252755}{2646}\log 2-\frac{78975}{1568}\log 3\,,$$ (72)

$$\Delta_{3}\equiv\pi\left[-\frac{93098188434443}{150214901760}+\frac{1712}{21}\left(\gamma+\log(\omega_{\mathrm{i}}^{1/3})\right)+\frac{80\pi^{2}}{3}+\frac{3424}{21}\log 2\right]\,,$$ (73)

$$\begin{split}&\Delta_{10/3}\equiv\frac{474387630222958367413}{168408930661171200}-\frac{6470582647}{110020680}\left(\gamma+\log(\omega_{\mathrm{i}}^{1/3})\right)-\frac{578223115\pi^{2}}{12192768}-\frac{53992839431}{220041360}\log 2\\ &\phantom{\Delta_{10/3}=}+\frac{5512455}{87808}\log 3\,,\end{split}$$ (74)

$$\begin{split}&\Delta_{s2/3}\equiv-\frac{5881840409277979019197}{245540220903987609600}-\frac{310713839464837\pi}{5069752934400\sqrt{6}}+\frac{8765086219\pi^{2}}{1975228416}-\frac{32}{81}\sqrt{\frac{2}{3}}\pi^{3}\\ &\phantom{\Delta_{s2/3}=}+\frac{(112087291999-7175443968\sqrt{6}\pi)\gamma}{17823350160}+\frac{43923447107}{17823350160}\log 2-\frac{1712\pi}{2835}\sqrt{\frac{2}{3}}\log\left(\frac{8}{3}\right)\\ &\phantom{\Delta_{s2/3}=}-\frac{5159694245689}{570347205120}\log 3-\frac{743}{56}\log(\omega_{\mathrm{i}}^{1/3})\,,\end{split}$$ (75)

$$\Delta_{s2}\equiv\frac{10747149910451}{26824089600}-\frac{856}{15}\left(\gamma+\log(\omega_{\mathrm{i}}^{1/3})\right)-\frac{112\pi^{2}}{3}-\frac{1712}{15}\log 2\,.$$ (76)

Let us now plug in numbers corresponding to LISA EMRIs into this formula to investigate its convergence (we will use the same numbers in Section V.3). We use a primary mass $10^{6}M_{\odot}$, initial frequency of the $l=2,m=2$ mode equal to $f_{m=2}=2\Omega^{\phi}_{\rm i}/(2\pi)=1\rm\,mHz$, and a secondary mass of $\mu=100M_{\odot}$ and thus $\epsilon=10^{-4}$. The spin is also picked as $s_{\parallel}=100M_{\odot}$, which corresponds to a maximally spinning and aligned secondary black hole. We evaluate each term separately and summarize the results in Table 1. In general alignment with the observations made by Burke et al. [51] and Isoyama et al. [69] for non-spinning secondaries, we see that even though the geodesic adiabatic terms are far from sub-radian accuracy at 5PN, the spin terms are suppressed by a mass ratio factor and have converged well below radians in this scenario. This supports the hybrid approach which we will use for the evolution of eccentric inspirals in the next Section.

Table 1: Contributions of various terms to the final phase a of a quasi-circular EMRI entering band at $1$ mHz at primary mass $M=10^{6}M_{\odot}$ and secondary mass $\mu=100M_{\odot}$ with a maximally spinning and aligned secondary.

PN Ord.	Geodesic $\Delta\phi$	Spin $\Delta\phi$
0	325 402	0
1	74 454.5	0
1.5	-158 135	0.629199
2	18 877.6	0
2.5	-47 387.3	-2.2748
3	6578.54	0.978466
3.5	2312.36	-0.0365018
4	2420.03	0.087333
4.5	-1172.24	-0.0000844637
5	688.9	0.00572313

V.2 Evolution equations for eccentric inspirals

To evolve the orbital parameters $p,e$, we must first use Eq. (36) to derive their (average) time derivatives. From the relation between the constants of motion $E$ and $J_{z}$ and the $p,e$, the evolution equations can be written as

$$\matrixquantity(\dot{p}\\ \dot{e})\equiv\matrixquantity(\derivative{p}{t}\\ \derivative{e}{t})=-\mathbb{J}^{-1}\matrixquantity(\mathcal{F}^{E}\\ \mathcal{F}^{J_{z}})\,,$$ (77)

where $\mathbb{J}$ is the Jacobian

which is known analytically and can be found in Appendix C. Eq. (77) can be expanded in the secondary spin as

$$\matrixquantity(\dot{p}\\ \dot{e})=-\mathbb{J}_{\rm(g)}^{-1}\matrixquantity(\mathcal{F}_{\rm(g)}^{E}\\ \mathcal{F}_{\rm(g)}^{J_{z}})\\ +\frac{s_{\parallel}}{M}\quantity(\mathbb{J}_{\rm(g)}^{-1}\delta\mathbb{J}\mathbb{J}_{\rm(g)}^{-1}\matrixquantity(\mathcal{F}_{\rm(g)}^{E}\\ \mathcal{F}_{\rm(g)}^{J_{z}})-\mathbb{J}_{\rm(g)}^{-1}\matrixquantity(\delta\mathcal{F}^{E}\\ \delta\mathcal{F}^{J_{z}})),$$ (80)

where we used the relation for the derivative of inverse matrix. The geodesic energy and angular momentum fluxes $\mathcal{F}^{E,J_{z}}_{\rm(g)}$ can be written using the geodesic evolution of $p$ and $e$ as

$$\begin{split}&\matrixquantity(\dot{p}\\ \dot{e})=\matrixquantity(\dot{p}_{\rm(g)}\\ \dot{e}_{\rm(g)})\\ &\phantom{\matrixquantity(\dot{p}\\ \dot{e})=}+\frac{s_{\parallel}}{M}\quantity(\mathbb{J}_{\rm(g)}^{-1}\delta\mathbb{J}\matrixquantity(\dot{p}_{\rm(g)}\\ \dot{e}_{\rm(g)})-\mathbb{J}_{\rm(g)}^{-1}\matrixquantity(\delta\mathcal{F}^{E}\\ \delta\mathcal{F}^{J_{z}})).\end{split}$$ (81)

Note that the first term is associated with the adiabatic term while the second term contributes only to the postadiabatic term. Therefore, the requirements for the accuracy of the first term are much higher than the requirements for the accuracy of the second term. Thus, we can use the PN expansion of the linear-in-spin parts of the energy and angular momentum fluxes. The geodesic evolution of $p$ and $e$ in the fully relativistic regime was calculated numerically and subsequently interpolated on a grid in the $p$-$e$ plane in [67]. Therefore, the evolution equations we use in this work read

$$\matrixquantity(\dot{p}\\ \dot{e})=\matrixquantity(\dot{p}_{\rm(g)}^{\rm num}\\ \dot{e}_{\rm(g)}^{\rm num})\\ +\frac{s_{\parallel}}{M}\quantity(\mathbb{J}_{\rm(g)}^{-1}\delta\mathbb{J}\matrixquantity(\dot{p}_{\rm(g)}^{\rm num}\\ \dot{e}_{\rm(g)}^{\rm num})-\mathbb{J}_{\rm(g)}^{-1}\matrixquantity(\delta\mathcal{F}^{E}_{\rm PN}\\ \delta\mathcal{F}^{J_{z}}_{\rm PN}))\,,$$ (84)

where the superscript “num” means fully relativistic results, subscript “PN” denotes the PN expansion and the Jacobian and its $s_{\parallel}$-derivative are fully relativistic as well because they can contain some nontrivial behavior near the last stable orbit. The explicit form of the matrix product $\mathbb{J}_{\rm(g)}^{-1}\delta\mathbb{J}$ can be found in Appendix C.

After the evolution of $p(t)$ and $e(t)$ is obtained, the inspiral waveform from two-timescale expansion can be calculated as

$$h=\frac{1}{r}\sum_{lmn}A_{lmn}(t)Y_{lm}(\theta)e^{-i\Phi_{mn}(t-r^{\ast})+im\phi}\,,$$ (85)

where the amplitude and phase read, respectively,

$$A_{lmn}(t) =\frac{2C^{+}_{lmn}(p(t),e(t))}{\omega_{mn}^{2}(p(t),e(t))}\,,$$ (86)

$$\Phi_{mn}(t) =m\Phi_{\phi}(t)+n\Phi_{r}(t)\,,$$ (87)

$$\Phi_{r,\phi}(t) =\int_{0}^{t}\Omega_{r,\phi}(p(u),e(u),s_{\parallel})\differential u\,.$$ (88)

V.3 LISA band inspirals

To verify the validity of the hybrid model (84) containing the PN expansions, in this Section we compare inspirals calculated using this model and a fully relativistic model for astrophysically relevant EMRIs that will be possible to detect with LISA. Similarly to Section V.1, the primary mass is chosen as $M=10^{6}M_{\odot}$. For this primary mass, the frequency of the dominant mode $l=m=2$, $n=0$ in the innermost stable circular orbit is $f_{20}=\omega_{20}/(2\pi)=2\Omega_{\phi}/(2\pi)=4.4$ mHz, which is close to the minimum of the LISA noise curve. The mass of the secondary and parallel spin is chosen as $\mu=s_{\parallel}=100M_{\odot}$. Therefore, the mass ratio is $\epsilon=10^{-4}$ and the secondary corresponds to a maximally spinning Kerr black hole. We evolved the inspirals in a range where the frequency is $2\Omega_{\phi}/(2\pi)\geq 1$ mHz.

The inspirals cannot be evolved all the way to the separatrix in our setup. This is caused by the fact that the grid, on which we interpolated the numerical fluxes, starts at finite distance from the separatrix, and, also, because some quantities linear in the secondary spin diverge there. Therefore, we need to choose a consistent condition to end the inspirals. In our setup, this condition reads

$$\frac{\derivative{\Omega_{r}}{t}}{\Omega_{r}^{2}}=10^{-2}\,,$$ (89)

which corresponds to a radial inverse adiabaticity parameter (we draw inspiration from a similar parameter defined in Ref. [27]). This quantity is small for adiabatic inspirals and grows near the separatrix where the two-timescale expansion breaks.

Figure 4: Adiabatic inspirals in the $p$-$e$ plane. The black line shows the separatrix $p=6+2e$. The two models are indistinguishable in this plot (see Figure 5 for phase differences).

In this setup, we found values of $p$ that satisfy the condition (89) for $e$ between $0.05$ and $0.4$ and evolved the inspirals backward using the fully relativistic model. The evolution was stopped when the condition $f=2\Omega_{\phi}<1$ mHz was reached. Then, we used the hybrid model (84) to evolve the inspirals from the end points of the previous calculation. In this way, we obtained two sets of evolutions of $p$ and $e$ with different models for comparison. The results in the $p$-$e$ plane are depicted in Fig. 4.

Figure 5: Absolute differences between the azimuthal (top) and radial (bottom) phases obtained from the PN model and fully relativistic model.

In the next step, we used the analytic formulas for the orbital frequencies $\Omega_{r}$ and $\Omega_{\phi}$ to calculate the phases $\Phi_{r}$ and $\Phi_{\phi}$ from Eq. (88). In Figure 5 we plot the absolute difference between the phases calculated with the hybrid and fully relativistic model. We can see that the phase differences are below unity.

Figure 6: Waveforms of the $+$ polarization of an inspiral ending at $e=0.1$ calculated with the hybrid model (blue) and with the fully relativistic (numerical) model (yellow). The inspiral is observed from the distance of 1 Gpc at the viewing angle $\theta=\pi/3$, $\phi=\pi/4$ in the source frame.

Figure 7: Waveforms of the $+$ polarization of an inspiral ending at $e=0.3$ calculated with the hybrid model (blue) and the fully relativistic (numerical) model (yellow). The distance and viewing angle are the same as in Figure 6.

To calculate both the evolution of the orbital parameters and the phases, we used NDSolve function in Mathematica with Adams’ method. The resulting time series $p(t)$, $e(t)$, and $\Phi_{r,\phi}(t)$ were then used to calculate the waveform using the FastEMRIWaveforms (FEW) package [59, 60, 61, 62]. The distance of the observer was chosen as 1 Gpc and the viewing angle was chosen as $\theta=\pi/3$, $\phi=\pi/4$ in the source frame. FEW calculates the waveforms (85) with the geodesic amplitudes $C^{{\rm(g)}+}_{lmn}$, which introduces $\order{s_{\parallel}/M}$ error in the amplitudes. However, this is not an issue since the requirements for the accuracy of the amplitudes are lower than the requirements for the accuracy of the phases [70]. In Figures 6 and 7 we show a comparison of the waveforms calculated with the two models for two inspirals ending at $e=0.1$ and $e=0.3$.

Figure 8: Mismatches between inspirals calculated with the hybrid model and the fully relativistic model for different final eccentricities $e_{f}$. The blue points show inspirals which end at $\dot{\Omega}_{r}/\Omega_{r}^{2}=10^{-2}$ while the yellow points show inspirals ending at $\dot{\Omega}_{r}/\Omega_{r}^{2}=10^{-3}$. Small differences between these cases indicate that the mismatches are almost independent of the ending criterion.

From the obtained waveforms we calculated the mismatch between the fully relativistic model and the hybrid model with PN linear-in-spin parts. Mismatch is defined from the overlap $\mathcal{O}$ as

$$\begin{split}\mathcal{M}(h_{1},h_{2})&\equiv 1-\mathcal{O}(h_{1},h_{2})\\ &=1-\frac{\langle h_{1},h_{2}\rangle}{\sqrt{\langle h_{1},h_{1}\rangle\langle h_{2},h_{2}\rangle}}\,,\end{split}$$ (90)

where $\langle\cdot,\cdot\rangle$ is the $L^{2}$ product. When the two waveforms are identical, the mismatch is zero. We plotted the mismatches for different final eccentricities in Figure 8. To test whether the ending criterion 89 influences the mismatches, we calculated the inspirals for two ending criterion, namely $\dot{\Omega}_{r}/\Omega_{r}^{2}=10^{-2}$ and $10^{-3}$ and compared them in Fig. 8. In this plot, we can see that the mismatch is consistent for the two ending criteria and it is lower than $10^{-2}$ for lower final eccentricities.

VI Discussion and outlooks

In the previous sections, we calculated the PN expansions of the trajectories of spinning bodies on eccentric orbits around Schwarzschild black holes. Then, we found the PN expansions of the energy and angular momentum fluxes from the aforementioned orbits. The linear-in-spin parts of the fluxes were then used in a hybrid model, where the subleading secondary-spin effects were analytically approximated by using the PN series. Mismatches between waveforms from the fully relativistic model and our hybrid model showed that for lower eccentricities the models are indistinguishable. This results shows that in some cases the linear-in-spin part of the fluxes can be approximated as an analytical PN series without the need to numerically calculate the fully relativistic contribution. However, to accurately assess the possible biases introduced by the model across the parameter space, a Fisher-matrix or Markov-chain Monte Carlo analysis such as those carried out by Burke et al. [51] and Piovano et al. [71] should be performed with this model.

Figure 8 shows that the mismatch is greater for inspirals with higher eccentricity. This could be improved by expanding to higher order in eccentricity or by finding exact (or arbitrary-order in eccentricity) formulas such as in [46, 47]. However, in Figure 4 we can see that the inspirals with higher eccentricity enter the LISA band in stronger field at lower $p$. Therefore, expanding to higher PN order may also improve accuracy. However, the computations at higher PN order increase in complexity. For example, modes with $n=-m$ and higher $l$ modes must be included starting from $5.5$PN and higher for the spin fluxes. Additionally, the horizon fluxes will be needed as well since they start at 4PN for the geodesic part and at 5.5PN for the linear-in-spin part.

Poor convergence of the PN series for higher eccentricities can be caused by the fact that the secondary body reaches a stronger field at the pericenter $r_{2}=Mp/(1+e)$ even for high $p$ (i.e. small $v$). However, the convergence of the series in $v$ is better than the convergence of the series in $x$ since at fixed $x$ the pericenter approaches zero when $e\rightarrow 1$ and the fluxes diverge there. This is connected to the cancellation of the divergent factors $(1-e^{2})^{-k/2}$ appearing in the $x$ series when it is reparametrized by $v$.

Nevertheless, other non-analytical terms of the type $(1-e^{2})^{k/2}$ with $k>0$ systematically appear in the flux series. The factorization of such terms on a case-to-case basis allowed us to resum the otherwise infinite $e$-series for a number of terms. However, at high PN orders this requires more and more terms in the $e$-series to be verified. Additionally, Figs 1 and 2 show that the higher-PN terms are clearly not as well converged as lower-order terms at $e^{10}$.

On the other hand, the requirement $2\Omega_{\phi}>10^{-3}\,\text{Hz}$ we impose for the starting point of inspirals works best only for circular or low-eccentricity inspirals; it may be too crude for the highly eccentric cases. This is because for higher eccentricities, higher $n$ modes are present, thus introducing higher-frequency harmonics into the spectrum, which then enter the LISA band earlier than our cutoff. Thus, a more sophisticated analysis of LISA mismatches of extended waveforms without such simplifications is needed.

Therefore, conclusions about highly eccentric inspirals should not be drawn from the results for quasicircular inspirals. We see this also in Table 1, where the contribution to the phase from the last, 5PN, term is of the order of $10^{-2}$. Such a truncation error would be sufficient for LISA waveforms, but this convergence property unfortunately does not generalize to eccentric inspirals. We can extrapolate our observations using secondary spin even to the hybrid model of Burke et al. [51], where 3PN approximations of second-order fluxes and conservative self-force were used in quasi-circular inspirals of nonspinning binaries with encouraging results (see also the earlier work of Isoyama et al. [69]). We do not expect these encouraging results to generalize to eccentric inspirals. What is more, we do not expect even 5PN-$e^{10}$ expansions of the second-order fluxes and conservative self-force to be sufficient for LISA parameter estimates of highly eccentric inspirals in hybrid models.

How could the results of this paper be further generalized or expanded? A possible extension would be to calculate the PN expansion of energy and angular momentum fluxes from generic orbits of spinning bodies in Kerr spacetime. This could be achieved by expanding the equations of motion obtained from the Hamilton-Jacobi equation [64] in a PN series and solving them order-by-order. We are already preparing a work in which we solve for the fundamental frequencies of motion of the spinning particle in Kerr in closed form and expand them in a PN series [72] (see also [73]), but the full trajectories pose more of a technical challenge.

Nonetheless, generic orbits of spinning test particles in Kerr are parametrized by one additional constant of motion, the Rüdiger (Carter-like) constant $K_{R}$. Hence, to calculate the inspirals, the evolution of this constant must first be derived, similarly to the evolution of $s_{\parallel}$ presented here in Section III. Until then, one can only evolve equatorial inspirals as done in [67]. Another loophole possibility to drive inspirals without the need for the evolution of $K_{R}$ turns out to be when one restricts to the inspirals of nearly spherical orbits with $e\approx 0$; we are also working on this topic [74].

Appendix A Linear-in-spin parts of the PN-expanded fluxes

In this Appendix, we present the results for the linear parts of the energy and angular momentum fluxes as series in the PN parameter $v=\sqrt{1/p}$ and eccentricity $e$ with some terms expressed in closed form.

The linear part of the energy flux reads

$$\begin{split}&\frac{\delta\mathcal{F}^{E}}{\mathcal{F}^{E}_{\rm N}\quantity(1-e^{2})^{3/2}}=\delta f_{3}(e)v^{3}+\delta f_{5}(e)v^{5}+\delta f_{6}(e)v^{6}+\delta f_{7}(e)v^{7}+\delta f_{8}(e)v^{8}+\\ &\phantom{\frac{\delta\mathcal{F}^{E}}{\mathcal{F}^{E}_{\rm N}\quantity(1-e^{2})^{3/2}}=}\quantity(\delta f_{9}^{0}(e)+\delta f_{9}^{\log(n)}(e)+\delta f_{9}^{\log(v)}(e)\quantity(\gamma-\frac{35\pi^{2}}{107}+\log(v)))v^{9}+\delta f_{10}v^{10}+\order{v^{11}},\end{split}$$ (91)

where

$$\mathcal{F}^{E}_{\rm N}=\frac{32}{5}\quantity(\frac{\mu}{M})^{2}v^{10}\,,$$ (92)

is the Newtonian flux from circular orbits and $\delta f_{i}(e)$ are functions of eccentricity which take the form

$$\delta f_{3} =-\quantity(\frac{25}{4}+\frac{151e^{2}}{4}+\frac{443e^{4}}{16}+\frac{355e^{6}}{192}),$$ (93)

$$\delta f_{5} =\frac{2403}{112}+\frac{16435e^{2}}{112}+\frac{6701e^{4}}{224}-\frac{78383e^{6}}{896}-\frac{108813e^{8}}{14336}+\frac{329e^{10}}{512}+\order{e^{12}},$$ (94)

$$\frac{\delta f_{6}}{\pi} =-\frac{187}{6}-\frac{32257e^{2}}{96}-\frac{67141e^{4}}{128}-\frac{9184435e^{6}}{55296}-\frac{663581e^{8}}{110592}-\frac{1761277e^{10}}{35389440}+\order{e^{12}},$$ (95)

$$\delta f_{7} =\frac{285211}{4536}+\frac{1027841e^{2}}{1134}+\frac{182196563e^{4}}{72576}+\frac{187711757e^{6}}{145152}-\frac{12275083e^{8}}{129024}+\frac{599257e^{10}}{258048}+\order{e^{12}},$$ (96)

$$\frac{\delta f_{8}}{\pi} =\frac{62471}{672}+\frac{1525357e^{2}}{1344}+\frac{65526409e^{4}}{43008}-\frac{110900285e^{6}}{774144}-\frac{15548354173e^{8}}{49545216}-\frac{1021275631e^{10}}{117964800}+\order{e^{12}},$$ (97)

$$\delta f_{9}^{0} =-\frac{29174232523}{34927200}-\frac{7994937281e^{2}}{436590}-\frac{389472520471e^{4}}{6652800}-\frac{1184172237779e^{6}}{22353408}-\frac{15468420403e^{8}}{698544}$$

$$\hphantom{=}-\frac{351984359281e^{10}}{30105600}+\order{e^{12}},$$ (98)

$$\frac{\delta f_{10}}{\pi} =\frac{86803}{216}+\frac{758931497e^{2}}{96768}+\frac{5808116575e^{4}}{193536}+\frac{2774132820133e^{6}}{83607552}+\frac{6683237111993e^{8}}{668860416}$$

$$\phantom{=}+\frac{524219611964239e^{10}}{1070176665600}+\order{e^{12}}.$$ (99)

The various logarithmic terms then given as

$$\delta f_{9}^{\log(n)} \equiv\delta f_{9}^{\log(2)}\log(2)+\delta f_{9}^{\log(3)}\log(3)+\delta f_{9}^{\log(5)}\log(5)+\delta f_{9}^{\log(7)}\log(7)$$ (100)

$$\delta f_{9}^{\log(v)} =\frac{15943}{105}+\frac{37985e^{2}}{14}+\frac{39804e^{4}}{5}+\frac{172163e^{6}}{28}+\frac{1088725e^{8}}{896}+\frac{257121e^{10}}{8960}$$ (101)

$$\begin{split}\delta f_{9}^{\log(2)}&=\frac{19153}{63}+\frac{243853e^{2}}{630}+\frac{44397296e^{4}}{315}-\frac{3518430389e^{6}}{2268}+\frac{4038595503577e^{8}}{362880}-\frac{1099970344876951e^{10}}{18144000}\\ &\phantom{=}+\order{e^{12}},\end{split}$$ (102)

$$\delta f_{9}^{\log(3)} =\frac{702027e^{2}}{140}-\frac{147191661e^{4}}{2240}+\frac{427378437e^{6}}{1120}+\frac{591272646381e^{8}}{573440}$$

$$\phantom{=}-\frac{1574559880088247e^{10}}{57344000}+\order{e^{12}},$$ (103)

$$\delta f_{9}^{\log(5)} =\frac{31138671875e^{6}}{72576}-\frac{50939732421875e^{8}}{9289728}+\frac{56782099609375e^{10}}{1769472}+\order{e^{12}},$$ (104)

$$\delta f_{9}^{\log(7)} =\frac{56067099797765e^{10}}{5308416}+\order{e^{12}}.$$ (105)

The term $\delta f_{5}$ can be resummed in eccentricity as

$$\delta f_{5} =\frac{1731}{112}+\frac{15399e^{2}}{112}+\frac{11811e^{4}}{224}-\frac{81743e^{6}}{896}-\frac{139613e^{8}}{14336}+6(1-e^{2})^{3/2}\quantity(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4})\,,$$ (106)

which is consistent with the results of Henry and Khalil [19]. We managed to resum also the term $\delta f_{7}$ in the form

$$\delta f_{7} =-\frac{473}{648}+\frac{848905e^{2}}{1296}+\frac{193342649e^{4}}{72576}+\frac{237963833e^{6}}{145152}-\frac{25289371e^{8}}{129024}-\frac{10496681e^{10}}{258048}$$

$$\phantom{=}+(1-e^{2})^{3/2}\quantity(\frac{1781}{28}+\frac{38839e^{2}}{112}+\frac{21935e^{4}}{64}+\frac{15179e^{6}}{448})\,.$$ (107)

The angular momentum fluxes can be expressed as

$$\begin{split}&\frac{\delta\mathcal{F}^{J_{z}}}{\mathcal{F}^{J_{z}}_{\rm N}\quantity(1-e^{2})^{3/2}}=\delta g_{3}(e)v^{3}+\delta g_{5}(e)v^{5}+\delta g_{6}(e)v^{6}+\delta g_{7}(e)v^{7}+\delta g_{8}(e)v^{8}+\\ &\phantom{\frac{\delta\mathcal{F}^{J_{z}}}{\mathcal{F}^{J_{z}}_{\rm N}\quantity(1-e^{2})^{3/2}}=}\quantity(\delta g_{9}^{0}+\delta g_{9}^{\log(n)}(e)+\delta g_{9}^{\log(v)}\quantity(\gamma-\frac{35\pi^{2}}{107}+\log(v)))v^{9}+\delta g_{10}v^{10}+\order{v^{11}},\end{split}$$ (108)

where

$$\mathcal{F}^{J_{z}}_{\rm N}=\frac{32}{5}\frac{\mu^{2}}{M}v^{7},$$ (109)

is the Newtonian flux from circular orbits and the functions $\delta g_{i}(e)$ read

$$\delta g_{3} =-\frac{19}{4}-\frac{683e^{2}}{48}-\frac{403e^{4}}{96},$$ (110)

$$\delta g_{5} =\frac{3559}{224}+\frac{30509e^{2}}{672}-\frac{135161e^{4}}{5376}-\frac{1251e^{6}}{64}+\frac{15e^{8}}{32}+\frac{99e^{10}}{512}+\order{e^{12}},$$ (111)

$$\frac{\delta g_{6}}{\pi} =-\frac{151}{6}-\frac{3671e^{2}}{24}-\frac{21607e^{4}}{192}-\frac{102199e^{6}}{13824}+\frac{88799e^{8}}{442368}-\frac{27779e^{10}}{44236800}+\order{e^{12}},$$ (112)

$$\delta g_{7} =\frac{1006711}{18144}+\frac{2535215e^{2}}{5184}+\frac{1709749e^{4}}{2304}+\frac{1212805e^{6}}{13824}-\frac{795497e^{8}}{18432}+\frac{38149e^{10}}{4096}+\order{e^{12}},$$ (113)

$$\frac{\delta g_{8}}{\pi} =\frac{100369}{1344}+\frac{296213e^{2}}{896}-\frac{7426709e^{4}}{28672}-\frac{620472319e^{6}}{1548288}-\frac{728183081e^{8}}{33030144}+\frac{42902045527e^{10}}{4954521600}+\order{e^{12}},$$ (114)

$$\begin{split}\delta g_{9}^{0}&=-\frac{19353142307}{27941760}-\frac{258040969517e^{2}}{27941760}-\frac{39004110703e^{4}}{2540160}-\frac{7526116434163e^{6}}{1117670400}-\frac{19004463748957e^{8}}{5960908800}\\ &\phantom{=}-\frac{953163710537e^{10}}{425779200}+\order{e^{12}},\end{split}$$ (115)

$$\begin{split}\frac{\delta g_{10}}{\pi}&=\frac{1068677}{3024}+\frac{6760111e^{2}}{1344}+\frac{113252009e^{4}}{9216}+\frac{56647762139e^{6}}{10450944}-\frac{104974004267e^{8}}{95551488}+\frac{174401161771e^{10}}{44590694400}\\ &\phantom{=}+\order{e^{12}}.\end{split}$$ (116)

The various logarithmic terms are then

$$\delta g_{9}^{\log(n)} \equiv\delta g_{9}^{\log(2)}\log(2)+\delta g_{9}^{\log(3)}\log(3)+\delta g_{9}^{\log(5)}\log(5)+\delta g_{9}^{\log(7)}\log(7),$$ (117)

$$\delta g_{9}^{\log(v)} =\frac{2675}{21}+\frac{347429e^{2}}{252}+\frac{1090651e^{4}}{504}+\frac{6934349e^{6}}{10080}+\frac{105823e^{8}}{4480},$$ (118)

$$\begin{split}\delta g_{9}^{\log(2)}&=\frac{80357}{315}-\frac{615143e^{2}}{1260}+\frac{181175503e^{4}}{2520}-\frac{63567501707e^{6}}{90720}+\frac{18959825303e^{8}}{4480}-\frac{3138557026727e^{10}}{162000}\\ &\phantom{=}+\order{e^{12}},\end{split}$$ (119)

$$\delta g_{9}^{\log(3)} =\frac{1794069e^{2}}{560}-\frac{21138813e^{4}}{560}+\frac{196801569e^{6}}{1024}+\frac{136487855331e^{8}}{573440}-\frac{1881487783393587e^{10}}{229376000}+\order{e^{12}},$$ (120)

$$\delta g_{9}^{\log(5)} =\frac{14503515625e^{6}}{82944}-\frac{877191015625e^{8}}{442368}+\frac{503789897265625e^{10}}{49545216}+O\left(e^{12}\right),$$ (121)

$$\delta g_{9}^{\log(7)} =\frac{8289498036460823e^{10}}{2654208000}+\order{e^{12}}.$$ (122)

The second and fourth term in eq. (108) can be again resummed as

$$\delta g_{5} =\frac{2215}{224}+\frac{33029e^{2}}{672}-\frac{104921e^{4}}{5376}-\frac{1401e^{6}}{64}+6(1-e^{2})^{3/2}\quantity(1+\frac{7}{8}e^{2}),$$ (123)

$$\begin{split}\delta g_{7}&=\frac{56743}{18144}+\frac{2111585e^{2}}{5184}+\frac{2036089e^{4}}{2304}+\frac{1977175e^{6}}{13824}-\frac{1549517e^{8}}{18432}\\ &\phantom{=}+(1-e^{2})^{3/2}\quantity(\frac{733}{14}+\frac{35897e^{2}}{224}+\frac{2215e^{4}}{28}).\end{split}$$ (124)

The linear parts of the fluxes as functions of alternate $x-e$ parametrization are given in the Supplemental material [58].

Appendix B Comparison with Phys. Rev. D 108, 104016 (2023)

In this Appendix, we show the derivation of the transformation between time eccentricity $e_{t}$ and Darwin eccentricity $e$ which is needed for the comparison between our results and the results of [19].

In the quasi-Keplerian parametrization and harmonic coordinates $(t_{\text{H}},r_{\text{H}},\phi_{\text{H}})=(t,r-M,\phi)$, the orbit is given as [19]

$$r_{H} =a_{r}(1-e_{r}\cos u)\,,$$ (125)

$$\Omega_{r}t_{H} =u-e_{t}\sin u+f_{v-u}(v-u)+f_{v}\sin v\,,$$ (126)

$$\frac{2\pi}{\Phi}\phi_{H} =v+g_{2v}\sin 2v+g_{3v}\sin 3v\,,$$ (127)

$$\tan\frac{v}{2} =\sqrt{\frac{1+e_{\phi}}{1-e_{\phi}}}\tan\frac{u}{2}\,.$$ (128)

where $a_{r}$ is the semi-major axis, $u$ is the eccentric anomaly, $\Phi$ is the total phase between two successive periastron passages, $v$ is the true anomaly and $f_{v-u}$, $f_{v}$, $g_{2v}$, and $g_{3v}$ are functions given in [19].

First, we find the relation between $(a_{r},e_{r})$ and $(p,e)$ parametrization from the expression for the turning points $r_{1,2}$ and $r^{\text{H}}_{1,2}$

$$\begin{split}&r_{\text{H}1,2}=a_{r}(1\pm e_{r})=r_{1,2}-M\\ &\phantom{r_{\text{H}1,2}}=Mp/(1\mp e)-M\,.\end{split}$$ (129)

The parameters $a_{r}$, $x=(M\Omega_{\phi})^{2/3}$, $e_{r}$, and $e_{t}$ are given in the supplemental material of [19] as functions of $\tilde{E}=-(E-Mc^{2})/\mu$ and $h=L/(GM\mu)$. By inverting the PN series to obtain $\tilde{E}$ and $h$, we were able to express the time eccentricity using the Darwin eccentricity $e$ and the PN parameter $x$ as

$$e_{t}^{2}/e^{2} =\quantity(e_{t}^{2}/e^{2})_{\rm(g)}+\frac{s_{\parallel}}{M}\delta(e_{t}^{2}/e^{2})\,,$$ (130)

$$\begin{split}\delta(e_{t}^{2}/e^{2})&=\frac{2x^{3/2}}{\sqrt{1-e^{2}}}+\frac{6(e^{2}-2+2\sqrt{1-e^{2}})x^{5/2}}{(1-e^{2})^{3/2}}\\ &\phantom{=}+\order{x^{7/2}}\,.\end{split}$$ (131)

The geodesic part $\quantity(e_{t}^{2}/e^{2})_{\rm(g)}$ can be found in Eq. (4.38) of [45].

Alternatively, one can solve the equation for $t$ as a function of the eccentric anomaly $u$ and collect all the terms that generate $e_{t}$, as was done in [49], however, this process is long and difficult and we leave it for future work.

Appendix C Evolution of the orbital parameters

In this Appendix we present the formulas for the evolution of the orbital parameters $p$ and $e$ used in the hybrid model in Eq. (84). The elements of the geodesic part of the inverse Jacobian

read

$$\partial_{E}p =\frac{-2p^{3/2}\sqrt{P_{3}P_{2}}}{P_{1}}\,,$$ (132)

$$\partial_{J_{z}}p =\frac{2(p-4)^{2}\sqrt{P_{3}}}{P_{1}}\,,$$ (133)

$$\partial_{E}e =\frac{(p-6-2e^{2})\sqrt{pP_{2}P_{3}}}{eP_{1}}\,,$$ (134)

$$\partial_{J_{z}}e =-\frac{(1-e^{2})((p-2)(p-6)+4e^{2})\sqrt{P_{3}}}{peP_{1}}\,,$$ (135)

where we introduced the polynomials

$$P_{1} =(p-6)^{2}-4e^{2}\,,$$ (136)

$$P_{2} =(p-2)^{2}-4e^{2}\,,$$ (137)

$$P_{3} =p-3-e^{2}\,.$$ (138)

Note that the polynomial $P_{1}$ vanishes at the separatrix $p=6+2e$, therefore the inverse Jacobian diverges there.

We can factor out some terms from the matrix product $\mathbb{J}_{\rm(g)}^{-1}\delta\mathbb{J}$ and express it in the form

$$\mathbb{J}_{\rm(g)}^{-1}\delta\mathbb{J}=\frac{1}{P_{1}P_{3}\sqrt{P_{2}p^{3}}}\mathbb{M}\,.$$ (139)

where

$$\mathbb{M}_{1,1} =-4e^{6}p+e^{4}(p(p(p+8)-36)+96)-4e^{2}\left(p\left(4p^{2}-30p+83\right)-48\right)-(p-6)(p-2)(p(4p-15)+24)\,,$$ (140)

$$\mathbb{M}_{1,2} =ep\left(-\left(\left(e^{2}-21\right)p^{3}\right)-8\left(e^{2}+25\right)p^{2}+4\left(e^{4}+10e^{2}+165\right)p-768\right),$$ (141)

$$\begin{split}\mathbb{M}_{2,1}&=\frac{1-e^{2}}{4ep}\big{(}16e^{6}(p-3)-4e^{4}((p-4)p(p+14)+60)+e^{2}(p(p(p(5p-12)-96)+336)-144)\\ &\phantom{=}+3(p-6)^{2}(p-2)^{2}\big{)}\,,\end{split}$$ (142)

$$\mathbb{M}_{2,2} =\left(e^{2}-1\right)\left(4e^{4}p-e^{2}\left(p^{3}+16p^{2}-120p+96\right)+2p^{4}-13p^{3}-24p^{2}+228p-288\right).$$ (143)