A 3PN Fourier Domain Waveform for Non-Spinning Binaries with Moderate Eccentricity

TL;DR

This paper develops a 3PN order Fourier domain gravitational waveform model for non-spinning eccentric binaries, improving accuracy and including effects like periastron precession, to aid future gravitational wave data analysis.

Contribution

It extends previous eccentric waveform models to 3PN order, incorporating higher-order effects and validating the model against numerical simulations.

Findings

Achieves 97-99% match with numerical models in key parameter regions.

Identifies where eccentric effects significantly impact gravitational wave signals.

Provides a validated, accurate model suitable for data analysis of eccentric binaries.

Abstract

While current gravitational wave observations with ground based detectors have been consistent with compact binaries in quasi-circular orbits, eccentric binaries may be detectable by ground-based and space-based instruments in the near future. Eccentricity significantly complicates the gravitational wave signal, and we currently lack fast and accurate models that are valid in the moderate to large eccentricity range. In a previous paper, we built a Fourier domain, eccentric waveform model at leading order in the post-Newtonian approximation, i.e. as an expansion in small velocities and weak fields. Here we extend this model to 3rd post-Newtonian order, incorporating the effects of periastron precession and higher post-Newtonian order effects that qualitatively change the waveform behavior. Our 3PN model combines the stationary phase approximation, a truncated sum of harmonics of combinations of two orbital frequencies (an azimuthal one and a radial one), and a bivariate expansion in the orbital separation and the eccentricity. We validate the model through comparisons with a fully-numerical, time domain post-Newtonian model, and find good agreement (matches between $97\% - 99\%$) in much of the parameter space. We estimate in what regions of parameter space eccentric effects are important by exploring the signal to noise ratio of eccentric corrections. We also examine the effects of higher post-Newtonian order terms in the waveform amplitude, and the agreement between different PN-consistent, numerical, time-domain models. In an effort to guide future improvements to the model, we gauge the error in our 3PN model incurred by the different analytic approximations used to construct it. This model is useful for preliminary data analysis investigations and it could allow for a phenomenological hybrid that incorporates eccentricity into an inspiral-merger-ringdown model.

Tables

Table 1. Table 1: Eccentricity order to which each of the coefficients in Eq. ( 65 ) is expanded past leading order for the TaylorF2e+ model and the TaylorF2e- model in parenthesis. We keep the same number of terms in l n ′ ​ ( e ) superscript subscript 𝑙 𝑛 ′ 𝑒 l_{n}^{\prime}(e) as in λ n ′ ​ ( e ) superscript subscript 𝜆 𝑛 ′ 𝑒 \lambda_{n}^{\prime}(e)

$n$	${\lambda }_n^{\prime }\left(e\right)$	$t_n^{\prime }\left(e\right)$
0	10 (4)	36 (12)
2	240/19 (8)	620/19 (202/19)
3	208/19 (4)	512/19 (170/19)
4	278/19 (126/19)	26 (164/19)
5	246/19 (132/19)	468/19 (182/19)
6	14 (94/19)	328/19 (100/19)

Equations

$${\tilde{h}}_{+,\times }\sim \sum _{j=1}^N{A}_{+,\times }^{(j)}(f)e^{i{\psi }_j(f)},$$

$$r$$

$$\varphi$$

$$\lambda$$

$$l$$

$$v$$

$$l={\int }^{t}n(t^{\prime })d{t}^{\prime }.$$

$$\frac{dn}{dt}$$

$$\frac{de}{dt}$$

$$r$$

$$\varphi$$

$$W$$

$$+\left(g_{4\varphi }+g_{6\varphi }\right)\sin 3\nu +i_{6\varphi }\sin 4\nu +h_{6\varphi }\sin 5\nu ,$$

$$\lambda$$

$$l$$

$$=u-e\sin u+\left(g_{4t}+g_{6t}\right)(\nu -u)$$

$$+\left(f_{4t}+f_{6t}\right)\sin \nu +i_{6t}\sin 2\nu +h_{6t}\sin 3\nu ,$$

$$v$$

$$y=\frac{(m(1+k)n{)}^{1/3}}{\sqrt{1-e^2}}=\mathcal{O}(v).$$

$$k$$

$$\displaystyle+\bigg{\{}6720-(20000-492\pi^{2})\eta+896\eta^{2}$$

$$+[18336-(22848-123{\pi }^2)\eta +5120{\eta }^2]e^2$$

$$+(2496-1760\eta +1040{\eta }^2)e^4$$

$$\displaystyle+\left[1920-768\eta+(3840-1536\eta)e^{2}\right]\sqrt{1-e^{2}}\bigg{\}}\frac{y^{6}}{128}.$$

$$e_t^{\mathrm{ADM}}$$

$$+[\frac{1}{2}+\left(\frac{16739}{1680}-\frac{21}{16}{\pi }^2\right)\eta -\frac{83}{24}{\eta }^2$$

$$+\left(\frac{1}{2}+\frac{249}{16}\eta -\frac{241}{24}{\eta }^2\right)e_t^2]y^6\}.$$

$$\frac{dl}{dt}$$

$$-\frac{1}{128}{y}^6[-192-(14624-492{\pi }^2)\eta +896{\eta }^2+e^2\left(5120{\eta }^2-\left(17856-123{\pi }^2\right)\eta +8544\right)$$

$$+e^4\left(1040{\eta }^2-1760\eta +2496\right)+\sqrt{1-e^2}\left(e^2(3840-1536\eta )-768\eta +1920\right)]\}$$

$$\frac{d\lambda }{dt}$$

$$m\frac{dy}{dt}$$

$$m\frac{de}{dt}$$

$$h_{+,\times }$$

$$+y^2{H}_{+,\times }^{(1)}+y^3{H}_{+,\times }^{(3/2)}+y^4{H}_{+,\times }^{(2)})$$

$$H_{+}^{(0)}$$

$$+\left(P_{C2S2}^0\cos 2W+P_{S2S2}^0\sin 2W\right)\sin 2\bar{\lambda }+P^0,$$

$$H_{\times }^{(0)}$$

$$+\left(X_{C2S2}^0\cos 2W+X_{S2S2}^0\sin 2W\right)\sin 2\bar{\lambda },$$

$$P_{C2C2}^0$$

$$P_{C2S2}^0$$

$$P^0$$

$$X_{C2C2}^0$$

$$X_{C2S2}^0$$

$$C$$

$$\mathcal{O}(h_1,h_2)=\frac{(h_1|h_2)}{\sqrt{(h_1|h_1)(h_2|h_2)}},$$

Full text

I Introduction

With the multiple detections of gravitational waves (GWs) from binary black hole and neutron star systems by ground based detectors Abbott et al. (2017a, b, c, d, 2016a, 2016b) and the imminent launch of space based detectors lis , we are now in a position to perform gravitational wave astronomy. With this comes the ability to answer important questions about the physical characteristics of the sources producing those waves, the environments in which they are formed, the validity of General Relativity, and the nature of the universe at large. To answer these questions, we rely heavily on our ability to model these sources – modeling which initially was focused on binaries in quasi-circular orbits. In more recent years, there has been much work to model the GWs incorporating the effects of orbital eccentricity.

This work has been motivated by astrophysical studies which suggest that, while quasi-circular binaries dominate potential ground based detector sources, there are a number of scenarios that can produce binaries in eccentric orbits while emitting GWs in the sensitivity band of these detectors. Recently, Rodriguez et al. incorporated PN effects into globular cluster simulations and as a result predicted a merger rate for LIGO/Virgo of 0.5 Gpc*-3* yr*-1* for binaries with detectable eccentricity above $0.1$ (with an overall merger rate between 10 and 213 Gpc*-3* yr*-1*) Rodriguez et al. (2018). Other promising eccentric formation scenarios include field triples Antonini et al. (2017), GW capture and triples near supermassive BHs O’Leary et al. (2009); Antonini et al. (2014). As a result, the detection of GWs from multiple eccentric sources could help us identify the environment in which they form Nishizawa et al. (2017).

Neglecting the effects of eccentricity in source modeling can lead to a loss in detection rate and parameter bias. Several studies have investigated the ability of quasi-circular template banks to recover eccentric signals and found significant loss in detection rate for binaries with eccentricities greater than about 0.1 Huerta and Brown (2013); Brown and Zimmerman (2010); Cokelaer and Pathak (2009). However, even for eccentricities smaller than 0.1, parameter biases can arise when the eccentricity is neglected in the model and the signal comes from a binary with eccentricity as small as about $10^{-3}$ Favata (2013).

