Eccentric self-forced inspirals into a rotating black hole

We develop the first model for extreme mass-ratio inspirals (EMRIs) into a rotating massive black hole driven by the gravitational self-force. Our model is based on an action angle formulation of the method of osculating geodesics for eccentric, equatorial (i.e., spin-aligned) motion in Kerr spacetime. The forcing terms are provided by an efficient spectral interpolation of the first-order gravitational self-force in the outgoing radiation gauge. We apply a near-identity (averaging) transformation to eliminate all dependence of the orbital phases from the equations of motion, while maintaining all secular effects of the first-order gravitational self-force at post-adiabatic order. This implies that the model can be evolved without having to resolve all $\mathcal{O}(10^5)$ orbit cycles of an EMRI, yielding an inspiral model that can be evaluated in less than a second for any mass-ratio. In the case of a non-rotating central black hole, we compare inspirals evolved using self-force data computed in the Lorenz and radiation gauges. We find that the two gauges generally produce differing inspirals with a deviation of comparable magnitude to the conservative self-force correction. This emphasizes the need for including the (currently unknown) dissipative second order self-force to obtain gauge independent, post-adiabatic waveforms.


Introduction
The detection of the first gravitational wave signal [1] ushered in a new era of astronomy, with ground-based observatories having now observed just over a hundred signals [2][3][4].The next generation of space-based detectors, such as the Laser Interferometer Space Antenna (LISA), will probe the previously inaccessible millihertz band of the gravitational wave spectrum, allowing for the detection of hitherto unseen gravitational wave sources.Among these sources are extreme mass-ratio inspirals (EMRIs), which consist of a stellar mass compact object (such as a black hole or neutron star) spiralling into a supermassive black hole.These systems are characterised by their extremely small mass ratio, typically between 10 −4 and 10 −7 .Unlike the signals detected by groundbased detectors, EMRIs will radiate in the LISA frequency band for up to hundreds of thousands of orbital cycles [5].They are also expected to be eccentric and precessing, potentially resulting in multi-year long waveforms with rich and complex morphology [6].These signals encode the spacetimes of supermassive black holes, promising exquisite parameter estimation and some of the most sensitive probes for new physics beyond general relativity [7].
The majority of EMRIs will have a very low instantaneous signal-to-noise ratio (SNR), and so the data must be processed with matched filtering techniques which will allow for the build up of the SNR over time [8].Such techniques require the development of theoretical waveform templates to compare against the data.To achieve LISA's science objectives, these templates need their phase to be accurate to within a fraction of a radian, even after hundreds of thousands of orbital cycles.They also need to be extensive across the large parameter space of possible EMRI configurations.Moreover, since many template evaluations would be needed, they should also be fast to compute, ideally in less than a second.
To meet this challenge, several so-called "kludge" models have been developed, which are both extensive and quick to compute [9][10][11][12].However, they also make use of non-relativistic approximations which cause them to fall short of the accuracy requirement, though they may still be sufficient for the detection of loud EMRI signals [13].Despite their shortcomings, these models are invaluable for testing data analysis techniques for LISA through the mock data challenges [14][15][16].In order to detect more EMRI signals, relativistic, "adiabatic" waveforms are required.The adiabatic trajectory can be calculated by balancing the fluxes of energy and angular momentum through null infinity the event horizon with the energy and angular momentum lost by the secondary throughout the inspiral.This has been implemented in a practical framework for non-spinning, eccentric EMRIs [17,18].While this represents a significant accuracy improvement over the kludge models, work remains to be done to extend this framework to generic inspirals into rotating black holes.
To detect all EMRI signals and enable precision parameter estimation, requires "post-adiabatic" waveforms.Producing such waveforms requires calculating the local force experienced by the secondary due to the presence of its own gravitational field, known as the gravitational self-force (GSF).This can be calculated perturbatively by expanding the field equations in powers of the small mass ratio of the binary, which makes this approach ideally suited to EMRIs.Post-adiabatic EMRI waveforms will require not only full knowledge of the the first-order self-force, but also the orbit averaged contribution of the second-order self-force [19,20].
At each instant, the self-force is a functional of the the past history of the secondary which can make it challenging to compute.One approach is to couple the field equations and the equations of motion and self-consistently solve both in a time-domain simulation.While this has been implemented for a toy model of a particle carrying a scalar charge orbiting a Schwarzschild black hole [21], numerical stability issues have so far stifled similar attempts for the gravitational case [22].Moreover, this approach is computationally very expensive, making it impractical for generating large numbers of templates.However, it does promise waveforms against which more efficient schemes should be tested.
An alternative method is to compute the self-force for a body moving along fixed geodesics of the background spacetime and then use that force to move to another geodesic at a later timestep.The periodic nature of these geodesic orbits allows for calculations in frequency domain leading to many efficient calculations of the first-order self-force in both Schwarzschild [23][24][25][26] and Kerr [27,28] spacetimes.Second-order calculations are also emerging, though at present these are restricted to quasi-circular inspirals into non-rotating black holes [29,30].These calculations can be repeated across the parameter space of bound geodesics and then interpolated in a preprocessing step, as has been done for eccentric inspirals into a Schwarzschild black hole [31,32].One of the goals of this work is to compute self-forced inspirals into rotating (Kerr) black holes.
With an interpolated model of the geodesic self-force in hand, one can restate the EMRI equations of motion in a more convenient form for numerical integration using the method of "osculating orbital elements" (or "osculating geodesics" when applied to the relativistic context).In this approach, the inspiral is described as a smooth evolution through neighbouring geodesics that are instantaneously tangent to the true inspiral.Formally, one identifies a set of constants of motion which uniquely identify a geodesic, known as "orbital elements".These constants are then promoted to functions of time which are governed by a set coupled ordinary differential equations that are derived from the "osculating conditions".These new equations of motion are then solved numerically to obtain the inspiral trajectory of the secondary.There are a number of equivalent formulations for these ODEs which have been derived for both Schwarzschild [33] and Kerr [34,35] inspirals.In this work, we make use of a formulation based on action angles of the geodesic motion that was sketched out in Ref. [34].
However, for all of these formulations, the resulting equations of motion depend on the orbital phases.This means that the numerical integrator will have to resolve features on the orbital timescale, requiring the use of many small time steps.Since a typical EMRI undergoes ∼ 10 4 −10 6 orbital cycles during a radiation reaction timescale, this results in computational times of minutes to hours for a single inspiral, depending on the orbital configuration.
Following Ref. [36] (hereafter Paper I), we overcome this problem by applying a near-identity (averaging) transformation (NIT) [37] to the self-forced equations of motion.These transformed equations have two important properties: (i) they no longer depend on the orbital phase, and (ii) they capture the long-term secular evolution of the original inspiral to the same order of approximation in the mass ratio as the original equations of motion.The first property means the transformed equation of motion can be numerically solved for any mass ratio in less than a second as the numerical integrator no longer needs to resolve the thousands of oscillations on the orbital timescale.
This approach has been applied to the case of eccentric inspirals into non-rotating black holes [36,38].In this work we apply the NITs to orbital motion in the Kerr spacetime.Combined with an interpolated model of first-order gravitational self-force data, these averaged equations of motion allow us to efficiently compute inspirals around a Kerr black hole which, for the first time, include all first order in the mass ratio effects.Since these inspirals are fast to compute, this approach can provide EMRI waveforms useful for practical data analysis when coupled to a fast waveform generation scheme, e.g., [17].
At this point, we emphasize that the inspirals we present do not reach the subradian accuracy required for EMRI data analysis as there are not yet any second order self-force calculations in this domain.To stress the importance of the second order contribution we examine the effects of driving the inspirals using first-order self-forces computed in two different gauges and demonstrate explicitly that without the inclusion of the second-order self-force the inspiral phase is not gauge invariant.Nonetheless, the framework we present can readily incorporate new self-force results as they become available.
This paper is laid out as follows.In Sec. 2, we review geodesic motion and introduce our action angle formulation of the osculating geodesic equations for generic Kerr orbits.We end this section by specialising to equatorial (i.e., spin aligned) orbits for the rest of this work.Using these equations, along with a model for the gravitational self-force, allows us to calculate eccentric, self-forced, inspirals in Kerr spacetime for the first time.However, these inspirals, henceforth referred to as "osculating geodesic" (OG) inspirals, are slow to compute, taking on the order of hours or days.
This motivates us to apply a near identity transformation, as developed in Paper I, which we summarize in Sec.3.1.In Sec.3.2 we explicitly derive the averaged equations of motion for the case of eccentric Kerr inspirals.Inspirals calculated with these equations of motion can be evaluated in less than a second, and are henceforth referred to as "NIT" inspirals.
For both OG and NIT inspirals, we require a model for the gravitational self-force.To this end, we review the calculation of the gravitational self-force in the radiation gauge in Sec.4.1 before outlining a procedure used to interpolate this data in Sec.4.2.Using this interpolated self force model we describe our numerical implementation for calculating the various terms in the NIT equations of motion in Sec. 5.
We then present the results of this implementation in Sec. 6.We start with a consistency check with energy and angular momentum fluxes and an interpretation of the various terms in the NIT equations of motion in Sec.6.1.We then compare OG and NIT inspirals to check that they agree at the appropriate order in the mass-ratio and assess the speed up the NIT provides.We then explore some post-adiabatic effects of the first-order self force by comparing NIT and adiabatic inspirals in Sec.6.3.To round off our results, in Sec.6.4 we make comparisons between inspiral trajectories around a Schwarzschild black hole calculated using two different first order self-force models: one calculated using outgoing radiation gauge self-force data and the other using Lorenz gauge self-force data.These comparisons indicate that trajectories calculated using only first order self-force data are gauge dependent.We end with some concluding remarks in Sec. 7.
Throughout this paper we work in geometric units such that the gravitational constant and the speed of light are both equal to one (i.e., G = c = 1).

