Spin effects in gravitational waveforms and fluxes for binaries on eccentric orbits to the third post-Newtonian order

Compact binaries can have non-negligible orbital eccentricities in the frequency band of ground-based gravitational-wave detectors, depending on their astrophysical formation channels. To accurately determine the parameters of such systems, waveform models need to incorporate eccentricity effects. In this paper, we consider an eccentric binary of spinning nonprecessing compact objects, and derive the energy and angular momentum fluxes at infinity, as well as the gravitational waveform modes to the third post-Newtonian order. The novel results of this paper include the next-to-leading order instantaneous spin-orbit and spin-spin contributions to the waveform modes, in addition to the hereditary (tail and memory) contributions to the modes and fluxes for eccentric orbits. The instantaneous contributions are derived for generic motion, while the hereditary contributions are computed in a small-eccentricity expansion, but we consider a resummation that makes them valid for large eccentricities. We employ a quasi-Keplerian parametrization of the motion using harmonic coordinates and the covariant spin-supplementary condition, which complements some results in the literature in other coordinates. Our results can be useful in improving the accuracy of waveform models for spinning binaries on eccentric orbits.

Most waveform models rely on the post-Newtonian (PN) approximation (small-velocity and weak-field expansion, i.e., v 2 /c 2 ∼ GM/c 2 r ≪ 1) to describe the binary dynamics.The PN conservative dynamics is known for generic orbits to 5.5PN order [44][45][46][47][48], though the nonlocal-in-time contributions need to be computed separately for bound and unbound orbits in an eccentricity expansion [49,50].The radiative sector (energy and angular momentum fluxes, and the waveform phase and amplitude) is generally more difficult to derive than the conservative dynamics, and currently the PN information is known at a lower order for eccentric than circular orbits.In the nonspinning case, and for quasi-circular inspirals, the fluxes and GW phase have recently been derived to 4.5PN, with the dominant (2,2) waveform mode to 4PN [51,52], and the higher modes at 3.5PN [53][54][55][56][57].However, for eccentric orbits, the radiative dynamics is known to 3PN; in particular, the energy flux was derived in Refs.[58,59] and angular momentum flux in Ref. [60], while the waveform modes were completed to 3PN in Refs.[61][62][63].
The spin contributions, for quasi-circular inspirals and nonprecessing spins, are fully known in the waveform modes and energy flux to 3.5PN [64][65][66][67][68][69][70][71][72][73].However, for eccentric orbits, the spin part is known in the modes only to 2PN [74,75], which includes the leading order (LO) spin-orbit (SO) and spin-spin (SS) effects.The instantaneous contributions to the energy and angular momentum fluxes are known to 3PN from Refs.[72,76,77], which includes next-to-leading (NLO) SO and SS contributions, but the hereditary contributions (tail and memory effects) are known for quasi-circular orbits only [67,70,78,79].The highest known order in the instantaneous contributions to the energy flux is the next-to-NLO SS at 4PN, which was recently derived in Ref. [78] for generic orbits and precessing spins.
The instantaneous contributions to the fluxes and waveform modes can be derived in a closed form for generic orbits, but the hereditary contributions have only been computed in a small-eccentricity expansion.A useful technique in the calculation of the hereditary contributions is the use of the quasi-Keplerian (QK) parametrization [80][81][82], which was derived for nonspinning binaries to 3PN in harmonic and Arnowitt-Deser-Misner (ADM) coordinates in Ref. [83], and to 4PN in ADM coordinates in Ref. [84].The 3.5PN SO and SS contributions were derived in Refs.[85,86] in ADM coordinates and using the Newton-Wigner spin-supplementary condition (SSC) [87,88].
In this paper, we complete the fluxes and waveform modes to 3PN for eccentric orbits and nonprecessing spins.The novel results in this paper are the following: 1.The spin contributions in the QK parametrization to 3PN in harmonic coordinates and using the covariant Tulczyjew-Dixon SSC [89,90].
2. The spin contributions to the tail part of the energy and angular momentum fluxes at infinity, computed in a small-eccentricity expansion to O(e 8 ).We also compute the orbit-averaged fluxes using the QK parametrization, resum the tail contribution to improve its validity for high eccentricities, and obtain the time evolution of the secular orbital elements from the fluxes.
3. The spin contributions to the waveform modes at the 2.5PN and 3PN orders, which includes the instantaneous and tail contributions, in addition to direct-current (DC) and oscillatory memory contributions.We compute all hereditary contributions to the modes in a small-eccentricity expansion to O(e 6 ).Note that even for circular orbits, our result for the spin part of the memory beyond 2PN is new.
This paper is structured as follows: in Sec.II, we summarize the approach we use for the calculations, which is the PN multipolar post-Minkowskian (PN-MPM) formalism, and explain some techniques for evaluating the hereditary contributions based on the QK parametrization.Sections III and IV include our results for the fluxes and waveform modes, respectively.We conclude in Sec.V and include some expressions for the QK parametrization in Appendix A. All the results presented in this paper are written in harmonic coordinates using the covariant SSC, and are provided as Mathematica files in the Supplemental Material [91].
For a binary with masses m 1 and m 2 , we assume m 1 ≥ m 2 and define where M is the total mass, µ is the reduced mass, ν is the symmetric mass ratio, and δ is the anti-symmetric mass ratio.
For a spinning binary with (constant-magnitude) spins S 1 and S 2 , the dimensionless spin magnitudes are denoted χ 1 ≡ c|S 1 |/(Gm 2 1 ) and χ 2 ≡ c|S 2 |/(Gm 2 2 ).The spin-quadrupole constants are denoted κ 1 and κ 2 , which equal one for black holes.To simplify our expressions, and to make it easier to specialize them for black holes, we define the following combinations of the spin magnitudes and quadrupole constants: ) κ+ ≡ Several quantities are used in the QK parametrization, as explained in Subsection II C, and we summarize them in Table I below.For a geometric picture for some of these quantities, see Fig. 2 of Ref. [86].

