Fourier domain gravitational waveforms for precessing eccentric binaries
Antoine Klein Yannick Boetzel Achamveedu Gopakumar Philippe Jetzer Lorenzo de Vittori
gr-qc astro-ph.CO







Department of Physics and Astronomy, The University of Mississippi, University, Mississippi 
38677, USA
CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade de 
Lisboa, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal
CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Bd Arago, 75014 Paris, 
France


Physik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich, 
Switzerland


Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research, 
Mumbai 400005, India


Physik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich, 
Switzerland


Physik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich, 
Switzerland



We build two families of inspiral waveforms for precessing binaries on eccentric orbits in the 
Fourier domain. To achieve this, we use a small eccentricity expansion of the waveform amplitudes 
in order to separate the periastron precession timescale from the orbital timescale, and use a 
shifted uniform asymptotics transformation to compute the Fourier transform in the presence of 
spin-induced precession. We show that the resulting waveforms can yield a median faithfulness above 
0.993 when compared to an equivalent time domain waveform with an initial eccentricity of e_0 
≈ 0.3. We also show that when the spins are large, using a circular waveform can potentially 
lead to significant biases in the recovery of the parameters, even when the system has fully 
circularized, particularly when the accumulated number of cycles is large. This is an effect of the 
residual eccentricity present when the objects forming the binary have nonvanishing spin 
components in the orbital plane.




 04.30.-w, 
 04.30.Tv 


Fourier domain gravitational waveforms for precessing eccentric binaries
    Lorenzo de Vittori
    December 30, 2023
========================================================================






§ INTRODUCTION


The recent discoveries of gravitational wave (GW) signals by the LIGO and Virgo Collaborations have
opened a new observation window on the Universe <cit.>, through which the potential for new discoveries in 
astrophysics is truly tremendous. So far, those events have been analyzed with the assumption that 
the systems that produced them were evolving on circular orbits. Indeed, it has been a well-known 
fact that the emission of gravitational waves by binary systems has the tendency to circularize 
their orbits <cit.>. Nevertheless, it has been argued that certain astrophysical 
scenarios could lead to stellar-origin black holes binaries having high initial 
eccentricities <cit.>, so they would 
still be measurable when the signal reaches the frequency window of the space-based GW detector 
LISA <cit.>. Furthermore, recent results have 
shown that eccentricity measurements by LIGO could be used to constrain stellar-mass black hole 
formation mechanisms <cit.>. It has been estimated that large biases in the recovery of the parameters of the 
first direct detection GW150914 could have occurred if the initial eccentricity in the detector was 
e_0 ≳ 0.05 <cit.> Supermassive black hole binaries could also have 
important eccentricities in the late inspiral, if triple systems are a significant ingredient of 
supermassive black hole evolution <cit.>. Furthermore, in some spin configurations, it has 
been shown that the eccentricity of the system never truly vanishes, but reaches a stationary value 
where it ceases to decrease through the emission of GWs <cit.>.

This has motivated the development of waveforms that include the effects of 
a nonzero eccentricity in GW binary signals. The first steps towards this goal rely on the 
derivation of quasi-Keplerian equations describing the orbits <cit.>, the 
derivation of the evolution equations for the orbital elements <cit.>, and the derivation of GW polarization 
amplitudes <cit.>. The effects of individual spins were later added to this 
approach <cit.>. 
Using these solutions, several waveforms have been developed. Yunes et al. <cit.> 
proposed an analytic eccentric waveform in the post-Newtonian (PN) postcircular approximation, by 
solving for the Fourier phase of a binary signal analytically at Newtonian order using a small 
eccentricity expansion. Cornish and Key <cit.> and 
Gopakumar and Schäfer <cit.> independently developed a numerical waveform in the 
time domain by solving the 1.5PN equations of motion numerically together with the spin-orbit 
precession equations, and using 1.5PN accurate amplitudes. Huerta et al. <cit.> 
expanded the analytical work of Yunes et al. by including the most important eccentricity-dependent 
terms up to 3.5PN order and at eighth order in the initial eccentricity for non-spinning systems. 
Tanay et al. <cit.> later computed the full 2PN Fourier phase for nonspinning systems 
at second order in the eccentricity. Moore et al. <cit.> then expanded this result to 
3PN order. Huerta et al. <cit.> and Hinder et al. <cit.> 
combined those results with numerical relativity to produce an eccentric inspiral-merger-ringdown 
waveform for nonspinning binaries. Recently, Hinderer and Babak <cit.> and Cao and 
Han <cit.> developed an eccentric waveform using a new approach in the effective one-body 
(EOB) formalism.

In this work, we further develop the formalism of post-Newtonian eccentric waveforms to include the 
effects of spin-induced precession in the Fourier domain. The advantage of Fourier domain waveforms 
over time domain ones is that they provide a much more computationally efficient way of computing a 
GW signal. Indeed, in order to produce a time domain waveform, one has to compute an equally spaced 
time series of the signal before computing its Fourier transform to use in detection or parameter 
estimation algorithms. The relevant timescale for this time series is the inverse of the maximum 
orbital frequency, which, being very short, makes this process computationally very expensive. 
Having a waveform available directly in the Fourier domain circumvents this problem and greatly 
reduces the computational cost of GW data analysis. In order to construct such a waveform, we solve 
the evolution equations for the orbital elements together with the orbit-averaged spin precession 
equations numerically at 3PN order, including spin effects at 2PN order. Using a quasi-Keplerian 
description of the orbit, we employ instantaneous nonspinning amplitudes to construct the resulting 
GW polarizations. We then use a shifted uniform asymptotics (SUA) technique <cit.> to 
compute the waveforms in the Fourier domain. The resulting waveform has the advantage that the 
phasing is computed without any expansion for small eccentricities and thus can be be very 
faithful compared to corresponding time domain waveforms for moderate to large eccentricities (e 
≲ 0.4). However, the amplitudes require a small-eccentricity expansion, and thus we do not 
expect the present waveforms to be faithful for 
arbitrarily large eccentricities.

In Sec. <ref>, we derive two different families of eccentric waveforms. Due to the 
similarity between the orbital timescale and the periastron-to-periastron timescale, we derive the 
first family by expanding the Fourier domain waveform into combined harmonics of the mean orbital 
phase and of the mean anomaly. We then derive the second family by further expanding the resulting 
Fourier phase and time-frequency relations for small differences between the two similar phases.

In Sec. <ref>, we describe simulations that we performed to compute the faithfulness 
between our Fourier domain waveforms and a corresponding time domain waveform, including a detailed 
summary of how these different waveforms are constructed. We also compare a circular waveform to 
probe which domain of the parameter space allows for such circular waveforms to be effectively used 
for parameter estimation of binary signals. 

We give concluding remarks in Sec. <ref>. Throughout this paper, we use geometric 
units where G=c=1.




§ WAVEFORM


In the presence of spins, the orbit of a binary system is, in general, not restricted to an orbital 
plane <cit.>. Indeed, interactions between the spins and the orbit cause them to 
precess. However, in the post-Newtonian regime, the timescale on which this precession occurs is 
well separated from the other timescales present in the problem. We can therefore approximate the 
spin-orbital precession to be occurring much more slowly than the orbit, which allows us to 
describe it using a so-called quasi-Keplerian parametrization inside an orbital plane that stays 
perpendicular to the orbital angular momentum as the latter precesses. A quasi-Keplerian 
parametrization of the orbit of a spinning binary system is known at 3PN order for the 
nonspinning part <cit.>, and at 2PN order for the 
spin-dependent part <cit.>. In this work, we restrict the quasi-Keplerian orbital 
description at 2PN for the computation of the polarization amplitudes. We can express the orbit at 
2PN order as


    r    = a ( 1 - e_r cos u) +f_r(v)  , 
    ϕ   = (1 + k) v + f_ϕ(v)  , 
    tanv/2   = √(1 + e_ϕ/1 - e_ϕ)tanu/2 , 
    
    	l    = u - e_t sin u + f_t(u,v)  , 
    l̇   = n  ,


where (r, ϕ) is a polar representation of the separation vector in the orbital plane, a is 
the semimajor axis; u is the eccentric anomaly; v is the true anomaly; l is the mean 
anomaly; n is the mean motion; e_r, e_ϕ and e_t are eccentricity parameters; and the 
functions f_i are general relativistic corrections given 
by <cit.>


    f_r(v) =    ∑_i=0^2 b_r,icos (2v - 2ψ_i)  , 
    
    	f_ϕ(v) =    ∑_k=2^3 a_ϕ,ksin(k v) 
       + ∑_k=1^2 ∑_i=0^2 b_ϕ,k,isin(k v - 2ψ_i)  ,
    
    	f_t(u,v) =     g_t (v - u) + a_tsin(v)  ,



where ψ_i is the angle between the periastron line and the projection of spin i onto the 
orbital plane (see Fig. <ref>), ψ_0 = (ψ_1 + ψ_2)/2 and the constants a_A, 
b_A and g_t are listed in Appendix <ref>. We complemented the spinning solution of 
 <cit.> by including quadrupole-monopole terms as described in Appendix 
 <ref>.



The orbital phase ϕ and the mean anomaly l can be decomposed as the sum of a linearly 
growing part and a periodic part <cit.>,


    ϕ   = λ + W_ϕ , 
    λ̇   = (1 + k)n  , 
    l̇   = n  , 
    
    	W_ϕ   = (1 + k) (v - l) + f_ϕ(v)  .


We choose to express our equations in terms of the post-Newtonian parameter y and the 
eccentricity parameter e defined by


    y    = [ M (1 + k) n ]^1/3/√(1 - e_t^2) , 
    
    	e    = e_t  .


The constants in the quasi-Keplerian parametrization are given in terms of these parameters in 
Appendix <ref>.

As the system emits gravitational waves, the orbital frequency and the eccentricity will evolve 
according to the following equations <cit.>


    M d y/dt   = ( 1 - e^2	)^3/2ν y^9 ( a_0 + ∑_n=2 a_n y^n 
    		)  , 
    
    	M de^2/dt   = - ( 1 - e^2 )^3/2ν y^8 ( b_0 + ∑_n=2 b_n y^n 
    		)  ,


