Adiabatic inspirals under electromagnetic radiation reaction on Kerr spacetime

A compact body in orbit about a black hole loses orbital energy and angular momentum through radiation-reaction processes, inspiralling towards the black hole until a final plunge. Here we consider a scenario with a charged compact body in which fluxes of electromagnetic radiation drive this inspiral. We calculate trajectories in the $(p,e)$ plane for inspirals in the equatorial plane of a rotating black hole within the adiabatic (orbital-averaged-dissipative) approximation. We make comparisons with a non-relativistic Keplerian approximation based on the Abraham-Lorentz force law, and with standard gravitational-wave driven scenarios. We find that EM-driven inspirals are less efficiently circularized (i.e. orbits remain more eccentric at the point of plunge) than their gravitational counterparts, and we quantify the effect of black hole spin.


I. INTRODUCTION
Since the celebrated first-light detection of 2015 [1], networks of gravitational-wave interferometers have observed numerous "chirp" signals originating from compact binaries during the final stages of inspiral and coalescence [2,3].The emission of gravitational waves, which is strongest at periapsis, removes eccentricity from a binary system, causing the orbits to become more circular as the inspiral progresses [4,5].For comparable-mass binaries -those accessible in the frequency band of groundbased detectors -the circularisation process is efficient, and by the point of detection the bodies are on nearly-circular orbits.Thus, the observed signals are regular "chirps" characterised by a small number of parameters (chiefly, the chirp mass) [6].Conversely, extreme mass-ratio inspirals (EMRIs) are less efficiently circularised, and the orbit of the secondary is expected to remain significantly eccentric and non-equatorial all the way up to the plunge [7].EMRIs are key targets for space-based detectors (i.e.LISA), which can probe much low frequencies due to the absence of seismic noise in space.
One approach to modelling EMRIs is by calculating the so-called self-force [8][9][10].The essential idea is that the emission of gravitational waves is associated with a back-reaction on the motion of the system that is described by a force acting on the secondary (of mass m) moving on the spacetime of the primary (typically, a black hole of mass M and spin a = J/M ).A complementary view is that the motion is geodesic in a certain effective spacetime [11].The self-force programme is welldeveloped at leading order in m/M [9,12], with first results at second order in the mass ratio now available [13][14][15].
Self-force, as a concept, has historical roots in electromagnetism: it builds upon Dirac's idea of separating the field into its symmetric/singular and radiative/regular parts, with the latter generating the self-force [8,11,16].Since 1997, the focus of the programme has been on calculating gravitational self force (GSF), because this is the predominant driver of astrophysical binaries.However, the electromagnetic self force (ESF) and radiation reaction (associated with the dissipative part of ESF) is interesting in its own right, and has received some attention in the literature of recent years [17][18][19][20][21][22][23][24].In this work, we seek to extend the ESF frontier by studying the evolution of the eccentric, equatorial orbit of a charged particle around a Kerr black hole for the first time (see Ref. [23] the circular orbit case).
There are at least two good reasons for considering the electromagnetic self force in binary systems.Firstly, it provides an opportunity to compare and contrast an electromagnetically-driven scenario with the standard gravitationally-driven scenario.A key difference is that electromagneticallydriven systems radiate principally in the dipole (ℓ = 1) mode, whereas gravitationally-driven systems radiate principally in the quadrupole (ℓ = 2) mode.A consequence is that the rate of the increase of the orbital frequency during the inspiral, f (t) ∝ (−t) −β (for r ≫ M ), is markedly different, with an index β = 3/8 in the gravitational case and an index β = 1/2 in the electromagnetic case [2,23].
A second consequence is that the orbit is less efficiently circularised by the ESF, as we investigate in more detail below.
Another good reason for considering ESF effects is to improve the modelling of astrophysical EMRIs.Firstly, compact bodies could sustain electromagnetic charges that could conceivably cause some dephasing over the final circa 10 4 -10 5 orbital cycles in the final year before plunge.Secondly, the presence of a dipole mode in the electromagnetic sector makes it a reasonable proxy for modelling EMRI evolution in beyond-GR theories.Theories with a radiative dipole that have been considered in the EMRI context include scalar-(vector)-tensor theories [25], massive-gravity theories [26], and Einstein-Maxwell-Dilaton theory [27].Recent work by Barsanti et al. [28,29] has shown that, in in scenarios with a new fundamental scalar field, the secondary experiences a scalar self-force of exactly the kind earlier considered as a toy model in the literature (e.g. in Refs.[30,31]).It seems plausible (but is not shown here) that electromagnetic self force may play the same role in extended theories with a radiative spin-one degree of freedom.
In the 1960s, the problem of a particle undergoing a gravitationally-driven inspiral into a black hole was first tackled using a Keplerian approximation to the underlying dynamics [4,5].The Keplerian approximation -which can be thought of as the leading term in a post-Newtonian expansion -gives a robust model for understanding the inspiral over much of its life.In this paper, we extend the Keplerian treatment of a binary black hole system to one where the particle's inspiral is driven by electromagnetism through the Abraham-Lorenz force, and then we proceed to a full-blooded radiation-reaction approach to calculate more accurate inspiral trajectories.
The paper is organised as follows.In section II A, we derive a Keplerian approximation to the problem of an electromagnetically-driven inspiral, and we recap the results in [4,5] for a gravitationally-driven inspiral.In Sec.II B, we summarise the calculation of the fluxes of energy and angular momentum radiated by particles on eccentric equatorial orbits in Kerr spacetime.The adiabatic inspiral model, in which the orbital parameters (p, e) evolve in proportion to the fluxes, is described in Sec.II C. In Sec.II D, we outline the application of Chebyshev interpolation to calculate fluxes across the (p, e) parameter space from values at node points.Various implementation details are covered in Sec.II E. A selection of numerical results are presented in Sec.III, including trajectories of inspirals in the (p, e) plane.The article concludes in Sec.IV with a recap of the key results, and an outlook on avenues for further exploration.