Clearly then, if we wish to recover and perform data analysis on eccentric signals we require models which incorporate the effects of orbital eccentricity, particularly efficient Fourier domain waveforms. The Post-Circular (PC) formalism provides a fully analytic (at Newtonian order) prescription to compute a frequency-domain waveform, while incorporating eccentric effects Yunes et al. (2009). This was extended to 2PN order in Tanay et al. (2016) and recently to 3PN order Tiwari and Gopakumar . In Moore et al. (2016), the authors used a circular GW amplitude and incorporated the leading-order eccentric effects in the binary dynamics while working at 3PN order. These models expand all relevant quantities in small eccentricity and relate the Fourier frequency to the orbital frequency in order to map the model to the Fourier domain. As a result, these models are computationally efficient, but likely only valid in the small eccentricity regime ($e_{0}\leq 0.3$).

In Reference Moore et al. (2018) (hereafter referred to as paper 1), we built a Fourier domain, eccentric waveform model at leading order in the post-Newtonian (PN) approximation 111The PN approximation is one in which the field equations are solved assuming small velocities and weak gravitational fields in an expansion in powers of $(\frac{v}{c})$, where $v$ is the orbital velocity and $c$ is the speed of light Blanchet (2006). By $n$PN order we mean an expansion to order $({v}/{c})^{2n}$.. This Newtonian order model schematically took the form

[TABLE]

where $A^{(j)}_{+,\times}$ are slowly-varying amplitudes and the $\psi_{j}$ are rapidly-varying Fourier phases. The multiple harmonics in the Fourier domain arise from the fact that an eccentric signal can be decomposed as a sum over harmonics of the mean orbital frequency in the time domain. We validated this waveform for large eccentricities ($e_{0}\leq 0.9$) using the match (a measure of agreement between two models) for advanced LIGO at design sensitivity against a fully-numerical, time-domain model. The key to our Newtonian Fourier-domain model’s accuracy and distinction from PC models is an analytic prescription in terms of the orbital eccentricity $e$ for the Fourier phases and the amplitudes, with $e$ then related to the fourier frequency though the stationary phase condition. A similar model to 1PN order is presented in Mikóczi et al. (2015).

We here extend the Newtonian model of paper 1 to 3PN order (developing a model similar to TaylorF2 Buonanno et al. (2009) but for eccentric binaries), and we validate it in a large region of initial-separation/initial-eccentricity parameter space. We follow the same procedure as in paper 1, with the most major complication arising from the inclusion of periastron precession due to PN effects. This gives rise to harmonics that are combinations of the two orbital frequencies in the eccentric problem – one related to the radial period and the other related to the azimuthal period. In addition, the Fourier phases do not admit exact solutions in terms of $e$ as was the case in paper 1, where those phases could be expressed exactly in terms of hypergeometric functions. However, in paper 1 we found that the small eccentricity expansion of those hypergeometric functions were still highly accurate, and so we here employ that same expansion.

Figure 1 shows the match between our 3PN Fourier domain model and a fully numerical time-domain PN waveform for a $(1.4,1.4)M_{\odot}$ binary inspiral for a number of initial separations and initial eccentricities. We see that our model is very faithful to the numerical one in a large region of the parameter space explored. The Fourier and time-domain models begin to disagree for very eccentric binaries that start at very small initial separation, since then their pericenter velocity becomes quite large. The results presented in this figure are representative of systems with other comparable masses. The match, however, deteriorates more rapidly for systems with disparate masses, suggesting that higher PN order models may be needed for these cases.

In order to develop our model and to better understand the features of the model, we produce several byproducts. First, we obtain a 3PN accurate decomposition of the time-domain waveform into an infinite sum of harmonics and we truncate it at a finite order. With this at hand, we then use numerical solutions to the orbital dynamics to validate the truncated model against another model that keeps all harmonics, so as to determine the loss of faithfulness due to harmonic truncation. Second, we investigate how different PN-consistent ways of solving the orbital dynamics lead to a loss of faithfulness in numerically-solved, time domain models due only to the truncation of the PN series. Third, we produce two Fourier domain waveforms, one with more terms kept in the eccentricity expansions and one with less terms kept, and perform a faithfulness analysis on each in order to understand the error incurred by our expansions in small eccentricity. Lastly, we swap the analytic components of our Fourier domain waveform with numerical solutions in order to get a sense of which of these analytic approximations leads to a dominant loss in faithfulness. Since each of these results isolate a particular approximation made in generating our Fourier domain model, we are able to pinpoint which approximation causes the most error, which can guide future advances aimed at increasing the accuracy of our model.

The organization of the paper is as follows. In Section II we review the PN dynamics of an eccentric orbit, and in Section III the time domain GW signal and its decomposition into harmonics of the two orbital frequencies. Section IV presents the GW signal in the Fourier domain, and provides the formalism that we employ to obtain our model. Section V provides details on the implementation of the Fourier domain model, and carries out a faithfulness study. In Section VI we conclude and point to future work. Throughout this work we use geometric units (c = G = 1).

II Review of Eccentric Orbital Dynamics

In this section we review the orbital dynamics of an eccentric binary. We first begin by reviewing the dynamics at Newtonian order. Then we present the dynamics including effects up to 3PN order in the quasi-Keplerian (QK) formalism presented in Damour et al. (2004). The most notable changes that occur when including higher PN effects is periastron precession (entering at 1PN order), which introduces a second frequency (and thus a second phase angle) into the dynamics, and higher order corrections to Kepler’s equation.

II.1 Newtonian Order

At Newtonian order the dynamics of an elliptical orbit restricted to a plane are described by

[TABLE]

Here $r$ is the magnitude of the relative separation vector, which we choose to be on the $x-y$ plane $\vec{r}=(r\cos\phi,r\sin\phi,0)$. The angle $\phi$ is the orbital phase, $e$ the orbital eccentricity, and the angle $l$ is the mean anomaly. The mean motion $n$ is related to the orbital period $P$ by $n=2\pi/P$, and the semi-major axis via Kepler’s third law $n^{2}a^{3}=m$ with $m$ the total mass of the system. We have written $\phi$ in terms of a linearly increasing phase angle $\lambda$, and an oscillatory phase angle $W$, which accounts for the fact that the rate of change of the orbital phase $\dot{\phi}$ is not constant throughout one orbit. In the absence of radiation reaction, the conservation of energy and angular momentum demands that the orbital eccentricity and semi-major axis remain fixed. To obtain $r(t)$ and $\phi(t)$, a separation and an orbital eccentricity must be specified and then these are used to solve Kepler’s equation [Eq. (2d)] for $u(t)$ – which can only be exactly solved numerically. This points to a clear complication of eccentric binary modeling – even in the Newtonian treatment and in the absence of radiation reaction, the explicit time domain dynamics require a numerical solution.

In GR, GWs carry away energy and angular momentum from the binary, and in response the orbital eccentricity and mean motion (or alternatively the separations associated with the ellipse) vary in time. Equations (2) remain the same with the exception of $l$, which obeys

[TABLE]

The time evolution of the eccentricity and mean motion were first derived by Peters and Matthew Peters and Mathews (1963) and is given by

[TABLE]

The appearance of the negative sign in $\dot{e}$ implies the well known result that radiation reaction causes a binary to circularize. In paper 1, we solved this coupled set of differential equations exactly for $n(e)$, which in turn allowed us to obtain $t(e)$ and $l(e)$ exactly in terms of hypergeometric functions and the AppelF1 function. We now turn to treat the effects of higher order PN effects in the orbital dynamics.

II.2 Inclusion of PN Effects

The quasi-Keplerian formalism (QK) of Damour et al. (2004) allows one to describe elliptic motion while including higher order PN effects. Equations (2) keep a similar form, with the notable difference being that the orbital angle $\lambda$ is related to an azimuthal frequency and $l$ is related to a radial one, a distinction caused by periastron precession. Including PN effects causes Eqs. (2) to take the form

[TABLE]

Note that now there are three different eccentricities: the time eccentricity $e_{t}$, the radial eccentricity $e_{r}$, and the phi-eccentricity $e_{\phi}$. The three eccentricities can be related to one another through PN accurate relations, which can be found in Memmesheimer et al. (2004); Damour et al. (2004). Throughout this work we parameterize the problem in terms of the time eccentricity, $e_{t}$, and therefore we will drop the $t$ subscript from hereon out. The functions $g_{n\phi},g_{nt},f_{n\phi},f_{nt},i_{n\phi},i_{nt},h_{n\phi},h_{nt}$ are $n/2$ PN coefficients that depend on the orbital energy, angular momentum and the masses of the binaries Memmesheimer et al. (2004).

We adopt a PN parameter $y$ defined by

[TABLE]