where the constants a_i and b_i are given at 3PN order for nonspinning systems and at 2PN 
order for spinning systems in Appendix <ref>. Here we also complemented the spinning 
solution of <cit.> by including quadrupole-monopole terms as described in 
Appendix <ref>. We found that the minimum value for the eccentricity e^2min 
found in <cit.> is unchanged by the addition of quadrupole-monopole effects, with

    e^2min   = 5 y^4/304σ(-1, 1, 0, 2, -2, 0) 
       = 5 y^4/304| s_⊥^(-)|^2  ,

where the 2PN spin-spin coupling σ can be found in Appendix <ref>, and

    s_⊥^(-)   = [ s_1 - L( L·s_1 ) ] 
    		- [ s_2 - L( L·s_2 ) ] ,

where L is a normal to the orbital plane.

Note that we found a typo in <cit.>, where the constant factor in e^2min should 
read 5/304 instead of 5/340. This minimum eccentricity depends on the spin orientations: it is 
maximal when the spins lie inside the orbital plane and are opposite to one another, and it vanishes 
when the projections of s_1 and s_2 onto the orbital plane are equal to each other. 
The maximum value it can take is independent of the mass ratio; it is e^2min = 5 y^4/304, 
which evaluates to emin≈ 0.021 at the ISCO defined by y = 6^-1/2, and it is 
multiplied by a factor (f/fISCO)^2/3 earlier in the inspiral. Note that this minimum 
eccentricity, being a spin effect, is unrelated to a similar effect observed in extreme mass-ratio 
inspirals around Schwarzschild black holes in <cit.>, and also unrelated to another 
effect due to orbital effects derived in <cit.>, which is of order e^2min∼ 
y^10 and is independent of the spins.

The 2PN orbit-averaged equations of precession are given by <cit.>


    M L   = - ( 1 - e^2 )^3/2 y^6 ( Ω_1 + Ω_2 ) 
    		 , 
    
    	Ms_1    = ( 1 - e^2 )^3/2μ_2 y^5 Ω_1  , 
    
    	Ms_2    = ( 1 - e^2 )^3/2μ_1 y^5 Ω_2  ,


where we defined the reduced spins s_i = S_i / m_i, the reduced individual masses 
μ_i = m_i/M, and the precession vectors Ω_i are given by

    Ω_i =    {[ 2 μ_i + 3/2μ_j - 3/2 yL·( s + (q_i - 1) s_i ) ] L
       + 1/2 y s_j }×s_i  ,

where i,j ∈{1, 2}, i ≠ j, and the q_i are quadrupole parameters defined in such a way that q_i=1 for black holes.

The gravitational waveform emitted by a binary system on such an orbit has been computed at 3PN 
order for nonspinning binaries, omitting tail effects <cit.>. The result has the 
following structure:


    h(t) =     F_+ (t) h_+ (t) + F_× (t) h_× (t)  ,
    
    	h_+,×(t) =    ∑_n ∈ℤ H_+,×^(n) (y, e, e cos u, e sin u) 	
    		e^in (ϕ + ϕ_T) ,


where F_+ and F_× are antenna pattern functions <cit.>, and the Thomas 
precession phase ϕ_T is given by <cit.>

    ϕ̇_T    = L·N/1 - (L·N)^2( L×N) ·L ,

with respect to a given sky location vector N.

In order to compute the Fourier transform of this signal, we need to separate the orbital timescale 
from the precessional one, and express the orbital timescale dependence in terms of linearly 
growing phases. To do so, we follow <cit.> and compute an inversion of the 
PN-accurate Kepler equation (<ref>) as

    u    = l + ∑_s=1^∞ A_s sin(sl)  ,

with the Fourier coefficients A_s given by

    A_s    = 2/s J_s(s e) + ∑_j=1^∞α_j [ J_s+j (s e) - 
    		 J_s-j(se)]  .

The PN-accurate constants α_j can be computed from <cit.> and are 
given in Eq. (18) of <cit.>. Similarly, we can find a Fourier expansion of the true 
anomaly v and the orbital phase ϕ in terms of the mean anomaly l:


    v    = l + ∑_s=1^∞ B_s sin(sl)  ,
    ϕ   = λ + ∑_s=1^∞ C_s sin(sl)  .


The Fourier coefficients A_s, B_s and C_s can be found up to y^4, e^5 in 
Appendix <ref>. Using this solution, we can then express


    e^i k u   = ∑_p ∈ℤϵ_p^ku  e^-i p l , 
     
    	e^i k v   = ∑_p ∈ℤϵ_p^kv  e^-i p l ,
    
    	e^i k ϕ   = e^i k λ∑_s ∈ℤ𝒫_p^kϕ  e^-i p l ,


where the coefficients ϵ_p^ku, ϵ_p^kv and 𝒫_p^kϕ are 
given as a Taylor expansion in both e and y. We refer to Eqs. (30), (34), (E11) 
of <cit.> for how to calculate these Fourier coefficients.

This small eccentricity expansion allows us to express the waveform as

    h_+,×(t)    = ∑_n∈ℤ∑_p∈ℤ H_+,×^(p,n) 
    		e^-i (n λ + p l) ,

where we included the Thomas phase ϕ_T into the amplitudes H_+,×^(p,n), which vary 
on the spin-precession timescale. To separate the periastron precession timescale from the orbital 
timescale, we define δλ = λ - l, such that


    M λ̇   = (1 - e^2 )^3/2 y^3  , 
    δλ̇   = k/1 + kλ̇ .


This new angle defines the periastron precession timescale, which is similar to the spin precession 
timescale since δλ̇/λ̇ = y^2.

Using this, we can then further simplify the waveform with


    h_+,×(t)    = ∑_n ∈ℤ H_+,×^(n) e^-i n λ ,	
    		
    
    	H_+, ×^(n)   = M ν y^2/d_L∑_m ∈ℤ G_+, ×^(m,n) 
    		e^-i m δλ e^-i(n+m)ϕ_T .


The amplitudes G_+,×^(m,n) are given in Appendix <ref> at order y^2, e 
[A Mathematica version of all amplitudes to order y^4, e^10 is available as 
supplemental material or upon request from 
mailto:boetzel@physik.uzh.chboetzel@physik.uzh.ch.].




 §.§ Fourier transform approximations


Before we compute an approximation of the Fourier transform of our signal, let us introduce 
two useful techniques.

Let us first assume that we have a signal of the form

    h(t)    = A(t) e^-i ϕ(t) ,

with ϕ̇(t) a positive and monotonically increasing function of time, and that we want to 
compute its Fourier transform

    h̃(f)    = ∫_-∞^∞ h(t) e^2π i f t dt  .


The stationary phase approximation (SPA) of this Fourier transform consists in Taylor expanding the 
amplitude A(t) and phase ϕ(t) around the stationary point t_0 defined by the relation 

    2 π f    = ϕ̇(t_0)  ,

keeping only the constant term in the expansion of the amplitude and up to the quadratic order in 
the expansion of the phase. We get

    h(t) ≈     A(t_0) exp{-i [ ϕ(t_0) 
       + ϕ̇(t_0) (t - t_0) + 1/2ϕ̈(t_0) (t - t_0)^2 ] } .

We can compute the Fourier transform of this approximate signal analytically, and we get

    h̃(f)    ≈√(2 π/| ϕ̈(t_0) |^2) A(t_0) e^i (2 π f t_0 - 
    		ϕ(t_0) - π/4) .

This approximation will be accurate if |Ȧ/A| ≪ |ϕ̈|^1/2 around the stationary point, and if the quadratic approximation is accurate around the stationary point. For a formal derivation, see e.g. <cit.>.

Let us now suppose that our signal is of the form

    h(t)    = A e^-i [ϕ_C + B sinβ] ,

with Ȧ/A = y^8, ϕ̇_C = y^3, B = y, β̇ = y^5, 
and that each additional time derivative adds a factor y^8 to the various quantities 
present in the signal, with y a small expansion parameter. This is the simplified form of a GW 
signal that we expect from a binary system undergoing spin-induced orbital precession, with y 
being a PN expansion parameter. The SPA cannot be directly used in this case, because the two terms 
in the second time derivative of the signal phase

    ϕ̈   = ϕ̈_C - B β̇^2 sinβ + y^13

are of the same PN order and can cancel each other. The shifted uniform asymptotics (SUA) 
method <cit.> offers an approximation of the Fourier transform of such a signal by first 
expanding the signal using Bessel functions as

    h(t)    = A ∑_k J_k(B) e^-i(ϕ_C + k β) ,

so its Fourier transform can be approximated by a series of SPA, since β̈≪ϕ̈_C. The Fourier transform then becomes

    h̃(f) ≈    ∑_k A √(2 π/ϕ̈_C + k β̈)
       ×exp[ i ( 2 π f t_k - ϕ_C - k β - π/4) 
    		 ]  ,

where the various functions are evaluated for each k ∈ℤ at the stationary time t_k defined by 

    2 π f    = ϕ̇_C(t_k) + k β̇(t_k)  .


The different stationary times can be related to each other by Taylor expanding their defining 
equations around t_0 and solving for the difference order by order:

    t_k - t_0    = - k β̇(t_0)/ϕ̈_C(t_0) + y^-4 .

By Taylor expanding Eq. (<ref>) around t_0, and keeping only the leading PN order 
amplitude and the phase accurate to order y^0, we obtain


    h̃(f) ≈    h̃_0(f) h̃corr(f)  , 
    h̃_0(f) =    √(2 π/| ϕ̈(t_0) |^2) A(t_0) e^i (2 π f t_0 - 
    		ϕ_C(t_0) - π/4) , 
    h̃corr(f) =    ∑_k J_k [B(t_0)] exp[ - k β(t_0) + 1/2 T^2 k^2 
    		β̇^2(t_0) ], 
    
    	T =    √(1/ϕ̈_C(t_0)) .


