Fast Self-forced Inspirals

We present a new, fast method for computing the inspiral trajectory and gravitational waves from extreme mass-ratio inspirals that can incorporate all known (and future) self-force results. Using near-identity (averaging) transformations we formulate equations of motion that do not explicitly depend upon the orbital phases of the inspiral, making them fast to evaluate, and whose solutions track the evolving constants of motion, orbital phases and waveform phase of a full self-force inspiral to $O(\eta)$, where $\eta$ is the (small) mass ratio. As a concrete example, we implement these equations for inspirals of non-spinning (Schwarzschild) binaries. Our code computes inspiral trajectories in milliseconds which is a speed up of 2-5 orders of magnitude (depending on the mass-ratio) over previous self-force inspiral models which take minutes to hours to evaluate. Computing two-year duration waveforms using our new model we find a mismatch better than $\sim 10^{-4}$ with respect to waveforms computed using the (slower) full self-force models. The speed of our new approach is comparable with kludge models but has the added benefit of easily incorporating self-force results which will, once known, allow the waveform phase to be tracked to sub-radian accuracy over an inspiral.


Introduction
Deducing the parameters of gravitational-wave sources requires accurate theoretical waveform templates. One challenging class of systems to model are extreme mass-ratio inspirals (EMRIs), a key source for future space-based detectors such as LISA [1]. These binary systems, composed of a stellar mass compact object in orbit about a massive 10 5 -10 7 M black hole, will radiate tens to hundreds of thousands of gravitational wave cycles whilst in the millihertz band of the detector [2]. These sources, unlike the compact binaries detected with ground-based detectors [3,4], will generally not have circularized resulting in a complicated waveform with a rich morphology [5].
For a typical EMRI we expect the gravitational-wave strain induced in the detector to have a very low instantaneous signal-to-noise ratio (SNR). Instead the data can be processed using matched filtering techniques which allow the build-up of the SNR over time. This approach involves comparing the signal against expected theoretical waveform templates. Ideally the waveform templates will have two important properties. First, to avoid significant loss of SNR, they need to be accurate, ideally with a phase accuracy of a fraction of a radian over hundreds of thousands of wave cycles. Second, due to the large parameter space of possible EMRI configurations, they need to be rapid to generate, ideally on a sub-second timescale.
The 'self-force' and 'kludge' modeling approaches have arisen to meet these requirements, focusing on either accuracy or speed of computation, respectively. There is some overlap between these two methods and both are based in black hole perturbation theory, which expands the Einstein field equations in powers of the mass ratio around an analytically known black hole solution.
The primary aim of the self-force approach is to reach the sub-radian accuracy goal. Obtaining this level of accuracy requires calculating the local radiation reaction force, or 'selfforce' [6]. The equations of motion and regularization procedures employed by this method are now well understood [7][8][9][10] and many concrete calculations have been made recently [11][12][13][14][15]. Depending on the orbital configuration, computing the self-force at an instance along a worldline takes minutes to days and even if the self-force along all possible worldlines is precomputed solving the equations of motion can take minutes or hours due to the need to resolve oscillations in the inspiral trajectory on the orbital timescale.
On the other hand, the primary goal of the 'kludge' approach is a rapid speed of computation [16][17][18]. The inspiral is computed by combining fits to (orbit-averaged) numerical flux data with post-Newtonian expansions. As the equations of motion only depend on orbit-averaged quantities there is no need to resolve the inspiral on the orbital timescale. This results in a very rapid computation of the inspiral, with the tradeoff that it does not capture the physics necessary to reach the sub-radian accuracy goal. Initially developed to scope out the data analysis task these models have, over the years, been improving in accuracy by incorporating ever more physics [19,20].
In this work we develop and implement a new framework for computing EMRI waveform templates that can both easily incorporate current, and future, self-force results, and also be evaluated on a timescale comparable to kludge models. We achieve this by applying nearidentity (averaging) transformations (NITs) to the self-forced equations of motion. Before we describe this technique let us first discuss why self-force inspirals are slow to evaluate.
Within the self-force approach the secondary is treated as a point particle and the inspiral trajectory is computed by calculating the (self-)force this particle experiences due to its interaction with its own metric perturbation. At each instant the self-force is a functional of the past, inspiralling, worldline because radiation that was emitted at an earlier time can backscatter off the spacetime curvature to interact with the particle later on. This dependence on the history of the particle is what makes the self-force challenging to calculate. A self-consistent inspiral can be computed by directly coupling the equations of motion and the field equations and this has been achieved for a toy model involving a scalar charge in orbit about a Schwarzschild black hole [21] (stability issues have so far prevented a similar calculation in the gravitational case [22]). This method is very slow to compute making it infeasible for generating banks of waveform templates (but still important as the gold-standard against which faster methods can be tested). An alternative approach is to approximate the self-force at each instant by the self-force for a particle moving along the unique geodesic tangent to the worldline at that instant [23][24][25][26]. The advantage of this approach is that the self-force can be computed and interpolated across a, now finite dimensional, parameter space in a preprocessing step. Using the geodesic self-force, rather than the self-force computed using the true inspiral, reduces the equations of motion to a finite dimensional phase space and it is to these equations of motion that we apply our method. The error in the gravitational wave phase induced by using the geodesic, rather then the true, self-force contributes at O(η 0 ). Formally, this contributes to the second-order in the mass-ratio dissipative corrections, but initial calculations in the scalarfield case suggest the coefficient of this correction term is small [27].
The crucial feature of the equations of motion that makes numerically finding their solution slow is that they depend explicitly upon the orbital phase(s). As a consequence, a numerical integrator must resolve features of the inspiral on the orbital timescale. With a typical EMRI undergoing on the order of 10 5 -10 6 orbits whilst in the detector band, resolving the orbital timescale results in inspiral calculations than take minutes to hours depending upon the massratio of the binary [25]. We circumvent this problem by transforming the equations of motion to a new set of variables via a near identity transformation. This transformation has two important properties: (i) the resulting equations of motion no longer depend explicitly on the orbital phase and (ii) the transformation is small (hence 'near identity') such that the solution to the transformed equations of motion remains always close to the solution to the original equations of motion. The first of these properties allows the transformed equations of motion to be numerically solved in milliseconds, rather than minutes or hours as for the original equations. The second property ensures that the resulting solution encapsulates all the self-force physics that the original, slow to compute, solution did. The explicit form of the NIT is derived by positing a general form, with undetermined functions, for a transformation which obeys the second property (that the transformation is 'small'), substituting into the original equations of motion, and changing variables before finally choosing the undetermined functions such that they cancel the dependence on the orbital phase in the equations of motion.
Near identity transformations are not new. They have a rich history being applied to dynamical systems, and in particular planetary dynamics, stretching back more than a century [28,29]. Sometimes called near-identity averaging transformations, the effect of the transformation is to average over the short-timescale physics to produce an equation of motion which captures the long-term secular evolution of the system without the need to resolve the shorter timescale. Averaged equations of motion such as this are ideally suited to the EMRI problem where the main concern is accurately tracking the long-term evolution of the waveform phase. Near identity transformations are closely related to two-timescale expansions, which have been applied in both the PN [30][31][32][33] and self-force regimes [6,34,35]. The two methods produce equivalent results, but sometimes one is easier to use. It also seems likely that there is a close relation with the dynamical renormalization group methods of [36] when applied to an expansion in the mass-ratio.
The existence of an averaging NIT depends only on minimal conditions on the form of the equations of motion which can always be achieved when the unperturbed system is integrable. When an averaging NIT exists it is normally not unique. In section 2 we derive a general form for the NITs and discuss different choices that can be made for their form. The NIT method is applicable to inspirals in Kerr spacetime away from orbital resonances, but as a first implementation we compute generic inspirals in Schwarzschild spacetime. Our particular choice of variables, implementation details and waveform generation approach are discussed in section 3. In section 4 we discuss our numerical results, showing that the solution to the transformed equations of motion remain close to the solution to the original self-force equation of motion and that the new equations of motion can be solved orders of magnitude more quickly than the originals. We also compute quadrupole waveforms from the NIT and full self-force inspirals and show that the overlap between the two is excellent. We finish with some concluding remarks in section 5. Throughout this article we use geometric units such that the speed of light and the gravitational constant are equal to unity.