Forced motion near a rotating black hole
In this section we describe the motion of a non-spinning compact object of mass µ moving in the Kerr spacetime under the influence of some arbitrary force.Later in this work, we will take this to be the self-force experienced by the secondary due to its interaction with its own metric perturbation.We denote the mass of the primary by M and parametrize its spin by a = |J|/M where J is the angular momentum of the black hole.The Kerr metric can then be written in modified Boyer-Lindquist coordinates, x α = {t, r, z = cos θ, φ}, as where If a force acts upon the secondary its motion can be described by the forced geodesic equation where u α = dx α /dτ is the secondary's four-velocity, ∇ β is the covariant derivative with respect to the Kerr metric, and a α is the secondary's four-acceleration.We seek to recast Eq. ( 3) into a form useful for applying the near-identity transformations.Before considering the forced equation it is useful to first examine the geodesic limit.

Geodesic motion and orbital parametrization
In the absence of any perturbing force, the secondary will follow a geodesic, i.e., The symmetries of the Kerr spacetime allow for the identification of integrals of motion P = {E, L, K} given by where K αβ is the Killing tensor, E is the orbital energy per unit rest mass µ, L is the z-component of the angular momentum divided by µ and K is the Carter constant divided by µ 2 [39].This definition of the Carter constant is related to another common definition of the Carter Constant, Q, by The geodesic equation can be written explicitly in terms of these integrals of motion [40]: dz dλ where r 1 > r 2 > r 3 > r 4 are the roots of the radial potential V r , z + > z − are the roots of the polar potential V z , and λ is Mino(-Carter) time that decouples the radial and polar equations [41].This time is related to the proper time of the particle, τ , by The two largest roots of V r correspond to the apoapsis and periapsis of the orbit, respectively.Explicit expressions for the other roots are derived in Ref. [42] and given for completeness in Appendix A. Rather than parametrize an orbit by the set {E, L, K} it is useful to instead use more geometric, quasi-Keplerian constants P = {p, e, x}.
Here p is the semi-latus rectum, e is the orbital eccentricity and x measures the orbital inclination.These are related to the radial and polar roots via The explicit relation between the integrals of motion {E, L, K} and {r 1 , r 2 , z − } can be found in, e.g., Appendix B of Ref. [43].We note that one advantage of using x over other other common choices for inclination is that x smoothly parametrizes orbits between prograde equatorial motion with x = 1 to retrograde equatorial motion with x = −1 orbits.Not all values of {p, e, x} correspond to bound geodesics and we denote the value of p at the last stable orbit by p LSO (a, e, x) [44,45].
In order to later apply the near-identity transformations it will be useful to employ action-angle formulation to parametrize the orbital motion [42,43,46].In this description the orbital phases q = {q r , q z } are such that the geodesic equations can written in the form Ṗj = 0 and qi = Υ i ( P ), (10a-b) where an overdot denotes derivative with respect to Mino time and the Υ i are the Mino time fundamental frequencies [42].Note that the right hand side of the qi equation depends only on P and not also on the orbital phases.Semi-analytic solutions to Eqs. (7) in terms of q can be found in Refs.[42,46] and we present the key equations in Appendix A.
At this point we note that it is also common in the literature [47,48] to express the radial and polar motion in terms of quasi-Keplerian angles, ψ and χ, via With this parametrization the evolution equations for ψ and χ depend on both P and q.This makes them inconvenient for deriving averaging transformations, though it is not an insurmountable challenge [35].For this reason we prefer the action-angle formulation.

Osculating Geodesics
We now wish to describe the forced motion of a body obeying Eq. (3).To do so, we make use of the method of osculating orbital elements, or osculating geodesics when applied to the relativistic context [33,34].We first identify a set of orbital elements that uniquely identify a given geodesic, such as the integrals of motion P along with the initial values of the orbital phases of the geodesic orbit q 0 and designate them as a set of "orbital elements" I = { P , q 0 }.For accelerated orbits, these orbital elements are promoted from constants to functions of Mino time.Note that now the orbital elements q 0 (λ) are different quantities from the values of q evaluated at λ = 0, i.e., q(0).We assume that the worldline and four-velocity of the secondary at each instant can be described as the worldline and four-velocity of a test body on a tangent geodesic, i.e., x α (λ) = x α G (I A (λ), λ) and u α (λ) = u α G (I A (λ), λ).From these assumptions, one can derive the osculating geodesic equations [33,34]: From these, evolution equations for each of the orbital elements can be calculated, which we have done in Appendix B and Appendix C.This was first done for generic Kerr orbits in Ref. [34] which lays out four different formulations of the equations.Three of these formulations use the quasi-Keplerian angles ψ and χ as the orbital phases, and two of these were numerically implemented and shown to agree with each other.We make use of the final formulation where one instead uses the geodesic actions angels q r and q z as the orbital phases.
We first find the evolution equations for the integrals of motion P .We do so by finding evolution equations for P and relating these to evolution equations for the roots r 1 , r 2 and z − .We derive these relations in Appendix B. From here, one can invert Eqs. ( 9) to find The orbital phases q still evolve with their respective Mino-time frequencies, but now pick up a correction due to the evolution of the initial values, i.e., To find the evolution equations for the initial values for the orbital phases, we can rearrange the first osculating condition (12a) and exploit the fact that the evolution of r is independent of q z , and the evolution z is independent of q r , to get where x i G is the geodesic expression for r or z given by Eqs.(A.3) and (A.4) respectively.Unfortunately, this expression is difficult to evaluate numerically at the orbital turning points where both ∂x i G /∂q i and the term in the parentheses go to zero.However, one can derive an alternative expression for this quantity that is regular at the turning points [34]: See Appendix C for details on the derivation.This expression instead has a singularity whenever ∂V i /∂x i G = 0. Thus, for our numerical implementation, we use Eq. ( 15) for the majority of the orbital cycle and switch to Eq. ( 16) in the vicinity of turning points.
Finally, we also require evolution equations for "extrinsic quantities" that don't show up on the right hand side of the equations of motion, but are necessary to compute the trajectory and the waveform.These are the time and azimuthal coordinates of the secondary which, as a set, we denote by S = {t, φ}.Their evolution is given by the geodesic equations for t and φ, i.e., equations (7c) and (7d).Putting it all together, the equations of motion take the form: These equations of motion are valid for generic inspirals under the influence of an unspecified perturbing force.We find that the action angle implementation produces inspiral trajectories that are identical to inspirals calculated using the "null tetrad" formulation described in Ref. [34].We have implemented both the action angle and null tetrad osculating element equations into a Mathematica package that will be publicly available on the Black Hole Perturbation Toolkit [49].We find numerically that the null tetrad formulation is more computationally efficient as it does not have any singular equations that necessitate switching between different formulations twice during each orbital period.As such, for direct comparisons between OG and NIT inspirals and waveforms we make use of the null tetrad formulation, but use the action angle formulation as the starting point for our averaging procedure.

Specialising to equatorial motion
For the rest of this work, we specialize to the case of eccentric inspirals in the equatorial plane (i.e., spin aligned) under the influence of the first-order ratio gravitational self force (GSF).This corresponds to setting x = ±1 for prograde and retrograde orbits, respectively.Due to symmetry, motion in the equatorial plane will stay in the equatorial plane, and thus ẋ = 0.As such, we only need to track the evolution of P = {p, e}.
Similarly, the equations of motion no longer depend on the polar phase q z , and so we only need to evolve the radial phase q = {q r }.The gravitational self force scales with the small mass ratio = µ/M , meaning that the secondary's four-acceleration can be expressed as a α = a α GSF .Factoring out this scaling, the equations of motion for equatorial inspirals become dp dλ = F p (a, p, e, q r ), (18a) dq r dλ = Υ r (a, p, e) + f r (a, p, e, q r ), (18c) dt dλ = s t (a, p, e, q r ), (18d) dφ dλ = s φ (a, p, e, q r ).(18e)

