Radiation-reaction force and multipolar waveforms for eccentric, spin-aligned binaries in the effective-one-body formalism

While most binary inspirals are expected to have circularized before they enter the LIGO/Virgo frequency band, a small fraction of those binaries could have non-negligible orbital eccentricity depending on their formation channel. Hence, it is important to accurately model eccentricity effects in waveform models used to detect those binaries, infer their properties, and shed light on their astrophysical environment. We develop a multipolar effective-one-body (EOB) eccentric waveform model for compact binaries whose components have spins aligned or anti-aligned with the orbital angular momentum. The waveform model contains eccentricity effects in the radiation-reaction force and gravitational modes through second post-Newtonian (PN) order, including tail effects, and spin-orbit and spin-spin couplings. We recast the PN-expanded, eccentric radiation-reaction force and modes in factorized form so that the newly derived terms can be directly included in the state-of-the-art, quasi-circular--orbit EOB model currently used in LIGO/Virgo analyses (i.e., the SEOBNRv4HM model).


I. INTRODUCTION
The observation of gravitational waves (GWs) by the LIGO-Virgo detectors [1,2] have corroborated the existence of binary black holes (BBHs) in our universe. But how and in which astrophysical environments these binaries form is not yet fully understood. However, the masses, spins (magnitude and orientation), and binary eccentricities inferred from GWs provide invaluable clues to determine BBH formation channels [3,4]. So far, the observed GWs are consistent with binary coalescences of negligible eccentricity, i.e., on quasi-circular orbits [5][6][7][8].
In general, binaries are expected to circularize [9,10] as they approach merger due to the emission of gravitational radiation. But depending on their astrophysical formation channel, a small fraction of binaries could have non-negligible orbital eccentricity, as they enter the frequency bands of current detectors. This can occur in dense stellar environments, such as globular clusters or galactic nuclei, where dynamic capture [11][12][13][14][15][16] or the Lidov-Kozai mechanism in hierarchical triples [17][18][19] can lead to eccentric binary inspirals at close separations.
In particular, Ref. [13] (and Ref. [14]) showed that ∼ 5% (or ∼ 10%) of all mergers in globular clusters enter the LIGO band with eccentricity e > 0.1. Binaries formed via dynamic capture in galactic nuclei are expected to have high eccentricities [16], with 92% having e > 0.1 and 50% − 85% having e > 0.8 at 10 Hz. For a BBH around a supermassive BH, the Lidov-Kozai mechanism can secularly drive the BBH to eccentricities near unity for some orientations [19]. Hence, inferring those eccentricities from GWs is important for understanding the origin and environment of BBHs. Interestingly, Ref. [8] pointed out that GW190521 [20] could be consistent with either an eccentric nonprecessing or a quasicircular precessing binary, which illustrates both the difficulties and prospects of further observations in the upcoming and future LIGO, Virgo and KAGRA runs [21].
Recent approaches to extend the EOB formalism to eccentric orbits include Ref. [78], which derived the RR force with eccentricity up to 2PN order, but without tail effects and for nonspinning BHs. More recently, Ref. [79] incorporated eccentricity effects in the RR force and in the (2, 2) waveform mode through 1.5PN order, including tail effects, using the Keplerian parametrization and phase variables that evolve only due to RR. References [80,81] extended the quasi-circular SEOBNRv1 [72] model to eccentric orbits, while Ref. [82] added eccentric corrections in the SEOBNRv4 [74,83] waveform model, notably in the (2, 2), (2,1), (3,3), (4,4) modes through 2PN order, including spin-orbit (SO) and spin-spin (SS) couplings 1 , but not tail effects. They employed these eccentric modes to construct a RR force for eccentric orbits, which however does not include the Schott terms. As argued in Ref. [78] and Sec. II below, these Schott terms are necessary for generic orbits to satisfy the flux-balance equations. Furthermore, Refs. [84,85] incorporated noncircular effects in the TEOBResumS SM [76,86] model at leading PN order in the azimuthal component of the RR force, and used a quasi-circular 2PN-expanded radial RR force without spin or tail effects. They included eccentric corrections at leading PN order to all modes m = 0 up to ℓ = |m| = 5.
In this paper, we develop a multipolar EOB waveform model for eccentric binaries with the compactobjects' spins aligned or antialigned (henceforth, for short aligned) with the orbital angular momentum. We derive the eccentric PN expressions for the RR force (including the Schott terms) and the gravitational modes up to ℓ = |m| = 6, including the m = 0 mode, through 2PN order, including tail effects, and SO and SS couplings. We recast our results for the RR force and modes in a form that can be directly incorporated in the state-ofthe-art, quasi-circular-orbit EOB model currently used in LIGO/Virgo analyses (SEOBNRv4HM [74,83]).
The paper is structured as follows. In Sec. II, we derive the RR force from the energy and angular momentum fluxes using the balance relations. We use the gauge freedom in the RR force to impose that it reduces to the relation used in SEOBNRv4HM in the quasi-circular-orbit limit. In Sec. III, we obtain initial conditions for eccentric orbits. In Sec. IV, we calculate all the gravitational waveform modes that contribute up to 2PN order relative to the leading order (LO) of the (2, 2) mode, i.e., up to the ℓ = |m| = 6 mode. These higher-order modes are even more important for eccentric orbits than for quasi-circular ones [87]. We conclude in Sec. V with a discussion of results and potential future work. Finally, Appendix A provides the coordinate transformation from harmonic to EOB coordinates, Appendix B includes a derivation of the LO spin-squared contribution to the angular momentum flux, Appendix C lists the spin contri- 1 Our results for those modes are mostly in agreement with Ref. [82] except for the SO part, where we disagree with their findings (their expressions contain two extra SO terms).
butions to the waveform modes in harmonic coordinates, Appendix D provides some relations for dynamic quantities in the Keplerian parametrization, and Appendix E includes the transformation to tortoise coordinates. We provide our results for the RR force and waveform modes as Mathematica files in the Supplemental Material [88].