After some manipulation, we can resum the Bessel functions in h̃corr(f) as

    h̃corr (f)    = ∑_p ≥ 0( - i T^2)^p/2^p p!∂_t^2p e^-i B sinβ ,

where the functions are evaluated at t=t_0. Truncating this series at some order p = kmax 
and using a stencil around t_0 to approximate the different order time derivatives, we obtain the 
SUA approximation

    h̃corr (f)    ≈∑_k=-kmax^kmax a_k,kmax e^-i B 
    	sinβ(t_0 + k T) ,

where the constants a_k,kmax satisfy the following linear system of equations:



    ∑_k=1^kmax a_k,kmax + 1/2 a_0, kmax = 1  , 
       ∑_k = 1^kmax a_k, kmaxk^2p/(2p)! = (-i)^p/2^p p! , 
    		   p ∈{1, …, kmax} , 
       a_-k,kmax = a_k,kmax .



To summarize, if we are able to separate the spin-precessional timescale effects from
a carrier phase ϕ_C that satisfies ϕ̇_C > 0 and ϕ̈_C > 0 as

    h(t)    = A(t) e^-i ϕ_C(t) ,

where all spin-precessional timescale effects are included in A(t), then we can 
write the SUA approximation of its Fourier transform:

    h̃(f) =    √(2π/ϕ̈_C) e^i(2π f t_0 - ϕ_C(t_0) - π/4)
       ×∑_k = -kmax^kmax a_k, kmax A(t_0 + k T)  ,

with the constants a_k,kmax satisfying the linear system of Eqs. (<ref>), and T 
= [ϕ̈_C(t_0)]^-1/2.




 §.§ Periastron precession effects


Let us first derive a waveform in the Fourier domain taht is valid for nonprecessing spins, and add the 
effects of spin-orbit precession later. Putting aside spin-orbit precession, the signal in the time 
domain can be expressed as in Eqs. (<ref>):

    h_+,×(t)    = M ν y^2/d_L∑_n ∈ℤ∑_m ∈ℤ 
    		G_+,×^(m,n) e^- i (n λ + m δλ) .

Using the SPA, we can approximate its Fourier transform by


    h̃_+,× (f) =    ∫ h_+,×(t) e^2 π i f t dt 
    
    		=    M ν y^2/d_L∑_n ≥ 1∑_m ∈ℤ√(2 π/n 
    		λ̈ + m δλ̈) G_+,×^(m,n)
       × e^i [2 π f t_n,m - n λ(t_n,m) - m δλ (t_n,m) - 
    		π/4] , 
    
    	2π f =     n λ̇(t_n,m) + m δλ̇(t_n,m)  ,


where each of the harmonics (n,m) has to be evaluated at a different time. It is worth noting 
here that we assumed that n λ̇ + m δλ̇ > 0, which is not necessarily 
true for every (n, m) pair during the whole inspiral. However, for this assumption to break down, 
one needs negative and sufficiently large m, since δλ̇ / λ̇ = 
y^2, and the corresponding amplitude will be suppressed by a factor e^m. We verified that 
ignoring this fact does not lead to high inaccuracies, at least for initial eccentricities e_0 
≲ 0.4.

In order to simplify the expression of the Fourier domain waveform and to improve its computational 
efficiency, we look for an expression of the following form:


    h̃_+,×(f)    = M ν y^2/d_L∑_n ≥ 1h̃_n,0(f) ∑_m ∈ℤh̃_n,mPP (f)  ,
    h̃_n,0(f)    = √(2 π/n λ̈) e^i[2 π f t_n - n λ(t_n) - 
    		π/4] , 
    
    	2 π f    = n λ̇(t_n)  ,


where h̃_n,0(f) is a waveform harmonic without any periastron precession effects and 
h̃_n,mPP(f) are corrections to it. In order to evaluate 
h̃_n,mPP (f), we define

    Δ t_n,m   = t_n,m - t_n  ,

and Taylor expand Eq. (<ref>) around t_n:

    2π f    = ∑_p ≥ 0Δ t_n,m^p/p!. d^p/dt^p(n 
    		λ̇ + m δλ̇) |_t=t_n .

We can use this together with Eq. (<ref>) to solve for the PN expansion of Δ 
t_n,m order by order, and we obtain

    Δ t_n,m   = ∑_p = 1^P 1/p!(-m/n)^p D^p-1( 
    		δλ̇^p/λ̈)  ,

where the differential operator D is given by

    D    = 1/λ̈d/dt ,

and every function of time is evaluated at t=t_n. We have checked that this expression remains 
valid at least up to P=6.

Using this, we can then Taylor expand the phase in Eq. (<ref>) around t_n to compute


    h̃_n,mPP (f)    = G_+, ×^(m,n) e^i ΔΨ_n,m ,
    ΔΨ_n,m   = - m δλ + n ∑_p = 2^P+11/p!( -m/n)^p D^p-2( δλ̇^p/λ̈)  ,


where all functions are once again evaluated at the stationary time t_n defined by Eq. (<ref>), and we checked that the latter equation is valid at least up to P=6. 
Equations (<ref>) and (<ref>) are PN expansions in the sense that each increasing 
 order in m is multiplied by a factor of PN order (δλ̇/λ̈) 
(d/d_t) = y^2, as both δλ̇ and λ̈ evolve on the radiation 
reaction timescale. This implies that the formal expansion in m in these two equations coincides with a PN expansion at order 2P beyond leading order.




 §.§ Complete waveform


We can now add spin precession by using a SUA transformation <cit.> instead of a SPA. We 
start by noting that we can express the waveform in the time domain by

    h(t)    = ∑_n,m𝒜_n,m(t) e^-i(nλ + mδλ) ,

where all spin-precession timescale effects are included in the amplitudes 

    𝒜_n,m(t) =    M ν y^2/d_L[ F_+(t) G_+^(n,m)(t) 
       + F_×(t) G_×^(n,m) (t) ] e^-i (n + m) ϕ_T .

This allows us to directly use a SUA transformation. If we restrict the amplitudes to y^N, 
e^M, we then obtain


    h̃(f) =    ∑_n = max(1, 2-N)^2+N∑_m = -M^Mh̃_n,m (f)  , 		
    		
    h̃_n,m(f) =    h̃_n,m^(0)(f) h̃_n,mSP(f)  , 
    h̃_n,m^(0)(f) =    √(2π) T_n,mexp[i ( 2 π f t_n,m
        -n λ (t_n,m) - m δλ (t_n,m) - π/4 ) ]  , 
    
    	2 π f =     n λ̇ (t_n,m) + m δλ̇(t_n,m)  , 
    
    	T_n, m =    [ n λ̈ (t_n,m) + m δλ̈ (t_n,m) 
    		]^-1/2 , 
    
    	h_n,mSP (f) =    ∑_k = - kmax^kmax a_k, kmax𝒜_n,m(t_n,m + k T_n,m)  ,


where the constants a_k, kmax satisfy the linear system of equations defined in 
Eq. (<ref>). This waveform is in the Fourier domain and consistently includes the effects 
of spin-induced precession and periastron precession. As we will see in the next section, it allows 
for large matches with time domain waveforms with eccentricities e ≲ 0.3 that we can 
consider as moderate in the modeling sense, because only the amplitudes, not the phasing, have been expanded for small 
eccentricities.

The waveform defined by Eq. (<ref>) suffers from the fact that it includes a double sum, and 
therefore its computational cost rises quickly as the precision of the amplitudes increases. 
However, in order to increase its computational efficiency, we can use a similar strategy as 
described in the previous subsection and expand the waveform in powers of 
δλ̇/λ̇.

First, we can approximate the SUA timescale in Eq. (<ref>) by T_n,m≈ T_n = 
T_n,0. Next, we can use Eqs. (<ref>) and (<ref>) to define Δ 
t_n,m and ΔΨ_n,m at order P:



    Δ t_n,m =    ∑_p = 1^P1/p!( -m/n)^p D^p-1( 
    		δλ̇^p/λ̈)  , 
    ΔΨ_n,m =     - m δλ + n ∑_p = 2^P+11/p!(- 
    		m/n)^p D^p-2( δλ̇^p/λ̈) .


We can use Eqs. (<ref>) and (<ref>) together with the chain rule

    d/dt f(y, e^2)    = ∂ f/∂ ydy/dt + ∂ f/∂ e^2de^2/dt ,

to get the necessary derivatives of λ and δλ as PN expanded functions. Thus, 
we can simplify the waveform as


    h̃(f) =    ∑_n = max(1, 2-N)^2+Nh̃_n^(0) (f) h̃_nPP(f) 
    		 , 
    h̃_n^(0)(f) =    √(2π)  T_n exp[i ( 2π f t_n - n λ (t_n) - 
    		π/4) ]  , 
    
    	2 π f =     n λ̇ (t_n)  , 
    
    	T_n =    [ n λ̈ (t_n) ]^-1/2 , 
    h̃_nPP (f) =    ∑_m = - M^M e^i ΔΨ_n,m
       ×∑_k = - kmax^kmax a_k, kmax𝒜_n,m (t_n + 
    		Δ t_n,m + k T_n)  .


Equation (<ref>) presents a further expanded waveform, and can possibly be made more efficient 
than the one defined by Eq. (<ref>), especially for amplitudes of high (N,M) order. Thus 
we get a family of Fourier domain gravitational waveforms for spin-precessing binaries on eccentric 
orbits characterized by the expansion orders (P, N, M).




§ COMPARISONS


We have run different sets of simulations in order to probe under what circumstances our waveforms 
defined in Eqs. (<ref>) and (<ref>) are sufficiently faithful to equivalent waveforms 
obtained in the time domain. For all the waveforms used in our comparisons, we use nonspinning 
amplitudes at 2PN order omitting tail terms <cit.>, and we use evolution equations for 
y and e^2 at 3PN nonspinning order and 2PN spinning order, including tail terms, as described 
in Appendix <ref>. For all Fourier domain waveforms, we use a SUA transformation as 
in <cit.> with kmax = 3.

