Transition from Inspiral to Plunge: A Complete Near-Extremal Trajectory and Associated Waveform

We extend the Ori and Thorne (OT) procedure to compute the transition from an adiabatic inspiral into a geodesic plunge for any spin, with emphasis on near-extremal ones. Our analysis revisits the validity of the approximations made in OT. In particular, we discuss possible effects coming from eccentricity and non-geodesic past-history of the orbital evolution. We find three different scaling regimes according to whether the mass ratio is much smaller, of the same order or much larger than the near extremal parameter describing how fast the primary black hole rotates. Eccentricity and non-geodesic past-history corrections are always sub-leading, indicating that the quasi-circular approximation applies throughout the transition regime. However, we show that the OT assumption that the energy and angular momentum evolve linearly with proper time must be modified in the near-extremal regime. Using our transition equations, we describe an algorithm to compute the full worldline in proper time for an extreme mass ratio inspiral (EMRI) and the resultant gravitational waveform in the high spin limit.


I. INTRODUCTION
The LIGO observation of the transient gravitational wave (GW) signal from the collision of two stellar mass black holes [1] in September 2015 spectacularly opened the new field of gravitational wave astronomy.By the end of the O2 observing run in August 2017, the LIGO/Virgo detectors had observed ten binary black hole mergers and a single binary neutron star inspiral [2].This handful of observations has already had a profound impact on our understanding of the astrophysics of compact objects and ruled out a number of modified theories of gravity [3][4][5][6][7].During the ongoing O3 observing run new events are being reported at the rate of one per week, so these constraints are rapidly improving.However, the masses of the objects being observed are all in the range of 1-100M , which is determined by the frequency sensitivity of the instruments [8].Black holes with much higher masses are expected to exist in the centres of most galaxies [9] and will be even stronger sources of GWs, but these waves will be at millihertz frequencies which are inaccessible to ground-based detectors due to the seismic noise background.
The launch of the Laser Interferometer Space Antennae (LISA) [10], scheduled for 2034, will open the millihertz band from 10 −4 -10 −1 Hz for the first time.Expected sources in this frequency band include massive black hole binaries, cosmic strings and extreme mass ratio inspirals (EMRIs).Detection of these sources, and * ollie.burke@aei.mpg.deestimation of their parameters, will rely on the comparison of accurate theoretical models of the expected gravitational waveforms to the observed data.Building these models for LISA is extremely challenging, in particular for EMRIs, which are expected to have a very rich structure and to be observed for hundreds of thousands of waveform cycles prior to merger with the central object [11].In this paper we focus on modelling of a particular class of EMRIs, in which the central black hole has very large angular momentum (spin).All of the LIGO observations to date are consistent with zero or small spin [2], but the massive black holes that will be probed by LISA are a different population.These black holes are observed in high accretion states as quasars, and accretion tends to spin the black holes up.Semianalytic models predict that the typical spins of these objects are a 0.95 [12].
The maximum spin of massive black holes is a quantity of fundamental interest for understanding the origin of black holes in the Universe.It was shown by Thorne [13] that the angular momentum of black holes being spun up through thin disc accretion saturates at a limit of a = 0.998 where an equilibrium is reached between spin up by accreted material and spin down by captured retrograde photons.Black holes with higher spin could in principle be formed directly in the early Universe and for sufficiently high mass these black holes can retain spins above the Thorne limit for a Hubble time [14].A direct observation of a system with spin above the Thorne limit would thus have profound implications for our understanding of the origin and growth of black holes.It is therefore important to understand how well observations of EMRIs can constrain the spin of near-extremal black holes and to determine this we first need to build accurate representations of the gravitational waves emitted by such systems.
The near extremal limit is also relevant for more theoretical considerations.Indeed, as the primary rotates faster, its Hawking's temperature decreases because the distance between the inner (r − ) and outer (r + ) horizons in Boyer-Lindquist (BL) coordinates reduces according to r± = 1 + 1 − a 2 = 1 + , where r = r/M and a the dimensionless Kerr spin parameter.The existence of a double pole in the function determining the black hole horizons in this limit is responsible for an enhancement of symmetry in the near horizon geometry of the Kerr black hole [15], a feature that remains true for any extremal black hole [16].This enhancement of symmetry from time translations to the conformal group has allowed several groups to analytically solve the master Teukolsky equation in the presence of the in spiraling probe particle leading to an analytic expression for the energy fluxes carried by the gravitational waves generated by this source [17][18][19][20][21][22][23][24][25].This provides a very exciting opportunity where analytic tools developed in the high energy theoretical physics community can provide accurate predictions to generate gravitational waveform templates.Future observations using such templates will be directly testing these theoretical predictions.It has already been shown that gravitational waveforms emitted by these sources contain unique qualitative features that provide a smoking gun for the existence of near-extremal systems [23].The amplitude of an EMRI waveform (averaged over a suitable amount of orbits) typically increases linearly in time for moderate spin a ≈ 0.9.It was shown in [23] that the amplitudes of these signals dampen in the high spin limit due to behaviour of the flux close to the horizon.There has been progress in modelling the inspiral from radial infinity to the innermost stable circular orbit (ISCO) [23], by integrating the geodesic equations in the near horizon geometry of the Kerr black hole [26] and exploiting the enhanced set of symmetries to compute the energy fluxes for more source trajectories [25].However, no one has focused on providing a model which encapsulates the inspiral and plunge in the limit of high spins.This is precisely what this paper seeks to do.
In this work, we build such a model for an EMRI comprised of a small compact object of mass µ gravitationally bound to a supermassive Kerr black hole of mass M and study the transition from an adiabatic inspiral into a geodesic plunge for any spin of the primary black hole.This transition to plunge was originally discussed by Ori and Thorne (OT) [27] for moderate values of the spin in the limit of small mass ratios.A similar but independent analysis conducted by Buonanno and Damour in [28] solved the problem for Schwarzschild black holes with arbitrary (reduced) mass ratio.The technical reason why high spins require a separate discussion is because of the existence of a second independent small parameter competing with the mass ratio η = µ/M 1.This new parameter is the near-extremal parameter = √ 1 − a 2 encoding the distance of the spin parameter a from its upper/lower bound, since Kerr black holes have spin parameters a ∈ [−1, +1].Since the dynamical equations describing the transition depend on the spin, the near extremal limit, i.e. → 0, modifies the original scaling discussed by OT.The transition to plunge for near-extremal EMRIs was previously considered in [29] and our work clarifies and extends those results in a number of ways.We point out the physical interpretation of the mathematical procedure used in that paper, identify a missing term in the nearextremal regime and incorporate recent analytic results for the near-extremal energy flux for the first time.
In this paper, we will first review the treatment of the transition regime given by OT in [27].We analyse their methodology and approximations and carefully estimate the scaling of terms that are being omitted.In each of [27,30,31] the notion of eccentricities and noncircular motion was ignored.We discuss the potential growth of eccentricities before and during the transition regime and find that corrections to our equations due to eccentric motion are sub-leading for any spin.We identify three separate transition regimes, each with a slightly different equation of motion: η , ∼ η and η.We then discuss a numerical algorithm to generate full inspiral trajectories in Boyer-Lindquist coordinates, alongside the corresponding evolution of the integrals of motion E(τ ) and L(τ ).Finally, we extend the waveform from the inspiral only results of [23] to include the plunge in the regime ∼ η.This paper is organised as follows.In section II, we review the properties of equatorial and circular orbits in the Kerr black hole and, in section II A, we review and compare the results describing gravitational fluxes emitted by circular EMRIs as a function of the spin.In section III we set-up the master transition equation of motion in general and, in subsection III C, we estimate corrections due to eccentricity and non-geodesic pasthistory of the orbital evolution.The transition equations of motion in the three different sxcaling regimes are described in subsections III D, III E and III F respectively.The numerical scheme to integrate our transition equations of motion for the ∼ η regime is presented in section IV A. We describe how to generate a nearextremal EMRI gravitational waveform encapsulating inspiral and plunge in subsections IV B and IV C. We finish with a summary of our main results in section V.
Notation: Any quantity carrying a tilde refers to a dimensionless quantity in units of the primary mass M, i.e. r = r/M , τ = τ /M , t = t/M , Ẽ = E/µ and L = L/M µ, but we keep a as the dimensionless Kerr spin parameter with no tilde.Dotted quantities (eg Ė) denote coordinate time derivatives of that quantity.
Finally, expressions A ∼ O(B) or, for brevity, A ∼ B stress that both A and B scale in the same way with the small parameters under consideration.We impose geometrized units by setting the constants G = c = 1.
Note added.In the final stages of this work, we became aware of overlapping results that were independently obtained in [32].