Notation
We use the metric signature (−, +, +, +), and use units in which c = G = 1, but write c explicitly in PN expansions.
We consider an aligned-spin binary with masses m 1 and m 2 , with m 1 ≥ m 2 , and we define the following constants: In the binary's center of mass, we introduce the canonical phase-space variables (R, φ, P R , P φ ), where R is the separation, φ the azimuthal angle, P R the radial momentum, and P φ the angular momentum. The total relative momentum P is given by P 2 = P 2 R + P 2 φ /R 2 . We use the rescaled dimensionless variables where the dimensionless quantities are denoted with either a hat or a lowercase letter. The energy and angular momentum fluxes far away from the binary are denoted by Φ E and Φ J respectively, and scale as follows: where quantities with a tilde are the physical dimensionful fluxes. The components of the RR force are denoted by F r and F φ , and are scaled similarly to Φ E and Φ J , respectively.

II. RADIATION REACTION FORCE
The RR force accounts for the energy and angular momentum losses by the system, and is added to the righthand side of the Hamilton equations of motion (EOMs) such thatṙ where the leading order of F r,φ is of order 1/c 5 (2.5 PN). From the EOMs, with ∂H/∂φ = 0, the time derivatives of energy and angular momentum are given bẏ The energy and angular momentum lost by the system are not equal to the energy and angular momentum fluxes, Φ E and Φ J , because of additional contributions to E and J due to interactions with the radiation field. The balance equations are modified by Schott terms, as in electrodynamics, that appear as total time derivatives in the balance equations [78] E system +Ė Schott + Φ E = 0, Substituting the expressions for the energy and angular momentum losses, we obtaiṅ The energy and angular momentum fluxes are gaugeindependent, but the RR force and Schott terms are gauge-dependent. This coordinate gauge freedom in the RR force was discussed by Iyer and Will in Refs. [89,90], and by Gopakumar et. al. in Ref. [91]. Bini and Damour showed in Ref. [78] how the gauge freedom in F is related to the freedom in defining the Schott terms. Note that while we only consider aligned spins in this paper, an extension to precessing spins is straightforward; the RR force F is added to the EOM for the total momentum p and a RR contribution is added to the spin evolution equations, such that The balance equations are then given bẏ withĖ system =ṙ · F , See, e.g., Refs. [92,93] for more details.
A. Summary of the approach used in this paper for the RR force The aim of this paper is to extend the quasi-circular RR force and gravitational modes employed in the SEOBNRv4HM waveform model to eccentric orbits. The Hamilton equations that describe the dynamics of the SEOBNRv4HM model use the following relations between the RR force and the energy flux for quasi-circular orbits, which are based on results from Refs. [46,94], with Ω being the (angular) orbital frequency. However, these two relations are only valid for quasi-circular orbits and are not consistent for generic orbits, since they use the circular-orbit relation Φ qc E = ΩΦ qc J and do not include the Schott terms.
Hence, the approach we use to obtain the RR force is to write a generic ansatz with unknown coefficients for the Schott terms, and calculate the RR force from the fluxes using the balance equations Then, we specify the free unknown coefficients in the Schott terms such that the force reduces to the conditions in Eq. (11) in the limit of quasi-circular orbits, i.e., since both p r andṗ r are zero for circular orbits. Finally, we factorize the RR force into the quasi-circular part used in SEOBNRv4HM times eccentric corrections where the quasi-circular parts are given by Eq. (11), and the eccentric corrections scale as F nc i ∼ 1 +ṗ r + p 2 r + . . . . In the following subsections, we provide the details of these steps.