Averaged equations of motion
In this section, we derive the near identity transform needed to produce averaged equations of motion for a very generic system. Because of the general nature, this discussion is quite abstract. Readers more interested in the results could skip ahead to the summary in section 2.8.

EMRI equations of motion
We start from the self-force corrected equations of motion for an extreme mass-ratio inspiral in first-order form, where is some small parameter, which we purposefully leave unspecified. The obvious choice would be the small mass-ratio η := m 2 /m 1 , but we could also take the symmetric mass-ratio (or something else). The P = {P 1 , . . . , P jmax } is some set of 'geodesic' constants of motion (i.e. quantities that do not change along a geodesic), which together specify a zeroth-order trajectory in phase space. These could be the actions, energy, angular momentum, eccentricity, angle between the secondary spin and total angular momentum, etc. This set can also include quantities that only acquire evolutionary terms at second order such as the primary mass and spin.
The q = {q 1 , . . . , q imax } are some set of 'phases' that specify where along a zeroth-order trajectory the system currently is. Together P and q should uniquely specify a point in phase space for the system. We require these phase to satisfy two properties: (1) All functions on the RHS are 2π periodic in these phases, (2) the zeroth-order term in their evolution equation (i.e. their 'frequencies', Ω), are independent of the phases q themselves. Such a choice is guaranteed to exist if the zeroth-order system is integrable (such as the equations of motion for a test gyroscope in Kerr spacetime), in which case action-angle variables will satisfy the required property [6]. However, we stress that any other choice that satisfies the required properties will work for us.
The S = {S 1 , . . . , S kmax } are a set of quantities that are extrinsic to the EMRI's dynamics in the sense that the RHS functions in the evolution equations do not depend on them. They may or may not be extrinsic to the binary itself. The most relevant examples here are the t and φ coordinates of the secondary. Due to symmetries of the background spacetime they cannot appear explicitly in the equations of motion. Besides these it can also include truly extrinsic quantities such as the center-of-mass velocity of the binary.
Finally, the over dots represent differentiation with respect to some 'time' parameter used for the evolution of the inspiral. This could be the background Boyer-Lindquist time coordinate, proper time, or something more abstract such as Mino time [37].

