Ready-to-use Fourier domain templates for compact binaries inspiraling along moderately eccentric orbits

I Introduction

Observations of GW events by the advanced LIGO and VIRGO GW interferometers are ushering in the era of GW astronomy Abbott et al. 2016a; Acernese et al. 2015. These GW events include merging black hole (BH) binaries and an inspiraling neutron star (NS) binary Abbott et al. 2016b; Abbott et al. 2016c; Abbott et al. 2017a; Abbott et al. 2017b; Abbott et al. 2017c; Abbott et al. 2017d; The LIGO Scientific Collaboration and the Virgo Collaboration 2018. Several scenarios that include long-lived (galactic) field binaries, star clusters, galactic nuclei and active galactic nuclei can produce these observed GW events Belczynski et al. 2016; Park et al. 2017; Hoang et al. 2018; McKernan et al. 2018; Kremer et al. 2018. Fortunately, it may be possible to extract valuable information about the astrophysical origins of GW events in the near future. This requires accurate GW measurements of the spin-orbit misalignment or the orbital eccentricities of these GW events Rodriguez et al. 2016a; Chen and Amaro-Seoane 2017; Nishizawa et al. 2017. Using both frequency and time domain inspiral-merger-ringdown (IMR) waveforms, residual orbital eccentricities of the first two GW events were restricted to be below $0.15$ when these binaries entered aLIGO frequency windowAbbott et al. 2016d; Huerta et al. 2017. Strictly speaking, the so far detected GW events do not exhibit any observational signatures of residual orbital eccentricities and are faithfully captured by IMR templates associated with compact binaries merging along quasi-circular orbits.

However, there exists a number of astrophysical scenarios that can produce GW events with non-negligible eccentricities in the frequency windows of ground-based GW detectors. Dense star clusters like the ubiquitous globular clusters are the most promising sites to form aLIGO relevant compact binaries with non-negligible orbital eccentricities Rodriguez et al. 2016b. A recent realistic modeling of globular clusters that involve general relativistic few body interactions provided non-negligible fraction of BH binaries with eccentricities $>0.1$ as they enter the aLIGO frequency window Samsing et al. 2014; Samsing and Ramirez-Ruiz 2017; Rodriguez et al. 2018a; Samsing 2018; Rodriguez et al. 2018b; Kremer et al. 2018. Additionally, there exists a number of other astrophysical scenarios that can force stellar mass compact binaries to merge with orbital eccentricities. This include GW induced merger during hyperbolic encounters between BHs in dense clusters O’Leary et al. 2016 and mergers influenced by Kozai effect in few body systems as explored in many detailed investigations (see Ref. Randall and Xianyu 2018 and references therein). Further, a very recent investigation pointed out that less frequent binary-binary encounters in dense star clusters can easily produce eccentric compact binary coalescence Zevin et al. 2018. These detailed investigations suggest that it may be reasonable to expect GW events with non-negligible orbital eccentricities in the coming years. Non-negligible orbital eccentricities may be helpful to improve the accuracy with a network of GW interferometers to constrain parameters of compact binary mergers Gondán et al. 2018; Gondán and Kocsis 2018. Moreover, massive BH binaries in eccentric orbits are of definite interest to maturing Pulsar Timing Arrays and the planned Laser Interferometer Space Antenna (LISA) Burke-Spolaor et al. 2018; Bonetti et al. 2018.

To derive our eccentric approximant, we extend the post-circular scheme of Ref. Yunes et al. 2009 to higher PN orders. This scheme involves expanding the Newtonian accurate $h_{\times}$ and $h_{+}$ as a power series in orbital eccentricity that requires analytic solution to the classic Kepler equation. We extend such a Newtonian approach by invoking a recent effort to solve analytically PN-accurate Kepler equation in the small eccentricity limit Boetzel et al. 2017. This detailed computation also provided analytic 1PN-accurate amplitude corrected expressions for $h_{\times}$ and $h_{+}$ as a sum over harmonics in certain mean anomaly $l$ of PN-accurate Keplerian type parametric solution Boetzel et al. 2017. Additionally, the above PN-accurate decomposition explicitly incorporated the effect of periastron advance on individual harmonics, numerically explored using PN description in Ref. Tessmer and Gopakumar 2007. We combine such 1PN-accurate amplitude corrected $h_{\times}$ and $h_{+}$ expressions that incorporated eccentricity contributions to sixth order at each PN order with the two beam pattern functions, $F_{\times}$ and $F_{+}$, to obtain fully analytic time domain GW response function $h(t)$. Our eccentric TaylorF2 approximant is obtained by applying the method of stationary phase approximation to such an analytic $h(t)=F_{+}h_{+}+F_{\times}h_{\times}$ expression.