B. EOB Hamiltonian and angular momentum
The EOB Hamiltonian is calculated from an effective Hamiltonian H eff via the energy map with H eff given in Refs. [50,51,56]. When calculating the RR force to 2PN, we only need to work with the PN expansion of the EOB Hamiltonian. The nonspinning part to 2PN order is given bŷ the LO (1.5PN) spin-orbit part and the LO (2PN) spin-spin part where C iES 2 are the spin quadrupole constants, which equal one for BHs. The orbital frequency expanded to 2PN is given by From the EOMṗ r = −∂Ĥ/∂r, we can obtain an expression for p φ (r,ṗ r , p r ) which we use to express the noncircular part of the RR force and modes in terms of p r andṗ r . It will also be useful below, when taking the circular-orbit limit, to have an expression for p φ as a function of r for circular orbits. Settingṗ r = 0 = p r in the previous equation yields

C. Energy and angular momentum fluxes
The energy and angular momentum fluxes for nonspinning binaries were derived to 3PN order in harmonic and Arnowitt-Deser-Misner (ADM) coordinates in Refs. [95][96][97]. The 2PN instantaneous part of the fluxes for nonspinning bodies is given in EOB coordinates in Appendix A of Ref. [78]. The leading order reads The hereditary contributions to the fluxes can be expressed as an infinite sum over Bessel functions [79,96] that can be evaluated numerically, or resummed analytically [31,98]. Here, we follow the method from Ref. [79] to obtain the LO tail part (1.5 PN) of the orbit-averaged fluxes in an eccentricity expansion and we extend their derivation to O(e 6 ), which yields 2 where x ≡ Ω 2/3 . The eccentricity e in these equations is defined using the Keplerian parametrization, which is given by where u p is the inverse semilatus rectum and χ is the relativistic anomaly.
Since we are not using the adiabatic approximation and are not working with orbit-averaged fluxes, we can obtain an approximate expression for the tail contribution to the fluxes by writing an ansatz in terms of (r, p r , p φ ) in a p r expansion of the form calculate the average of that ansatz in terms of (e, x) (see Appendix D), and then match it to the average flux in Eq. (24) to determine the unknowns c n . This yields The LO (1.5PN) SO fluxes for generic orbits and generic spins were derived in Refs. [92,100]. (The nextto-leading-order (NLO) SO energy flux was derived in Ref. [101].) It should be noted that Ref. [100] used the Tulczyjew-Dixon (covariant) spin supplementary condition (SSC) [102][103][104][105], while Ref. [92] used the Newton-Wigner (NW), or canonical, SSC [106,107]. In this paper, we use the NW SSC since we are working in a canonical Hamiltonian formulation of the spinning two-body dynamics [108,109]. Changing the velocities in Eq. (17) where I ij is the mass quadrupole moment, and b = 2r 0 e −11/12 with r 0 a gauge parameter.
of Ref. [92] to momenta, which involves spin-orbit terms, the aligned-spin fluxes reduce to For the LO (2PN) SS contributions, the LO spin 1 -spin 2 energy and angular momentum fluxes in harmonic coordinates were derived in Refs. [93,100], while the spinsquared energy flux was derived in Refs. [110,111], and we obtain in Appendix B the spin-squared angular momentum flux. Transforming from harmonic to EOB coordinates, using the transformations in Appendix A, we get the following SS contributions to the fluxes for alignedspins: The total 2PN energy and angular momentum fluxes are the sum of all the above contributions, i.e.,

D. Ansatz for the Schott terms
As an ansatz for the Schott terms E Schott and J Schott , we consider Note that this ansatz for E Schott is more general than the one used in Eq. (4.4) of Ref. [78], since we found that such an ansatz is needed for the RR force to satisfy the conditions in Eq. (11).
For the LO tail, we use the ansatz while for the LO SO part, and for the SS part, The total energy and angular momentum Schott terms are the sum of the above contributions, i.e.
Note that when taking the time derivative of these Schott terms using the EOMs, the LO nonspinning part contributes to the LO SO and SS parts of the RR force.

E. Solving for the eccentric-orbits RR force
Using the fluxes and the Schott terms, the RR force can be calculated from the balance equations (12), which fix some of the unknowns in the ansatz for the Schott terms. The remaining unknowns can be determined by requiring that the RR force satisfies the conditions (13) in the circular-orbits limit.