Near-Identity Transformations
Near identity (averaging) transformations (NITs) are well known technique in applied mathematics and celestial mechanics [37].This technique involves making small transformations to the equations of motion, such that the short timescale physics is averaged out, while retaining information about the long term evolution of a system.This is well suited to modelling EMRIs, where we do not require perfect knowledge of the trajectory on the orbital timescale, so long as we can accurately track its evolution on the radiation reaction timescale.We note that there can be a relationship between gauge transformations and NITs, as explored in Ref. [50], the NITs we perform in this work are distinct from the choice of gauge [51].In Paper I [36], these transformations were derived for EMRIs in general and then applied to the OG equations of motion around a Schwarzschild black hole.We briefly review the work of Paper I for a general EMRI before applying these transformations to the specific case of eccentric self-forced inspirals in Kerr.

Near identity averaging transformations for generic EMRI systems
The NIT variables, Pj , qi and Sk , are related to the OG variables P j , q i and S k via i ( P , q) + 2 X (2) Here, the transformation functions i , and are required to be smooth, periodic functions of the orbital phases q.At leading order, Eqs.(19) are identity transformations for P k and q i but not for S k due to the presence of a zeroth order transformation term Z (0) k .The inverse transformations can be found for P k and q i by requiring that their composition with the transformations in Eqs.(19) must give the identity transformation.Expanding order by order in , this gives us j ( P, q) ∂ qk X (1) i ( P, q) i ( P, q) ∂ Pj Y (1) i ( P, q) ∂ qk X (1) To proceed it is useful to decompose various functions into Fourier series where we use the convention: where N is the number of orbital phases.Based on this, we can split the function into an averaged piece and an oscillating piece Using the above transformations along with the equations of motion, and working order by order in , we can chose values for the transformation functions Y (1) j , X i , X i , Z k and Z k such that the resulting equations of motion for Pj , qi and Sk take the following form Crucially, these equations of motion are now independent of the orbital phases q.
Deriving the relationship between the transformed forcing functions ( i and s(0\1) k ) to the original forcing functions is quite an involved process with several freedoms and choices, each with their own merits and drawbacks.This is discussed at length in [36], so for brevity we will summarize the results and the particular choices we have made in this work.
The transformed forcing functions are related to the original functions by In deriving these equations of motion, we have constrained the oscillating pieces of the NIT transformation functions to be X( 1) This equation is satisfied by the oscillating pieces for the analytic solutions for the geodesic motion of t and φ in Eqs.(A.2), i.e., We have chosen the averaged pieces such that Y In order to generate waveforms, one only needs to know the transformations in Eq. ( 19) to zeroth order in the mass ratio, i.e., Furthermore, to be able to directly compare between OG and NIT inspirals, we will need to match their initial conditions to sufficient accuracy.To maintain an overall phase difference of O( ) in the course of an inspiral, this requires known the transformation of the P j 's (19a) to linear order in , while it is sufficient to know the rest of Eq. ( 19) to zeroth order.In particular, we will never need an explicit expression for X , which would require solving a PDE to obtain.

Averaged Equations of Motion for Eccentric Equatorial Kerr Inspirals
We now apply the near identity averaging transformation procedure to the equations of motion for equatorial Kerr inspirals to obtain: The leading order terms in each equation of motion are simply the original function averaged over a single geodesic orbit, i.e., F where Υ t and Υ φ are the Mino-time t and φ fundamental frequencies.The remaining terms are more complicated and require Fourier decomposing the original functions and their derivatives with respect to to the orbital elements (p, e).To express the result, we define the operator With this in hand, the remaining terms in the equations of motion are found to be

Gravitational self-force for eccentric Kerr inspirals
In order to drive the inspiral we need to rapidly evaluate the gravitational self-force (GSF) given any values of (p, e, q r ).Codes that compute the GSF generally take minutes to hours to compute the force along a geodesic for a given (p, e) value and so it is not practical to directly couple the equations of motion to a GSF code.Instead, it is common practice to calculate the GSF on a discrete set of points across the parameter space and then build an interpolation or fitted model that smoothly connects the GSF data.The following subsections describe our approach.