II. ANALYSIS A. Keplerian approximation
For a Keplerian orbit, the position vector is x(t) = r(t) cos ϕ(t) i + r(t) sin ϕ(t) j, where Here p = a(1 − e 2 ) is the semi-latus rectum, a is the semi-major axis, and e is the eccentricity of the orbit.In the absence of other forces, the specific orbital energy E and angular momentum L are conserved, and the orbital parameters (p, e) are constant.In the adiabatic approximation, the radiative losses cause E and L to change slowly (relative to the orbital timescale), and ⟨ Ė⟩ and ⟨ L⟩ are found by taking an average over one orbit.

Electromagnetic-driven inspirals
Invoking the Abraham-Lorentz [32,33] force law, where v = ẋ and a = v, the (instantaneous) rate of loss of energy and angular momentum are given by Ė = −v • F and L = −x × F, respectively.Taking an average over one orbit, Here we have adopted units in which GM = 1, m = 1, c = 1 and 4πϵ 0 = 1.Consequently, and so This ODE has an elementary solution, where p i and e i are initial values.Hence for a Schwarzschild black hole, the eccentricity at plunge, e f , may be estimated from the iterative formula The characteristic timescale for a circular-orbit inspiral is, from Eq. (4a), T = p 3 i /(4q 2 ).

Gravitational-wave-driven inspirals
For comparison with the above, we here quote the corresponding formulae for inspirals driven by gravitational radiation, also in the Keplerian approximation (see Refs. [4,5]), and for m ≪ M .The rate of change of orbital parameters is where where c 0 is a constant of integration.By comparing the index 4/3 in Eq. ( 7) with the index 12/19 in Eq. ( 12), it is clear that the emission of gravitational waves will "circularize" the orbit more efficiently than the emission of electromagnetic waves (at least within this approximation).