II. FLUXES AND WAVEFORM IN THE PN-MPM FORMALISM
In this section, we briefly recall some aspects of the PN-MPM formalism.Then, we derive the QK parametrization, and finally describe how to evaluate the hereditary integrals that are required for the 3PN waveform and fluxes.
The PN-MPM formalism is an iterative algorithm that allows computing, to a given PN order, the radiated fluxes and the waveform through a multipolar decomposition of the system [92].The relevant quantities are parametrized by a set of symmetric trace-free multipole moments {U L , V L }, called "radiative" because they are defined at infinity on an asymptotically flat spacetime.The PN-MPM formalism allows linking the radiative multipole moments to the "canonical" multipole moments {M L , S L } parametrizing the source and computed in harmonic coordinates.
Both sets of multipole moments differ by nonlinear interactions of the GW field.Notably, hereditary effects, such as tail and memory integrals, appear at the considered PN order.Fortunately for us, the canonical moments were derived in Refs.[53,55,73] for arbitrary motion up to the 3.5PN order.Thus, we can obtain the fluxes and waveform modes from those moments, which is straightforward for the instantaneous contributions, but for the hereditary contributions, we first derive the QK parametrization, then express the canonical moments within this parametrization to evaluate the hereditary integrals.

A. Fluxes and waveform modes in terms of radiative multipole moments
Let us consider an asymptotically flat spacetime with a background Minkowski metric η µν .We define h µν ≡ √ −gg µν − η µν being the deviation to the gothic metric, where g is the determinant of g µν .Within the harmonic gauge, which imposes ∂ ν h µν = 0, the Einstein field equations read where τ µν is the usual Landau-Lifshitz pseudo stress-energy tensor.There exists a coordinate system (T, X), called radiative coordinates, in which the solution metric, in the transverse and traceless (TT) gauge, reads [93] h where ε ijk is the Levi-Civita symbol, P ijkl = P i(k P l)j − 1 2 P ij P kl is the usual TT projector, with P ij = δ ij − N i N j being the projector orthogonal to the unit vector N , which is the direction to the observer.The relation (2.2) defines the radiative multipole moments {U L , V L }.After introducing the orthonormal triad (P , Q, N ), we can define the usual GW polarizations ) Next, the amplitude of the GW can also be linked to the radiative moments.After a mode decomposition and projecting the multipoles from the STF to a spherical-harmonics basis, the complex amplitude h reads where with R being the distance between the source and the observer, and the STF tensors α ℓm L are defined by The explicit expression for α ℓm L is given by, e.g., Eq. (4.7) of Ref. [56].Note that in the case of planar orbits, the mass-type and current-type multipoles do not mix within a given (ℓ, m) mode.More specifically for even ℓ + m, only the mass-type moments U L contribute, while for odd ℓ + m, only the current-type moments V L contribute.
Let us now focus on the radiated energy and angular momentum.The radiated energy at future null infinity can be written in terms of the TT metric at the lowest order in R as [93] (2.7) Using Eq. (2.2), we can write the energy flux in terms of the radiative multipole moments, and a similar procedure can be done to express the angular momentum flux [93,94].In the end, both fluxes read L U (1) (2.8b) We recall that we are interested in the spin contributions to the fluxes and modes to the 3PN order.The spin contributions to the mass-type multipoles arise at 1.5PN, while they contribute to the current-type multipoles at 0.5PN.Thus, restricting ourselves to the 3PN order for the spin contributions, the energy and angular momentum fluxes are obtained only from the mass and current-type multipoles for ℓ = 2 and ℓ = 3.In the next subsection, we explain how the radiative multipole moments are derived.