Note that in the zero eccentricity limit, $y$ corresponds exactly to the orbital speed $v$. At Newtonian order, the parameter $y$ is related to the dimensionless semilatus rectum, $p$, by $y=1/\sqrt{p}$. This PN parameter is well-behaved even for highly eccentric orbits, since one can think of it as effectively requiring the semi-latus rectum to be large instead of requiring the semi-major axis to be large, which could result in large pericenter velocities.

The effect of periastron precession is encoded in the parameter $k=\delta\Phi/(2\pi)$, where $\delta\Phi$ is the advance of periastron over a period. In terms of $y$ at 3PN order

[TABLE]

We can use this expression for $k$ to relate the radial orbital angular frequency, $n$, which corresponds to the periastron to periastron period to the azimuthal orbital angular frequency, $\omega=(1+k)n$, which is related to the period to return to the same angle in the orbit.

The functional dependence of the time eccentricity on the PN expansion parameter $y$ depends on the choice of gauge. Two commonly employed gauges are the Arnowitt-Deser-Misner (ADM) and modified harmonic (MH) ones, which differ at 2PN order; throughout this work, we use the MH gauge. In order to express results in either gauge one needs to relate the time eccentricity in one gauge to another. These relations can be found in Moore et al. (2016), but we also list them here:

[TABLE]

Much like in the Newtonian case, the effect of the emission of GWs is to cause the orbital frequency (or alternatively $y$) and the eccentricity to vary with time. To fully solve for an eccentric orbit including PN effects, one has to solve the following set of coupled ordinary differential equations:

[TABLE]

where the $a_{n}$ and $b_{n}$ can be found in modified harmonic gauge in Appendix C of Klein et al. (2018) and $\eta=m_{1}m_{2}/m^{2}$ is the symmetric mass ratio. The equations for $\dot{y}$ and $\dot{e}$ depend on certain enhancement functions $\phi_{y}$, $\phi_{e}$, $\psi_{y}$, $\psi_{e}$, $\zeta_{y}$, $\zeta_{e}$, $\kappa_{y}$, and $\kappa_{e}$  Klein et al. (2018), which we have expanded to $\mathcal{O}(e^{50})$. Once the solution to this set of ordinary differential equations is found, one can insert them into Eqs. (5) to obtain the evolution of the different orbital quantities (while Kepler’s equation must still be solved numerically).

III Eccentric GWs in the Time Domain

In this section we present the GW signal in the time domain. We begin by analyzing the effect of PN corrections to the amplitude of the waveform. This analysis suggests that we can faithfully model signals while restricting to leading PN order in the amplitude. We then obtain and analyze the decomposition of this signal (keeping only leading order PN corrections in the amplitude) into harmonics of the azimuthal and radial orbital frequencies. We validate this decomposition by computing the overlap (a measure of the agreement between two models) between the waveform decomposed into a number of harmonics and one which has not been decomposed into harmonics. Lastly, we discuss different PN consistent ways of solving Eq. (9) and the implications this has for modeling eccentric systems. Throughout we provide details on how we obtain the fully numerical model with which we will compare and validate our frequency domain model.

III.1 PN Effects in the Amplitude

To 2PN order in the amplitude, the GW polarizations read:

[TABLE]

The lengthy expressions for $H^{(n)}_{+,\times}$ can be found in Gopakumar and Iyer (2002) in ADM gauge and using a PN expansion parameter related to $mn$. The corrections needed to express them in MH and in terms of $y$ can be found in Appendix A. For reference we list the 0PN polarizations here (which are unaffected by gauge differences or the change of PN parameter between $mn$ and $y$):

[TABLE]

where we have defined

[TABLE]

and where the different orbital parameters entering the polarizations are described in Sec. II, while

[TABLE]

with angles $\beta$ and $i$ polar angles that describe the polarization axes.

Before we can attempt to obtain a Fourier domain model for the GW signal, we must first decide on how many PN corrections to keep in the amplitude of the polarization in Eq. (III.1). In order to do so, we will make use of the overlap, i.e. the normalized inner product in Fourier space, weighted by the spectral noise of the detector. Between two signals $h_{1}$ and $h_{2}$, the overlap is defined by

[TABLE]

where the inner product is defined as

[TABLE]

and $S_{n}(f)$ is the spectral noise density of the detector. An overlap of unity between two models indicates perfect agreement, while an overlap of -1 indicates models perfectly out of phase, with a vanishing overlap indicating perfect disagreement.

In order to gauge the importance of inclusion of higher order PN corrections to the amplitude we calculate the overlap between two models. One model consists of numerically solving for the orbital parameters (i.e. $u,~{}W,~{}l,~{}\lambda,~{}e,~{}\text{and }~{}y$) using their evolution equations at 3PN order, and keeping all corrections up to 2PN order in the amplitude. The other consists of again numerically solving for the orbital parameters using their evolution equations at 3PN order, but varying the PN order to which the amplitude of the polarizations in Eq. (III.1) is kept. We also explore the effect of solving for $u(t)$ using the Newtonian version of Kepler’s equation.

Before presenting the results of the overlap, let us provide some technical details of how these models are constructed. We begin by numerically solving (using the Burlisch-Stoer method implemented in C++) the coupled set of differential equations $({\dot{\lambda},\dot{l},\dot{e},\dot{y}})$ provided by Eqs (9). We stop the evolution when $y>(3(1+e))^{-1}$, a choice that is explained in detail in Section IV.2, but we see that in the circular limit this corresponds to an orbital velocity of a third the speed of light. With the time evolution of these parameters in hand, we must then solve the 3PN Kepler’s equation. This is done by adapting Mikkola’s method Mikkola (1987), a fast and accurate numerical scheme to solve Kepler’s equation at Newtonian (0PN) order. To solve Kepler’s equation at higher PN order [Eq. (5e)], we first use Mikkola’s method on the Newtonian equation, and then use the result as an initial guess to solve Eq. (5e) using Mikkola’s method again (see Tessmer and Gopakumar (2007); Tanay et al. (2016) for more details). To obtain $W$, we plug in the solution for $u$, $y$, and $e$ into Eq. (5f) to obtain the true anomaly $v$, and then, all of these are substituted into Eq. (5) to obtain $W$.

The numerical solutions for the orbital elements are plugged into the GW polarization to form $h(t)$ (where we choose to sample all time domain quantities at 8192 Hz). This time series is then zero padded on either side such that its length is the nearest integer power of 2, and we take the discrete Fourier transform using the FFTW library implemented in C++. The overlaps are computed from $f_{1}=1$Hz to $f_{2}=4096$Hz. We use an interpolation of a set of values of the noise curve corresponding to stationary Gaussian noise in aLIGO at design sensitivity Abbott et al. (2017e). The integrals appearing in the overlap are evaluated as sums with frequency resolution determined by the frequency resolution of the discrete Fourier transform, which depends on the length of signal, and as a result varies between different initial conditions.

Figure 2 presents density contours for the value of the overlap between various templates as a function of the initial eccentricity and dimensionless semi-latus rectum. In the equal mass case, the overlaps are not affected past the third decimal place if we keep the amplitude to Newtonian order. We also find that keeping Kepler’s equation at Newtonian order when solving for $u(t)$ only affects the overlaps in the fifth decimal place. This suggests we can keep the amplitude and Kepler’s equation to Newtonian order, which considerably simplifies the decomposition of $h(t)$ into a sum over harmonics of $l$ and $\lambda$, which in turn is crucial to analytically approximate the Fourier transform. However, half order contributions vanish in the equal mass limit and thus the overlap is affected in the 2nd digit for unequal mass binaries. This amount of loss in overlap may be considerable, depending on the data analysis application one has in mind, and it is not due to solving Kepler’s equation at Newtonian order, but rather due to neglecting higher order PN effects in the amplitude. We will return to this loss of accuracy for non-equal mass systems later.

III.2 Harmonic Decomposition

Following paper 1, a crucial step in obtaining a Fourier domain waveform is to decompose the time domain signal into a sum over harmonics of the angles related to the mean orbital frequencies $n$ and $\omega$. We carry out such a decomposition here and then verify that it faithfully represents the un-decomposed signal using overlaps. The results of Section III.1, however, indicate that we may retain only the leading order (Newtonian) contributions in the GW amplitude and in Kepler’s equation, at least for comparable mass binaries, and thus, we will restrict attention to such a model.

The results of this section rely heavily on the work to analytically solve 3PN accurate Kepler’s equation by Boetzel et. al. in Boetzel et al. (2017), which provides a useful decomposition of harmonics of $W$ into an infinite sum of harmonics of $l$:

[TABLE]

Re-expressing the polarizations in Eq. (11) in terms of complex exponentials and inserting the above decomposition yields

[TABLE]

with

[TABLE]

where c.c. stands for the complex conjugate, and where $F_{+}$ and $F_{\times}$ are the antenna functions of the detector