Near identity transform
Our objective is to rewrite (1) in a form where the right hand side is completely independent of the phases q. For this we use a tool that has a long history in the study of dynamical systems-and planetary dynamics in particular-the near identity transform (NIT). This type of transform has previously appeared, though not necessarily by this name, in various studies of EMRIs [38][39][40][41]. The presentation here closely follows that of Kevorkian and Cole [42]. Focusing on the intrinsic quantities first (the extrinsic quantities S will be dealt with in section 2.7), the idea is to introduce a small transformation of our phase space coordinates, where we require the X to be smooth periodic functions of the phases q. Consequently, the difference between the tilded and untilded variables will always be O( ), anywhere in the phase space.
The inverse transformation is easily derived by requiring that that the composition with the original transformation is the identity and working order by order in , (3b)

Transformed equations of motion
By taking the time derivative of the NIT (2), substituting the EMRI equations of motion (1) and inverse NIT (3), and expanding in powers of we obtain the NIT transformed equations of motionsṖ withF and Here all functions on the right hand side are evaluated at P and q, and we have adopted the convention that all repeated roman indices are summed over.

Cancellation of oscillating terms at O( )
To proceed it is useful to distinguish between oscillating and average pieces of functions. For this we recall that all functions appearing on the RHS of the equations of motion are 2π periodic in all the phases. Consequently, we can decompose them into Fourier modes. If A( P, q) is such a function then we write its Fourier expansion, Based on this we can define decomposition of A in an average and an oscillatory part with Using this notation the expression for F Consequently, we can eliminate the oscillatory part of F (1) j by choosing the the oscillatory part of Y Obviously, this choice is only possible if κ · Ω = 0 for all κ such that F (1) j, κ = 0. The surfaces in orbital phase space that fail to satisfy this condition are known as orbital resonances. Evolving through these points requires a separate treatment [38,43]. For this work, we will assume no resonances occur along the inspiral. This is true generically if the primary black hole has no spin, or for equatorial or spherical inspirals into a spinning black hole. In all three cases the coefficients of the offending terms vanish by the virtue that the forcing terms depend on at most one orbital phase. However, for generic inspirals (featuring both eccentricity and inclination) into a spinning black hole resonances will have to be dealt with.

With the above choice for
Consequently, (in the absence of resonances) we can eliminate the oscillatory part of f (1) i by choosing the oscillatory part of X (1) i such that

Cancellation of oscillating terms at O( 2 )
With the choice for Y (1) j above the oscillatory part of the expression for F (2) j becomes, where we introduced the additional notation {·} to denote the oscillatory part of a product of functions. Consequently, when not at a resonance (i.e. κ · Ω = 0), we can eliminate the the oscillatory part by choosing the the oscillatory part of Y We continue in similar fashion with the oscillatory part of f (2) i . With the previous choices this reduces to, where we have left the Y (n) i , and X (1) i unexpanded if the explicit choices do not lead to a simplification. Consequently, when κ · Ω = 0 we can cancel the oscillatory part by choosing X (2) i such that,

Freedom in average pieces of transformation
With the oscillatory pieces removed by the choices in the previous sections, the remaining average parts of the tilded forcing terms become, and We have thus achieved our primary goal: effective equations of motion where the forcing terms do not depend on the phases. Beyond this there is still considerable freedom due to the unconstrained average parts of the near identity transformation X (n) and Y (n) . Various choices can significantly simplify the equations of motion. We discuss some possibilities in the following subsections.

No average terms in NIT.
The easiest choice is to simply not include any average terms in the near identity transformation, i.e. set X (n) = Y (n) = 0. Unlike some of the options below this choice is available regardless of the particular details of the original equations of motion. With this choice the NIT'd forcing functions become, (1) . The equations (26) are already quite simple, except for the expression for f (2) . We can improve on this by noting that X (1) appears in the forcing functions only through f (2) . Hence we can eliminate f (2) by solving a set of uncoupled first order PDEs, Although it may not be possible to provide an explicit solution, it is clear that solutions to these PDEs will exist. Given a numerical realization of the RHS, numerical integration of these equations should be straightforward. Combining this choice with Y (1) = Y (2) = X (2) = 0, we obtain the fairly simple expressions This is the choice we will use in the practical implementation in section 3. We further note that if one would continue the NIT to higher orders in , this choice can be made at arbitrary order to eliminate f (n) using the freedom in X (n−1) .