II. PRELIMINARIES
In Boyer-Lindquist (BL) coordinates (r, φ, θ, t), the motion of a point particle with mass µ in a Kerr black hole on the equatorial plane (θ = π/2) is given by [33] dr dτ where the largest root of ∆ = r2 − 2r + a 2 corresponds to the outer horizon r+ r+ = 1 + 1 − a 2 , τ = τ /M denotes proper time in units of the Kerr black hole mass M and a is the dimensionless spin parameter a ∈ [−1, 1].The particle is on a prograde (retrograde) orbit if it follows the same (opposite) direction as the rotation of the primary hole.Prograde (retrograde) orbits correspond to a > 0 (a < 0) while keeping the azimuthal component of the angular momentum L > 0. Since retrograde orbits do not reach the near horizon geometry of the primary black hole in the near-extremal limit (see appendix B), while prograde orbits do, we only consider the latter from here on.
For an equatorial orbit to be circular, the BL radial coordinate r must be constant and to be stable, the latter must be at a minimum of the potential V eff in (2) so that G = ∂G ∂ r = 0, and These conditions determine the energy Ẽ and angular momentum L of these orbits to be [34] Substituting ( 5)-( 6) into ( 3)-( 4) gives rise to whose ratio defines the angular velocity Ω of the particle Equatorial circular orbits are also known to satisfy the identity [35] where we stress the equality holds for any circular orbit labelled by (r, Ẽ, L).Differentiating (10) with respect to r, we can derive further equalities satisfied for any such orbits.The ones below will play a role in our analysis later on.The innermost stable circular orbit (ISCO) is the marginal circular stable orbit satisfying The last equality, describing marginality, allows to solve for its radius as a function of the spin [36] This is the last radii before plunging into the horizon occurs.In appendix A, we derive general formulas for (13) and higher order derivatives of G(r, Ẽ, L), which are valid for any spin a > 0, when evaluated at ISCO that will be relevant in the rest of this work.
For near-extremal Kerr black holes, it is natural to introduce the near extremal parameter to mathematically capture the large spin limit a → 1.For these black holes, the ISCO location can be expanded in and the physical parameters of this marginal orbit reduce to Notice ) for prograde orbits, whereas it is O(1) for retrograde ones, as mentioned below (B2).This further justifies our interest in prograde orbits in the near extremal limit.

A. Gravitational Wave Flux
For equatorial orbits, the motion of a (point) particle on the Kerr spacetime background generates gravitational waves carrying energy and angular momentum, either escaping towards infinity or being absorbed by the horizon of the primary hole.
Due to energy and angular momentum conservation, the orbit averaged rates of change Ė and L satisfy Ė = − ĖGW and L = − LGW , where the averaged gravitational wave dissipative fluxes are determined in terms of the time-averaged t and φ components of the dissipative self-force f diss normalised by the t−component of the four velocity u t .These quantities are discussed in more depth in the next subsection III A (see appendix C, in particular the discussion around (C7), for a derivation of relations such as (20)).Notice we split these fluxes into their horizon and asymptotic infinity contributions in the right hand side.We stress here this flux balance law only holds for adiabatically evolving binaries forcing the small mass ratio limit η → 0. The (orbit-averaged and dissipative) gravitational wave fluxes are determined by solving the Teukolsky equation in the presence of the point particle source [37][38][39][40][41]. From hereon, to avoid cumbersome notation, we shall drop the angled-brackets to denote time-averaging and simply write, for example, Ė = Ė.In [42], Finn and Thorne (F&T) parameterise the energy flux as the (Peters and Mathews [43]) leading order post newtonian correction with an extra general relativistic correction Ė factor These fluxes are spin dependent and are typically computed through numerical means.See the tables in [42] for some of the values of these relativistic corrections.
As we increase the spin of the black hole, the two roots r± = 1 ± of the function ∆ determining the outer and inner horizons of the rotating black hole coincide in the extremal limit = 0.In this limit, the geometry close to the horizon of the black hole, which can be isolated using the change of coordinates with λ → 0, has an enhancement of symmetry from R × U(1), i.e. time translations and rotational symmetry, to SL(2, R) × U(1).The resulting near horizon geometry is warped AdS 2 over a 2-sphere.The enhanced SL(2, R), the isometry group of AdS 2 , includes the scaling symmetry ρ → cρ and T → T /c.This was already observed in the original work [15] and it is true for any extremal black hole [16].Larger symmetry in physics implies larger kinematic constraints which can provide further analytic control over the given problem, in this case, the calculation of the gravitational wave fluxes (20).It is precisely the emergence of this conformal group (SL(2, R)) and the use of asymptotic expansion matching methods that allowed to find analytic expressions for these energy and angular momentum fluxes for equatorial circular orbits close to the horizon [17-20, 22, 23, 25].This body of work led to the simple relationship for the flux given in [23] The quantities CH and C∞ are constants representing how much wave emission goes towards the horizon and infinity respectively.These constants are given analytically in equations ( 76) and (77) of [20].Numerically evaluating them and summing the contribution of the first |m| ≤ l = 30 modes gives CH = 0.987 and C∞ = −0.133.(23) with exact flux data found in the BHPT.We fix the radial coordinate at r = risco and change the spin parameter a.This data can also be found in [20].
The flux (23) is only reliable when working with near extremal black holes and interested in near horizon physics.This fact can be checked by comparing the exact fluxes (21), using exact results found in the black hole perturbation toolkit (BHPT), with the near extremal approximation (23).This comparison is shown in figure 1. Fixing the radial coordinate to r = risco and varying the spin parameter a, we observe in Table (I) that as a → 1, the NHEK flux (23) converges towards the exact value computed using the BHPT.Furthermore, fixing the spin parameter to a = 1 − 10 −9 , as in figure (1), the NHEK flux (23) provides a nearly-perfect agreement up to a coordinate radii r ≈ 1.012.The reason for the (extremely small) discrepancy at the ISCO is because Eq.( 23) is only valid for → 0 and we consider ≈ 10 −5 .Thus we can use (21) to build a trajectory throughout the adiabatic inspiral regime.Then, as we near the ISCO, we can use the powerful analytic result given by Eq. (23).Using Eq.( 23) allows for a more analytic treatment of the analysis of the transition regime.

III. THE TRANSITION EQUATION OF MOTION
In this section we revisit the earlier work by OT [27] describing how a small body following an initial equatorial circular orbit around the large black hole inspirals and eventually transitions into a plunging trajectory falling into the black hole.
Our discussion is organised as follows.First, we analyse in section III A the effects arising from the radial self-force in the vicinity of the ISCO on the dynamics of this small body, justifying the starting point in OT.Second, assuming the dissipative fluxes of energy and angular momentum for quasi-circular and equatorial orbits are still related as in circular orbits [35,44] we derive in section III B the transition equation for arbitrary black hole spins without the OT assumption that both energy Ẽ and angular momenta L evolve linearly in proper time τ .Third, given the quasi-circular nature of our assumed orbits, we argue in section III C there can be corrections to (24) of the form whose scaling behaviour on the trajectory of the small body is determined.Finally, in sections III D-III F, we discuss in great detail the existence of three different scaling regimes in our transition equation, depending on the black hole spin, paying special attention to the nearextremal ones which contain new physics.We show the corrections due to (25) are subleading in all the regimes.