[TABLE]

with $\theta$ and $\Phi$ the sky angles associated with the orientation of the detector relative to the source, and $\psi$ the polarization angle that defines the rotation from the wave’s polarization basis to that defined by the arms of the detector.

Inspection of Eqs. (11) and  (12) reveals that harmonic decompositions of three additional expressions are required:

[TABLE]

The decomposition of $\chi_{0}$ at 0PN order, which is consistent with our 0PN treatment of Kepler’s equation, is very straightforward and we obtain (using Eq.(B13) of Pierro et al. (2001))

[TABLE]

with $G_{0}=0$ and $G_{s\neq 0}=J_{s}(se)$, where $J_{s}(x)$ are Bessel functions of the first kind. With a little more work, we further obtain

[TABLE]

with the definitions

[TABLE]

where the $A_{j}^{s}$ are given by Eq. (42b) of Boetzel et al. (2017). Lastly, for $\chi_{2}$ we have

[TABLE]

where we have defined

[TABLE]

We can combine all of these results to find the harmonic decomposition of the waveform. First, we find that

[TABLE]

where we have defined $M_{s}=I_{s}+2L_{s}$, and we note that $M_{s}=M_{-s}$. With this, we then have the full decomposition of $H^{(0)}_{+,\times}(t)$ into harmonics of $l$ and $\lambda$:

[TABLE]

where we have defined

[TABLE]

The $N_{j}$ coefficients scale as $e^{|j|}$ to leading PN order. In this work, we represent $N_{j}$ as the first 20 terms in a low-eccentricity expansion past leading order and at 3PN in the contributions from $W$ (keeping only 0PN contributions from the solution to Kepler’s equation) with $j=[-15,15]$. We provide each of these coefficients to 16 digits of precision in the supplemental material.

Putting the above together to form the GW strain we then have

[TABLE]

where $Q$ is the complex function of the antenna functions and sky location

[TABLE]

and $\mathcal{A}=\frac{m\eta}{R}y^{2}(1-e^{2})$. Notice that in the circular limit the first term in Eq. (33) is not present, as it scales as $G_{s}={\cal{O}}(e^{|s|})$.

We are now in a position to check if the low eccentricity expansion of the $N_{j}$ and $G_{s}$ can be used to faithfully reconstruct the signal. In order to do so, we compute the overlap between the harmonically decomposed waveform keeping $j=[-15,15]$ and $s=[-15,15]$ harmonics in Eq. (33), and a waveform where trigonometric functions of $u$ and $W$ have not been harmonically decomposed. Figure 3 shows these overlaps for a $(10,10)M_{\odot}$, $(10,1.4)M_{\odot}$, and $(1.4,1.4)M_{\odot}$ binary. We see that for a majority of the parameter space considered, the overlap is above 0.98.

For very high initial eccentricities ($e_{0}>0.7)$ and small dimensionless initial separations ($p_{0}<50$), however, the overlap decreases sharply (below the minimum value shown) to values as low as  0.93. This is due to a combination of spectral truncation (i.e. not enough harmonics have been kept in Eq. (33) to faithfully reconstruct the full signal) and inaccuracy of the harmonic amplitudes, as they have been expanded in low eccentricity. Since Bessel functions appear all throughout the harmonic amplitudes with the form $J_{s}(se)$, and expansions of Bessel functions are known to converge very slowly as their argument is near their index, we conclude that the faithful modeling of this small subset of the parameter space would require a considerably higher number of terms in the eccentric expansions of the amplitudes, as well as more harmonics in the sum of Eq. (33).

III.3 Qualitative Analysis of Eccentric Signals

Let us now develop a sense of the relative strength of the harmonics and the structure of the Fourier domain signal before moving on to explicitly deriving such a Fourier-domain waveform. We achieve this through a few different methods. First, we examine the structure of the Fourier response without including the effects of radiation reaction. This reveals the relative strength of the different harmonics and the mapping between the orbital frequencies, $n$ and $\omega$, and the Fourier frequency, $f$. Second we examine the structure of the Fourier response with radiation reaction using both the Fourier transform and a time frequency representation of the signal. Lastly we investigate the relative strength of the $N_{j}$ by plotting their values as function of $j$ and $e$. Knowing which amplitudes are the strongest is key if one wishes to make the model more efficient by careful truncation or selection of which harmonics to include in Eq. (33).

Let us first inspect the Fourier response of a steady state system (i.e. in the absence of radiation reaction so that $\dot{e}=\dot{y}=0$). As a result the orbital angles $l$ and $\lambda$ are given by $l=n(t-t_{0})$ and $\lambda=\omega(t-t_{0})$. The Fourier transform of the signal as written in Eq.(33) is then trivial: the sum of the product of the harmonic amplitude $N_{j}$ or $G_{s}$ and delta functions centered at $sn/(2\pi)$ or $\pm(jn+2\omega)/(2\pi)$, where $s$ and $j$ are the indices of the sums in Eq. (33). We proceed by taking the discrete Fourier transform of the numerical un-decomposed time domain signal in order to get a sense of the distribution of strength of the harmonics, as well as to verify that the result is consistent with Eq. (33).

The left panel of Fig. 4 shows the amplitude of the discrete Fourier transform of the numerically evolved (at 3PN order in the dynamics, keeping only leading PN order contributions in the amplitude) time-domain signal as written in Eq. (11). As predicted, we see spectral lines appearing at Fourier frequencies that are combinations of the two orbital frequencies. In the case of the signals with modest initial eccentricities, we have labeled a number of frequency positions of the different harmonics. The subscript on the index $j$ indicates whether the plus or minus has been taken in the location of the $j$th harmonic given by $\pm(jn+2\omega)/(2\pi)$. In the high eccentricity case, there are too many closely spaced harmonics to clearly indicate which is which. The harmonics of index $j_{+}$ tend to be dominant, followed by the harmonics of $s$, and the harmonics of index $j_{-}$ are clearly subdominant. We choose all sky angles to be $3\pi/7$, and have verified that the results do not qualitatively change much with the choice of sky angle (other than the fact that the harmonics of $s$ vanish in the $\iota\rightarrow 0$ limit).

When we allow the binary to inspiral due to radiation reaction, the orbital eccentricity and the two orbital frequencies change with time. As a result, the Fourier transform is not a simple sum of spectral lines any longer. Since the amplitude is slowly varying compared to the complex exponentials, we can expect the spectral lines of the left panels of Fig. 4 to sweep with the orbital frequencies as the latter increase, forcing different harmonics to overlap and to interference. The right panel of Fig. 4 shows the Fourier spectrum when the effects of radiation reaction are included. We see that the amplitude oscillates rapidly due to harmonic interference as expected.

Since the amplitude is so rapidly oscillating when radiation-reaction is included, it is more enlightening to look at a time-frequency representation of the GW signal. Figure 5 shows the Q-transform of the GW signals considered in the right panel of Fig. 4. A Q-transform is one in which the basis functions are sine-Gaussians with quality factor Q, such that a low (high) value of Q localizes the basis more in time (frequency). With this representation of the signal, we can clearly see the behavior of the different harmonics. Each harmonic takes a distinct track in time-frequency space, as a result of the orbital frequencies increasing in time. We also see that harmonics, besides the $j=0$ one (the only non-vanishing harmonic in the circular limit), decay with time and very rapidly with frequency. The lower left panel shows the time frequency evolution for the high eccentricity case ($e_{0}=0.9$). Although at late times the harmonic structure is clear, this is not so at early times, which is an artifact of the particular choice of Q we made in this figure (a choice of $Q=150$ does well at capturing the oscillatory nature of the late signal since it is more localized in frequency). At early times the signal consists of a series of bursts. In the lower right panel we have plotted the Q-transform of some of the early signal with $Q=20$, to appreciate how much more localized in time the signal is at early times.

Figure 6 presents the value of $N_{j}$ as a function of the index $j$ and the eccentricity $e$ at which they are being evaluated, keeping 20 terms in eccentricity past leading (left) and 16 terms past leading (right). We have set $\eta=0.25$ in these results, but verified that setting $\eta=0.1$ changes the result negligibly. We have set $y$ such that the results correspond to a dimensionless semilatus rectum of 50. In the low eccentricity limit, the $j=0$ harmonic dominates as expected as it is the only non-vanishing harmonic in the circular limit. As the eccentricity is increased, the harmonics become more comparable in strength with the $j>-1$ harmonics clearly dominating. For large eccentricities, there is considerable strength in the larger negative $j$ harmonics. We see that at such high eccentricity there are visibly different features between the two panels suggesting that the eccentricity expansions are beginning to break down.

III.4 Different Time Domain Approximants