B. Radiative multipole moments
The radiative moments are composed of two types of contributions: instantaneous ones, which depend on the local time and can be computed for arbitrary motion; and hereditary contributions, which depend on the entire history of the binary and require specifying the orbit.Furthermore, as we see later, an eccentricity expansion within the QK parametrization is required to evaluate the hereditary integrals analytically.In this subsection, we restrict ourselves to the moments that affect the spin terms at the 3PN order.An exhaustive list of nonlinearities coming from the difference between radiative and canonical moments can be found in, e.g., Ref. [55].
The general structure of the radiative moments in terms of the canonical moments at 3PN order read

Instantaneous contributions
The instantaneous part of the radiative multipole moments are generally given by (2.10a) where δU inst L and δV inst L are local-in-time, nonlinear interactions of the gravitational field.The full expressions of these interactions for the 3.5PN waveform are provided in Ref. [55].At the considered order of this paper, the spin contributions of these interactions are given by [73] ) where the subscript S refers to the spin part of the interaction.The rest of δU inst L and δV inst L do not contribute to the spin terms for the waveform at 3PN.Both the linear and nonlinear parts of the radiative moments can be computed for arbitrary motion, i.e., without the need for specifying the orbit.

Tail contributions
The tail integrals contributing to the 3PN spin terms in the GW fluxes and modes come from those of the mass quadrupole, current quadrupole, and current octupole.Their expressions read [95] where M is the ADM mass, T R is the retarded time, and b is an arbitrary gauge parameter that can be eliminated after a phase redefinition, as we discuss in Subsection IV D. Contrary to the instantaneous contributions, these integrals require specifying the binary's orbit to be evaluated analytically.We detail the procedure in Subsection II D.