Leading order
At leading order, calculating the RR force with the ansatz in Eqs. (31) and expanding in p r gives (36) Requiring that the 1/p r term is zero, leads to the solution . (38) Requiring that the first term in that series expansion is zero gives the solution With that solution, we obtain the LO RR force This force satisfies the conditions in Eq. (13) for circular orbits since

1PN
Following the same steps as above, we obtain the following solution for the unknowns at 1PN: with 3 arbitrary coefficients out of 9 coefficients at that order. To simplify the resulting expressions for the RR force, we choose to set all arbitrary coefficients to zero, i.e. α 3 = β 4 = β 5 = 0, which yields the following 1PN contribution to the RR force:

LO tail
Solving for the unknowns at the LO tail contribution leads to the solution with either λ 1 or λ 4 arbitrary. Choosing λ 1 = 0, the tail contribution to the RR forces becomes

LO spin-spin
At LO SS, we obtain the unique solution With that solution, we get F. Factorizing the RR force into circular and noncircular parts The total RR force is the sum of the contributions calculated in the previous section, i.e., We have checked that our gauge-dependent RR force agrees with that in Refs. [78,92,93] by using the balance equations. Denoting the RR force from those references by (F r ,F φ ) with corresponding Schott terms (Ē Schott ,J Schott ), Eqs. (7) lead tȱ Then, by writing an ansatz for (Ē Schott ,J Schott ) with unknown coefficients, we checked that a solution exists, implying that (F r , F φ ) and (F r ,F φ ) are related via a coordinate transformation.
To implement our results in the SEOBNRv4HM model, we factorize the RR force into a quasi-circular part times eccentric corrections as in Eqs. (14) and (11), which read and for the quasi-circular part we use the unexpanded force used in SEOBNRv4HM, in which the energy flux has the following PN expansion in terms of the orbital velocity v Ω ≡ Ω 1/3 : This leads to the eccentric part The full 2PN expressions are provided in the Supplemental Material [88].
In these eccentric corrections to the RR force, we useḋ p r instead of p φ because it improves the agreement of our model with SEOBNRv4HM in the quasi-circular orbit limit, in which p r = 0 =ṗ r leading to F ecc r,φ = 1. However, havingṗ r on the right-hand side of the EOM for p r would complicate solving the system of differential equations (4). Therefore, when evolving the EOMs, we simply replaceṗ r in the RR force with the derivative of the Hamiltonian with respect to r calculated numerically, i.e. p r → −∂Ĥ EOB /∂r. SEOBNR waveform models use p r * (the conjugate momentum to the tortoise radial coordinate r * ) instead of p r since it improves stability of the EOMs near the EOB event horizon [69,112]. The two momenta are related by Eq. (E3). In Appendix E, we also obtain Eq. (E8) for the transformation betweenṗ r andṗ r * .

III. INITIAL CONDITIONS
Having determined the RR force that enters the EOMs, we need to specify the initial conditions to be used in evolving the system of equations. In this section, we first review how the initial conditions are implemented in SEOBNRv4HM for quasi-circular orbits [46,94], and then discuss a simple extension for eccentric orbits.

A. Initial conditions for quasi-circular orbits
Let us recapitulate the initial conditions for quasicircular/spherical orbits in the SEOBNRv4HM model as derived in Refs. [46,94]. We start by specifying an initial orbital frequency Ω 0 , with initial orbital phase φ 0 = 0, and solve for the initial values of r and p φ , while neglecting RR, p r ≈ 0. The initial condition for p r is then obtained by solving for p r , after calculating [ṙ] 0 using the result from adiabatic evolution [46] [ṙ] 0 = dp φ /dt dp whereĖ is the circular-orbits energy flux, and the derivative dE/dr = dH/dr can be determined using the following equations for circular orbits: p r = 0, dp r = 0, d ∂H ∂r = 0.
This leads to d dr which can be solved for dp φ /dr to obtain dp φ dr Plugging that solution into dH/dr = (∂H/∂p φ )dp φ /dr yields the result in Eq. (4.14) of Ref. [94], which reads and hence The complete procedure to obtain the initial conditions for the orbital phase-space is now as follows. Given Ω 0 , masses, and spins, we numerically solve the relations in Eq. (60) for the initial values r 0 and p φ 0 , choosing φ 0 = 0 and assuming p r ≈ 0. Using these values, we numerically solve Eq. (68) for the initial value p r0 .