We use as a reference time domain waveform h_R obtained in the following way:

	
  * Equations (<ref>)-(<ref>) are solved numerically together with 
		Eqs. (<ref>), (<ref>), and (<ref>) in order to 
		yield solutions for y(t), e^2(t), L(t), s_1(t), s_2(t), 
		λ(t), ϕ_T(t), and δλ(t), from an initial time t_0 to a maximum 
		time tmax defined by the ISCO-like condition M λ̇(tmax) = 
		6^-3/2.
	
  * Time is equally sampled between t_0 and tmax using a sampling time Δ t = 
		2π/[24 λ̇(tmax)], in order to ensure that the first 12 waveform 
		harmonics fall below the Nyquist frequency. Equations (<ref>)-(<ref>) are solved 
		at each step to get the orbital phase ϕ and the eccentric anomaly u. 
		Equation (<ref>) is inverted numerically to yield u( l = λ - δλ, e = 
		e_t).
	
  * A waveform signal is constructed using Eqs. (<ref>) and (<ref>), and 	
		the solutions for y(t), e(t), ϕ(t), ϕ_T(t), and u(t). The antenna pattern 
		functions are chosen in the low-frequency limit, for a static detector <cit.>. 
		The waveform amplitudes are included at 2PN order, with the omission of spin effects and 
		tail terms.
	
  * A Tukey window is introduced in order to reduce spectral leakage and a discrete Fourier 	
		transform of the signal is taken to yield the waveform in the Fourier domain.


We compare different waveforms to the reference one

	
  * A nonexpanded eccentric one (NE_M) defined by Eq. (<ref>) and N=4, M ∈{0, …, 6}, i.e. with amplitudes at N/2=2PN order and amplitudes expanded at Mth 
	order in e.
	
  * An expanded eccentric one (EE_M,P) defined by Eq. (<ref>) and N=4, M ∈{0, …, 6}, P ∈{0, …, 3}, i.e. with amplitudes at N/2=2PN order, 
	amplitudes expanded at Mth order in e, and the waveform expanded at Pth order in 
	δλ̇ / λ̇ as in Eqs. (<ref>).
	
  * A circular one (C) with amplitudes at 2PN order, taken from <cit.>.

Note that Eqs. (<ref>) imply that the waveforms NE_0 and EE_0,P are identical 
for any P.

In order to make our comparisons, we compute the faithfulness F = maxℳ, defined by 
the match ℳ maximized over some of the waveform parameters, with


    ℳ   = (h, h_R)/√((h,h)(h_R,h_R)) , 
    
    	(a,b)    = ∫_fmin^fmaxã(f) b̃^*(f) df  ,


where we chose a white detector noise in order to make as few assumptions about the detector as 
possible. For the eccentric waveforms, since they use the same phasing as the reference one, we do 
not maximize over any parameters, while for the circular waveform, we maximize the match over time 
and orbital phase shifts to obtain the faithfulness. We compare the faithfulness obtained this 
way to a fiducial value of F = 0.993, corresponding to a faithfulness level at which we can 
estimate that the errors in the recovered parameters due to mismodeling are smaller than the 
statistical errors coming from the detector noise, for D=10 intrinsic parameters and a 
signal-to-noise ratio (SNR) of ρ = 25 <cit.>. The relation between the 
faithfulness and the SNR at which the mismodeling error becomes likely to exceed the statistical 
error in a GW detection is <cit.>

    F    ≈ 1 - D-1/2 ρ^2 .


We ran two different sets of simulations: one to study systems in the late inspiral as 
observed by the LIGO/Virgo network and by LISA in the case of massive black hole binaries 
[denoted by (Xa)], and the other to study systems in the early inspiral as observed by LISA for 
stellar-origin black hole binaries [denoted by (Xb)].
We made six different runs in order to probe the faithfulness of our waveforms as a function of 
the starting eccentricity in different situations:

	
  (I) We randomize the initial eccentricity with a log-flat distribution 10^-5 < e_0 < 	
		0.5 and the spin magnitudes with a flat distribution 0 < χ_i < 1.
	
  (II) We randomize the initial eccentricity with a log-flat distribution 10^-5 < e_0 < 
		0.5 and spin magnitudes with a flat distribution 0 < χ_i < 0.1.
	
  (III) We randomize the initial eccentricity with a log-flat distribution 10^-5 < e_0 
		< 0.5 and spin magnitudes set to the maximum value χ_i = 1.
	
  (IV) We start with zero initial eccentricity and spin magnitudes with a flat 	
		distribution 0 < χ_i < 1.
	
  (V) We start with zero initial eccentricity and spin magnitudes with a flat 	
		distribution 0 < χ_i < 0.1.
	
  (VI) We start with zero initial eccentricity and spin magnitudes set to the maximum 
		value χ_i = 1.

We thus have twelve runs: (Ia)–(VIa) in the late inspiral case and (Ib)–(VIb) in the early inspiral 
case.

To get the binary parameters used in our runs, we randomize all vector directions with a flat 
distribution on the sphere. Since the distance does not affect the match ℳ in 
Eq. (<ref>), we fix it at some fiducial value. We randomize the initial orbital phase and 
the initial periastron-ascending node angle ϕ_e (see Fig. <ref>) with a flat 
distribution in [0, 2π]. Whenever the randomized initial eccentricity is lower than the minimal 
value given in Eq. (<ref>), we set e_0 = emin. Note that cases (IV) to (VI) 
correspond to fully circularized binaries, but Eq. (<ref>) prevents them from having truly 
zero eccentricity unless the reduced spins have exactly equal support in the orbital plane. In each 
late inspiral run, we start our simulations with an initial eccentricity e_0 and at a frequency 
M λ̇start = 6^-3/2/10, and stop at M λ̇stop = 6^-3/2. We 
also randomize the mass ratio q = m_2/m_1 with a log-flat distribution between 1 and 1/30, 
and use a fixed total mass M = 100 M_⊙, taking advantage of the white detector noise.
In each early inspiral run, we randomize the two masses with a log-flat distribution 10 M_⊙ < 
m_i < 100 M_⊙. We then use the Newtonian time-frequency relation and the initial eccentricity 
to determine the starting frequency such that the system will evolve to have an orbital frequency 
of fend = 1 Hz after T = 4 yr:


    ystart   = [ 5 M yend^8/5M + 32 ν T (1 - e_0^2)^3/2(8 + 7 e_0^2) yend^8]^1/8 , 
    
    	yend   = [ 2 π M fend/(1- e_0^2 )^3/2]^1/3 .


We then let the system evolve and stop after four years, and set the maximum frequency 
fmax=1 Hz in Eq. (<ref>).




 §.§ Late inspiral systems




We present in Fig. <ref> the results from late inspiral run (Ia), with starting 
eccentricity 10^-5 < e_0 < 0.5 and spin magnitudes 0 < χ_i < 1. In it, we compare the mean 
faithfulness as a function of the initial eccentricity for different waveforms. The top panel 
shows a comparison between the results for the circular waveform C, the nonexpanded eccentric 
waveform NE_6, and the expanded eccentric waveforms EE_6,P, 0 ≤ P ≤ 3, and the bottom 
panel shows a comparison between the expanded eccentric waveforms EE_M,2, 0 ≤ M ≤ 
6. We can see in the top panel that the circular waveform stays above the fiducial faithfulness 
only for initial eccentricities of e_0 ≲ 0.008. Furthermore, the results for the
expanded eccentric waveform become very close to the nonexpanded version starting at second order in 
δλ̇/λ̇, and leads to a faithfulness above the fiducial threshold for 
eccentricities below e_0 ≲ 0.3. On the bottom panel, we can see the effects of the 
expansion of the waveform amplitudes for small eccentricities. We can see that the largest starting 
eccentricity for which the median faithfulness stays above the threshold increases with increasing 
order in the expansion. Furthermore, we can see that below a certain starting eccentricity 
depending on the specific order, increasing the expansion order has no effect on the faithfulness, 
as the errors due to this approximation become subdominant.



We present in Fig. <ref> the results from late inspiral run (IIa), with starting 
eccentricity 10^-5 < e_0 < 0.5 and spin magnitudes 0 < χ_i < 0.1. These results are very similar 
to the results of run (Ia), but due to the reduced spin magnitudes the starting eccentricities 
reach smaller values. On the top panel, we can see that below a starting eccentricity of e_0 
≲ 10^-3, the loss of faithfulness using circular waveforms with respect to our eccentric 
models becomes negligible.



We present in Fig. <ref> the results from late inspiral run (IIIa), with starting 
eccentricity 10^-5 < e_0 < 0.5 and spin magnitudes χ_i = 1. The results are similar to the ones shown in Figs. <ref> and <ref>, but the increased magnitudes of the spins slightly reduce the performance of the circular waveform.
Comparing this figure to Figs. <ref> and <ref>, we can conclude that 
the value of the spin magnitudes has little effect on the faithfulness, other than on  the limiting residual eccentricity. 


 
We present in Fig. <ref> the results from late inspiral run (IVa), with starting 
eccentricity e_0 = emin and spin magnitudes 0 < χ_i < 1. We can see here an 
effect due to the residual eccentricity. Indeed, the circular waveform performs poorly in some 
cases, even when the binaries are fully circularized. In our simulations, 7% of the 
faithfulness for the circular waveform was below the threshold line, while virtually no 
faithfulness were found below it for waveforms that used eccentric phasing, even with the lowest 
order amplitudes. While this does not represent a large proportion of binaries, this number will 
only increase when considering binaries with higher SNRs.


 
We present in Fig. <ref> the results from late inspiral run (Va), with starting 
eccentricity e_0 = emin and spin magnitudes 0 < χ_i < 0.1. Comparing with the 
results shown in Fig. <ref>, we can see that assuming lower spins prevents the 
circular waveforms from having faithfulness below the threshold line. Thus, eccentricity effects in the inspiral can be safely ignored when only the last part of it is visible.
This further shows that the starting eccentricity is the most important factor to influencing the 
accuracy of our waveforms in the late inspiral.



 
We present in Fig. <ref> the results from late inspiral run (VIa), with starting 
eccentricity e_0 = emin and spin magnitudes χ_i = 1. The results here are similar to the ones shown in Fig. <ref>, but more pronounced. The proportion of binaries for which the circular waveform has a faithfulness lying below the threshold line increases to 25 %, indicating that the inclusion of eccentricity effects might be important even for fully circularized binaries in the last stages of their inspiral when their spins are large.




 §.§ Early inspiral systems