Elimination of post-adiabatic dissipative terms using Y (n) .
In the same spirit as the previous option we can try to eliminate F (2) using Y (1) . This again requires solving a set of first order PDEs, which are now coupled, Solutions to these equations should still exist, at the very least locally/numerically. Moreover, if one would continue the NIT to higher orders in , this choice can be made at arbitrary order to eliminate F (n) using the freedom in Y (n−1) . Together with the option of eliminating all f (n) terms with n 2 using the freedom in X (n−1) , this means that in principle (and provided there are no non-perturbative-e.g. e −α/ -terms) we can find NIT'd equations of motion that are linear in , Note that whilst the equations of motion now appear simpler, unless ∂Ωi ∂Pj = 0, the solutions for Y (1) j will appear explicitly in the expression for f (1) i . Furthermore, even if ∂Ωi ∂Pj = 0 the solutions for Y (1) j will appear explicitly in the expressions for the extrinsic parameters-see equation (42) in the section on the treatment of the extrinsic parameters.
Consequently, if there exists a left-inverse for the matrix ∂Ωi This choice eliminates all f (n) terms, yielding the following forcing functions (with some abuse of notation we write ∂Pj ∂Ωi for the left-inverse of ∂Ωi ∂Pj ) f However, the existence of a left-inverse of ∂Ωi ∂Pj is not always guaranteed. Section 3 shows a trivial way in which this can happen. Namely, if one chooses the 'time' parameter along the trajectory such that one of the Ω i is constant as a function of P, then the rank of ∂Ωi ∂Pj is smaller then i max and no left-inverse exists. Barring that particularly pathological situation, we normally have less (intrinsic) phases than 'constants of motion' (i.e. i max < j max ), because-due to symmetries of the background-some phases conjugate to the actions will be extrinsic to the local dynamics. Consequently, we should generically expect the left-inverse of ∂Ωi ∂Pj to exist. However, one may still worry that the left-inverse of ∂Ωi ∂Pj may fail to exist on local nodes in the parameter space. One particular reason to worry about this, is the occurrence of isofrequency pairs of orbits-pairs of physically distinct orbits with the same orbital frequencies Ω r , Ω θ , and Ω φ -in some regions of orbital parameter space, but not others. On the boundary between two such regions ∂Ωi ∂Pj will become singular. In [44], the existence of isofrequency orbits in Kerr spacetime was shown when the frequencies are measured w.r.t. coordinate time (no such pairings seem to exist for Mino time frequencies [38]). However, unless there are external perturbations that break axisymmetry, the φ-phase is extrinsic to the local dynamics.
So we would only need the matrix ∂Ωi ∂Pj two have rank 2 when restricted to i ∈ {r, θ} in order for the choice in this subsection to be available. We have not proven so, but this seems likely to be satisfied.
A nice aspect of obtaining the forcing functions in the form (32) is that it allows one to directly read off the successive terms in the post-adiabatic (PA) expansion of the inspiral from the F (n) terms. To evolve the orbit at adiabatic (0PA, n = 1) order, we just need the average changes of the constants of motion. At 1PA, in addition, we need the local first order self-force correction and the average changes of the constants of motion at second order. In this way, NIT reproduces the results from the two-timescale expansion of [6].

Evolution of extrinsic quantities
We now turn our attention to the evolution the quantities extrinsic to dynamics, S. Since, by definition, these quantities do not appear explicitly in the equations of motion we only need their equations of motions up to terms of order , By substituting the inverse NIT (3) and re-expanding in we can write this as an equation involving only the NIT'd variables P and q, where all functions on the RHS are now understood to be functions of P and q.
We would like to recast these equations in an 'averaged' form that is independent of the dynamic phases q. To this end we introduce a new set of transformed extrinsic coordinates S , defined by the transformation, Note that this is not a near-identity transform due to the inclusion of the Z (0) k term at zeroth order. This means that for the production of waveforms it will be necessary to know the details of this transformation.
By taking the time derivative of (35) and substituting the equations of motion for S we obtain equations of motion for S , We can eliminate the oscillatory parts of the forcing functions s for the oscillatory parts of the transformation, Z (n) k . Solutions for both equations clearly exist. Solving the first equation is akin to solving equations of motion at the test body level, which in many cases can be done analytically. The second equation would have to be solved numerically. However, in practice it is sufficient to know that it exists, since s (1) k will only explicitly appear in the second order forcing term for S (1) k . Consequently, it will only enter the waveform at order O( ), and can thus be neglected.
The remaining forcing functions depend only on P and are given by, In principle, it is possible to eliminate the first order forcing term s k completely by solving a first order linear partial differential equation for Z (0) k . However, since Z (0) k will appear explicitly in any construction of the waveform, this is of little utility. Instead it is much easier to just set Z (0) k = 0.

Summary of NIT results
Using a set of averaging transformation we have recast the the small mass-ratio expanded equations of motion for a compact binary (1) in an orbit averaged form that is independent of the phases,Ṗ The forcing functions are given bỹ whereY To recover the original variables ( P, q, S)-which are needed to construct the generated waveform-we need to apply the inverse transformation at leading order where Z (0) k is found by solving (39), preferably analytically.
This analysis independently confirms an important result from the two-timescale analysis of the same system [6]: to evolve the dynamics of the system over an O( −1 ) time making an error in the phases of at most O( ), one needs the first order corrections to the equations of motion and the average dissipative corrections at second order.
Finally, we stress an important caveat: the transformation above is only possible if the κ· Ω stay bounded everywhere along the inspiral. In other words, this procedure works only in the absence of orbital resonances. If resonances do occur, a different analysis is needed in the vicinity of the resonance [38,43].

Schwarzschild case
The previous section was purposefully very abstract so that it is applicable to, e.g. generic inspirals into a rotating black hole away from orbital resonances. In this section we apply NITs to a concrete evolution problem: the evolution of a non-spinning extreme mass-ratio inspiral under the gravitational self-force.