A. The self-force
This subsection shows both that quasi-circular and equatorial orbits have vanishing dissipative effects and the conservative piece of the radial self-force can be neglected close to the ISCO.
Consider the radial geodesic equation, Eq. ( 2).Differentiating it with respect to proper time, one obtains It is shown in Appendix C this is equivalent to Hence, the terms on the left hand side of eq. ( 26) correspond to the usual ones for geodesic motion, whereas the one on the right hand side can be understood as the perturbing force f r diss exerted on the particle driving energy d Ẽ/dτ and angular momentum d L/dτ loss due to gravitational wave emission.
In [45], Mino recognised that the forcing term for general geodesic motion can perturbatively be split into a radiative reactive dissipative and a conservative piece at first order in the mass ratio η More details on this splitting can be found in [46].For circular orbits, as considered by OT, the dissipative fluxes of energy and angular momentum are related by [35,44] Figure 1: These plots show the deviation between using the exact results for the flux (21) and the near extremal approximation given in (23).Notice that, to keep the error < 5%, we require r 1.01.For each of these plots, we used a spin parameter a = 1 − 10 −9 .Similar plots can be found in [20].
Hence, the dissipative part of the self-force f r (1)diss vanishes, leading to which is precisely Eq.(3.10) of OT in [27].A gauge invariant way to quantify the leading order in η effect on the trajectory due to f r (1)cons is to study the orbital velocity Ω (see [46] for a review).This generates a shifted orbital velocity Ωshifted isco with respect to the Kerr orbital velocity Ωisco at ISCO given by [47,48] ( with the quantity C Ω(a) discussed in depth and independently (numerically) calculated in both [48,49].According to [48], C Ω(a) ∈ (1.24, 1.39).Hence, C Ω(a) ∼ O(1) for all spins, and since Ωisco ∼ O(1), it follows that for η 1 It is further shown in [48] that C Ω(a) → 1 + 1/2 √ 3 as a → 1 in an averaged sense1 .Using Eq.( 19) in the high spin limit and η 1, equation (31) becomes This implies that the change in the orbital velocity at the ISCO due to conservative self-force effects is an O(η) quantity.
Since Ωshifted isco is related to the shifted Boyer-Lindquist radial coordinate at the ISCO by it follows, using eq.( 33), that the "radial thickness" differs by an O(η) quantity when including the conservative self-force effects.
It will be shown in this paper that there are three different transition regimes depending on the ratio of and η.The "radial thickness" of the transition in each regime scales according to Thus, the effect of the conservative piece of the self-force is subleading in all regimes.For this reason, like in the original OT analysis, we shall ignore these effects.

B. Transition Equation -Generalities
To discuss the evolution of the orbit, we pursue the following strategy : we assume the corrections in (25) are subleading, and once the scaling behaviour of the different dynamical regimes is identified, we double check the consistency of our original assumption.
To evolve the orbit, OT used the circular flux relationship (29) and additionally assumed that the energy Ẽ and angular momentum L evolve linearly in proper time τ throughout the transition regime In our analysis of the transition, we will not assume a strict equality in Eq. (36).Instead, we will keep track of the evolution of Ω−1 isco ( Ẽ − Ẽisco ) − ( L − Lisco ), as also considered in [29].
OT proposed to analyse the transition to the plunging geodesic by expanding (26) around the ISCO trajectory (r isco , Ẽisco , Lisco ), since the latter provides the natural starting point for the plunging trajectory for equatorial and circular orbits.It is physically natural to introduce the new variables δE, δL and to study the inspiral evolution of the small body perturbatively around the primary.The presence of Ωisco is for technical convenience.
Instead of expanding ( 26), we find it more convenient to expand (2).Our conclusions do not depend on this choice.The latter is given by Since G(r, Ẽ, L) is quadratic in Ẽ and L, we have ignored the derivatives Plugging the perturbative variables (37), using the definition of the coefficients (A8) and the results in (A9)-(A11), one can rewrite the general transition equation as where Γ is defined by Notice that Γ ∝ δE − δL at leading order in R. Hence, it encodes the deviations from the OT approximation (36).
The time evolution of δE − δL near risco is controlled by the fluxes and the angular velocity.Throughout a quasi-circular inspiral far from ISCO, the compact object inspirals on a sequence of circular geodesics defined by the constants of motion Ẽ(r circ ) = Ẽcirc and L(r circ ) = Lcirc , as given in Eq (5) and Eq.( 6) respectively.The evolution of the constants of motion is linked through Eq. ( 29) above, which simply states that circular geodesics evolve into circular geodesics.It can be shown that solutions to the Teukolsky equation for circular orbits obey this condition [35,50].For circular evolutions we therefore see that where we expanded Ω(r) to first order in R and approximated d L/dτ ≈ (d L/dτ ) isco = −ηκ for κ constant defined by Thus we deduce that δE − δL ∼ ηRτ for circular inspirals close to risco .We shall see that these corrections are indeed subleading in the regime considered by OT [27].However, they will not be negligible for near-extremal black holes.

C. Corrections arising from deviations from adiabatic nearly-circular inspiral
Given our assumption that the orbit is nearly circular when it reaches the transition regime, one expects corrections to the relation (29) between the fluxes of energy and angular momentum satisfied for an exactly equatorial circular adiabatic inspiral.We discuss below two possible physical effects giving rise to such corrections : eccentricity and the non-geodesic past-history of the orbital evolution.These will give rise to the corrections (25).
Eccentricity can lead to corrections to the transition equation which we will discuss further below, but eccentricity corrections to the fluxes tend to be suppressed during the transition regime.This is because the transition, for an arbitrary eccentric inspiral, corresponds to the orbit passing over the maximum of the effective potential given by Eq.( 2).The radial velocity throughout the transition regime is therefore always small, while the angular velocity remains O(1).Hence the orbit looks very much like a circular orbit, even if it is technically eccentric or even plunging.For nearly-circular transitions, the orbit is passing over a point of inflection of the effective potential and corrections to this approximately-circular assumption are even smaller.
Corrections from non-geodesic past-history enter because the self-force acting on the small object at a particular time is generated by the intersection of the particle world line with gravitational perturbations generated by the orbital motion in the immediate past [51].The self-force acting on the orbit when it is at a particular radius will therefore have corrections that depend on how far, in radius, the orbit has moved over the relevant past-history.The latter is determined by the dominant, azimuthal, timescale, and is an O(1) quantity, when expressed in coordinate time 2 .The orbital radius therefore changes by an amount of O( ṙ) over the relevant past-history.This is the scaling of the fractional change in the fluxes, and since Ė ∼ O(η be considered O(1) and so the scaling of this correction is η ṙ.This is the first type of correction in Eq.( 25).In the adiabatic inspiral phase, these corrections are O(η 2 ) and form part of the second-order component of the selfforce.However, in the transition phase these corrections can be larger.
We have argued above that eccentricity corrections to the fluxes should be suppressed in the transition regime.We now make this more concrete.Eccentricity corrections to the fluxes enter as fractional corrections of O(e 2 ), since corrections to the orbit at linear order in eccentricity are oscillatory and average to zero over a complete orbit [35].The corrections to the coordinate time fluxes thus scale like η 2 3 e 2 (which is ηe 2 in the OT regime discussed in Section III D).This is the second type of correction in Eq. (25).If these corrections are to be small relative to the non-geodesic past-history corrections, we need e 2 < ṙ.In the transition zone we will see that the proper time scales like R −1/2 , where R = r − risco is the distance from the ISCO, regardless of the spin of the primary.For non near-extremal black holes, i.e., those with η , proper time and coordinate time scale in the same way and the scaling of ṙ is therefore the same as that of R 3/2 .The constraint we obtain on eccentricity is therefore e < R 3/4 .However, there is also a geometric constraint, which is that the variation in the orbital radius due to eccentricity should be small compared to the variation due to radiation reaction through the transition zone.The latter is the scaling of R, while the former is a quantity of O(e), so we deduce an additional constraint e < R < R 3/4 , the latter inequality following from the fact that R is a small quantity throughout the transition.We deduce that the geometrical constraint is stronger than the flux-correction constraint in the regime η .In the near-extremal regime, η , d t/dτ ∼ −2/3 and so the constraint on the eccentricity changes to e < 1/3 R 3/4 if these corrections are to be subleading.This is then more stringent than the geometric constraint.However, in this regime we will see below that eccentricity cannot grow until deep inside the transition zone, so even the more stringent constraint is easily satisfied.
Eccentricity during the transition can arise either from the presence of residual eccentricity prior to the start of the transition zone, or due to the excitation of eccentricity during the transition.The latter manifests itself as additional terms in the transition equation, the existence of which we will check for carefully in our analysis.To understand the former, we need to analyse the growth of eccentricity during the adiabatic inspiral.We will assume that at the beginning of the inspiral the orbit is nearly circular.It was shown in [35] that, for small eccentricity, the evolution of eccentricity under radiation reaction takes the form ė = f (r 0 )e, where r0 is the mean orbital radius and e is an eccentricity defined such that the orbital apoapsis is at r = r0 (1 + e).For large r0 , f (r 0 ) < 0 and so the eccentricity decreases.In this regime any small eccentricity that is excited by small perturbations arising due to inspiral evolution or other effects is damped away and does not grow.However, for all spins a < 1, as the innermost stable circular orbit (or separatrix) is approached the sign of f (r 0 ) changes and is greater than zero in the vicinity of the ISCO.This means that orbits near to the separatrix are unstable to eccentricity growth.We would therefore expect any eccentricity that is excited to begin to grow.
Denoting ṽ2 = 1/r 0 , Kennefick [35] showed that the evolution of the orbital parameters, for small eccentricity, was governed by equations of the form where Both Γ and Ė0 are components of the self-force, which can be evaluated by solving the Teukolsky equation.The quantity Γ is in fact a linear combination of quantities that are time derivatives and so the above equation takes the same form for any choice of time coordinate with respect to which to evaluate the fluxes.Kennefick's analysis used coordinate time and so we make the same choice in the following discussion.An explicit expression for Γ is given in [35] and the quantity Ė0 is the energy flux given in (20).Numerical calculations show that these are finite quantities of O(1) throughout parameter space.The first term in the eccentricity evolution equation vanishes for evolution driven by gravitational radiation reaction, while the quantity h(ṽ) is singular at the ISCO.Therefore, close to ISCO the eccentricity evolution takes the form ė e ≈ j(ṽ)h(ṽ) Ė0 Notice that the expression is entirely geodesic and independent of the energy flux Ė0 .For non-extremal spin, both j(ṽ) and h(ṽ) have simple poles at r = risco and there is a simple zero in the term (1−6ṽ 2 +8aṽ 3 −3a 2 ṽ4 ) in the numerator.Therefore as ISCO is approached the eccentricity evolves as with R = r − risco as before, and e 0 denotes the eccentricity when R = R 0 and r0 risco .The exponent k(a) is given by k(a) = H(ṽ isco )/D(ṽ isco ) where D(ṽ) = 2ṽ 2 (1 − 2ṽ 2 + aṽ 4 ) and ṽ2 isco = 1/r isco .We find that k(a) = 1/4 for all a < 1.The behaviour for near-extremal black holes is slightly different, which we will discuss further below.
For extremal black holes the various factors in the expression for d ln e/dr 0 have repeated roots at the ISCO.To understand the behaviour for near-extremal black holes we therefore need to do an expansion in both R and .This takes the form The terms omitted from both the numerator and denominator above are O(1) in .The ratio a 0 /b 1 = −1/4, agreeing with the result for k(a) found above.However, for R, the behaviour is not dominated by this term, but by the terms from a 6 in the numerator and from b 7 in the denominator.The leading order behaviour in this regime is therefore This is also exponential, but we find the ratio a 6 /b 7 = 3/2, i.e., it is greater than zero and therefore the eccentricity decreases exponentially until we reach the regime R ∼ .This is the statement that the critical curve, where the sign of the eccentricity evolution changes, is in the near-horizon region, which is consistent with results in [40].We conclude that for near-extremal black holes, eccentricity can only grow once the inspiraling object is already very close to the ISCO, which is typically already inside the transition zone.
To complete this discussion we need to determine the scaling of the initial eccentricity e 0 .If the orbit is truly circular then the eccentricity remains zero, so there must be some mechanism to excite an initial eccentricity which can then grow.Eccentricity can be excited by other physical processes, such as the presence of perturbing material, e.g., dust, or gravitational interactions with third bodies.Those processes are important, but in the pure-vacuum case eccentricity could still in principle be excited by the evolution under radiation reaction.We argued earlier that corrections to the fluxes far from the horizon scale like η ṙ which is η 2 during the adiabatic inspiral.These corrections mean that the first term in Eq. ( 45) is no longer exactly zero.Setting that term to η 2 we find an evolution equation of the form de 2 /d t ∼ η 2 .After a few orbits the eccentricity is then O(η).This eccentricity induced by second order corrections to the evolution is damped by the process described above, until we reach the critical curve where it grows, eventually exponentially near the ISCO.This suggests appropriate initial conditions are e 0 ∼ η and R 0 ∼ O(1). 3 We note that this mechanism could also excite eccentricity during the transition zone itself, but this would be of order e 2 ∼ η ṙ 2 3 and hence no larger than the non-geodesic past-history corrections described above.If eccentricity grew coherently throughout the transition zone, the eccentricity induced by this process would be no larger than e 2 ∼ η ṙ 2 3 T , where T is the coordinate time elapsed through the transition zone, which is typically smaller than the eccentricity grown during adiabatic inspiral prior to the start of the transition zone.
To summarise, we expect corrections to the evolution equations that arise from higher-order terms in the flux to scale like η 2 3 ṙ (which is η ṙ in the OT regime discussed in Section III D), and we expect residual eccentricity in the transition zone to be no more than e ∼ ηR −k(a) .In the non-near extreme case, these eccentricity corrections will be important when e > R, which implies R < η 1/(1+k(a)) .In the near-extreme case, the corrections only become important when R ∼ , so we simply need to check that this is well inside the transition zone.In the analysis that follows we will evaluate 3 A natural continuation of this argument would be to say that the second-order self-force induced corrections continue to drive eccentricity growth, over the whole of the inspiral, lasting a coordinate time ∼ η −1 , leading to a final eccentricity of O(η 1/2 ), which can be larger than the eccentricity grown through the mechanism discussed here.However, this assumes that the eccentricity grows coherently and monotonically.In practice, once the eccentricity is O(η), the radial motion due to eccentricity becomes larger than the amount the radius evolves over the relevant past-history that determines the self-force and so the argument that the latter is the dominant contribution to corrections no longer applies.Knowledge of the second-order selfforce would be required to fully explore the further evolution of the eccentricity and this is not currently available.However, we expect that the growth of initial eccentricities of O(η) through the instability mechanism will be the dominant contributor to the residual eccentricity in the transition zone.
the scaling of these terms and show that they are subdominant for inspirals into near-extremal black holes.

D. Ori and Thorne regime
Consider non-extremal black holes, i.e. rotating black holes where the extremality parameter is not close to zero so that η . In this regime of spins and according to the discussion below (A13)-(A18), all the coefficients controlling the general transition equation ( 40) and (41) are O(1).This is the regime originally discussed in [27].
Omitting coefficients of order one, the dominant contributions to the transition equation are where we also omitted any further terms from ( 40) and (41) since they are subleading.Looking for a scaling solution R ∼ η p and τ ∼ η q , it follows, using equation ( 36) that δL ∼ η 1+q .Requiring all dominant terms to have the same scaling fixes p = 2/5 and q = −1/5, so that Notice the overall scaling of the transition equation is (dr/dτ ) 2 ∼ η 6/5 .The remaining question is whether the dominant term in Γ ∼ δE − δL is subleading or not.From (42), it follows δE −δL ∼ η 6/5 in this regime, suggesting the change of variables This allows to write the schematic transition equation as Terms in Eqs. ( 40) and ( 41) that have been dropped can be seen to scale like the above terms multiplied by additional powers of R or δL.Since both R and δL are small quantities in the transition zone, these terms are sub-leading and we can ignore them.
The above scaling analysis proves the dominant terms in (40) in the regime η are captured by where we neglected all subleading corrections, kept the same original notation as in OT [27] for the coefficients and the dominant contribution to (41) Keeping all coefficients of order one, the natural scaled variables to introduce are where Plugging this into (55), one obtains where we defined C 0 = α 4/5 (κβ) −6/5 .From now on, we ignore the subleading corrections.
The analogue of the acceleration equation ( 26) reduces to This depends on the time evolution of the circularity deviation parameter Y , whose dominant contribution is derived in (42).Inserting the re-scaled variables (59) into Eq.(42) leads to a transition equation Evaluating (11) at ISCO, we find the term in square brackets equals and so the last term vanishes.This was inevitable, since this term is precisely the term that arises from the dissipative part of f r from Eq.( 28), as identified earlier.
The leading order evolution of Y is driven by maintaining the circularity of the orbit and so with this condition we expect the radial self-force corrections to be subleading.
The resulting transition equation of motion in the regime of low spins η is and Y is evolved through the ODE (62).We note that the transition equation does not depend on Y in this regime.Corrections to this equation arising from evolution of Y enter at an order η 2/5 higher than leading and so are subdominant.As discussed earlier the evolution of Y is related to deviations from the linear-in-propertime evolution of energy and angular momentum and so the fact that these corrections do not enter the transition equation for η demonstrate that the linear evolution assumed by OT is appropriate in this regime.
Let us check the self-consistency of the transition equation ( 64) by verifying that all neglected corrections to it are indeed smaller when evaluated on the scaling regime (59).First, as discussed in section III A, the corrections to the orbit due to the conservative piece in the self-force are order O(η), see (35).This is indeed smaller than the "radial thickness" R ∼ η 2/5 in (59).Second, corrections due to η ṙ, appearing in (25), are O(η 8/5 ).Hence, these corrections are O(η 3/5 ) smaller than the dominant δL and δE scaling in (59) 4  Finally, corrections arising from eccentricity are subleading provided e < r − risco , as discussed in Section III C. In the non-extremal case we therefore need e < η 2 5 , due to (47).This yields the constraint We saw previously that k = 1/4 for all spins a < 1, which satisfies this bound.We deduce that eccentricity corrections are subdominant in the non-near-extremal regime. 4Using the more conservative assumption that the averaging timescale is determined by the period of radial oscillations, which scales with T ∼ η − 1 5 , the corrections are still suppressed by a factor of η 2 5 .Third, corrections to dΓ /dτ arising from non-geodesic past history corrections to the fluxes scale like η 8/5 and those arising from additional terms in the expansion of the azimuthal frequency as a function of radius scale as η 9/5 , which are both subdominant to the leading η 7/5 scaling, albeit only by a factor of η 1/5 .

E. General Transition Equation of Motion -Near-Extremal
Let us consider rapidly rotating black holes with spin parametrized by a = √ 1 − 2 for 1, as in (15).The discussion below equations (A13)-(A18) allows to identify the a priori dominant contributions to the transition equation (40) as (65) Since the functional dependence of the above equation does not depend on η, we learn the η scaling should be the same as before if we keep the R 3 and R δL terms.Hence, we are left to determine any possible scaling.Proceeding as before, we look for scalings of the form R ∼ η 2/5 p and τ ∼ η −1/5 q .We learn from equation (36) that δL ∼ η 4/5 q .Requiring these dominant terms to scale in the same way determines p = 4/15 and q = −2/15, so that Hence, if η ∼ , the term R 2 δL scales like the velocity squared (dr/dτ ) 2 ∼ η 6/5 4/5 ∼ 2 and must be kept in the transition equation, whereas the term δL 2 4/3 is O( 2/3 ) smaller and, consequently, subdominant.

Introducing the finite variable
the general transition equation in the η ∼ regime reduces to dR dT Notice the radial velocity throughout the transition regime scales like dr/dτ ∼ η 3/5 2/5 ∼ η in the regime ∼ η.This is as in the adiabatic regime, but smaller than in the OT regime where dr/dτ ∼ η 3/5 .
As a self-consistency check, we can write the radial geodesic equation using the change of variables (66) and (67) It is apparent that the dominant terms are the i = 3 and m = 2 terms, all others being subleading.

(73)
Since η ∼ , it follows R ∼ 2/3 .Hence, the near ISCO expansion corresponds to the near horizon geometry of the primary black hole since, in Boyer-Lindquist coordinates, |r isco − r+ | ∼ 2/3 .As a result, we will be able to use the (leading order and analytic) expression for the energy flux due to gravitational radiation in (23).This allows to compute κ in (60) in this regime as Ignoring subleading terms, the general transition equation (70) reduces to with C 1 = γ(α βκ) −3/5 κ (76) Γ = −4/5 η −6/5 α 4/5 ( βκ) −6/5 Γ . (77) Notice the appearance of the new term proportional to T X 2 , compared to the OT regime, is due to the regime η ∼ .Taking a further T derivative, we find the analogue of the acceleration equation (26) in this regime The time evolution of Γ in (72) has two contributions : one proportional to dY /dT , which can be computed using ( 42) and a second one proportional to Y (dX/dT ).Altogether yields where we used (11) to simplify the first term.The latter cancels the −2X term in (78).Using the dominant contribution to the identity (12) evaluated at ISCO, the second term cancels the C 1 X 2 term in (78).Finally, the third term gives a non-trivial contribution to the acceleration equation with constant defined by and evolution equation for Y such that (82) In our treatment of the OT regime (non near-extremal spins), the terms in Eq. (80) were neglected since they scaled with η 2/5 and were subdominant.In the nearextremal case, one can clearly see that the XT and Y term are comparable to the (rescaled) radial acceleration provided η ∼ .As such, they must be included in the analysis.Our final transition equation of motion differs from Eq.( 43) in [29], which correctly included the Y term but missed the cross term XT , which is the same order as the terms being retained.Our analysis improves on [29] in two additional ways.Firstly, Y was introduced in [29] as a mathematical construct to ensure conservation of the four-velocity norm.The evolution equation for Y was derived by forcing the equation of motion obtained from differentiation of the kinetic energy equation, Eq. ( 2), to agree with that obtained by expansion of the left-hand-side of the acceleration equation, Eq. ( 26).This is equivalent to setting the radial self-force term to zero, which is equivalent to imposing the circular-to-circular condition.This physical interpretation of the procedure was not made clear in [29], nor the interpretation of Y as representing departures from the linear-in-proper-time evolution.Secondly, the scaling of the flux given in Eq. ( 23) was not known at that time and this was left as an unspecified power of .Now that we know this scaling we can do a more complete analysis of the near-extremal regime.
The quantities above can be computed in the nearextremal limit, → 0, Equations ( 80) and ( 82) are a coupled set of ODEs which will link the adiabatic inspiral to a plunging geodesic.As in the previous section we now consider the size of corrections to the transition equation.Corrections to the circular flux-balance law in the geodesic part of the transition equation scale like η ṙ 2 3 according to (25).These are O( ) smaller than the terms kept in this regime 5 .Similarly, corrections to the linear-in-propertime angular momentum evolution enter through corrections to δE and δL and scale like dr/dτ times terms that are being retained.These are therefore subdominant since dr/dτ ∼ η 3/5 2/5 1.These corrections also contribute additional terms through corrections to the radial self-force part of the transition equation.These are of order η • ∂G/∂ Ẽ and η • ∂G/∂ L, which scale like η 2/3 and so are a factor of (η/ ) 1/5 1/3 smaller than the leading order terms in the transition equation and are therefore sub-dominant.
Eccentricity corrections enter like fractional e 2 corrections to the fluxes, and are only more important than the corrections described above if e > R ∼ η 2/5 or e 2 > ṙ.In the near-extremal regime ṙ ∼ 2/3 dr/dτ ∼ η 5/3 and so eccentricity corrections become important when e > η 5/6 .However, as shown in Section III C, for near-extremal inspirals eccentricity can only grow once 5 Using the conservative assumption about the averaging timescale, these corrections are sub-leading by a factor of η r − risco ∼ O( ).In the transition zone r − risco ∼ (η/ ) 2/5 2/3 and so eccentricity has not started to grow when the transition zone is reached.Residual eccentricity from the adiabatic inspiral would be O(η) and eccentricity excited during the transition would be O(η 4/5 8/15 ) (or O(η 7/10 7/15 ) if it was coherently excited throughout the transition).These are smaller than the threshold η 5/6 at which the eccentricity corrections become more important than the non-geodesic past history corrections, which we have already shown to be sub-leading.

F. General Transition Equation -Very
Near-Extremal The final regime concerns very rapidly rotating black holes, where η.Using the results in appendix A, one can identify the a priori dominant contributions to the transition equation ( 40) and ( 41) to be (ignoring coefficients of O( 1)) (84) It is natural to expect that terms involving some explicit factors of should be sub-leading in this regime.Assuming a scaling solution of the form R ∼ η α and τ ∼ η β , we learn using (36) that δL ∼ η β+1 .Imposing the dominant terms R 3 and R 2 δL scale like (dR/dτ ) 2 yields the scaling solutions α = 2/3 and β = −1/3, so that As a consistency check, notice the terms 2/3 RδL ∼ η 2 ( /η) 2/3 and 4/3 δL 2 ∼ η 8/3 ( /η) Notice the radial velocity throughout the transition regime scales as dr/dτ ∼ η, as it were throughout the adiabatic inspiral regime and in the near-extremal case [see sec.(IIIE)].Thus the radial motion is fastest throughout the transition regime when the primary is of moderate spin: η .As a consistency check, we can substitute the scalings (85) and (86) into the general transition equation ( 40) Clearly the dominant terms occur when both i = 3 and m = 2 with the rest being subleading.
The above scaling analysis proves the dominant terms in (40) in the regime η are captured by where α and γ are given in Eq.( 71) with Keeping all coefficients of order one, the natural rescaled variables in this regime are In these variables, the radial velocity equation (89) can be expressed as with and κ as in (60).Taking a further derivative with respect to T yields the acceleration equation Using (42) together with (91), one finds that (95) Plugging this back in (94) and using the dominant contribution to the identity (12), the K 1 X 2 term cancels and one is left with together with the evolution equation for Y (T ) given by where In the limit → 0, the constants K 1 and K 2 approach the values As argued in previous sections, corrections to the circular flux-balance law contribute terms to the transition equation which scale like dr/dτ ∼ O(η) times terms that are being retained and like η 2/3 .Corrections to the linear-in-time angular momentum evolution enter with the same scaling as the former.The retained terms in the transition equation scale like η 4/3 in the very near extremal regime and so these corrections are both subleading.Eccentricity corrections enter like fractional e 2 corrections to the fluxes, but, as in the near-extremal case, eccentricity cannot grow until the transition zone has already been reached, and so these corrections are no larger than O(η 5/3 2/3 ) and are also sub-leading.We conclude this subsection by noting that the transition equation of motion (96) is perfectly well behaved in the limit → 0 and can therefore be used to compute an inspiral into a maximally spinning black hole with a = 1.In this case, the horizon coincides with the ISCO in BL coordinates.However, the proper distance is ∆ = − M 3 ln (see Fig. 2 in [52] together with explanations in [52,53] and more recently in appendix A of [54]).Hence, we terminate the integration of the ODE (96) at r = r+ , since our numerics are specific to BL coordinates.The presence of the horizon manifests itself in the transformation from proper time to coordinate time, which will be discussed, for non-extremal inspirals, in the next sub-section.

A. Numerical Integration
We now seek to compute a full worldline r(τ ) for ∞ > r ≥ r+ .Out of the three regimes just discussed, we restrict ourselves to the ∼ η one.This is because the η regime has already been considered in the literature [27,[29][30][31] and the η regime has been argued to be inaccessible throughout the transition regime in [32].The latter conclusion follows from the observation that the waves emanating from the secondary produce a spin down effect on the primary leading to a maximum attainable spin with ∼ η.For ∼ η, we try to find the solution to which deviates off the past adiabatic inspiral and evolves into a geodesic plunge.The constants in Eq.( 98) are given by Eq.( 83).We can derive an equation for an adiabatic inspiral in proper time by using the quasi-circular approximation.Using our far-horizon expression for the energy flux defined by Eq.( 21) with both equations ( 8) and ( 5), one derives This equation diverges at the ISCO which is a break down of the quasi-circular approximation.We shall use Eq.(98) to smoothly transition from the adiabatic inspiral Eq.(99) into a geodesic plunge to the horizon.We used a cubic spline to interpolate values for the relativistic correction Ė(r) using exact flux data found in the BHPT.We then numerically integrate Eq.(99) by stepping forwards in proper time until L(τ ) − Lisco ∼ η 4/5 −2/15 .We feel this criteria is suitable for turning on the transition equation of motion since our model for the flux is well represented during the transition regime.When this criteria is met we can be sure that our model for flux evolution throughout the transition regime is correct to leading order.Once this is satisfied, we stop integrating our adiabatic inspiral solution and begin integrating our transition equation of motion (98).Since we do not terminate our adiabatic inspiral solution at the ISCO, we do not know the precise proper time where the particle crosses the ISCO.As such, the variable T is not a good choice of variable to integrate on the right hand side of (98).Instead, we substitute T for δL from Eq.(73) into our transition equation of motion, then (100) We use initial conditions determined by the end of the adiabatic inspiral Eq.(99) at some time τinit .
(101) Where Lcirc (r init ) corresponds to the circular angular momenta evaluated at the end of the inspiral, rinit .Using this prescription, we are able to integrate the coupled ODEs Eq.( 100) with initial conditions (101) to obtain Fig. (4).The transition solution smoothly deviates away from the adiabatic inspiral (blue curve), passes through the ISCO and reaches the horizon where the solution terminates.The plot on the right shows the full worldline in proper time r(τ ) where the inspiral starts at r = 1.006 and terminates at the horizon.This method ensures that r(τ ) is both continuous and once differentiable everywhere.
Also, by our choice of integrating (100) using the variable δL, we ensure continuity but not differentiability in L throughout the full inspiral.We note here that Apte and Hughes in [31] also found discontinuities in their evolution of both L and Ẽ and added corrections to ensure both (first order) differentiability and continuity at τinit .We consider a correction of the form We have discussed previously that the leading order term in L(τ ) − Lisco scales proportionally to η 4/5 −2/15 .So we choose to add a constant offset ∆ Lcor ∼ η 6/5 2/15 to the angular momenta evolution to ensure continuity in the Ltrans evolution.
To calculate the evolution in Ẽ, one computes Ẽcirc given by Eq.( 5) during the adiabatic inspiral regime.Then, during the transition regime, one integrates from the flux balance law Ė = Ω(r) L. The correction to ∆ Ẽcor is chosen to ensure continuity with the end of the inspiral energy given by Eq.( 6) as previously discussed after equation (102).
Notice here that this ensures that the energy flux obeys Ė = Ω(r) Lisco and is thus not constant.This ensures that we are still granted a full cancellation of the dissipative part of the forcing term f r in Eq.( 28).This will yield a continuous evolution Ẽ at the matching point with a discontinuous first derivative.At this point we will have a full trajectory r(τ ) with (continuous) integrals of motion in proper time Ẽ(τ ) and L(τ ).In each of [27,30,31], the authors compute three separate worldlines in proper time; Adiabatic inspiral, transition, geodesic plunge.Apte et al in [31], provide an algorithm in which they freeze the constants of motion Ẽ and L when the extra terms in Eq.(64) exceed the leading order terms X 2 and T by 5%.As one would expect, as one ventures farther from the ISCO, the Taylor expansion used to derive these transition equations of motion will break down.As such, it is very natural for each of the aforementioned authors to compute a geodesic plunge to complete their worldlines in proper time r(τ ).Simply because, for moderate spins (non near-extremal), |r + − risco | ∼ O(1) η 2/5 .For nearextreme black holes the ISCO is close to the horizon in Boyer Lindquist coordinates |r isco − r+ | ∼ 2/3 .The scaling of the near-extremal transition zone is also 2/3 and so the horizon is reached while the object is still in the transition zone.We therefore do not expect to need to add a geodesic plunge to compute full near-extremal inspirals.To verify this we numerically calculate the extra terms in (100), which are (104) We compare the solution to (100) when these terms are omitted or included in Figure 3.The difference is at most 1% even the horizon r+ .We conclude that we can use the solution from (100) throughout the plunging regime, for r(τ ) ∈ [r + , risco ].It would be useful in the future to compare our results with the analytic geodesic plunges found in [25].

B. Worldline in Boyer-Lindquist Coordinates
In the previous section, we computed the full worldline comprised of inspiral, transition and plunge parametrized as r(τ ).We now intend to do the same but in coordinate time so that our worldline is in Boyer-Lindquist coordinates ( t, r( t), θ = π/2, φ( t)).Loosely speaking, this is the time measured from Earth (at radial infinity) so is useful for observable purposes.
For the quasi-circular inspiral solution, we simply integrate the circular relation relating coordinate time to proper time via Eq.( 8 where T (r, Ẽ, L, a) is given by 4 and tinsp is defined through t(τ init ).Throughout the transition regime, we use the model for both Ẽ(τ ) and L(τ ) given by Eq.( 103) and Eq.( 102).This will yield the r( t) throughout the transition regime.Combining these results yield a full trajectory from radial infinity to the horizon in coordinate time r( t).
To then calculate the orbital velocity dφ/d t = Ω in coordinate time we substitute r( t) found previously into Eq.( 9).This now gives Ω( t) valid throughout the adiabatic inspiral regime.Using our solutions for Ẽ(τ ) and L(τ ) defined through Eq.(102) and Eq.(103) and r( t) throughout the transition regime, we calculate where ∆(r) = r2 −2r+a 2 .This algorithm will provide a worldline in coordinate time r( t) which will be used for our waveforms.We stress here that r( t) is continuous and (once) differentiable.

C. Near-Extremal Waveform
Following [42], the root mean square (rms) amplitude of gravitational waves emitted towards infinity at harmonic m is given by h o,m = h 2 +,m + h 2 ×,m .The plus and cross each represent individual transverse-traceless polarisations of the gravitational wave strain h.The amplitudes are averaged • over the direction and over the period of the waves.Furthermore, the rms amplitude generated by a particle on an equatorial circular orbit in the limit η → 0 is related to the outgoing radiation flux in harmonic m by for distance D and outgoing fluxes defined by where the amplitude A m equals and Ė∞,m is the relativistic correction to Ė∞,m at each harmonic m.
An EMRI signal is a superposition of infinitely many harmonics of the fundamental frequency Ω with 2π fm = m • Ω. Recall that the total emission of radiation through gravitational waves is related to the outgoing and ingoing flux by where ĖH,m is the ingoing flux (towards the horizon) including the contribution from all l for each harmonic m.Using the exact results from the BHPT for a spin parameter of a = 1−10 −9 , we constructed a cubic spline for each outgoing flux Ė∞,m .Our results are plotted in figure 5.It is clear that including the higher order modes become increasingly important as the spin parameter increases towards unity.This has already been observed in [25].Hence, for near-extremal systems, only using the m = 2 harmonic is not an accurate representation of the EMRI signal in general.Figure 5 suggests that truncating the sum at the eleventh harmonic in the outgoing flux (111) is a good approximation to model a near-extremal waveform encapsulating quasi-circular inspiral, transition, and then plunge with suitable accuracy.The remaining difference from modes with m > 11 contributes a small difference in the amplitude of the waveform, but gravitational wave detectors are much less sensitive to amplitude corrections than corrections to the phase.The phase evolution is determined by the total flux, in which we are including all modes.We therefore believe that the approximate waveform with 11 modes is a sufficiently accurate model for parameter estimation studies and will use this henceforth Once the ISCO is reached, we smoothly extrapolate each of the fluxes Ė∞,m → 0, as r → r+ .This is a similar approach to that found in Taracchini et al in [55].Using (112) and the results obtained in this paper, we plot a near-extremal waveform, including the transition from inspiral to plunge, in Fig. (6).We notice that the waveform in Figure 6 exhibits the usual dampening before the ISCO is reached as seen by Gralla et al in [23].This, qualitatively, is a unique feature to near-extreme EMRIs as a gravitational wave source.

V. CONCLUSION
This paper has presented a solution to the problem of the transition from inspiral to plunge, for any primary spin, for EMRIs on circular and equatorial orbits.This work has extended the treatment of Ori & Thorne [27] which was the first analysis of this problem but did not apply to systems with near-extremal spins.This work also extended the analysis of [29] which did consider near extremal spins, by providing a better physical interpretation of the procedure, identifying a missing term in the analysis and updating the treatment to use recent calculations of the near-extremal energy flux.We have also carefully identified the scaling of the various higher order terms arising from effects such as eccentricity and non-geodesic past-history to carefully demonstrate that these are all sub-dominant.Previous treatments have assumed that the quasi-circular assumption holds throughout the inspiral, but without rigorous justification.We have demonstrated that initial eccentricities excited during the adiabatic inspiral regime grow by the time the transition regime is reached, but are still sufficiently small to be sub-dominant.We have shown that corrections to the flux balance law (29) arising from eccentricity and from the non-geodesic pasthistory of the orbital evolution are also sub-dominant, if only marginally, but there are non-trivial deviations from the linear-in-proper-time evolution of energy and angular momentum in (36) that was assumed in OT.These deviations are encoded in the evolution of the parameter Γ through the transition regime.
Based on these arguments, we have derived a transition equation for each of the three scaling regimes: η , η ∼ and η and described a numerical scheme to generate a full inspiral trajectory in coordinate time, from radial infinity to the horizon.For near-extremal black holes, we found that there was no need to attach a geodesic plunge onto the transition solution as the inspiraling object reaches the horizon while still within the transition regime.Finally, we used these inspiral trajectories to construct a near-extremal waveform exhibiting the transition and plunging dynamics using results from the BHPT [56].
The OT procedure is straightforward, but with surprisingly rich phenomenology.Through semi-analytic means, one is able to derive an equation which describes the dynamics within the vicinity of the ISCO.However, in practice, the OT theory has several shortcomings.The point at which the transition solution is taken to start has a significant influence on the time it takes the particle to reach the horizon and so the OT procedure does not define a unique worldline given a particular set of parameters for the source.This is clearly not physical behaviour.We argued in section IV A that if the switch from the adiabatic inspiral to the transition equation is made when the constraint δL ∼ η 4/5 −2/15 is satisfied, the solution will be almost unique.This was verified numerically and we found it leads to plunge times consistent within ±0.5M .This very same problem was found in [31] but they saw no effect in their waveform analysis.For the η case, there is a further degree of freedom as to when to attach the geodesic plunge.To do so, one must "freeze" the integrals of motion Ẽ, L at the end of the transition regime and integrate the Kerr geodesic equations forward in coordinate time.Attaching the geodesic plunge is discussed, at length, in [31] but does not have a unique solution.Care must be taken as to when the transition solution and the geodesic plunge is attached or comparatively different radial trajectories will be produced.Fortunately for the ∼ η case, there is no need to attach a geodesic plunge as shown earlier in section IV A.
Another issue with the OT method is that it can lead to discontinuities in the constants of motion Ẽ(τ ) and L(τ ) if the OT equations are integrated backwards from the ISCO rather than forwards from the point of the switch from the adiabatic inspiral to transition regime.Discontinuities in the constants of motion lead to discontinuities in the coordinate time trajectories and in the waveforms which must be avoided if these waveforms are to give physically reasonable results in parameter estimation studies.Our solution, which was to integrate forward not backwards, yields continuous, but not first order differentiable, trajectories.The procedure described in [31] provides both.For parameter estimation studies we only require continuity of Ẽ and L and first order differentiability of r( t) and so our procedure should be sufficient, although this should be examined more carefully.
There are natural extensions of this work.First, the waveforms constructed in this paper can be used to carry out a parameter estimation study to understand how well the parameters of near-extremal EMRIs can be measured with observations by LISA.Of particular interest is how well the spin can be determined, since the identification of an object that definitely has spin above the Thorne limit would be of profound significance.Second, our waveforms are missing the quasinormal mode ringdown contribution.Hence, it would be very interesting to generate a full waveform taking these into account, together with the plunging dynamics discussed here.Details on how to construct the waveform including this effect were discussed in [57].Finally, it would also be of interest to extend the current analysis to inspirals that are not circular and equatorial.The extension of OT was first performed in [30] who attempted to solve the problem for generic orbits; both eccentric and inclined orbits.Sundararajan's treatment was corrected by [31] in the case of arbitrary inclina-Figure 6: (Top Plot) Here we plot the root mean square gravitational waveform for both inspiral, transition and plunge using the first eleven harmonics.Notice the smooth evolution of h( t).We terminate evolution of the waveform close to the plunge r = r+ + δ for suitably chosen 0 < δ 1, otherwise the waveform will continue to decay for infinite coordinate time.This is obvious since the (point-like) particle as observed from infinity will never reach the horizon.In this example, we considered a = 1 − 10 −9 and η = 10 −5 so that we are in the ∼ η regime.(Bottom Plot) We plot a zoomed in version of the top plot to show the reader the the smooth evolution of our adiabatic inspiral waves into the transition waves.The faded black dot indicates the moment the transition solution is turned on.
tion.Hence, no one, as of yet, has considered the transition from inspiral to plunge in the case of eccentric orbits and inclined orbits.These orbits are expected for EMRIs formed through standard astrophysical channels [11].The extension to eccentric orbits will require more careful modelling of the self-force and the use of the (eccentricity-dependent) separatrix in place of the ISCO among other complications.A model of the transition for inspirals on generic orbits into black holes arbitrary spin will be invaluable for the analysis of future LISA EMRI observations and is an important future topic of study.

ACKNOWLEDGMENTS
This work made use of data hosted as part of the Black Hole Perturbation Toolkit at bhptoolkit.org.We wish to thank Maarten van de Meent for his comments on the final stages of this manuscript.We also give our thanks to Maarten for pointing out useful literature concerning the conservative part of the self-force and giving direction on why we could neglect it.We also wish to thank Niels Warburton for his guidance using the toolkit and Peter Zimmerman for his useful discussions regarding the behaviour of near-extremal Kerr Black holes.Finally, we thank Geoffrey Compère for comments on an earlier version of this manuscript.where δ n0 = 1 for n = 0 and zero otherwise.Finally, evaluating these derivatives at (r isco , Ẽisco , Lisco ) and using the properties (A2)-(A7), we can derive the exact results Let us study the behaviour of these derivatives for near extremal black holes, i.e. in the limit → 0 as introduced in section II.Remember risco is given by risco → 1 + 2 1/3 2/3 + 7 4 • 2 1/3 4/3 + O( 2 ) .(A12) Using this expansion together with a = √ 1 − 2 , we can evaluate the leading terms of all previous derivatives to be into Eqs.(6), ( 5), ( 14) and ( 9  The coeffcients in (B6) can be approximated for → 0 under the retrograde condition Eq.(B1) and Γ in (B6) defined through equation (41).Notice that none of the coefficients in our transition equation of motion depend on the extremality parameter .This gives us no reason to introduce any scalings on r, τ and δL like we did for prograde orbits around rapdily rotating black holes.As such, let us introduce similar scalings to OT R = η 2/5 α −3/5 (βκ) for K 0 = α 4/5 (βκ) −6/5 .Which is precisely the equation of motion for the transition regime derived by OT in [27].Although the quantities α, β and κ present in the change of coordinates are different, the physics and ultimate end goal are the same.As a result, we stop our analysis of retrograde orbits here since we feel that this problem has already been solved by the community for smaller spin values a ≥ −0.999.We conclude that, for near-extremal retrograde orbits, there is nothing new to learn about the transition regime.It can be solved in the matter of OT in [27].We do remark that the quantity κ can no longer be computed using the near-extremal formula defined by ĖGW = ( CH + C∞ )(r− r+ )/r + .This is because the transition region is far from the horizon of the primary hole [see figure (1)].Instead we have to use the numerical quantity Various results are tabulated (including retrograde orbits) in [42].The downside of this equation is that it can only be evaluated numerically.

Appendix C: Osculating Elements Equations
The proper time derivative of the radial geodesic equation ( 2 The purpose of this appendix is to review how this equation is equivalent to the radial component of a forced geodesic equation where xµ = (r, t, θ, φ), ∇ ν = ∇ xν , u µ = dx µ /dτ is the four-velocity of the particle and f µ a forcing term driving deviations from geodesic motion.
To show the equivalence between (C1) and (C2), we use the osculating elements formulation [58].Since this method does not take into account conservative effects arising from the self-force [46], the component f r in this appendix will only account for the dissipative piece in (28).Its conservative piece is treated in more detail in the main text (See section III A).
Since the four velocity u α is normalised, it follows f α is normal to it by proper time differentiation Evaluating (C2) along the radial direction, solving (C3) for f r and plugging it into (C2) yields Ẽ dτ , dτ .

Figure 3 :
Figure 3: Solution to (100) (black dashed line) and difference in solution when including the higher-order corrections given in Eq. (104) (red solid line).The numerical difference is small throughout the transition regime reaching a maximum at plunge of ∼ 1%.

Figure 4 :
Figure 4: The red curve shows the orbital velocity Ω and the black curve shows the trajectory in coordinate time r( t).Notice the smooth evolution of both r( t) and Ω during the start of the transition (green dashed curve).This smooth evolution continues through the ISCO (blue dashed curve) and evolves towards the horizon (black dashed curve).

Figure 5 :
Figure 5: Comparison of the total energy flux at infinity (black curve) including different harmonic Ė∞,m contributions.Note that at r ≈ 1.3, the m = 2 harmonic energy flux Ė∞,2 contributes ∼ 32% of the total energy flux, whereas including the first 11 harmonics (violet curve) contributes ∼ 98% at the least.
Notice here that the expansion in is no longer increasing in powers of 2/3 and now in2.Also notice that |r isco − r+ | ∼ O(1) rather than of order 2/3 like in the case of near-extremal prograde orbits.Like we have done previously, we consider the Kerr radial velocity expanded around the ISCO dR dτ