B. Initial conditions for eccentric orbits
Since eccentricity is a gauge-dependent concept, we do not need to calculate accurate initial conditions for eccentric orbits in a specific gauge. Instead, we can choose a measure for eccentricity that can be adjusted to be as convenient as possible for numerical implementation. The only strict requirement is that for zero eccentricity e = 0 one recovers the quasi-circular case. Hence, we can start with very accurate initial conditions for quasicircular orbits and perturb them for eccentric orbits.
We choose to specify an initial orbital frequency Ω 0 and an initial eccentricity e 0 using the Keplerian parametrization 1/r = u p (1 + e cos χ). We also assume that the orbit starts with φ 0 = 0 at periastron (χ = 0), where p r = 0 in absence of RR, which simplifies calculating the initial conditions for r and p φ . An advantage of starting at periastron instead of apastron is that the specified initial frequency is then the maximum orbital frequency (over the first orbit), and can be used to estimate the frequency at which the binary enters a GW detector's frequency band.
To obtain r 0 and p φ 0 , we solve Eq. (60) with a nonzerȯ p r , i.e., with p r ≈ 0, and [ṗ r ] 0 given as a 2PN expansion in terms of p φ and e. For quasi-circular orbits, these equations reduce exactly to Eqs. (60) sinceṗ r ∝ e.
To obtain the PN expansion forṗ r at periastron, we first invert the Hamiltonian at the turning points r ± = 1/(u p (1 ± e)) with p r = 0 and solve for the energy and angular momentum as functions of e and u p , which are given by Eqs. (D2). Then, we invert p φ (e, u p ) to obtain Eq. (D3) for u p (p φ , e) and insert it into the PN expansion forṗ r = −∂H/∂r at periastron (r = 1/[u p (1 + e)]). This yields The initial condition can now be obtained in analogy to the quasi-circular case: given Ω 0 , e 0 , masses, and spins, we obtain r 0 and p φ 0 from Eqs. (69) and (70) (assuming p r ≈ 0), p r0 then follows from Eq. (68) as before, and φ 0 = 0 by convention. We can keep using the circularorbits energy flux in Eq. (68), instead of replacingĖ with ΩF φ for eccentric orbits, because the difference on the orbital dynamics is negligible since it involves p r (which at periastron is numerically much smaller than r or p φ ).
To assess the accuracy of the initial conditions for eccentric orbits, we compare the specified value for the eccentricity in the Keplerian parametrization with the value calculated from the orbital frequency (or the frequency of the (2, 2) mode) at periastron Ω p and apastron Ω a , which is given by [113] e and to calculate it, we follow the steps explained after Eq. (2.8) of Ref. [41]. Since we evaluate e Ω by evolving the binary over one orbit (including RR), it only holds approximately that Ω p ≈ Ω a in the quasi-circular case. That is, e Ω does not vanish exactly for quasi-circular orbits, in contrast to our e 0 . Table I shows the value of the eccentricity e Ω calculated from the orbital frequency compared to the specified eccentricity e 0 , and we see good agreement between the two measures of eccentricity.

IV. GRAVITATIONAL WAVEFORM MODES
In this section, we obtain the 2PN waveform modes to 2PN order including LO tail effects, and SO and SS couplings for aligned spins. The instantaneous nonspinning part of the modes was derived for eccentric orbits in Refs. [114,115] in harmonic coordinates and we convert their results to EOB coordinates. For the LO tail part, we extend the results of Ref. [79] to O(e 6 ) and to higher modes using the Keplerian parametrization, and then convert those results to an expansion in p r andṗ r . The spin contributions to the modes were derived to 2PN order for circular orbits in Ref. [116], and here we derive them for eccentric orbits.
The GW spherical harmonic modes h ℓm are the expansion of the complex polarization waveform The modes h ℓm can be calculated directly from the radiative multipole moments via [117][118][119] where D L is the luminosity distance of the source, and the radiative multipole moments are related to the symmetric trace-free (STF) moments U L and V L by whereȲ ℓm L is the complex conjugate of the STF tensors relating the unit vectors N L (which point from the source to the detector) to the spherical harmonics basis Y ℓm (Θ, Φ) such that and we can express the unit vector N in terms of the angles Θ and Φ as N = sin Θ cos Φê x + sin Θ sin Φê y + cos Θê z .
For planar binaries, nonspinning or with aligned spins, it was shown in Ref. [119] that the modes can be determined using the mass-type multipole moments for even ℓ + m, or the current-type multipole moments for odd ℓ + m, i.e., We define H ℓm such that which makes the LO part of H 22 = x for circular orbits. Note that different conventions for the phase origin contribute a factor of (−i) m to the modes [120].
In this paper, we compute the modes to 2PN order beyond the leading order of the (2, 2) mode, which means we consider modes up to the ℓ = 6, m = even modes. To 2PN order, the instantaneous contributions to the radiative multipole moments coincide with the source multipole moments. Including the hereditary terms that contribute to 2PN, the radiative multipole moments are given by [115,117,118] where the constants b i are gauge parameters that will be eliminated via a phase shift as was done in Ref. [117]. The source multipole moments for nonspinning binaries are given in, e.g., Refs. [115,117], while the spin contributions to the source moments are given in Refs. [116,121].