Memory contributions
The memory contributions to the radiative moments can be derived within the PN-MPM formalism.They take the form of integrals over the past of products of multipoles without logarithms.Formally, these integrals are of order O(c −5 ); for example, the leading order memory term in the radiative mass-quadrupole moment reads j⟩a (τ ) dτ + O(7). (2.13) These memory contributions have been derived consistently for each multipole entering the 3.5PN waveform in Ref. [55].Furthermore, recent work [96] tackled the computation of the so-called tail-of-memory effects that formally appear at 4PN.
However, as shown in Subsection II D, the explicit evaluation of these integrals can lower the PN order of the memory contributions to the modes.Thus, in order to be consistent to a given PN order in the waveform amplitude, one would need to push the MPM procedure to higher orders, which is difficult to do using this method.Fortunately, one can bypass this difficulty by computing the memory contributions to the waveform modes using another method, allowing us to be consistent in the present case to 3PN.Following Refs.[95,[97][98][99], one can extract from the flux, Eq. (2.7), the memory contributions to the amplitude modes.They read With this relation at hand, one can express the flux in terms of the time derivatives of the modes instead of the radiative multipole moments.Then, after performing the angular integration, as explained in detail in Ref. [97], Eq. (2.14) becomes an integral over time of the following source where the value of the angular integral G ss ′ s ′′ ll ′ l ′′ mm ′ m ′′ can be found in Appendix A of Ref. [97].We then use the modes computed in Subsections IV A and IV B, take their time derivative and insert them in Eq. (2.15) which gives the source of the memory integrals.In Subsection II D, we explain how to evaluate the elementary integrals generated by those sources.

C. Quasi-Keplerian parametrization
We parametrize the motion of an eccentric-orbit binary using a Keplerian-like parametrization.The periastron advance implies that we cannot take a purely Keplerian parametrization as the motion is no longer elliptic.Instead, we choose the quasi-Keplerian parametrization introduced in Ref. [80], and extended to spin in Refs.[85,86,100].
The starting point of this computation is the conserved energy E and angular momentum J , which were derived in harmonic coordinates to 3PN for spinning binaries in Refs.[72,101].The spin definition in these papers, which is the one we adopt here, is chosen to use the covariant (Tulczyjew-Dixon) SSC.Therefore, we rederive the spin contributions to the QK parametrization using the covariant SSC and harmonic coordinates, as opposed to the Newton-Wigner SSC and ADM coordinates used in Refs.[85,86].
We restrict ourselves to the nonprecessing case, which implies that the motion is planar.When considering spin effects, it is better to work with the orbital angular momentum L = J − (S 1 + S 2 )/c, which is also conserved for nonprecessing spins.We use the orthonormal unit vectors ( l, n, λ), where n is the unit vector in the radial direction, l is the direction of the orbital angular momentum, and λ ≡ l × n.Then, we express the conserved energy and angular momentum, computed in [72,101], in terms of ṙ and φ and obtain, at leading order as an example, where E is (minus) the reduced binding energy and h is the reduced angular momentum, i.e., (2.17) We can invert the above relations for E and h to obtain an expression for ṙ and φ in terms of E, h and s ≡ 1/r by replacing iteratively the values of ṙ and φ appearing at higher PN orders.We find ṙ2 = P(s) , (2.18a) where P and Q are two polynomials of s whose coefficients are explicit expressions of E and h.Note that, at the 3PN order for nonspinning bodies in harmonic coordinates, ṙ2 is not a polynomial of s due to the appearance of ln(r) in the conserved energy and angular momentum.Thus, one has to change the coordinate system to remove these logarithms in order to apply the following procedure [83].However, to derive the spin contributions in the waveform to the 3PN order, we do not need to include the 3PN nonspinning contributions.Hence, we derive the QK parametrization including 2PN nonspinning and 3PN spinning contributions without changing coordinates.
To obtain the QK parametrization, we followed the well-explained method in Ref. [83] which we briefly summarize here.The first step to integrate Eqs.(2.18) is to remark that for bound orbits, s goes through a maximum and a minimum s + and s − corresponding to the minimal and maximal value of the separation r − and r + , respectively, for quasi-elliptic orbits [102].Using this property, one can write where R is also a polynomial.The above relation fixes s + and s − .Furthermore, we impose the radial motion to follow the ansatz where the semi-major axis a r and radial eccentricity e r read a r = 1 2 Thus factorizing P as in Eq. (2.19) directly fixes a r , e r and R(s).At Newtonian order, we recover the Keplerian parametrization of the orbits which means that R(s) is a constant.
We next focus on the integration of Eq. (2.18a).We define the total radial period P as the time integral over one orbit This integral is straightforward to compute; we can then express P in terms of a r and e r , and then the conserved quantities.Then, we compute the time dependency of the QK parametrization using Finally, we define the mean motion n ≡ 2π/P and mean anomaly l = n(t − t 0 ).The results of the two integrals above can be expressed in terms of the eccentric and true anomalies u and v, defined in Eq. (2.26) below.
We treat Eq. (2.18b) in the same way as Eq.(2.18a).We define the advance of periastron angle as the angular integral over one orbit which is computed the same way as P .Then, the phase variable is computed as a function of u and v through Finally, we get an expression for 2π Φ (ϕ − ϕ 0 ) which constitutes the last piece of the QK parametrization that we adopt.In harmonic gauge, our result reads where we have introduced the time eccentricity e t and the phase eccentricity e ϕ .This choice of parametrization is done so that at Newtonian order, we recover the Keplerian parametrization [80].The phase eccentricity e ϕ is fixed by imposing that ϕ does not contain any term proportional to sin(v).
The explicit expressions for the main quantities are displayed in Appendix A. In the nonspinning case, we recover the results in Ref. [83] to 2PN, at which the structure of the 3PN spin contribution is similar to the 2PN nonspinning contribution.The 3.5PN SO and SS contributions are known from Refs.[85,86] in ADM coordinates and the NW SSC.Our results differ from those references because of using harmonic coordinates and the covariant SSC, but we checked that the gauge-invariant quantities, such as n and Φ when expressed in terms of conserved quantities, are in agreement.
With this parametrization at hand, we can express (r, ṙ, ϕ, φ) in terms of (x, e t , u), where x is a gauge-independent PN parameter defined as where Ω = Kn is the orbital frequency.First, we express the dynamical quantities in the QK parametrization in terms of x and e t , then make use of the chain rule which gives, e.g., for the separation and we derive ṙ(x, e t , u) using Eq.(2.26d) by eliminating v for u.The same steps can be used to derive φ(x, e t , u).
Note that the above results for the QK parametrization correspond to the conservative part of the dynamics; the radiation-reaction part, coming from the emission of GWs, is derived in Subsection III E by taking into account the energy and angular momentum radiated by the system.For now, the expressions of the conservative part are exact in the sense that no eccentricity expansion has been performed.One can also choose to express the variables of the problem in terms of (x, e t , l), which requires performing a small eccentricity expansion as we see in the following subsection.

D. Evaluating hereditary contributions using the QK parametrization
As mentioned, there are two types of hereditary integrals to compute: the tail integrals containing logarithms, and the memory integrals.The latter enter the even ℓ + m modes and can be split into two classes of contributions: the DC memory, affecting only the (ℓ, 0) modes, and the oscillatory memory.For both the memory and tail contributions, we compute the integrands for generic orbits, i.e., in terms of (r, ṙ, ϕ, φ).Then, we express them in terms of (x, e t , l) in a small-eccentricity expansion, in order to evaluate the integrals analytically.

Tail integrals
To express (r, ṙ, ϕ, φ) as functions of (x, e t , l), we invert Eq. (2.26b) to obtain u as a function of l.This relation reads (see, e.g., Refs.[62,103]) with where J k are Bessel functions and α k reads with )/e ϕ .Note that the value of α k is specific to this computation and has a more general form for higher PN orders in the QK parametrizations.Once u(l) is known, one can insert it into the expressions for (r, ṙ, ϕ, φ) in terms of (x, e t , u), and perform an eccentricity expansion to get them as functions of (x, e t , l).
An important feature of this parametrization, and especially of ϕ(x, e t , l), is the fact that the phase variable can be split into two parts: a linear-in-l part given by λ ≡ Kl, and an oscillatory part W (l) (2.32) After expressing ϕ as a function of l using Eqs.(2.26c) and (2.29), and performing an eccentricity expansion, we identify W (l) = ϕ(l) − λ.This split is useful in evaluating the hereditary integrals.
With this parametrization at hand, we compute the integrands in Eqs.(2.12) for generic orbits, then express them in terms of (x, e t , l).After simplifying the expressions, having performed an eccentricity expansion, we evaluate the following elementary integrals: where ω ∈ R * , k ∈ N, γ E is the Euler constant and H k is the k th harmonic number.Note that the case ω = 0 does not appear in this computation.