Equations of motion
Our first task will be to find a set of equations of motion in the form of (1). For this we employ the method of osculating geodesics [23]. At each point in time, the trajectory of the secondary is described by a tangent geodesic in the background Schwarzschild spacetime generated by the primary. To describe such a geodesic we need two constants of motion and one phase. For the two constants of motion we use the semi-latus rectum p and eccentricity e. These are defined following Darwin [45,46] using the periapsis and apoapsis distance r min and r max , p := 2r min r max (r min + r max )M , e := r max − r min r max + r min (52) where M is the mass of the massive black hole. As the phase we use the relativistic anomaly ξ (also introduced by Darwin [45,46]) defined by the relation To fully describe the trajectory of the secondary we also need two quantities extrinsic to the dynamics; the coordinate values of t and φ. The osculating geodesics evolution equations in these coordinates were provided by Pound and Poisson [23] and take the form where η is the mass ratio m 2 /m 1 and the 'time' parameter along the trajectory, χ, is defined such that when η = 0, dξ/dχ = 1. The full details of the functions f ξ , F p/e , and ω t,φ are given in appendix A. The equations (54) are of the form (1) with q = {ξ}, P = { p, e}, S = {t, φ}, and = η. We can thus follow the procedure of section 2 (using the choices of section 2.6.2) to a find an averaged version of the equations of motion, where T r (p,e) and Φ r ( p, e) are the radial period and total accumulated φ over such a period of a Schwarzschild geodesic described by ( p, e), and the averaged forcing functions are given bỹ The full details of the NIT need to achieve this form are given in appendix B, where we also give equation (56) written in terms of Fourier coefficients of the original forcing functions, a form which is particularly useful for practical implementation. In that appendix we also give the analytic formula for T r and Φ r .

Implementation
Constructing the functions, F (1/2) p/e ,f ξ/t/φ on the right-hand side of the NIT equations of motion (55) requires knowledge of both the self-force and its derivatives with respect to the ( p, e) orbital elements. At present there are no self-force codes that directly compute these derivatives. Instead, we employ an analytic model for the self-force with numerically fitted coefficients in the range 0 e 0.2 and 6 + 2e < p 12 from [24]. The analytic nature of the model makes it straightforward to take derivatives of the self-force with respect to p and e. A number of pre-processing, or offline, steps are applied to the full self-force model to construct the NIT inspiral model. These offline steps only need to be computed once. The inspiral trajectory can then be rapidly evaluated in an online step for a given mass-ratio and initial parameters ( p 0 , e 0 ). The steps that can be precomputed offline are: (i) [Offline] Compute the gravitational self-force along geodesic orbits at many thousands of points in the ( p, e) parameter space using codes such as those presented in [11][12][13]. This step can takes days running on hundreds of processors and produces gigabytes of data. Once all the data is in hand it can be interpolated using a global [24] or local [25] fit to produce the rapidly evaluated functions F p/e and f ξ . (ii) [Offline] Compute the coefficients in the Fourier expansion (8) of the functions F p/e , f ξ/t/φ on a grid of points in the ( p, e) parameter space. We choose to use a grid with regular spacing in p and e as it simplifies the construction of the two-dimensional interpolatants in step (iii). The decomposition into Fourier modes is performed using the efficient FFTW C-library [47]. With a spacing of ∆p = 0.05 and ∆e = 0.002 this step takes ∼2.5 min. This step is the first step in the calculation presented in this work as the prior step was carried out in [24] and the fit made publicly available.

(iii) [Offline] Compute the averaged forcing functions
ξ/t/φ at each point in the ( p, e) parameter space using the Fourier form of equation (56) given in equation (B.20a), and save the output to disk. This step takes less than 2 s and the stored data takes up ∼2 megabytes of disk space. (iv) [Offline] Interpolate the grid of data for each of F (1/2) p/e ,f (1) ξ/t/φ . In our implementation we use cubic spline interpolation from the GNU Scientific Library [48]. This step takes ∼35 milliseconds.
All the times quoted above are computed on a single core of a 2.5 GHz MacBook laptop. The online steps that can be computed rapidly for each set of initial conditions are: (v) [Online] Compute an inspiral using equations (55). In our implementation we solve the ODEs using an adaptive Runge-Kutta algorithm from the GNU Scientific Library [48]. (vi) [Online] With the inspiral trajectory in hand, the waveform can be computed as outlined in the next subsection.
We discuss in the results section below the computation time of the online steps. An implementation of the above steps in the C++ programming language is publicly available as part of the the Black Hole Perturbation Toolkit [49]. The code is licensed under the open-source GNU Public Licence (GPL).