We present in Fig. <ref> the results from early inspiral run (Ib), with starting 
eccentricity 10^-5 < e_0 < 0.5 and spin magnitudes 0 < χ_i < 1. We can see that, in this 
case, using circular waveforms will likely result in large biases even when the starting 
eccentricity is below 10^-3. The large number of orbital cycles accumulated is such that the 
small difference in the frequency evolution induces very low faithfulness even for very low 
eccentricities. On the other hand, the eccentric waveforms perform better than in the late inspiral 
case. In the top panel, we can see that the low-order EE_6,0 waveform stays above the 
faithfulness threshold for e_0 ≲ 0.05, while the high-order one EE_6,2 is above the 
threshold for the whole parameter space that we investigated. In the bottom panel, we can see that 
the waveform with circular amplitudes EE_0,2 stays above the threshold for e_0 ≲ 0.1, 
while the waveforms EE_M,2, M ≥ 2 do so for e_0 ≲ 0.3.
 



We present in Fig. <ref> the results from early inspiral run (IIb), with starting 
eccentricity 10^-5 < e_0 < 0.5 and spin magnitudes 0 < χ_i < 0.1. We can see that, for 
initial eccentricities e_0 ≳ 10^-4, circular waveforms yield a faithfulness below F = 
0.9. Thus, even if they are slowly spinning, the use of circular waveforms for parameter estimation 
for such binaries is likely to yield important biases. Using eccentric waveforms for early inspiral 
systems is therefore crucial in order to ensure accurate parameter recovery, even with initial 
eccentricities as low as e_0 ∼ 10^-4. In the bottom panel, similarly to run (Ib), we can see 
that the waveform with circular amplitudes EE_0,2 stays above the threshold for e_0 ≲ 
0.1, while the waveforms EE_M,2, M ≥ 2 do so for e_0 ≲ 0.3.



We present in Fig. <ref> the results from early inspiral run (IIIb), with starting 
eccentricity 10^-5 < e_0 < 0.5 and spin magnitudes χ_i = 1. While the results for the eccentric waveforms are similar to the ones shown in Figs. <ref> and <ref>, the circular waveform never reached a median faithfulness above the threshold line above an initial eccentricity of e_0 = 3 × 10^-5. This indicates that highly spinning systems in the early inspiral will require the use of an eccentric model irrespective of their initial eccentricity.



We present in Fig. <ref> the results from early inspiral run (IVb), with starting 
eccentricity e_0 = emin and spin magnitudes 0 < χ_i < 1. We can see that for these 
systems, including the eccentricity in the phasing is important, but the order used in other 
effects matters very little. Indeed, the faithfulness distributions for the two eccentric waveforms 
EE_0,0 and EE_6,2 are indistinguishable and have support almost exclusively above 
the faithfulness threshold, whereas the faithfulness distribution for the circular waveforms has 
46 % support below the threshold.



We present in Fig. <ref> the results from early inspiral run (Vb), with starting 
eccentricity e_0 = emin and spin magnitudes 0 < χ_i < 0.1. We can see that for these 
systems, circular waveforms have a faithfulness distribution almost identical to those of eccentric 
waveforms, indicating that when the spins are small and the binaries have fully circularized, the 
use of circular waveforms may be sufficient for unbiased parameter estimation.



We present in Fig. <ref> the results from early inspiral run (VIb), with starting 
eccentricity e_0 = emin and spin magnitudes χ_i = 1. We observe in this figure that the faithfulness distribution for the circular waveforms has 94% support below the threshold line, indicating that for highly spinning binaries, the use of eccentric waveforms will be crucial for unbiased parameter estimation. However, the distributions for the two eccentric waveforms 
EE_0,0 and EE_6,2 are indistinguishable also in this case, indicating that the precision of 
the waveform amplitude is of little importance. Thus, eccentricity and spins will be important 
to include in the analysis of stellar-origin black hole binaries with LISA to account for the 
possibility of high spins, even if the binaries have fully circularized. However, using accurate 
amplitudes might be unnecessary for those sources.



Comparing the different results, we find that the median faithfulness for each waveform is 
mainly influenced by the initial eccentricity, and the stage in the inspiral that they find 
themselves in. We summarize in Table <ref> the initial eccentricity below which the 
median faithfulness falls above the threshold line for a few of the waveforms compared in our 
simulations. Interestingly, we find that waveform EE_0,0 performs slightly better than waveform 
EE_6,0. We find the same to be true by comparing EE_0,0 to any waveform EE_M,0 or 
EE_M,1 with M>0. We thus remark that in order for the inclusion of beyond-circular effects in 
the amplitudes to increase the accuracy of the waveform, one also needs to include periastron 
precession effects at least at second order.




§ CONCLUSION


We have constructed two families of Fourier domain waveforms for spin-precessing binaries on 
eccentric orbits. These include phasing at the third nonspinning post-Newtonian order, including 
leading-order spin-orbit and spin-spin interactions, as well as instantaneous amplitudes at 
second post-Newtonian order as small eccentricity expansions. In this work, we have used amplitudes 
up to e^6, but the extension to higher orders in the eccentricity would be trivial though 
lengthy. Through comparisons with a complete time domain waveform at consistent post-Newtonian 
order, we find that our new waveforms faithfully reproduce their Fourier transform for initial 
eccentricities up to e_0 ∼ 0.3 for systems in the late inspiral, and at least up to e_0 ∼ 
0.4 for systems in the early inspiral such as stellar-origin black hole binaries as observed by 
LISA. 

Comparing results, we find that using circular waveforms would likely lead to significant biases in 
parameter recovery, even for fully circularized binaries with a signal-to-noise ratio around 25, 
provided they are highly spinning. Indeed, a 2PN spin effect prevents the eccentricity of a binary 
system from vanishing completely unless the projections of the reduced spins in the orbital plane are 
exactly equal to each other.
We find that the use of circular waveforms can cause biases if fully circularized systems 
with large spin magnitudes and random orientations are observed in the late inspiral, but not if 
the spin magnitudes are small. This situation is made worse if binary systems are observed in the 
early inspiral, and we expect large biases with circular waveforms irrespective of the initial 
eccentricity for highly spinning systems, even if they are fully circularized. However, if the 
spins are sufficiently small and the binaries have circularized below an eccentricity of 10^-4 
when the observations start, we expect the use of circular waveforms to be appropriate for 
parameter estimation. 
Overall, we expect circular waveforms to be safe to use for parameter estimation in the late inspiral if the initial eccentricity falls below 10^-2 and in the early inspiral when it falls below 10^-4, but we would recommend the use of eccentric phasing in the waveform to describe highly spinning systems, even if they have fully circularized.

Those waveforms provide a step towards the inclusion of the eccentricity in gravitational wave data 
analysis such as that performed by the LIGO/Virgo community. We argue from the simulations described in 
this paper that the inclusion of spins and eccentricity might be of importance for reducing 
potential biases in the parameter recovery of binaries, even when they are fully circularized. 
While circular templates might be appropriate to describe slowly spinning systems, it can be 
important to include in the modeling of highly spinning systems. It is worth noting that the 
faithfulness measurements described in this work are not suitable to estimate the loss of events 
due to mismodeling, or the measurability of binary parameters, including the initial eccentricity. 
We leave those questions open for future work.

Some assumptions made in this work, particularly the neglect of orbital timescale effects in the 
spin-orbit precession dynamics, have to be more closely investigated. Furthermore, the inclusion of 
the merger and ringdown signals in our waveforms is also very important work for the future, and 
will have to be taken into account in the construction of waveform templates to use in current and 
future detectors. The waveform that we have presented in this work, while useful to describe 
inspiral-dominated signals such as stellar-origin black hole binaries in LISA or neutron star 
binaries in the LIGO/Virgo network, is inspiral-only and therefore cannot be used alone in the 
characterization of merger-dominated signals such as black hole binaries as observed by the 
LIGO/Virgo network.




We thank Katerina Chatziioannou and Eliu Huerta for useful comments. We thank our referees for useful comments.
A. K. is supported by NSF CAREER Grant No. PHY-1055103, by FCT Contract No. IF/00797/2014/CP1214/CT0012 
under the IF2014 Programme, and by H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. This work was 
supported by the Centre National d'Études Spatiales. Y. B. and L. d. V. are supported by the 
Swiss National Science Foundation. This work has made use of the Horizon Cluster, hosted by the 
Institut d’Astrophysique de Paris. We thank Stéphane Rouberol for smoothly running this cluster 
for us.








§ QUADRUPOLE-MONOPOLE EFFECTS


The 2PN part of the quasi-Keplerian parametrization found in <cit.> is based upon the 
reduced Lagrangian

    ℒ/ν   = r^2/2 + 1/r + 1/r^3s_1 ·s_2 - 3/r^3(r·s_1 ) ( r·s_2 )  ,

where the reduced spins s_i = S_i / m_i. The quadrupole-monopole part of the reduced 
Lagrangian is <cit.>

    ℒQM/ν   = 1/2r^3∑_i=1^2 q_i [s_i^2 - 3 ( 
    		r·s_i )^2 ]  ,

where the quadrupole parameter q_i is defined in such a way that q_i = 1 for black holes. The 
total Lagrangian can then be written as

    ℒ/ν   = r^2/2 + 1/r + 1/2r^3s^2 - 
    		3/2r^3(r·s)^2 + 1/2r^3∑_i=1^2 ( 
    		q_i - 1 ) [s_i^2 - 3 ( r·s_i )^2 ]  ,

where s = s_1 + s_2.