Memory integrals
Regarding the memory contributions, we start from Eq. (2.15) and we need to compute ḣℓm .We computed them in two different ways: the first is to take the derivative of the radiative moments and project them on the spherical harmonic basis defined above; the second is to start from the modes themselves and differentiate them using [63] (2.34) Both approaches yield the same results.Once the source of the time integrals is known, we get two types of elementary integrals: oscillatory and non-oscillatory.
The oscillatory memory integrals are of the form which was derived in Appendix C of Ref. [63] in the case of nonspinning binaries.After generalizing this method to spins, we found that at the considered order for the spin terms, the form of the integral remains the same, namely This is due to the fact that we do not need to include post-adiabatic corrections for the computation of the spin terms to 3PN order in the waveform modes.However, if one were to compute the memory integrals for the 3.5PN waveform, one would need to include the radiation reaction and thus I r,s would need to be generalized.
A crucial point to note is that when r = −s in Eq. (2.36), the PN order of I r,s is lowered compared to the source of this integral, since K = 1 + O(2), which complicates the power counting when trying to be consistent to a given PN order.The oscillatory memory effects only enter the even ℓ + m modes, and we derived them up to the ℓ = 7 modes.The results are displayed in Subsection IV C and in the Supplemental Material.
The second type of memory integrals are the non-oscillating (or DC) integrals, which take the form (2.37) These effects only enter the (ℓ, 0) modes.For the purpose of this computation, we need to introduce an initial eccentricity e 0 , and substitute the time variable in terms of the eccentricity Note that in the practical computations, we used for ė the value of its orbit average, since we are interested in the secular contribution to these integrals.In Subsection III E, we derive the secular evolution of the orbital elements and notably ⟨ ė⟩.Then, by knowing the expression of ⟨ ė⟩ as a function of (x, e) and the explicit value of x(e) and we can compute J p,q .

III. 3PN ENERGY AND ANGULAR MOMENTUM FLUXES FOR NONPRECESSING SPINS
In this section, we present the 3PN energy and angular momentum fluxes, which are related to the radiative moments via Eqs.(2.8).We split the fluxes into instantaneous and tail contributions, such that where the instantaneous contributions are valid for generic motion, while the tail contributions are computed for bound orbits in a small-eccentricity expansion to O(e 8 ).

A. Instantaneous contributions to the fluxes
Using the instantaneous part of the radiative moments from Eqs. (2.10) and plugging them in Eqs.(2.8), we obtain the instantaneous contributions to the fluxes to 3PN order for the spin terms.The 3PN nonspinning part of the energy flux is given by Eqs.(5.2) of Ref. [59], while the angular momentum flux is given by Eqs.(3.4) of Ref. [60].We only computed the nonspinning part to 2PN, and found agreement with the results of those references.
For the energy flux, we obtain the following SO and SS parts: This result agrees with Refs.[72,77,78].However, the flux is defined in Ref. [78] in terms of the source moments instead of the radiative moments, which means the 3PN SO contribution was not taken into account, while we find agreement on the overlapping terms.Nonetheless, these terms vanish when taking an orbit average.
For the angular momentum flux, we obtain the following SO and SS parts:  where this result agrees with Ref. [77], except for the 3PN SO contribution, which was not accounted for in that paper.

B. Tail contributions to the fluxes
To compute the tail contribution to the fluxes, we use the radiative moments from Eqs. (2.12) computed in a small-eccentricity expansion, as explained in Subsection II D, then plug the moments in Eqs.(2.8).
To the 3PN order, the only spin terms in the tail part of the fluxes are SO contributions at 3PN order.The nonspinning tail contributions, starting at 1.5PN, were derived to 3PN in Refs.[58,60,63].We write here the SO part of the tail to O(e 2 ) and provide the full expressions to O(e 8 ) in the Supplemental Material.Henceforth, we use e instead of e t to ease the notation.
For the energy flux, we obtain  We checked that our results agree in the circular-orbit limit with Ref. [78].

