Self-forced inspirals with spin-orbit precession

We develop the first model for extreme mass-ratio inspirals (EMRIs) with misaligned angular momentum and primary spin, and zero eccentricity -- also known as quasi-spherical inspirals -- evolving under the influence of the first-order in mass ratio gravitational self-force. The forcing terms are provided by an efficient spectral interpolation of the first-order gravitational self-force in the outgoing radiation gauge. In order to speed up the calculation of the inspiral 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. The resulting solutions are defined with respect to `Mino time' so we perform a second averaging transformation so the inspiral is parametrized in terms of Boyer-Lindquist time, which is more convenient of LISA data analysis. We also perform a similar analysis using the two-timescale expansion and find that using either approach yields self-forced inspirals that can be evolved to sub-radian accuracy in less than a second. The dominant contribution to the inspiral phase comes from the adiabatic contributions and so we further refine our self-force model using information from gravitational wave flux calculations. The significant dephasing we observe between the lower and higher accuracy models highlights the importance of accurately capturing adiabatic contributions to the phase evolution.


I. INTRODUCTION
With the advent of spaced-based gravitational wave (GW) detectors such as Laser Interferometer Space Antenna (LISA) [1], TianQin [2], and Taiji [3] in the 2030s, comes the need to model sources in the millihertz frequency range.A prime source for these future detectors are extreme mass ratio inspirals (EMRIs).These systems consist of a massive black hole (MBH) primary with mass M and a stellar-mass compact object secondary with mass µ (e.g., a steller-mass black hole or neutron star).These space based detectors will be to be sensitive to EMRIs with mass ratios ϵ = µ/M ∼ 10 −4 −10 −7 .Gravitational waves from these systems are expected to stay in the LISA band for months to years with the secondary completing ∼ 1/ϵ = 10 4 − 10 7 orbits spiralling around the primary before merging under the influence gravitational radiation reaction.Accurately modelling these systems through such a large number of orbital cycles to sub-radian accuracy with methods efficient enough for parameter estimation remains an open challenge, but one that will grant precise parameter estimation for MBHs along with some of the most sensitive probes for new physics beyond general relativity [4][5][6].
In this work, we look at the subclass of these systems which have zero eccentricity and misaligned orbital angular momentum and primary spin.In the strong field regime this will cause the orbital plane to rapidly precess around the primary's spin, ergodically filling (part of) a sphere, if we were to ignore orbital evolution due the emission of gravitational waves.Consequently, such systems are also known as quasi-spherical inspirals.One probable EMRI formation channel is the tidal disruption of a comparable mass black hole binary by a nearby MBH which will result in a circular and arbitrarily inclined inspiral by the time it enters the LISA band [7][8][9].While the primary EMRI formation channel of direct capture is expected to produce LISA band inspirals with eccentricity and inclination, radiation reaction will cause an EMRI to circularize (with the exception of right before the last stable orbit) [10] and cause its inclination to slowly increase [11].Thus an EMRI with any amount of inclination but low eccentricity is likely to evolve into a quasi-spherical inspiral before plunge.
The leading order behaviour of quasi-spherical EMRIs have been modelled by calculating the asymptotic fluxes at infinity and through the horizon of the primary and relating these to the averaged rate of change of energy, angular momentum, and Carter constant of the system [12][13][14].This same approach has recently been generalized to inspirals with inclination and eccentricity [14].While such inspiral calculations are numerically efficient and maybe accurate enough for detection of loud EMRI signals, the resulting adiabatic (0PA) inspirals do not reach the sub-radian accuracy requirement needed to reach all of LISA's science goals [15].
For that, one must understand the sub-leading order in mass ratio behaviour and compute post-adiabatic (1PA) inspirals.This will require knowledge of the local force on the secondary due to the presence of its own gravitational field, i.e., its gravitational self-force (GSF).This is computed by expanding the metric of the binary as the Kerr metric of the primary plus perturbations which are of expressed as a power series of the small mass ratio ϵ [16][17][18].Using a two-timescale analysis, it has been determined that to reach sub-radian accuracy across the inspiral, one will need to know all first-order in the mass ratio effects of the GSF as well as the orbit averaged dissipative contri-bution from the second-order in the mass ratio GSF [19].While second-order results for Schwarzschild spacetime are now emerging [20][21][22][23][24], calculations for inspirals in Kerr spacetime are likely still a few years away.As such, we do not include them in this work, but they can be incorporated into our computations once they are known.However, we do have a code that computes first-order GSF in Kerr spacetime for both eccentric equatorial orbits [25] and generic orbits [26] in the outgoing radiation gauge.We have adapted our code to compute the first ever GSF results for spherical orbits.
In order to drive inspirals, we require a model for the GSF can be rapidly evaluated for any point in the parameter space.This is typically done by calculating the GSF at many points and then fitting or interpolating the resulting data [27,28].We opt to interpolate our self-force data using Chebyshev polynomials, a technique which has proven to be very effective for EMRI calculations [29,30].We grid the parameter space on a 18 × 9 Chebyshev nodes for prograde orbits where the primary spin is 0.9M and for inclination angles up to 45 • with respect to the equatorial plane, and interpolate the data using Chebyshev polynomials to sub-percent accuracy.
With an interpolated model in hand, we can then compute the inspiral trajectory using the method of osculating geodesics (OG) [31,32].We note that the OG equations for generic Kerr inspirals derived in either Ref. [32] or Ref. [29] become singular in the spherical case.This is not a physical singularity and in this work we derive OG equations in this limit that are finite (see Appendix A).In this paper we explicitly assume the eccentricity is always zero so that the resulting inspiral model depends only on spherical-geodesic first-order GSF data.Neglecting terms proportional to eccentricity in our model incurs an error at 1PA order in the waveform phase but this is acceptable as we exclude other terms that contribute at this order, such as the currently unknown second-order GSF.The main goal for this work is not to produce a complete 1PA inspiral model, but to demonstrate how such a model can made computationally efficient without significant loss of accuracy.
Using OGs equations can take minutes to hours just to compute a single inspiral which is prohibitively slow for LISA data analysis.This is due to the OG equations' dependence on the rapidly oscillating orbital phases, which forces a numerical solver to take small time steps to resolve all ∼ 10 4 − 10 7 orbital cycles of a typical EMRI.To overcome this, we employ the technique of near identity (averaging) transformations (NITs) to find equations of motion that average out the dependence on the orbital phases while faithfully capturing the long term behaviour of the binary.This framework was first sketched out for generic EMRI systems and applied to eccentric Schwarzschild inspirals in Ref. [33], hereafter paper I, and then extended to equatorial, eccentric Kerr inspirals in Ref. [29], hereafter paper II.In both cases, the resulting inspirals retained sub-radian accuracy when compared to the OG equations and could be computed in less than a second.However, in both of these works, the resulting inspirals were not parametrized in terms of Boyer-Lindquist coordinate time t (the time of an asymptotic observer), which is much more convenient for data analysis as this is directly related to the time at the detector.Overcoming this required an additional interpolation step which increases the overall computation time.Following a procedure presented in Ref. [18], one can parametrize the averaged equations of motion in terms of t by performing an additional averaging transformation.We apply this improved NIT to the case of quasi-spherical Kerr inspirals and find that the inspiral calculation retains sub-radian accuracy and sub-second computation time, whilst being parametrized in a form that is more convenient for practical waveform production.
Much of the work done to model EMRIs to postadiabatic order makes use of the two-timescale framework [18,19,34] which has been long known to operate in an almost equivalent way to NITs [35] when applied to the equations of motion.Inspired by Ref. [18], we perform a two-timescale expansion (TTE) on the NIT equations of motion and show that the TTE equations of motion also produces fast inspirals that are accurate to post-adiabatic order but with the drawback of doubling the number of equations to numerically solve.On the other hand these equations are independent of the mass ratio and so, after being solved once (in under a second), inspirals for any mass ratio can be generated in a fraction of a second.This could provide extra efficiency for an appropriately designed search or parameter estimation algorithm.The TTE equations of motion can also be used in conjunction with the two timescale expansion of the field equations to produce complete post-adiabatic waveforms [23].
Unlike the OG equations, both the NIT and TTE equations of motion cleanly separate the adiabatic and postadiabatic contributions.As such, it is straightforward to increase the overall accuracy of the inspiral by calculating the adiabatic contributions from computationally cheaper, and generally more accurate, gravitational wave flux calculations and interpolating them to higher precision.Since this is the leading order contribution to the inspiral evolution, it will need to be known with a fractional accuracy of ∼ ϵ [28,36].We achieve this using a 18 × 19 Chebyshev grid which covers all inclinations and reaches a fractional accuracy in the adiabatic pieces of ≲ 10 −5 .Using these higher precision interpolants for the adiabatic pieces along with the GSF model for the postadiabatic pieces results in an improvement to the phase ranging from 10 − 10 4 radians, highlighting just how important high precision flux calculations are to both adiabatic and post-adiabatic inspirals [14,37].
With our more accurate quasi-spherical inspiral model, we examine the effects of the self-force on quasi-spherical Kerr inspirals.We find that the evolution is consistent with what was known at adiabatic order and that the post-adiabatic contributions of the first-order GSF become more important as the inclination of the orbit in-creases.
This paper is organized as follows.In Section II, we introduce spherical geodesic orbits before giving an overview of the method of osculating geodesics.In Section III, we briefly review the gravitational self-force approach, discuss the modifications made to the Kerr GSF code of Ref. [26] for spherical orbits and outline our interpolation scheme for both the gravitational wave fluxes and the GSF.In Section IV, we review the technique of NITs for generic EMRI systems before describing an additional averaging transformation used to obtain averaged equations of motion parametrized by t.We then perform a two timescale expansion on these averaged equations of motion before applying these frameworks to the case of quasi-spherical inspirals.In Section V, we outline our numerical implementation for calculating and using these averaged equations of motion.We demonstrate the results of this implementation in Section V where we start by comparing inspirals and waveforms calculated using either OG, NIT or TTE equations of motion.We then show the effect of using higher precision adiabatic pieces from flux calculations before finally exploring the effects of the self-force across a subsection of the spherical Kerr parameter space.We conclude with some discussion in Section VII.We detail how one can derive OG equations in the quasi-spherical limit in Appendix A and how the phases used in the t parametrized equations of motion relate to the waveform phases in Appendix B.
Throughout this paper we work in geometrical units where G = c = 1.

II. QUASI-SPHERICAL INSPIRALS AROUND A ROTATING BLACK HOLE
In this section we describe the motion of a non-spinning compact object of mass µ on a spherical orbit of Boyer-Lindquist radius r p , with arbitrary inclination z p = cos θ p in the Kerr spacetime under the influence of some arbitrary force.We model the secondary as a point particle with a position given in modified Boyer-Lindquist coordinates by x α p = {t p , r p , z p , ϕ p }.Note that hereafter, we use a subscript 'p' to denote a quantity evaluated at the particle's location.Later in this work, we will take this perturbing force 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 as where ∆ := r 2 + a 2 − 2M r, Σ := r 2 + a 2 z 2 , and ϖ := √ r 2 + a 2 .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. ( 2) into a form that is useful for applying the near-identity transformations.Before considering the forced equation it is useful to first examine the geodesic limit.

A. 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 given by where E is the orbital energy per unit rest mass µ, L z is the z-component of the angular momentum divided by µ, Q is the Carter constant divided by µ 2 [38], and ).This definition of the Carter constant is related to another common definition of the Carter constant, K, by The geodesic equation can be written explicitly in terms of these integrals of motion [38,39]: where B = Eϖ 2 − aL z , the radial potential V r = 0 for spherical and circular orbits, and z + > z − are the roots of the polar potential V z .We make use of (Carter-)Mino time λ which is related to the secondary's proper time τ , by dτ /dλ = Σ p , as this decouples the radial and polar equations [40].While this is not strictly necessary in the spherical case as there is no radial motion, using λ as our time parameter allows us to take advantage of the analytic solutions for generic Kerr geodesics derived in Ref. [41], and allows for a more uniform treatment with paper II and the generic case in the future.Rather than parametrize an orbit using the integrals of motion ⃗ J = {E, L z , K}, we find it more convenient to work with the constants ⃗ P = {r p , e = 0, x}.Here x is a measure of the orbital inclination given by x = cos θ inc .The inclination angle θ inc is related to θ min (the minimum value of θ p measured with respect to the primary's spin axis) by θ inc = π/2 − sgn(L z )θ min .Explicit relationships between the integrals of motion ⃗ J in terms of ⃗ P for spherical orbits are too lengthy to be expressed here, but they were first derived in Ref. [12] and can be found in the KerrGeodesics package as part of the Black Hole Perturbation Toolkit [42].
There are other common choices for inclination in the literature such as the inclination angle ι [11][12][13]43] given by cos ι = L z / L z 2 + Q.Both inclination angles smoothly parametrize all inclinations between prograde equatorial motion where x = 1 = cos ι to retrograde equatorial motion where x = −1 = cos ι.However, we opt to use x as it has a simple relation to the roots of the polar potential V z via It is worth noting that not all values of {r p , x} correspond to bound geodesics and we denote the position of the innermost stable spherical orbit (ISSO) as r ISSO (a, x) [10,44].
In order to later apply the near-identity (averaging) transformations, it will be useful to employ action-angle formulation to parametrize the geodesic motion [41,45,46].In this description the orbital phase ⃗ q = {q z } is such that the geodesic equations can written in the form where Υ z is the Mino time fundamental polar frequency.This is known analytically [41] and is given by where k z = Az 2 − /z 2 + , K is the complete elliptic integral of the first kind with K(m) = F (π/2|m) and F is the incomplete elliptic integral of the first kind given by One can also analytically express the polar coordinate z in terms of r, x and q z via where sn is the Jacobi elliptic sine function given by sn(u|m) = sin(am(u|m)) and the Jacobi amplitude am is the inverse function of F .At this point we note that it is also common in the literature [47,48] to express this polar coordinate in terms of a quasi-Keplerian angle χ z , via However, the rate of change of χ z depends on χ z itself, which will be inconvenient in later sections as we look to derive averaged equations of motion that have no dependence on the orbital phases.

B. Osculating Geodesics for quasi-spherical inspirals
To go beyond the geodesic orbits and describe quasispherical inspirals, we make use of the method of osculating orbital elements (or osculating geodesics) [31,32] to recast the equations of motion of a body under the influence of a force obeying Eq. ( 2).
We chose a set of parameters that uniquely specify a geodesic orbit, 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 }1 .For spherical geodesics, we chose ⃗ I = (r p , x, q z,0 ).For accelerated orbits, these orbital elements are promoted from constants to functions of time.
Here we are implicitly assuming that the eccentricity e stays zero throughout the inspiral of an initially spherical system.While this has been shown to be valid at 0PA order [11,49], at 1PA this becomes a subtle issue [50,51].Whether the eccentricity remains zero throughout the post-adiabatic evolution is intimately related to the exact setup for the calculation of other post-adiabatic corrections such as the second-order fluxes [21].For the sake of simplicity, and to have an inspiral model that depends only on spherical-geodesic first-order GSF data, in this paper we explicitly assume the eccentricity is always zero.In doing so, we are neglecting terms that contribute at the same order as the currently unknown second-order GSF.
To produce waveforms, we will also require the evolution of certain "extrinsic quantities" ⃗ S = {t p , ϕ p } which are the secondary's Boyer-Lindquist time and azimuthal coordinates respectively.The combined osculating geodesic evolution equations take the following form: where f i = dq 0,i /dλ.For generic Kerr inspirals, this was implemented in Ref. [32] using quasi-Keplerian angles (e.g., χ) as the orbital phases.In order to derive averaged equations of motion, we implemented an alternative version in Paper II [29] which makes use of geodesic action angles ⃗ q for the orbital phases and the integrals of motion.In Appendix A, we derive the explicit expressions for spherical inspirals.Combining these modified evolution equations with the rest of the evolution equations derived in Paper II [29] gives us equations of motion for quasi-spherical inspirals under the influence of the gravitational self-force with the form dr p dλ = ϵF (1)  r (r p , x, q z ) + ϵ 2 F (2) r (r p , x, q z ) + O(ϵ 3 ), (15a) dx dλ = ϵF (1)  x (r p , x, q z ) + ϵ 2 F (2)  x (r p , x, q z ) + O(ϵ 3 ), (15b)

III. GRAVITATIONAL SELF-FORCE FOR QUASI-SPHERICAL KERR INSPIRALS
The gravitational self-force approach consists of expanding the metric of the binary around the metric of the primary, i.e., g µν = ḡµν + ϵh µν + . . .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.Using matched asymptotic expansions, one can then derive how the interaction between these metric perturbations and their source effects the motion of the secondary [17].

A. Gravitational self-force for spherical orbits
This work represents the first calculation of the gravitational self-force for spherical orbits in Kerr spacetime.While the gravitational self-force had been previously calculated for generic (eccentric and inclined) orbits in Kerr in [26], building on previous developments in [25,[52][53][54], and the scalar self-force for spherical orbits had been calculated in [55], the case of the gravitational self-force for spherical orbits has not been explored before now.
To calculate the gravitational self-force for spherical orbits, we follow the same method as was used in the generic orbit case in [26], which can be adapted without major adjustments.We summarize our method here.We start by solving the Teukolsky equation for the Weyl scalar ψ 4 in the frequency domain using a numerical implementation [56,57] of the semi-analytical Mano-Suzuki-Takasugi method [58,59].From ψ 4 , we obtain the Hertz potential by algebraically inverting the fourth order differential equation relating it to ψ 4 [54,60] mode-by-mode.From the modes of the Hertz potential we can reconstruct the modes of the local metric perturbation in the outgoing radiation gauge [61][62][63], from which, in turn, we obtain the modes of the gravitational self-force.Initially these modes are expressed in spin-weighted spheroidal harmonics and their derivatives.To facilitate calculating the regular part of the GSF modes using mode-sum regularization [64][65][66][67], the spheroidal harmonics are first projected to spin-weighted spherical harmonics using the method introduced in [12], which in turn are projected onto scalar spherical harmonics [25,26,54].The mode-sum is accelerated by a tail fitting procedure leading to an exponential convergence of the mode-sum [54].
The above metric reconstruction procedure captures most of the metric perturbation, except for a perturbation of the background Kerr spacetime's mass and spin, and a purely gauge contribution.The mass and spin completions are directly related to the energy and angular momentum of the orbit [68,69].The gauge completion is more involved.The version of the outgoing radiation gauge used for the GSF calculation, the no-string gauge [53], has a distributional singularity on sphere of the orbit.Recently, [70] has proposed a generalization of our method that produces the metric perturbation in a gauge that is sufficiently smooth to allow it to be used as a source for a second order calculation.We here strive for a less lofty goal and will add the minimal gauge completion needed to ensure orbit averaged observables like the frequency-shifts are invariant.This involves fixing a 2-parameter gauge freedom representing a rescaling of time and time dependent coordinate rotations inside the orbit.In [70] a nice derivation of the needed coefficients was given.However, for this work we used an implementation [71] of a generalization of the method used for the mass and spin completion in Ref. [68].The gauge vector representing this freedom is given by, with the gauge fixing parameters α and β given by where the numerators of the integrands, I α and I β are given in Appendix C, and ∆ p and Σ p denote ∆ and Σ evaluated at the particle location.
It is useful to split the self-force into dissipative (time anti-symmetric) and conservative (time-symmetric) contributions [67].The dissipative pieces cause the orbit to shrink until the secondary plunges into the primary.They also have a small effect on the orbital inclination, causing the orbit to become more inclined over time [11][12][13]43].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 [10,14,37,72,73].The conservative pieces have more subtle effects on the inspiral, such as altering the rate of periapsis advance, the rate of nodal precession, and the location of the innermost stable circular orbit [74][75][76][77][78][79][80].
Computing 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 second-order self-force [19].There are as yet no calculations of the latter in Kerr spacetime, so we will make do with only the first-order gravitational wave fluxes and the first-order self-force in this work, though this means that the resulting "post-adiabatic" inspirals and waveforms will be gauge dependent until the missing second-order contributions are added [29].
In order to drive the OG equations of motion to calculate the inspiral, we require a model for the GSF that can be rapidly evaluated at each time step during the inspiral.This is typically done by tiling the parameter space with with GSF data which is then either fitted to a model or interpolated.While this has been done in a variety of ways for eccentric orbits in Schwarzschild [27,28] and Kerr [29], we will now describe the first interpolated GSF model for quasi-spherical Kerr inspirals.

B. Interpolated gravitational wave fluxes for quasi-spherical Kerr inspirals
To introduce our interpolation procedure, we will first interpolate the energy and angular momentum fluxes for spherical orbits.Since flux calculations, which can be obtained directly from ψ 4 , are significantly cheaper than calculating the GSF, it is much more feasible to densely tile a large section of the parameter space with flux data.This in turn will result in more accurate interpolation of the leading-order, adiabatic effects.It will also allow us to carry out consistency checks on our GSF model by comparing the fluxes with the orbit averaged GSF.
We start by fixing the value of the spin parameter of the primary, which we choose to be a = 0.9M .This reduces our parameter space to two parameters; the orbital radius r p and the inclination x.We then define a parameter y using r p and the position of the innermost stable spherical orbit r ISSO .We choose y to be Tiling the parameter space with y instead of r p will concentrate more points near the ISSO where the fluxes and the GSF experience the most variation.We let y range from 0 to 1 and x range from −1 to 1 thus covering all inclinations for both prograde and retrograde orbits.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 fluxes on a 18 × 19 grid of Chebyshev nodes, with the y values given by the roots of the 18th order Chebyshev polynomial and the x values given by the roots of the 19th order Chebyshev polynomial.The size of this grid was chosen empirically until the interpolant reached the desired accuracy across the parameter space.
At each point on this grid, we calculate the energy and angular momentum flux both at infinity and at the event horizon.From these, one can calculate the leading order orbit averaged rate of change of the energy and angular momentum of the secondary via the following balance laws [40,[81][82][83][84]: where the orbit average of a quantity A is given by We could interpolate the rates of change of energy and angular momentum, but we find it more convenient to work with the orbital elements r p and x.As such, we find their rates of change via the chain rule: = ϵΓ (1)  x + O(ϵ 2 ).
The partial derivatives can be found using the analytic expressions for E(r p , x) and L z (r p , x) from Ref. [12] to construct the Jacobian This can then be inverted to give To improve the accuracy of our interpolation, we rescale the data for Γ  x , times a term that is zero for the limiting cases of either the separatrix or the equatorial plane respectively.Finally, we take a discrete cosine transformation of the data on our Chebyshev grid to obtain the Chebyshev polynomial coefficients C ij rp/x .Summing these coefficients together with Chebyshev polynomials gives us the following interpolants for Γ (1) r and Γ (1) x : Γ (1)  x = 1 Using the largest coefficient for i = 17 and j = 18 to estimate the relative error, we infer that these interpolants should have a relative error of ∼ 10 −6 .
To test the accuracy of this interpolation, we compare the interpolants to data on a grid that were not used for the interpolation while using the relations p − (∂E/∂x) Γ (1) x and x .In Fig. 1, we see that the interpolants match the energy and angular momentum fluxes to a relative error of ≲ 10 −5 with the exception of points very close to the ISSO.In principle this model interpolates out to r p = ∞ and the leading PN scaling should give it the right behaviour in the weak field.We would expect the model to retain a comparable level of accuracy at large r p , but it has not been tested beyond r p ∼ 30 as, in this work, we are more interested in the strong field dynamics.

C. Interpolated gravitational self-force model for quasi-spherical Kerr inspirals
We now use a similar interpolation scheme to create a model for the gravitational self-force that is continuous throughout the parameter space and fast to evaluate which can be used with the osculating geodesic equations to describe quasi-spherical Kerr inspirals.However, given the cost of our GSF code -the computation of the GSF on a single geodesic can take hundreds of CPU hours -we restrict ourselves to a 2D slice of the EMRI parameter space using only a few hundred points.Once again, we restrict a = 0.9M and let y range from y min = 0.1 to y max = 1, but instead of x, we opt to tile in z 2 − = 1 − x 2 and let z 2 −,min = 0 to z 2 −,max = 0.5.
This allows us to cover inclination angles ranging from 0 • ≤ θ inc ≤ 45 • , in a way that is easy to rescale for Chebyshev interpolation.We define parameters u and v which cover this parameter space as they range from (−1, 1) u := y − (y min + y max )/2 (y min − y max )/2 and, (27a) We then calculate the GSF on a 18 × 9 grid of Chebyshev nodes, with the u values given by the roots of the 18th order polynomial and the v values given by the roots of the 9th order polynomial.At each point on our grid, we Fourier decompose each component of the force with respect to the polar action angle q z .As inclination increases more and more Fourier modes of the GSF become relevant.However, the number of relevant modes stays finite even for polar orbits.To build a Chebyshev interpolant we need to resolve all Fourier modes at all grid points.Resolving the higher order harmonics at low inclination grid points is numerically challenging and therefore computationally expensive.This is the main reason for limiting the range of inclinations to 0 • ≤ θ inc ≤ 45 • as this limits the number of Fourier modes that need to be resolved.
To smooth the behaviour of the force near the separatrix and improve the accuracy of our interpolation, we then multiply the data for each Fourier coefficient by a factor of (1−y) 2 .Next, we use Chebyshev polynomials to interpolate each Fourier coefficient across the (u, v) grid.We then sum the modes to reconstruct our interpolated gravitational self-force model: We note that this choice of rescaling forces each component to become singular at the ISSO, and while the components of the GSF change rapidly as one approaches the ISSO, we still expect them to be finite at the ISSO.A greater understanding of the analytic structure of the GSF in this region would greatly improve this and any future interpolated GSF models.
We also 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.Since the OG equations are derived assuming this condition to be true [29,32], we project the force so that this condition is always satisfied, i.e., To verify the accuracy of our interpolated model, we employ the flux balance laws to compare the local energy and angular momentum lost by the secondary against the energy and angular momentum fluxes radiated at infinity and down the horizon: where Υ and Υ (0) τ = ⟨Σ p ⟩ are the fundamental Mino time frequencies for the secondary's Boyer-Lindquist time coordinate and proper time with analytic expressions found in Ref. [41] and appendix C of Ref. [46] respectively.
Fig. 2 shows the relative error in the fluxes across the parameter space.The white circles indicate the points used to generate the interpolated GSF model, whereas the green crosses indicate the points in the parameter space where the fluxes were calculated and thus where the comparisons were made.We find that the relative error in the fluxes only rises to ∼ 10 −2 when close to the ISSO, and it is typically of the order ∼ 10 −3 − 10 −5 .This is comparable to previous methods which required tens of thousands of points to achieve the same level of accuracy [27,28].While this would not be sufficiently accurate for LISA data analysis which would require fluxes accurate to ∼ ϵ, such an interpolating error may be permissible for the post-adiabatic contributions of the first order GSF [28,36].We will go on to improve this model by incorporating the higher precision flux interpolants in Sec.VI C.
We can now use this model in conjunction with the OG equations of motion (15) to calculate inspiral trajectories.However, we find these trajectories take minutes to hours to compute, due to the need to resolve hundreds of thousands of orbital cycles.We will now look to leverage averaging transformations which can remove the dependence of the orbital phases from the equations of motion while retaining the accuracy required to produce post-adiabatic EMRI waveforms.

IV. AVERAGING TECHNIQUES FOR EMRI EQUATIONS OF MOTION
NITs are a well known technique in applied mathematics and celestial mechanics [35].This approach involves making small transformations to the equations of motion, such that the short timescale physics is averaged over, while retaining information about the long-term evolution of a system.In paper I [33], these transformations were derived for a generic EMRI system which we review in Section IV A below.
As we have seen above, for Kerr geodesics it is helpful to parametrize the motion by Mino time, λ.How-ever, for data analysis one requires the final waveforms to be parametrized by the time at the detector.Naively, transforming back to Boyer-Lindquist time, t, requires either numerically interpolating λ(t), or adding orbital timescale oscillations back into our NIT equations of motion.The former additional step is unwelcome in the pursuit of fast waveform generation, and the latter would remove all the benefit of the averaging transformation.Fortunately, it was noted in Ref. [18] that these oscillations can be averaged out if one performs an additional averaging transformation.We outline the details of this procedure in Section IV B.
NITs are not the only averaging procedure that can be used to accelerate the inspiral calculation.In Section IV C we present a two-timescale expansion of the equations of motion.One advantage of this approach is that the adiabatic and post-adiabatic inspirals can be computed offline.Given a particular mass ratio, the true inspiral can then be computed from interpolating functions as a very fast online step.It may be advantageous to design data analysis pipelines to leverage this additional speed-up.
We conclude this section by presenting averaged equations of motion for the case of quasi-spherical Kerr inspirals in Section IV D.

A. Review of near-identity 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 Pj = P j + ϵY (2) j ( ⃗ P , ⃗ q) + O(ϵ 3 ), (31a) qi = q i + ϵX (1) i ( ⃗ P , ⃗ q) + ϵ 2 X (2) i ( ⃗ P , ⃗ q) + O(ϵ 3 ), (31b) Here, the transformation functions i , and Z (n) k are required to be smooth, periodic functions of the orbital phases ⃗ q.At leading order, Eqs.(31) 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 k .The inverse transformations can be found for P k and q i by requiring that their composition with the transformations in Eqs.(31) must give the identity transformation.To first order in ϵ, this gives us It will be 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 ⟨A⟩ ( ⃗ P ) and an oscillating piece Using the above transformations along with the equations of motion, and working order by order in ϵ, we can choose values for the transformation functions i , X i , Z k and Z (1) k such that the resulting equations of motion for Pj , qi and Sk take the following form i Crucially, these equations of motion are now independent of the orbital phases ⃗ q.Deriving the relationship between the transformed forcing functions, α ), and the original forcing functions is quite an involved process with several freedoms and choices, each with its own merits and drawbacks.This is discussed at length in paper I, 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 Υ Note that since we currently lack second order in mass ratio contributions, we will set F (2) j = 0.In deriving these equations of motion, we have constrained the oscillating pieces of the NIT transformation functions to be X( 1) k is found by solving This equation is satisfied by the oscillating pieces for the analytic solutions for the geodesic motion of t and ϕ, We have chosen the averaged pieces Y = 0 for simplicity, though one could make other choices like in Ref. [18].In order to generate waveforms, one only needs to know the transformations in Eq. ( 31) 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(ϵ) throughout an inspiral, this requires the transformation of the P j 's (31a) to linear order in ϵ, while it is sufficient to know the rest of Eq. ( 31) to zeroth order.

B. Averaging transformations for motion parametrized by Boyer-Lindquist coordinate time
Solving the above equations results in solutions for ⃗ P , ⃗ q and ⃗ S as functions of Mino time λ.While this would include t(λ), the transformation to λ(t) is non-trivial, and in practice it is done via interpolation which can be costly for long inspirals.It would be significantly more convenient for the solutions to be functions of t from the start so that one can produce waveforms for data analysis without this postprocessing step.This can be accomplished for the OG equations by simply using the chain rule: Notice that we have one less equation of motion to solve.However, using the same approach to the NIT equations of motion results in ϕ ( ⃗ P ) .
To obtain the equations of motion for P j and φ i , we take the time derivative of Eq. ( 44), substitute in the expression for the NIT equations of motion, and then use the inverse transformation of Eq. ( 44) to ensure that all functions are expressed in terms of ⃗ P and ⃗ q and expanding order by order in ϵ.We then choose the oscillatory functions ∆t, Φ i , Π j and Π (2) j to cancel out any oscillatory terms that appear at each order in ϵ.
This results in averaged equations of motion that take the following form: (2) These equations of motion are related to the Mino time averaged equations of motion, Eqs.(35), where the adiabatic terms are given by Γ (1) and the post-adiabatic terms are given by Γ (2) t Ω (1) i .

(47b)
This constrains the oscillating pieces of our transformation to be ∆t Π(1) Φ (1) t Ω (0) t We are free to chose the averaged pieces of Π j , and we make the simplification that Π  and Ω (1) t Ω (0) α .

(49b)
What is most useful about these equations of motion is that their solutions ⃗ P(t) and ⃗ φ(t) is exactly what is required to feed into waveform generating schemes.We show the equivalence between ⃗ φ(t) and the relationship derived in Ref. [85] between the solutions to the original NIT equations of motion and the waveform phases Φ mn in Appendix B.

C. Two-timescale expansion
There is a related way of obtaining the above solutions via the TTE.We exploit the difference between the timescales of the system by defining T := ϵt as the slow time which governs the long-term, secular behaviour of the system and defining t as the fast time of the system that governs the short-term, orbital dynamics and treating these two times independent variables.
As such, we expand the transformed variables as Applying this expansion to the t parametrized NIT equations of motion, one finds that the equations of motion for the two-timescale expanded variables takes the form: (1) j dT = Γ (2) dφ There is a trade-off for solving these equations of motion.We now have to solve a system of coupled differential equations that is twice the size and thus is more expensive to solve numerically, but the solutions are independent of ϵ and one can construct a solution for any given value of ϵ using Eqs.50.Thus, if one wants to compute multiple inspirals with varying mass ratios, the TTE can be more efficient overall.However, there is also an issue where the inspiral will stop earlier, as the variables P (0) j typically reach values which correspond to the ISSO before P j do.In this regime, one should instead employ a transition to plunge as this is where the adiabaticity assumptions of the OG equations and the two-timescale expansion is expected to break down.As such we will use both the NIT and the TTE equations of motion to produce waveforms and compare them to waveforms generated using the OG equations to assess which is the more practical framework for producing post-adiabatic EMRI waveforms.

D. Averaged equations of motion for quasi-spherical Kerr inspirals
We now specialize the above results to the case of quasi-spherical inspirals into a Kerr black hole.The NIT equations of motion parameterized by Mino time take the form: x (a, rp , x), (52b) z (a, rp , x) + ϵΥ (1)  z (a, rp x), (52c) t (a, rp , x), (52d) The leading order terms in each equation of motion are simply the original functions, derived in Appendix A and paper II, averaged over a single geodesic orbit, i.e., F where Υ (0) t and Υ (0) ϕ are the Mino time t and ϕ fundamental frequencies which are known analytically [41].The remaining terms are more complicated and require Fourier decomposing the original functions and their derivatives with respect to the orbital elements (r p , x).To express the result, for any function we define the operator With this in hand, the remaining terms in the equations of motion are found to be x = N (F (1)  x ), Υ (1) Combining these results with Eqs.(45), one can find the NIT equations of motion parametrized by Boyer-Lindquist time t for the phases ⃗ φ = {φ z , φ ϕ } and orbital elements ⃗ P = {r φ , x φ } in the form The leading order terms in these equations are given by The sub-leading terms are given by t Γ (1)   x , (58b) Finally, using the two-timescale expansion the adiabatic equations of motion are given by x , (59a-b) The post-adiabatic contributions to the equations of motion are given by (1) z dT = Ω (1)  p + r (1)   φ ϕ + r (1)   φ Solutions for post-adiabatic inspirals can be obtained by solving Eqs. ( 59) and ( 60) simultaneously and using Eqs.(50) along with a value for the mass ratio ϵ to recover r φ (t), x φ (t), φ z (t) and φ z (t).
With the averaged equations for quasi-spherical Kerr inspirals in hand, we will now outline our numerical implementation for rapidly computing quasi-spherical selfforced inspirals.

V. IMPLEMENTATION
Combining the interpolated GSF model along with our action-angle formulation of the OG equations provides us with all the information required to calculate the NIT and TTE equations of motion.We first evaluate and interpolate the various terms in the NIT and TTE equations of motion across the parameter space.While this offline process can be expensive, it only needs to be completed once.Once completed, the online process of calculating self-forced inspirals can be completed in less than a second.This procedure is very similar to Papers I & II, though now we also export interpolating functions for the partial derivatives needed for the TTE equations of motion.

A. Offline Steps
The offline calculation consists of the following steps.
1. We start by selecting a grid which covers the parameter space.We choose y values between 0.099 and 0.999 in 451 equally spaced steps and z 2 − values from 0.002 to 0.5 in 250 equally spaced steps (giving a total of 112,750 points). 2   2. For each point in the parameter space (a, y, z 2 − ) we evaluate the functions t\ϕ along with their derivatives with respect to r p and x for 49 equally spaced values of q z from 0 to 2π.
3. We then perform a fast Fourier transform on the output data to obtain the Fourier coefficients of the forcing functions and their analytical derivatives with respect to r p and x.
4. With these, we then use Eqs.(53a-c), (53d-e), ( 55) and (55a-d) to construct r\x /∂x, and Υ t\z\ϕ at that point in parameter space.All other terms needed for the NIT and TTE can be derived from these terms or are already known analytically.
5. We also use Eqs.(37) to construct the Fourier coefficients of the first-order transformation functions Y (1) r\x .These are needed when comparing NIT and OG/TTE inspirals, to ensure that the initial conditions are comparable.Otherwise this step can be skipped.
6.We then repeat this procedure across the parameter space for each point in our grid.

Finally we interpolate the values for
r\x /∂p, and Υ t\z\ϕ along with the coefficients of Y (1) r\x across this grid using Hermite interpolation and store the interpolants for future use. 2 Evaluating the NIT functions is computationally cheap so using a dense grid does not significantly increase the computational burden.Using an equally spaced 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 so that the interpolation error is a negligible source of error for our comparisons between the OG, NIT and TTE inspirals, though a less dense grid may also achieve this.
We implemented the above algorithm in Mathematica 12.2 and find the calculation takes about 5 hours to complete when parallelized across 20 CPU cores.Since these offline steps need only be completed once, this is a comparatively small price to pay.

B. Online Steps
By contrast, the online steps are required for every inspiral calculation, but are computationally inexpensive.The online steps for computing a NIT or TTE inspiral are as follows.
1. We load in the interpolants for F (1\2) r\x and Υ t\z\ϕ , and define the NIT equations of motion.For TTE inspirals, one also needs to load in the interpolated derivatives ∂ F (1) r\x /∂r p and ∂ F (1) r\x /∂x and then define the TTE equations of motion.

In order to facilitate comparisons between OG,
NIT, and TTE inspirals, we also load interpolants of the Fourier coefficients of Y (1) r/x and Eqs. ( 31), ( 44) and (48b) to construct first order near-identity transformations.
3. We state the initial conditions of the inspiral (r p (0), x(0), q z (0)) and use the NIT to leading order in the mass ratio to transform these into initial conditions for the NIT/TTE equations of motion, i.e., (r φ (0), x φ (0), φ z (0)).
4. We then evolve the NIT or TTE equations of motion using an ODE solver (in our 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(ϵ).

C. Waveform Generation
With a trajectory in hand, can now generate the gravitational waveform.Ideally, this would be done using interpolated Teukolsky fluxes as implemented by, e.g., the FastEMRIWaveform package for eccentric Schwarzschild inspirals [86,87].Unfortunately, such a model is not currently available for quasi-spherical Kerr inspirals, and so we generate all of our waveforms using the "semirelativistic" approximation used by the numerical kludge models [88].In this approach the Boyer-Lindquist coordinates of the particle are mapped to flat spacetime coordinates which are fed into the quadrupole formula that is used to produce the final waveform strain h = h + + ih × .This approximation fairs surprisingly well compared to Teukolsky snapshot waveforms, even in the strong field regime.What is most important for our purposes is that all waveforms are calculated using the same kludge formula so that any differences in the waveforms are due to differences in the calculations of the inspiral trajectory.For our implementation, we sample every t step = 2M and calculate the waveform strain at each time step to produce our numerical waveforms.
In order to calculate the fractional overlap O(h 1 , h 2 ) between two time-domain waveform strains h 1 and h 2 , we first define the noise weighted inner product of these two waveforms as where h(f ) denotes the Fourier transform of the time domain waveform h(t) and h * is the complex conjugate of h.We take the power spectral density (PSD) of the detector noise S n (f ) to be a flat noise curve.From this, one can calculate the fractional waveform overlap O via This overlap ranges from 1, when the two waveforms are identical, to 0, when the waveforms are perfectly orthogonal.When dealing with waveform overlaps close to 1, it is often useful to talk in terms of the fractional waveform mismatch M = 1 − O.In practice, we make use of the WaveformMatch function from the SimulationTools package to calculate waveform overlaps [89].

A. OG vs NIT and TTE inspirals
To test the accuracy of our NIT and TTE implementations, we compare inspirals calculated using the OG equations of motion against inspirals calculated using the NIT or TTE equations of motion.We use a case study of a typical system EMRI with a primary of mass M = 10 6 M ⊙ .We chose the initial conditions for this inspiral to have an inclination of x = 0.75 and radius r p = 7.75M .This highly inclined, strong field inspiral provides a good test of our numerical implementations.We see similar results for other initial conditions.We also set the initial phases q z (0) = ϕ(0) = 0 for simplicity.
Figure 3 demonstrates the difference in the trajectories produced by each method.We have purposely chosen an unusually large mass ratio for an EMRI of ϵ = 10 −2 FIG. 3. Trajectory through (rp, θinc) space for an inspiral with ϵ = 10 −2 , a = 0.9M , and initial conditions rp(0) = 7.75M, x(0) = 0.75.We use such a large mass ratio to highlight the orbital timescale oscillations one encounters when using the OG equations.Using the NIT or TTE equations of motion averages out these oscillations and results in almost identical inspirals.We also see that the TTE equations of motion break down further from the ISSO than the OG or NIT equations of motion.so that the orbital timescale oscillations are clearly visible on the plot.We see that for all trajectories the orbital separation decreases substantially, while the orbit becomes slightly more inclined.This is consistent with previous results for adiabatic quasi-spherical inspirals [11][12][13].The OG trajectory oscillates on the orbital timescale, while the NIT and TTE trajectories faithfully capture the average evolution of the trajectory.While the NIT and TTE trajectories may appear at first glance to be identical, it is important to note that the TTE equations of motion break down sooner than the NIT equations of motion.This is due to the evolution of r (0) φ and x (0) φ reaching the location of the ISSO sooner than r φ and x φ .This is not as important of an issue as one might expect, as this close to the ISSO one should swap over to a "transition to plunge" scheme [90][91][92][93].
Figure 4 shows the difference in the phases between the OG trajectory and the NIT and TTE trajectories for a year long inspiral with a primary of mass M = 10 6 M ⊙ and mass ratio of ϵ = 10 −5 .In both cases, we find that the difference in the phases stays below ∼ 10 −3 throughout the inspiral, spiking only when the inspiral approaches the ISSO where the adiabatic assumption implicit in the OG, NIT and TTE equations of motion starts to break down.Even then, the difference in the phases is substantially lower than subradian accuracy requirement needed for LISA data analysis.We also find that the growth in the error over time is most closely correlated with the interpolation error for the terms in the NIT and TTE equations of motion, and so interpolating on a denser grid should reduce this error even more.We note that formally both the NIT and the TTE should induce an error in both the phases and the orbital elements that scales linearly with the mass ratio.To ensure that our implementation is converging correctly, we evolve inspirals from initial conditions r p (0) = 5M and x(0) = 0.75 until they reach r p = 4M with different values of the mass ratio.Since the OG quantities will have oscillations of O(ϵ), sampling the error at only the final time step would make it difficult to determine the convergence with the mass ratio.As such we calculate the error at several time-steps in the last three orbital cycles and report the largest error.The results of this test can be seen in Fig. 5.We see that for larger mass ratios the errors converge linearly, as expected.However, for mass ratios ≲ 10 −4 , the formal error in the NIT and TTE is no longer the dominant source of error for the evolution of the phases.Instead the error is dominated by either interpolation error or the error in the ODE solver when solving the OG equations as we have found empirically that this error is reduced by increasing the number of grid points used for interpolation and decreasing the relative tolerance of the integrator.While both of these sources can be further suppressed at the cost of more expensive offline and online steps respectively, we show in Sec.VI B that it is not necessary for producing trajectories accurate enough for LISA data science.Now that we have established our averaging procedure does not change the accuracy of our model, we demonstrate the speedup that one enjoys from using either the NIT or TTE equations of motion instead of using the OG equations of motion in Table.I.For each of these calculations, the initial conditions are set to r p (0) = 7.75M and FIG. 5.The absolute difference in the quantities of a prograde inspiral with a = 0.9M and an initial inclination of x0 = 0.75 after evolving from rp = 5M to rp = 4M when using different evolution equations.The solid lines show the difference between OG and NIT equations of motion and the dashed lines show the difference between using the OG and TTE equations of motion.We see that the differences generally follow the black ϵ curve and so converge linearly with the mass ratio until we reach mass ratios ≤ 10 −3 where the error in the phases becomes dominated by interpolation and numerical error.TABLE I. Computational time required to evolve an inspiral from its initial conditions of rp(0) = 7.75M and x(0) = 0.75 to the last stable orbit for different values of the mass ratio, as calculated in Mathematica 13 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.The computation time for the TTE inspiral is 0.53s, which is slightly longer than any of the NIT inspirals.However, if we consider the TTE calculation as an offline step, one can then immediately recover the solution for any value of ϵ in a matter of milliseconds.
x(0) = 0.75, and the inspirals are evolved to just before the ISSO, when y(r p , x) = 0.998.All calculations were done using machine precision arithmetic and an accuracy goal of 7 digits for Mathematica's NDSolve function.We find that using the OG inspiral calculation takes longer as the mass ratio gets smaller since the solver will have to resolve many more orbital cycles before the inspiral reaches the ISSO.However, the NIT inspirals take roughly the same amount of time regardless of the mass ratio.The TTE inspiral takes slightly longer than the NIT inspiral, as one has to solve for twice as many equations.However once this step is complete, inspirals with the same initial conditions but any mass ratio can be computed extremely rapidly by evaluating the resulting interpolating functions.

B. Waveform mismatches
To generate waveforms, we use a semi-relativistic approximation to produce quadrupole waveforms which can be seen in Figs.6a and 6b [88].These figures show the first and last 3.5 hours of the waveforms produced by the OG, NIT and TTE trajectories respectively, sampled once every 10s.The source is viewed edge on, i.e., the detector is located at a latitude Θ = π/2 and an azimuth of Φ = 0 with respect to the source.From the plots it is clear that the two waveforms overlap significantly.The waveform mismatch between the OG waveforms and both the NIT and TTE waveforms is ∼ 2.5×10 −8 .This means that they meet the indistinguishability criteria [94][95][96] from the OG waveforms for signal-to-noise ratios (SNRs) of up to at least 4500.
Intuitively, one would expect the error in our averaging procedure to scale with the size of the orbital timescale oscillations which themselves scale with the mass ratio and the orbital inclination.As such, we wish determine the section of the parameter space where the difference in the waveforms between the OG and the NIT/TTE inspirals would be small enough not to matter for LISA data science.To do this, we first fix the mass of the primary to be 10 6 M ⊙ and create a function which uses adiabatic inspirals and root finding to numerically compute the initial orbital separation for an inspiral which will take one year to reach the ISSO for a given inclination and mass ratio.We then create a grid of mass ratios ϵ = {1, 10 −1/2 , 10 −1 , 10 −3/2 , 10 −2 , 10 −5/2 , 10 −3 , 10 −7/2 , 10 −4 } and inclinations x = {0.74,0.79, 0.84, 0.89, 0.94, 0.99}.For each point on this grid, we calculate a year-long inspiral using the OG equations of motion, the NIT equations of motion, and the TTE equations of motion and then generate a waveform from each.
The waveform mismatch between the OG and NIT or TTE waveforms is displayed in Fig. 7. From these plots we see that there does not seem to be a substantial difference in terms of accuracy from using either the NIT or TTE equations of motion.While there may be some correlation between mismatch and inclination, this does not appear to be a strong effect for θ inc ≤ 45 • .The strongest effect on the waveform mismatch comes from the mass ratio.It is worth noting that 3×10 −2 is a commonly chosen maximum mismatch for a waveform template bank that corresponds to a 90 %-ideal observed event rate [97].
Our results suggest that one could, in principle, produce such a template bank of quasi-spherical inspirals using NIT or TTE waveforms even for mass ratios as large as ϵ ≈ 0.5.In practice, we are still missing important postadiabatic contributions such as second-order effects and  contributions from the spin of the secondary; thus such a template bank would have substantial systematic biases.Also note that, formally, when all post-adiabatic corrections are included, the OG, NIT, and TTE inspirals have the same accuracy in terms of powers of the mass ratio.Consequently, we cannot say a priori which will be more faithful to Nature.

C. Using Higher Precision Fluxes
Now that we can compute fast and accurate inspirals, it is worth recalling that the relative accuracy of our interpolated GSF model is currently too low for production-level waveforms.This is primarily due to the harsh relative accuracy requirements for the adiabatic pieces of ≲ ϵ for subradian accuracy in the phases.This can be improved by incorporating information from the asymptotic fluxes, which can be interpolated to a much higher accuracy across the parameter space due to the much cheaper cost of flux calculations compared to GSF calculations.This is why our interpolated flux model is accurate to ∼ 10 −6 while our interpolated GSF model is only accurate to ∼ 10 −3 .Such information can be incorporated into the GSF model itself to improve the accuracy of the OG inspirals as well as the resulting NIT and TTE inspirals [28].However since the NIT and TTE equations of motion are naturally split into adiabatic and post-adiabatic pieces, one can calculate Γ (1) j directly from the fluxes via Eqs.( 21) and (23), interpolate them to higher precision than using a GSF model, and then substitute these improved adiabatic terms into the averaged equations of motion.
To test the difference this would make to the overall accuracy of our post-adiabatic inspirals, we looked at the final error in the φ ϕ phase when evolved from one point in the parameter space to the ISSO when using the flux+GSF model verses using only the GSF model for inspirals with ϵ = 10 −5 .From Fig. 8, we see that im-   We do not see a substantial difference in accuracy between using either NIT or TTE waveforms.We also see that the mismatch remains smaller than 0.03 for mass ratios as large ϵ ≈ 0.5.
proving the adiabatic pieces results in a phase difference that can range for from tens to tens of thousands of radians for multi-year long inspirals and this difference gets larger as one moves away from the ISSO.
As such, using high precision fluxes for the adiabatic contributions to the equations of motion is vital for obtaining accurate adiabatic and post-adiabatic inspirals.

D. Impact of the self-force on spherical inspirals
Since this is the first time that the full first-order self-force has been computed for quasi-spherical Kerr inspirals, we examine the impact that the adiabatic and post-adiabatic contributions have on the inspiral in Fig. 9.The blue curves show typical adiabatic trajectories through {r p , θ inc } space.From this we see that the self-force causes the orbital radius to decrease over time, but also cause the inclination angle to increase slightly throughout the inspiral.This is consistent with previous work on adiabatic quasi-spherical inspirals [11,13].The post-adiabatic contributions to the solutions for r p and θ inc do not change this trend in any significant way since their contribution to the orbital elements is O(ϵ).However, the post-adiabatic contribution to the orbital phases is of order O(ϵ 0 ) and is thus very significant.This is demonstrated by the dashed contours on Fig. 9, which indicate the final value of the post-adiabatic piece of the azimuthal phase φ (1) ϕ , when evolved from a given point in the parameter space to the ISSO.Since it is the azimuthal phase that most strongly impacts the final waveform phase, this gives a good estimate of how many radians one can expect the waveform to be out of phase in the absence of post-adiabatic contributions.It is worth restating that this does not include the second-order GSF contributions to φ (1) ϕ as these are not currently known.when evolved from that point in the parameter space to the ISSO.

VII. DISCUSSION
This paper presents the first calculations of the firstorder gravitational self-force in the radiation gauge for spherical orbits in Kerr spacetime by utilizing a modified version of the code found in Refs.[25,26,54].The main new element in this calculation is the inclusion of the gauge completion piece that puts the result in the correct asymptotically flat gauge.
To produce a continuous model this data is interpolated using Fourier decomposition and Chebyshev interpolation.Using only 162 points we obtain a model with sub-percent accuracy for inclinations up to 0 • ≤ θ inc ≤ 45 • .This same method also allows interpolation of the orbit averaged rate of change of energy and angular momentum from the asymptotic fluxes for all inclinations (both prograde and retrograde) with an accuracy of ∼ 10 −6 using only 342 points.This is a significant improvement over other interpolation methods found in the literature which require at least an order of magnitude more points to achieve a comparable level of accuracy [27,28].This could be further improved with a better choice when rescaling the data before interpolation, ideally choosing a function informed by the leading-order PN contributions in the weak fields and/or the analytic structure of the GSF near the ISSO.So far this model is only valid for orbits where the primary's spin is a = and so interpolating over the other values of spin is left for future work.However this work, along with paper II and Ref. [30], shows that the Chebyshev interpolation methods are a promising approach for interpolating information from expensive GSF and flux codes across the vast four-dimensional parameter space of generic Kerr inspirals.
Using our interpolated GSF model with the osculating geodesics (OG) formulation outlined in Ref. [32] and paper II along with the spherical limit derived in Appendix A, allows us to calculate the first ever quasi-spherical self-forced inspirals around a Kerr black hole.Our Mathematica implementation of this will be made publicly available on the Black Hole Perturbation Toolkit [42].However, for a binary with a small mass ratio ϵ, numerically evolving these inspirals can take minutes to hours due to the need to resolve the ∼ 1/ϵ orbital oscillations.
We overcome this by employing the technique of nearidentity averaging transformations (NITs), as outlined in papers I and II, which produce equations of motion that capture the correct long-term secular evolution of the binary but can also be rapidly solved numerically.Following the methodology of Ref. [18], we improve upon this formulation by employing a second averaging transformation such that the solutions to our equations of motion are parametrized by Boyer-Linquist coordinate time instead of Mino time.Since Boyer-linquist time can be related to the time at the detector, this improved NIT procedure is much more convenient for generating waveforms and for data analysis.We also employ a two-timescale expansion (TTE) of the NIT equations of motion which factors out the dependence of the mass ratio, at the cost of doubling the number of equations to be solved.
The quantities calculated from either the NIT or TTE equations of motion remain close to the OG evolution variables throughout the inspiral to the expected order in the mass ratio, scaling similarly as the error expected from omitted higher order post-adiabatic corrections.These averaging schemes work particularly well for quasi-spherical inspirals as we find that the mismatch between waveforms calculated from NIT or TTE inspirals and waveforms calculated from OG inspirals is ≤ 10 −3 , even for binaries with ϵ ∼ 10 −1 .
With our efficient inspiral trajectories, we investigate the effect of first-order self-force on inspirals across the spherical Kerr parameter space.We find that the orbital separation decreases with time while the orbit becomes more inclined over time, which is consistent with the findings of adiabatic evolutions of quasi-spherical Kerr inspirals.We also find that neglecting the conservative effects can result in an O(ϵ 0 ) radian dephasing in the azimuthal phase, with this effect becoming larger with the inclination angle of the orbit.
We note that both NIT and TTE equations of motion allow us to further improve the accuracy of our waveforms by replacing the adiabatic terms with higher accuracy interpolants calculated from the asymptotic energy and angular momentum fluxes.Making this improvement can result in phase differences ranging from tens to tens of thousands of radians for multi-year-long EMRIs.This highlights the necessity of efficient gravitational wave flux codes for both adiabatic and post-adiabatic EMRI wave-forms [14].
This work only incorporates first-order GSF results.Post-adiabatic waveforms will require second-order results not only to reach O(ϵ 0 ) accuracy in the phases [19] but also, as was noted in paper II, to produce gaugeinvariant waveforms [29].As such, the inspirals and waveforms produced here are only representative of the outgoing radiation gauge.
Thus, when second order effects are available for Kerr orbits, they should be folded into our framework.However, care should be taken as current second order calculations are produced assuming two-timescale expanded equations of motion and not assuming osculating geodesics and then performing the two-timescale expansion, as we have done.
We are also missing any post-adiabatic contributions due to the spin of the secondary [98][99][100][101][102][103][104][105][106].Incorporating the linear-in-spin contribution to the energy and angular momentum fluxes into our averaged equations should be straightforward.However, incorporating the conservative effects from the Mathisson-Papapetrou-Dixon equations will require careful consideration and will be the subject of future work.
We plan to extend this framework to the case of generic Kerr inspirals, but there are two challenges standing in our way.The first is the computational cost of generic Kerr GSF code coupled with the four-dimensional parameter space makes interpolating a GSF model impractical for now.The second challenge is the presence of transient self-force resonances where our NIT equations are formally singular which will force us to apply an alternate averaging procedure in their vicinity [107].There has been a lot of work covering the effects of transient resonances on EMRI trajectories due to the self-force [108][109][110][111][112][113][114] or an external third body [115,116].However, evolving through self-force resonances while incorporating all self-force effects and retaining the subradian accuracy requirement remains an open problem.
Finally, we note we have used the semi-relativistic quadrupole formula to generate the waveforms from the OG, NIT and TTE inspirals.This is sufficient for this work as we only wish to compare the difference in the waveforms caused by different inspiral calculations.However LISA data analysis will require fully relativistic waveform amplitudes such as those currently in the FastEMRIWaveforms (FEW) package for Schwarzschild inspirals [86].Currently FEW only uses adiabatic inspirals, but this can be improved by employing either our NIT or TTE equations of motion.Once the waveform amplitudes have been interpolated for Kerr inspirals, they can be combined immediately with the implementation presented in this work.

( 1 )
x by a factor of r 3 p (1 − y) and r 11/2 p (1 − x 2 ) respectively.This scaling comes from the leading order PN terms for Γ FIG. 1. Relative error in the rate of energy (a) and angular momentum (b) loss between the interpolated fluxes and the fluxes calculated with a grid of verification points that were not used for the interpolation.The white dots indicate the locations of the data points that were interpolated to produce the model and the green crosses indicate the positions of verification data.The relative error is always < 10 −5 with the exception of orbits that are very close to the ISSO.
FIG. 2. Relative error in the rate of energy (a) and angular momentum (b) loss between the asymptotic fluxes and the interpolated GSF model.The white dots indicate the locations of the data points the model is interpolating.The green crosses indicate the locations of the flux calculations and thus where the comparisons are made.The relative error only rises to ∼ 10 −2 when close to the ISSO, and is otherwise typically of the order ∼ 10 −3 − 10 −5 .

t
d⃗ q) = 0, we get the further simplification f

FIG. 4 .
FIG.4.The difference in the orbital phases for a quasispherical Kerr inspiral with M = 10 6 M⊙, ϵ = 10 −5 , and a = 0.9M when using the OG equations versus the NIT or TTE equations of motion.In both cases, the difference stays small throughout the inspiral, only becoming large as the secondary approaches the ISSO when the adiabatic assumption breaks down.
FIG.6."Semi-relativistic" quadrupole waveform strains for an EMRI with M = 10 6 M⊙, ϵ = 10 −5 , rp(0) = 7.75M , and x(0) = 0.75 as viewed edge on (i.e., with the detector at a latitude of Θ = π/2 with respect to the source's frame), sampled once every 10s.The waveform strain is normalised by the luminosity distance from the source to the detector D and the mass of the secondary µ.Panel (a) shows the first 3.5 hours, and panel (b) shows the last 3.5 hours of this year-long waveform.The blue curve is the waveform generated from the OG inspiral, the dashed yellow curve is the waveform generated from the NIT inspiral, and the purple dotted curve is the waveform generated from the TTE inspiral.
Mismatch between OG and NIT waveforms.
Mismatch between OG and TTE waveforms.

FIG. 7 .
FIG. 7. Mismatch between year long OG and NIT/TTE waveforms as a function of orbital inclination and mass ratio.We do not see a substantial difference in accuracy between using either NIT or TTE waveforms.We also see that the mismatch remains smaller than 0.03 for mass ratios as large ϵ ≈ 0.5.

FIG. 8 .
FIG.8.Absolute difference in φ ϕ for a NIT inspiral evolved from a point in the parameter space to the ISSO with a typical EMRI mass ratio of ϵ = 10 −5 , depending on whether one incorporates high precision fluxes or not.Using the high precision fluxes results in a more faithful inspiral, with phase differences with respect to the lower accuracy model in the range ∼ 10 − 10 4 radians.As such, interpolating the adiabatic contributions to the averaged equations of motion to high precision is vitally important for accurate adiabatic and post-adiabatic EMRI waveforms.