B. Flux calculation on Kerr spacetime
The approximations derived in the previous direction are valid in the far-field, p ≫ M .For a more accurate study, it is necessary to compute the fluxes for orbits in the strong-field region, all the way up to the plunge.This is done by applying the Teukolsky formalism [34][35][36][37][38][39].

C. Adiabatic inspiral model
In the adiabatic approximation, we assume that the particle's worldline is, at each moment, locally tangent to an equatorial eccentric geodesic with orbital parameters p and e (the semi-latus rectum and eccentricity, respectively).These orbital parameters change smoothly with time -and slowly in comparison with the orbital timescale -as the system radiates away energy and angular momentum [40].Thus each inspiral corresponds to a trajectory in the (p, e)-plane.For a particle on an eccentric orbit, the radius of orbit r 0 is related to the orbit's semi-latus rectum p and eccentricity e by where χ is a monotonically increasing function of proper time τ .
The energy and angular momentum of a particle on an equatorial (geodesic) orbit can be expressed as functions of the orbital parameters, as E = E(p, e) and L = L(p, e).The form of these functions is specified in Appendix B. Thus, we can relate Ė and L to ṗ and ė by i.e., a system of differential equations that can be written in matrix form and inverted giving where the partial derivatives are calculated from Eq. (B1) and Eq.(B2).Thus Ė and L determine how e and p change.If we specify Ė and L, along with a set of initial conditions, we can solve this system of equations with an ODE integrator, such as Mathematica's NDSolve[].
Using a flux-balancing argument, in the adiabatic approximation the loss of orbital energy (orbital angular momentum) matches the averaged radiated energy flux (angular momentum flux).Hence the quantities Ė and L are given by the corresponding fluxes produced: This completes the scheme for being able to calculate inspiral trajectories in the (p, e) domain.

D. Chebyshev interpolation
To model an inspiral, we need fluxes across a domain in (p, e) for a given spin parameter a.The numerical evaluation of Φ (E/L) ∞/h in Sec.II B comes with a computational cost.We wish to be able to calculate these at any arbitrary point on the (p, e) plane rapidly, so as to be able to calculate multiple trajectories.An efficient way to achieve this is by using interpolation with Chebyshev nodes to generate a polynomial that approximates the fluxes, as in Ref. [41].The resulting interpolation polynomial minimises the effect of Runge's phenomenon (i.e. it minimises spurious oscillations near the boundaries of the interpolated region).
A Chebyshev polynomial of the first kind is a polynomial T n (x) such that T n (cos θ) = cos(nθ).The roots of these polynomials x k are known as Chebysehev nodes and are given by x k = cos 2k−1 2n π where n is the degree of the polynomial and k is an integer between 1 and n.For our application, we wish to formulate a Chebyshev interpolation function for a function of two variables, f (x, y).We wish to obtain an approximation of the form: where u and v are scaled versions of our variables x, y such that u, v ∈ [−1, 1] and f ij are the associated Chebyshev coefficients.The procedure for calculating these coefficients is: 1. Set up an n × m grid of Chebyshev nodes (x i , y j ) on the domain of interest.
2. At each node calculate the value of the function f (x i , y j ), and plug into Eq.( 22) giving: This is a set of n × m equations for n × m unknowns.
3. Solve this system of equations using Mathematica's NSolve[] function to obtain f ij .

E. Implementation details
These methods were implemented in Mathematica.A 25×20 grid of nodes was chosen to perform the interpolation.We define a parameter y = p separ (a, e)/p.This transforms the pe-plane into the ye-plane.We do this transformation to only calculate nodes beyond the Separatrix, the line in the pe-plane that separates orbits from plunging trajectories.We then map this plane onto the [−1, 1] 2 domain upon which the Chebyshev polynomials are defined using the following transformations: where the subscripts indicate the minimum and maximum values we want to consider on the plane.
We choose e min = 0 and e max = 0.6, and y max = 1, corresponding to the point on the separatrix, and y min = p separ (a, 0.5)/150, which is a value of y that roughly corresponds with p = 150 as a maximum p value.We then define the Chebyshev nodes on this new uv-plane and invert the transformations to obtain the nodes in the pe-plane, these are illustrated in Fig. 1.
In principle, the total fluxes in Eq. (13a) and Eq.(13b) are formed from infinite sums over ℓmn modes.These partial fluxes are typically monotonically decreasing in magnitude with |n|, and diminish in an exponential fashion.Thus when calculating the final flux we can truncate the sum to some N max and L max values without significant loss of accuracy.Furthermore, there is a symmetry that can be exploited to increase calculation speed two-fold.For a given flux mode F ℓmn we have that F ℓmn = F ℓ−m−n , and that F ℓ00 = 0. Applying this symmetry to the sum over flux modes we get that the total flux Φ will be given by: Each of the four fluxes were calculated at each node point for a/M = 0, 0.25, 0.5, 0.75, 0.9, and q = 1.
Once calculated, we plotted each set of flux points against its p-value.This is illustrated in Fig. 2. We can fit a straight line to each of these, the coefficient of which can be used as a smoothing parameter on the flux to give a flatter function before applying the Chebyshev interpolation.Numerically, we obtain smoothing parameters (s E∞ , s Eh , s L∞ , s Lh ) = (4.16,6.53, 2.68, 5.14) for a = 0, which are close to the indices anticipated from the post-Newtonian expansions.We then define the function we will interpolate as f (p, e) = p s (i) Φ (i) (p, e), and then the final interpolated fluxes Φ  With all this in place we can now calculate inspiral trajectories by feeding Eqs.(21) with the interpolated fluxes, and solving the ODE system (20) numerically using NDSolve[].

III. RESULTS
Figure 3 shows the calculated inspiral trajectory in the (p, e) plane for initial values e i = 0.4, p i = 100 in the Schwarzschild case (a = 0).In the far field (p ≫ M ) there is excellent agreement between the Keplerian approximation and the numerical results obtained from the flux-balancing argument.As the particle gets close to the black hole, some deviation begins to occur, with the trajectory first dipping below the eccentricity predicted by the Keplerian approximation, before an uptick in eccentricity in the final stages of the inspiral before the plunge.The uptick in eccentricity is familiar from the gravitational case [41].
As expected, the inspiral speeds up as the particle gets closer to the black hole and the radiated fluxes increase in magnitude.This progression is visible in Fig. 3 by the number of orbits before plunge (for the case q/m = 1); the calculation of this quantity is described in Appendix C.  Figure 4 shows electromagnetic inspiral trajectories for various initial eccentricities.In every case, the orbit becomes more circular as the inspiral progresses.However, circularization due to EM emission is less efficient than circularization due to gravitational-wave emission, and significant eccentricity is retained at the point of plunge.
Figure 5 shows inspiral trajectories for rotating (Kerr) black holes.The trajectories for a ̸ = 0 are qualitatively similar in nature to the Schwarzschild trajectory, again exhibiting a slightly faster loss of eccentricity than predicted by the Keplerian approximation 7, before an uptick in eccentricity in the approach to the plunge phase (i.e.near the separatrix).The position of the separatrix depends on the spin of the black hole a/M .In the rapidly-spinning case, the co-rotating orbits make a closer approach to the black hole before plunging than the counter-rotating orbits.
As shown above, the Keplerian approximation (7) describes the inspiral accurately at large p, and robustly up to p ∼ 10, but it fails to capture the uptick in eccentricity before the plunge.Consequently, Eq. ( 8) does not give a particularly accurate prediction of the eccentricity parameter at the plunge.Figure 6 shows ∆e (defined in the caption), which characterises the additional eccentricity acquired in the uptick.As shown, the effect of the uptick diminishes somewhat as the  7).The labels on the plot indicate the number of orbits remaining until plunge, in the case q/m = 1 (see Appendix C).
spin of the black hole increases.
Inspirals in astrophysics are driven principally by the emission of gravitational waves, rather than electromagnetic waves.This motivates us to consider inspirals that are driven by both electromagnetic and gravitational interactions.We define the total mixed flux as: where µ is a mixing parameter taking values on [0, 1], with µ = 0 corresponding to gravitational flux only and µ = 1 corresponding to electromagnetic flux only.Employing the mixed flux to calculate inspirals driven by the combined action of electromagnetism and gravity leads to the results shown in Fig. 7.Here we see that the gravitationally-driven inspiral loses eccentricity faster than the electromagnetically-driven inspiral, as expected.All the trajectories show an uptick in the parameter e in the approach to the separatrix.FIG.4: Inspiral trajectories in the (p, e) plane for a selection of initial values of e i starting at p i = 100 (Schwarzschild case a = 0), comparing the adiabatic approximation (solid lines) with the Keplerian approximation (dotted lines).The inset shows the uptick in eccentricity as the trajectory approaches the separatrix.

IV. DISCUSSION AND CONCLUSIONS
In the preceding sections, we have examined the adiabatic evolution of the orbital parameters of binary inspiral driven by radiation reaction in electromagnetism.Figures 3,4,5,6,7 show trajectories in the (p, e) plane, computed by interpolating the fluxes calculated through Teukolsky formalism (see Sec. II B).We find that the leading-order Keplerian approximation, Eq. ( 7), gives a robust characterisation in the weak-field regime (Fig. 3 and 4).Near the plunge, there is a characteristic uptick in the eccentricity parameter e (see Figs. 4 and 6; a similar uptick is also observed in the gravitationally-driven case [42]).The spin of the central black hole has some effect on the evolution of the inspiral in the strong-field regime, most notably in Fig. 5 by shifting the innermost stable circular orbit that marks the transition to plunge (shifting inwards for prograde orbits, and outwards for retrograde orbits).As expected, the electromagnetic radiation-reaction causes the orbit to become less eccentric over time, but this circularization process is observed to be less efficient than in the gravitationally-driven case (Fig. 7), which is consistent with the Keplerian approximation (see Eq. ( 7) and Eq. ( 12)).We expect this observation to hold true for any field theory with a radiative dipole mode.The inset shows the uptick in eccentricity in the approach to the separatrices (black solid lines).
A few natural extensions of this work present themselves.Firstly, it would be straightforward to include the interpolated EM flux model in the modular Fast EMRI Waveforms (FEW) framework [43].This would allow for fully Bayesian Markov Chain Monte Carlo-based parameter estimation studies in the time- [12] or frequency-domain [44] that could determine the smallest measurable charge on the secondary.In fact, there has been recent work in this direction in Ref. [24] using circular orbits, based on an assessment of the dephasing of gravitational waves over the lifetime of the binary.Extending such efforts to eccentric orbits is a natural next step.
If residual charge were present in the binary, then it possible that the primary could also be charged, leading to a modified inspiral due to (i) direct EM interactions between the primary and secondary, and (ii) modifications of the spacetime due to the charge of the primary.Studying these effects in full would require the study of coupled EM and gravitational perturbations of the Kerr-Newman spacetime, which is a significant undertaking (see e.g.[45]).
A second natural extension would be to consider non-equatorial orbits.Very recently the Teukolsky package in the Black Hole Perturbation Toolkit was extended by one of us to compute the EM perturbation for a particle moving on generic bound orbit of a Kerr black hole.Thus, computing the EM flux for spherical (fixed Boyer-Lindquist radius) or generic bound orbits is a straightforward, if computational demanding, task.Calculating the evolution of non-equatorial orbits adds a new complexity, as now the evolution of the Carter constant must also be obtained.
For spherical inspirals the zero-eccentricity constraint allows an evolution equation of the Carter constant to be obtained in terms of the asymptotic fluxes [46].For generic inspirals (eccentric and inclined) no such balance law can be obtained in general.Fortunately, so long as the binary is evolving adiabatically, Mino showed how the correct radiation reaction force for scalar, EM, or gravitational perturbations can be derived [47].This result was then used by others to provide explicit, practical evolution equations for the Carter constant in the scalar [48] and gravitational [49,50] cases.To the best of our knowledge, work remains to be done to extend these results to the electromagnetic case.
Thirdly, there is clearly scope for obtaining post-Newtonian expansions of fluxes and/or the local self-force.From these, closed-form approximations for p(e) would follow, to go substantially beyond the Keplerian approximation of Sec.II A. In the gravitational case, there has been impressive progress in calculating fluxes (e.g.Refs.[51][52][53]), the self-force acting on the particle (e.g.Ref. [54]),  7) and (12).The mixing parameter µ is defined in Eq. (26).
and certain gauge-invariant scalars (see e.g.[55,56]) at very high orders in v/c, at leading order in m/M .The technical machinery for expanding solutions of the Teukolsky equations should transfer straightforwardly from the s = 2 case to the s = 1 case, and so the prospects for such calculations appear to be good.
One limitation of the adiabatic approximation, employed here, is that conservative effects of the self-force are neglected.At leading order in r −1 the conservative part of the self-force is [17] F cons ≈ q 2 4πϵ 0 c 2 • GM r 3 r.The conservative component oscillates over one orbit, unlike the dissipative force which drives the secular change in E and L. Including the conservative force is necessary for an accurate calculation of the orbital phase, and thus dephasing.An obvious next step is to calculate the local self-force on the charged particle using e.g. the ℓ-mode regularization scheme (as was done for Schwarzschild circular orbits in Ref. [22], Schwarzschild eccentric orbits in Ref. [21], and Kerr circular orbits in Ref. [23]).This would enable the calculation of the shift in the innermost stable circular orbit (ISCO) under the conservative piece of the electromagnetic self-force; its gravitational counterpart is already well characterised [57,58].
Finally, in this work we have restricted our attention to binaries where the primary spin |a| ≤ 0.9M .Astrophysically, we might expect spins higher spins up to |a| = 0.998M .Beyond this, the near-horizon, near-extremal (|a| > 0.999) is also interesting as new physical effect occur in this regime.Previous work has both analytically and numerically calculated how the flux behaves in this case for an orbiting scalar and gravitational charge [59,60].The calculations could be extended to the EM case.
= f (p, e)/p s (i) , where (i) indicates the type of flux.Full range of nodes.

FIG. 2 :
FIG. 2: The absolute magnitude of fluxes as a function of semi-latus rectum p, for eccentricities e ∈ [0, 0.6].Each line is increasing in e from left to right.Shown are the energy and angular momentum flux at infinity (black and red) and at the horizon (blue and green).

FIG. 3 :
FIG.3: Inspiral trajectory in the (p, e) plane for p i = 100, e i = 0.4 and a = 0.The vertical dotted line is the separatrix, marking the onset of the plunge phase.The horizontal dotted line is the Keplerian approximation, Eq. (7).The labels on the plot indicate the number of orbits remaining until plunge, in the case q/m = 1 (see Appendix C).

9 FIG. 5 :
FIG.5: Inspiral trajectories in the (p, e) plane for initial values p i = 100M , e i = 0.5, for Kerr black holes of varying spin parameter a/M .The Keplerian approximation(7) is shown as a dashed line.The inset shows the uptick in eccentricity in the approach to the separatrices (black solid lines).

9 FIG. 6 :
FIG.6:The difference ∆e ≡ e f − e Kf evaluated on the separatrix, where e f is evaluated on the adiabatic inspiral trajectory, and e Kf is evaluated with the Keplerian inspiral approximation(8).This is shown as a function of the initial value of eccentricity e i at p i = 100, for black holes of spin a/M ∈ {0, 0.5, 0.9} .