Gravitational self-force
The gravitational self-force approach is reviewed extensively elsewhere [52,53] and so we just give a brief overview of the calculations that we employ.The GSF approach starts by expanding the metric of the binary around the metric of the primary, i.e., µν + . . .where ḡµν is the Kerr metric and the h (n) are the n-th order perturbations to the spacetime due to the presence of the secondary.The interaction between these metric perturbations and the motion of the secondary can be derived through a matched asymptotic expansion analysis [53].In this work we use only first-order (in the mass ratio) results, as second order results are still emerging [30].The first order metric perturbation generated by a compact object can be found by solving the linearized Einstein Field Equations with a point particle source moving on a geodesic of ḡµν .The matched asymptotic expansion analysis identifies a regular part of this (divergent) metric perturbation that is responsible for the backreaction on the compact object [53][54][55][56].The result is a (self-)force that appears on the right-hand side of Eq. ( 3) computed from derivatives of the regular metric perturbation [53].
Solving the perturbation equations requires picking a gauge, and the resulting selfforce is gauge dependent [57].The self-force was first computed in the Lorenz gauge where the procedure for obtaining the regular part was best understood.Numerical calculations of the Lorenz gauge self-force have been made in both the frequency- [25,26,58] and time-domains [22][23][24] All these results have been for motion in Schwarzschild spacetime, with one exception in Kerr [59].
Calculating perturbations of the Kerr spacetime is hampered by the lack of separability of the linearized Einstein Field Equations on this background.This difficulty can be circumvented by using the Teukolsky formalism for describing perturbations to the Weyl scalars [60], which is fully separable in the frequency domain.From the Weyl scalars, the metric perturbation can be reconstructed in a radiation gauge [61][62][63].There has also been recent progress understanding how to reconstruct the metric in the Lorenz gauge [64].Regularization of the metric perturbation in radiation gauges is more subtle [65], but self-force calculations in the radiation gauge are now routine [27,28,66,67].
In this work we primarily use the self-force computed using the code of Ref. [27,28,67].This code uses the Mano-Suzuki-Takasugi methods [68][69][70][71] to compute the perturbations to the Weyl scalars in the frequency domain.The metric is then reconstructed into an outgoing radiation gauge (including mass and angular momentum perturbations [72,73] and gauge completion contributions [74]).The metric perturbation is then projected onto a basis of spherical harmonics before regularization is carried out using the mode-sum approach [65,75].Depending on the eccentricity of the orbit the code must compute the metric perturbation by summing over thousands to tens of thousands of Fourier harmonic modes.With the current Mathematica implementation, the self-force for an, e.g., a = 0.9M, p = 3.375, e = 0.5 orbit takes approximately 90 CPU hours to compute.
The self-force can be split into dissipative (time anti-symmetric) and conservative (time-symmetric) contributions [76].
The dissipative pieces causes the orbit to shrink until the secondary plunges into the primary.It also generally causes the orbit to circularize, with the exception being just before the transition to plunge where the orbit gains eccentricity [44,[77][78][79].To produce adiabatic waveforms, we only require knowledge of the orbit averaged dissipative pieces of the first-order self-force.These can be related, via balance laws, to the fluxes of GWs to infinity and down the event horizon.Since calculating fluxes avoids regularization of the metric perturbation, adiabatic inspirals are typically calculated via flux balance laws [44,[78][79][80][81][82][83].The conservative pieces have more subtle effects on the inspiral, such as altering the rate of periapsis advance and the location of the innermost stable circular orbit [50,[84][85][86][87][88][89].
To compute post-adiabatic inspirals requires knowledge of both the dissipative and conservatives pieces of the first-order self-force and the orbit average piece of the secondorder self-force [19].There are as yet no calculations of the latter so we will make do with just the first-order self-force information in this work.This will allow us to explore some of the effects of the conservative self-force on equatorial Kerr inspirals for the first time.

Interpolation method
Driving inspirals requires a model for the self-force that can be rapidly evaluated at each instant during the inspiral.To achieve this we tile the parameter space with GSF data which we can then interpolate.While there have been many different approaches for doing this, we has been implemented for eccentric Schwarzschild inspirals [31,32].In these works the data was interpolated using standard cubic spline methods, which required computing the self-force at tens of thousands of points in the parameter space.While this might not pose much of a problem for the 2D parameter space of eccentric, Schwarzschild inspirals, these approaches would not scale well to the 4D parameter space for generic Kerr inspirals.Motivated by this, as well as the computational expense of the eccentric Kerr self-force code, we build an interpolation model based on Chebyshev polynomials that is accurate to percent level across a 2D slice of the EMRI parameter space using only a few hundred points.
We start by fixing the value of the spin parameter of the primary, which we choose to be a = 0.9M for Kerr inspirals or a = 0 for Schwarzschild inspirals and set the inclination x to be either 1 or −1 for prograde orbits or retrograde orbits respectively.This reduces our parameter space to two parameters; the semilatus rectum p and the eccentricity e.We then define a parameter y using the p and the position of the last stable orbit p LSO .For Kerr orbits, we chose y to be With this parametrization we found that the accuracy of the Chebyshev interpolation is limited by the appearance of cusps at the LSO in the data.To ameliorate their impact we instead used a parameter y given by for later runs in Schwarzschild spacetime.In either case, tiling the parameter space in y instead of p will concentrate more points near the separatrix where the self force varies the most.We let y range from y min = 0 (0.01 for Schwarzschild) to y max = 1 and e range from e min = 0 to e max = 0.5 for Kerr and e min = 0 to e max = 0.3 for Schwarzschild.We define parameters u and v which cover this parameter space as they range from (−1, 1) This parametrization is convenient when using Chebyshev polynomials of the first kind, where the order n polynomial is defined by T n (cos ϕ) := cos(nϕ).The Chebyshev nodes are the roots these polynomials, and the location of the kth root of nth polynomial is given by We then calculate the GSF on a 15 × 7 grid of Chebyshev nodes, with the u values given by the roots of the 15th order polynomial and the v values given by the roots of the 7th order polynomial.At each point on our grid, we Fourier decompose each component of the force with respect to the radial action angle q r .We then multiply the data for each Fourier coefficient by a factor of (1 − y)/(1 − e 2 ), as we find that this smooths the behaviour of the force near the separatrix and improves the accuracy of our interpolation.Next, we use Chebyshev polynomials to interpolate each Fourier coefficient across the (u, v) grid.We then resum the modes to reconstruct our interpolated gravitational self force model: where Using this procedure forces each component to become singular at the last stable orbit.
While the GSF changes rapidly as one approaches the last stable orbit, we do not expect the components of the self force to diverge at the LSO.Understanding the analytic structure of the self-force in this region would likely improve future interpolation models.
We note that the GSF should satisfy the orthogonality condition with the geodesic four-velocity, i.e., a α u α = 0. Interpolation will bring with it a certain amount of error which can cause this condition to be violated.We find empirically that we can reduce this interpolation error by projecting the force so that this condition is always satisfied, i.e., This procedure allows us to create a smooth, continuous model for the gravitational self force with relative errors less than 5 × 10 −3 in the strong field -see Fig. 1.The variation in the accuracy of the model is primarily a by-product of how close a given test point (green cross) is to the data points (white dots) used to create the model.We note that this level of precision would not be sufficient for production grade waveforms for LISA, as we would need the relative error of the orbit averaged dissipative self-force to be less than ∼ 1/ , whereas the oscillatory pieces of the self-force only need to be interpolated to an accuracy of a fraction of a percent [32].Our present interpolation model already likely reaches the latter criteria and a future hybrid method that combines flux and self-force data, similar to the one constructed in Ref. [32], can likely reach the overall accuracy goal.Nonetheless, our present model is more than sufficient to test our averaging procedure and to explore the effects of the GSF for eccentric Kerr inspirals.This will now be treated as the underlying forcing model for both the OG and NIT inspirals.

Implementation
Combining the above model with our action angle formulation of the osculating geodesic equations provides us with everything required to calculate the NIT equations of motion.We first evaluate and interpolate the various terms in the NIT equations of motion across the parameter space.This offline process is costly but it only needs needs to be completed once.By contrast, the online steps are computationally cheap, which allows us to rapidly compute eccentric self-forced inspirals into a Kerr black hole.

Offline Steps
To make the offline calculation we complete the following steps.
(i) We start by selecting a grid to evaluate the NIT functions upon.We chose y values between 0.2 and 0.998 in 320 equally spaced steps and e values from 0.001 to 0.5 in 500 equally spaced steps (160,000 points) in the case of Kerr, or use the same spacing in y but only grid in e from 0.001 to 0.3 in 300 equally spaced steps (96,000 points) in Schwarzschild.‡ ‡ Evaluating the NIT functions is computationally cheap so using a dense grid does not significantly increase the computational burden.Using a dense grid also allows us to use Mathematica's default Hermite polynomial interpolation method for convenience of implementation.The grid spacing is chosen to be sufficiently dense that interpolation error is a negligible source of error for our comparisons between the OG and NIT inspirals, though a less dense grid may also achieve this.(ii) For each point in the parameter space (a, y, e) we evaluate the functions F p\e , f r and s t\φ along with their derivatives with respect to p and e for 30 equally spaced values of q r from 0 to 2π.
(iii) We then perform a fast Fourier transform on the output data to obtain the Fourier coefficients of the forcing functions and their derivatives.
(iv) With these, we then use Eqs.( 32), ( 33), (34) and (35) to construct r and s t\φ for that point in the parameter space.(v) We also use Eqs.(26) and (27) to construct the Fourier coefficients of the first order transformation functions Y p\e and X r .(vi) We then repeat this procedure across the parameter space for each point in our grid.
(vii) Finally we interpolate the values for r and s(1) t\φ along with the coefficients of p\e and X (1) r across this grid using Hermite interpolation and store the interpolants for future use.
We implemented the above algorithm in Mathematica 12.2 and find, parallelized across 20 CPU cores takes, the calculation takes about one day to complete.This is a small price to pay, since these offline steps need only be completed once.

Online Steps
The online steps are required for every inspiral calculation, but are comparatively inexpensive.The online steps for computing an NIT inspiral are as follows.(ii) In order to make comparisons between NIT and OG inspirals, we also load interpolants of the Fourier coefficients of Y p/e and X(1) r and Eq. ( 19) to construct first order near-identity transformations.§ (iii) We state the initial conditions of the inspiral (p 0 , e 0 , q r,0 ) and use the NIT to leading order in the mass ratio to transform these into initial conditions for the NIT equations of motion, i.e., (p 0 , ẽ0 , qr0 ).
(iv) We then evolve the NIT equations of motion using an ODE solver (in this case Mathematica's NDSolve).
As with the offline steps we implement the online steps in Mathematica.Note that steps (ii) and (iii) are only necessary because we want to make direct comparisons between NIT and OG inspirals with the same initial conditions.In general, the difference between the NIT and OG variables will always be O( ), and so performing the NIT transformation or inverse transformation to greater than zeroth order in mass ratio will not be necessary when producing waveforms to post adiabatic order, i.e. with phases accurate to O( ).

Results
In this section we present the results from the NIT equations of motion.We first perform some consistency checks in Sec.6.1.We then show that our NIT and OG inspirals agree to the relevant order in the mass ratio in Sec.6.2.Here we also compute, for the first time, self-forced inspirals in Kerr spacetime.With our fast NIT model we then explore the impact of the conservative effects of the first-order GSF as calculated in radiation gauge for Kerr inspirals in Sec.6.3.Finally, in Sec.6.4, we compare Schwarzschild inspirals calculated using a radiation gauge GSF model and a Lorenz gauge GSF model.

Consistency checks
Before computing inspirals, we perform a series of consistency checks on the NIT equations of motion.A useful feature of the NIT is how it separates adiabatic and post-adiabatic effects of the gravitational self-force.At first order in the mass ratio, this corresponds to the dissipative and conservative pieces respectively.We note that when we substitute a α → a α diss , we find that We compared these to the energy and angular momentum fluxes at infinity tabulated in the Black Hole Perturbation Toolkit [49] and generated with a variant of the Gremlin code [80,81] and found that the balance laws were upheld up to relative errors < 10 −3 throughout the parameter space which is consistent with the interpolation error of our self-force model.From all of this, we can infer the significance of each of the terms in Eq. ( 31): Υ r , Υ t and Υ φ capture the background geodesic motion, F p and F e capture the adiabatic effects due to the first order dissipative self-force, f (1) r , s(1) t , and s(1) φ capture the postadiabatic effects due to the first order conservative self-force, and e capture the interplay between the first order dissipative and conservative self-force, as well as the effect of the orbit averaged contribution from the second order self-force.1, e 0 = 0.48).We show the inspiral computed using the osculating geodesic equations, the NIT equations of motion and the inverse NIT to first order in .The insets zoom into the start and end of the inspiral to reveal the small orbital timescale oscillations.The NIT averages through these oscillations, and when using the inverse NIT to add the oscillations back on, we see that the NIT trajectory remains almost perfectly in phase with the OG trajectory throughout the inspiral.

Comparison between OG and NIT inspirals
In order to test the accuracy of our implementation, we compare inspirals calculated using the OG equations of motion found in Ref. [34] to those calculated using the nearidentity transformed equations of motion.To demonstrate these results, we choose a binary with a primary of mass M = 10 6 M and a secondary of mass µ = 10M for a typical EMRI mass ratio of = 10 −5 .To push our procedure to the limit, we chose the initial conditions of our prograde inspiral to be deep in the strong field and highly eccentric with p 0 = 7.1 and e 0 = 0.48 such that the resulting inspiral would take approximately 1 year to plunge.We also set q r,0 = t 0 = φ 0 = 0 for simplicity.
Figure 2 shows the evolution of p and e over time.The trajectories calculated with the OG equations of motion have order oscillations on the orbital timescale which requires the numerical integrator to take small time steps to accurately resolve.The NIT trajectory does not have these oscillations so the numerical integrator can take much larger steps and still faithfully track the averaged trajectory throughout the entire inspiral.The inverse NIT given in Eq. ( 20) through O( ) can be used to add the oscillations back on to the NIT trajectory.We find that while this is unnecessary for computing accurate waveforms, it demonstrates that the NIT trajectory remains in The difference in the orbital phase and extrinsic quantities for a equatorial Kerr inspiral with = 10 −5 and a = 0.9M calculated using the OG and NIT equations of motion with initial conditions p 0 = 7.1, e 0 = 0.48.We find that the differences remain small throughout the inspiral, only becoming large as the secondary approaches the last stable orbit where the adiabatic approximation breaks down.
phase with the OG trajectory -see the insets of Fig. 2.
The accuracy of our NIT model is further demonstrated by Fig. 3 which shows the absolute difference in the orbital phase q r and the extrinsic quantities t and φ between the NIT and OG evolutions.Over the course of the year long inspiral, |t − ( t − Z Finally, we test the effect the NIT procedure has on the waveform.In principle, we could use our averaged equations of motion in conjunction with the FastEMRIWaveforms (FEW) framework to rapidly compute waveforms with relativistic amplitudes.However, currently, the FEW framework only has amplitude data for Schwarzschild inspirals.As such, we make use of the same procedure as the Numerical Kludge [10] by mapping the Boyer-Lindquist coordinates {t, r, θ, φ} to flat space coordinates and using the quadrupole formula to generate the waveform.The resulting waveforms are only an approximation to the true waveforms, but since both inspiral trajectories are being fed through the same waveform generation scheme this should not bias the results when finding the difference in the waveform as a result of using the NIT trajectory instead of the OG trajectory.
From Fig. 4, we can see that the waveforms generated by each evolution  Two snapshots of the dominant (l, m) = (2, 2) mode of the quadrupole waveform for our prograde, equatorial Kerr inspiral with (a, , p 0 , e 0 ) = (0.9M, 10 −5 , 7.1, 0.48).These snapshots correspond to the first and last hours of the inspiral.This shows that the waveform generated using the NIT trajectory almost perfectly overlaps with the waveform generated using the OG trajectory.It also demonstrates how dramatically an EMRI waveform evolves throughout the inspiral.scheme,sampled every t = 1M ≈ 5s, are almost identical by eye.We can further quantify this by calculating the waveform mismatch using the WaveformMatch function from the SimulationTools [90] Mathematica package and assuming a flat noise curve.From Fig. 5, we see that the mismatch remains below 5 × 10 −8 throughout the inspiral.At this level of mismatch the two waveforms would be completely indistinguishable for EMRIs with signal-to-noise ratio (SNR) of upto (at least) 3000 [91][92][93].
Next, the difference between the OG and NIT quantities should scale linearly with the mass ratio.This is illustrated in Fig. 6, where starting with initial conditions p 0 = 4 and e 0 = 0.2 we evolved the inspiral until it reached p = 3 for mass ratios ranging from 10 −1 to 10 −5 .While working with only machine precision arithmetic we found that for smaller mass ratios the numerical error of the solver of the OG inspiral became dominant The mismatch between the semi-relativistic quadrupole waveforms between inspirals calculated using the OG equations with (a, , p 0 , e 0 ) = (0.9M, 10 −5 , 7.1, 0.48) and the adiabatic EOM matched initial conditions, the adiabatic EOM calculated with matched initial frequencies, and the near-identity transformed EOM.We also mark the mismatch that would be indistinguishable for signals with SNR = 100. .The overlap between OG and NIT waveforms for year-long, prograde, a = 0.9M , equatorial Kerr inspirals as a function of the mass ratio and initial eccentricity.The difference between the two waveforms is less than the accuracy benchmark of 0.97 for mass ratios ≤ 3 × 10 −3 , but not for mass ratios larger than this.While increasing eccentricity does have an effect on the overlap, this effect is not as strong as the effect observed in Fig. 9 of Ref. [38].
over the difference with the NIT.To rectify this, we increased the working precision of our solver to 30 significant digits and found that the difference does, in fact, scale linearly with the mass ratio.This requirement for higher precision only affected the OG solver, the NIT equations of motion can be solved with machine precision arithmetic without introducing any significant error.Since the difference between OG and NIT quantities scales with the mass ratio, it is natural to ask how large can the mass ratio be before the NIT and OG waveforms differ enough to affect data analysis.Following the procedure outlined in Ref. [38], we used our fast NIT inspiral code along with a root-finding algorithm to find the initial value of p that corresponds to a year long inspiral for a given value of the mass ratio and initial eccentricity, and assuming a primary mass of 10 6 M .We use these initial conditions to calculate the overlap between year-long NIT and OG waveforms.This calculation is repeated with mass ratios = {1, 3, 5, 7, 9} × 10 −3 and initial eccentricities e 0 ranging from 0.05 to 0.45 in equally spaced steps of 0.05.The result of this analysis can be seen  1. Computational time required to evolve an inspiral from its initial conditions of p 0 = 7.1 and e 0 = 0.48 to the last stable orbit for different values of the mass ratio, as calculated in Mathematica 12.2 on an Intel Core i7 @ 2.2GHz.The computational time for the OG inspiral scales inversely with the mass ratio, whereas the computational time for NIT inspirals is independent of the mass ratio.This demonstrates how the smaller the mass ratio of the inspiral, the greater speed-up one obtains from using the NIT equations of motion.
in Fig. 7.This demonstrates that NIT and OG waveforms have overlaps larger than the benchmark of 0.97 [94] for mass ratios less than ≈ 3 × 10 −3 , but these overlaps decrease substantially for mass ratios larger than this.We also see that the overlap generally decreases as the initial eccentricity increases, though this effect is not as strong as the effect demonstrated by a similar analysis in Ref. [38] for NITs applied to highly eccentric inspirals in Schwarzschild.They also found that the mismatch between NIT and OG waveforms became substantial for mass ratios larger than 2 × 10 −4 .These differences between the two analyses are most likely the result of our inspirals being deeper in the strong field and driven by a self-force computed in a different gauge (Ref.[38] uses the Lorenz gauge self-force).Such mismatches should not be an issue for EMRI data analysis as EMRIs have mass ratios that range from 10 −7 to 10 −4 .However, these mismatches become significant for intermediate mass ratio inspirals, with mass ratios between 10 −4 to 10 −1 .Since both the OG and NIT equations of motion are formally valid to the same order in the mass ratio, it is not clear a priori which of the two would be closer to the true inspiral.When completed at 1-post-adiabtatic (1PA) order the two sets of equations represent different resummations of the 1PA equations of motion, differing only in their higher order (2+) PA terms.The fact that we are seeing a significant difference between these two resummations for intermediate mass ratios suggests that such higher order PA terms might become relevant.However, in this case it might just be signalling the importance of the missing orbit-averaged dissipative self-force term at 1PA order.Finally, we note that using the NIT equations of motion produces a substantial speed-up over using the OG equations.From Table 1, we see the typical computation time for an inspiral starting at p 0 = 7.1 and e 0 = 0.48 and evolved until the inspiral reaches the last stable orbit for different values of the mass ratio.We see that as we decrease the mass ratio by an order of magnitude, the OG inspiral takes roughly an order of magnitude longer to compute, as it would have to resolve an order of magnitude more orbital cycles before reaching last stable orbit.The NIT inspirals all take roughly the same amount of time to evolve to the last stable orbit, regardless of the mass ratio.
Using our current Mathematica implementation, the NIT inspirals can be computed in less than a second.This time could be further reduced tens of milliseconds if one uses a compiled language such as C/C++, as was done in Paper I [36].We see that using the NIT equations of motion is most advantageous for long inspirals with small mass ratios.Another benefit of using the NIT is that the inspiral requires taking fewer time steps, which results in less numerical error, making it easier achieve a given target accuracy.
The only disadvantage of our formulation is that our final trajectory is parametrized in terms Mino time λ, whereas LISA data analysis applications will need waveforms parametrized Boyer-Lindquist retarded time t.Since our formulation also outputs t(λ), we can numerically invert this to get λ(t) which allows us to resolve this issue at the cost of additional computation time.This was also a problem with the NIT formulation in Schwarzschild where the final trajectory is outputted as a function of the quasi-Keplerian angle χ [36,38].This problem might be circumvented entirely by performing an additional transformation to our NIT equations of motion which would produce averaged equations of motion parametrized by t as outlined in [35].
Since we are now satisfied that our formulation can produce fast and accurate selfforce driven trajectories, we can now use this procedure to explore the phenomenology of eccentric, equatorial Kerr inspirals.

Impact of adiabatic and post-adiabatic effects
With the ability to generate fast and accurate inspirals, we can survey the physics of equatorial Kerr inspirals and examine how this differs from the Schwarzschild case.From Fig. 8a, we see the familiar effect of gravitational radiation reaction on the semilatus rectum, p, and eccentricity, ẽ, whereby p and ẽ both decrease over the inspiral with ẽ growing a little as the last stable orbit is approached [44,78,79].As the inspiral approaches the last stable orbit adiabaticity breaks down and the inspiral undergoes a transition to plunge [95][96][97][98].As such, we stop our inspirals just before the last stable orbit.Our results are the first inspirals to include conservative self-force corrections to the equations of motion in Kerr spacetime.The initial phase q r,0 only evolves secularly when conservative self-force corrections are present and so we use this as a measure of the influence of these corrections [31].This is illustrated by the dashed orange curves in Fig. 8a, which mark the number of radians qr,0 will evolve from a given pair of initial conditions (p 0 , ẽ0 ) until the last stable orbit.For retrograde Kerr (and Schwarzschild orbits in Fig. 10), we find that qr,0 increases throughout the inspiral, whereas for prograde Kerr qr,0 decreases during the inspiral before increasing slightly just before plunge.This is consistent with the change of sign in the correction to the rate of periapsis advance induced by the conservative self force as a function of spin in the circular orbit limit [88] -see Appendix D for further details.
As discussed in Sec.6.1, one can readily calculate adiabatic inspirals using the NIT equations of motion by simply neglecting the post-adiabatic terms.However, when trying to determine how post-adiabatic corrections effect the inspiral, one must be   (p, ẽ) space for prograde and retrograde equatorial Kerr inspirals with = 10 −5 and a = 0.9M .From these plots, we see the familiar behaviour of EMRIs losing eccentricity as the compact object approaches the primary and then gaining eccentricity just before crossing the separatrix (dashed black line).The dashed orange curves are contours that mark the number of radians qr,0 will evolve from a given point until plunge.The conservative self-force for retrograde orbits has a similar effect to the non-spinning case as it causes qr,0 to increase throughout the inspiral.In the prograde case, qr,0 decreases for most of the inspiral and then slightly increases shortly before plunge.mindful of how one matches up an adiabatic inspiral with its post-adiabatic counterpart.Following the argument found in Refs.[31] and [32], matching the initial conditions (p 0 , ẽ0 ) results in an error in the orbital phases that grows linearly in t as the conservative self-force changes the orbital frequencies [86].Instead, one should instead match the Boyer-Lindquist time fundamental frequencies Ω r and Ω φ .For an adiabatic inspiral, these are directly related to the Mino-time fundamental frequencies via Ω Ad r\φ = Υ r\φ Υt [42].To calculate these frequencies as perturbed by the conservative self-force, one can either follow the method outlined in Ref. [32], or one can calculate them directly from the NIT equations of motion: We find that both approaches give the same result up to an error that scales as 2 .With this in hand, we can now choose a value for our initial conditions (p SF 0 , ẽSF 0 ) for our self-forced inspiral, and then root find for initial conditions (p Ad 0 , ẽAd 0 ) that satisfy the simultaneous equations Difference in φ as a function of t between an adiabatic and a first order self-forced inspiral when either matching initial conditions or matching the initial Boyer-Lindquist frequencies.The self-forced inspiral has initial conditions (p 0 , ẽ0 ) = (7.1,0.48) with mass ratio = 10 −5 .Matching initial conditions results in an error that grows linearly with t, while matching frequencies produces an error that is initially constant and then grows quadratically with t.
Using this procedure to match the initial frequencies we find that the linear-in-t growth of the difference in the orbital phases is removed and the phase difference grows quadratically in t as expected -see Fig. 9.

Comparing inspirals driven using radiation Gauge and Lorenz gauge self-force in Schwarzschild spacetime
We now turn our attention to the special case of Schwarzschild (a = 0), where we now have interpolated GSF models calculated in two different gauges.In addition to our outgoing radiation gauge self-force model, we make use of an interpolated Lorenz gauge self-force from Ref. [31], which is valid in the domain 6 ≤ p ≤ 12 and 0 ≤ e ≤ 0.2.We apply the same NIT procedure to inspirals driven by this force model, and find agreement with inspirals calculated in Paper I, up to the precision of the numerical solver.
To assess the accuracy of the dissipative self-force, we calculate the orbit averaged energy and angular momentum fluxes, and find that they agree with values from the literature with a relative error less than 10 −3 for both models across the parameter space.To assess the accuracy of the conservative self-force, we calculate the periapsis  .Sample Schwarzschild trajectories through (p, ẽ) space using either a radiation gauge or a Lorenz gauge model, accompanied by contours denoting the change in qr,0 (in radians) by the end of the inspiral if the inspiral had started in that point of the parameter space.While there are only slight differences in the (p, ẽ) trajectories, there is a stark difference in the evolution of qr,0 induced by each model.advance in the circular orbit limit as outlined in [85] using the formula found in [88].We find that both models show good agreement with the literature across the Lorenz gauge model's domain of validity, with the Lorenz gauge model producing errors less than 10 −3 and the radiation gauge model producing relative errors less than 10 −2 .
While we find good agreement between the two results for gauge invariant quantities, we see from Fig. 10 that the inspirals experience dramatically different conservative effects, depending on the gauge used.While in both cases, the conservative self-force acts against geodesic periapsis advance, we see that the evolution of q r,0 depends heavily on the gauge involved, while the trajectories through p and e space are less affected.This is to be expected as the leading order averaged rates of change of p and e are related to the gauge invariant asymptotic fluxes, while the change in q r,0 is induced entirely by the (gauge dependent) conservative self-force [27].
Just as when comparing adiabatic and self-forced inspirals it is important to match the initial frequencies (rather than the initial (p, e) values).We note that for the Lorenz gauge model, we must account for the fact that the perturbed time coordinate, t, is not asymptotically flat [99].We can define an asymptotically flat time coordinate for Lorenz gauge inspirals via the following rescaling where α is given by We make use of a code provided to us by S. Akcay to numerically calculate this quantity for Lorenz gauge values of p and e [25,100].Equation ( 46) means the perturbed Boyer-Lindquist frequencies must also be rescaled by: r(LG) and Ω φ(LG) In the radiation gauge model, the corresponding subtleties have been dealt with by including the gauge completion corrections, so the frequencies can be calculated using Eq. ( 44) as before.Thus, we can choose a value for p(LG) This allows us to make comparisons between inspirals driven by self-force models calculated in different gauges.We use an inspiral driven by the Lorenz-gauge force model with initial conditions (p 0 , e 0 ) = (11, 0.18), mass ratio = 10 −5 as our reference inspiral which should last just over two and a half years for a 10 6 M primary.In Fig. 11, we see the difference in the phase of the waveform Φ as a function of time between the Lorenz gauge NIT inspiral, and a number of reference models.We make use of the relations between the NIT quantities and the waveform phases derived in Ref. [38] to find We then feed the solutions for {p(t), ẽ(t), Φ r (t), Φ φ (t)} into the FastEMRIWaveforms package to generate these eccentric Schwarzschild waveforms [17].Finally, we make use of the SimulationTools Mathematica package to calculate the mismatches and decompose the waveforms into a single evolving amplitude A(t) and phase Φ(t).This allows us to find the difference in the waveform phase ∆Φ(t) between the Lorenz gauge inspiral and the other inspiral calculations.We use this as our point of comparison as the waveform phase is an observable and thus a gauge invariant quantity.
We note that in each case, we see constant error which gives way to quadratic growth with t just as in Fig. 9.As we discussed in Sec.6.3, this shows that the initial frequencies were correctly matched.From the blue curve, we see that the NIT radiation gauge inspiral quickly goes out of phase with the Lorenz gauge NIT inspiral, resulting in a very large mismatch of 0.93.We found that the largest source of error here is due to interpolation error for in the adiabatic pieces of the NIT.Since these are related to the gauge invariant fluxes, these pieces should be identical in both models.As such, we can rectify this error by using the Lorenz gauge functions for the adiabatic pieces and Φ = 0, and sampled every ∆t = 1M ∼ 5s.We also show the mismatch (MM) between the waveforms in each case.By matching the initial frequencies, we compare an inspiral calculated using a radiation gauge self-force model, an adiabatic inspiral, an inspiral with the adiabatic pieces of the Lorenz gauge model and conservative pieces from the radiation gauge model, and a Lorenz gauge model with a 10 % relative error added to each conservative piece.In all cases the difference grows quadratically in time.This plot suggests that post-adiabatic waveforms calculated using only the first-order self-force differ significantly depending on the gauge used.
and continue to use the radiation gauge functions for the conservative pieces of the NIT equations of motion.The improvement is evident in the green curve, which shows much better agreement with the Lorenz gauge NIT inspiral, with the mismatch falling to 0.83.However, it is only slightly better than matching an adiabatic inspiral (orange curve) using Eq. ( 45) resulting in a mismatch of 0.86.Both radiation gauge and adiabatic inspirals go out of phase by almost 100 radians by the time they reach the last stable orbit.
In order to rule out the possibility of interpolation error of the conservative effects being the primary cause of this difference, we repeat the Lorenz gauge inspiral, but this time we manually add a relative error of δ = 0.1 to all of the conservative pieces of both the NIT equation of motion and our matching procedure for the initial conditions, e.g., qr = Υ r + f (1) r etc.We note that this is an order of magnitude larger than the 10 −2 error produced by the radiation gauge model when calculating the gauge invariant quasi-circular periapsis advance.From the red curve, we see that manually adding a constant 10% relative error results in phase difference and a mismatch (0.54) which is significantly smaller than what we observe between the two self-forced inspirals.This gives us confidence that this difference is not dominated by numerical error.From these investigations, we infer that the trajectories driven using only the first order self-force are gauge dependent, and thus, so too are their waveforms.Since post-adiabatic waveforms are an observable quantity, this leads us to conclude that incorporating the orbit-averaged dissipative second-order self-force will be necessary to obtain gauge invariant, post-adiabatic waveforms.Moreover, since the difference between the Radiation and Lorenz gauge self-forced inspirals is of the same magnitude as the difference with the adiabatic inspiral, we further conclude that the impact of the orbit-averaged dissipative second-order self-force must be of a similar magnitude in at least one of the two gauges.

Discussion
In this paper, we present the first self-forced inspirals in Kerr spacetime.We computed the self-force in the radiation gauge using the code of Ref. [27] and interpolated it over a region of the parameter space of eccentric, equatorial orbits using Chebyshev interpolation.Our model achieves sub-percent accuracy for the self-force across the two dimensional parameter space using only 105 points is a substantial improvement over cubic spline interpolation which would require O( 103 ) points to achieve a comparable level of accuracy.So far we have applied our method to strong-field regions of the parameter space for three values of the primary's spin (a = 0, ±0.9M ).It remains as future work to interpolate over the spin of the primary, however, the Chebyshev interpolation method appears to be a promising approach to tiling data from expensive gravitational self-force codes across the 4-dimensional generic Kerr parameter space.This method could be further improved with the aid of a detailed of the study of the analytic structure of the GSF near the last stable orbit.
With an interpolated self-force model in hand, we computed inspirals using an action-angle formulation of the method of osculating geodesics (OG).This approach is sketched in Ref. [34] and in we implement these equations of motion for generic (eccentric and inclined) inspirals about a Kerr black hole.Our Mathematica implementation will be made publicly available on the Black Hole Perturbation Toolkit [49].For a binary with (small) mass ratio, , numerically solving the osculating geodesic equations takes minutes to hours due to the need to resolve the ∼ 1/ oscillations in the orbital elements.To overcome this, we follow Paper I [36] and apply near-identity (averaging) transformations (NIT) which produce equations of motion that capture the correct longterm secular evolution of the binary but can also be rapidly numerically solved.
As a test of this formulation, we applied it to our eccentric, equatorial self-forced inspirals.We showed that our NIT'd quantities remain close to the original evolution variables throughout the inspiral at the expected order in the mass-ratio.When the mass ratio is greater than 1 : 300, we find the difference between year-long NIT and OG inspirals becomes significant for data analysis, reinforcing the findings of Ref. [38].Note, however, that a priori it is not known which (the NIT or OG inspiral) is closer to the true inspiral, since both are accurate to the same order in the mass-ratio.
With our efficient NIT model of eccentric, equatorial inspirals we explored the effects of the gravitational self force.We find that prograde inspirals around a rapidly rotating black hole generally experience an additional periastron advance on top of the periastron advance induced by geodesic motion.This is in contrast to the "periastron retreat" experienced by retrograde inspirals and inspirals around non-rotating black holes [101].
The NIT equations of motion make it convenient to compare inspirals both with and without post-adiabatic effects included and we confirmed that without post-adiabatic effects, the orbital phases of a typical EMRI will incur an error of order O( 0 ).Moreover, by comparing inspirals under the influence of self-force models calculated in different gauges, we find that the resulting trajectories are gauge dependent.This difference due to gauge causes a de-phasing that is comparable in magnitude to not including any post-adiabatic effects.This suggests that in order to obtain gauge invariant postadiabatic waveforms, one must also include second order self-force results.Second order self-force calculations are presently made using a two-timescale framework [30,102].For the equations of motion, this framework is related to the NITs [35,37], but more work is required in order to explicitly transform the two-timescale results into forcing terms that could be used in the framework presented in this paper.Many of the averaging techniques developed in this work are also useful for the two-timescale approach [35].
For complete post-adiabatic waveforms, one would also need to include the spin of the secondary.Inspirals incorporating the leading conservative spin induced effects around a non-rotating primary have been calculated [103] and the effect of (anti-)aligned secondary spin on the energy and angular momentum fluxes have recently been computed for eccentric orbits [104].These effects can readily be incorporated into the NIT framework.
Another natural extension of our work is to non-equatorial orbital motion.We already have results for spherical inspirals that will be the topic of an upcoming paper.After that we plan to tackle generic orbits, but there are two major barriers to this.The first is the larger parameter space which will make calculating the self-force extremely expensive [83].Our Chebyshev interpolation method should help to reduce the number of points in the parameter space where the self-force needs to be calculated.The second barrier is the presence of orbital resonances [105][106][107][108][109].
Near these resonances the NITs break down and an alternative averaging procedure over the resonance timescale needs to be applied [35].While a lot is known about the effects of resonances on EMRI trajectories [105][106][107][108][109][110], rapidly computing inspiral trajectories while incorporating all resonant effects remains an open challenge.
Finally, we note that in this paper we use the leading-order quadrupole formula to generate the waveforms from the OG and NIT inspirals.This is sufficient for our purposes where we wish to compare waveforms from OG and NIT inspirals but for LISA data analysis we will want to use relativistic waveform amplitudes.These were recently efficiently interpolated in Ref. [17] for Schwarzschild inspirals.That work used orbit-averaged fluxes to drive the inspirals but it would be straightforward to use a NIT inspiral instead.Once the waveform amplitudes have been interpolated for Kerr inspirals it could be combined immediately with the implementation presented in this work.in Ref. [42] and then presented in a simplified form in Ref. [46].The radial and polar solutions to the geodesic equations are given by where sn is Jacobi elliptic sine function, K is complete elliptic integral of the first kind, and

Evolution Equations for the Integrals of Motion
Our goal for this section is to derive evolution equations for the integrals of motion P = {p, e, x}.To do so, we must first consider how a different set of integrals of motion P = {E, L, K} evolve in terms of the covariant components of the particle's 4-acceleration {a t , a r , a z , a φ }.
Using the second osculating geodesic equation (12b) along with definitions of E and L in Eq. ( 5), we can obtain the evolution equations for E and L: To find the the evolution of K, we note that the contravariant components of the Killing tensor can be written as [34] where l and n are null vectors with components Taking the derivative of K from Eq. (5c) with respect to proper time gives us Expanding this out explicitly while making use of the orthogonality condition, g αβ u α a β = 0 and converting to Mino time gives us: Using the above equations, we can express how the roots {r 1 , r 2 , z − } evolve with Mino time by exploiting the same trick as in Appendix A.2 and A.3 of [34].First, we note that, using the chain rule, we can express the rate of change of r 1 or r 2 as We then find expressions for ∂r 1,2 /∂P j be differentiating V r (r) with respect to P j .
∂V r ∂P j r=r We then note that the coefficients proceeding ∂r 1,2 /∂P j are also obtained by differentiating V r with respect to r and then evaluating at r 1,2 , i.e., where we have defined F (r) := (r) We can use similar steps as above to find the evolution of z − .Again, the chain rule tells us that the evolution of z − follows We can then take the derivative of these equations with respect to λ and use the chain rule to obtain Eqs.(13).

Appendix C. Evolution of the Orbital Phases
While it is straightforward to obtain Eq. ( 15) from the first osculating geodesic equation (12a) this equation is difficult to evaluate numerically at the turning points of the orbit, i.e., when ∂x i G /∂q i = 0.As such, following the procedure described in Ref. [34], we derive an equivalent evolution equation for the initial phases which is finite at the turning points.
We start by considering the definition of the geodesic potentials: If we add together the derivative of both sides of this expression with respect to P j and then multiply both sides by Ṗj , one obtains: We note that the square bracket term will vanish for geodesics.For perturbed orbits this means that this term will be proportional to the component of the four-acceleration a i scaled by a factor of Σ 2 to compensate for taking derivatives with respect to λ instead of τ .When evaluating this expression, we make use of the osculating condition x i (λ) = x i G (λ).This leads is to the simplification Combining these results with equation (C.5) gives us our final expression for the evolution of the initial phases which is regular at the turning points, as expressed in Eq. 16.
Appendix D. Self-force corrections to the periapsis advance around a spinning black hole The periapsis advance is a observable quantity that has been used to compare models of compact binary dynamics [101].The effect of the gravitational self-force on this observable for quasi-circular EMRIs around a rotating primary was explored in Ref. [88].
One important insight, that is present in the supplemental material of that work, is the effect that the spin of the primary and orbital radius has on the self-force correction to the rate of periapsis advance.For completeness, we highlight this result in this appendix.For quasi-circular inspirals the relation between the dimensionless quantity W = Ω 2 r /Ω 2 φ and Ω φ is an important benchmark for comparing between different calculational approaches to the two-body problem.[88,101,111].The linear in mass ratio correction to the quantity is defined via W ( ; a, Ω φ ) = W (0; a, Ω φ ) + ρ(a, Ω φ ) + O( 2), (D.1) where W (0; a, Ω φ ) is the background value for the periapsis advance, and ρ(a, Ω φ ) is the correction induced by the first-order gravitational self-force.Fig. D1 demonstrates how ρ varies as a function of orbital radius r and the spin of the primary a.We plot the ratio r ISCO /r, where r ISCO is the radius of the innermost stable circular orbit (ISCO).This ratio is convenient for plotting the results as goes from 1 at the ISCO for all spin values and asymptotically approaches zero as r grows large.As one would expect, the plot demonstrates that this correction grows larger  .The linear in mass ratio correction, ρ, to the periapsis advance, W , as a function of distance from the innermost stable circular orbit (ISCO) r ISCO /r and the spin of the primary, a.The contours show that ρ grows larger as the radius of the inspiral approaches the ISCO.ρ is positive for most of the parameter space, including for all retrograde orbits, implying that in these regions the self-force acts against the geodesic periapsis advance.However, for all prograde orbits, there is a region (in blue) where this correction is negative meaning the self-force instead increases the rate of periapsis advance.The region where this occurs grows larger as a increases The black crosses mark the location of the underlying data from Ref. [88] used to calculate the contour plot.
as the radius of the inspiral approaches the ISCO.This correction is positive for all retrograde orbits and in the strong field for prograde orbits.This means that self-force typically acts against the periapsis advance caused by the background geodesic motion, resulting in a reduction of the observed periapsis advance of the binary.However, for positive spins and at large radii, there is a region of the parameter space (in blue) where this correction is negative, meaning that the self-force increases the observed rate of periapsis advance compared to the background geodesic motion.The larger the spin, the smaller the radii at which this effect occurs.As such, this effect is most prominent

Figure 1 .
Figure1.The relative error of the components of the interpolated gravitational self force model for prograde equatorial orbits with a = 0.9M .The white dots represent the data points that were interpolated.The green crosses represent the data set that the model was tested against.The black dashed line represents the location of the last stable orbit.The relative error was calculated using the normalised L 2 error over a single orbital cycle.
t\φ , define the NIT equations of motion.

Figure 2 .
Figure2.The trajectory through (p, e) space for an inspiral with = 10 −5 , a = 0.9M , and initial conditions (p 0 = 7.1, e 0 = 0.48).We show the inspiral computed using the osculating geodesic equations, the NIT equations of motion and the inverse NIT to first order in .The insets zoom into the start and end of the inspiral to reveal the small orbital timescale oscillations.The NIT averages through these oscillations, and when using the inverse NIT to add the oscillations back on, we see that the NIT trajectory remains almost perfectly in phase with the OG trajectory throughout the inspiral.

Figure 3 .
Figure 3.The difference in the orbital phase and extrinsic quantities for a equatorial Kerr inspiral with = 10 −5 and a = 0.9M calculated using the OG and NIT equations of motion with initial conditions p 0 = 7.1, e 0 = 0.48.We find that the differences remain small throughout the inspiral, only becoming large as the secondary approaches the last stable orbit where the adiabatic approximation breaks down.
Last hour of the waveform.

Figure 4 .
Figure 4.Two snapshots of the dominant (l, m) = (2, 2) mode of the quadrupole waveform for our prograde, equatorial Kerr inspiral with (a, , p 0 , e 0 ) = (0.9M, 10 −5 , 7.1, 0.48).These snapshots correspond to the first and last hours of the inspiral.This shows that the waveform generated using the NIT trajectory almost perfectly overlaps with the waveform generated using the OG trajectory.It also demonstrates how dramatically an EMRI waveform evolves throughout the inspiral.

Figure 5 .
Figure 5.The mismatch between the semi-relativistic quadrupole waveforms between inspirals calculated using the OG equations with (a, , p 0 , e 0 ) = (0.9M, 10 −5 , 7.1, 0.48) and the adiabatic EOM matched initial conditions, the adiabatic EOM calculated with matched initial frequencies, and the near-identity transformed EOM.We also mark the mismatch that would be indistinguishable for signals with SNR = 100.

Figure 6 .Figure 7
Figure 6.The absolute difference in the quantities of a prograde inspiral with a = 0.9M and e 0 = 0.2 after evolving from p = 4 to p = 3 using either the OG or NIT equations of motion.We observe that all the differences follow the solid, black reference curve, as expected.

Figure 8 .
Figure 8. Sample trajectories through (p, ẽ) space for prograde and retrograde equatorial Kerr inspirals with = 10 −5 and a = 0.9M .From these plots, we see the familiar behaviour of EMRIs losing eccentricity as the compact object approaches the primary and then gaining eccentricity just before crossing the separatrix (dashed black line).The dashed orange curves are contours that mark the number of radians qr,0 will evolve from a given point until plunge.The conservative self-force for retrograde orbits has a similar effect to the non-spinning case as it causes qr,0 to increase throughout the inspiral.In the prograde case, qr,0 decreases for most of the inspiral and then slightly increases shortly before plunge.

Figure 9 .
Figure 9. Difference in φ as a function of t between an adiabatic and a first order self-forced inspiral when either matching initial conditions or matching the initial Boyer-Lindquist frequencies.The self-forced inspiral has initial conditions (p 0 , ẽ0 ) = (7.1,0.48) with mass ratio = 10 −5 .Matching initial conditions results in an error that grows linearly with t, while matching frequencies produces an error that is initially constant and then grows quadratically with t.

Figure 10
Figure 10.Sample Schwarzschild trajectories through (p, ẽ) space using either a radiation gauge or a Lorenz gauge model, accompanied by contours denoting the change in qr,0 (in radians) by the end of the inspiral if the inspiral had started in that point of the parameter space.While there are only slight differences in the (p, ẽ) trajectories, there is a stark difference in the evolution of qr,0 induced by each model.

Figure 11 .
Figure 11.The difference in the waveform phase Φ for various inspirals as a function of t when compared to NIT inspiral driven by a Lorenz gauge self-force model, with initial conditions (a, p 0 , e 0 ) = (0, 11, 0.18), mass ratio = 10 −5 , viewing angles Θ = π/4 and Φ = 0, and sampled every ∆t = 1M ∼ 5s.We also show the mismatch (MM) between the waveforms in each case.By matching the initial frequencies, we compare an inspiral calculated using a radiation gauge self-force model, an adiabatic inspiral, an inspiral with the adiabatic pieces of the Lorenz gauge model and conservative pieces from the radiation gauge model, and a Lorenz gauge model with a 10 % relative error added to each conservative piece.In all cases the difference grows quadratically in time.This plot suggests that post-adiabatic waveforms calculated using only the first-order self-force differ significantly depending on the gauge used.

2 − − z 2 15 )= 2r 1 r 2 M (r 1 + r 2 ) , e = r 1 − r 2 r 1 + r 2 ,
∂Vz ∂P j z − along with the second definition V z in equation (7b) to find an expression for ∂z − ∂P j ∂V z ∂P j z − = −2z − (βz the first definition of V z in terms of {E, L, Q} gives us the following explicit expressions for ∂Vz∂P j z − ∂V z ∂E z − = 2a 2 Ez 2 − (1 − z 2 − ), ∂V z ∂L z − = −2Lz 2 − ,and∂V z ∂Q z − = 1 − z 2 − .(B.14a-c)Combining the results from equations (B.12), (B.13) and (B.14) gives us 2z − (z + − βz − ) (1 − z 2 − )Since we have expressions for the evolution of {E, L, K}, we can derive and expression for the evolution of Q by taking the derivative of equation (6) with respect to λ: results and simplifying yields our final expression for the evolution of z −dz − dλ = 1 2z − (z 2 + − βz 2 − ) x 2 dK dλ − 2 L − ax 2 E dL dλ − ax 2 dE dλ , (B.17)where we've used Eq.(9c) to tidy up the final expression.Now that we know how {r 1 , r 2 , z − } evolve, determining the evolution of {p, e, x} is straightforward since we can convert from one set to the other using the relations p and x = ± 1 − z 2 − .(B.18a-c) Figure D1.The linear in mass ratio correction, ρ, to the periapsis advance, W , as a function of distance from the innermost stable circular orbit (ISCO) r ISCO /r and the spin of the primary, a.The contours show that ρ grows larger as the radius of the inspiral approaches the ISCO.ρ is positive for most of the parameter space, including for all retrograde orbits, implying that in these regions the self-force acts against the geodesic periapsis advance.However, for all prograde orbits, there is a region (in blue) where this correction is negative meaning the self-force instead increases the rate of periapsis advance.The region where this occurs grows larger as a increases The black crosses mark the location of the underlying data from Ref.[88] used to calculate the contour plot.