C. Orbit-averaged fluxes
To compute the orbit average of the instantaneous fluxes, we express them using the QK parametrization in terms of (x, e, u), then evaluate the integral where P is the radial period and ⟨...⟩ stands for the average over one orbit.To compute this integral, we made use of the following formulas, with k and n ∈ N: ) is the usual normalized hypergeometric function.Note that in the case of k = 0, this integral reduces to the well known result where P n is the Legendre polynomials.
The nonspinning part of the orbit-averaged fluxes (instantaneous and tail contributions) were computed in Refs.[58][59][60]63], and we find agreement with their results to 2PN order.Our result for the 3PN spin contributions to the orbit-averaged instantaneous energy flux reads and for the angular momentum flux ) 15 (e 2 − 1) 9/2   The orbit average of the tail contribution to the fluxes is easier to compute since the only dependence on l is through the exponential, which averages to zero, leading to The eccentricity expansion in the tail part of the fluxes implies that it is valid for small eccentricities only.However, a simple resummation of the tail can make it accurate for high eccentricities; from Eqs. (3.10) and (3.11), we see that each PN order in the instantaneous part of the fluxes contains a factor 1/(1 − e 2 ) k , for some power k depending on the PN order, meaning that it diverges as e → 1.Such a divergence should also occur for the tail part of the fluxes, as argued in Ref. [104], based on the known closed-form quadrupole formula and how it relates to the LO tail [105].Therefore, pulling out a factor of 1/(1 − e 2 ) k from the tail is expected to improve the series expansion for high eccentricities. 1Reference [104] performed such a resummation of the tail and compared the analytical resummed expressions with numerical gravitational-self-force calculations and with analytical expansions to very high orders in eccentricities.Thereby, confirming the accuracy of the resummation that pulls a factor (1 − e 2 ) −5 at 1.5PN, (1 − e 2 ) −6 at 2.5PN, and (1 − e 2 ) −13/2 at 3PN in the energy flux.
Following this procedure, we advocate for using the following resummed form of the tail in the energy flux: Similarly, for angular-momentum flux, the resummation is given by ), plotted on a log scale and normalized by the circular-orbit limit.The lower panels show the relative errors between the numerical results of Ref. [60] and the analytical results.We see that the resummed analytical expressions are much closer to the numerical results than the eccentricity-expanded expressions, even for eccentricities near unity.As a check of this resummation, we compare the LO nonspinning part of Eqs.(3.14) and (3.15) with the numerical results for the tail integrals (2.12), which were calculated numerically in Refs.[58,60] to 3PN in the nonspinning case.Using the numerical values from Appendix C of Ref. [60], which are only given to four significant digits, we plot in Fig. 1 the LO fluxes, and compare them to the eccentricity-expanded and resummed analytical expressions above.We see that the eccentricity-expanded fluxes provide a decent agreement with the numerical results for small eccentricities, but are several orders of magnitude different for large eccentricities.However, the resummed fluxes show excellent agreement between the analytical and numerical results for all eccentricities, comparable to the numerical errors of Ref. [60] in the case of the resummed O(e 8 ) expansion, and less than 1% for the O(e 4 ) expansion.

⟨G
We also checked that the higher-order nonspinning contributions are similarly close to their numerical solutions after factoring (1 − e2 ) −6 at 2.5PN and (1 − e 2 ) −13/2 at 3PN in the energy flux, and factoring (1 − e 2 ) −9/2 at 2.5PN and (1 − e 2 ) −5 at 3PN in the angular momentum flux.While it would be good to compare the resummed spin part of the tail to numerical results, we consider such a calculation to be outside the scope of this paper, since it requires using a different method 2 from what we followed here to compute the fluxes.Despite not having numerical results for the spin contribution, we expect that, similarly to the nonspinning case, the resummation given in Eqs.(3.14) and (3.15) is also significantly more accurate than a pure eccentricity expansion.
As a different approach for resuming the hereditary contributions to the fluxes, Ref. [109] obtained superasymptotic and hyperasymptotic series that are accurate to better than 10 −8 , compared to the numerical PN results, for arbitrary eccentricities.Therefore, if an extremely high accuracy is required, one could use the resummation of Ref. [109], but we expect that the simple and physically-motivated resummation considered above is sufficient for most eccentric-orbit waveform models.

E. Evolution of the secular orbital elements
The balance equations relate the orbit-averaged fluxes to the energy and angular momentum losses, such that where the vectors J i and G i can be replaced by their magnitudes in the case of planar orbits.
Using the balance equations and the orbit-averaged energy and angular momentum fluxes from the previous subsection, we compute the secular evolution of the orbital elements x, a r , and e, due to radiation reaction.For example, and taking the orbit average leads to The evolution of other secularly varying quantities, such as e, a r and K, can be similarly related to the fluxes.
The 3PN nonspinning contributions were obtained in Refs.[60,63] in ADM and modified-harmonic coordinates, and we find agreement with their results to 2PN order.For the 3PN spin contributions in ⟨ ẋ⟩, we obtain where the SO term of order x 8 comes from the tail contribution to the fluxes while the rest are instantaneous contributions.
For ⟨ ė⟩, we obtain which is proportional to e.We included the results for ȧr and K in the Supplemental Material.
From the expressions for ẋ(x, e) and ė(x, e), one can solve for x(e), which is needed for the derivation of the memory contribution to the modes, as discussed in Subsection II D above.This is done by first obtaining dx/de = ẋ/ ė, then solving this differential equation order by order, in a PN and an eccentricity expansion, starting from an initial frequency x 0 and eccentricity e 0 .For example, at leading order in eccentricity, and at leading PN order for the nonspinning and spin contributions, we get The full expression is provided in the Supplemental Material.
We stress that all the results obtained in this section are for the radiated fluxes at infinity, but we did not account for the horizon absorption, which starts at 2.5PN for the spin contributions and at 4PN for the nonspinning contributions, relative to the leading order of the flux at infinity.For arbitrary mass ratios and quasi-circular inspirals, the horizon flux was derived in Ref. [110] at leading order, in Refs.[111,112] to 5PN for nonspinning and to 3.5PN for spinning binaries, and in Refs.[113][114][115] to 4PN.It was also derived to high PN orders in the test-mass limit for quasi-circular and eccentric orbits in Refs.[116][117][118][119][120][121][122].Once the energy and angular-momentum horizon fluxes are derived for eccentric orbits and arbitrary mass ratios, one can combine them with our results for the fluxes at infinity, and obtain their contributions to ẋ(x, e) and ė(x, e) (see Ref. [123] for recent work in that direction).