To obtain analytic expressions for several Fourier phases at their associated stationary points of $h(t)$, we require additional PN-accurate expressions. This involves deriving 3PN-accurate expression for the time eccentricity $e_{t}$, present in the 3PN-accurate Kepler Equation Memmesheimer et al. 2004, as a bivariate expansion in terms of orbital angular frequency $\omega$, its initial value $\omega_{0}$ and $e_{0}$, the value of $e_{t}$ at $\omega_{0}$. This lengthy computation extends to 3PN order, the idea of certain asymptotic eccentricity invariant at the quadrupolar order, introduced in Ref. Królak et al. 1995, and extended to 2PN in Ref. Tanay et al. 2016. In fact, we adapted the approach of Ref. Tanay et al. 2016 by employing the appropriately modified 3PN-accurate $d\omega/dt$ and $de_{t}/dt$ expressions of Refs. Arun et al. 2009; Klein et al. 2018 to obtain 3PN-accurate bivariate expression for $e_{t}$. A careful synthesis of the above listed PN-accurate expressions lead to a fully analytic frequency domain TaylorF2 approximant that included 1PN-accurate amplitude corrections and 3PN-accurate Fourier phases. An additional feature of our approximant is the inclusion of periastron advance effects to 3PN order. To explore GW data analysis implications of these features, we perform preliminary match computations Damour et al. 1998. We conclude that the influences of periastron advance are non-negligible for moderately eccentric binaries, especially in the aLIGO frequency window. This observation should be relevant while constructing IMR waveform family for compact binaries merging along moderate eccentric orbits.

II Post-circular extensions to circular inspiral templates

We begin by reviewing two key efforts to include the effects of orbital eccentricity onto the circular inspiral templates Królak et al. 1995; Yunes et al. 2009. This involves listing in Sec. II.1 the steps that are crucial to compute analytic frequency domain GW response function with quadrupolar amplitudes and PN-accurate Fourier phase in some detail. Various lengthy expressions, extracted from Refs. Boetzel et al. 2017; Arun et al. 2009; Klein et al. 2018, are listed in Sec. II.2 that will be crucial to compute the time domain response function for eccentric binaries while incorporating effects of periastron advance, higher order radiation reaction and amplitude corrections.

Following Thorne 1987, we may express the GW interferometric response function as

where $F_{\times,+}\left(\theta_{S},\phi_{S},\psi_{S}\right)$ are the two detector antenna patterns. These quantities depend on $\phi_{S},\theta_{S}$, the right ascension and declination of the source, and certain polarization angle $\psi_{S}$ Thorne 1987. For eccentric inspirals, the explicit expressions for the quadrupolar order GW polarization states, $h_{\times}$ and $h_{+}$, are given by Eqs. (3.1) of Ref. Yunes et al. 2009. It is rather straightforward to express these Newtonian accurate expressions as a sum over harmonics in terms of the mean anomaly $l$. The resulting expressions read

where $D_{L}$ denotes the luminosity distance while the symmetric mass ratio $\eta$ of a binary consisting of individual masses $m_{1}$ and $m_{2}$ is defined to be $\eta=(m_{1}\,m_{2})/m^{2}$ while the total mass $m=m_{1}+m_{2}$. Further, we use the commonly used dimensionless PN expansion parameter $x=\,\left(\frac{G\,m\,\omega}{c^{3}}\right)^{2/3}$ where $G$, $c$ and $\omega$ are the gravitational constant, the speed of light in vacuum and the orbital angular frequency, respectively. The Newtonian accurate amplitudes, $C^{(j)}_{+,\times}$ and $S^{(j)}_{+,\times}$, are written as power series in orbital eccentricity $e_{t}$ whose coefficients involve trigonometric functions of the two angles $\iota,\beta$ that specify the line of sight vector in a certain inertial frame. The derivation of these expressions is detailed in Ref. Yunes et al. 2009 and the required inputs are obtained by adapting a standard analytic approach to solve the classical Kepler equation in terms of the Bessel functionsColwell 1993.

With the help of Eqs. (1) and (2), we obtain interferometric strain for GWs from eccentric binaries as

where $\alpha_{j}={\rm sign}(\Gamma_{j})\sqrt{\Gamma_{j}^{2}+\Sigma_{j}^{2}}$ and $\phi_{j}=\tan^{-1}{\left(-\frac{\Sigma_{j}}{\Gamma_{j}}\right)}$. The two new functions, $\Gamma_{j}$ and $\Sigma_{j}$, are defined as $\Gamma_{j}=F_{+}\,C^{(j)}_{+}+F_{\times}\,C^{(j)}_{\times}$ and $\Sigma_{j}=F_{+}\,S^{(j)}_{+}+F_{\times}\,S^{(j)}_{\times}$, respectively as in Ref. Yunes et al. 2009. We impose the effects of GW emission on the above strain by specifying how $e_{t}$ and $\omega=2\,\pi\,F$, $F$ being the orbital frequency, vary in time. In Ref. Yunes et al. 2009, the temporal evolutions of $\omega$ and $e_{t}$ are governed by the following Newtonian (or quadrupolar) equations that were adapted from Refs. Peters and Mathews 1963; Peters 1964; Junker and Schaefer 1992.

It is customary to solve these two coupled differential equations numerically to obtain $\omega(t)$ and $e_{t}(t)$ and hence temporally evolving $h(t)$. Interestingly, earlier efforts provided certain analytic way for obtaining temporal evolution for $\omega(t)$ and $e_{t}(t)$ that mainly involves the usage of hypergeometric functions Moore et al. 2018; Mikóczi et al. 2012; Pierro et al. 2002; Pierro et al. 2001

However, it is possible to obtain analytic frequency domain counterpart of the above $h(t)$ as demonstrated in Ref. Królak et al. 1995; Yunes et al. 2009. This traditional approach involves the method of SPA, detailed in Ref. Bender and Orszag 1999, to compute analytically the Fourier Transform of $h(t)$. This was essentially demonstrated at the leading order in initial eccentricity $e_{0}$ in Ref. Królak et al. 1995 and later extended to ${\cal O}(e_{0}^{8})$ in Ref. Yunes et al. 2009. Following Refs. Królak et al. 1995; Yunes et al. 2009, we write

In the approach of stationary phase approximation, the crucial Fourier phase is given by

where $\tau$ stands for $F/\dot{F}$. Note that one needs to evaluate the above integrals at appropriate stationary points $t_{0}$, defined by $F(t_{0})=f/j$.

where $\chi$ is defined as $\omega/\omega_{0}=F/F_{0}$. We note that the above result was first obtained in Ref. Królak et al. 1995 which influenced them to introduce the idea of an asymptotic eccentric invariant . This relation allows us to write $\tau$ in terms of $\omega,\omega_{0}$ and $e_{0}$ as

It is now straightforward to compute analytically the indefinite integral for $\Psi_{j}$, namely

where $\phi_{c}$ and $t_{c}$ are the orbital phase at coalescence and the time of coalescence, respectively. Note that $\chi$ now stands for $f/f_{0}$ due to the use of the stationary phase condition. Additionally, we have re-scaled $F_{0}\rightarrow f_{0}/j$ to ensure that $e_{t}(f_{0})=e_{0}$ while employing the above expression for $e_{t}$, given by Eq. (8). Indeed, our expression is consistent with Eq. (4.28) of Ref. Yunes et al. 2009 that employs the chirp mass to characterize the binary. A number of extensions to the above result is available in the literature. In fact, Ref. Yunes et al. 2009 computed the higher order corrections to $e_{t}$ in terms of $e_{0}$ up to ${\cal O}({e_{0}^{7}})$ and extended $\Psi_{j}$ to ${\cal O}({e_{0}^{8}})$. Its PN extension, available in Ref. Tanay et al. 2016, provided 2PN corrections for $\Psi_{j}$ that incorporated eccentricity corrections, accurate to ${\cal O}({e_{0}^{6}})$ at every PN order while Ref. Moore et al. 2016 computed 3PN-accurate $\Psi_{j}$ that included leading order $e_{0}$ contributions.

A crucial ingredient to such PN extensions is the derivation of PN-accurate $e_{t}$ expression in terms of $e_{0},\chi$ and $x$. In what follows, we summarize the steps that are required to obtain 1PN-accurate expression for $e_{t}$ (see Ref. Tanay et al. 2016 for details). The starting point of such a derivation is the 1PN-accurate differential equations for $\omega$ and $e_{t}$, obtainable from Eqs. (3.12) in Ref. Tanay et al. 2016. With these inputs, it is fairly straightforward to obtain the following 1PN accurate expression for $d\omega/\omega$ that includes only the leading order $e_{t}$ contributions as

The fact that $\omega$ term appears only at the 1PN order allows us to use the earlier derived Newtonian accurate $\omega=\omega_{0}\,\left(e_{0}/e_{t}\right)^{18/19}$ relation to replace $\omega$ on the right hand side of the above equation. This leads to

where $x_{0}=\left(G\,m\,\omega_{0}/c^{3}\right)^{2/3}$. We can integrate this equation to obtain $\ln\omega-\ln\omega_{0}$ in terms of $e_{t},e_{0}$ and $\omega_{0}$. The exponential of the resulting expression and its bivariate expansion in terms of $x_{0}$ and $e_{t}$ result in

We invert the above equation to obtain $e_{t}$ in terms of $e_{0}$ and $x_{0}$ after invoking the Newtonian accurate relation $e_{t}=e_{0}\,\chi^{-19/18}$ to replace the $e_{t}$ terms associated with the $x_{0}$ term. This inversion and the associated bivariate expansion in terms of $e_{0}$ and $x_{0}$ require that $e_{0}\ll 1$ and $x_{0}\ll 1$. The resulting $e_{t}$ expression reads

To obtain $e_{t}$ as a bivariate expansion in terms of the regular PN parameter $x$ and $e_{0}$, we employ the fact that $x/x_{0}=\chi^{2/3}$ and this results in

We are now in a position to obtain 1PN-accurate $\Psi_{j}$ expression that includes ${\cal O}(e_{0}^{2})$ contributions both at the Newtonian and 1PN orders with the help of 1PN-accurate $\tau=\omega/\dot{\omega}$ expression that is accurate to ${\cal O}(e_{t}^{2})$ terms. A straightforward computation leads to the desired $\Psi_{j}$ expression which reads

II.2 Analytic PN-accurate amplitude corrected time domain eccentric GW templates

where $\phi=(1+k)\,l$, $\phi^{\prime}=k\,l$ and $k$ provides the rate of periastron advance per orbit Damour et al. 2004. Further, we let $c_{i}=\cos\iota$, $s_{i}=\sin\iota$, $c_{2\beta}=\cos 2\beta$ and $s_{2\beta}=\sin 2\beta$. Note that crucial ingredients to obtain above analytic expressions include developing approaches to solve PN-accurate Kepler equation and adapting them to derive PN-accurate relations to connect true and eccentric anomalies, detailed in Ref. Boetzel et al. 2017. A close inspection of the above two equations with Eqs. (3.3) and (3.4) of Ref. Yunes et al. 2009 reveals that the arguments of cosine and sine functions in above expressions involve $\phi^{\prime}=k\,l$ and its multiples in addition to the usual orbital phase $\phi$ and its multiples. These additional $\phi^{\prime}$ contributions are clearly due to the periastron advance. It turns out that these additional angular contributions are sufficient to provide the numerically inferred side bands in the power spectrum of eccentric binaries due to the presence of $k$ Tessmer and Gopakumar 2007. This is why we explicitly included $e_{t}^{4}$ contributions to the above $h_{\times,+}$ expressions as these contributions are required to reveal the underlying side band structure of waveforms due to the influence of periastron advance.

We re-write the above expressions for $h_{\times,+}^{0}$ in a more compact form to explicitly show how various harmonics are affected by the advance of periastron. The resulting expressions read

where we denoted the coefficient of $\cos(j\,\phi-(j\pm n)\phi^{\prime})$ harmonic at the quadrupolar (Newtonian) order for the $+$ polarization by $C_{+}^{j,\pm n}(0)$ while the coefficient of $\sin(j\,\phi-(j\pm n)\phi^{\prime})$ is indicated by $S_{+}^{j,\pm n}(0)$. We adopt a rather heavy notation as it is amenable to higher PN order contributions which will be tackled below. In this convention, we represent the coefficient of $\cos(j\,\phi-(j\pm n)\phi^{\prime})$ that appears in the 1PN contributions to $\times$ polarization state by $C_{\times}^{j,\pm n}(1)$. It should be obvious that $j$ stands for the harmonic variable while $n$ provides a measure of the shift that each harmonic experiences due to periastron advance. A close comparison of Eqs. (18) and (19) reveals that these coefficients are functions of $\iota,\beta$ and contain powers of $e_{t}$. Moreover, the arguments of cosine and sine functions clearly show that the eccentricity induced higher harmonics are not mere multiples of $\omega=N(1+k)$, where $N$ is the PN-accurate mean motion. Clearly, this is due to the presence of non-vanishing $\phi^{\prime}$ contributions due to periastron advance. Interestingly, the plus polarization state does provide harmonics which are integer multiples of $N$. It is not difficult to show that these Newtonian like terms arise from specific cosine functions with arguments $j\phi-j\phi^{\prime}$, as evident from Eqs. (19). Further, it is possible to show that these contributions arise from $e_{t}\,\cos u\,s_{i}^{2}/(1-e_{t}\,\cos u)$ contributions to $H_{+}^{0}$, given by Eq. (F2a) in Ref. Boetzel et al. 2017 and therefore not influenced by the periastron advance. Interestingly, similar conclusions were obtained in Ref. Tessmer and Gopakumar 2007.

With the above inputs, we write the time-domain GW detector response function for eccentric inspirals as

where the amplitudes of the cosine and sine functions are denoted by rather complicated symbols $\Gamma_{j,\pm n}^{(0)}$ and $\Sigma_{j,\pm n}^{(0)}$. The definition of $h(t)=\,F_{+}\,h_{+}(t)\,+\,F_{\times}\,h_{\times}(t)$ ensures that $\Gamma_{j,\pm n}^{(0)}=F_{+}\,C_{+}^{j,\pm n}(0)+F_{\times}\,C_{\times}^{j,\pm n}(0)$ while $\Sigma_{j,\pm n}^{(0)}=F_{+}\,S_{+}^{j,\pm n}(0)+F_{\times}\,S_{\times}^{j,\pm n}(0)$. We list in Appendix A, the lengthy expressions for these quantities in terms of $\iota,\beta$ and eccentricity contributions, accurate to ${\cal O}(e_{t}^{4})$. We display up to $\mathcal{O}(e_{t}^{4})$ contributions to demonstrate the full harmonic structure of the quadrupolar order GW polarization states. It turns out that $\Sigma_{j,0}^{(0)}$ contributions are zero by construction. This is mainly because the un-shifted harmonics only appear with the cosine terms, present in the $+$ polarization state. Invoking familiar trigonometric identities, we simplify the above equation and obtain

where we introduce two new multi-index symbols $\alpha_{j,\pm n}^{(0)}$ and $\bar{\phi}_{j,\pm n}^{(0)}$ to ensure that detector strain can be written in terms of only cosine functions. Influenced by Ref. Yunes et al. 2009, these symbols are defined as $\alpha_{j,\pm n}^{(0)}={\rm sign}\left(\Gamma_{j,\pm n}^{(0)}\right)\sqrt{\left(\Gamma_{j,\pm n}^{(0)}\right)^{2}+\left(\Sigma_{j,\pm n}^{(0)}\right)^{2}}$ and $\bar{\phi}_{j,\pm n}^{(0)}=\tan^{-1}\left(-\frac{\Sigma_{j,\pm n}^{(0)}}{\Gamma_{j,\pm n}^{(0)}}\right)$. We do not list explicit expressions for these quantities that are accurate to $\mathcal{O}(e_{t}^{4})$ in eccentricity corrections as they can be easily obtained from our Eqs. (59) and (60).

A close inspection of above equations reveal that they provide GW response function for compact binaries moving along precessing eccentric orbits. To obtain temporally evolving $h(t)$ associated with compact binaries inspiraling along precessing eccentric orbits, we need to specify how $\phi,\phi^{\prime},\omega$ and $e_{t}$ vary in time due to GW emission. We adapt the phasing formalism, detailed in Refs. Damour et al. 2004; Tanay et al. 2016, to provide differential equations for these variables. And, for the time being, we will concentrate on the secular evolution of these variables. In other words, we will neglect GW induced quasi-periodic variations to orbital elements and angles, detailed in Ref. Damour et al. 2004. The 3PN-accurate secular evolution to $\phi$ and $\phi^{\prime}$ in the modified harmonic gauge that are accurate to ${\cal O}(e_{t}^{6})$ are given by

The explicit $1.5,\,2,\,2.5$ and $3$PN order contributions to $d\omega/dt$ and $de_{t}/dt$ that incorporates all the $\mathcal{O}(e_{t}^{6})$ corrections are provided in the Appendix B. The differential equations for $d\omega/dt$ and $de_{t}/dt$ are extracted from expressions, available in Refs. Arun et al. 2009; Klein et al. 2018 and are in the modified harmonic gauge. These papers provided above 3PN accurate expressions as sum of certain ‘instantaneous’ and ‘tail’ contributions

The 3PN-accurate instantaneous contributions that depend only on the binary dynamics at the usual retarded time while the hereditary contributions are sensitive to the binary dynamics at all epochs prior to the usual retarded time Blanchet et al. 1995. The instantaneous contributions to $d\omega/dt$ are extracted from Eqs. (6.14),(6.15a),(6.15),(C6) and (C7) of Ref. Arun et al. 2009 while for $de_{t}/dt$ such contributions originate from Eqs. (6.16),(6.19a),(6.19b),(C10) and (C11) in Ref. Arun et al. 2009. It should be obvious that we have Taylor expanded these equations around $e_{t}=0$ to obtain eccentricity contributions accurate to ${\cal O}(e_{t}^{6})$. The hereditary contributions to $d\omega/dt$ and $de_{t}/dt$ are adapted from Eqs. (6.24c) and (6.26) of Ref. Arun et al. 2009 and they depend on a number of eccentricity enhancement functions. We employ such enhancement functions provided in Ref. Klein et al. 2018 for our computations. We now have all the inputs to obtain the restricted time-domain $h(t)$ to model GWs from non-spinning compact binaries inspiraling along precessing moderately eccentric orbits. To obtain such time domain templates, we numerically solve the above listed differential equations for $\omega,e_{t},\phi$ and $\phi^{\prime}$ and impose their temporal evolution in the quadrupolar order GW response function, given by Eq. (II.2). We now move onto describe how we extend the quadrupolar order GW response function.

It should be obvious that we require a prescription to compute analytically PN-accurate amplitude corrected GW polarization states to improve the above listed quadrupolar order GW response function. Therefore, we adapt 1PN-accurate amplitude corrected and fully analytic expressions for $h_{\times,+}$, available in Ref. Boetzel et al. 2017, to compute GW response function for eccentric inspirals that incorporates PN contributions even to its amplitudes. We list below certain ingredients that will be crucial to write down analytic $h(t)$ that incorporates 1PN-accurate amplitude corrections to $h_{\times,+}$ while consistently keeping eccentricity contributions up to ${\cal O}(e_{t}^{6})$. We begin by displaying Eqs. (44) and (45) of Ref. Boetzel et al. 2017 as a single sum which reads

Various PN order amplitude contributions take the following form

where $\delta=(m_{1}-m_{2})/(m_{1}+m_{2})$ and we let $m_{1}$ to be the heavier of the two binary components. We do not list explicitly very lengthy expressions for these amplitudes. However, they can be easily extracted from the attached Mathematica notebook. The derivation of above lengthy expressions include developing analytic approaches to solve PN-accurate Kepler equation and PN-accurate relations connecting true and eccentric anomalies, detailed in Ref. Boetzel et al. 2017. Indeed, we have verified that these expressions reduce to their circular counterparts, provided in Ref. Blanchet et al. 1996.

The associated GW detector strain for eccentric binaries is given by

A further simplification is possible which requires, as expected, additional multi-index functions

A cursory look at the above equation may give the impression that the summation indices in various sums are terminated in an arbitrary manner. Interestingly, we find a possible way to predict the maximum value that j index can take in each of the above summations. This is related to the argument of $\phi^{\prime}$ in each of these cosine series. We infer that the argument of $\phi^{\prime}$ can take a maximum value of six as we are restricting eccentricity contributions to sixth order in $e_{t}$. This ensures that $j$ index can take maximum values of $8,6$ and $4$ at the Newtonian order in the above expression. In other words, $j_{\rm max}$ in the above expression is given such that $j_{\rm max}\pm n=6$ where $\pm n$ value arises from the the argument of $\phi^{\prime}$ variable in various summations. It is easy to see that the above relation holds true even at 0.5 and 1PN orders and it provides a natural check on the structure of these higher order PN contributions to $h(t)$.

We begin by listing the expanded version of our quadrupolar order $h(t)$, namely Eq. (II.2) with $\mathcal{O}(e_{t}^{4})$ eccentricity contributions as

Clearly, we see three distinct square brackets that contain three cosine functions with explicitly time dependent arguments, namely $j\phi-(j-2)\phi^{\prime}$, $j\phi-j\phi^{\prime}$ and $j\phi-(j+2)\phi^{\prime}$. Note that $\alpha_{j,\pm n}^{(0)}$ and $\bar{\phi}_{j,\pm n}^{(0)}$ experience implicit temporal evolution due to the GW emission induced variations to $\omega$ and $e_{t}$. The main reason for displaying the above equation is to show explicitly how the periastron advance, defined by $\phi^{\prime}$, influences the harmonic structure of $h(t)$ in comparison with Eq. (4.21) of Ref. Yunes et al. 2009 or our Eq. (3).

where $l>0$ and as expected $S(t)$ should be a product of slowly varying amplitude $s(t)$ and a rapidly varying cosine function with argument $l\phi(t)$. Due to the virtue of Riemann-Lebesgue lemma, as noted in Ref. Bender and Orszag 1999, the Fourier transform of $S(t)$ becomes

It is not difficult to gather that the argument of the exponential function vanishes at the stationary point $t_{0}$ such that $l\dot{\phi}(t_{0})=2\pi f$. This allows us to invoke the approach of SPA to obtain the asymptotic behaviour of $S_{f}(f)$ by the following expression:

Note that $F(t)=\dot{\phi}(t)/2\pi$ and therefore its value at the stationary point should be $F(t_{0})=f/l$. Interestingly, a rather identical computation can be done to obtain the Fourier transform of a similar sinusoidal time series to be $i\,S_{f}(f)$.

To make operational the above expression for $S_{f}(f)$, we require an explicit expression for the above defined Fourier phase at the stationary point $t_{0}$, namely

This is done by defining $\tau=F/\dot{F}$ such that $\phi(F)$ and $t(F)$ become

where $\phi_{c}$ and $t_{c}$ are the orbital phase and time at coalescence. In the present context, $\tau$ is defined using our 3PN-accurate expression for $\dot{\omega}$ given by Eq. (25). Additionally, we require 3PN-accurate $e_{t}(\omega,\omega_{0},e_{0})$ expression, namely 3PN extension of Eq. (16), for computing these integrals analytically. The expression for $\Psi[F(t_{0})]$ obtained using Eq. (38) and (39) in (37) may be written as

where $\dot{\phi}=N(1+k)$ and this by definition is $\omega$. The treatment of $\dot{\phi}^{\prime}$ requires PN approximation as $\dot{\phi}^{\prime}$ equals $k\,N$ ( this is because $\phi^{\prime}=k\,l$). We need to express $k\,N$ in terms of $\omega$ and this leads to $\dot{\phi}^{\prime}=\omega\,k/(1+k)$ as $\omega=N(1+k)$. For computing Fourier phase analytically, we express $\dot{\phi}^{\prime}$ as $\omega\,k^{(6)}_{(3)}$, where $k^{(6)}_{(3)}$ stands for the 3PN-accurate expression for $k/(1+k)$ that incorporates $e_{t}$ contributions accurate to ${\cal O}(e_{t}^{6})$. The resulting expression reads

With the help of these inputs, the stationary points $t^{\pm n}$, where $\dot{\Psi}^{\pm n}(t^{\pm n})$ vanish, are given by

In other words, the stationary phase condition is given by

Rewriting $\Psi^{\pm n}(t)\coloneqq-2\pi ft+j\phi-(j\pm n)\phi^{\prime}$ using relation between $\phi^{\prime}$ and $\phi$ $\left(\phi^{\prime}=k^{(6)}_{(3)}\phi\right)$ gives $\Psi^{\pm n}(t)\coloneqq-2\pi ft+\left(j-(j\pm n)k^{(6)}_{(3)}\right)\phi$. We are now in a position to obtain analytic PN-accurate expressions for Fourier phases, associated with these stationary points. With Eq. (38) and (39), our Eq. (40) becomes

Note that $n$ takes values $0$ and $2$ as we are dealing with quadrupolar order GW response function given by Eq.(III.1). However, $n$ varies from $0$ to $4$ if the underlying GW response function contains 1PN-accurate amplitude corrections that include at each PN order eccentricity corrections accurate to ${\cal O}(e_{t}^{6})$. Further, we do not display here 3PN-accurate expression for $\tau$ that includes the leading order $e_{t}$ corrections, listed as Eqs. (6.7a) and (6.7b) in Ref. Moore et al. 2016. However, we do list below the explicit 3PN-accurate $\Psi_{j}^{\pm n}[F(t^{\pm n})]$ that incorporates leading order $e_{0}$ contributions at each PN order:

We now employ fully the final result of SPA, namely Eq. (36), to compute the Fourier transform of Eq. (III.1). This gives us

where the Fourier amplitudes $\xi_{j,\pm n}^{(0)}$ are now given by

In the above expression, the Fourier amplitudes are given by

The explicit expressions for $e_{t}$ and $\Psi^{n}_{j}(f)$ that incorporate the next to leading order $e_{0}$ corrections at each PN order, as noted earlier, are listed in the Appendix C.

We move on to contrast our approach with other attempts in the literature. The Sec. VI of Ref. Yunes et al. 2009 indeed sketched a road map to include PN corrections to their Newtonian waveform family. This road map included a suggestion to incorporate the effect of periastron advance into their quadrupolar order GW polarization states, influenced by Ref. Damour et al. 2004. Their suggestion involves splitting the orbital phase evolution into two parts where one part remains linear in the mean anomaly $l$ while the other part is periodic in $l$. These considerations influenced them to re-write our Eq. (3) essentially to be

where $k^{(6)}_{(1)}$ stands for the 1PN accurate expression for $k$, given by $3\,x/(1-e_{t}^{2})$, expanded to the sixth order in $e_{t}$ (see our Eq.(III.1)). It is not difficult to see that the associated SPA based Fourier phase takes the following form:

III.2 Preliminary GW data analysis implications

We employ the familiar match computations to probe basic GW data analysis implications of our PN-accurate inspiral templates. Following Ref. Damour et al. 1998, the match ${\cal M}(h_{s},h_{t})$ between members of two waveform classes, namely signal $h_{s}$ and template $h_{t}$, is computed by maximizing a certain overlap integral $\mathcal{O}(h_{s},h_{t})$ with respect to the kinematic variables of the template waveform. In other words,

where $t_{0}$ and $\phi_{0}$ are the detector arrival time and the associated arrival phase of our template. The overlap integral involves the interferometer-specific normalized inner product between members of $h_{s}$ and $h_{t}$ families; it reads

We require additional steps to operationalize our inspiral templates while performing the ${\cal M}$ computations. Clearly, these waveform families should only be implemented within the physically allowed frequency intervals. This is to ensure that the many higher harmonics present in these waveform families do not cross the above listed upper frequency limit. Influenced by Ref. (Yunes et al. 2009), we invoke the Unit Step function ($\Theta$) to operationalize our inspiral templates. This step function allows us to appropriately terminate the waveform as $\Theta(y)=1$ for $y\geq 0$ and zero otherwise. The structure of our quadrupolar amplitude inspiral family, given by Eqs. (III.1), compels us to invoke $\Theta$ functions such that

We qualify the implications of our match estimates on GW data analysis by considering the threshold $\mathcal{M}(h_{s},h_{t})\geq 0.97$, denoted in the presentation of results in Figs. 1, 2 and 3 by solid black lines. This limit corresponds to a loss of less than $10\%$ of all signals in the matched filter searches. In regions of parameter space where the computed matches are high, i.e., $\mathcal{M}\geq 0.97$, waveform models are generally considered both effectual templates for the detection of fiducial GW signals and reasonably faithful in the estimation of GW source parameters Damour et al. 1998. However, even if $\mathcal{M}$ larger than $0.97$, certain errors in the model waveform (due to unmodeled effects of, e.g., eccentricity) may become distinguishable from noise at high signal-to-noise ratio (SNR) and can affect the accuracy of source parameter estimation. Negligible systematic errors in parameter estimation – despite differences between the true signal waveform and the template model – can be guaranteed only if $(h_{s}-h_{t},h_{s}-h_{t})<1$, the so-called indistinguishability criterion Creighton and Anderson 2011. In other words, such systematic errors in the estimated source parameters may become significant when the mismatch $1-\mathcal{M}_{c}\geq 1/{\rm SNR}^{2}$ and clearly depend on the amplitude of the signal. In the following analysis, we let the signal-to-noise ratio of our fiducial GW signals be SNR$\,=30$ (corresponding to the SNR of the binary neutron star inspiral GW170817) and probe the distinguishability of certain effects in our model waveforms for inspiraling eccentric binaries. In the inset plots of Figs. 1 and 2, we zoom into those regions of parameter space where we can expect waveform uncertainties to become indistinguishable from noise for SNR$\,=30$; the corresponding distinguishable limit $\mathcal{M}_{c}$ is represented by the dashed black lines.

IV Conclusions

We have provided fully analytic PN-accurate Fourier domain gravitational waveforms for compact binaries inspiraling along precessing moderately eccentric orbits. Our inspiral approximant contains 1PN-accurate amplitude corrections and its Fourier phase incorporates the effects of 3PN-accurate periastron advance and GW emission. Additionally, the eccentricity effects are accurate to sixth order in $e_{0}$ at each PN order. We infer from our analytic waveform expression that the orbital eccentricity induced higher harmonics are no longer integer multiples of orbital frequency due to the influence of periastron advance. This substantiates and extends what is detailed in Ref. Boetzel et al. 2017 for compact binaries inspiraling along PN-accurate precessing eccentric orbits. Preliminary GW data analysis implications of our waveforms are probed with the help of the usual match computations.

In what follows, we provide a step-by-step summary of our effort.

We start from our Eqs. (18) and (19) that provide quadrupolar order GW polarization states from compact binaries in PN-accurate eccentric orbits as a sum over various harmonics.

With above inputs, we compute the time domain GW detector response function and express it as a summation of several cosine functions whose arguments are sum of integer multiples of $\phi$ and $\phi^{\prime}$ associated with the orbital and periastron motions. Amplitudes of these functions are expressed in terms of $\omega$, $e_{t}$ and the angles that specify the antenna patterns $F_{\times},F_{+}$ and the direction of the orbital angular momentum vector. The quadrupolar version of $h(t)$ that explicitly incorporates the next to leading order $e_{t}$ corrections is given by Eq. (II.2) and associated expressions like Eqs. (59) and (60). Its 1PN extension is symbolically provided by Eq. (II.2) and the accompanying Mathematica file provide the explicit expressions for various PN coefficients while incorporating ${\cal O}(e_{t}^{6})$ corrections.

We also provide a prescription to obtain temporally evolving $h(t)$ for compact binaries inspiraling due to 3PN accurate GW emission along precessing 3PN accurate orbits of moderate eccentricities. This involves imposing temporal evolution for $\omega,e_{t},\phi^{\prime}$ and $\phi$ with the help of PN accurate differential equations. The relative 3PN accurate equations for $\omega$ and $e_{t}$ are due to the emission of GWs, as evident from our Eqs. (25) and (26). The conservative 3PN accurate differential equation for $\phi^{\prime}$ arises essentially due to periastron advance as evident from Eq. (24). The differential equation for $\phi$ is kinematical in nature as $d\phi/dt\equiv\omega$.