Thus, a quasi-Keplerian description of the orbit including quadrupole-monopole terms can be found 
by adding the 2PN terms of <cit.>, using the substitutions 
(s_1 →s/√(2), s_2 →s/√(2)), 
(s_1 →s_1 √((q_1 - 1)/2), s_2 →s_1 √((q_1 - 1)/2)), 
and (s_1 →s_2 √((q_2 - 1)/2), s_2 →s_2 √((q_2 - 1)/2)). It reads


    r    = a ( 1 - e_r cos u) + f_r(v)  , 
    ϕ   = (1 + k) v + f_ϕ,1(v) + f_ϕ,2(v)  , 
    tanv/2   = √(1 + e_ϕ/1 - e_ϕ)tanu/2 ,
    
    	l    = u - e sin u  , 
    l̇   = n  ,


with


    a    = 1/(1 - e^2) y^2[ 1 + 1/2(1 + e^2 ) γ_1 
    		y^4 ]  , 
    
    	e_r^2    = e^2 [ 1 + ( 1 - e^2 ) γ_1 y^4 ]  , 
    
    	k    = 3/2γ_1 y^4  , 
    
    	e_ϕ^2    = e^2 [ 1 + 2( 1 - e^2 ) γ_1 y^4 ]  , 
    
    	n    = (1 - e^2 )^3/2 y ^3 (1 - 3/2γ_1 y^4 )  , 
    γ_1    = 1/2{ 3 (L·s)^2 - s^2 + 
    		∑_i=1^2 (q_i - 1) [ 3 (L·s_i )^2 - 
    		s_i^2 ] } , 
    
    	f_r(v)    = - y^2/4[ | L×s|^2 cos (2v - 2ψ) + 
    		∑_i=1^2 ( q_i - 1 ) | L×s_i |^2 cos (2v - 
    		2ψ_i) ]  , 
    
    	f_ϕ,1(v)    = - y^4 e/2[ | L×s|^2 sin (v - 
    		2ψ) + ∑_i=1^2 ( q_i - 1 ) | L×s_i |^2 sin 
    		(v - 2ψ_i) ]  , 
    
    	f_ϕ,2(v)    = - y^4/8[ | L×s|^2 sin (2v 
    		- 2ψ) + ∑_i=1^2 ( q_i - 1 ) | L×s_i |^2 
    		sin (2v - 2ψ_i) ]  ,


where ψ is the angle subtended by the total reduced spin s and the periastron line, 
ψ_i is the angle subtended by the individual reduced spin s_i and the periastron line 
(see Fig. <ref>), and the periastron line is defined by the equation v = u = l = 2 p 
π, p ∈ℤ. We can then use this representation of the orbit together with the orbit 
averaged evolution equations for the energy and orbital angular momentum computed 
in <cit.> to find


    M dy/dt =    ν(1 - e^2 )^3/2{(32/5 + 28/5 
    		e^2 ) y^9 + σ[ - 84/5 - 228/5 e^2 - 33/5 e^4, 
    		242/5 + 654/5 e^2 + 381/20 e^4, -447/10 e^2 - 93/10 
    		e^4, 
       88/5 - 16 q + (48 - 216/5 q ) e^2 + ( 69/10 - 
    		63/10 q ) e^4, - 244/5 + 48 q + (-132 + 648/5 q ) 
    		e^2 + (- 96/5 + 189/10 q ) e^4, 
       (1 - q) ( 447/10 e^2 + 93/10 e^4 ) ] y^13} , 
    
    	M de^2/dt =     - ν( 1 - e^2 )^3/2{(608/15 e^2 + 
    		242/15 e^4 ) y^8 + σ[ 2/3 - 1961/15 e^2 - 
    		2527/12 e^4 - 157/8 e^6, -2/3 + 5623/15 e^2 + 
    		2393/4 e^4 
       + 447/8 e^6, -5527/30 e^2 - 10117/30 e^4 - 5507/160 e^6,
    		- 4/3 + ( 682/5 - 1876/15 q ) e^2 + (1337/6 
    		- 595/3 q ) e^4 + (83/4 - 37/2 q ) e^6, 
       4/3 + ( - 5618/15 + 1876/5 q ) e^2 + ( 
    		-1203/2 + 595 q ) e^4 + ( - 225/4 + 111/4 q) e^6 
    		,
       ( 2764/15 - 921/5 q ) e^2 + (1687/5 - 
    		5056/15 q ) e^4 + ( 551/16 - 172/5 q ) e^6 ] 
    		y^12} ,


where

    σ(a, b, c, a_1 + a_2 q, b_1 + b_2 q, c_1 + c_2 q) =     a s^2 + b ( L·s)^2 + c | L×s|^2 cos 2 ψ
       + ∑_i=1^2 [ (a_1 + a_2 q_i ) s_i^2 + ( b_1 + b_2 q_i 
    		) ( L·s_i )^2 + (c_1 + c_2 q_i ) | 
    		L×s_i |^2 cos 2 ψ_i ]  .

We thus find the residual eccentricity found in <cit.> unchanged by quadrupole-monopole 
effects.




§ QUASI-KEPLERIAN PARAMETRIZATION


A full quasi-Keplerian parametrization of the orbit at 2PN order in harmonic coordinates is <cit.>


    r    = a ( 1 - e_r cos u) + f_r(v)  , 
    ϕ   = (1 + k) v + f_ϕ(v)  , 
    tanv/2   = √(1 + e_ϕ/1 - e_ϕ)tanu/2 ,
    
    	l    = u - e sin u + f_t(u,v)  ,
    l̇   = n  ,