IV. 3PN WAVEFORM MODES FOR NONPRECESSING SPINS
In this section, we present our results for the 3PN spin contributions to the waveform modes, which are related to the radiative moments via Eqs.(2.5).We factor the modes such that and split H ℓm into instantaneous, tail, DC memory, and oscillatory memory contributions, i.e., where the instantaneous contributions are valid for generic motion, while the hereditary contributions are computed for bound orbits in a small-eccentricity expansion to O(e 6 ).To the order considered in this paper, there are no postadiabatic contributions to the waveform, which start in the nonspinning part at 2.5PN, as computed for eccentric orbits in Ref. [63], and start at 3.5PN for the spin contributions.
In this section, we write explicit expressions for the ℓ = 2 modes only and write eccentricity expansions to O(e 2 ).The full expressions for all modes up to the (6,6) mode, and to O(e 6 ) for the hereditary contributions, are included in the Supplemental Material [91].

A. Instantaneous contributions
The instantaneous contributions to the modes are computed using the instantaneous part of the radiative moments from Eqs. (2.10) and plugging them in Eq. (2.5), which yields the spin part of the modes to 3PN.
The 3PN nonspinning instantaneous contributions were derived in Ref. [61], and we find agreement with the results of that reference to 2PN order.For the spin contributions, our 3PN results agree with, e.g., Ref. [73] in the circular-orbit limit, and agree at 2PN for general orbits with Refs.[74,75].
For the (2,2) mode, we obtain the following instantaneous spin contributions Note that the (2,1) mode is the only mode containing cubic-in-spin terms for the 3PN waveform.Finally, we find for the (2,0) mode κ− + χ (2) Let us stress that the (ℓ, 0) modes do not vanish for eccentric orbits, unlike the case for circular orbits.The (2,0) mode in particular could be important for highly-eccentric binaries, since it is proportional to the eccentricity and starts at the same PN order as the (2,2) mode.

B. Tail contributions
The tail contributions to the modes are computed using the tail part of the radiative moments from Eqs. (2.12) and plugging them in Eq. (2.5), which yields the spin part of the modes to 3PN.The 3PN nonspinning tail contributions were derived in Ref. [62], and we checked that our results for the 2.5PN tail are in agreement.For the spin terms, we get agreement in the circular-orbit limit with Ref. [73].
To 3PN order, there are SO tail contributions in the (2,0), (2,1), (2,2), (3,0), and (3,2) modes, starting at 2.5PN for the (2,1) mode and at 3PN for the other modes.No SS tail contributions enter at this order, since they start at 3.5PN in the (2,2) mode.We write here the SO part of the tail to O(e 2 ) for the ℓ = 2 modes, and provide the full expressions to O(e 6 ) in the Supplemental Material.
In addition to the oscillatory memory, the (even ℓ, m = 0) modes contain a non-oscillatory (or DC) memory, which starts in the (2,0) mode at Newtonian order for the nonspinning part and at Our results for the memory agree with Ref. [63] in the nonspinning case, which we computed to 3PN for the oscillatory memory as a check, and to 2PN for the DC memory.For the spin part in the circular-orbit limit, we find agreement with the 2PN results of Ref. [79], which were computed using a different method, by relating the memory to the Moreschi supermomentum, which is equivalent to the energy flux at null infinity.For the 2.5PN and 3PN spin contributions, our results for the integrand of the DC memory, Eq. (2.15), agree in the circular-orbit and test-mass limits with Ref. [124], in which the authors computed the memory contribution analytically and numerically using gravitational-self-force methods.There is a difference between Ref. [124] and our results at 2.5PN for the integrated DC memory in the (2,0) and (4,0) modes, since we did not account for the horizon-flux contribution in ė(x, e), which is used when computing the DC memory integral (2.38).
To estimate the effect of the 2.5PN and 3PN spin terms derived here on the memory, compared to the lower orders in the circular-orbit case, we plot in Fig. 2 the (2,0) and (4,0) modes, scaled by their LO, as functions of x at each PN order for the spin part, while keeping the full 3PN nonspinning part.We consider two equal-mass configurations: one with spin magnitudes 0.9 aligned with the direction of orbital angular momentum, and the other with the same spin magnitudes but anti-aligned with the orbital angular momentum.We see that for large aligned spins, the NLO SO term at 2.5PN has a significant effect on the DC memory, while the 3PN SO and SS terms have a smaller effect.For anti-aligned spins, and also for small spin magnitudes, the effect of PN orders beyond the leading order is relatively small.Even though we did not include the horizon flux contribution when computing the DC memory, we checked that including it for circular orbits, using the results of Ref. [110], leads to a negligible effect on the plots in Fig. 2 for equal masses.