Given the increase in the number of relevant ODEs and parameters in the eccentric problem, there are a number of different and new PN-consistent ways to solve the orbital dynamics beyond what is possible in the circular limit. In order to introduce these, let us begin by reviewing the circular TaylorT4 approximant (for a more in-depth discussion of the different circular PN approximants, see Sec. III of Buonanno et al. (2009)). In the circular TaylorT4 approximant, one numerically solves the set of ODEs:

[TABLE]

The numerical solutions to these are then substituted into a circular GW polarization and a DFT can then be taken to obtain a Fourier-domain waveform.

In the eccentric case, however, there are several different ways of solving the relevant orbital dynamics presented in Eq.(9), and each of them is consistent in a PN sense. Given this, we introduce here 3 different time-domain approximants that extend the circular TaylorT4 approximant to eccentric inspirals: TaylorT4t, TaylorT4y, and TaylorT4e. Each of these is constructed by numerically solving a set of ODEs, with the difference being the independent variable $t$, $y$, or $e$ used to parameterized the solution.

The TaylorT4t approximant is obtained by using time $t$ as the independent variable when numerically solving the set of ODEs $(\dot{l},\dot{\lambda},\dot{e},\dot{y})$ for $[\lambda(t),~{}l(t),~{}e(t),y(t)]$. The resulting solutions are then inserted into the appropriate expressions for the full orbital dynamics (Eq. (5)) to obtain $\phi(t)$, and the result of this is then inserted into the expression for $h(t)$ in Eq. (11). Finally, we take a DFT to obtain a Fourier domain waveform.

The TaylorT4e approximant is obtained by using the eccentricity $e$ as the independent variable when solving the orbital dynamics equations. This implies we must first express the ODEs of TaylorT4t in terms of the eccentricity and then numerically solve $(d\lambda/de,dl/de,dt/de,dy/de)$ for $[\lambda(e),l(e),t(e),y(e)]$. This re-writing is obtained via the chain rule and the subsequent re-expansion of the resulting expressions to the appropriate PN order. For example, to obtain $dy/de$, we use $dy/de=\dot{y}/\dot{e}$, where the ratio is re-expanded to 3PN order in $y$. The numerical solution to these ODEs can then be expressed as functions of time by numerically inverting $t(e)$ for $e(t)$. With these substitutions done, one finds $h(t)$ through Eq. (11) and applies a DFT to find the Fourier-domain representation.

The TaylorT4y approximant is obtained by using the PN variable $y$ as the independent variable when numerically solving the set of ODEs $(d\lambda/dy,dl/dy,dt/dy,de/dy)$ for $[\lambda(y),l(y),t(y),e(y)]$. Again in each of these ODEs, any ratio of PN expressions arising due to the use of the chain rule is re-expanded to the appropriate PN order. With solutions for the orbital elements in terms of $y$ in hand, one then numerically inverts $t(y)$ for $y(t)$, and this is used to construct $h(t)$ via Eq. (11) before taking a DFT to find the Fourier-domain representation.

Each of these different ways to solve for the orbital elements is equally valid in the PN approximation, but each solution is slightly different due to the truncation of this series in the ODEs and the subsequent use of numerical solvers. Since our Fourier domain waveform will later rely on the differential equations of the TaylorT4e approximant, it is important to know how faithful TaylorT4e is relative to the TaylorT4t approximant. Such an analysis will reveal the PN limitations of our frequency domain waveform model. We choose to compare TaylorT4e to TaylorT4t because, at least in the circular limit, the latter have been shown to be the most faithful against EOB waveforms Buonanno et al. (2009).

Figure 7 shows the overlap maximized over a time offset and overall phase offset between the TaylorT4e and the TaylorT4t approximants for a $(10,10)M_{\odot}$, a $(10,1.4)M_{\odot}$, and a $(1.4,1.4)M_{\odot}$ binary for a variety of different initial conditions ($e_{0}$ and $p_{0}$). In the equal-mass cases, we see that as the two models agree with each other in most of the initial parameter space, with disagreements only arising at sufficiently large initial eccentricity or small initial separation. To be clear, there has been no small eccentricity expansion or approximations between these two models, with the only source of disagreement coming from the truncation of the PN series when reformulating the ODEs in terms of $t$ or $e$. Both models are obtained purely numerically and without any harmonic decomposition.

In the unequal mass case, the disagreement between the TaylorT4e and TaylorT4t approximants is much more pronounced. This is because of the well-known fact that the PN expansion is poorly behaved for unequal masses Yunes and Berti (2008), but this feature seems to become enhanced for even moderate initial eccentricities. Since the Fourier-domain model we will develop later solves for the orbital dynamics as formulated in the TaylorT4e approximants, we expect that it will also disagree with TaylorT4t to the level displayed in Fig. 7.

IV Eccentric GWs in the Fourier Domain

In this section we begin by reviewing the Fourier domain model that was derived in paper 1. Following the prescription laid out there, we extend the model to 3PN order in the Fourier phase. Lastly we describe how other frequency domain approximants could be derived in analogy to the different time domain approximants discussed in Section III.4.

IV.1 Review of Newtonian model

In paper 1 we computed the Fourier response via the SPA of a time domain signal whose time domain harmonic decomposition took the general form:

[TABLE]

Here the amplitudes $A_{j}$ are functions of eccentricity and vary on a slower timescale (the radiation reaction timescale) than the orbital angle $l$ (which varies on the orbital timescale). This makes it appropriate for the application of the stationary phase approximation (SPA) which is described in detail in Bender and Orszag (1999). The time domain signal in Eq. (36) is composed entirely of harmonics of $l$ because at Newtonian order $\lambda=l$, a result of periastron precession effects entering at 1PN order. Application of the SPA to the above signal yields

[TABLE]

where the Fourier phases are given by

[TABLE]

The stationary phase condition relates the Fourier frequency, $f$, to time domain quantities through

[TABLE]

Clearly then, it is necessary to express the functions $t^{\ast}_{j}$ and $l$ that appear in the Fourier phases in terms of some time domain quantity to evaluate the model at any given frequency. In paper 1, we chose to use the orbital eccentricity as that time domain quantity, because the orbital dynamics at Newtonian order can be solved exactly in terms of this quantity, i.e. one is able to exactly solve the Newtonian order $dn/de$ equation to obtain $n(e)$. This solution can then be used in the integrals that define the phase functions

[TABLE]

which admit exact solutions in the form of hypergeometric functions. Evaluation of these functions, however, is computationally expensive, and thus, in paper 1 we used a simple Taylor expansion of those hypergeometric functions about small eccentricity, which we found was sufficient to faithfully model the Fourier phases, even to eccentricities as high as $\sim 0.9$. The last step to evaluate the model at a given frequency was to numerically invert the stationary phase condition to find $e(n)$. The stationary phase condition can be more suggestively written as

[TABLE]

and this is a relatively cheap inversion because the solution is clearly independent of $j$ and, in the Newtonian limit, also system independent.

The key result which allows for the model to be readily extended to 3PN order is that the Taylor expansions in orbital eccentricity of the hypergeometric functions appearing in the phase were generally sufficient to faithfully represent the signal even for large eccentricities. Only the inversion of $n(e)$ behaves very poorly for even moderate eccentricities. Thus, the route for extension to higher PN is clear:

• (i)

solve for the frequency evolution as a function of eccentricity in an expansion in both a PN parameter and the orbital eccentricity,

• (ii)

use this solution to express the integrands appearing in the phase functions as expansions of those same parameters and integrate the series in eccentricity,

• (iii)

numerically invert the stationary phase condition to relate frequency to the eccentricity.

IV.2 Eccentric 3PN Fourier Domain Model

Let us now derive the Fourier domain model to 3PN order, which we will call TaylorF2e as it closely follows the quasi-circular TaylorF2 model, and corresponds to solving the orbital dynamics where the independent variable is $e$.

Following the general scheme laid out in paper 1, we begin by applying the SPA to the harmonic decomposition, which was derived in Eq. (33) of Sec.III.2. We reiterate the harmonic decomposition here, explicitly writing out the complex conjugate term:

[TABLE]

Applying the SPA to each of the oscillating terms separately (see Appendix A of paper 1 for a brief description of the SPA as applied to individual harmonics of the orbital phases), we find

[TABLE]

for the first term,

[TABLE]

for the second term and

[TABLE]

for the last term.

The stationary phase condition reveals which terms in the sum contribute to the Fourier transform. In the case of harmonics of $\dot{l}$, only positive $s$ terms contribute as $\dot{l}>0$. However, in Eqs.(44) and (45), the appearance of $(j\ddot{l}+2\ddot{\lambda})$ in the denominator could lead to catastrophes (instances where the amplitude diverges at a subset of frequencies) in the event that this term vanishes. To determine which terms could possibly lead to these catastrophes, we re-express the term appearing in the denominator of the amplitudes as