Waveform generation
In most approaches the method for computing the waveform is independent of the method used to compute the inspiral trajectory. Given an inspiral trajectory there are a number of ways to construct the associated gravitational waveform. The most robust method, but also the most computational expensive, is to use the trajectory as a source in a time-domain perturbation code, such as [11] (Schwarzschild) or [50,51] (Kerr). Computing tens to hundreds of thousands of waveform cycles using this method is infeasible, but for a smaller number of cycles this approach is an important benchmark for the methods outlined below.
One alternative method is to stitch together a sequence of so-called 'snapshot' waveforms. Each snapshot is the waveform associated with a particle moving along a bound geodesic. The periodic nature of bound geodesics means these snapshots can be rapidly computed using frequency-domain perturbation codes. These snapshots can be precomputed and interpolated across the parameter space in an offline step. The waveform for a given inspiral can then be constructed by smoothly moving from one snapshot to the next. This method has been implemented in e.g. [25].
Another commonly used waveform generation algorithm is the 'semi-relativistic approximation' [52] often used by kludge methods [16][17][18]. In this approach the Schwarzschild (or Boyer-Lindquist) coordinates of the inspiral trajectory are mapped to flat-space coordinates. The waveform is then constructed using the quadrupole formula (possibly with octupolar corrections). Despite the black hole to flat space coordinate map this method has been shown to produce surprisingly accurate results in the strong-field when compared to snapshot waveforms [17].
For our purposes it does not matter which waveform generation scheme we use so long as we use the same method with the full self-force and NIT inspiral to allow for a fair compariso n. Thus, we opt to use the semi-relativistic approximation in this work as it is the simplest to implement. Details of this method can be found in, e.g. [17,53].

Results
The key feature of a NIT inspiral is that it can be evaluated rapidly and at the same time the evolving constants of motion, phases and extrinsic parameters remain within O(η) of an inspiral computed using the full self-force. In this section we present numerical results which demonstrate these two properties of the NIT inspiral. We also show that the waveform computed using the NIT inspiral trajectory is an excellent match with respect to the waveform computed using the full inspiral trajectory.
First, let us demonstrate the accuracy of the NIT inspiral method. Figure 1 gives an example of the evolution of ( p, e) and (p,ẽ) for the full self-force and NIT inspirals, respectively. We compute the full inspiral trajectory using an osculating element prescription [23] coupled to an interpolated self-force model [24]. The resulting full self-force inspiral trajectory clearly shows oscillations on the orbital timescale which, as discussed in the introduction, is what slows down the numerical computation of the trajectory. To compute the corresponding NIT inspiral we first transform the initial conditions ( p 0 , e 0 ) using the first-order NIT equation (B.1a) up to O(η) to get (p 0 ,ẽ 0 ). We then numerically solve for the NIT inspiral using equations (55). The NIT inspiral trajectory is then a smooth curve with no oscillations that runs through the 'average' of the oscillating full self-force inspiral trajectory. The accuracy of the NIT inspiral trajectory can be illustrated by applying the inverse NIT transformation, equation (3), through O(η) and comparing to the full self-force trajectory. The inset of figure 1 shows close agreement between the two inspiral trajectories. This comparison improves as the mass ratio is made smaller (we used a relatively large mass-ratio of η = 10 −3 for figure 1 to make the oscillations in the full self-force inspiral clear).
The evolution of the phase and extrinsic parameters {ξ, t, ϕ} show similarly excellent agreement between the NIT and full self-force inspirals. Figure 2 shows sample results for an inspiral with η = 10 −5 . We find the difference in the phase, |ξ − ξ|, remains less than 10 −3 over the entire inspiral excluding the last few orbits where, with the onset of the plunge, the adiabatic approximation breaks down and with it the effectiveness of the NIT. To compare NIT'd extrinsic parameters, {t,φ}, with {t, ϕ} one must first restore the O(η 0 ) oscillatory terms, Z φ } implies the NIT waveform is a good approximation to the full self-force waveform. To quantify this we compute waveforms using the kludge quadrupole approximation described in section 3.3 and calculate the mismatch between the two waveforms, minimizing over phase and time shifts. To do this we use equation (4) of [54] assuming a flat noise spectral density for the detector (in practise we compute the mismatch integral using the WaveformMatch function from the SimulationTools Mathematica package [55]). For our sample inspiral with ( p 0 , e 0 ) = (11, 0.18) we find that the waveform mismatch is always less than 5 × 10 −4 for mass ratios in the range 10 −6 η 10 −4 over durations of 2 months to 2 years-see table 1. In figure 3 we show the full self-force and NIT inspiral trajectories, waveforms and waveform mismatch for an EMRI with M = 10 6 M and η = 10 −5 .
Having demonstrated that NIT inspirals and waveforms faithfully approximate the full self-force results, we now show the rapid speed at which NIT inspirals can be computed. In order to make a fair and detailed comparison between the two methods (and other methods for computing EMRI waveforms) it is worth considering the individual steps in the calculation and their computational time. The three steps are (i) compute the phase space trajectory, (ii) Figure 1. Inspiral trajectory for a binary with η = 10 −3 and initial parameters ( p 0 , e 0 ) = (11, 0.18). This relatively large mass ratio was chosen to make the oscillations in the self-forced inspiral clear. The oscillating (blue) curve shows the trajectory of ( p, e) for the self-force inspiral. The smooth (orange) curve shows the trajectory of (p,ẽ) for the NIT inspiral. The solid (black) line shows the location of the separatrix between bound and plunging orbits. The inset figure shows a zoom in of the region inside the black rectangle. In the inset the dotted (red) curve shows the result of applying the inverse NIT, equation (3), through O(η) to the NIT trajectory. The inverse NIT trajectory and the self-force trajectory are in good agreement at this late stage of the inspiral. This agreement improves further for smaller mass ratios. compute the physical trajectory and (iii) compute the waveform. Let us examine each one in turn.
The full self-force equations of motion, (54), depend upon the orbital phase and so the numerical integrator must take many small steps in order to resolve oscillations on the orbital timescale. Consequently, computing the phase space trajectory { p, e, ξ} using the full selfforce method takes tens of seconds to hours depending on the initial conditions and the massratio (smaller mass-ratio binaries evolve more slowly and so accumulate more orbits before plunge). As the NIT equations of motion, (55), do not depend on the orbital phase they can be numerically integrated in milliseconds which, depending on the mass ratio, is 2-5 orders of magnitude faster than the full self-force method. In table 2 we give the computation time of the phase space inspirals and the speed up between the full self-force and NIT methods. The Difference in the full self-force extrinsic parameters t, φ and phase ξ and their NIT equivalents for a binary with initial parameters ( p 0 , e 0 ) = (11, 0.18) and mass ratio η = 10 −5 . All the variables oscillate on the orbital timescale which is the origin of the noisy features in the curves. Apart from close to plunge, where the NIT breaks down, we find the full self-force and NIT inspiral are in excellent agreement with |ξ − ξ| 10 −3 , Waveform mismatch between the self-forced and NIT inspirals for an inspiral with initial parameters ( p 0 , e 0 ) = (11, 0.18) and M = 10 6 M . No data is shown for the 6 months and 2 years columns for η = 10 −4 as this inspiral plunges after ∼4 months. The dominant source of error in these results comes from interpolating the inspiral trajectory when computing the waveform. The mismatch can be further reduced, at the expense of computation time, by more densely sampling the NIT inspiral trajectory or using a higher-order interpolation method. millisecond computation time of the NIT model is comparable to kludge methods, but with the benefit of including self-force corrections.
In the full self-force method the computation of the physical trajectory is normally performed simultaneously with solving for the phase space trajectory as the equations for { p, e, ξ, t, φ} form a hierarchically coupled set of equations (r is trivially computed using equation (53) and without loss of generality θ = π/2). The addition of the {t, φ} equations adds little to the computation time as their righthand side is cheap to evaluate. Consequently, the second column of table 2 is also indicative of the time to compute the phase space and physical trajectory simultaneously. Computing the physical trajectory using the NIT method is a two-step process. First {t,φ} are solved for simultaneously with the phase space variables {p,ẽ,ξ}. As with the full self-force method, this adds little computation time and the tilded variables are computed in milliseconds. To compute the physical trajectory, accurate to O(η), we need to add the oscillatory O(η 0 ) terms, Z (0) t/φ . These are given analytically in terms of elliptic integrals in equations (B.4) and (B.5) and are quick to evaluate. The total time required to compute the physical trajectory in the NIT prescription strongly depends on the sampling rate and duration of the desired waveform. For example, for a 2 month duration equally sampled at 5 s intervals (this equates to ∼10 6 samples) computing the physical trajectory takes ∼0.2 s. For kludge models the time to compute the physical inspiral depends on the model being used. Using the EMRI Kludge Suite [56] implementation of the various kludges, the Analytic Kludge [16] and Analytic Augmented Kludge [18] take around ∼0.2 s for the same duration and sample rate. The Numerical Kludge takes ∼4 s, which is longer as it directly solves the t and φ equations of motion (54d) and (54e).
Similar to the computation of the physical trajectory, the waveform computation time depends upon the duration and sample rate. Use the semi-relativistic approximation briefly described in section 3.3 with the same 2 month duration and 5 s sample rate we find the waveform takes ∼1 s to compute. This time is the same for both the full self-force and NIT inspiral and is also comparable with kludge methods.