with


    a =    1/( 1 - e^2 ) y^2{ 1 + [ -1 + ν/3 + ( 3 - 
    		ν/3) e^2 ] y^2 + β( 2/3 + 2 e^2, 1 + e^2) 
    		y^3 + [5 + 11/4ν + ν^2/9
       + (21/2 -73/6ν - 2/9ν^2 ) e^2 
    		+ ( 1 + 5/12ν + ν^2/9) e^4 	+ (1 -	e^2)^3/2(-5 + 2 ν) + γ_1/2(1 + e^2)] y^4  , 
    
    	e_r^2 =     e^2 { 1 + (1 - e^2 ) {(8 - 3ν) y^2 + β(4, 2) y^3 + [	32 - 467/12ν + 4 ν^2 + (- 40 + 		
    		371/12ν - 4 ν^2 ) e^2 
        + √(1 - e^2)(15 - 6 ν) + γ_1] y^4 }} ,
    
    	k =     3 y^2 +β( 4, 3) y^3 + [ 27/2 - 7 ν + ( 51/4- 
    		13/2ν) e^2 + 3/2γ_1] y^4  ,
    
    	e_ϕ^2 =     e^2 {1 + (1 - e^2 ) {(8 - 2ν) y^2 + β(4, 4 ) y^3 + [ 42 - 113/12ν + 11/12ν^2 + (- 40 
    		+ 1043/48ν - 89/48ν^2 ) e^2 
       + √(1 - e^2)(15 - 6 ν) + 2γ_1 ] y^4 }} , 
    
    	n =    (1 - e^2)^3/2 y^3 { 1 - 3 y^2 - β( 4, 3 ) y^3 + [ 
    		-9/2 + 7 ν + ( -51/4 + 13/2ν) e^2 - 3/2γ_1] y^4 } ,


where


    β(a, b)    = -[ (a μ_1 + b μ_2 )s_1 + (bμ_1 + aμ_2 ) 
    		s_2 ] ·L , 
    γ_1    = 1/2{ 3 (L·s)^2 - s^2 + 
    		∑_i=1^2 (q_i - 1) [ 3 (L·s_i )^2 - 
    		s_i^2]} .


The functions f_r, f_ϕ, f_t, and f_n are given by


    f_r(v)    = ∑_i=0^2 b_r,icos (2v - 2ψ_i)  , 
    
    	f_ϕ(v)    = ∑_k=2^3 a_ϕ,ksin(k v) + ∑_k=1^2 ∑_i=0^2 b_ϕ,k,isin(kv-2ψ_i)  , 
    
    	f_t(u,v)    = g_t(u-v) + a_t sin(v)  ,


with


    b_r,i   = - y^2/4 F_i | L×s_i |^2  , 
     
    	a_ϕ,2   = e^2 ( ν/8 - 3/8ν^2 ) y^4  , 
    
    	a_ϕ, 3   = -e^33/32ν^2 y^4  ,
    
    	b_ϕ,1,i   = - e/2 F_i | L×s_i |^2 y^4  , 
    
    	b_ϕ,2,i   = - 1/8 F_i | L×s_i |^2 y^4  ,	
    
    	g_t    = (1 - e^2 )^3/2( 15/2 - 3 ν) y^4  , 
    
    	a_t    = e (1 - e^2 )^3/2( -ν/2 - ν^2/8) y^4  ,


where we defined, for convenience, s_0 = s, ψ_0 = ψ, F_0 = 1, F_1 = q_1 - 1, 
and F_2 = q_2 - 1.




§ EVOLUTION EQUATIONS


The evolution equations of y and e are given at 3PN order by 
 <cit.>


    M d y/dt   = ( 1 - e^2	)^3/2ν y^9 ( a_0 + ∑_n=2^6 a_n y^n 
    		 )  , 
    
    	 M de^2/dt   = - ( 1 - e^2 )^3/2ν y^8 ( b_0 + ∑_n=2^6 b_n y^n 
    		)  ,


where


    a_0 =    32/5 + 28/5 e^2  , 
    
    	a_2 =     -1486/105 - 88/5ν + ( 12296/105 - 5258/45ν) e^2 + ( 3007/84 - 244/9ν) e^4  , 
    
    	a_3 =    128 π/5ϕ_y + β( 904/15 + 2224/15 e^2 + 
    	99/5 e^4, 40 + 1916/15 e^2 + 314/15e^4 )  , 
    
    	a_4 =    34103/2835 + 13661/315ν + 944/45ν^2 + ( 
    		-256723/945 - 173587/315ν + 147443/270ν^2 ) e^2 + 
    		( 2095517/7560 - 589507/504ν + 34679/45ν^2 ) e^4 
    		
       + (53881/2520 - 7357/90ν + 9392/135ν^2 ) e^6 
    		+ e^2/1 - e^2( 85/6 + 1445/6ν) + 1 - √(1 
    		- e^2)/√(1 - e^2)[ 16 - 32/5ν + (266 - 532/5ν) 
    		e^2 
       + ( -859/2 + 859/5ν) e^4 + (-65 + 26 ν) e^6 
    		] + σ[ - 84/5 - 228/5 e^2 - 33/5 e^4, 242/5 
    		+ 654/5 e^2 + 381/20 e^4, - 447/10 e^2 - 93/10 e^4, 
    		
       88/5 - 16 q + (48 - 216/5 q ) e^2 + ( 69/10 - 
    		63/10 q ) e^4, - 244/5 + 48 q + (-132 + 648/5 q ) 
    		e^2 + (- 96/5 + 189/10 q ) e^4, 
       (1 - q) ( 447/10 e^2 + 93/10 e^4 ) ]  , 
    
    	a_5 =    π( -4159/105ψ_y - 756/5νζ_y )  , 
    
    	a_6 =    16447322263/21829500 - 54784/525γ_E + 512/15π^2 + 
    		(-56198689/34020 + 902/15π^2 ) ν + 541/140ν^2 - 
    		1121/81ν^3 + [ 247611308999/87318000
       - 392048/525γ_E + 3664/15π^2 + (-2828420479/680400 + 
    		477/4π^2 ) ν + 1070903/315ν^2 - 392945/324ν^3 
    		] e^2 + [-236637777001/58212000
       - 93304/175γ_E + 872/5π^2 + ( 2963572847/453600 - 
    		53131/960π^2 ) ν + 44123941/6048ν^2 - 2198212/405ν^3 ] e^4 + [ -28913792717/6468000
       - 4922/175γ_E + 46/5π^2 + ( 107275139/30240 - 
    		369/80π^2 ) ν + 5155951/1512ν^2 - 44338/15ν^3 
    		] e^6 + (-243511057/887040 + 4179523/15120ν
       + 83701/3780ν^2 - 1876/15ν^3 ) e^8 + e^2/1 - e^2[ 
    		91284763/378000 + ( 19505077/5040 - 595/8π^2 ) ν - 
    		48569/12ν^2 - 730168/23625 (1 +	√(1- e^2))] 
    		
       + 1 - √(1 - e^2)/√(1 - e^2){-1425319/3375 + ( 
    		9874/315 - 41/30π^2 ) ν + 632/15ν^2 + [ 
    		2385427/1050 + ( -274234/45 + 4223/240π^2 ) ν + 
    		70946/45ν^2 ] e^2 
       + [8364697/4200 + ( 1900517/630 - 32267/960π^2 
    		) ν - 47443/90ν^2 ] e^4 + [-167385119/25200 + ( 
    		4272491/504 - 123/160π^2 ) ν - 43607/18ν^2 ] 
    		e^6 
       + ( -65279/168 + 510361/1260ν - 5623/45ν^2 ) e^8
    		} + 1284/175κ_y 
       + ( 54784/525 + 392048/525 e^2 + 93304/175 e^4 + 
    		4922/175 e^6 ) log[ 1 + √(1 - e^2)/8 y (1 - e^2 
    		)^3/2]  ,
    
    	b_0 =    608/15 e^2 + 242/15 e^4  , 
    
    	b_2 =    ( -1878/35 - 8168/45ν) e^2 + ( 59834/105 - 
    		7753/15ν) e^4 + ( 13929/140 - 3328/45ν) 
    		e^6  , 
    
    	b_3 =    788 π e^2/3ϕ_e + β( 19688/45 e^2 +	28256/45 
    	e^4 + 263/5 e^6, 1448/5 e^2 + 1618/3 e^4 + 167/3 e^6 ) 
    		 , 
    
    	b_4 =    ( -952397/945 + 5937/7ν + 1504/5ν^2 ) e^2 + 
    		( -3113989/1260 - 388419/140ν + 64433/20ν^2 ) e^4 
    		+ (4656611/1512 - 13057267/2520ν
       + 127411/45ν^2 ) e^6 + ( 420727/1680 - 362071/1260ν + 1642/9ν^2 ) e^8 + √(1 - e^2)[ (2672/3 - 
    		5344/15ν) e^2 + ( 2321 - 4642/5ν) e^4 
       + ( 565/3 - 226/3ν) e^6] + σ[ 2/3 
    		-1961/15 e^2 - 2527/12e^4 - 157/8 e^6, - 2/3 + 
    		5623/15 e^2 + 2393/4 e^4 + 447/8 e^6, -5527/30 e^2 - 
    		10117/30 e^4 
       - 5507/160 e^6, -4/3 + ( 682/5 - 1876/15 q ) 
    		e^2 + (1337/6 - 595/3 q ) e^4 + ( 83/4 - 
    		37/2 q ) e^6,	4/3 + ( - 5618/15 + 1876/5 q 
    		) e^2 
       + (- 1203/2 + 595 q ) e^4 + ( - 225/4 + 111/4 q 
    		) e^6 , ( 2764/15 - 921/5 q ) e^2 + (1687/5 - 
    		5056/15 q ) e^4 + ( 551/16 - 172/5 q ) e^6 ] 
    		 , 
    
    	b_5 =    π( -55691/105ψ_e - 610144/315νζ_e ) e^2  , 
    
    	b_6 =    [ 61655211971/4365900 - 2633056/1575γ_E + 24608/45π^2 + ( 43386337/56700 + 1017/5π^2 ) ν - 
    		4148897/1260ν^2 - 61001/243ν^3 ] e^2 
       + [64020009407/21829500 - 9525568/1575γ_E + 89024/45π^2 + (770214901/12600 - 15727/96π^2 ) ν - 
    		80915371/7560ν^2 - 86910509/9720ν^3 ] e^4 
       + [-1167012417073/58212000 - 4588588/1575γ_E + 42884/45π^2 + ( 8799500893/453600 - 295559/960π^2 ) ν + 
    		351962207/10080ν^2 - 2223241/90ν^3 ] e^6 
       + [120660628321/12936000 - 20437/175γ_E + 191/5π^2 + 
    		(- 91818931/5040 - 6519/320π^2 ) ν + 		
    		2495471/126ν^2 - 11792069/1215ν^3 ] e^8 
       + (302322169/887040 - 1921387/5040ν + 41179/108ν^2 - 
    		386792/1215ν^3 ) e^10 + √(1 - e^2){[ 
    		-22713049/7875 + ( -11053982/945 + 8323/90π^2 ) ν
       +108664/45ν^2 ] e^2 +[178791374/7875 + ( 
    		-38295557/630 + 94177/480π^2 ) ν + 681989/45ν^2 
    		] e^4 + [5321445613/189000
       + ( -26478311/756 + 2501/1440π^2 ) ν + 450212/45ν^2 ] e^6 + [186961/168 - 289691/252ν + 3197/9ν^2 
    		] e^8 } + 1460336/23625( 1 - √(1 - e^2)) 
       + 428/1575 e^2 κ_e + ( 2633056/1575 e^2 + 9525568/1575 
    		e^4 + 4588588/1575 e^6 + 20437/175 e^8 ) log[ 1 + √(1 
    		- e^2)/8 y (1 - e^2 )^3/2]  ,


with the tail terms given, in terms of the functions found in <cit.>, by


    ϕ_y    = (1 - e^2 )^7/2ϕ̃ ,
    ϕ_e    = 192 (1 - e^2 )^9/2/985 e^2( √(1- e^2)ϕ - 
    		ϕ̃)  , 
    ψ_y    = ( 1 - e^2 )^9/2( -8064/4159√(1 - e^2)ϕ + 
    		4032/4159ϕ̃ + 8191/4159ψ̃)  , 
    ζ_y    = (1 - e^2 )^7/2[ 160 (1 - e^2 )^3/2/567ϕ 
    		+ ( -176/567 + 80/567 e^2 ) ϕ̃ + 583 (1 - 
    		e^2 )/567ζ̃]  , 
    ψ_e    = 16382 (1 - e^2)^9/2/55691 e^2[ ( 9408/8191 - 
    		14784/8191 e^2 ) √(1 - e^2)ϕ+(-9408/8191 + 
    		4032/8191 e^2 ) ϕ̃ + (1 - e^2 ) ( √(1 - e^2)ψ - ψ̃)]  , 
    ζ_e    = 12243 (1 - e^2)^9/2/76268 e^2[-16 (1 - e^2 
    		)^3/2/53ϕ	+ ( 16/53 - 80/583 e^2 ) ϕ̃ 
    		+ (1 - e^2 ) ( √(1 - e^2)ζ - ζ̃)]  , 
    κ_y    = - 934088 (1 - e^2 )^5/33705( κ̃ - F̃)	 , 
    κ_e    = -5604528 (1 - e^2 )^6/3745 e^2[ √(1 - e^2)( 
    		κ - F ) - (κ̃ - F̃) ]  .


We chose to only include in the 3PN enhancement functions κ_i the terms proportional to 
log n, as the other ones are in finite number and can be combined with nontail terms. Using the 
formalism developed in <cit.>, we give them here at tenth order in the 
eccentricity:


    ϕ_y =     1 + 97/32 e^2 + 49/128 e^4 - 49/18432 e^6 - 
    		109/147456 e^8 - 2567/58982400 e^10 + e^12 , 
    ϕ_e =     1 + 5969/3940 e^2 + 24217/189120 e^4 + 623/4538880 e^6 - 
    		96811/363110400 e^8 - 5971/4357324800 e^10 + e^12 , 
    ψ_y =     1 - 207671/8318 e^2 - 8382869/266176 e^4 - 8437609/4791168 
    		e^6 + 10075915/306634752 e^8 - 38077159/15331737600 e^10 + e^12 , 
    ζ_y =     1 + 113002/11907 e^2 + 6035543/762048 e^4 + 253177/571536 
    		e^6 - 850489/877879296 e^8 - 1888651/10973491200 e^10 + e^12 ,
    ψ_e =     1 - 9904271/891056 e^2 - 101704075/10692672 e^4 - 
    		217413779/513248256 e^6 + 35703577/6843310080 e^8 - 
    		3311197679/9854366515200 e^10 + e^12 , 
    ζ_e =     1 + 11228233/2440576 e^2 + 37095275/14643456 e^4 + 
    		151238443/1405771776 e^6 - 118111/611205120 e^8 - 		
    		407523451/26990818099300 e^10 + e^12 , 
    κ_y =     244 log 2 ( e^2 - 18881/1098 e^4 + 6159821/39528 e^6 - 
    		16811095/19764 e^8 + 446132351/123525 e^10) - 243 log 3 (e^2 
    		- 39/4 e^4 + 2735/64e^6 
       + 25959/512 e^8 - 638032239/409600 e^10) - 48828125 log 
    		5/5184( e^6 - 83/8 e^8 + 12637/256 e^10) - 
    		4747561509943 log 7/33177600 e^10 + e^12 , 
    κ_e =     6536 log 2 (1 - 22314/817 e^2 + 7170067/19608e^4 - 
    		10943033/4128 e^6 + 230370959/15480 e^8 - 866124466133/8823600 
    		e^10) 
       - 6561 log 3 ( 1 - 49/4 e^2 + 4369/64 e^4 + 214449/512 e^6 - 
    		623830739/81920 e^8 + 76513915569/1638400 e^10) 
       - 48828125 log 5/64( e^4 - 293/24 e^6 + 159007/2304 e^8 - 
    		6631171/27648 e^10) - 4747561509943 log 7/245760( e^8 - 
    		259/20 e^10) + e^12 .



It can be noted that those enhancement functions converge much more quickly than the ones presented 
in <cit.>. Indeed, because of the inclusion of factors of √(1 - e^2) in 
them, the enhancement functions seem to converge in the parabolic limit e → 1. We believe it to 
be related to the fact that the PN parameter y we used here is related to the Newtonian orbital 
angular momentum and thus is finite and nonzero in this limit. In contrast, the PN parameter 
(Mω)^1/3 is related to the energy and thus vanishes in this limit. In that case, in order 
for the tail effects to stay nonzero, the enhancement functions are forced to diverge.




§ TRUE AND ECCENTRIC ANOMALY EXPANSION


The Fourier coefficients of the eccentric anomaly, true anomaly and orbital phase are given to 
order y^4, e^5 by


    A_1 =     e-e^3/8+e^5/192+y^4 (e^3 (105/8-51η/64 
    		-19 η ^2/64)+e (-15/2+9 η/8 + 		
    		η^2/8)+ e^5 (-735/128-489 η/512+111 
    		η^2/512))  , 
    
    	A_2 =    e^2/2-e^4/6+y^4 (e^4 (75/4-5 η/32 
    		-47 η ^2/96)+e^2 (-75/8+15 η/16+3 
    		η^2/16))  , 
    
    	A_3 =    3 e^3/8-27 e^5/128+y^4 (e^5 (6825/256+705 η/1024 -775 η ^2/1024)+e^3 (-95/8+49 η/64 
    		+17 η ^2/64))  , 
    
    	A_4 =    e^4/3+e^4 y^4 (-975/64+35 η/64+71 η^2/192)  , 
    
    	A_5 =    125 e^5/384+e^5 y^4 (-5049/256+1167 η/5120+523 
    		η ^2/1024)  , 
    
    	B_1 =     2 e-e^3/4+5 e^5/96+y^2 (e (4-η )+e^5 ( 17/48 - 
    		17 η/192)+e^3 (-7/2+7 η/8)) 		
    		
       +y^4 (e^3 (-39/4+197 η/32-7 η ^2/32)+e 
    		(13-31 η/12+η ^2/12)+e^5 (-401/192-569 
    		η/144+107 η ^2/576))  , 
    
    	B_2 =    5 e^2/4-11 e^4/24+y^2 (e^2 (4-η )+e^4 	( -14/3 + 
    		7 η/6)) 
       + y^4 (e^4 (-217/24+3187 η/288-313 η^2/288) + e^2 (91/8-101 η/24+11 η ^2/24)) 
    		 , 
    
    	B_3 =    13 e^3/12-43 e^5/64+y^2 (e^3 (9/2-9 η/8) + e^5 (-207/32+207 η/128)) 
       +y^4 (e^5 (-1011/128+9471 η/512-1161 η^2/512) + e^3 (43/4-51 η/8+7 η ^2/8))  , 
    		
    
    	B_4 =    103 e^4/96+e^4 y^2 (16/3-4 η/3)+e^4 y^4 ( 
    		1969/192-2711 η/288+205 η ^2/144)  , 
    
    	B_5 =    1097 e^5/960+e^5 y^2 (625/96-625 η/384)+e^5 y^4 
    		(3641/384-315821 η/23040+10039 η ^2/4608)  ,
    
    	C_1 =     2 e-e^3/4+5 e^5/96+y^2 (e (10-η)+e^5 ( 49/96 - 
    		17 η/192)+e^3 (-17/4+7 η/8)) 
    		
       + y^4 (e^5 (-317/96+2497 η/2304-1579 η^2/2304) + e (52-235 η/12+η^2/12)+e^3 ( 13/8 
    		- 231 η/32+17 η^2/32)) 
    
    	C_2 =    5 e^2/4-11 e^4/24+y^2 (e^2 (31/4-η)+e^4 
    		(-145/24+7 η/6)) 
       +y^4 (e^2 (323/8-163 η/12+η ^2/12)+e^4 
    		(-331/24+η/18+25 η ^2/36)) 
    
    	C_3 =    13 e^3/12-43 e^5/64+y^2 (e^3 (31/4-9 η/8)+e^5 (-543/64+207 η/128)) 
       +y^4 (e^3 (313/8-1205 η/96+η ^2/32)+e^5 
    		(-749/32+7217 η/1536+567 η ^2/512)) 
    		
    
    	C_4 =    103 e^4/96+e^4 y^2 (821/96-4 η/3)+e^4 y^4 
    		(1975/48-118 η/9-11 η ^2/144) 
    
    	C_5 =    1097 e^5/960+e^5 y^2 (9541/960-625 η/384)+e^5 
    		y^4 (10813/240-338987 η/23040-1211 η ^2/4608)  .






§ WAVEFORM AMPLITUDES EXPANSION


The amplitudes G_+,×^(m,n) in Eq. (<ref>) are given to order 
y^2, e for n < 0, with C_i = cos i = L·N and S_i =sin i, by


    G_+^(1,-1)   = e (1/2-C_i^2/2)+e y^2 
    		(-35/16+9 C_i^2/4-C_i^4/16+η(-5/48-C_i^2/12+3 C_i^4/16))  ,
    		
     
    	G_+^(1,-2)   = e y δ(-3 S_i/2+1/2 C_i^2 
    		S_i)  ,
     
    	G_+^(1,-3)   = e (-9/4-9 C_i^2/4)+e y^2 
    		(-9/8-63 C_i^2/8+η(-75/8+87 
    		C_i^2/8))  ,
     
    	G_+^(1,-4)   = e y δ(4 S_i+4 C_i^2 S_i)  ,
    		
     
    	G_+^(1,-5)   = e y^2 (-625 S_i^2/96-625/96 C_i^2 
    		S_i^2+η(625 S_i^2/32+625/32 C_i^2 
    		S_i^2))  ,
     
    	G_+^(0,-1)   = y δ(-5 S_i/8-1/8 C_i^2 
    		S_i)  ,
     
    	G_+^(0,-2)   = -1-C_i^2+y^2 (19/6+3 
    		C_i^2/2-C_i^4/3+η(-19/6+11 
    		C_i^2/6+C_i^4))  ,
     
    	G_+^(0,-3)   = y δ(9 S_i/8+9/8 C_i^2 
    		S_i)  ,
     
    	G_+^(0,-4)   = y^2 (-4 S_i^2/3-4/3 C_i^2 
    		S_i^2+η(4 S_i^2+4 C_i^2 S_i^2)) 
    		 ,
     
    	G_+^(-1,-1)   = e (3/4+3 C_i^2/4)+e y^2 
    		(143/24+53 C_i^2/8+C_i^4/12+η(7/8-9 C_i^2/8-C_i^4/4)) 
    		 ,
     
    	G_+^(-1,-2)   = e y δ(-3 S_i/2-3/2 C_i^2 
    		S_i)  ,
     
    	G_+^(-1,-3)   = e y^2 (81 S_i^2/32+81/32 C_i^2 
    		S_i^2+η(-243 S_i^2/32-243/32 C_i^2 
    		S_i^2))  ,
     
    	G_×^(1,-1)   = i e y^2 (-1/4 C_i S_i^2+3/4η C_i S_i^2)  ,
     
    	G_×^(1,-2)   = i e y δ C_i S_i ,
     
    	G_×^(1,-3)   = i (9 e C_i/2+e y^2 (9 
    		C_i-27/8 C_i S_i^2+η(-3 
    		C_i/2+81/8 C_i S_i^2))) 
    		 ,
     
    	G_×^(1,-4)   = -8 i e y δ C_i S_i ,
     
    	G_×^(1,-5)   = i e y^2 (625/48 C_i S_i^2-625/16η C_i S_i^2)  ,
     
    	G_×^(0,-1)   = 3/4 i y δ C_i S_i ,
     
    	G_×^(0,-2)   = i (2 C_i+y^2 (-13 C_i/3-4/3 
    		C_i S_i^2+η(C_i/3+4 C_i 
    		S_i^2)))  ,
     
    	G_×^(0,-3)   = -9/4 i y δ C_i S_i ,
     
    	G_×^(0,-4)   = i y^2 (8/3 C_i S_i^2-8 η C_i 
    		S_i^2)  ,
     
    	G_×^(-1,-1)   = i (-3 e C_i/2+e y^2 (-38 
    		C_i/3+11/24 C_i S_i^2+η(C_i/2-11/8 C_i S_i^2))) 
    		 ,
     
    	G_×^(-1,-2)   = 3 i e y δ C_i S_i ,
     
    	G_×^(-1,-3)   = i e y^2 (-81/16 C_i S_i^2+243/16η C_i S_i^2)  .


Note that G_+,×^(m,-n) = G̅_+,×^(-m,n).