[TABLE]

We can rearrange the above to find the $j$ for which there are catastrophes,

[TABLE]

where

[TABLE]

Are there catastrophes in the SPA? Figure 8 shows $\chi(e)$, which clearly varies from $2/3$ at $e=0$ to about $0.05$ at $e=0.9$. If we bound the 1PN pericenter precession quantity $k$ to $0<k<1/2$, we see that it is possible to have a catastrophe in the $j=-3$ case when the binary is sufficiently close (i.e. when $k\sim 1/2$) and circular ($\chi(e)\sim 2/3$). However, since accuracy in the amplitude is less important than accuracy in the phase for GW data analysis, we proceed by keeping only the leading PN order corrections in $\ddot{l}$ and $\ddot{\lambda}$. At leading PN order, $\ddot{l}=\ddot{\lambda}$, and clearly, there are only catastrophes when $j=-2$. But this harmonic corresponds to a low-frequency contribution that has a negligible amplitude contribution, as we can see from the stationary phase condition and Figure 6. In fact, neglecting this harmonic altogether affects the matches shown in Fig. 3 negligibly, so henceforth we neglect the $j=-2$ harmonic, which ensures our model never hits catastrophes.

But exactly which $j$ indices contribute to the Fourier transform? Inspection of the stationary phase conditions (Eqs. (45c) and (44c)) does not clearly reveal the answer to this question. In order to find the $j$ indices that do contribute, let us begin by rearranging the stationary phase conditions to set a condition on $j$ for which $f$ is positive:

[TABLE]

where we have used $\dot{l}=n$ and $\dot{\lambda}=n(1+k)$. We see that as long as $k\leq 1/2$ then the stationary point corresponding to the “-” superscript is satisfied when $j>-2$ and that corresponding to the “+” superscript when $j<-3$.

Clearly then, in order to understand which indices $j$ contribute to the Fourier transform, we must first determine which values of $k$ are reasonable within the PN approximation. To do so, let us compare two constraints that we might impose on our waveforms:

• (i)

that the periastron velocity be small $v_{p}<1/3$,

• (ii)

that the rate of periastron advance be slow $k<1/2$.

Each of these constrains the final separation of the binary, but since our results are typically quoted in terms of initial dimensionless semi-latus rectum, $p$, it is natural to consider what lower limit on $p$ each of the above constrains implies. Constraining the periastron velocity ($v_{p}<1/3$) leads to $p>9(1+e)^{2}$ [or equivalently $y<(3(1+e))^{-1}$], while constraining the advance of periastron ($k<1/2$) leads to $p>6$. Since the periastron velocity constraint is more stringent, henceforth we use it to stop our waveform evolution. Semi-latus recta below this number would incur into a regime were the PN approximation may be inaccurate Yunes and Berti (2008).

Restricting $v_{p}<1/3$ ensures that $k<1/2$, and thus, we see that the stationary phase condition of Eq. (44c) is satisfied only when $j>2$ , while that of Eq. (45c) is satisfied when $j<-2$ (with the understanding that we have dropped the $j=-2$ harmonic). In practice, we have found that the stationary phase condition for the $j=-3$ harmonic can become double valued, instead of adding more complication by treating this harmonic with more care, we simply neglect it as we have verified that it does not affect the results of Fig. 3. Putting everything together then we have a frequency response in the SPA given by:

[TABLE]

where

[TABLE]

Following paper 1 and with this in hand, we now seek solutions to express the functions appearing in the phase ($t,l,\lambda$) in terms of the orbital eccentricity.

Let us begin by seeking a relation between the orbital frequency and the eccentricity. Since either of the two orbital frequencies ($n$ or $\omega$) can be expressed in terms of the PN parameter $y$, this amounts to seeking the mapping $y(e)$. Composing $dy/de$ using the chain rule (i.e. $dy/de=\dot{y}/\dot{e}$) and expanding to 3PN order in $y$ yields

[TABLE]

whose PN solution is

[TABLE]

for certain functions $d_{n}(e)$. Inserting Eq. (54) into Eq. (53) and expanding to 3PN order (now an expansion in the parameter $y_{0}$) gives a set of coupled ordinary differential equations for the $d_{n}(e)$ functions, namely

[TABLE]

where we have split the 3PN contribution into a component that is proportional to the natural logarithm of $y_{0}$ and one that is not proportional to it, such that the 3PN part of the solution to $y(e)$ takes the form $(d_{6}+d_{6l}\ln{y_{0}})y_{0}^{6}$.

Given that

[TABLE]

we see that solving any of Eqs. (55) amounts to solving a differential equation of the form

[TABLE]

for some functions $g(e)$. The above differential equation has the solution

[TABLE]

for some constant $C$. The constants arising from solving for any of the coefficients $d_{n}$ are determined by requiring that at each PN order $y(e_{0})=y_{0}$. Clearly, in solving for $d_{0}$, $g(e)=0$ and the solution is simply

[TABLE]

The solutions for the higher PN coefficients $d_{n}$ are not so simple as the form that $g(e)$ takes becomes a more complicated function of eccentricity. In order to solve the equations, we then expand Eq. (IV.2) in eccentricity, which allows us to integrate the series analytically. The $d_{n}(e)$ functions and the constants of integration $C_{n}$ to order ${\cal{O}}(e_{0}^{40})$ past leading are provided in the supplemental material.

With the solution of $y(e)$ in hand, we are now in a position to solve for the phase functions $t(e),\lambda(e),$ and $l(e)$. We begin by substituting the solution above for $y(e)$ into the integrands for each of the phase functions

[TABLE]

In order to carry out the integration above we expand the integrand to 3PN in $y_{0}$. We then expand the coefficients at each order in $y_{0}$ to a given order in $e$. Note that there is no restriction on the order to which we expand any PN coefficient (coefficients of $y_{0}$ in each of the integrands) in eccentricity other than the order to which we have obtained $y(e)$. We could, for example, keep $\mathcal{O}(e^{100})$ terms in each coefficient, but with the inclusion of so many terms the model would become very computationally inefficient to evaluate. On the other hand, if we keep too few eccentric corrections, the model becomes inaccurate. In Section V.1 we introduce a scheme to pick the order in eccentricity to which we keep the expansions of each $y_{0}$ coefficient, and then show in Section V.2 how the eccentricity expansion affects the validity of the waveform for systems with different initial eccentricities.

After carrying out the integration, we have the phase functions in terms of the orbital eccentricity,

[TABLE]

where the functions $t^{(n)},l^{(n)}$ and $\lambda^{(n)}$ are series in eccentricity (and also depend on the constants of integration appearing in $y(e)$), which are provided in the supplemental material. With these solutions in hand, the phases in Eq. (IV.2) are functions of only the orbital eccentricity. The final step then is to numerically invert the appropriate stationary phase condition to relate the Fourier frequency to the orbital eccentricity. The specifics of the evaluation of TaylorF2e are described in Sec. V.1. Section V.2 shows that this waveform model is as faithful to TaylorT4t as the TaylorT4e model is, giving good agreement for even moderate eccentricities.

IV.3 Alternative Frequency Domain Models

In Section III.4 we laid out different possible time domain models based on the independent variable appearing in the set of ODEs that one solves numerically to specify the orbital dynamics of the system. In the frequency domain model derived in the previous section, we solved the set of ODEs with the eccentricity $e$ as the independent variable, and thus, the Fourier domain TaylorF2e model corresponds to the SPA version of the TaylorT4e time domain approximant. In this section, we wish to explore schematically how one would go about obtaining other frequency domain approximants related to the other two time domain approximants TaylorT4t and TaylorT4y.

TaylorF2y would be obtained by first solving the set of ODEs $(dl/dy,d\lambda/dy,dt/dy,de/dy)$ analytically within the PN approximation (and possibly in a low-eccentricity expansion). With this at hand, one would then invert $y(e)$ as given in Eq. (54) to obtain $e(y)$, which would be a series in both $y$ and $y_{0}$. This solution would then be plugged into the integrands that define the phase functions and integrated to give $l(y)$, $\lambda(y)$, and $t(y)$. The different stationary phase conditions could be numerically inverted to relate $y$ to the Fourier frequency $f$; one could attempt to invert the stationary phase condition perturbatively as well, but the error incurred would have to be investigated with the waveform in hand. We expect, however, that the resulting TaylorF2y approximant would have only limited validity in $e$, because even at Newtonian order the inversion of $y(e)$ is transcendental and very poorly behaved even for eccentricities as low as $0.3$ (see Appendix C of paper 1 where we investigated this inversion at length at Newtonian order).

TaylorF2t would be obtained by first solving the set of ODEs $(\dot{l},\dot{\lambda},\dot{e},\dot{y})$ analytically within the PN approximation (and possibly in a low-eccentricity expansion) to find both $e(t)$ and $y(t)$. With this at hand, one would express the integrals for $l$ and $\lambda$ in terms of only the $t$ variable, but to obtain both, one would first have to invert the solutions found in TaylorF2e and TaylorF2y to obtain $t(y)$ and $t(e)$, or else come up with some scheme to solve $\dot{e}$ and $\dot{y}$ analytically. Upon substitution of these solutions into $\dot{l}$ and $\dot{\lambda}$, one would integrate to obtain $l(t)$ and $\lambda(t)$ analytically. Again, the stationary phase condition would be numerically solved to find $t(f)$ for any given $j$ or $s$ harmonic index, or alternatively one could try to perturbatively invert the stationary phase condition, but this could incur too much error and would have to be investigated with the approximant in hand.

V Implementation and Validation

In this section we validate the TaylorF2e model introduced in in Sec. IV against the TaylorT4e model. To quantify the agreement between two models $h_{1}$ and $h_{2}$ it is customary to introduce the match

[TABLE]

which is nothing but the overlap between two waveforms maximized over the extrinsic parameters $t_{c}$, $l_{c}$, and $\lambda_{c}$. Note that this differs from the match in the quasicircular case, which is only maximized over $\phi_{c}$ and $t_{c}$, because in the eccentric case there is an extra phase angle associated with pericenter precession. We employ the implicit maximization scheme developed in Section 3 of Paper 1 in order to maximize over these parameters.

V.1 Implementation

Let us begin by clearly describing how theTaylorF2e model is generated and implemented. As explained in Sec. IV, the TaylorF2e model can be deployed in slightly different ways, which in essence depend on the following choices:

• (i)

The order of the eccentricity expansion of each of the coefficients appearing in the phase functions of Eqs. (63).

• (ii)

The technique by which we invert the stationary phase condition to relate the orbital eccentricity to the Fourier frequency.

• (iii)

The method of obtaining the eccentricity dependence on the orbital frequency.

Let us begin by addressing item (i). The expansion of the phase functions in small eccentricity introduces an error in our waveform model, and thus, we wish to determine the optimal expansion order such that the resulting waveform is accurate enough to lead to sufficiently high matches against a TaylorT4e model. Let us then begin by introducing the following notation for our ODEs $(dl/de,d\lambda/de,dt/de)$ that lead to Eqs. (63):

[TABLE]

In each of the PN coefficients appearing above $[\lambda_{i}^{\prime}(e),l_{i}^{\prime}(e),t_{i}^{\prime}(e)]$, we have substituted in our solution for $y(e)$, which is valid to order $\mathcal{O}(e^{40})$ past leading, but we have not yet expanded these coefficients in eccentricity. We now wish to determine the order in eccentricity to which we should expand each of these coefficients.

Let us then introduce the following measure of error in each of the PN coefficients

[TABLE]

where the overhead bar stands for the Taylor expansion in small eccentricity to a given order. This is analogous to the Newtonian-order measure of error introduced in paper I

[TABLE]

where $\psi_{\rm exact}$ was our exact (in eccentricity) phase and $\psi_{\rm err}$ had a fixed amount of error introduced through $\epsilon$ via the above equation. We found that for moderately eccentric systems, the match began to drop for values of $\epsilon\sim 10^{-5}$, as shown in Table III of paper I.

Returning to our measures of error in Eq. (66), we must now choose a tolerance for $\gamma_{n}^{\lambda,l,t}$ that is acceptable for our model, i.e. that leads to a sufficiently small loss of match. Guided by the error investigation of paper 1, one may be tempted to set a uniform tolerance of $10^{-5}$. One may expect, however, that more error could be tolerated in the coefficients appearing at higher PN order, since they are multiplied by the small parameter $y_{0}$. Thus, at each PN order we impose a maximum error tolerance $\gamma_{n}^{\lambda,l,t}(e)\leq\epsilon_{n}$, where

[TABLE]

for some constant $\epsilon$. The factor of $(y_{0})$ allows for more error to be accumulated by coefficients at higher PN order, since their contribution to the phase is scaled by the small parameter $y_{0}$.

The tolerance condition described above clearly depends on $y_{0}$ (or $p_{0}$), $\eta$, $e_{0}$, and the constant $\epsilon$. We have investigated two different constraints, one in which we set $\epsilon=10^{-5.5}$, $\eta=0.25$, $p_{0}=12$, and $e_{0}=0.7$, and one in which we set $\epsilon=10^{-5}$, $\eta=0.25$, $p_{0}=18$, and $e_{0}=0.4$. The first set of conditions forces us to keep more terms in eccentricity than the second set, which we may then consider a conservative choice. With these choices, we then generate two waveforms, one with more terms kept in eccentricity (TaylorF2e+) and with less terms kept (TaylorF2e-), and then we investigate the effects in the match. Figure 9 shows $\gamma_{n}^{\lambda,l,t}$ as a function of eccentricity for the conservative choice. We see that $\gamma^{\lambda,l,t}_{n}$ is at least nearly bounded by the threshold imposed by Eq. (68) (the spikes correspond to zero-crossings on a logarithmic scale). In Table 1 we list the maximum order past leading to which each of our coefficients are expanded in eccentricity for TaylorF2e+ and TaylorF2e-. Note that while many terms have been kept in each of the phase functions, one can collect the resulting expression in like powers of eccentricity. For example, after such a collection, there are only 26 terms in eccentricity to be evaluated for $\lambda(e)$ in the TaylorF2e- case.

To address item (ii), let us first recall the different stationary point inversions that we must achieve to relate $f$ to $e$. The equations that must be inverted are

[TABLE]

where the sign in the first equation depends on the index $j$ as detailed in Sec. IV. The second equation is particularly simple to solve since it only requires the inversion of $n(e)$. That is, defining the inverse function $\kappa$ such that $\kappa[n(e)]=e$, one can relate $f$ to $e$ for the $s$ indices via

[TABLE]

Of course, the inverse function $\kappa$ must be obtained numerically, but this formulation has the advantage that once $\kappa$ is obtained, it can be used for any index $s$. The inversion of the equation for the $j$ harmonics is more complicated. First, let us rewrite the stationary condition as

[TABLE]

We then employ the secant method Press et al. (2007) to numerically solve Eq. (71), where we recognize that the term $kn$ is of higher PN order than $(j\pm 2)n$ and we leverage this to specify the approximate initial guess $e_{j}(f)^{\rm guess}=\kappa[2\pi f/(j\pm 2)]$. Indeed, this initial guess makes the secant method typically converge after 4 iterations to an accuracy of $10^{-8}$. We expect that this error tolerance could be relaxed.

Let us conclude this subsection by discussing item (iii), which in any case is required to invert the stationary phase conditions, as it needs a solution for $y(e)$. We obtain this solution by solving $\dot{y}$ and $\dot{e}$ and interpolating these solutions to compose $n(e)$ and $\omega(e)$. Given that we have obtained $y(e)$ numerically, this is then in turn used in the various harmonic amplitudes, which enter the SPA (in $\mathcal{A}$, $N_{j}$, and $\ddot{l}$). While somewhat computationally costly, we estimate that there are considerable speedups to be gained through numerical techniques and careful truncation of the sum appearing in Eq. (IV.2). We leave a thorough investigation of computational cost to future work. Alternatively we could use the solution found in Eq. (54), which is purely analytical. However, in the interest of being consistent with paper 1, we chose here to solve $\dot{y}$ and $\dot{e}$ numerically instead.

V.2 Validation

In this subsection we validate the TaylorF2e model by computing the match between it and the TaylorT4t model, which was introduced in Sec. III. Figure 10 shows the match for a $(10,10)M_{\odot}$, $(10,1.4)M_{\odot}$, and a $(1.4,1.4)M_{\odot}$ binary, with TaylorF2e+ on the left and TaylorF2e- on the right. In these match calculations there are three main sources of loss in match:

• (i)

the truncation of the summations in Eq. (IV.2) (explored in Fig. 3),

• (ii)

the PN disagreement between models obtained through solving different PN valid sets of ODEs (explored in Fig. 7), and

• (iii)

the truncation of the eccentricity expansions in the phase functions.

Let us discuss each of these potential sources of error in turn.

In order to address item (i), let us return to Fig.  3, which shows the overlap between two numerically calculated time-domain models (both TaylorT4t) with one model kept to all harmonics and the other restricted to $j=[-15,15]$ and $s=[1,15]$. Comparison of that figure to Fig. 10 reveals that spectral truncation is not the main source of loss in match between TaylorT4t and TaylorF2e. The loss in match due to spectral truncation (as shown in Fig. 3) occurs for much greater initial eccentricities and smaller initial semi-latus recta than the loss in match that is seen between TaylorF2e and TaylorT4t (as shown in Fig. 10).

When we compute the match between TaylorF2e+ and TaylorT4t we nearly exactly recover the results of taking the match between TaylorT4t and TaylorT4e. This is not a surprising result as TaylorF2e is built by solving the ODEs formulated with $e$ as the independent variable, as is done in TaylorT4e. This suggests that, at least when keeping many terms in eccentricity, the dominant source of loss in match is due to the PN differences between the ODEs governing the orbital dynamics, when re-expressing them in terms of a different independent variable. This result suggests that a more faithful model would require PN ODEs at higher PN order for the orbital dynamics.

Let us now return to item (iii). When comparing the left and right panels of Fig. 10, we see that for the $(10,10)M_{\odot}$ and $(1.4,10)M_{\odot}$ systems, there is nearly no change in the match whether we use TaylorF2e+ or TaylorF2e-. However, in the $(1.4,1.4)M_{\odot}$ case, there is significant loss in match as the initial eccentricity exceeds $0.5$. This suggests that for less massive systems, it is important to keep more terms in the eccentricity expansions. This conclusion is also supported by Fig. 11, which shows again the match between TaylorF2e+/- (left/right) and TaylorT4t but this time for a $(5,5)M_{\odot}$ system. In this case, the loss in match due to truncation of the eccentricity expansions becomes considerable for close initial separations and initial eccentricities around $0.6$.

The matches appearing in Fig. 10 misleadingly suggest that the model is unfaithful for a significant amount of the parameter space, but this is only an artifact of the particular parameter space region we chose to calculate the match over. For initial semi-latus recta that are larger than those considered here, the model is very faithful with matches larger than 0.98. Let us then study if eccentric corrections to the model matter at larger initial separations than those considered in Fig. 10. In order to roughly determine which combinations of mass, initial eccentricity, and initial semi-latus recta lead to significant eccentric corrections, we employ a signal-to-noise ratio (SNR) measure, defined by

[TABLE]

Figure 12 shows the SNR normalized by the largest SNR present in the parameter space explored for the three systems considered in Fig. 10. We see that particularly in the less massive case there is significant SNR for systems with larger initial semi-latus recta than those considered in Fig. 10. This suggests that our TaylorF2e model is useful for systems in a larger parameter space than considered in Fig. 10. We leave a more thorough investigation of the region of parameter space in initial eccentricity and semi-latus recta where eccentric effects matter to future work.

V.3 Error Analysis in Different Components of TaylorF2e

In this section we wish to investigate which of the major analytic components which enter the TaylorF2e contribute to the most error in match. Namely we wish to know which of the functions which appear in the phase ($\lambda(e)$, $l(e)$, and $t(e)$) give the most loss in match, so that if one wished to improve the model they would have some guidance as to where to start. In order to do so we obtain numerical solutions for the phase functions by numerically solving $\dot{\lambda}$, $\dot{l}$, $\dot{y}$, and $\dot{e}$. With these solutions in hand we invert $e(t)$ to obtain $t(e)$ and use this to map $\lambda$ and $l$ to the eccentricity domain which allows them to be used in our TaylorF2e model.

In Figure 13 we show the value of the match for a $(10,10)M_{\odot}$ system where we have kept the phase functions all numerical, or just one of any of them numerical. We see, as expected, when all of the phase functions are numerical we recover the result of Figure 3. The error in match is due to the truncation of harmonics in Eq.(IV.2). The best increase in match (besides keeping all phase functions numerical) occurs when $l(e)$ is computed numerically. Interestingly, the match when $\lambda(e)$ is kept numerical is slightly worse than had the analytic expression for $\lambda(e)$ been used. This is because there is cancellation in error between the analytic expression for $t(e)$ and $\lambda(e)$.

VI Conclusions & Future Work

We have extended the eccentric Fourier-domain Newtonian model of paper 1 to 3PN order. This model combines a harmonic decomposition of the time domain signal, the stationary phase approximation, and PN theory. We have thoroughly validated the model and shown that it is faithful to very high initial orbital eccentricities ($e_{0}\sim 0.7$), but of course this depends on the particular system and the initial separation at which the initial eccentricity is defined. We expect that the model is valid to even higher eccentricities in relevant regions of ($e_{0}$, $p_{0}$) parameter space.

In addition to this main result, we have validated the harmonic decomposition of the eccentric time-domain signal in terms of a sum of harmonics of the two orbital frequencies with amplitudes that are expanded in small eccentricity. We also investigated the importance of higher PN corrections in the overall amplitude of the time domain model, and we investigated the agreement of different PN consistent time domain approximants. We found that the harmonic decomposition presented here is valid for very large initial eccentricities ($e_{0}\sim 0.9$), and that in the equal mass case the effects of higher PN corrections in the amplitude are negligible. Interestingly, we found that different PN consistent ways of solving for the orbital dynamics lead to time-domain approximants that disagree at moderate initial eccentricities and relatively close initial separation – an effect which is enhanced in the unequal mass limit. However, we expect that this result only excludes faithful modeling in a small region of the relevant parameter space.

Given our model’s 3PN accuracy, we believe it is sufficient for preliminary data analysis applications. However, first we must address the efficiency of the model. In this work we sought to validate the model, but we did not implement it in the fastest possible way. We expect that through a careful evaluation of the Fourier phase, the inversion of the stationary phase condition, and neglecting certain low SNR harmonics, the model can retain its accuracy and become very efficient.

With such an efficient and 3PN accurate eccentric model we could seek to answer some important questions about eccentric inspirals using Bayesian techniques. We plan to investigate the region of $(M,\eta,p_{0},e_{0})$ parameter space where eccentric corrections to the signal are important (i.e. where in this parameter space we could differentiate between a quasi-circular signal and an eccentric one). We also plan to investigate the ability to measure parameters and any biases of the model.

We also expect that, using our model, we can constrain certain modified theories of gravity. If such a modified theory of gravity modifies the eccentricity evolution of a binary, these corrections could be incorporated into the model presented here. With these corrections due to a modified theory present in the model, we could investigate the ability of detectors to constrain these corrections and test General relativity perhaps more stringently than possible with a quasi-circular model.

Lastly, in order to fully model an eccentric signal we require the incorporation of effects due to the merger and the ringdown. The model presented here could be used to create an inspiral-merger-ringdown (IMR) type hybrid model. The only IMR model which incorporates orbital eccentricity is that of Huerta et al. (2018). It will be useful to have several IMR type models with which to extract signals in the future for validation purposes, as is the case for quasi-circular coalescences.

Acknowledgments

B. M. was supported by the Joan L. Dalton Memorial Fellowship in Astronomy from Montclair State University. B. M. and N. Y.  also acknowledge support from NSF Grant No. PHY-1759615, as well as NASA grants NNX16AB98G and 80NSSC17M0041. We thank Antoine Klein for useful conversations and providing us with the expressions for the tail terms. We also thank Travis Robson and Alejandro Cárdenas-Avendaño for useful conversations. Computational efforts were performed on the Hyalite High Performance Computing system, which is supported by University Information Technology at Montana State University.

Appendix A 2PN Polarizations in Modified Harmonic Gauge

Here we provide the necessary corrections in order to obtain Eq.(III.1) using the expressions in Gopakumar and Iyer (2002). In the latter reference, the authors list the amplitude corrections as

[TABLE]

where $\xi=mn$ and we have introduced an overbar to differentiate these amplitudes from those listed in Eq.(III.1). Because gauge differences arise at 2PN order, only the 2PN amplitude contains differences due to gauge. First, we must re-express the above in terms of $y$. At 2PN order, $y$ and $\xi$ are related by

[TABLE]

Inserting the above into the expressions for $h_{+,\times}$ in Eq. (A), expanding to 2PN order in $y$, and collecting like terms in $y$ such that it takes the form

[TABLE]

gives the unbarred amplitudes as a function of the barred ones

[TABLE]

In order to express these results in MH gauge we must substitute Eq.(II.2) into the expression for $\bar{H}^{(0)}_{+,\times}$. This will produce 2PN corrections that when added to $\bar{H}^{(2),\mbox{\tiny ADM}}_{+,\times}$ yields $\bar{H}^{(2),\mbox{\tiny MH}}_{+,\times}$. We have then

[TABLE]

where the $\Lambda_{+,\times}$ are given by

[TABLE]

We these expressions, we now have all that we need to express the polarization in MH gauge and in terms of the parameter $y$.