Discussion
In this paper we leveraged the existing machinery of near-identity transformations to obtain equations of motion for small mass-ratio binary systems that are independent of the orbital timescale degrees of freedom. The result is a system of equations that can be evolved at speeds similar to previously developed kludge models, while (in principle) accounting for all physics coming from a systematic expansion of the dynamics in the small mass-ratio-e.g. gravitational self-forces or corrections due to secondary spin or higher multipoles. As a proof of principle we implemented these equations using the self-forced evolution model of [24]. The results show a speed-up of the phase space evolution of 2-5 orders of magnitude compared to evolution of the 'full' self-force dynamics, while the phase difference between the two models stays O(η). Comparing two-year duration waveforms produced from phase space evolutions from both models, we find that mismatches stay 10 −4 .
We must however stress that the implemented model should be viewed as a proof of concept. It does not provide an evolution model that is faithful up to O(η) errors in the phases. The main issue is that the model (like that of [24] on which it is based) is missing the secondorder forcing terms in the osculating geodesics evolution equations. Contributions to these functions include the second order self-force and the post-geodesic corrections to the first order self-force, neither of which have currently been calculated. Calculation of the second order self-force is currently a topic of significant effort [57][58][59][60][61][62][63][64][65][66]. Preliminary investigations of the post-geodesic corrections to the first order self-force (using scalar toy models) suggest that the contribution may be negligibly small [27]. If not, we note that we are only interested in the corrections to the orbit averaged 'fluxes'. In principle, one should be able to calculate this by comparing the time domain fluxes of adiabatic inspirals with the geodesic equivalent fluxes. Once calculations of these contributions become available (which will come at significant one-time offline computational cost), it will be trivial to include them in the NIT averaged inspiral model. Since this model already contains contributions to the second order averaged forcing functions from the oscillatory part of the first order self-force, we expect no additional online computation cost to include these effects in the inspiral calculation.
Another point to note is that our evolution model is applicable only during the inspiral phase of the binary evolution. It will breakdown as the last stable orbit is approached and adiabaticity is lost. This is not an issue for EMRIs detectable by LISA as the plunge, merger, and ringdown phases account for a very small fraction of the SNR compared to the inspiral. For binaries with more comparable mass components the loss of adiabaticity near the last stable orbit will limit the applicability of our approach though how far our model can be pushed towards comparable mass binaries remains to be quantified.
We also note that the implementation here is based on the dataset from [24], which covers only a part of the expected EMRI parameter space for non-spinning binaries. Self-force data for most of the non-spinning parameter space was published in [25]. The choice for using the older dataset of [24] was motivated by the fact that it included global analytical fits to the data. Consequently, this dataset lent it self well to calculating the phase space derivatives needed for some of the NIT averaged forcing functions. With some care it should also be possible to obtain the phase space derivatives numerically from the data of [25]. An alternative approach would be to calculate the phase space derivatives directly when the self-force is computed on a geodesic. This would involve calculating the spacetime derivatives of the self-force and some straight-forward algebra. Both are well within the technological capabilities of the stateof-the-art self-force calculations. This approach may be particularly appealing for filling the EMRI parameter space for spinning binaries. That task will be computationally much more expensive than the non-spinning case, both due to the higher dimensionality of the parameter space and the higher computational costs for calculating the self-force on generic orbits [15].
Although our proof of principle was applied to non-spinning binaries, the averaging NIT can easily accommodate the addition of spin, both to the primary and secondary. When adding spin the equations of motion will involve more than one intrinsic phase. The general derivation Table 2. Comparison of the phase space trajectory computation time between the full self-force and NIT methods for a variety of mass ratios. All the inspirals start with initial parameters ( p 0 , e 0 ) = (11, 0.18), or their NIT'd equivalent, and continue to plunge. The time to compute the full inspiral depends on the mass ratio and takes seconds to hours, owing to the need to resolve oscillations on the orbital timescale. By constrast, a NIT inspiral takes milliseconds to compute for any mass ratio. This results in a speed-up of two to five orders of magnitude, depending on the mass ratio. of section 2 has shown that generically this does not pose any issues. The exception is when the frequencies of the different phases form resonant ratios, in which case the NIT breaks down. Unfortunately, such resonances will appear generically in EMRI evolutions [38,43,[67][68][69][70]. Resonance will therefore need to be dealt with separately. The simplest approach would be to simply switch back to the full evolution equations just before hitting the resonance, evolving through the resonance, and switching back to the NIT averaged equations. This will undoubtedly work, but will come at a significant computational cost. We can already do a lot better by rather than switching back to the full equations of motion, using a NIT to eliminate all non-resonant oscillating terms as in [38]. By definition the resonant terms will only vary slowly in the vicinity of the resonance limiting the computational cost. However, the best solution would be obtained if one could implement the effects of the resonance as an instantaneous jump on the orbital parameters. The results of [43] and [38] suggest that this may be possible. In any case, however, the inclusion of resonances in fast-evolution models will require further consideration in future work. Finally, we note that in this work we employed a simple model to produce a gravitational waveform from the inspiral dynamics. For our current purposes this was sufficient, as we were using the same waveform generation scheme for both evolutionary models that were compared. However, for application in EMRI data analysis a more realistic, and fast to compute, waveform model will be needed. One approach would be to utilize the two-timescale expansion of the waveform at infinity. At leading order this will be given by a function h( P(t), q(t), S(t)) [6], which can be rewritten as function h ( P (t), q(t), S (t)). It is worth investigating whether an efficient numerical surrogate for this function can be build from the waveforms generated by particles on geodesic orbits. This will be pursued in future work.
(A.1) can be rewritten as a set of first order equations for the geodesic elements ( p, e, ξ, t, φ) [23], This can be solved analytically in terms of elliptic functions [45,46], where F(ϕ|k), E(ϕ|k), and Π(h; ϕ|k) are elliptic functions of first, second, and third kind (following the conventions for the arguments used in Mathematica), and we introduced the short-hand k r := 4e p − 6 + 2e . (B.6) Note that although the expressions for Z (0) t/φ contain explicit linear terms ξ, these are canceled by secular contributions from the elliptic functions, and as a whole the Z (0) t/φ are purely oscillatory functions of ξ.
Finally T r (p,e) and Φ r ( p, e) are the radial period and the accumulated φ over one such period or a geodesic with semi-latus rectum p and eccentricity e, The O(η) terms in the transformation are given by ∂T r ∂e dξ, (B.12) ∂Φ r ∂e F e dξ dξ, (B.13) where X (1) satisfies the first-order PDE (no explicit solution is needed anywhere), and the primitive · dξ of a purely oscillatory is chosen to be purely oscillatory. That is, given a Fourier decomposition A κ e iκξ , (B.15) its primitive is given by Finally, the second order terms in the transformation are given by To conclude this appendix we give explicit expressions for the averaged forcing functions in terms of the Fourier coefficients of the original forcing functions. Expanding the original forcing terms using equation (8), the averaged forcing functions in (56) are given by: