The observation of gravitational waves (GWs) with the LIGO-Virgo [1][2][3] and LIGO-Virgo-KAGRA (LVK) [4] collaborations marked a new era in gravitational physics, uncovering unique properties of stellar-mass black holes (BHs) and neutron stars. As future data acquisition becomes characterized by increased sensitivity, it is necessary to improve the precision and accuracy of waveform models used for matched-filtering and parameter-estimation pipelines. In particular, modeling waveforms from eccentric and precessing-spin binaries will become increasingly important. This is further motivated by the fact that upcoming observational runs [5] and future detectors, like the Einstein Telescope [6], Cosmic Explorer [7] and LISA [8], will increase the number of detections by a factor ∼10 3 and be able to probe binary's subpopulations, in lower frequency bands and for smaller mass ratios, which can exhibit larger eccentricities [9][10][11][12].
While eccentricity decreases toward merger due to the energy and angular momentum loss caused by the emission of GWs [13,14], the residual eccentricity can help constrain different binary-formation scenarios and thus the origin of GW sources [15][16][17][18]. Indeed, the eccentricity is indicative of binaries formed through dynamical formation channels, which could occur in dense stellar environments, like globular clusters, where the three-body Kozai-Lidov mechanism [19,20] or dynamic capture [18,[21][22][23][24] play a role. Efforts are currently underway to detect signs of orbital eccentricity in the GW signals observed by the LVK detectors [25][26][27][28][29][30][31][32][33][34][35].
Among the different methods to solve the two-body problem in general relativity, the effective-one-body (EOB) approach [36,37] is a framework that provides accurate and fast waveforms for quasicircular (QC) binaries [38][39][40][41][42][43][44][45][46][47][48][49], due to a strong synergy between analytical approximation methods and numerical relativity (NR) results.
In recent years, generalizations of EOB models to eccentric inspirals have been developed [50][51][52][53][54][55][56][57][58][59][60][61]. In particular, Ref. [60] derived the second-order post-Newtonian (PN) expressions for the radiation-reaction (RR) force and gravitational-waveform modes for eccentric inspirals. They include tail effects, in addition to spin-orbit (SO) and spin-spin (SS) couplings. Reference [61] introduced the SEOBNRv4EHM model: an extension of the QC SEOBNRv4HM model [44] to eccentric orbits, where the authors showed an EOB/NR unfaithfulness less than 1% when comparing with the 28 eccentric NR simulations that were publicly available at the time from the Simulating eXtreme Spacetimes (SXS) Collaboration [62,63]. 1 However, in the work of Ref. [61], the dynamics of the binary is modeled through the use of the QC RR force from the SEOBNRv4HM model, and the eccentric corrections are considered only when computing the eccentric waveforms modes introduced in Ref. [60]. By contrast, the more recent waveform model SEOBNRv5EHM employs a RR force with eccentric corrections [64,65].
Among the other examples of eccentric EOB models, we mention TEOBResumS [48,66], which has been extended to eccentric orbits, after investigating several prescriptions for incorporating eccentricity effects, in Refs. [56][57][58][59][67][68][69][70][71][72][73]. The latest version of their model, known as TEOBResumS-Dalì, includes eccentric 2PN information, and is characterized by factorizing the leading PN order of the waveform modes and azimuthal component of the RR force, which include high-order time derivatives of the radial separation and orbital frequency. The eccentric 2PN radial component of the RR force is adapted from Ref. [50] and is Padé resummed. When compared against the 28 SXS publicly available eccentric NR simulations, the TEOBResumS-Dalì model shows unfaithfulness around ∼0.1%, although there are some NR waveforms for which the unfaithfulness is close to, or above ∼1% [74]. However, note that these unfaithfulness values cannot be directly compared against the ones computed in Ref. [61] for the SEOBNRv4EHM model. This is because Refs. [61,74] employ different prescriptions for determining the EOB waveform that corresponds to a given NR simulation.
A direct comparison of these models employing the same prescription is presented in Ref. [65].
Several studies [38,58,68,71,[75][76][77][78][79] showed the importance of augmenting EOB waveform models for the plunge-merger and ringdown with insights from BHperturbation theory. The common approach is to consider a test mass (TM) orbiting or scattering off a Kerr BH and use this system as a laboratory to test and provide benchmarks to the models in the comparable-mass case. Among these works, Refs. [58,68] assessed different EOB eccentric RR force prescriptions (either with distinct RR gauge choices or different factorizations of the RR force). In particular, Ref. [68] also provides an analysis for a proxy to the SEOBNRv4HM QC RR force and of a resummed version of the 2PN eccentric RR force introduced in Ref. [60]. (We will discuss some comparisons with their results in Sec. III A.)
Here, we aim to extend past analyses in different ways. We consider the TM limit of a 3PN-eccentric RR force, recently computed by some of the authors of this work in Ref. [64]. This force is computed employing the same procedure of Ref. [60], but it considers a different gauge choice for the leading-order of the RR force, which avoids a 2.5PN modification (relative to the leading-order) of the QC orbital phase when transforming between harmonic and EOB coordinates. The force that we consider contains the full nonspinning contributions up to 3PN order, and spin contributions only up to 2PN order. For the spin contributions, we employ the 1.5PN spin-orbit (SO) and 2PN spin-spin (SS) parts to the RR force as in Ref. [60], but taking into account the leading-order gauge choice of Ref. [64]. We resum the RR force in two ways: by extracting the QC RR force of SEOBNRv5HM [49] as a multiplicative and an additive term. In the TM limit, we remark that this QC RR force differs from the SEOBNRv4HM one by new higher-order PN contributions in the waveform modes used for the flux computation, as explained in Ref. [49]. We study the effects of the individual PN eccentric contributions to the RR force in different gravitational-field regimes of the parameter space. The analysis is performed by comparing the analytical fluxes, computed from the eccentric RR force in the TM limit, against numerical fluxes that are computed by solving the Teukolsky equation [80] through the use of a frequency-domain (FD) [81] and a time-domain (TD) [82][83][84] code. Both fluxes are computed on equatorial geodesics of the Kerr metric and extracted at future null infinity. By comparing the fluxes, we test the two resummations of the eccentric RR force, as we explain in detail in Sec. II. In our study, we focus on equatorial bound orbits of Kerr in the weak and strong gravitational-field regimes, and also explore the RR force for hyperbolic encounters, by comparing the fluxes for some Schwarzschild hyperbolic geodesics with fixed energy.
The article is structured as follows. In Sec. II, we introduce the methodology of our analysis. In particular, 1 The unfaithfulness is a metric that quantifies the disparity between two waveforms as observed by a GW detector (a lower unfaithfulness indicates greater similarity between the waveforms). The definition of the unfaithfulness is presented in, e.g., Ref. [61]. To determine the accuracy of the SEOBNRv4EHM model, Ref. [61] employed an optimization over eccentricity and starting frequency to find the best-fitting EOB waveform corresponding to a given eccentric NR waveform. Thus, the SEOBNRv4EHM model showed an unfaithfulness against NR waveforms less than 1% for systems with eccentricities below ∼0.3.
Sec. II A describes the EOB model we use in the TM limit and how we compute eccentric-planar geodesics in the Kerr metric. In Sec. II B, we describe how to compute the numerical fluxes by solving the Teukolsky equation numerically, while Sec. II C shows how the analytical fluxes are computed from the RR force through the use of the balance equations. In Sec. III, we provide the main results of our analysis. In particular, Sec. III A gives an overview of the fluxes comparison, and Secs. III B and III C show the results for bound orbits in the Schwarzschild and Kerr spacetime, respectively. Section III D focuses on hyperbolic encounters in Schwarzschild spacetime. Finally, Sec. IV summarizes the results, and points out future steps. In the Appendices, we provide supplemental information. Notably, in Appendix A we summarize how the nonspinning 3PN terms of the eccentric RR force are derived and we provide the full TM expressions of the eccentric corrections to the QC RR force, together with the expressions of the Schott terms, which are necessary to compute the instantaneous fluxes. Appendix B provides details on how the instantaneous fluxes are computed from the FD Teukolsky-equation code.
We adopt natural units G ¼ c ¼ 1 and consider a nonspinning TM of mass μ ¼ νM orbiting a Kerr BH of mass M with dimensionless spin a ¼ J=M 2 . The Kerr metric is expressed in Boyer-Lindquist coordinates fT; R; θ; φg and we restrict our analysis to the equatorial plane, θ ¼ π=2. The dynamics of the TM is described by canonical coordinates fR; φ; P R ; P φ g. Throughout this article, we consider scaled dimensionless variables
The Hamiltonian H, and RR force F ¼ ðF r ; F φ Þ are scaled by the TM μ.
In this work, we assess the analytical EOB eccentric RR force of Ref. [60], extended to 3PN in the nonspinning part in Ref. [64], by comparing it with numerical results. In particular, as we explain in Sec. II C, we compare numerical fluxes obtained by solving the Teukolsky equation against the analytical fluxes that are connected to the RR force through the energy and angular-momentum balance equations. These fluxes are computed on Kerr equatorial geodesics of a TM. In the following two sections, we describe how the orbits are computed and the methodology used to derive the numerical and analytical fluxes.
To describe the dynamics of a TM orbiting a Kerr BH in the equatorial plane, we work within the EOB framework [36,37] and consider the Kerr Hamiltonian restricted to equatorial orbits (θ ¼ π=2, p θ ¼ 0):
with quantities Λ and Δ being
Instead of the radial momentum p r we consider p r à , which is the momentum conjugate to the tortoise radial coordinate r à . The tortoise coordinate is related to the Boyer-Lindquist coordinate r by
This is a general practice [76,85] that is done to improve the numerical stability of the dynamical evolution, since p r diverges at the horizon while p r à does not. The evolution of the dynamics is provided by the Hamilton equations
where the dot symbol represents a total derivative with respect to the scaled coordinate time in Eq. ( 1), Ω is the orbital frequency, scaled by the total mass, and F ¼ ðF r ; F φ Þ corresponds to the RR force connected to the emission of GWs for generic equatorial orbits.
The RR force components we consider, F r and F φ , are two resummed versions of the RR force originally computed in Ref. [60] and here extended to 3PN order in the nonspinning eccentric sector [64]. These two resummations are given by
where F QC r and F QC φ are the radial and azimuthal components of the RR force using the QC prescription, while F ecc;mult r;φ and F ecc;add r;φ are the eccentric corrections. We refer to the two resummations as the multiplicative implementation, given in Eq. (6a) and the additive implementation, given in Eq. (6b). The complete expressions in the TM limit of the eccentric corrections are in Appendix A. Note that we restrict our attention to resummations that reduce to SEOBNRv5HM in the QC limit. This precludes considering resummations of the type introduced in [56], which also change the QC radiation reaction.
The QC RR force F QC ¼ ðF QC r ; F QC φ Þ is calculated using the prescription and PN information of the SEOBNRv5HM waveform model [49], defined by the expressions
where d L is the luminosity distance of the binary to the observer and h F lm are the GW modes in factorized form [38,39,76,86], given by:
Here, ϵ is the parity of the multipolar waveform mode, such that ϵ ¼ 0 for even l þ m, and ϵ ¼ 1 for odd l þ m. The leading term in Eq. ( 8), h ðN;ϵÞ lm is the Newtonian contribution
where Y l-ϵ;-m ðθ; ϕÞ are the scalar spherical harmonics, n ðϵÞ lm and c lþϵ ðνÞ are functions given in Eqs. ( 28) and (29) of Ref. [49], and v Ω is given by
Note that the SEOBNRv5HM model employs the variable
instead of v Ω in Eq. ( 9), which simplifies to v φ ¼ Ω½r 3=2 þ a 2=3 in the TM limit. However, the generalization of Eq. ( 11) to eccentric orbits is not straightforward due to the p r ¼ 0 requirement, which would affect the PN expansion of the eccentric modes unless a correction factor is applied. Thus, by choosing v Ω instead of v φ , one obtains a prescription that is easily applicable in generic orbits.
Furthermore, this choice does not have a significant impact on the accuracy of the underlying QC model SEOBNRv5HM [65]. Therefore, the eccentric 3PN RR force in Eqs. (5c) and (5d) is computed by employing v Ω . The function ŜðϵÞ eff is the effective source term which is given by
The factor T lm resums the leading order logarithms of tail effects and corresponds to
where Γ is the Euler gamma function, k ¼ mΩ in the TM limit and r 0 ¼ 2= ffiffi ffi e p . The remaining part of the factorized modes ( 8) is expressed as an amplitude f lm and a phase δ lm , which are computed such that the expansion of h lm agrees with the PN expanded modes. We point the reader to Appendix B of Ref. [49] for the explicit expressions of the different f lm and δ lm terms.
The TM limit of Eqs. ( 7) and ( 8) is obtained by setting the mass ratio ν to zero in the expressions, except for the leading ν in the Newtonian prefactor (9).
In our analysis, we consider equatorial planar geodesics of the Kerr background; hence, we consider F r ¼ F φ ¼ 0 in Eqs. (5c) and (5d), when evolving the dynamics. We characterize the planar orbits through the parameters fp; e; ag, which correspond to the semilatus rectum, the eccentricity and the spin of the Kerr BH, respectively. We adopt the Keplerian parametrization, where for the definitions of p and e we have:
where r a and r p are the radial separation at the apocenter and at the pericenter, respectively. As we show in Sec. II C, after evolving Eq. ( 5) without RR forces (i.e., for geodesics), we compute the analytical fluxes by evaluating the RR force in Eq. ( 6) on the geodesics. Hence, to avoid any possible confusion to the reader we stress that whenever the QC expressions in Eqs. (7) and (8) are employed, they are evaluated on the geodesic although they are quantities constructed assuming QC trajectories.
The core of our analysis relies on the computation of the energy flux Φ E and angular-momentum flux Φ J radiated by the TM to future null infinity. These fluxes are computed numerically by solving the Teukolsky master equation [80], which in Boyer-Lindquist coordinates reads -
This equation describes the evolution of scalar, vector, and tensor perturbations of a Kerr BH. The function Δ, the spin parameter a and the coordinates correspond to the same quantities defined in the previous sections. The parameter s is the spin weight of the field. In particular, when s ¼ AE2 the equation describes degrees of freedom of gravity that radiate, and for s ¼ -2 it is Ψ ¼ ðria cos θÞ 4 ψ 4 , where ψ 4 is the Weyl curvature scalar that describes outgoing GWs.
A system composed of a TM orbiting a Kerr BH can be interpreted as a perturbed Kerr metric, and within this interpretation the source term T in the right-hand side of Eq. ( 15) describes a TM moving in the Kerr spacetime. The details on how the source term T of the TM is constructed and how Eq. ( 15) is numerically solved are beyond the scope of this section. We mention the fact that the source term T of a TM orbiting a Kerr BH is constructed from Dirac-delta functions of the variables r and θ, as well as first and second derivatives of the delta functions in these variables. These terms are sourced at the location of the TM, hence the source T depends on the trajectory that the TM follows in the Kerr spacetime. The details can be found in Refs. [81,83,84]. In this analysis, the trajectories used to source the term T are the geodesics introduced at the end of Sec. II A, constructed evolving Eq. (5).
To solve Eq. ( 15), we make use of two different codes: when considering bound orbits we adopt the FD code of Ref. [81], while when considering unbound orbits we employ the TD code developed in Ref. [82][83][84]. At future null infinity, the Weyl scalar ψ 4 and the waveform strain h ¼ h þih × are related by the expression
and following standard practices, the waveform is decomposed in spin-weighted spherical harmonics
The fluxes Φ E and Φ J are computed from the modes through the expressions
In Appendix B, we give closed form expressions for constructing the instantaneous fluxes from FD Teukolsky solutions. In this work, we truncate the summation over the modes at l ¼ 8. More details on the numerical errors of the FD and TD Teukolsky codes employed in our analysis can be found in Refs. [77,[81][82][83][84].
In order to assess the EOB eccentric RR force of Eq. ( 6) by comparing with numerical results, it is necessary to compute the analytical fluxes. The connection between the RR force and the fluxes is given by the balance equations
These equations relate the time-dependent fluxes at future null infinity, Φ E and Φ J , to the change in the energy and angular momentum of the system, Ėsystem and Jsystem , together with two other terms that appear as total time derivatives, ĖSchott and JSchott , known as Schott terms. These two terms take into account the contributions to the fluxes due to the interaction of the system with the radiation field, as originally pointed out in the context of electromagnetism in Ref. [87] and they were introduced in the context of the EOB framework in Ref. [50].
The connection between the fluxes Φ E=J , the RR force F and the Schott terms is made explicit by first considering the Hamilton equations (5), which lead to
where the fact that the Hamiltonian in Eq. ( 2) does not depend on the azimuthal angle φ is exploited. By plugging these expressions in Eq. ( 19), one gets
From these equations we are able to compute the fluxes from the RR force in Eq. ( 6), providing a way to test its prescriptions. This is done by first computing Kerr geodesics as explained in Sec. II A and then evaluating all the quantities involved in the right-hand side of Eq. ( 21) on the geodesics. We point the reader to Appendix A for the expressions of the PN time derivatives of the Schott terms, ĖSchott and JSchott , as functions of the EOB dynamical variables, as used in this work.In these expressions, which are written for the first time in this work, we follow the gauge choice employed in Ref. [64] by fixing the gauge constants as α ¼ -16=3 and β ¼ -13=2.
Note that in the above, we have not included any effects due to the central Kerr BH absorbing GWs, changing its mass and spin. If included, these effects would alter the relationship between the RR forces F r=φ and fluxes Φ E=J in Eq. (21). In practice, the effects of absorption are typically several orders of magnitude smaller than the fluxes to infinity see e.g., [88] (but can become order 10% for extreme orbits close to the horizon of a nearly extremal BH [78,89,90]). As such, the absorption fluxes are mostly relevant when their effects can accumulate over a large number of orbits in an inspiral. For the single orbit comparisons in this work, their impact would be minimal. Nonetheless, to ensure an apples-to-apples comparison, we also include only the fluxes to infinity in the numerical Teukolsky fluxes.
In the following, we present the comparison between the analytical and numerical fluxes computed as explained in Sec. II. In particular, we focus on studying the different PN contributions in the eccentric part of the analytical fluxes.
Throughout this section, we label the fluxes that include PN corrections due to eccentricity at 1PN, 2PN and 3PN order (including SO and SS corrections up to 2PN order), as "ECC1PN," "ECC2PN" and "ECC3PN," while the fluxes that simply use the QC prescription of the fluxes evaluated on an eccentric trajectory are labeled "QC." Note that the fluxes labeled "ECCnPN" still contain all available PN orders in the QC part of the flux.
We start by comparing the numerical and analytical instantaneous fluxes. Before doing so, we emphasize that the Schott terms in Eq. ( 21), are known only as PN expansions up to 3PN order, while the RR force terms have been resummed in the EOB formalism. Consequently, there is a limit to what can be learned when comparing the instantaneous fluxes to the numerical-Teukolsky fluxes in the strong-field regime, as it is unclear whether any particular disagreement is due to an inaccuracy of the RR forces or a possible degradation of the PN approximation in the Schott terms. This ambiguity is not present when considering orbit-averaged fluxes, as we mention in Sec. III B. Nonetheless, we believe it is instructive to look at the comparison of the instantaneous fluxes as this helps build an intuitive picture of our analysis.
The comparison of the numerical and analytical fluxes is illustrated in Fig. 1. We consider an orbit with semilatus rectum p ¼ 13, eccentricity e ¼ 0.5 and spin a ¼ 0. The analytical and numerical fluxes are plotted over a radial period on the bound geodesic. In the left panel, the orbit is shown, while in the top panels the instantaneous fluxes between two consecutive apocenters are plotted. The black curves correspond to the numerical fluxes obtained from FIG. 1. Fluxes for an eccentric geodesic with p ¼ 13, e ¼ 0.5 and a ¼ 0. On the left panel, the planar orbit is shown, highlighting in blue the trajectory over one radial period. The top panels show the angular momentum Φ J and energy Φ E fluxes at infinity: black curves are the Teukolsky fluxes, while colored curves are the analytical fluxes. "QC" refer to the QC fluxes while "ECCnPN" refer to fluxes computed considering eccentric corrections at nPN order. The bottom panels show the relative difference with respect to the Teukolsky fluxes.
the FD Teukolsky code, while the colored lines are the analytical fluxes. We show the QC fluxes (blue curves), which are computed with the QC RR force in Eq. ( 7), and the eccentric fluxes ECC1PN (red curves), ECC2PN (orange curves) and ECC3PN (green curves). Finally, the bottom panels of Fig. 1 show the fractional difference between the analytic and Teukolsky fluxes. The analytical fluxes are computed considering the multiplicative implementation (6). The additive corrections provide similar results to Fig. 1, and are not shown here.
For the geodesic considered in Fig. 1, we do not find a clear improvement of the fluxes at different PN orders along the orbit. Around the apocenter passage (for t ≈ 600 and t ≈ 1100 where r ≈ 25) we observe that the ECC1PN fluxes better approximate the numerical ones, whereas the ECC2PN and ECC3PN fluxes are close, but do not improve the approximation. By removing the eccentric PN tail terms, we find that the 1.5PN tail contribution is degrading the accuracy of the instantaneous results. However, in Sec. III B, we will discuss that these terms are essential to recover the correct scaling in the weak-gravitational field regime when considering the orbit-averaged fluxes. This confirms what we anticipated at the beginning of this section: one PN order may be worse for the instantaneous fluxes, but better for the averaged fluxes, and thus, we cannot conclude from those results if the RR force is accurate since we cannot disentangle it from the effect of the PN-expanded Schott terms.
Near the pericenter passage (for t ≈ 850 where r ≈ 9) we find that the ECC3PN energy flux improves the other PN orders, while for the angular-momentum flux it is the QC curve that best approximates the numerical flux. This close agreement of the QC fluxes near the pericenter passage was already pointed out in Ref. [68] for some orbital configurations, and explained as a numerical coincidence. This will become more apparent in Sec. III B, as we show the behavior of this agreement as a function of the orbital separation.
The bottom panels of Fig. 1 show an asymmetry of the relative differences with respect to the pericenter passage, which was also observed in Ref. [58]. After convincing ourselves that this behavior is not due to any numerical artifact, we find that it comes from an asymmetry of the Teukolsky fluxes with respect to the pericenter. We conclude that this asymmetry arises from contributions that are not modeled by the EOB fluxes, but that are present in the numerical fluxes. By inspecting different orbits and observing the same pattern, especially for orbits that lay in the weak-field regime, we suggest that this asymmetry comes from delayed contributions coming from higher-order-PN tail terms that are not present in the analytical fluxes.
We notice that there are qualitative differences between our results and Fig. 1 of Ref. [68], despite the two figures ostensibly plotting the same quantities for the same orbit with both the QC and ECC2PN (labeled "QC2PN" in Ref. [68]); in our version, the fluxes are much closer to the numerical values. The discrepancy in the QC flux may arise due to the proxy used in Ref. [68] for the QC flux not faithfully reproducing the QC flux of SEOBNRv4HM (and notably the QC flux of SEOBNRv5HM, which is used in this paper) when applied to eccentric orbits. The discrepancy between the ECC2PN fluxes, on the other hand, is due to the fact that in Ref. [68], the authors use the Schott terms of Ref. [50] when computing the fluxes through Eq. ( 21). However, this is not compatible with the gauge of the eccentric RR force of Ref. [60], which is used to compute the ECC2PN (QC2PN) fluxes. There is no freedom in using Schott terms in a different gauge, and the ones of Ref. [60] must be considered when computing the instantaneous ECC2PN fluxes.
Figure 1 indicated that the 1PN eccentric corrections of the multiplicative implementation of the RR force (6) better approximate the numerical flux than the higher order eccentric PN corrections near the pericenter of the orbit. The additive implementation (6b) (not plotted in Fig. 1) shows similar behavior. This motivates a further investigation, because, in principle, one would expect the higher order PN corrections to improve the approximation, at least in the weak-field regime.
To better assess this, we start by restricting to the Schwarzschild case and we analyze what happens when considering the eccentric force in the weak-field scenario. For the study, we consider the averaged fluxes over one radial period T r , given by:
The averaging eliminates the ambiguity due to the Schott terms contribution in the instantaneous-fluxes expressions (21), as the Schott terms appear as total time derivatives. We consider a set of bound orbits defined by eccentricity e ¼ 0.5, spin a ¼ 0 and decreasing semilatera recta p ¼ f960; 480; 240; 120; 60; 30; 15g. This choice allows us to explore the weak-field regime in order to check whether the flux residuals, obtained by subtracting the eccentric corrections (6a) and (6b) from the numerical fluxes, have the expected PN scaling.
In Fig. 2, we show the fractional difference of the averaged numerical fluxes with respect to the averaged analytical fluxes
computed on these orbits for different values of the gaugeinvariant quantity x defined by
For the different PN corrections, we follow the same nomenclature used in Fig. 1, and we consider both implementations: solid lines correspond to the multiplicative RR force from Eq. (6a) while dotted lines correspond to the additive RR force from Eq. (6b). We observe that the eccentric corrections to the RR force provide a consistent improvement over the QC force. For orbits that are in the weakest regimes, x ≈ 8 × 10 -4 , including the ECC3PN corrections improves the agreement with the numerical flux by a factor 10 8 over the QC flux. We also find that the curves follow the general expected scaling at low x. The fractional difference between the QC analytical and the Teukolsky flux approaches a constant in the weak field regime, indicating that the QC prescription already needs corrections at the leading "Newtonian" order. The slopes of fractional residuals, after subtracting the ECCnPN corrections in the log-log plot of Fig. 2, are compatible with the expected x -ðnþ1=2Þ behavior of a residual that starts at ðn þ 1=2ÞPN order.
Since for the instantaneous flux comparison in Sec. III A the inclusion of the PN tail corrections did not lead to an unequivocal improvement of the flux, we also considered the effects of omitting the PN tail terms on the orbitaveraged flux. The gray lines in Fig. 2 correspond to the ECC3PN fluxes computed without the tail part, showing that the inclusion of these terms is essential in order to recover the expected PN scaling, corroborating the hypothesis that degradations of higher-order terms in the instantaneous fluxes may come from the Schott terms and disappear when orbit averaging.
When moving into the stronger-field regimes (i.e., for higher values of x), the discrepancy between analytical and numerical averaged fluxes increases, as expected due to the PN nature of the analytical fluxes.
In order to assess the RR force over a wider region of the parameter space, we consider more orbits with different eccentricities spanning the milder and stronger field regimes. Figure 3 shows the fractional differences of the fluxes evaluated at geodesics with eccentricities e ¼ f0.1; 0.3; 0.5; 0.7g, spin a ¼ 0 and semilatera recta p ¼ fp LSO þ 0.025;p LSO þ 0.05;p LSO þ 0.1;7;8;9;10;13; 15g,foot_0 where p LSO corresponds to the semilatus rectum of the last stable orbit (LSO) with the corresponding e. For Schwarzschild geodesics this is given by
The strong-field residuals in Fig. 3 show a much less organized picture than their weak-field counterparts in Fig. 2. There is not always a clear order-by-order improvement from adding higher PN eccentric corrections. This signifies (the start of) the breakdown of the convergence of the PN series in this regime. Notably, the QC and ECC1PN fluxes seemingly outperform the higher-order corrections for larger x. In fact-contrary to naive expectation-the residuals from the QC and ECC1PN fluxes actually decrease, for both the angular-momentum and energy fluxes, up to certain values of x. This behavior is evident starting from values of x ≈ 0.06 for the orbits with e ¼ f0.1; 0.3; 0.5g and of x ≈ 0.04 for the ones with e ¼ 0.7. This decreasing trend is interrupted at larger values of x, where the QC and ECC1PN residuals of the fluxes start to FIG. 2. Absolute value of the fractional differences defined in Eq. ( 23) of the averaged fluxes computed for orbits in the weak-field regime. The x-axis corresponds to the gauge-invariant variable x ¼ h φi 2=3 . The orbits are characterized by the parameters e ¼ 0.5, a ¼ 0 and p ¼ f960; 480; 240; 120; 60; 30; 15g. Dotted lines correspond to the additive implementation (add) while solid ones correspond to the multiplicative (mult). Empty circle dots represent negative values of the fractional differences, highlighting zero crossing in the log-log plots.
monotonically increase up to the closest orbits to the LSO with p ¼ p LSO þ 0.025. We find that this is connected to a change in sign of the fractional differences between the analytical/numerical fluxes. Since the absolute fractional differences are plotted in Fig. 3, here we represent this change of sign by using different dots: filled dots for positive values of the relative difference (23), empty dots for negative values. This change of sign signifies that in the weak field the average QC fluxes overestimate the average numerical flux, while they underestimate it in the stronger field. Consequently, there are some "goldilocks" configurations for which the QC (and ECC1PN) fluxes are "just right," and by coincidence produce the correct flux, corroborating what was found in Ref. [69] for the QC fluxes.
We do not observe a significant difference between the multiplicative and additive implementations in the nonspinning case. However, in Sec. III C we show that this changes when considering spin.
We also find that larger eccentricity values impact the accuracy of the analytical averaged fluxes. More specifically, we observe that for similar values of the x parameter, the relative differences are more prominent for higher eccentricities. This trend is due to the fact that the pericenter of the orbits-where most radiation is generated-is pushed more and more in the strong-field regime, where in turn, the PN expressions are less reliable.
Figure 3 allows us to gauge the overall performance of the flux approximations: in the worst case scenario (i.e., for orbits with e ¼ 0.7 and close to the LSO), we observe a difference of ≈5%, when considering the ECC3PN fluxes.
So far, we have shown a comparison between analytical and numerical fluxes for eccentric (bound) orbits of a TM moving on geodesics around a Schwarzschild BH. In this section, we provide results that assess the RR force corrections in Eq. ( 6) when we allow the central BH to have a spin aligned or antialigned (henceforth, for simplicity aligned) with the orbital angular momentum (i.e., we consider eccentric-equatorial orbits around a Kerr BH).
To test the spinning case, we consider orbits with a fixed value of the x parameter defined in Eq. ( 24). This choice allows us to identify in a gauge-invariant manner different gravitational-field regimes. The orbits we analyze have eccentricity values e ¼ f0.0; 0.1; 0.3; 0.5; 0.7g and spins of the central BH a ¼ f-0.99; -0.9; -0.7; -0.5; 0; 0.5; 0.7; 0.9; 0.99g.
In Fig. 4, we show contour plots of the fractional difference of the averaged fluxes (22) with respect to the eccentricity and the spin values for orbits with fixed x ¼ 0.06. We choose this value for x since it allows computing stable geodesics all over the subset of the parameter space fe; ag we want to investigate. We consider the 3PN order in the eccentric sector, ECC3PN, for both implementations (6a) and (6b). FIG. 3. Absolute value of the fractional differences of the averaged fluxes computed on orbits in mild/strong-field regimes. The orbits with highest x (last points on the right) are characterized by a semilatus rectum p ¼ p LSO þ 0.025. The nomenclature is the same as in Fig. 2.
We observe that, regardless of the spin values, for small eccentricities (e ≤ 0.3), the fractional difference is less than 0.5% for both resummations. When considering higher values of e, the differences increase. This trend corroborates what we mentioned at the end of Sec. III B, notably, for fixed values of x, larger eccentricities degrade the accuracy of the analytical-averaged fluxes, because the pericenter of the orbit is pushed more in the strong-field regime. This degradation increases when considering prograde orbits with high-spin values, for which the pericenter is even closer to the central BH. For these orbits, we find that the multiplicative implementation provides a mismatch of the fluxes less than 5%, for the considered portion of parameter space. In comparison, the additive implementation provides relative differences that are less than 15%. In particular, we find that for more extreme orbits (i.e., for the ones with high eccentricity and spin), the multiplicative implementation improves over the additive one by a 10% difference.
As a final remark, we mention an overall general improvement of the multiplicative ECC3PN factorization with respect to the additive ECC3PN factorization for regimes with x ≤ 0.06 over the parameter space ðe; aÞ for both the energy and angular-momentum fluxes. This can be seen in Fig. 7 in Appendix A 4 where we show the behavior of the fluxes at different PN orders, for both factorizations, from weak-field regimes up to x ¼ 0.06. The considered orbits use the particular configuration of eccentricity and spin of Fig. 4, e ¼ 0.5 and prograde spin a ¼ 0.9, with semilatera recta p ¼ f960; 480; 240; 120; 60; 30; 15g. Figure 7 also highlights (as in Fig. 2) the convergence of the different PN truncations of the eccentric fluxes pointing out the correctness of the procedure (and of the expressions) used to derive the EOB eccentric RR force in Eq. ( 6) in the general spinning scenario. We performed this same test for other eccentricity and spin configurations ðe; aÞ of Fig. 4, finding similar results. FIG. 4. Contour plots of the absolute value of the fractional differences of the ECC3PN averaged fluxes with respect to the numerical one. They are computed on orbits that span the parameter space (e, a) with fixed value x ¼ 0.06. The upper panels show the multiplicative implementation while the lower ones show the additive implementation. The multiplicative implementation shows a worst case scenario of 5% relative difference for higher eccentric spinning prograde orbits.
Finally, we push the comparison of the eccentric analytical and numerical fluxes to even higher eccentricities to the e > 1 regime (i.e., hyperbolic-scattering orbits). In this section, we restrict our attention to nonspinning Schwarzschild BHs.
We consider hyperbolic orbits with fixed energy E 0 ¼ 1.005, and we vary the angular momentum J 0 . Figure 5 shows the five different orbits we consider. To produce them, we examine different values of the angular momenta by first computing the critical value J c through Eq. ( 6) of Ref. [91]. This critical angular momentum represents the smallest value a TM with fixed energy must have to still scatter back to infinity without plunging into the central BH. For the chosen value of E 0 , it is J c ¼ 4.0397. We consider different values for J 0 ¼ fJ c þ 2; J c þ 0.7; J c þ 0.2; J c þ 0.07; J c þ 0.02g, which provide orbits with increasing scattering angle and decreasing pericenter distance. As in the bound case, we compute the analytical and numerical fluxes on these orbits. The numerical fluxes are obtained through the TD Teukolsky code of Refs. [82][83][84] as mentioned in Sec. II B, while the analytical fluxes are computed similarly to the bound-orbits case (i.e., by evaluating Eq. ( 21) on the hyperbolic geodesics). In the following, we consider only the multiplicative implementation (6), since in Fig. 4 we observe it improves over the additive one.
In Table I, we compare the analytical/numerical total energy and angular momentum emitted by the TM on the examined orbits. These quantities are evaluated by integrating the fluxes on each orbit and truncating at different PN orders. As in the case of the averaged fluxes for bound orbits (see Sec. III B), we also find that the eccentric corrections improve the agreement with the numerical fluxes with respect to the QC case. We observe an improvement of the PN series in the eccentric sector, especially for orbits with smaller J 0 . In almost all the examined cases, the 3PN corrections provide the best approximation of the total radiated energy and angular momentum with respect to the numerical flux. However, we find that the fractional differences for the orbits with smaller J 0 (J 0 ¼ 4.0597 and J 0 ¼ 4.1097) are greater than 10%. This is expected, because these orbits have angular momentum close to the critical one J c and are characterized by strong-field regimes, where the PN expansions start to breakdown. For the orbits with larger angular momentum, characterized by weaker-field regimes, the residuals are smaller (≤ 8%) but we do not observe a large improvement of the PN series when increasing the order in general. Moreover, for the trajectory with higher angular momentum (J 0 ¼ 6.0397), we observe a slight deterioration of the PN series, with the eccentric 1PN corrections providing the best approximation. We believe this is because even in FIG. 5. The hyperbolic orbits of a TM around a Schwarzschild BH considered in our analysis. They all have energy E 0 ¼ 1.005 and critical angular momentum J c ¼ 4.0397. We consider five orbits which have total angular momenta
Comparison of the analytical/numerical total emitted energy and angular momentum for the different hyperbolic orbits shown in Fig. 5. The analytical fluxes are computed with the multiplicative implementation (6). For each orbit with energy E 0 and angular momentum J 0 , we show the values of the pericenter distance r p , the numerical total emitted energy E Teuk and angularmomentum J Teuk , and the fractional differences between the analytical emitted energy/angular-momentum, truncated at different PN orders in the eccentric sector, and the total numerical fluxes.
.005 4.0597 4.39 2.605 × 10 -1 0.3304 0.2095 0.2122 0.1330 2.805 0.2121 0.1948 0.1637 0.0998 1.005 4.1097 4.83 1.450 × 10 -1 0.3316 0.1571 0.1650 0.0964 1.818 0.1701 0.1433 0.1135 0.0593 1.005 4.2397 5.64 6.680 × 10 -2 0.3523 0.1056 0.1204 0.0807 1.062 0.1318 0.0927 0.0711 0.0412 1.005 4.7397 8.14 1.381 × 10 -2 0.4282 0.0473 0.0673 0.0631 3.836 × 10 -1 0.0964 0.0394 0.0357 0.0370 1.005 6.0397 14.67 1.477 × 10 -3 0.5225 0.0065 0.0093 0.0103 9.761 × 10 -2 0.0911 0.0159 0.0224 0.0280 the far weak field these orbits still have high (relativistic) velocities.
To provide an intuitive picture of this last point we consider again the instantaneous fluxes. Figure 6 shows the time-dependent fluxes for the hyperbolic encounters with the highest (J 0 ¼ 4.0597) and smallest (J 0 ¼ 6.0397) scattering angle. From the plots, we observe that for the orbit with the highest scattering angle, in the mild/ weak-field regime part of the orbits (when r > 20), the analytical fluxes differ by > 80% from the numerical ones. This is particularly evident for the ECC2PN and ECC3PN fluxes, for which the fractional difference is > 100%. We observe a similar pattern for the orbit with the smallest scattering angle (J 0 ¼ 6.0397) (i.e., in the weak-field regime), a large mismatch between the numerical and the analytical curves arises. This large fractional difference we observe for r > 20 is caused by the fact that, for hyperbolic orbits, considering the weak-gravitational-field regime does not necessarily correspond to considering small velocities (v ≪ 1) of the TM. This large-velocity regime deteriorates the convergence of the PN expansion, as shown in the last line of Table I.
In this work, we tested an EOB eccentric RR force computed with the same procedure as Ref. [60] with a different gauge choice for the leading PN order. The nonspinning 3PN part of this RR force is derived in Ref. [64] while the SO and SS contributions, at 1.5PN and 2PN respectively, are similar to the ones of Ref. [60] but computed taking into account the different gauge choice of the leading order. We considered two possible resummations of this force: a multiplicative (6a) and an additive (6b) one. The assessment of the resummations is performed by considering a TM orbiting the equatorial plane of a Kerr BH. We computed the analytical energy and angular-momentum fluxes, which are linked to the RR FIG. 6. Instantaneous fluxes for hyperbolic orbits with E 0 ¼ 1.005 and J 0 ¼ f4.0597; 6.0397g. Black curves are the numerical fluxes, blue-dotted ones represent QC fluxes, red curves are the fluxes containing 1PN eccentric corrections while orange and green curves are the eccentric 2PN and 3PN corrected fluxes. The portion of orbit we consider is highlighted in light-blue in the orbit-plot and it is for r < 30.
force through the balance equations (21), on eccentric geodesics of the Kerr metric, and we compared them with the numerical fluxes computed through the use of a FD [81] and a TD [82][83][84] Teukolsky codes. We focused our analysis on the orbit-averaged fluxes for bound orbits and on the emitted energy and angular momentum for hyperbolic orbits. This is because when integrating the fluxes over an orbit, we are sure that any degradation coming from the PN-expanded Schott terms is integrated out and we can test the RR force resummations.
When considering bound orbits in the weak-field limit, we recovered the expected scaling of the different PN truncations of the fractional difference of the averaged fluxes, indicating that the procedure used to compute the RR force and its resummed expressions are both correct. We also showed that the eccentric corrections to the RR force are necessary to improve the QC fluxes in this regime. This improvement over the QC fluxes is of the order of 10 8 when considering the ECC3PN fluxes for orbits with x ≈ 10 -3 . For nonspinning bound orbits in stronger-field regimes, computed up to the LSO, we found that the fractional differences of the averaged fluxes are < 5% for eccentricity up to e ¼ 0.7. This result holds for both the multiplicative and additive resummations.
For the more general case of bound equatorial geodesics of a Kerr BH, characterized by x ≤ 0.06, e ¼ ½0; 0.7 and a ¼ ½-0.99; 0.99, we found that for e < 0.3 the fluxes discrepancies are < 0.5% for both resummations. For larger eccentricities up to e ¼ 0.7, we showed that the multiplicative implementation provides discrepancies with the numerical fluxes that are < 5%, while for the additive implementation, we found larger fractional differences, up to ≈15% values. This indicates that, in the considered part of the parameter space, the multiplicative resummation ( 6) is a better approximation of the fluxes up to regimes with x ¼ 0.06. This is especially true for highly eccentric (e ¼ 0.7) and prograde orbits with large Kerr spin (a ¼ 0.99).
Finally, we concluded our analysis by considering geodesic hyperbolic orbits with fixed energy E ¼ 1.005 and different angular momenta. We found that the 3PN eccentric corrections of the multiplicative resummation improve the total emitted energy and angular momentum for hyperbolic encounters with angular momenta closer to the critical one. However, for these orbits, the relative differences at 3PN are larger than 8%, due to the strongerfield regimes characterizing them. For the unbound orbits in weaker-field regimes, we found smaller values of the fractional differences (< 7%) of the emitted energy and angular momentum. However, we observed a deterioration of the convergence of the PN series, with the 3PN corrections providing higher discrepancies. We pointed out this may be due to the fact that for unbound orbits the weak-gravitational-field regime does not necessarily imply small velocities, impacting the convergence of the eccentric PN orders.
The results we found for the unbound case indicate that the RR force should be further improved for large eccentricities and high velocities. There are possible different strategies to tackle this last point: one possibility is to determine whether different parametrizations improve the results we find. Another strategy could be to perform a further resummation of the eccentric corrections to the RR force F ecc r=φ in Eq. (6). A different approach would be to numerically inform the eccentric corrections by introducing terms that fit the numerical fluxes in the parameter space. We believe all these strategies are promising for improving the EOB eccentric fluxes we studied in this work, and we leave them to future works.
More specifically, we employ the complete 3PN expressions of the EOB fluxes presented in Ref. [64]. Then, the free unknown coefficients (gauge constants) appearing in the Ansätze for the Schott terms are determined by imposing the regularity of the radial component of the RR force in the QC limit (this means no 1=p r terms in F r ), and by choosing a specific gauge. Following Refs. [60,64] we employ a gauge that satisfies
With this choice, the QC limit of the resulting RR force satisfies the gauge used in SEOBNRv5HM, which is given by Eq. ( 7). At leading PN order, the components of the RR force, in a generic gauge, are given by
where p 2 ¼ p 2 r þ p 2 φ =r 2 , and fα; βg are gauge constants, representing the gauge freedom in defining the RR force [93,94]. 3 In particular, the gauge in Eq. (A2) fixes the value of α to be -16=3, but leaves β unspecified. As shown in Ref. [64], β needs to equal -13=2 to avoid a 2.5PN modification (relative to the leading order) of the QC orbital phase when transforming between harmonic and EOB coordinates. For the expressions of this work, we use the aforementioned values of α and β, to be consistent with Ref. [64]. As a consequence of this, the expressions used here do not follow the RR gauge employed in Ref. [60].
The next step consists in factorizing the RR force in terms of a QC part and an eccentric correction. In this way, given the PN expressions of F φ and F r , and the QC part given by Eq. ( 7), we determine the different PN expressions for F ecc φ and F ecc r which are specified in Eq. ( 6). The results valid for the TM limit are explicitly shown in Appendixes A 1 and A 2. Additionally, the expressions for the time derivatives of the Schott terms in the TM limit are shown in Appendix A 3.
In the following subsections, we write the expressions in terms of the variables fr; p r à ; v 0 g, where r is the radial separation, p r à is the radial momentum conjugate to the tortoise coordinate r à , and v 0 is given by
The expressions of the eccentric corrections for the multiplicative implementation ( 6) are given by
where the expressions of the different PN orders are
3 Different values for the gauge constants α and β have been considered in the literature. For example: the choice α ¼ -1, β ¼ 0 corresponds to the harmonic RR gauge [95]; the values α ¼ 5=3, β ¼ 3 correspond to the Arnowitt-Deser-Misner (ADM) RR gauge [96,97]; and the values α ¼ 0, β ¼ 2 correspond to a choice employed in Ref. [50] for a "minimal decomposition" of the energy and angular momentum fluxes.
TESTING ECCENTRIC CORRECTIONS TO THE RADIATION-… PHYS. REV. D 111, 044036 (2025) 044036-15
r à 24192r 5 v 10 0 -868811p 4 r à 12096r 6 v 10 0 -3011119p 2 r à 13608r 7 v 10 0 -804565 6804r 8 v 10 0 -52459p 6 r à 217728r 7 v 14 0 -63503p 4 r à 48384r 8 v 14 0 -85591p 2 r à 36288r 9 v 14 0 -151855 108864r 10 v 14 0 þ 2465p 8 r à 18144r 7 v 16 0 -930955p 6 r à 217728r 8 v 16 0 -1604525p 4 r à 72576r 9 v 16 0 -102055p 2 r à 3024r 10 v 16 0 -840635 54432r 11 v 16 0 -95p 8 r à 972r 10 v 22 0 -1415p 6 r à 1944r 11 v 22 0 -655p 4 r à 324r 12 v 22 0 -1205p 2 r à 486r 13 v 22 0 -275 243r 14 v 22 0 ;
where γ ≈ 0.577 is the Euler-gamma constant. While for the azimuthal component of the RR force we have
with
123157v 2 0 ln 2 630r 2 -80 9 π 2 rv 8 0 -856 63 γrv 8 0 -86801rv 8 0 4410 -856 21 rv 8 0 ln v 0 -856 ln 2 21 rv 8 0 :
Here we provide the additive implementation (6b) expressions. For the radial component we have
where the different terms are
TESTING ECCENTRIC CORRECTIONS TO THE RADIATION-… PHYS. REV. D 111, 044036 (2025)
The azimuthal component is
while all the different PN orders are given by
5r 5 -8 15 r 4 v 18 0 -224v 4 0 15r 3 -4v 6 0 15r 2 -12 5 rv 12 0 -256v 10 0 15 ;
Here we provide the expressions of the total time derivative of the Schott terms introduced in Eqs. (19a) and (19b). They are expressed as functions of the variables fr; p r à ; v 0 g. For the Schott term contribution to the energy balance equation we have
where
10395r 8 -1483r 7 v 30 0 6930 þ 7537561v 6 0 38500r 5 -142096v 6 0 ln r 1575r 5 þ 1417r 4 v 24 0 770 -1133672317v 12 0 2079000r 2 þ 29104v 12 0 ln r 225r 2 þ 182383141rv 18 0 693000 -6848 175 rv 18 0 ln r;
7 -247384p 2 r à ln r 175r 7 þ 468018p 2 r à ln 3 175r 7 -409168p 2 r à ln 2 525r 7 -6208π 2 p 2 r à v 6 0 45r 4 þ 664256γp 2 r à v 6 0 1575r 4 þ 19507708p 2 r à v 6 0 55125r 4 -332128p 2 r à v 6 0 ln r 525r 4 þ 624024p 2 r à v 6 0 ln 3 175r 4 -609472p 2 r à v 6 0 ln 2 225r 4 -512π 2 15r 8 þ 54784γ 525r 8 þ 3736352 18375r 8 -27392 ln r 175r 8 þ 109568 ln 2 525r 8 -16π 2 v 6 0 45r 5 þ 1712γv 6 0 1575r 5 -8511719v 6 0 55125r 5 -856v 6 0 ln r 525r 5 þ 156006v 6 0 ln 3 175r 5 -279056v 6 0 ln 2 315r 5 þ 1552π 2 v 12 0 45r 2 -166064γv 12 0 1575r 2 -2697337v 12 0 55125r 2 þ 83032v 12 0 ln r 525r 2 -156006v 12 0 ln 3 175r 2 þ 152368v 12 0 ln 2 225r 2
; ðA14iÞ while for the Schott terms related to the angular momentum flux we have
with J0PN Sch ¼ 256p 2 r à v 3 0 15r þ 128v 3 0 15r 2 -128rv 9 0 15 ; ðA16aÞ J1PN Sch ¼ -64p 4 r à 15r 4 v 3 0 -64p 4 r à v 3 0 15r -32p 2 r à 3r 5 v 3 0 -6227p 2 r à v 3 0 105r 2 -32 5 p 2 r à rv 9 0 -64 15r 6 v 3 0 -32 15 r 3 v 15 0 -219v 3 0 35r 3 þ 443v 9 0 35 ; ðA16bÞ J1.5PN;Tail Sch ¼ -10πp 2 r à r 4 -10π 3r 5 þ 10πv 6 0 3r 2 ; ðA16cÞ J1.5PN;SO Sch ¼ a -16p 4 r à r 3 -448p 2 r à 5r 4 -24p 2 r à v 6 0 5r -152 5r 5 þ 116v 6 0 5r 2 þ 36 5 rv 12 0 ;
83032p 2 r à v 3 0 ln r 105r 4 þ 156006p 2 r à v 3 0 ln 3 35r 4 -1552π 2 v 3 0 45r 5 -152368p 2 r à v 3 0 ln 2 45r 4 þ 166064γv 3 0 1575r 5 þ 2697337v 3 0 55125r 5 -83032v 3 0 ln r 525r 5 þ 156006v 3 0 ln 3 175r 5 -152368v 3 0 ln 2 225r 5 þ 1552π 2 v 9 0 45r 2 -166064γv 9 0 1575r 2 -2697337v 9 0 55125r 2 þ 83032v 9 0 ln r 525r 2 -156006v 9 0 ln 3 175r 2 þ 152368v 9 0 ln 2 225r 2 :
We remark that these total time derivatives of the Schott terms have to be computed on the geodesics and plugged in the balance equations.
In Fig. 7, we show a similar plot to the one in Fig. 2 for eccentric orbits with e ¼ 0.5 around a Kerr BH with spin a ¼ 0.9. We test the weak-field regime scaling of the different PN order of the eccentric corrections to the RR force in Eq. ( 6) by studying the fractional differences of the orbit-averaged analytical and numerical fluxes. As in the nonspinning case, we recover the correct scaling of the different PN truncations, which is expected to be x -ðnþ1=2Þ , with n being the PN order at which we perform the truncation. However, in the spinning case, our analysis has an important caveat: the analytical fluxes we consider contain SO and SS spin corrections at 2PN. This means that the ECC3PN fluxes do not contain the spin contributions at 2.5PN and 3PN. As a consequence the ECC3PN lines in Fig. 7 exhibit the same scaling (slope) as the ECC2PN, but they still improve the ECC2PN lines exhibiting a shift down. In the plot we highlight with a vertical gray curve the orbit with x ¼ 0.06 which is present in the contour plot of Fig. 4. This plot extends to regimes x ≤ 0.06 what we found in Fig. 4: the multiplicative ECC3PN corrections in Eq. ( 6) perform better than the additive ECC3PN in Eq. ( 6), especially for eccentric orbits with prograde spin.
The frequency domain solutions to the Teukolsky equation at future null infinity for a particle traveling on an eccentric geodesic take the form ψ 4 ¼ X lmn Z lmn -2 S lmω mn ðcos θÞe imϕ e -iω mn t ; ðB1Þ
where the -2 S lmω mn spin-weighted spheroidal harmonics of weight -2, and the mode frequencies are ω mn ¼ mΩ ϕ þ nΩ r with Ω ϕ and Ω r the azimuthal and radial frequencies. To build the instantaneous fluxes using Eq. ( 18) we first project the spheroidal harmonics
where the coefficients ð -2 b mn Þ l l are obtained using the algorithm of Ref. [98]. This allows us to define the spinweighted spherical harmonic coefficients of ψ 4 ,
Combining the above expressions with Eq. ( 18), and applying the orthogonality relations for products of spinweighted harmonics when integrated over the sphere, we obtain the following closed-form expressions for the instantaneous flux,
ω mn ω mn 0 e -iðn-n 0 ÞΩ r t ; and
mn 0 e -iðn-n 0 ÞΩ r t : ðB4bÞ FIG. 7. Absolute value of the fractional differences defined in Eq. ( 23) of the averaged fluxes computed on orbits in a weak-field regime. The x-axis corresponds to the gauge-invariant variable x ¼ h φi 2=3 . The orbits are characterized by the parameters e ¼ 0.5, a ¼ 0.9 and p ¼ f960; 480; 240; 120; 60; 30; 15; 12.797g. Dotted lines correspond to the additive implementation (add) while solid lines correspond to the multiplicative (mult) implementation. Empty circle dots represent negative values of the fractional differences, highlighting zero crossing in the log-log plots. We highlight with a dashed vertical gray line the orbit with x ¼ 0.06, which is also shown in the contour plot of Fig. 4.
For the orbits with e ¼ f0.5; 0.7g, we consider the set p ¼ fp LSO þ 0.025; p LSO þ 0.05; p LSO þ 0.1; 8; 9; 10; 13; 15g, since p LSO ≥ 7 for these orbits. We also limit the eccentricity to ≤ 0.7 since frequency-domain Teukolsky codes become computationally expensive for high eccentricities.