A number of extensions are possible. Influenced by Refs. Königsdörffer and Gopakumar 2005; Kidder 1995, we are incorporating the effects of leading order aligned spin-orbit and spin-spin interactions into these waveforms. It will be interesting to explore data analysis implications of our present waveforms. A possible avenue is to explore the astrophysical implications of using PN-accurate periastron advance contributions that depend both on $m$ and $\eta$, influenced by Refs. Mikóczi et al. 2012; Seto 2001. There are on-going efforts to construct analytic IMR templates to model eccentric compact binary coalescence Huerta et al. 2017; Huerta et al. 2018. The present waveform family will be relevant to construct IMR templates for moderately eccentric compact binary mergers which can be used to extract orbital eccentricity and periastron advance as done in Ref Mroué et al. 2010. Efforts are on-going to obtain various constructs, using elements of our post-circular Fourier domain approximant, that should allow us to make comparisons with brand new PN-accurate frequency domain waveform family, developed in Refs. Moore et al. 2018; Moore and Yunes 2019 for moderate eccentricities.

V Acknowledgements

We thank Yannick Boetzel for helpful discussions, suggestions and providing us the lengthy eccentricity enhancement functions. We are grateful to Marc Favata and Blake Moore for their helpful comments. M. H. acknowledges support from Swiss National Science Foundation (SNSF) grant IZCOZ0_177057. We have used software packages from PyCBC Nitz et al. 2019 and Matplotlib Hunter 2007 to compute and plot match estimates.

Appendix A Γj,±n(0)\Gamma^{(0)}_{j,\pm n} and Σj,±n(0)\Sigma^{(0)}_{j,\pm n} coefficients

We list the $\Gamma^{(0)}_{j,\pm n}$ and $\Sigma^{(0)}_{j,\pm n}$ coefficients appearing in Eq.(21). The relevant $\Gamma_{j,\pm n}^{(0)}$ expressions read

The $\Sigma_{j,\pm n}^{(0)}$ counterparts of above expressions read

Appendix B 3PN accurate d​ωd​t\frac{d\omega}{dt} and d​etd​t\frac{de_{t}}{dt}

We give here the 3PN accurate expressions for temporal evolution of $\omega$ and $e_{t}$ for obtaining $h(t)$ associated with compact binaries inspiraling along precessing eccentric orbits. 1PN accurate $\frac{d\omega}{dt}$ and $\frac{de_{t}}{dt}$ with $\mathcal{O}(e_{t}^{6})$ eccentricity corrections are given by Eq.(25) and Eq.(26) respectively. The 1.5PN - 3PN contributions to $\frac{d\omega}{dt}$ appearing in Eq. 25 with $\mathcal{O}(e_{t}^{6})$ corrections are,

where $\gamma$ stands for the Euler-Mascheroni constant. The 1.5PN - 3PN contributions to $\frac{de_{t}}{dt}$ appearing in Eq. 26 with $\mathcal{O}(e_{t}^{6})$ corrections are,

Appendix C 3PN accurate analytic expressions for ete_{t} and Ψj±n\Psi_{j}^{\pm n}

We display explicit expressions for 3PN-accurate $e_{t}$ and Fourier phases that incorporate next to leading order $e_{0}$ corrections at each PN order. These expressions along with Eqs. (III.1), (48), (31), (59) and (60) are required to make operational the fully analytic frequency domain quadrupolar order GW response function for eccentric inspirals that includes ${\cal O}(e_{0}^{4})$ corrections at every PN order. We begin by listing explicit expression for the 3PN accurate $e_{t}$ in terms of $e_{0},\chi$ and $x$. The underlying computation is detailed in Ref. Tanay et al. 2016 and requires 3PN-accurate expressions for $\dot{\omega}$ and $\dot{e}_{t}$, given by Eqs. (25) and (26). The fully 3PN accurate $e_{t}$ expression that accounts for all the ${\cal O}(e_{0}^{3})$ contributions read

The coefficients $\mathcal{E}_{m}$ with next to leading order eccentricity corrections ${\cal O}{\left(e_{0}^{3}\right)}$ at each PN order can be listed as,

Due to the lengthy nature of 3PN order terms in $e_{t}$, we split it in two parts as

The explicit form of these two contributions are

We have pursued careful checking of our results with what is available in Ref. Tanay et al. 2016 and observed a slight typo in the $\mathcal{O}(e_{0}^{5})$ contributions for the $e_{t}$ expression (Eq. (A6e)) of Ref. Tanay et al. 2016. The $\eta$ independent term present in the coefficient of $\chi^{-119/18}$ should be $16952610560003855/162260186038272$ instead of $16633441088056655/162260186038272$. Note that the above $e_{t}$ expression is required while computing the Fourier amplitudes $\xi_{j}$. Additionally, it is a crucial ingredient while computing analytic expression for our Fourier phases $\Psi_{j}$. It should be obvious that its frequency dependence is encapsulated in $\chi=F/F_{0}$ and the PN expansion parameter $x=\left(G\,m\,2\,\pi\,F/c^{3}\right)^{2/3}$.

Various PN coefficients $\mathcal{P}_{m}$ with next to leading order eccentricity contributions are given by

For the ease of presentation, we split the 3PN contributions to $\Psi_{j}^{n}$ in to three parts

Various contributions to $\mathcal{P}_{6}$ are given by,

References