A. Instantaneous nonspinning contributions
The instantaneous contributions to the modes for nonspinning binaries in eccentric orbits were derived in Ref. [114] to 2PN, and in Ref. [115] to 3PN. The results of Ref. [115] are in harmonic coordinates and in terms of the variables (r, φ,ṙ,φ). Hence, we can simply transform their results from harmonic to EOB coordinates using the transformations in Appendix A. For the (2, 2) mode we obtain The expressions for the other modes that contribute to 2PN, i.e., up to ℓ = |m| = 6, are provided as a Mathematica file in the Supplemental Material [88]. Note that the (ℓ, 0) modes are zero for circular orbits but not for eccentric orbits. For example, the LO part of the (2, 0) mode is given bŷ which is zero for circular orbits since p 2 = 1/r + . . . .

B. Hereditary contributions
The hereditary contributions to the modes can be calculated analytically in an eccentricity expansion, as was done in Ref. [79] for the (2, 2) mode to O(e 2 ), and in Ref. [122] for all modes to 3PN order and to O(e 6 ). The results of Ref. [122] use the quasi-Keplerian parametrization, while here we use the Keplerian parametrization following the method developed in Ref. [79], which is based on results from Refs. [96,117,123], to derive the leading order tail effects that contribute to the modes up to 2PN order and to O(e 6 ). (See Ref. [79] for a discussion of the advantages of the Keplerian parametrization over the quasi-Keplerian parametrization.) We finally convert the eccentricity-expanded tail contributions to an expansion in p r andṗ r .

1
, and the unit vectors n L are related to spherical harmonics via Eq. (76), leading to with the coefficients (for equatorial orbits) Decomposing the phase φ into an oscillatory part and a linearly growing part φ = φ 0 + ω φ t + ∆φ r allows expressing the oscillatory part ∆φ r as a Fourier series expansion. Hence, with the functions J ℓm defined by where ψ r and ψ φ are the radial and azimuthal angle variables associated with the frequencies ω r = dψ r /dt and ω φ = dψ φ /dt. The coefficients J ℓmk are given by where, in the last line, we assume φ 0 = 0. The function P denotes the conservative part ofχ and is related to the radial angle ψ r via dψ r /dχ = ω r /P, with P = (1 + e cos χ) 2 u 3/2 p at LO (see Ref. [79] for more details).
Thus, the Newtonian mass multipole moments can be expressed as where the azimuthal angle is related to the radial angle by ψ φ = φ − ∆φ r with ∆φ r = χ − ψ r at LO. This allows us to write the LO tail contribution to the mass-type radiative moments as where Ω mk ≡ mω φ + kω r and The exponential e −ikψr can be expressed in terms of χ and e by integrating dψ r = ω r dχ/P, with ω r = (u p − e 2 u p ) 3/2 , leading to Eq. (3.41) of Ref. [79], which reads e −ikψr = e ik e √ 1−e 2 sin χ 1+e cos χ Hence, the modes with ℓ = 2 and m even are given by while the modes with ℓ = 3 and m odd are given by To obtain analytical expressions for the modes, we expand the above equations in eccentricity, where the infinite sum over k can be stopped at the order of the expansion in e. The result of that expansion is complicated, but we can perform a phase redefinition in the leading order instantaneous part 3 of the form φ → φ + x 3/2 δ φ , and absorb in δ φ all terms that are not proportional to π 3/2 . This modifies the phase at 4PN relative order, which we can ignore when working to 2PN order. (See Ref. [117] for more details.) The result for the (2,2) and the (3, 1) modê We checked that our results agree with those of Ref. [122] after converting between the quasi-Keplerian and Keplerian parametrization, and performing a phase shift.
To express the modes in terms of (r, p r , p φ ) instead of (x, e, χ), we use the following leading order relations: As explained above, it is advantageous to replace p 2 φ withṗ r usingṗ r = (p 2 φ − r)/r 3 and expand in both p r andṗ r (since p r andṗ r are both of order e) to obtain  The Newtonian order current quadrupole moment is given by [96] whereê i z is the unit vector in the z-direction. The term e i z n j can be expressed in terms of Y ij 21 , as was done in Ref. [124], by defining the complex vector which leads to Since, for equatorial orbits, n i = cos φê i x + sin φê i y and λ i = − sin φê i x + cos φê i y , we obtain Hence, e i z n j = Re e −iφ e i z ζ j = −2 Since V ij is contracted withȲ ℓm ij in Eq. (74), and Y ℓm ijȲ ij ℓm = 0, only the term with Y ij 21 in the above equation contributes to the modes. Thus, we only need to consider the following part of the current quadrupole Then, we follow the same steps as in the previous subsection. Decomposing the phase into φ = ω φ t + ∆φ leads with and Thus, the current quadrupole source moment can be expressed as (112) and the current quadrupole radiative moment leading to the (2, 1) mode which is in agreement with the results of Ref. [122]. In terms of (r, p r , p φ ,ṗ r ), we obtain The spin contributions to the modes were derived for circular orbits in Refs. [116,120,125]. To derive the spin part of the modes to 2PN for eccentric orbits, we use the source moments from Refs. [116,121], which are in harmonic coordinates and in terms of the covariant SSC. Dif-ferentiating the source moments to obtain the radiative moments (82), and plugging them into Eq. (73), we obtain the modes listed in Appendix C. Transforming from harmonic to EOB coordinates, and from the covariant to the NW SSC using the transformations in Appendix A, we obtain the following spin contributions to the modes: , The circular-orbit limit of these modes, when expressed in terms of the orbital frequency, agrees with the results of Refs. [116,125]. The spin contributions to the (2, 2), (2, 1), and (3, 3) modes for eccentric orbits were calculated in Ref. [82]; however, we find a small disagreement with their results for the SO part. 4