D. Phase redefinition for the full waveform
The waveform modes expressed as in Eq. (4.1), in terms of the orbital phase ϕ, depend on the arbitrary gauge constant b.The logarithmic dependence on b can be eliminated by a phase redefinition, or a shift in the coordinate time, as was done for circular orbits in Refs.[125][126][127].This procedure was generalized for eccentric orbits in Ref. [62] (see Sec. V.C. there), and we follow the same steps here, except for not including the post-adiabatic corrections, which are not needed for the 3PN spin contributions.
First, one redefines the mean anomaly l as The top panels are for the (2,0) mode, while the bottom ones are for the (4,0) mode.The left panels are for spins aligned with L, while the right panels are for spins anti-aligned with L, and all plots are for mass ratio q = 1.We see that the NLO SO term at 2.5PN can have a significant contribution to the memory, particularly for large aligned spins, while the 3PN spin terms have a smaller effect.
where M = M (1 − νx/2) + O(x 2 ) is the ADM mass, and x ′ 0 is related to b via Then, one defines the new phase ψ as where W (l) is the oscillatory part of the orbital phase ϕ = Kl + W (l).
In terms of ψ, the waveform modes can be expressed as where the nonspinning part H ψ,nS 22 is given at 3PN by Eq. ( 76) of Ref. [62].The full expressions to O(e 6 ) for all modes that have tail contributions are provided in the Supplemental Material.

V. CONCLUSIONS
In this paper, we completed the spin contributions in the radiative sector at 3PN order for eccentric orbits.We derived the energy and angular-momentum fluxes, as well as the waveform modes.In the fluxes, our results for the tail contributions complement the instantaneous ones derived in Refs.[72,77,78], and the 3PN nonspinning results of Refs.[58][59][60]63].We also computed the orbit-averaged fluxes and from them, the time evolution of the secular orbital elements, namely the frequency, eccentricity, semi-major axis, and periastron advance.For the modes, we obtained the instantaneous, tail, and memory (oscillatory and DC) contributions, which extends the 3PN nonspinning results of Refs.[61][62][63] and the 2PN spin results of Refs.[74,75].
To compute the hereditary contributions to the fluxes and modes, we first derived the quasi-Keplerian (QK) parametrization at 3PN for nonprecessing spins in harmonic coordinates using the covariant SSC, which complements the ADM-coordinates results of Refs.[85,86].Then, using the QK parametrization, we computed the hereditary contributions to the fluxes and modes in an eccentricity expansion for bound orbits, to O(e 8 ) in the fluxes and to O(e 6 ) in the modes.Thus, while the instantaneous contributions are valid for general planar orbits, the hereditary contributions are valid for small eccentricities.However, resuming the eccentricity expansion can lead to very good agreement with numerical calculations, as shown in Fig. 1 for the orbit-averaged fluxes.A similar resummation can also be performed to the hereditary contributions in the waveform modes.
The expressions for (r, ṙ, ϕ, φ) in terms of (x, e t , u), in addition to their eccentricity expansions in terms of (x, e t , l), are too lengthy to write here; we provide them in the Supplemental Material [91].

L
denote the instantaneous contributions, while U tail L and V tail L are the tail contributions.The memory contribution U mem L only appears in the mass-type moments.

1 FIG. 1 .
FIG.1.The leading-order nonspinning tail part of the energy flux (left panel) and angular-momentum flux (right panel), for eccentricity expansions to O(e 4 ) and O(e 8 ), plotted on a log scale and normalized by the circular-orbit limit.The lower panels show the relative errors between the numerical results of Ref.[60] and the analytical results.We see that the resummed analytical expressions are much closer to the numerical results than the eccentricity-expanded expressions, even for eccentricities near unity.

TABLE I .
Definition of the main quantities used in the QK parametrization. ) D. Resummation of the tail contributions ) − e 2 ) factor in the resummed expressions is such that it is consistent with the powers of those factors in the instantaneous part in Eqs.(3.10) and(3.11),and the rest of the expression is such that the O(e 8 ) expansion agrees with Eqs.(3.12) and (3.13).