D. Factorized modes
The quasi-circular waveform modes used in SEOBNRv4HM are factorized as follows [70,83,126,127]: where h N,qc ℓm is the Newtonian part of the mode,Ŝ qc eff is an effective source term given bŷ T qc ℓm resums the infinite number of "leading logarithms" entering the tail effects, δ ℓm contains the part of the tail not included in T qc ℓm , and f ℓm contains PN corrections such that the expansion of h F,qc ℓm agrees with the known PN expansion of the modes. See Refs. [70,83] for more details and for expressions of these terms.
which is likely due to the coordinate/SSC transformations detailed in Appendix A.
We include the eccentric corrections in the factorized modes as follows: (120) where the effective source term is given bŷ T ecc ℓm contains the eccentric corrections to the hereditary contributions, δ ℓm is the same as in the quasi-circular case, and f ecc ℓm contains the eccentric corrections to the instantaneous contributions (both spinning and nonspinning, including the Newtonian part). For example, for the leading order of the (2, 2) mode, we obtain f ecc 22 = 1 2(r 2ṗ r + 1) 1/3 2 + r 2ṗ r − rp 2 r − 2 r 2ṗ r + 1 For the tail part, we simplified the results of Sec. IV B and eliminated the gauge parameter by using a phase shift, which led to the circular part of the tail contribution to the (2, 2) mode simply being 2πv 5 Ω ; however, this phase redefinition is not done in SEOBNRv4HM, and the corresponding expression reads v 5 Ω (2π + 12i log (2ǫv Ω ) − 17i/3 + 12iγ E /3). Therefore, when including the eccentric corrections in T ecc ℓm , we assume that the phase redefinition was done only for the eccentric part and keep using the same circular part as in SEOBNRv4HM. In addition, since we expanded the tail part in eccentricity to O(e 6 ), when factorizing the modes as in Eq. (120) and writing the quasi-circular part in terms of frequency, we reexpand T ecc lm in eccentricity (or p r anḋ p r ). For example, for the (2, 2) mode, we obtain The full expressions for T ecc lm and f ecc lm are provided in the Supplemental Material [88].

V. CONCLUSIONS
Extending the waveform models used today in GW astronomy from quasi-circular to eccentric orbits is important for future observations with LIGO, Virgo and KAGRA detectors [21], and with new facilities on the ground (Cosmic Explorer and Einstein Telescope), and in space (LISA). In fact, sources with non-negligible eccentricity might come into reach of observations soon and should routinely be included in searches and parameter inference. While this presents a challenge for waveform modeling and data analysis, it also offers the unique opportunity to unveil the formation channels of compact binaries and probe their environment (through eccentricity measurements). In this paper, we constructed an EOB waveform model for eccentric binaries. For this purpose, we obtained analytical results for the RR force and waveform modes to 2PN order, including the leading-order tail effects, and SO and SS couplings for aligned spins.
In particular, we first derived the RR force for eccentric orbits in PN expanded form, and then we recast it in a form that it can be directly incorporated in the quasicircular RR force employed in the SEOBNRv4HM [74,83] model, currently used in LIGO/Virgo analyses [1]. We then obtained initial conditions for the binary evolution which generalize those from Ref. [94] to eccentric orbits, and which allow starting the binary's evolution from a specified initial frequency at periastron and an initial eccentricity (in the Keplerian parametrization). We also calculated all the waveform modes that contribute up to 2PN order relative to the leading order of the (2, 2) mode. It should be noted that the (ℓ, 0) modes are proportional to the eccentricity and are hence important for eccentric orbits, especially the (2, 0) mode since it starts at the same PN order as the (2, 2) mode. Also the gravitational modes were rewritten in a factorized form to be straightforwardly implemented in the SEOBNRv4HM model.
Our results for the RR force and modes are valid for moderate to high eccentricities during the inspiral phase, since we do not use an eccentricity expansion except for the tail part, which is known analytically as an infinite series expansion. We provided expressions for the tail part in an expansion to O(e 6 ), but we checked that expanding to O(e 10 ) produces negligible difference on the waveform even for high eccentricities ( 0.9). If results for e close to 1 are needed, one could calculate the series expansion for the tail part numerically, or use analytical resummation methods as was done in Refs. [31,98].
We are currently incorporating the eccentric RR force and gravitational modes of this paper in the inspiralmerger-ringdown quasi-circular-orbit SEOBNRv4HM waveform model (SEOBNRv4EHM [128]) and validating it against NR simulations with eccentricity. We leave to future work the extension of the model to higher PN orders and the inclusion of spin precession.

ADM to EOB transformation
To find the canonical transformation from the ADM Hamiltonian with LO SO and SS using the NW SSC (see e.g. Refs. [129][130][131]) and the 2PN expansion of the EOB Hamiltonian of Ref. [50], we write an ansatz with unknown coefficients for the generating function G(x, p), perform the following transformation on the ADM Hamiltonian [45]: and match it to the EOB Hamiltonian to solve for the unknowns.
The result for the generating function is given by which has no LO SO terms since the ADM and EOB Hamiltonian are the same at that order. This generating function yields

Harmonic to EOB transformation
The transformation from harmonic to ADM coordinates is given by Eq. (E1) of Ref. [78], which is independent of spin since the ADM and harmonic coordinates agree at LO SO and SS. Using that equation together with Eq. (A3), we obtain the following transformation from harmonic to EOB coordinates: form the resulting modes to the NW SSC, we use the center-of-mass shift [100] x and the spin transformation [132] S cov where the spin transformation is only required for the NLO SO part of the 2PN (2, 1) mode. For the scalars (r, φ,ṙ,φ, χ 1 , χ 2 ), we obtain the transformations r cov = r − νrφ 2c 3 (χ 1 + χ 2 ), φ cov = φ + νṙ 2c 3 r (χ 1 + χ 2 ), r cov =ṙ + νṙφ 2c 3 (χ 1 + χ 2 ), φ cov =φ − ν 2c 3 r 3 1 + rṙ 2 − r 3φ2 (χ 1 + χ 2 ), Appendix B: Angular momentum flux at leading-order spin-squared In this appendix, we derive the angular momentum flux at leading spin-squared (S 2 i ) order. Here, we use unscaled variables in harmonic coordinates, but we drop the subscript 'h' to simplify the notation. We denote the orbital angular momentum L = µr × v, the relative position r = x 1 − x 2 , and relative velocity v = dr/dt.
Since the spin evolution equations start at 1PN order, we can assumeṠ 1 = 0 =Ṡ 2 for the calculation of the LO fluxes.
The source multipole moments needed are the spin quadrupole I ij and the current quadrupole J ij , which are given by [100,110,121] where the indices in angle brackets denote a symmetric trace-free part.
To transform from the coordinates of the two bodies x i 1 and x i 2 to the center-of-mass relative coordinates x i = x i 1 − x i 2 , we use [133] x i 1 = where The energy and angular momentum fluxes in terms of the multipole moments, to the order needed for the LO fluxes, are then calculated from [100,134] This yields the LO SO and S 1 S 2 fluxes derived in Refs. [92,93,100], in addition to the S 2 i energy flux from Ref. [110]. For the S 2 i angular momentum flux, we obtain L µr S 2 1 − n · S 1 (v × S 1 ) + v · S 1 (n × S 1 ) + 2m 2 2 C 1ES 2 5c 4 M r 4 L µr S 2 1 −30ṙ 2 + 12v 2 + 24 M r + L µr (n · S 1 ) 2 210ṙ 2 − 60v 2 − 90 M r This is in agreement with the recent results of Ref. [135], although our expression appears simpler because of using the individual spins S i and masses m i , instead of different combinations of them.
In the Keplerian parametrization, r = 1 u p (1 + e cos χ) , where u p is the inverse semilatus rectum and χ is the relativistic anomaly. Inverting the Hamiltonian at the turning points r ± = 1/(u p (1 ± e)) and solving for the