Accuracy of Estimating Highly Eccentric Binary Black Hole Parameters with Gravitational-wave Detections

The Advanced Laser Interferometer Gravitational Wave Observatory¹ (aLIGO) detectors (Aasi et al. 2015) and Advanced Virgo² (AdV; Acernese et al. 2015) have made the first six detections of gravitational waves (GWs) from approximately circular inspiraling binaries (Abbott et al. 2016c, 2016d, 2017a, 2017b, 2017c, 2017d), and opened a new window through which to observe the universe. These advanced GW detectors, together with upcoming instruments KAGRA³ (Somiya 2012) and LIGO-India⁴ (Iyer et al. 2011; Abbott et al. 2016f), are expected to continue to detect GW sources in the upcoming years (Abbott et al. 2016f). Orbital eccentricity was ignored in the analysis of the detected GW sources, but a preliminary upper limit was claimed to be e ≲ 0.1 at 10 Hz (Abbott et al. 2016a, 2016e, 2017a, 2017e). In this paper we estimate the future potential of the aLIGO–AdV–KAGRA network of advanced GW detectors to measure the orbital eccentricity and other physical parameters of initially highly eccentric sources.

Initially highly eccentric black hole (BH) binaries are inspiraling systems that have orbital eccentricities beyond ${e}_{0}\geqslant 0.9$ when their peak GW frequency (Wen 2003) enters the sensitive frequency band of advanced Earth-based GW detectors. The orbital eccentricity decreases in the inspiral phase from this value until the last stable orbit (LSO; Peters 1964). Such systems can form in multiple ways, including single–single encounters due to GW emission (Kocsis et al. 2006b; O'Leary et al. 2009; Gondán et al. 2017) in dense, high-velocity-dispersion environments; dynamical multibody interactions (Gültekin et al. 2006; O'Leary et al. 2006; Kushnir et al. 2013; Amaro-Seoane & Chen 2016; Antonini & Rasio 2016; Rodriguez et al. 2017); and the secular Kozai–Lidov mechanism (Wen 2003; Thompson 2011; Aarseth 2012; Antonini & Perets 2012; Antognini et al. 2014; Antonini et al. 2014, 2016; Breivik et al. 2016; Rodriguez et al. 2016a, 2016b; VanLandingham et al. 2016; Hoang et al. 2017; Petrovich & Antonini 2017; Silsbee & Tremaine 2017; Randall & Xianyu 2018) in hierarchical triples or binary–single interaction (Samsing et al. 2014, 2017; Samsing 2017; Samsing & Ramirez-Ruiz 2017). Eccentric BH binaries offer promising new detection candidates.

Previous parameter estimation studies of stellar-mass compact binaries have mostly focused on circular binaries (see Finn 1992; Finn & Chernoff 1993; Marković 1993; Cutler & Flanagan 1994; Jaranowski & Krolak 1994; Kokkotas et al. 1994; Królak et al. 1995; Poisson & Will 1995 for the first papers and Chatziioannou et al. 2014; Favata 2014; Mandel et al. 2014; O'Shaughnessy et al. 2014; Rodriguez et al. 2014; Berry et al. 2015; Canizares et al. 2015; Farr et al. 2016; Miller et al. 2015; Moore et al. 2016; Veitch et al. 2015; Lange et al. 2017; Vitale et al. 2017 for recent developments) due to their predicted high detection rates. The current detections constrain the merger rate density of BH–BH mergers in the universe to $12\mbox{--}213\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ (Abbott et al. 2017a), which corresponds to a detection rate between 400 and 7000 yr⁻¹ for a typical $2\,\mathrm{Gpc}$ detection range for aLIGO's design sensitivity. See Abadie et al. (2010) for a partial list of historical compact binary coalescence rate predictions and Dominik et al. (2013), Kinugawa et al. (2014), Abbott et al. (2016a, 2016b, 2016h, 2016i, 2016g), Belczynski et al. (2016), Rodriguez et al. (2016a), Bartos et al. (2017), McKernan et al. (2017), Hoang et al. (2017), Stone et al. (2017), and references therein for recent rate estimates.

However, several theoretical studies have shown that the detection rates of highly eccentric BH binaries may be non-negligible. For sources formed by GW emission in galactic nuclei (GNs), the expected aLIGO detection rate at design sensitivity may be higher than ≈100 yr⁻¹ if the BH mass function extends to masses above 25 ${M}_{\odot }$ (O'Leary et al. 2009; Kocsis & Levin 2012). Recently, such heavy BHs have been observed in several LIGO/VIRGO detections (Abbott et al. 2016b, 2017a, 2017c). In addition, the expected merger rate densities in the Kozai–Lidov channel are $1\mbox{--}1.5\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ for BH binaries forming in nuclear star clusters without supermassive BHs (SMBHs) through multibody interactions (Antonini & Rasio 2016) and $0.14\mbox{--}6.1\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ in isolated triple systems (Silsbee & Tremaine 2017). A merger rate density of order $1\mbox{--}5\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ is expected for BH binaries forming via the Kozai–Lidov mechanism in globular clusters (Antonini et al. 2014, 2016; Rodriguez et al. 2016a) and in GNs (Antonini & Perets 2012; Hoang et al. 2017), and non-spherical nuclear star clusters may produce BH binary merger rates of up to $15\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ (Petrovich & Antonini 2017). Smaller size GNs with intermediate-mass BHs may produce higher rates (VanLandingham et al. 2016). Binary–single gravitational interactions may greatly increase the rates (Samsing et al. 2014, 2017; Samsing & Ramirez-Ruiz 2017; Samsing 2017). In a companion paper (Gondán et al. 2017), we have shown that GW capture sources in GNs, which appear to be circular to within e < 0.2 near the LSO, may be highly eccentric at the beginning of the detected waveform at 10 Hz, and that heavier BH binaries are expected to be systematically more eccentric in this channel. The ongoing development of detectors toward their design sensitivity at low frequencies may open the possibility of detecting eccentricity in such systems.

In this paper, we determine the expected accuracy with which a network of ground-based interferometric GW detectors may determine the physical parameters that describe highly eccentric BH binaries in comparison to circular sources. We investigate how signal-to-noise ratios (S/Ns) and parameter measurement errors depend on the initial orbital parameters, particularly the initial pericenter distance and eccentricity. We examine if it is possible to measure the initial binary parameters (initial eccentricity and pericenter distance) at formation for sources that form in the GW frequency band of the instrument.

Previous GW parameter estimation accuracy studies for eccentric waveforms were carried out for extreme mass ratio (EMRI) sources around SMBHs for LISA (Barack & Cutler 2004; Cornish & Key 2010; Porter & Sesana 2010; Mikóczi et al. 2012; Nishizawa et al. 2016) and for low-eccentricity stellar-mass compact binaries for Earth-based GW detector network (Sun et al. 2015). The premerger localization accuracy of eccentric neutron star (NS) binary systems was determined by Kyutoku & Seto (2014), and the source localization accuracy was investigated for low-eccentricity binaries by Ma et al. (2017).

The parameter space of an eccentric spinning binary waveform is generally very large, with 17 dimensions (Vecchio 2004; Cornish & Key 2010). Therefore, state-of-the-art methods such as Monte Carlo Markov Chain calculations (see O'Shaughnessy et al. 2014 and references therein) are numerically prohibitively expensive to explore the full range of source parameters for a large set of binaries. For Gaussian noise and a large S/N, the posterior distribution function of the measured parameters is generally well approximated by a multidimensional Gaussian, and the parameter measurement errors can be estimated accurately and very efficiently using the Fisher matrix method (Finn & Chernoff 1993; Cutler & Flanagan 1994; Cutler & Vallisneri 2007). Using this technique, we determine the physical parameters' measurement accuracy.

We restrict this first study to waveforms introduced by Moreno-Garrido et al. (1994, 1995), which account for part of the leading-order post-Newtonian (PN) correction, the GR pericenter precession (hereafter simply precession), and ignore other first-order PN and higher-order corrections including those due to spins. Future extensions of this work should include higher-order PN and merger waveforms (see Levin et al. 2011; Csizmadia et al. 2012; East et al. 2013 for waveform generators, Damour et al. 2004; Memmesheimer et al. 2004; Königsdörffer & Gopakumar 2005, 2006; Yunes et al. 2009; Tessmer & Schäfer 2010, 2011; Huerta et al. 2014; Mikóczi et al. 2015; Moore et al. 2016; Tanay et al. 2016; Boetzel et al. 2017; Cao & Han 2017; Hinderer & Babak 2017; Huerta et al. 2017a, 2017b; Loutrel & Yunes 2017 for analytic waveform models, and Hinder et al. 2008; East et al. 2012, 2015, 2016; Gold et al. 2012; Gold & Brügmann 2013; Paschalidis et al. 2015; Lewis et al. 2017 for waveforms of numerical relativity simulations of eccentric compact binary inspirals). In this paper, we focus on BH–BH binaries, but the method is also applicable to NS–NS and NS–BH binaries on highly eccentric orbits as long as tidal interactions and matter exchange among the components are negligible (see Gold et al. 2012; East et al. 2015, 2016; Radice et al. 2016 and references therein).⁵

Once a large number of GW sources is detected, the correlations between the orbital eccentricity, binary total mass, reduced mass, and spins may be distinct among the different astrophysical mechanisms leading to BH mergers (O'Leary et al. 2009; Cholis et al. 2016; Rodriguez et al. 2016a, 2016b; Chatterjee et al. 2017; Gondán et al. 2017; Samsing & Ramirez-Ruiz 2017; Silsbee & Tremaine 2017; Kocsis et al. 2017). Therefore, detections of eccentric BH binaries have the potential to constrain GW source populations.

However, detecting eccentric sources and recovering their physical parameters are very challenging. So far, three search methods have been developed to find the signals of stellar-mass eccentric BH binaries in data streams of GW detectors (Tai et al. 2014; Coughlin et al. 2015; Tiwari et al. 2016). All three methods achieve substantially better sensitivity for eccentric BH binary signals than existing localized burst searches or chirp-like template-based search methods. Once a source is detected, different algorithms are used to recover its physical parameters. For compact binary coalescences, BAYESTAR (Singer & Price 2016) is an online fast sky-localization algorithm that produces probability sky maps, LALInference (Veitch et al. 2015) is an offline full parameter estimation algorithm, and gstlal (Cannon et al. 2012; Privitera et al. 2014) is a low-latency binary BH parameter estimation algorithm. All three algorithms use waveform models of compact binaries on circular orbits. In addition, for short-duration GW "bursts" with poorly modeled or unknown waveforms, Coherent WaveBurst (Klimenko et al. 2016), BayesWave (Cornish & Littenberg 2015), and LALInferenceBurst (Veitch et al. 2015) pipelines produce reconstructed waveforms with minimal assumptions on the waveform morphology. The development of algorithms recovering the parameters of compact binaries on eccentric orbits are currently underway. These algorithms will play an important role in the astrophysical interpretation of eccentric sources.

This paper is organized as follows. In Section 2, we summarize the basic formulae describing the time-domain and frequency-domain eccentric waveform model. In Section 3, we outline the properties of the advanced detectors we use in the analysis. In Section 4, we describe the signal parameter measurement estimation method. In Section 5, we discuss which parameters of an eccentric binary can be measured through the binary's waveform. We present our main results in Section 6, and compare or results with previous papers. Finally, we summarize our conclusions in Section 7. Several details of our methodology are included in the appendices. In Appendix A, we consider the values of the source parameters in the circular limit. Next, in Appendix B, we introduce the geometric conventions we use to describe how the GWs interact with ground-based detectors. In Appendix C, we discuss the applicability of assumptions ignoring the Earth's rotation around its axis and the Earth's motion around the Sun. In Appendix D, we derive numerically effective formulae to reduce the computational cost of numerical calculations of the S/N and the Fisher matrix. In Appendix E, we present numerical comparisons to validate our codes for both precessing and non-precessing waveforms.

We use G = 1 = c units when referring to the initial orbital parameters and when determining the phases of waveforms. We work in the observer frame assuming a binary at cosmological redshift z. In this frame, all of the formulae have redshifted mass parameters ${m}_{z}=(1+z)m$.⁶

In this section, we summarize the basic formulae describing the time-domain (Section 2.1) and frequency-domain (Section 2.2) eccentric waveform models including precession in the leading quadrupole-order radiation approximation using the Fourier–Bessel decomposition. Note that we ignore the radiation of higher multipole orders, which are typically subdominant at least in cases where the initial pericenter distance is not close to a grazing or zoom–whirl configuration and the initial velocity is much less than the speed of light (Davis et al. 1972; Berti et al. 2010; Healy et al. 2016).

2.1. The Waveform in the Time Domain

We adopt the waveform model of Moreno-Garrido et al. (1994) and Moreno-Garrido et al. (1995), which describes the quadrupole waveform emitted by a spinless binary on a Keplerian orbit undergoing slow precession. For a fixed semimajor axis a and orbital eccentricity e, the two polarization states of a GW, h₊ and h_×, with component masses m_A and m_B, and at luminosity distance D_L, can be given in the observer's time domain as (Moreno-Garrido et al. 1995)

$$\begin{eqnarray}{h}_{+}(t) & = & -\displaystyle \frac{h\,{\sin }^{2}{\rm{\Theta }}}{2}\displaystyle \sum _{n=1}^{\infty }{A}_{n}\cos {{\rm{\Phi }}}_{n}(t)\\ & & +\displaystyle \frac{h(1+{\cos }^{2}{\rm{\Theta }})}{2}\displaystyle \sum _{n=1}^{\infty }({B}_{n}^{+}\cos {{\rm{\Phi }}}_{n}^{-}(t)-{B}_{n}^{-}\cos {{\rm{\Phi }}}_{n}^{+}(t)),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{h}_{\times }(t)=-h\,\cos {\rm{\Theta }}\sum _{n=1}^{\infty }({B}_{n}^{-}\sin {{\rm{\Phi }}}_{n}^{+}(t)+{B}_{n}^{+}\sin {{\rm{\Phi }}}_{n}^{-}(t)),\end{eqnarray} \tag{ 2 }$$

where Θ is the angle between the orbital plane and the line of sight to the observer, ${{\rm{\Phi }}}_{n}^{\pm }$ describes the orbital phase given below for the $n\mathrm{th}$ harmonic,

$$\begin{eqnarray}&&h=\displaystyle \frac{4\,{M}_{\mathrm{tot},z}\,{\mu }_{z}}{a\,{D}_{L}},\end{eqnarray} \tag{ 3 }$$

where ${M}_{\mathrm{tot},z}=({m}_{A}+{m}_{B})(1+z)$ is the redshifted total binary mass at cosmological redshift z, and ${\mu }_{z}=(1+z){m}_{A}{m}_{B}({m}_{A}+{{m}_{B})}^{-1}$ is the redshifted reduced mass. We can express the luminosity distance for a flat ΛCDM cosmology as a function of z as

$$\begin{eqnarray}&&{D}_{L}=\displaystyle \frac{(1+z)c}{{H}_{0}}{\int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{\sqrt{{{\rm{\Omega }}}_{{\rm{M}}}{(1+z^{\prime} )}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}},\end{eqnarray} \tag{ 4 }$$

where ${H}_{0}=68\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ is the Hubble constant, and ${{\rm{\Omega }}}_{{\rm{M}}}=0.304$ and Ω_Λ = 0.696 are the density parameters for matter and dark energy, respectively (Planck Collaboration et al. 2014a, 2014b).

The A_n and ${B}_{n}^{\pm }$ prefactors in Equations (1) and (2) are the linear combinations of the Bessel function of the first kind, J_n(x),

$$\begin{eqnarray}&&{A}_{n}={J}_{n}({ne}),\quad {B}_{n}^{\pm }=\displaystyle \frac{{S}_{n}\pm {C}_{n}}{2},\end{eqnarray} \tag{ 5 }$$

where e is the orbital eccentricity,

$$\begin{eqnarray}&&{S}_{n}=-\displaystyle \frac{{(1-{e}^{2})}^{1/2}}{e}\displaystyle \frac{2}{n}{J}_{n}^{{\prime} }({ne})+\displaystyle \frac{{(1-{e}^{2})}^{3/2}}{{e}^{2}}2{{nJ}}_{n}({ne}),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{C}_{n}=-\displaystyle \frac{2-{e}^{2}}{{e}^{2}}{J}_{n}({ne})+\displaystyle \frac{2(1-{e}^{2})}{e}{J}_{n}^{{\prime} }({ne}),\end{eqnarray} \tag{ 7 }$$

and ${J}_{n}^{{\prime} }({ne})$ is the first derivative of J_n(ne) with respect to e, which satisfies

$$\begin{eqnarray}&&{J}_{n}^{{\prime} }({ne})=\displaystyle \frac{n}{2}[{J}_{n-1}({ne})-{J}_{n+1}({ne})].\end{eqnarray} \tag{ 8 }$$

We will also need the second derivative of J_n(ne) with respect to e when calculating the Fisher matrix, thus we introduce $J{{\prime\prime} }_{n}({ne})$ as

$$\begin{eqnarray}&&J{{\prime\prime} }_{n}=\displaystyle \frac{{n}^{2}}{4}[{J}_{n-2}({ne})-2{J}_{n}({ne})+{J}_{n+2}({ne})]\quad \mathrm{if}\,n\geqslant 2,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&J{{\prime\prime} }_{1}=-\displaystyle \frac{{J}_{1}(e)}{2}-\displaystyle \frac{1}{4}[{J}_{1}(e)-{J}_{3}(e)].\end{eqnarray} \tag{ 10 }$$

The phase functions ${{\rm{\Phi }}}_{n}(t)$ and ${{\rm{\Phi }}}_{n}^{\pm }(t)$ in Equations (1) and (2) are

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{n}(t)={{\rm{\Phi }}}_{c}-2\pi n{\int }_{t}^{{t}_{c}}\nu (t^{\prime} ){dt}^{\prime} ,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{n}^{\pm }(t)={{\rm{\Phi }}}_{n}(t)\pm 2\gamma (t),\end{eqnarray} \tag{ 12 }$$

where the second term on the right-hand side is n times the mean anomaly expressed by the time integral of the redshifted Keplerian mean orbital frequency ν, ${{\rm{\Phi }}}_{c}$ is the phase extrapolated to coalescence time $t={t}_{c}$, and γ is the azimuthal angle of the pericenter relative to the x-axis of the coordinate system defined by the orbital plane. The redshifted Keplerian mean orbital frequency may be expressed with the dimensionless pericenter distance

$$\begin{eqnarray}&&{\rho }_{{\rm{p}}}=\displaystyle \frac{a(1-e)}{{M}_{\mathrm{tot},z}}\end{eqnarray} \tag{ 13 }$$

(where a is the semimajor axis in the observer frame) as

$$\begin{eqnarray}&&\nu (e,{\rho }_{{\rm{p}}})=\displaystyle \frac{{(1-e)}^{3/2}}{2\pi {\rho }_{{\rm{p}}}^{3/2}{M}_{\mathrm{tot},z}}.\end{eqnarray} \tag{ 14 }$$

For an inspiraling binary, both the eccentricity and the Keplerian orbital frequency evolve in time. Assuming quadrupole radiation and the adiabatic evolution of orbital parameters, the equations of the time evolution of e and ν, as seen at some cosmological redshift, can be given to leading order as (Peters 1964)

$$\begin{eqnarray}&&\dot{e}=-\displaystyle \frac{304}{15}\displaystyle \frac{e{{ \mathcal M }}_{z}^{5/3}{(2\pi \nu )}^{8/3}}{{(1-{e}^{2})}^{5/2}}\left(1+\displaystyle \frac{121}{304}{e}^{2}\right),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\dot{\nu }=\displaystyle \frac{48}{5\pi }\displaystyle \frac{{{ \mathcal M }}_{z}^{5/3}{(2\pi \nu )}^{11/3}}{{(1-{e}^{2})}^{7/2}}\left(1+\displaystyle \frac{73}{24}{e}^{2}+\displaystyle \frac{37}{96}{e}^{4}\right),\end{eqnarray} \tag{ 16 }$$

where ${{ \mathcal M }}_{z}={\mu }_{z}^{3/5}{M}_{\mathrm{tot},z}^{2/5}$ is the redshifted chirp mass, and the overdot denotes a redshifted time derivative, $\dot{x}\equiv {dx}/{dt}$. The fraction of the two equations

$$\begin{eqnarray}&&\displaystyle \frac{d\nu }{{de}}=-\displaystyle \frac{18}{19}\displaystyle \frac{\nu }{e}\displaystyle \frac{\left(1+\tfrac{73}{24}{e}^{2}+\tfrac{37}{96}{e}^{4}\right)}{(1-{e}^{2})\left(1+\tfrac{121}{304}{e}^{2}\right)}\end{eqnarray} \tag{ 17 }$$

may be integrated as (Peters 1964; Mikóczi et al. 2012)

$$\begin{eqnarray}&&\nu (e)={c}_{0}/H(e),\end{eqnarray} \tag{ 18 }$$

where we define

$$\begin{eqnarray}&&H(e)={e}^{18/19}{(1-{e}^{2})}^{-3/2}{\left(1+\displaystyle \frac{121}{304}{e}^{2}\right)}^{\tfrac{1305}{2299}},\end{eqnarray} \tag{ 19 }$$

and ${c}_{0}$ is an integration constant set by the initial condition $\nu ({e}_{0},{\rho }_{{\rm{p}}0})={\nu }_{0}$ or the conditions at the LSO, $\nu ({e}_{\mathrm{LSO}},{\rho }_{\mathrm{pLSO}})={\nu }_{\mathrm{LSO}}$ (see Equations (21) and (23) below). Equation (18) shows that the product ${c}_{0}=\nu (e)H(e)$ is conserved during the evolution. Similarly, it is straightforward to determine the evolution of the dimensionless pericenter distance,

$$\begin{eqnarray}{\rho }_{{\rm{p}}}(e) & = & \displaystyle \frac{{c}_{1}{e}^{12/19}}{{M}_{\mathrm{tot},z}^{2/3}(1+e)}{\left(1+\displaystyle \frac{121}{304}{e}^{2}\right)}^{\tfrac{870}{2299}}\\ & = & {\rho }_{{\rm{p}}0}\displaystyle \frac{{e}^{12/19}{(1+e)}^{-1}{[1+(121/304){e}^{2}]}^{\tfrac{870}{2299}}}{{e}_{0}^{12/19}{(1+{e}_{0})}^{-1}{[1+(121/304){e}_{0}^{2}]}^{\tfrac{870}{2299}}}\\ & = & {\rho }_{\mathrm{pLSO}}\displaystyle \frac{{e}^{12/19}{(1+e)}^{-1}{[1+(121/304){e}^{2}]}^{\tfrac{870}{2299}}}{{e}_{\mathrm{LSO}}^{12/19}{(1+{e}_{\mathrm{LSO}})}^{-1}{[1+(121/304){e}_{\mathrm{LSO}}^{2}]}^{\tfrac{870}{2299}}}\end{eqnarray} \tag{ 20 }$$

(Peters 1964), where ${c}_{1}={(2\pi {c}_{0})}^{-2/3}$, and in the second and third lines we expressed the evolution with the initial condition ${\rho }_{{\rm{p}}0}$ and e₀, or the "final condition" at the LSO, which satisfies

$$\begin{eqnarray}&&{\rho }_{\mathrm{pLSO}}={\rho }_{{\rm{p}}}({e}_{\mathrm{LSO}})=\displaystyle \frac{6+2{e}_{\mathrm{LSO}}}{1+{e}_{\mathrm{LSO}}}\end{eqnarray} \tag{ 21 }$$

in the leading-order approximation in the test-mass geodesic zero-spin limit (Cutler et al. 1994). This shows that the evolution may be parameterized with the single parameter e_LSO, or with the two parameters ${\rho }_{{\rm{p}}0}$ and e₀. Note that for any e, the orbital frequency depends only on the single parameter ${c}_{0}$, which is set uniquely by ${e}_{\mathrm{LSO}}$ and ${M}_{\mathrm{tot},z}$ as

$$\begin{eqnarray}&&{c}_{0}=\displaystyle \frac{1}{2\pi {M}_{\mathrm{tot},z}}\displaystyle \frac{{(1-{e}_{\mathrm{LSO}})}^{3/2}H({e}_{\mathrm{LSO}})}{{[{\rho }_{\mathrm{pLSO}}({e}_{\mathrm{LSO}})]}^{3/2}}.\end{eqnarray} \tag{ 22 }$$

We restrict our interest to the repeated burst (O'Leary et al. 2009; Kocsis & Levin 2012) and eccentric inspiral phases of the waveform model between $0\lt {e}_{\mathrm{LSO}}\leqslant e\leqslant {e}_{0}\leqslant 1$ and ignore the merger and ringdown phases in this analysis. The repeated burst phase starts when the binary is formed with initial eccentricity e₀ > 0.9 and initial dimensionless pericenter distance ${\rho }_{{\rm{p}}0}$, and the eccentric inspiral phase ends when the binary reaches the LSO with eccentricity e_LSO. Note that during the evolution, e and ρ_p both shrink strictly monotonically in time.

Let us also note for further use that the Keplerian redshifted orbital frequency at the end of the assumed eccentric inspiral waveform (i.e., at the LSO) is given by Equations (14) and (21) as

$$\begin{eqnarray}&&{\nu }_{\mathrm{LSO}}=\nu ({e}_{\mathrm{LSO}})=\displaystyle \frac{1}{2\pi {M}_{\mathrm{tot},z}}{\left(\displaystyle \frac{1-{e}_{\mathrm{LSO}}^{2}}{6+2{e}_{\mathrm{LSO}}}\right)}^{3/2}.\end{eqnarray} \tag{ 23 }$$

Precession leads to a time-dependent γ in Equation (12). Using the analysis in Mikóczi et al. (2012), we adopt pericenter precession from classical relativistic motion, and assume that the adiabatic evolution of the orbital parameters is governed by Equations (15) and (16). The angle of precession for a single eccentric orbit in the test-particle geodesic approximation around a Schwarzschild BH is

$$\begin{eqnarray}&&{\rm{\Delta }}\gamma =\displaystyle \frac{6\pi {M}_{\mathrm{tot}}}{a(1-{e}^{2})}.\end{eqnarray} \tag{ 24 }$$

Using an adiabatic approximation, we approximate the redshifted precession rate to be constant during the orbit with

$$\begin{eqnarray}&&\dot{\gamma }\approx \displaystyle \frac{{\rm{\Delta }}\gamma }{T}=\displaystyle \frac{3{(2\pi \nu )}^{5/3}{M}_{\mathrm{tot},z}^{2/3}}{1-{e}^{2}}.\end{eqnarray} \tag{ 25 }$$

The phase functions given by Equations (11) and (12) can be calculated from Equations (15) and (18) as⁷

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{n}(t)={{\rm{\Phi }}}_{c}+2\pi n{\int }_{0}^{e(t)}\displaystyle \frac{\nu (e^{\prime} )}{\dot{e}(\nu (e^{\prime} ),e^{\prime} )}{de}^{\prime} \end{eqnarray} \tag{ 26 }$$

(Cutler & Flanagan 1994). The phase functions, which arise due to precession, ${{\rm{\Phi }}}_{n}^{+}(t)$ and ${{\rm{\Phi }}}_{n}^{-}(t)$, follow from Equations (12), (18), (25), and (26),

$$\begin{eqnarray}{{\rm{\Phi }}}_{n}^{\pm }(t) & = & {{\rm{\Phi }}}_{n}(t)\pm 2{\gamma }_{c}\mp 2{\displaystyle \int }_{t}^{{t}_{c}}\dot{\gamma }(t^{\prime} ){dt}^{\prime} \\ & = & {{\rm{\Phi }}}_{c}\pm 2{\gamma }_{c}+{\displaystyle \int }_{0}^{e(t)}\displaystyle \frac{2\pi n\nu (e^{\prime} )\pm 2\dot{\gamma }(\nu (e^{\prime} ),e^{\prime} )}{\dot{e}(\nu (e^{\prime} ),e^{\prime} )}{de}^{\prime} ,\end{eqnarray} \tag{ 27 }$$

where ${\gamma }_{c}$ is the angle of periapsis extrapolated to coalescence. Note that ${\rho }_{{\rm{p}}0}$, t_c, ${\gamma }_{c}$, Φ_c, and ${e}_{\mathrm{LSO}}$ are free parameters of the waveform. Alternatively, we may use the corresponding initial values e₀, t₀, γ₀, Φ₀, and ${\rho }_{{\rm{p}}0}$.

2.2. The Waveform in the Frequency Domain

Since the expressions defining the S/N and the Fisher matrix are both given in Fourier space (Section 4), we construct the Fourier transforms of the waveform⁸

$$\begin{eqnarray}&&{\tilde{h}}_{+,\times }(f)={\int }_{-\infty }^{\infty }{h}_{+,\times }(t){e}^{2\pi {if}}{dt},\end{eqnarray} \tag{ 28 }$$

where ${h}_{+}(t)$ and ${h}_{\times }(t)$ are given in Equations (1) and (2) as an infinite sum over orbital harmonics n. In the stationary phase approximation, each frequency harmonic splits into a triplet due to precession (see Equation (25)), ${\boldsymbol{f}}\equiv ({f}_{n},{f}_{n}^{\pm })$ (Moreno-Garrido et al. 1995; Mikóczi et al. 2012), where

$$\begin{eqnarray}&&{f}_{n}=n\nu ,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{f}_{n}^{\pm }\,=\,n\nu \pm \displaystyle \frac{\dot{\gamma }}{\pi },\end{eqnarray} \tag{ 30 }$$

and the Fourier transform simplifies to

$$\begin{eqnarray}{\tilde{h}}_{+}({\boldsymbol{f}}) & = & -\displaystyle \frac{{h}_{0}}{4}{\sin }^{2}{\rm{\Theta }}\displaystyle \sum _{n=1}^{\infty }{A}_{n}{{\rm{\Lambda }}}_{n}{e}^{i({{\rm{\Psi }}}_{n}-\pi /4)}\\ & & -\displaystyle \frac{{h}_{0}}{4}(1+{\cos }^{2}{\rm{\Theta }})\displaystyle \sum _{n=1}^{\infty }{B}_{n}^{+}{{\rm{\Lambda }}}_{n}^{-}{e}^{i({{\rm{\Psi }}}_{n}^{-}-\pi /4)}\\ & & +\displaystyle \frac{{h}_{0}}{4}(1+{\cos }^{2}{\rm{\Theta }})\displaystyle \sum _{n=1}^{\infty }{B}_{n}^{-}{{\rm{\Lambda }}}_{n}^{+}{e}^{i({{\rm{\Psi }}}_{n}^{+}-\pi /4)},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}{\tilde{h}}_{\times }({\boldsymbol{f}}) & = & -\displaystyle \frac{{h}_{0}}{2}\,\cos {\rm{\Theta }}\displaystyle \sum _{n=1}^{\infty }{B}_{n}^{-}{{\rm{\Lambda }}}_{n}^{+}{e}^{i({{\rm{\Psi }}}_{n}^{+}+\pi /4)}\\ & & -\displaystyle \frac{{h}_{0}}{2}\,\cos {\rm{\Theta }}\displaystyle \sum _{n=1}^{\infty }{B}_{n}^{+}{{\rm{\Lambda }}}_{n}^{-}{e}^{i({{\rm{\Psi }}}_{n}^{-}+\pi /4)},\end{eqnarray} \tag{ 32 }$$

where

$$\begin{eqnarray}&&{h}_{0}=\displaystyle \frac{4{{ \mathcal M }}_{z}^{5/3}{(2\pi \nu )}^{2/3}}{{D}_{L}},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{n}=\displaystyle \frac{1}{\sqrt{| n\dot{\nu }| }},\quad {{\rm{\Lambda }}}_{n}^{\pm }=\displaystyle \frac{1}{\sqrt{| n\dot{\nu }\pm \ddot{\gamma }/\pi | }},\end{eqnarray} \tag{ 34 }$$

$\dot{\nu }$ and $\ddot{\gamma }$ are given by Equations (16), (17), and (25), and ${{ \mathcal M }}_{z}=(1+z){ \mathcal M }$. The ${{\rm{\Psi }}}_{n}$ and ${{\rm{\Psi }}}_{n}^{\pm }$ phases and their first (${\dot{{\rm{\Psi }}}}_{n},{\dot{{\rm{\Psi }}}}_{n}^{\pm }$) and second (${\ddot{{\rm{\Psi }}}}_{n},{\ddot{{\rm{\Psi }}}}_{n}^{\pm }$) derivatives with respect to redshifted time t are (Mikóczi et al. 2012)

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{n}(e,{f}_{n})=2\pi {f}_{n}{t}_{n}-{{\rm{\Phi }}}_{n},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{n}^{\pm }(e,{f}_{n}^{\pm })=2\pi {f}_{n}^{\pm }{t}_{n}^{\pm }-{{\rm{\Phi }}}_{n}^{\pm }.\end{eqnarray} \tag{ 36 }$$

Here, the $({t}_{n},{t}_{n}^{\pm })$ parameters of the stationary phase approximation specify the times at which the orbital frequency satisfies Equations (29) and (30) for a given $({f}_{n},{f}_{n}^{\pm })$; see Appendices A and B in Mikóczi et al. (2012) for details. In Equations (35) and (36), $({{\rm{\Phi }}}_{n},{{\rm{\Phi }}}_{n}^{\pm })$ are to be substituted from Equations (26) and (27). We eliminate $\nu (t)$ and e(t) for $({f}_{n},{f}_{n}^{\pm })$ using Equations (18), (22), and (25) together with Equations (29) and (30) to obtain the frequency-domain waveform.⁹ The result depends on constant parameters t_c, Φ_c, γ_c, e_LSO, ${{ \mathcal M }}_{z}$, and ${M}_{\mathrm{tot},z}$. Further, we note that if the precessing eccentric BH binary forms with ${\rho }_{{\rm{p}}0}$ and e₀, then the frequency-domain waveform is truncated at the corresponding minimum frequency $({f}_{n,\min }$, ${f}_{n,\min }^{\pm })=(n{\nu }_{0},n{\nu }_{0}\pm {\dot{\gamma }}_{0}/\pi )$. Furthermore, the waveform model becomes invalid after reaching the LSO (with ${\rho }_{{\rm{p}}\mathrm{LSO}}$ and e_LSO), which corresponds to a maximum frequency for each harmonic $({f}_{n,\max }$, ${f}_{n,\max }^{\pm })=(n{\nu }_{\mathrm{LSO}},n{\nu }_{\mathrm{LSO}}\pm {\dot{\gamma }}_{\mathrm{LSO}}/\pi )$, where this model is applicable. If we truncate the waveform at these maximum frequencies, this respectively introduces an explicit $({\rho }_{{\rm{p}}0},{e}_{0})$ and $({\rho }_{{\rm{p}}\mathrm{LSO}},{e}_{\mathrm{LSO}})$ parameter dependence in the waveform model. This is shown in Equation (104) in Appendix D. Examples of the frequency-domain waveforms are shown in Kocsis & Levin (2012).

In principle, the number of spectral harmonics of an eccentric binary system is infinite. Note, however, that a large fraction of the signal power is accumulated in a finite number of harmonics. Therefore, in order to reduce the necessary computation time, we truncate n at ${n}_{\max }({e}_{0})$ (O'Leary et al. 2009; Mikóczi et al. 2012),

$$\begin{eqnarray}&&{n}_{\max }({e}_{0})=\left\{5\displaystyle \frac{{(1+{e}_{0})}^{1/2}}{{(1-{e}_{0})}^{3/2}}\right\},\end{eqnarray} \tag{ 37 }$$

which accounts for 99% of the signal power (Turner 1977). Here, the bracket {} denotes the floor function. In Appendix D, we discuss other technical details to optimize the calculation of the S/N and the Fisher matrix.

To test our calculations, we examine the limiting cases of no precession $(\dot{\gamma }\to 0)$ and circular orbits ($e\to 0$), respectively. In Appendix D, we discuss numerical and analytical tricks to optimize the calculation and discuss the results for the precessing (Prec) and precession-free (NoPrec) waveform model (i.e., $\dot{\gamma }\equiv 0$).

Here we summarize the GW detectors and the assumed properties of the detector noise in our analysis.

The aLIGO and AdV detectors completed their first two observing runs and made the first six detections of GWs (Abbott et al. 2016c, 2016d, 2017a, 2017b, 2017c, 2017d). Two additional GW detectors are planned to join the network of aLIGO and AdV: (i) the Japanese KAGRA is under construction with baseline operations beginning in 2018 (Somiya 2012), while (ii) the proposed LIGO-India is expected to become operational in 2022 (Iyer et al. 2011; Abbott et al. 2016f). LIGO-India was approved by the government of India and a study has already suggested the site location and orientations of the arms of the detector based on scientific figures of merit (Raffai et al. 2013). These parameters, however, have not been finalized yet, and because of this we omit LIGO-India from the analysis.

Due to the expected similarities in the design sensitivities of the aLIGO, AdV, and KAGRA detectors within the frequency range of BH inspiral waveforms, for simplicity we adopt the design sensitivities of the two aLIGO (Abbott et al. 2016f) for the AdV (Abbott et al. 2016f) and KAGRA (Somiya 2012) detectors. Table 1 gives the locations and orientations of these detectors, which we used to calculate the response functions. For each detector, we define the detector's orientation angle, ψ, as the angle measured clockwise from north between the x-arm of the detector (see Appendix B for the geometric conventions of detectors) and the meridian that passes through the position of the detector.

Table 1.  Locations and Orientations of the Considered GW Detectors in the Coordinate System Defined in Appendix B

Detector	East Long.	North Lat.	Orientation ψ
LIGO H	−1194	465	−36°
LIGO L	−908	306	−108°
VIRGO	105	436	20°
KAGRA	1373	364	65°

Note. LIGO H marks the advanced LIGO detector in Hanford, WA, and LIGO L marks the advanced LIGO detector in Livingston, LA.

Download table as:  ASCIITypeset image

We assume that the noise in each detector is a stationary colored Gaussian with zero mean, and that it is uncorrelated between the different detectors. In reality, the detector noise arises from a combination of instrumental, environmental, and anthropomorphic sources that are difficult to characterize precisely (Aasi et al. 2012, 2015; Aso et al. 2013), and non-Gaussian noise transients (glitches) may arise as well (Blackburn et al. 2008). However, there are existing techniques to identify and remove glitches from GW strain channels and to reduce the level of these artifacts (Littenberg & Cornish 2010; Prestegard et al. 2012; Biswas et al. 2013; Powell et al. 2015; Bose et al. 2016; Torres-Forné et al. 2016; George et al. 2017; Mukund et al. 2017; Powell et al. 2017; Shen et al. 2017). Furthermore, correlated noise between widely separated detectors can arise from the so-called Schumann resonances (predicted in Schumann 1952a, 1952b and observed soon thereafter (Schumann & König 1954; Balser & Wagner 1960)), as well as from other EM phenomena such as solar storms, currents in the van Allen belt (Rycroft 2006), and anthropogenic emission (see Shvets et al. 2010; Thrane et al. 2013, and references therein). Note, however, that Schumann resonances mostly affect the stochastic GW background searches (Thrane et al. 2013), and a strategy against such a noise artifact already exists (Thrane et al. 2014). Our simplifying assumptions on uncorrelated Gaussian noise are therefore partly justified.

In this section, we provide a brief overview of the Fisher matrix method to estimate the measurement errors of the physical parameters characterizing a precessing eccentric BH binary source, and refer the reader to Finn (1992) and Cutler & Flanagan (1994) for further details.

The output of a GW detector, s(t), is a combination of a signal, h(t), and a noise term, $n(t),$ i.e.,

$$\begin{eqnarray}&&s(t)=h(t)+n(t).\end{eqnarray} \tag{ 38 }$$

We assume the noise of a detector to be stationary, zero mean, and Gaussian, where the different Fourier components of the noise are uncorrelated, i.e.,

$$\begin{eqnarray}&&\langle \tilde{n}(f){\tilde{n}}^{* }(f^{\prime} )\rangle =\displaystyle \frac{1}{2}\delta (f-f^{\prime} ){S}_{n}(f)\,\end{eqnarray} \tag{ 39 }$$

(Nissanke et al. 2010), where $\langle \cdot \rangle$ denotes the average, S_n(f) is the one-sided noise power spectral density of the detector, and the ${}^{* }$ superscript denotes complex conjugate. With these assumptions, the probability for the noise to have some realization ${n}_{0}(t)$ is given as

$$\begin{eqnarray}&&p(n\equiv {n}_{0})\propto {e}^{-({n}_{0}| {n}_{0})/2}\,\end{eqnarray} \tag{ 40 }$$

(Finn 1992), where p(n) is the probability distribution function of the noise to assume a value n, and $(\ldots | \ldots )$ denotes the following inner product between any two functions of frequency, e.g., x(f) and y(f):

$$\begin{eqnarray}&&(x| y)\equiv 4{\int }_{0}^{\infty }\displaystyle \frac{\tilde{x}(f)\,{\tilde{y}}^{* }(f)}{{S}_{n}(f)}\,{df}.\end{eqnarray} \tag{ 41 }$$

The optimal S/N is given by the standard expression

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}=\sqrt{(h| h)}=\sqrt{4{\int }_{0}^{\infty }\displaystyle \frac{| \tilde{h}(f){| }^{2}}{{S}_{n}(f)}{df}}.\end{eqnarray} \tag{ 42 }$$

Here, the signal waveform, h(f), depends on the parameter set $\{{\lambda }_{p}\,| \,p\in \{1\ldots P\}\}$, which characterizes the source. For a large S/N, the parameter estimation errors ${\boldsymbol{\Delta }}{\boldsymbol{\lambda }}\,=\{{\rm{\Delta }}{\lambda }_{p}\,| \,p\in \{1,\,\ldots ,\,P\}\}$ defined as the measured value minus the true value have the Gaussian probability distribution for a given signal

$$\begin{eqnarray}&&p({\boldsymbol{\Delta }}{\boldsymbol{\lambda }})={ \mathcal N }\,\exp \left(-\displaystyle \frac{1}{2}{{\rm{\Gamma }}}_{{ij}}{\rm{\Delta }}{\lambda }_{i}{\rm{\Delta }}{\lambda }_{j}\right)\end{eqnarray} \tag{ 43 }$$

(Finn 1992), where ${ \mathcal N }$ is a normalization constant. In Equation (43), we assume summation over repeated indices, and Γ_ij is the Fisher information matrix defined as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}\equiv ({\partial }_{i}h| {\partial }_{j}h)=4{\int }_{0}^{\infty }\displaystyle \frac{{\mathfrak{R}}({\partial }_{i}{\tilde{h}}^{* }(f){\partial }_{j}\tilde{h}(f))}{{S}_{n}(f)}\,{df},\end{eqnarray} \tag{ 44 }$$

where ${\partial }_{i}h=\partial h/\partial {\lambda }_{i}$ and ${\mathfrak{R}}$ labels the real part.

Following Cutler & Flanagan (1994), we define the combined S/N of a network of detectors (${\rm{S}}/{{\rm{N}}}_{\mathrm{tot}}$) as an uncorrelated superposition of individual S/Ns,

$$\begin{eqnarray}&&{({\rm{S}}/{\rm{N}})}_{\mathrm{tot}}^{2}=\displaystyle \sum _{k=1}^{{N}_{\det }}{({\rm{S}}/{\rm{N}})}_{k}^{2},\end{eqnarray} \tag{ 45 }$$

where the number of detectors in the network is denoted by N_det, and ${({\rm{S}}/{\rm{N}})}_{k}$ denotes the S/N in the $k\mathrm{th}$ detector.

Similarly, for uncorrelated Gaussian noise, the Fisher matrix of a network of detectors is the sum of the Fisher matrices of individual detectors,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij},\mathrm{tot}}=\sum _{k=1}^{{N}_{\det }}{{\rm{\Gamma }}}_{{ij},k}.\end{eqnarray} \tag{ 46 }$$

The covariance matrix is defined with the inverse of the Fisher matrix,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{ij}}={({{\rm{\Gamma }}}_{{ij},\mathrm{tot}})}^{-1}=\langle {\rm{\Delta }}{\lambda }_{i}{\rm{\Delta }}{\lambda }_{j}\rangle ,\end{eqnarray} \tag{ 47 }$$

where the angle brackets denote an average over the probability distribution function in Equation (43). The root mean square parameter measurement error σ_i in the parameters λ_i marginalized over all other parameters is

$$\begin{eqnarray}&&{\sigma }_{i}=\langle {({\rm{\Delta }}{\lambda }_{i})}^{2}{\rangle }^{1/2}=\sqrt{{{\rm{\Sigma }}}_{{ii}}}.\end{eqnarray} \tag{ 48 }$$

The off-diagonal elements of Σ_ij give the cross-correlation coefficients between the parameters λ_i and λ_j.

Parameters can be measured independently if the corresponding Fisher matrix ${{\rm{\Gamma }}}_{{ij},\mathrm{tot}}$ is nonsingular. Otherwise, if the Fisher matrix is singular, then the eigenvector(s) corresponding to the zero eigenvalue(s) of the Fisher matrix represents the linear combination(s) of the parameters, which cannot be measured by the network.

We derive efficient formulae to compute the S/N and the Fisher matrix in Appendix D.

In this section, we identify the parameters of a precessing eccentric binary that can be extracted from the detected waveform for the signal model introduced in Section 2. We set the parameters in our calculations and measure their errors as follows:

• 1.  
  ${D}_{L}:$ We set ${D}_{L}=100\,\mathrm{Mpc}$, and measure its relative error $\langle {\rm{\Delta }}{D}_{{\rm{L}}}^{2}{\rangle }^{1/2}/{D}_{{\rm{L}}}=\langle {({\rm{\Delta }}\mathrm{ln}{D}_{{\rm{L}}})}^{2}{\rangle }^{1/2}$. This choice is arbitrary, smaller than the nearest circular BH–BH merger detection to date, 340 ± 140 Mpc (Abbott et al. 2017b). The Fisher matrix method gives accurate results for the parameter measurement errors for high S/N. For moderately larger distances, the errors scale as $\propto {D}_{{\rm{L}}}$.
• 2.  
  ${\theta }_{N}$ and ϕ_N: We generate an isotropic random sample of the sky position angles θ_N and ϕ_N by drawing $\cos {\theta }_{N}$ and ϕ_N from a uniform distribution between [−1, 1] and [0, 2π], and calculate the parameter estimation covariance for each sample. The errors of the sky position is described by a localization ellipse. We characterize the sky-localization accuracy either by the corresponding proper angular length of the semimajor and semiminor axes of the sky-localization error ellipsoid given by Lang & Hughes (2006), $({a}_{N},{b}_{N})$, or its proper solid angle ${\rm{\Delta }}{{\rm{\Omega }}}_{N}=\sqrt{\pi {a}_{N}{b}_{N}}$. The calculated results are valid if ${a}_{N}\ll 1$ radian and ${b}_{N}\ll 1$ radian.
• 3.  
  ${\theta }_{L}$ and ϕ_L: We draw the angular momentum vector direction angles from an isotropic distribution and construct their error ellipsoids or solid angles similar to that given for θ_N and ϕ_N.
• 4.  
  m_A and m_B: We fix the fiducial component masses to ${m}_{A}={m}_{B}=30\,{M}_{\odot }$, consistent with the first discovered source GW150914 (Abbott et al. 2016d). Such high-mass sources are expected in GNs since mass segregation helps to increase their numbers relative to the lower-mass binaries, and the S/N is also higher for these binaries (O'Leary et al. 2009). Since we ignore additional PN corrections of the GW phase, we restrict the measurement error estimation to ${{ \mathcal M }}_{z}$ for calculations evaluated for comparison in which we ignore precession. However, generally, the assumed precessing eccentric waveform model depends on two independent combinations of component masses: ${{ \mathcal M }}_{z}$ sets the inspiral rate, and ${M}_{\mathrm{tot},z}$ sets both the apsidal precession rate and the final frequency at the LSO. We calculate the relative errors for both of these mass parameters and for the precessing eccentric waveform model $\langle {\rm{\Delta }}{{ \mathcal M }}_{z}^{2}{\rangle }^{1/2}/{{ \mathcal M }}_{z}=\langle {({\rm{\Delta }}\mathrm{ln}{{ \mathcal M }}_{z})}^{2}{\rangle }^{1/2}$, and similarly for ${M}_{\mathrm{tot},z}$.
• 5.  
  t_c, Φ_c, and γ₀: These parameters only enter in the complex phase of the waveform through Ψ_n and ${{\rm{\Psi }}}_{n}^{\pm }$ (see Equations (35) and (36)), but do not affect the S/N. Since these parameters are responsible for an overall phase shift of the waveform, we do not randomize their values but assume the fiducial value ${t}_{c}={{\rm{\Phi }}}_{c}={\gamma }_{c}=0$ for each binary in the Monte Carlo sample.
• 6.  
  e_LSO: The adopted eccentric inspiral waveform model depends explicitly on the final eccentricity at the LSO; see Equation (21). This quantity parameterizes the evolutionary path of the binary during its eccentric inspiral in the $({\rho }_{{\rm{p}}},e)$ plane as shown in Equations (20) and (21); see also Figure 3 in Kocsis & Levin (2012) for illustration. In fact, any segment of the evolutionary path ${\rho }_{{\rm{p}}}(e)$ specifies the value of e_LSO uniquely. Conversely, e_LSO specifies ${\rho }_{{\rm{p}}}(e)$, which sets a constraint on the possible values of $({\rho }_{{\rm{p}}0},{e}_{0})$, if the PN binary inspiral model is extrapolated backwards in time. Indeed, in some cases, this is the only indirect information we may have on the formation parameters $({\rho }_{{\rm{p}}0},{e}_{0})$. In particular, ${e}_{0}\leqslant 1$ puts an upper bound on ${\rho }_{{\rm{p}}0}$ for a given e_LSO.¹⁰
• 7.  
  e₀: We choose several e₀ values from the highly eccentric (e₀ ≥ 0.9) limit when discussing the e₀ dependence of the measurement errors (see Section 6.2). However, we restrict to e₀ = 0.9 for calculations of a large survey of binaries.¹¹
• 8.  
  ${\rho }_{{\rm{p}}0}$: We examine two values for the dimensionless initial pericenter distance ${\rho }_{{\rm{p}}0}=\{10,20\}$, and the circular limit corresponds to ${\rho }_{{\rm{p}}0}\to \infty$ (O'Leary et al. 2009). These values are likely for sources that form through the GW capture mechanism in high-velocity dispersion environments such as in GN, as shown in O'Leary et al. (2009) and Gondán et al. (2017), or the core-collapsed regions of star clusters without a central massive BH (Kocsis et al. 2006b; Antonini & Rasio 2016).

If the peak GW frequency of the initial orbit is large enough to be in the detectors' sensitive frequency band, then ${\rho }_{{\rm{p}}0}$ and e₀ are directly measurable due to the truncation of the time-domain waveform for times when e < e₀ and ${\rho }_{{\rm{p}}}\gt {\rho }_{{\rm{p}}0}$. In the opposite case, only a lower limit may be given for ${\rho }_{{\rm{p}}0}$, which corresponds to ${e}_{0}\to 1$ (Kocsis & Levin 2012).

In summary, we use the following free parameters in the Fisher matrix analysis,

$$\begin{eqnarray}&&{{\boldsymbol{\lambda }}}_{\mathrm{Prec}}=\{\mathrm{ln}({D}_{L}),\mathrm{ln}({{ \mathcal M }}_{z}),\mathrm{ln}({M}_{\mathrm{tot},z}),{\theta }_{N},{\phi }_{N},{\theta }_{L},{\phi }_{L},{e}_{0},\\ &&\qquad {e}_{\mathrm{LSO}},{t}_{c},{{\rm{\Phi }}}_{c},{\gamma }_{c}\}.\end{eqnarray} \tag{ 49 }$$

Given these parameters, the other parameters' marginalized measurement errors may be determined by linear combinations of the covariance matrix based on Equation (48). For example, ${\rho }_{{\rm{p}}0}$ is given by e₀ and e_LSO using Equations (20) and (21). Its measurement error is

$$\begin{eqnarray}\langle {({\rm{\Delta }}{\rho }_{{\rm{p}}0})}^{2}\rangle & = & {\left(\displaystyle \frac{\partial {\rho }_{{\rm{p}}}({e}_{0},{e}_{\mathrm{LSO}})}{\partial {e}_{0}}\right)}^{2}\langle {({\rm{\Delta }}{e}_{0})}^{2}\rangle \\ & & +{\left(\displaystyle \frac{\partial {\rho }_{{\rm{p}}}({e}_{0},{e}_{\mathrm{LSO}})}{\partial {e}_{\mathrm{LSO}}}\right)}^{2}\langle {({\rm{\Delta }}{e}_{\mathrm{LSO}})}^{2}\rangle \\ & & +2\displaystyle \frac{\partial {\rho }_{{\rm{p}}}({e}_{0},{e}_{\mathrm{LSO}})}{\partial {e}_{0}}\displaystyle \frac{\partial {\rho }_{{\rm{p}}}({e}_{0},{e}_{\mathrm{LSO}})}{\partial {e}_{\mathrm{LSO}}}\langle {\rm{\Delta }}{e}_{0}{\rm{\Delta }}{e}_{\mathrm{LSO}}\rangle .\end{eqnarray} \tag{ 50 }$$

The parameter estimation errors of individual component masses or the mass ratio can be estimated similarly using ${\rm{\Delta }}{M}_{\mathrm{tot}}$ and ${\rm{\Delta }}{ \mathcal M }$ after inverting ${M}_{\mathrm{tot}}({m}_{a},{m}_{b})$ and ${ \mathcal M }({m}_{a},{m}_{b})$.

The measurement errors depend on the sky position of the source with respect to the detectors and on the relative orientation of the binary. We generate random Monte Carlo samples of ∼4500 binaries by drawing from isotropic distributions of the sky position and of the binary-orientation normal vector. We present the results for the ${\rm{S}}/{{\rm{N}}}_{\mathrm{tot}}$ for detecting precessing highly eccentric BH binaries with the GW detector network described in Table 1 and for the expected parameter measurement errors.

6.1. S/N Distributions

Figure 1 displays the distribution of the S/N_tot for precessing highly eccentric BH binaries detected with the detector network described in Table 1, assuming binary parameters ${m}_{A}={m}_{B}\,=30\,{M}_{\odot }$, e₀ = 0.9, ${\rho }_{{\rm{p}}0}=\{10,20\}$, and similar binaries in the circular limit (see Appendix A for details). Generally, similar to the results of O'Leary et al. (2009; see Figure 11 therein, which corresponds to a single aLIGO detector), the S/N_tot is systematically higher for binaries with ${\rho }_{{\rm{p}}0}=10$ than for binaries with ${\rho }_{{\rm{p}}0}=20$. We find that increasing the initial eccentricity from e₀ = 0.9 to 0.97 for a fixed ${\rho }_{{\rm{p}}0}$ does not change the S/N_tot significantly (see Table 3 below), hence we expect that the distribution of S/N_tot for fixed ρ_p0 converges in the ${e}_{0}\to 1$ limit. This is expected since Figure 10 in O'Leary et al. (2009) shows that a low value of the S/N_tot accumulates near e₀ ≈ 1 for low to moderately high BH masses.

Zoom In Zoom Out Reset image size

Figure 2 shows that the S/N_tot increases rapidly with ρ_p0 for low ρ_p0, has a maximum between ρ_p0 ∼ 9 and 20, and decreases for higher ρ_p0 approaching the circular binary limit for ${\rho }_{{\rm{p}}0}\to \infty$. These findings may be understood qualitatively as follows. Within ${\rho }_{{\rm{p}}0}\leqslant 40.9{({M}_{\mathrm{tot},z}/20{M}_{\odot })}^{-2/3}$, the binary forms with a characteristic frequency above 10 Hz in the detector band (see Equation (59) in Gondán et al. 2017). The rapid decrease of the S/N_tot for decreasing ${\rho }_{{\rm{p}}0}\lt 9$ is due to the fact that we ignore the GWs of the first hyperbolic encounter (Kocsis et al. 2006b). The decrease of the S/N_tot at high ρ_p0 is due to the fact that part of the GW spectrum falls outside of the detectors' sensitive frequency band. For very large ρ_p,0, the binary becomes circular by the time it enters the detectors' sensitive frequency band, and a significant fraction of the S/N_tot accumulates only in the n = 2 harmonic (Figure 9), which explains the flat asymptotics for high ρ_p0. As ρ_p0 decreases from high to moderate values, higher harmonics start to contribute to the S/N_tot (Figure 9), which explains the increase of the S/N_tot. A combination of these arguments leads to the peak of the S/N_tot at an intermediate ρ_p0 value seen in Figure 2. However, note that in addition to ignoring the initial hyperbolic encounter and the final coalescence/ringdown segments of the signal, our waveform model also ignores contributions of spherical moments beyond the quadrupole order and deviations from a precessing Keplerian orbit (Davis et al. 1972; Berti et al. 2010; Healy et al. 2016). This approximation may not be valid for low ρ_p0 (particularly for ρ_p0 ≤ 10). The S/N _tot is expected to be underestimated in this region in Figure 5.

Zoom In Zoom Out Reset image size

6.2. Parameter Measurement Errors

We present the measurement accuracy for the final eccentricity at the LSO for parameters grouped as

$$\begin{eqnarray}&&{{\boldsymbol{\lambda }}}_{\mathrm{slow}}=\{\mathrm{ln}({D}_{L}),{\theta }_{N},{\phi }_{N},{\theta }_{L},{\phi }_{L}\},\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\lambda }}}_{\mathrm{fast}}=\{{{\rm{\Phi }}}_{c},{t}_{c},\mathrm{ln}({{ \mathcal M }}_{z}),\mathrm{ln}({M}_{\mathrm{tot},z}),{e}_{\mathrm{LSO}},{\gamma }_{c}\}\,\ \end{eqnarray} \tag{ 52 }$$

(Kocsis et al. 2007). The ${{\boldsymbol{\lambda }}}_{\mathrm{fast}}$ fast parameters are related to the high-frequency GW phase, while the ${{\boldsymbol{\lambda }}}_{\mathrm{slow}}$ slow parameters appear only in the slowly varying amplitude of the GW signal. Slow parameters are mostly determined from a comparison of the GW signals measured by the different detectors in the network. For the polar angles $({\theta }_{N},{\phi }_{N})$ describing the source direction, we calculate the minor and major axes (${a}_{N},{b}_{N}$) of the corresponding 2D sky location error ellipse and its area (${{\rm{\Omega }}}_{N}=\pi {a}_{N}{b}_{N}$), and we do the same for the binary-orientation error ellipse $({a}_{L},{b}_{L})$ and its area (${{\rm{\Omega }}}_{N}=\pi {a}_{L}{b}_{L}$).

Figures 3 and 4 show the distribution of the measurement errors for a randomly chosen source sky position and binary orientation for ${\rho }_{{\rm{p}}0}=10$ and 20, and for the circular limit ${\rho }_{{\rm{p}}0}\to \infty$ (see Appendix A), while Table 2 shows the 10%, 50%, and 90% quantiles of the error distributions. Compared to a ρ_p0 = 20 binary, a ρ_p0 = 10 binary is more eccentric throughout its evolution, which leads to a higher S/N_tot, and most of its measurement errors are smaller. There are, however, exceptions to this finding: the fast parameters such as the mass parameters and the eccentricity have higher errors for ρ_p0 = 10 than for ρ_p0 = 20 (see the discussion below).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  The 10%, 50%, and 90% Quantile of Measurement Errors for Parameters of 30 ${M}_{\odot }$–30 ${M}_{\odot }$ Precessing Eccentric BH Binaries with Initial Eccentricity e₀ = 0.9 and Dimensionless Pericenter Distance ${\rho }_{{\rm{p}}0}=10$ and 20, and Circular Binaries at Distance ${D}_{L}=100\,\mathrm{Mpc}$, Random Source Sky Location, and Orientation Using the Detector Network in Table 1

ρp0	10	10	10	20	20	20	Circular	Circular	Circular
Quantile	10%	50%	90%	10%	50%	90%	10%	50%	90%
${\rm{\Delta }}(\mathrm{ln}{D}_{L})$	1.13(–2)	2.49(−2)	7.21(−2)	1.89(−2)	4.36(−2)	0.154	8.81(−3)	5.86(−2)	0.311
${\rm{\Delta }}{{\rm{\Omega }}}_{N}\,[\mathrm{sr}]$	1.13(−5)	1.46(−4)	7.23(−4)	2.21(−4)	7.28(−4)	2.89(−3)	5.65(−5)	1.58(−3)	6.67(−3)
${\rm{\Delta }}{{\rm{\Omega }}}_{L}\,[\mathrm{sr}]$	7.89(−5)	3.85(−3)	0.13	2.21(−5)	1.13(−2)	0.85	4.36(−4)	2.21(−2)	9.6
${a}_{N}\,[\deg ]$	0.14	0.62	2.42	0.57	1.36	4.07	0.28	2.01	5.16
${b}_{N}\,[\deg ]$	0.11	0.22	0.48	0.26	0.52	0.98	0.41	0.76	1.44
${a}_{L}\,[\deg ]$	0.96	2.78	20.04	1.44	4.65	45.48	0.78	6.98	2.02( + 2)
${b}_{L}\,[\deg ]$	0.55	1.33	7.51	0.84	2.46	20.58	0.32	3.41	50.63
${\rm{\Delta }}{{\rm{\Phi }}}_{c}\,[\mathrm{rad}]$	9.27(−2)	0.23	0.71	0.29	0.66	2.44	0.36	0.72	57.91
${\rm{\Delta }}{t}_{c}\,[\mathrm{ms}]$	4.32(−2)	8.41(−2)	0.181	9.28(−2)	0.167	0.311	9.41(−2)	1.582	3.011
${\rm{\Delta }}(\mathrm{ln}{{ \mathcal M }}_{z})$	3.53(−5)	6.17(−5)	1.71(−4)	1.04(−5)	1.81(−5)	3.42(−5)	1.36(−3)	2.34(−3)	4.27(−3)
${\rm{\Delta }}(\mathrm{ln}{M}_{\mathrm{tot},z})$	5.42(−4)	9.51(−4)	2.43(−3)	2.82(−4)	4.81(−4)	9.18(−4)	5.88(−3)	1.13(−2)	1.81(−2)
${\rm{\Delta }}{e}_{\mathrm{LSO}}$	1.18(−4)	2.16(−4)	5.83(−4)	2.39(−5)	3.19(−5)	5.88(−5)	⋯	⋯	⋯
Δ e0	1.44(−3)	2.16(−3)	3.95(−3)	1.72(−3)	2.91(−3)	5.79(−3)	⋯	⋯	⋯
${\rm{\Delta }}{\rho }_{{\rm{p}}0}$	6.64(−3)	1.08(−2)	2.28(−2)	1.33(−2)	2.29(−2)	4.58(−2)	⋯	⋯	⋯
${\rm{\Delta }}{\gamma }_{c}\,[\mathrm{rad}]$	0.04	0.11	0.47	0.14	0.34	1.51	⋯	⋯	⋯

Note. Here, $({{\rm{\Omega }}}_{N},{a}_{N},{b}_{N})$ are, respectively, the area, and the semimajor and semiminor axes of the 2D error ellipse corresponding to the source's sky direction, and similarly for $({{\rm{\Omega }}}_{L},{a}_{L},{b}_{L})$ describing the source's orbital plane orientation (i.e., angular momentum vector direction). Note that ${e}_{\mathrm{LSO}}=0.187$ and 0.059 for ${\rho }_{{\rm{p}}0}=10$ and 20, respectively, if e₀ = 0.9 is assumed. In the circular limit, the binary forms outside of the sensitive frequency band of the detector network, and ${\rm{\Delta }}{e}_{0}\to \infty$ and ${\rm{\Delta }}{\rho }_{{\rm{p}}0}\to \infty$. We adopt the following notation in the table: $1.13(-2)=1.13\times {10}^{-2}$.

Download table as:  ASCIITypeset image

Many of the binaries in GNs¹² form with very high e₀, close to unity, in single–single encounters due to GW emissions (O'Leary et al. 2009; Gondán et al. 2017). However, similar to the finding that the S/N_tot does not increase significantly for 1 > e₀ ≥ 0.9, we find that ${{\boldsymbol{\lambda }}}_{\mathrm{slow}}$ parameter errors do not improve due to the early very eccentric evolutionary period beyond e > 0.9 (repeated burst phase) compared to waveforms with e₀ = 0.9 as shown in Table 3. However, some ${{\boldsymbol{\lambda }}}_{\mathrm{fast}}$ parameters' measurement errors improve more significantly with e₀ > 0.9. In particular, the measurement errors of the mass parameters (${{ \mathcal M }}_{z},{M}_{\mathrm{tot},z}$) improve by a factor of ∼2, and the measurement error of e_LSO improves by ∼50% if increasing e₀ from 0.9 for 0.97. This difference is due to the fact that eccentricity modifies the GW phase significantly, which affects the determination of the ${{\boldsymbol{\lambda }}}_{\mathrm{fast}}$ parameters only.

Table 3.  Measurement Errors for Parameters of 30 ${M}_{\odot }$–30 ${M}_{\odot }$ Precessing Eccentric BH Binaries with Initial Eccentricities e₀ = 0.9, 0.95, and 0.97 for Initial Dimensionless Pericenter ${\rho }_{{\rm{p}}0}=10$ and 20, Luminosity Distance D_L = 100 Mpc, and Arbitrarily Fixed Source Direction $({\theta }_{N},{\phi }_{N})=(\pi /2,\pi /3)$ and Orientation $({\theta }_{L},{\phi }_{L})=(\pi /4,\pi /5)$ for the Detector Network Specified in Table 1

${\rho }_{{\rm{p}}0}$	10	10	10	20	20	20
e0	0.9	0.95	0.97	0.9	0.95	0.97
eLSO	0.1872	0.1932	0.1956	5.89(−2)	6.07(−2)	6.15(−2)
${\rm{S}}/{{\rm{N}}}_{\mathrm{tot}}$	251.1	257.7	260.4	179.4	182.1	182.9
${\rm{\Delta }}(\mathrm{ln}{D}_{L})$	3.92(−2)	3.79(−2)	3.76(−2)	6.41(−2)	6.31(−2)	6.27(−2)
${\rm{\Delta }}{{\rm{\Omega }}}_{N}\,[\mathrm{sr}]$	8.37(−5)	7.78(−5)	7.53(−5)	4.06(−4)	3.91(−4)	3.85(−4)
${\rm{\Delta }}{{\rm{\Omega }}}_{L}\,[\mathrm{sr}]$	9.68(−3)	9.07(−3)	8.84(−3)	2.74(−2)	2.65(−2)	2.61(−2)
${a}_{N}\,[\deg ]$	0.501	0.484	0.477	0.994	0.977	0.970
${b}_{N}\,[\deg ]$	0.174	0.168	0.165	0.427	0.418	0.415
${a}_{L}\,[\deg ]$	4.31	4.18	4.13	6.96	6.85	6.81
${b}_{L}\,[\deg ]$	2.35	2.26	2.23	4.12	4.04	4.01
${\rm{\Delta }}{{\rm{\Phi }}}_{c}\,[\mathrm{rad}]$	0.331	0.313	0.308	1.102	1.049	1.031
${\rm{\Delta }}{t}_{c}\,[\mathrm{ms}]$	9.68(−2)	9.08(−2)	8.88(−2)	0.198	0.192	0.189
${\rm{\Delta }}(\mathrm{ln}{{ \mathcal M }}_{z})$	4.61(−5)	1.89(−5)	1.39(−5)	1.38(−5)	6.03(−6)	4.33(−6)
${\rm{\Delta }}(\mathrm{ln}{M}_{\mathrm{tot},z})$	7.09(−4)	5.62(−4)	4.85(−4)	3.67(−4)	2.87(−4)	2.38(−4)
${\rm{\Delta }}{e}_{\mathrm{LSO}}$	1.61(−4)	1.28(−4)	1.11(−4)	2.43(−5)	1.93(−5)	1.62(−5)
Δ e0	1.95(−3)	1.74(−3)	1.69(−3)	2.25(−3)	2.17(−3)	2.14(−3)
${\rm{\Delta }}{\rho }_{{\rm{p}}0}$	8.33(−3)	7.45(−3)	7.11(−3)	1.71(−2)	1.61(−2)	1.56(−2)
${\rm{\Delta }}{\gamma }_{c}\,[\mathrm{rad}]$	0.169	0.162	0.159	0.556	0.529	0.520

Download table as:  ASCIITypeset image

We calculate the parameter measurement errors for precessing highly eccentric BH binaries as a function of ρ_p0 for some arbitrarily fixed binary direction and orientation. For one such binary direction and orientation, Figure 5 shows the ρ_p0 dependence of ${\rm{\Delta }}{D}_{{\rm{L}}}/{D}_{{\rm{L}}}$, ${\rm{\Delta }}{{ \mathcal M }}_{z}/{{ \mathcal M }}_{z}$, ${\rm{\Delta }}{e}_{0}$, ${\rm{\Delta }}{\rho }_{{\rm{p}}0}$, ${\rm{\Delta }}{e}_{\mathrm{LSO}}$, semimajor and semiminor axes of the sky position error ellipse $({a}_{N},{b}_{N})$, and semimajor and semiminor axes of the error ellipse for the binary orbital plane normal vector direction $({a}_{L},{b}_{L})$. Note that we find similar trends with ${\rho }_{{\rm{p}}0}$ for other random choices of binary direction and orientation. We find that measurement errors systematically decrease with decreasing ρ_p0 for precessing highly eccentric binaries relative to similar binaries in the circular limit (Appendix A), and the errors have a minimum in the range 8 < ρ_p0 < 80 and deteriorate rapidly for ρ_p0 < 8. The latter is due to the rapid decrease of the S/N_tot in that range. The ρ_p0 dependence of ${\rm{\Delta }}{D}_{L}/{D}_{L}$ and the principal axes of the sky-position and binary-orientation error ellipses $({a}_{N},{b}_{N})$ and $({a}_{L},{b}_{L})$ (i.e., the quantities derived from the slow parameters; see Section 5) are qualitatively similar to that of 1/(S/N_tot) in the complete range of ρ_p0 (i.e., they decrease rapidly with ρ_p0 for low ρ_p0, have a minimum at moderate ρ_p0, and converge asymptotically to the value of precessing highly eccentric binaries in the circular limit); see Figure 6 for details. However, Figure 5 shows that the chirp mass errors have a minimum at a much higher ρ_p0, i.e., between 50–60 and 20–40 for 10 ${M}_{\odot }$–10 ${M}_{\odot }$ and 30 ${M}_{\odot }$–30 ${M}_{\odot }$ precessing highly eccentric BH binaries, respectively. The main reason for the different behavior of the chirp mass from the distance and angular errors is the chirp mass is a fast parameter, while the distance and angular parameters are slow parameters. Slow parameters are insensitive to the GW phase perturbations and depend on the GW amplitude, which is set by the S/N_tot. The S/N of the early part of the waveform near the low-frequency noise wall of the detector is small. However, fast parameters depend sensitively on the GW phase, and the GW phase accumulates mostly at low frequencies, since the residence time (i.e., $\nu /\dot{\nu }$) is largest at low orbital frequencies. Thus, the fast parameters' errors are minimized for binaries that form with e₀ ∼ 1 with a ρ_p0 value for which the GW characteristic frequency is near the detectors' minimum frequency. The peak of the spectrum is initially at f_min if ${\rho }_{{\rm{p}}0}=40.9{[{M}_{\mathrm{tot}}/{(20{M}_{\odot })]}^{-2/3}({f}_{\min }/10\mathrm{Hz})}^{-2/3}$ (Gondán et al. 2017). A slightly lower value of ${f}_{\min }\sim 7\,\mathrm{Hz}$ leads to values that represent the minimum of the fast parameters. This also leads to the result observed in Figure 3 that these parameters have higher errors for ρ_p0 = 10 than for ρ_p0 = 20.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 5, note that the Δe_LSO errors are relatively small for relatively high ρ_p0 up to ρ_p0 ∼ 200. At high ρ_p0, the orbital eccentricity approaches zero when it enters the aLIGO band, and Δe_LSO increases. We note that the posterior probability distribution function of e_LSO is well-defined even in the circular limit ${\rho }_{{\rm{p}}0}\to \infty$, and Δe_LSO is finite for a given confidence region. However, the Fisher matrix algorithm becomes invalid in this regime as the signal is not approximated well by its linear Taylor expansion with respect to the Δe_LSO parameter, since its first e_LSO derivative vanishes in the circular limit. Therefore, the true asymptotic value of Δe_LSO for high ρ_p0 cannot be recovered with the Fisher matrix technique used in this paper. Further, note that Δe₀ and Δρ_p0 also increase rapidly with ρ_p0 for high ρ_p0. This is due to the fact that for these parameters, the binary forms with a pericenter frequency smaller than the minimum frequency of the detector network, and the information on e₀ and ρ_p0 is limited to higher harmonics with small power. Thus, these parameters indeed have a very high error and become indeterminate in the circular limit. The fact that the relative error of e₀ and ρ_p0 can be less than ∼5% percent in the range 5 < ρ_p0 < 50 (6 < ρ_p0 < 100) for 30 ${M}_{\odot }$–30 ${M}_{\odot }$ (10 ${M}_{\odot }$–10 ${M}_{\odot }$) precessing highly eccentric BH binaries implies that the GW detections might have the potential to constrain the formation environment of these systems (O'Leary et al. 2009; Cholis et al. 2016; Rodriguez et al. 2016a, 2016b; Chatterjee et al. 2017; Gondán et al. 2017; Kocsis et al. 2017; Samsing & Ramirez-Ruiz 2017; Silsbee & Tremaine 2017).

Furthermore, we found from numerical investigations that Δe₀ does not correlate significantly with other parameters' errors, which is due to the fact that e₀ is measured from the truncation of the signal for e > e₀ at the start of the waveform, while other parameters of a precessing eccentric binary are measured from the inspiral rate (Section 5). However, Δρ_p0 behaves differently from Δe₀ in this regard, which is due to the fact that ρ_p0 is determined by e_LSO in Equation (50), and e_LSO depends on the mass parameters.

6.3. Comparison with Previous Results

In this paper, we have determined the S/N and the expected accuracy with which the aLIGO–AdV–KAGRA detector network may determine the parameters that describe highly eccentric BH binaries, and investigated how these quantities depend on the initial pericenter distance ρ_p0 and initial eccentricity e₀. There are some previous studies that also made similar investigations for eccentric compact binaries with significant differences (Yunes et al. 2009; Kyutoku & Seto 2014; Sun et al. 2015; Ma et al. 2017). They considered different detector networks, applied different waveform models, and used different definitions for e₀ and ρ_p0. As a consequence, only a qualitative comparison is possible with those results, which we discuss in this section. At the end of this section, we compare our results for the measurement errors in the circular limit with those presented in previous studies.

We first compare our results with a previous study for the ρ_p0 dependence of the S/N_tot. Our result for the ρ_p0 dependence of the S/N_tot (Figure 2) is qualitatively in agreement with the result of Figure 2 in Kyutoku & Seto (2014), i.e., the S/N_tot increases rapidly with ρ_p0 for low ρ_p0, peaks at a moderate ρ_p0, and converges asymptotically to the value of highly eccentric binaries in the circular limit for high ρ_p0.

In order to compare our results for the e₀ dependence of the S/N_tot with Yunes et al. (2009) and Sun et al. (2015), we also set the lower bound of the advanced GW detectors' sensitive frequency band to ${f}_{\min }\,=\,20\,\mathrm{Hz}$. We define ${e}_{20\mathrm{Hz}}$ to be the eccentricity at which the peak GW frequency of the binary defined in Wen (2003) is ${f}_{\mathrm{GW}}=20\,\mathrm{Hz}$, and evaluate ρ_p0 corresponding to ${e}_{0}={e}_{20\mathrm{Hz}}$ from Equation (37) in Wen (2003) as

$$\begin{eqnarray}&&{\rho }_{20\mathrm{Hz}}={[{(1+{e}_{20\mathrm{Hz}})}^{0.3046}{f}_{\mathrm{GW}}\pi {M}_{\mathrm{tot},z}]}^{-2/3}.\end{eqnarray} \tag{ 53 }$$

We recalculate the distribution of the S/N_tot for 10 ${M}_{\odot }$–10 ${M}_{\odot }$ precessing eccentric compact binaries with ${e}_{20\mathrm{Hz}}\,=(0.1,0.2,0.3,0.4)$. The top panel of Figure 7 shows that the S/N_tot is roughly the same for different ${e}_{20\mathrm{Hz}}$, which is consistent with the results presented in Figure 2 in Sun et al. (2015) and in the left panel of Figure 8 in Yunes et al. (2009). Moreover, we find that the S/N_tot increases weakly with ${e}_{20\mathrm{Hz}}$, which is in agreement with the results in the left panel of Figure 8 in Yunes et al. (2009). Note that this result disagrees with Table 5 in Sun et al. (2015). The S/N_tot does not depend significantly on ${e}_{20\mathrm{Hz}}$ in the range of [0.1,0.4] for 10 ${M}_{\odot }$–10 ${M}_{\odot }$ for ${\rho }_{{\rm{p}}0}={\rho }_{20\mathrm{Hz}}\sim 28$ as seen in the bottom panel of Figure 7.¹³ For 10 ${M}_{\odot }$–10 ${M}_{\odot }$ binaries with e₀ < 0.4 and ρ_p0 ∼ 28, we find that ${\rho }_{20\mathrm{Hz}}$ is high enough to fall into the range of ρ_p0, where the S/N_tot depends on e₀ at the ∼10% level for ${e}_{20\mathrm{Hz}}\lt 0.4$. The influence of ${e}_{20\mathrm{Hz}}$ on the S/N_tot increases with M_tot since in this case ${\rho }_{20\mathrm{Hz}}$ is lower as shown by Equation (53). Thus, the influence of ${e}_{20\mathrm{Hz}}$ on the distribution of S/N_tot is more significant for higher-mass low-eccentricity binaries. Examples for this characteristic of the S/N_tot are seen in Figure 8 in Yunes et al. (2009).

Zoom In Zoom Out Reset image size

Finally, we compare our results with previous studies for the e₀ dependence of the measurement errors of the parameters describing eccentric binaries. Sun et al. (2015) and Ma et al. (2017) set the initial orbital parameters to be ${e}_{20\mathrm{Hz}}$ (${e}_{10\mathrm{Hz}}$) and ${\rho }_{20\mathrm{Hz}}$ (${\rho }_{10\mathrm{Hz}}$). For various values of ${e}_{20\mathrm{Hz}}$ (${e}_{10\mathrm{Hz}}$) in the range [0.1, 0.2, 0.3, 0.4] and the corresponding values of ${\rho }_{20\mathrm{Hz}}$ (${\rho }_{10\mathrm{Hz}}$), they determined the measurement accuracies for various parameters of eccentric binaries. Since they applied different waveform models and different parameters describing the eccentric binaries, we resort to a qualitative comparison. We repeated the analysis of Figure 5 for e₀ = (0.1, 0.2, 0.3, 0.4) and determined the measurement error of parameters as a function of ρ_p0 as shown in Figure 8. We find qualitative agreement with Sun et al. (2015) and Ma et al. (2017). The measurement accuracies of the parameters increase strictly monotonically with e₀.

Zoom In Zoom Out Reset image size

Previous papers have investigated the M_tot dependency of the measurement errors for $\{{t}_{c},{{\rm{\Phi }}}_{c},\mathrm{ln}({{ \mathcal M }}_{z}),\mathrm{ln}(\eta )\}$ by using different PN order waveform models for non-spinning inspiraling binaries for a fixed S/N in a single aLIGO-type detector. Previous results showed that the measurement accuracy of these parameters decreases with increasing M_tot for 2.8 ${M}_{\odot }$ ≤ M_tot ≤ 20 ${M}_{\odot }$, provided that the S/N accumulated in one GW detector is fixed; see Table 1 in Arun et al. (2005) and references therein. Therefore, we determined the measurement errors in the circular limit for $\{{t}_{c},{{\rm{\Phi }}}_{c},\mathrm{ln}({{ \mathcal M }}_{z}),\mathrm{ln}(\eta )\}$ for a qualitative comparison.¹⁴ To calculate the measurement error of the ln (η) parameter, we use the fact that $\eta ={({{ \mathcal M }}_{z}{M}_{\mathrm{tot},z}^{-1})}^{5/3}$ and so

$$\begin{eqnarray}\displaystyle \frac{\langle {\rm{\Delta }}{\eta }^{2}\rangle }{{\eta }^{2}} & = & \displaystyle \frac{25}{9}\displaystyle \frac{\langle {\rm{\Delta }}{{ \mathcal M }}_{z}^{2}\rangle }{{{ \mathcal M }}_{z}^{2}}+\displaystyle \frac{25}{9}\displaystyle \frac{\langle {\rm{\Delta }}{M}_{\mathrm{tot},z}^{2}\rangle }{{M}_{\mathrm{tot},{\rm{z}}}^{2}}\\ & & +\displaystyle \frac{50}{9}\displaystyle \frac{\langle {\rm{\Delta }}{{ \mathcal M }}_{z}{\rm{\Delta }}{M}_{\mathrm{tot},z}\rangle }{{{ \mathcal M }}_{z}{M}_{\mathrm{tot},z}}.\end{eqnarray} \tag{ 54 }$$

In agreement with the 1PN-order case in Arun et al. (2005), we find that Δt_c, ΔΦ_c, ${\rm{\Delta }}{{ \mathcal M }}_{z}/{{ \mathcal M }}_{z}$, and Δη/η increase with M_tot for fixed S/N_tot (Table 4). Such a qualitative agreement is expected since the adopted precessing eccentric waveform approximates the full 1PN waveform in its most important features, and the M_tot-dependent trends of error distributions do not depend on the number of detectors or on the sky position or angular momentum unit vectors of the source.

Table 4.  Errors in t_c (ms) and Φ_c (rad) and the Relative Errors in ${{ \mathcal M }}_{z}$ and η in the Circular Limit for Equal-mass Binaries for a Specific Sky Position and Inclination ${\theta }_{N}=\pi /2$, ${\phi }_{N}=\pi /3$, ${\theta }_{L}=\pi /4$, and ${\phi }_{L}=\pi /5$

${m}_{A}-{m}_{B}$	${\rm{\Delta }}{t}_{c}$	${\rm{\Delta }}{{\rm{\Phi }}}_{c}$	${\rm{\Delta }}{{ \mathcal M }}_{z}/{{ \mathcal M }}_{z}$	${\rm{\Delta }}\eta /\eta$
10 ${M}_{\odot }$–10 ${M}_{\odot }$	0.30	0.32	3.6 × 10−4	3.4 × 10−3
15 ${M}_{\odot }$–15 ${M}_{\odot }$	0.55	0.39	8.7 × 10−4	6.5 × 10−3
20 ${M}_{\odot }$–20 ${M}_{\odot }$	0.94	0.49	1.5 × 10−3	10−2
25 ${M}_{\odot }$–25 ${M}_{\odot }$	1.38	0.57	2.2 × 10−3	1.3 × 10−2
30 ${M}_{\odot }$–30 ${M}_{\odot }$	1.84	0.62	2.8 × 10−3	1.5 × 10−2

Note. In each case, we have assumed detection with the detector network introduced in Table 1, and the errors correspond to a fixed ${\rm{S}}/{{\rm{N}}}_{{\rm{t}}{\rm{o}}{\rm{t}}}=100$. We find similar trends for other random choices of binary directions and orientations (not shown).

Download table as:  ASCIITypeset image

We carried out a Fisher-matrix-type study to determine the accuracy with which the parameters of highly eccentric BH binaries may be measured using the aLIGO–AdV–KAGRA GW detector network. Eccentricity changes the GWs of binaries compared to circular binaries in several ways. In the time domain, the gravitational waveform of eccentric binaries is quasiperiodic but not sinusoidal. Relativistic precession adds a slow amplitude modulation to the waveform for each polarization. Eccentricity also changes the inspiral rate at which the binary separation and period shrink. We take all of these effects into account using the stationary phase approximation (Moreno-Garrido et al. 1994; Mikóczi et al. 2012). In contrast to circular binaries, the waveform of eccentric binaries includes several prominent orbital frequency harmonics, general relativistic precession causes each harmonic to split into three frequencies for both GW polarizations respectively, and the eccentric inspiral creates a spectrum, which is different from the $\tilde{h}\propto {f}^{-7/6}$ waveform of circular inspiral sources for each harmonic. These features in the waveform make it possible to accurately determine the eccentricity and angle of periapsis, and the modulated inspiral rate improves the measurement accuracy of mass parameters for eccentric inspirals.

The main parameters that describe eccentric inspiraling binaries are the initial pericenter distance ρ_p0 when the eccentricity is close to unity and the final eccentricity at the LSO e_LSO. These parameters are systematically different for different formation channels. Thus, their measurements may have important implications on the astrophysical origin of the sources (O'Leary et al. 2009; Cholis et al. 2016; Rodriguez et al. 2016a, 2016b; Chatterjee et al. 2017; Gondán et al. 2017; Samsing & Ramirez-Ruiz 2017; Silsbee & Tremaine 2017). Based on a survey with 10 ${M}_{\odot }$–10 ${M}_{\odot }$ and 30 ${M}_{\odot }$–30 ${M}_{\odot }$ precessing highly eccentric BH binaries at 100 Mpc using the planned aLIGO–AdV–KAGRA detector network, our results can be summarized as follows:

• 1.  
  The S/N_tot improves by a factor of ∼1–1.7 (depending on ρ_p0, the component masses, and the sky-position and binary-orientation angles; see Figure 1) for 30 ${M}_{\odot }$–30 ${M}_{\odot }$ precessing highly eccentric BH binaries compared to similar binaries in the circular limit¹⁵ with the same masses and distance. The volume in the universe for a fixed maximum S/N is $\langle {({\rm{S}}/{\rm{N}})}^{3}\rangle \sim 5\times (2\times )$ larger for eccentric inspiraling binaries with ρ_p0 = 10 (ρ_p0 = 20) than for similar binaries in the circular limit.
• 2.  
  We determined how the parameters' measurement accuracies depend on the initial dimensionless pericenter distance (ρ_p0) for precessing highly eccentric BH binaries. The smallest errors are obtained for small ρ_p0 < 10 values for the sky position and angular momentum and ρ_p0 < 20 for the luminosity distance and ρ_p0. However, the errors for fast parameters, which are sensitive to the GW phase, like the chirp mass, the initial eccentricity, and the eccentricity at the LSO, improve most significantly for a higher ρ_p0 between 10 and 80 (Figure 5).
• 3.  
  The parameter estimation errors can improve significantly for highly eccentric precessing BH binaries compared to similar binaries in the circular limit by a factor of (depending on ρ_p0, the component masses, and the sky-position and binary-orientation angles; see Figures 3 and 5)

  • i.  
    ∼1–200 for the mass errors,
  • ii.  
    ∼1–4.5 for the semimajor and semiminor axes of the sky-localization ellipse,
  • iii.  
    ∼1–2 for the distance errors,
  • iv.  
    ∼1–3 for the semimajor and semiminor axes of the error ellipse for the binary orientation.
• 4.  
  For initially highly eccentric BH binaries at D_L = 100 Mpc, the measurement errors for the parameters specific to precessing highly eccentric BH binary sources may be as low as of order (depending on ρ_p0, the component masses, and the sky-position and binary-orientation angles; see Figures 4 and 5)

  • i.  
    10⁻⁵ for the final eccentricity errors at LSO,
  • ii.  
    10⁻⁴ for the initial eccentricity errors,
  • iii.  
    10⁻³ for the initial pericenter distance.

For initially moderately eccentric to low-eccentricity binaries, the parameter measurement errors and S/Ns improve by a smaller amount for low to intermediate ρ_p0 (Figure 8).

Note that the eccentricity errors are remarkably low, which is not surprising given that eccentricity is encoded in several measurable features of the waveform, including the orbital harmonics, the splitting of each frequency harmonic into triplets, the frequency evolution of harmonics (the "chirp"), the frequency evolution of the distance between the spectral triplets, the low-frequency cutoff of the signal at the initial pericenter frequency, and the eccentricity dependence of the LSO where the inspiral transitions into a rapid coalescence.

However, there are several factors that may significantly increase the measurement errors in more typical cases. First, more typical sources are expected to be at much larger distances than 100 Mpc. Assuming crudely that the eccentricity errors scale with D_L, the median measurement errors for a 30 ${M}_{\odot }$–30 ${M}_{\odot }$ precessing highly eccentric BH binary at ∼410 Mpc (similar to GW150914) are expected to be ${\rm{\Delta }}{e}_{\mathrm{LSO}}\sim 8.8\times {10}^{-4}$ ($1.3\times {10}^{-4}$) for the final eccentricity at the LSO if ${\rho }_{{\rm{p}}0}=10$ (20) for the design sensitivity of the aLIGO/AdV/KAGRA instruments. In these cases, the expected median initial eccentricity error is ${\rm{\Delta }}{e}_{0}\sim 8.9\times {10}^{-3}$ (1.2 × 10⁻²), and the median initial pericenter distance is ${\rm{\Delta }}{\rho }_{{\rm{p}}0}\sim 4.4\times {10}^{-2}$ (9.4 × 10⁻²; see Table 2).

Another important simplifying assumption, which may have skewed the errors to lower values, was to ignore higher-order PN corrections that depend on the spin of the merging objects. The spins of the two binary components introduce six additional parameters, which may become partially degenerate with all other parameters, thereby increasing their errors. On the other hand, spin precession breaks degeneracies between the binary orientation and other slow parameters (Lang & Hughes 2006; Kocsis et al. 2007; Chatziioannou et al. 2014). However, the eccentricity-induced orbital harmonics enter at the Newtonian order, and the frequency triplets due to GR precession enter at the low 1PN order. Therefore, the eccentricity-related spectral features are already dominant at the early stages of the inspiral when higher-order PN corrections are negligible. For this reason, the estimated Δe₀ and Δρ_p0 errors are expected to be robust. On the other hand, most of the S/N accumulates at late times for stellar-mass BH binaries, where the high-order PN corrections are significant. We leave an estimate of the parameter estimation errors for spinning binaries to future work.

Finally, an important simplification is the Fisher matrix method itself, which is valid only if the waveform model is a faithful representation of the GW signal and the noise is Gaussian and sufficiently small that the S/N is sufficiently large that the parameter error region is an ellipsoid in parameter space and when the parameter derivate of the waveform is non-vanishing. For smaller S/N, the parameter error region geometry is more complex and the uncertainties are generally higher (see Cornish & Littenberg 2015 and references therein). The parameter estimation errors may also be affected by theoretical uncertainties of the waveform model (Cutler & Vallisneri 2007), which may be especially important for highly eccentric binaries, where the PN expansion is slowly convergent (Kocsis & Levin 2012).

Such low-eccentricity errors may give the aLIGO–AdV–KAGRA GW detector network the capability to distinguish among different astrophysical formation channels. In a companion study (Gondán et al. 2017), we illustrate the expected distribution of eccentricities and other physical parameters for single–single GW capture sources in GNs. Similar studies for other astrophysical formation channels are underway.

Future multiwaveband searches of eccentric inspiral sources with LISA and aLIGO–AdV–KAGRA (Kocsis & Levin 2012; Sesana 2016) have good prospects for even more accurate measurements of the physical parameters well beyond the level reported here. In that case, the GW frequency range is much wider. Since eccentricity decreases due to GW emission, eccentricity may be expected to be much higher at lower frequencies in the LISA band. This leads to a much larger total GW phase shift caused by eccentricity. A better measurement of relativistic precession may more efficiently break degeneracies between mass and other parameters. The modulation caused by the orbit of the instrument around the Sun and Earth's spin can help break degeneracies among the source direction, orientation, and other parameters. Accounting for eccentricity for third-generation Earth-based (e.g., Einstein telescope), deci-Hertz to mHz space-based instruments will be essential (Chen & Amaro-Seoane 2017).

We thank the anonymous referee for constructive comments, which helped improve the quality of the paper. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No. 638435 (GalNUC) and by the Hungarian National Research, Development, and Innovation Office grant NKFIH KH-125675. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607761. The calculations were carried out on the NIIF HPC cluster at the University of Debrecen, Hungary.

The circular limit follows from the limit ${e}_{0}\to 0$ for arbitrary ρ_p0 or ${\rho }_{{\rm{p}}0}\to \infty$ for arbitrary e₀ ≤ 1 (see Figure 5 in O'Leary et al. 2009). In practice, we set e₀ = 10⁻⁴, ρ_p0 = 1000, and omit the parameters from the Fisher matrix that become degenerate or unconstrained in the circular limit: γ_c, e₀, and e_LSO. Thus, the parameters in the circular limit are

$$\begin{eqnarray}{{\boldsymbol{\lambda }}}_{\mathrm{Prec},\mathrm{circ}} & = & \{{t}_{c},{{\rm{\Phi }}}_{c},\mathrm{ln}({D}_{L}),\mathrm{ln}({{ \mathcal M }}_{z}),\\ & & \mathrm{ln}({M}_{\mathrm{tot},z}),{\theta }_{N},{\phi }_{N},{\theta }_{L},{\phi }_{L}\},\end{eqnarray} \tag{ 55 }$$

where M_tot,z arises due to precession even in the circular limit (see Equation (25)). In the ${\rho }_{{\rm{p}}0}\to \infty$ limit, ${\rm{\Delta }}{\gamma }_{0}\to \infty$ and ${\rm{\Delta }}{e}_{0}\to \infty$; however, the Fisher matrix algorithm becomes invalid for Δe_LSO in this limit (Section 6.2).

We have also examined how the measurement errors depend on the total mass of the binary in the circular limit, and qualitatively compared our results to those of previous parameter estimation studies in Section 6.3. We have found an excellent match between our results and the results of previous studies.

Here we describe the measured signal of individual ground-based detectors, and introduce the adopted coordinate system.

We define the Cartesian coordinate system with basis vectors ${\boldsymbol{i}}$, ${\boldsymbol{j}}$, ${\boldsymbol{k}}$ and a spherical coordinate system (θ, ϕ) fixed relative to the center of the Earth, such that ${\boldsymbol{k}}$ and θ = 0 is along the north geographic pole and ϕ = 0 is along the prime meridian. We denote the unit vector pointing from the center of the Earth to the binary's sky position as ${\boldsymbol{N}}$, and define ${\boldsymbol{L}}$ to be the normal vector parallel to the binary's orbital angular momentum,

$$\begin{eqnarray}&&{\boldsymbol{N}}=\sin {\theta }_{N}\cos {\phi }_{N}\,{\boldsymbol{i}}+\sin {\theta }_{N}\sin {\phi }_{N}\,{\boldsymbol{j}}+\cos {\theta }_{N}\,{\boldsymbol{k}},\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&{\boldsymbol{L}}=\sin {\theta }_{L}\cos {\phi }_{L}\,{\boldsymbol{i}}+\sin {\theta }_{L}\sin {\phi }_{L}\,{\boldsymbol{j}}+\cos {\theta }_{L}\,{\boldsymbol{k}}.\end{eqnarray} \tag{ 57 }$$

We denote the unit vectors parallel to the arms of the kth detector as ${\boldsymbol{x}}$_k and ${\boldsymbol{y}}$_k, and set ${\boldsymbol{z}}$_k = ${\boldsymbol{x}}$_k × ${\boldsymbol{y}}$_k. As ${\boldsymbol{x}}$_k and ${\boldsymbol{y}}$_k are parallel to the surface of the Earth for all detectors, ${\boldsymbol{z}}$_k points from the center of the Earth toward the geographical location of the kth detector. Let the coordinates (θ_k, ϕ_k) denote the location of the kth detector, thus the unit vectors along the arms can be expressed as

$$\begin{eqnarray}{{\boldsymbol{x}}}_{k} & = & (\cos {\psi }_{k}\sin {\phi }_{k}-\sin {\psi }_{k}\cos {\phi }_{k}\cos {\theta }_{k}){\boldsymbol{i}}\\ & & +(-\cos {\psi }_{k}\cos {\phi }_{k}-\sin {\psi }_{k}\sin {\phi }_{k}\cos {\theta }_{k}){\boldsymbol{j}}\\ & & +(\sin {\psi }_{k}\sin {\theta }_{k}){\boldsymbol{k}},\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}{{\boldsymbol{y}}}_{k} & = & (-\sin {\psi }_{k}\sin {\phi }_{k}-\cos {\psi }_{k}\cos {\phi }_{k}\cos {\theta }_{k}){\boldsymbol{i}}\\ & & +(\sin {\psi }_{k}\cos {\phi }_{k}-\cos {\psi }_{k}\sin {\phi }_{k}\cos {\theta }_{k}){\boldsymbol{j}}\\ & & +(\cos {\psi }_{k}\sin {\theta }_{k}){\boldsymbol{k}},\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&{{\boldsymbol{z}}}_{k}=\sin {\theta }_{k}\cos {\phi }_{k}\,{\boldsymbol{i}}+\sin {\theta }_{k}\sin {\phi }_{k}\,{\boldsymbol{j}}+\cos {\theta }_{k}\,{\boldsymbol{k}}\,\end{eqnarray} \tag{ 60 }$$

(Creighton & Anderson 2011), where the orientation angle of the $k\mathrm{th}$ detector, ψ_k, is defined in Section 3.

These vectors define the response tensor for the kth detector:

$$\begin{eqnarray}&&{D}_{k}^{{ij}}=\displaystyle \frac{1}{2}({x}_{k}^{i}{x}_{k}^{j}-{y}_{k}^{i}{y}_{k}^{j})\,\end{eqnarray} \tag{ 61 }$$

(Finn & Chernoff 1993), where x_kⁱ and y_kⁱ are the ith Cartesian components of ${\boldsymbol{x}}$_k and ${\boldsymbol{y}}$_k.

We adopt the basis vectors following the conventions of previous studies (Finn & Chernoff 1993; Cutler & Flanagan 1994; Anderson et al. 2001; Dalal et al. 2006; Nissanke et al. 2010),

$$\begin{eqnarray}&&{\boldsymbol{X}}=\displaystyle \frac{{\boldsymbol{N}}\times {\boldsymbol{L}}}{| {\boldsymbol{N}}\times {\boldsymbol{L}}| },\quad {\boldsymbol{Y}}=\displaystyle \frac{{\boldsymbol{X}}\times {\boldsymbol{N}}}{| {\boldsymbol{X}}\times {\boldsymbol{N}}| }\,\end{eqnarray} \tag{ 62 }$$

with preferred polarization basis tensors

$$\begin{eqnarray}&&{e}_{{ij}}^{+}={X}_{i}{X}_{j}-{Y}_{i}{Y}_{j},\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&{e}_{{ij}}^{\times }={X}_{i}{Y}_{j}+{Y}_{i}{X}_{j},\end{eqnarray} \tag{ 64 }$$

where i and j are Cartesian components. Thus, the transverse-traceless metric perturbation describing the GW is written as

$$\begin{eqnarray}&&{h}_{{ij}}={h}_{+}{e}_{{ij}}^{+}+{h}_{\times }{e}_{{ij}}^{\times },\end{eqnarray} \tag{ 65 }$$

where h₊ and h_× are given in Equations (1) and (2).

The response of the kth detector to a GW with frequency f can be given in the time domain by

$$\begin{eqnarray}&&{h}_{k}={e}^{i{\rm{\Delta }}{{\rm{\Phi }}}_{k}}{D}_{k}^{{ij}}{h}_{{ij}}={e}^{i{\rm{\Delta }}{{\rm{\Phi }}}_{k}}({h}_{+}{F}_{+,k}+{h}_{\times }{F}_{\times ,k})\,\end{eqnarray} \tag{ 66 }$$

(Nissanke et al. 2010), where ${{\boldsymbol{r}}}_{k}={R}_{\oplus }{{\boldsymbol{z}}}_{k}$ is the position of the kth detector, the factor $-{\boldsymbol{N}}\cdot {{\boldsymbol{r}}}_{k}$ measures the light travel time between the kth detector and the coordinate origin, and thus the factor ${\rm{\Delta }}{{\rm{\Phi }}}_{k}=-2\pi f{\boldsymbol{N}}\cdot {{\boldsymbol{r}}}_{k}$ measures the phase shift between the kth detector and the coordinate origin. In Equation (66), F_+,k and ${F}_{\times ,k}$ are the antenna factors

$$\begin{eqnarray}&&{F}_{+,k}={e}_{{ij}}^{+}{D}_{k}^{{ij}},\quad {F}_{\times ,k}={e}_{{ij}}^{\times }{D}_{k}^{{ij}}.\end{eqnarray} \tag{ 67 }$$

In our calculations, Earth is taken to be a sphere with a radius of ${R}_{\oplus }=6,370\,\mathrm{km}$.

If the time that the GW signal spends in the detectors' sensitive frequency band is negligible compared to the rotation period of the Earth, then the measured waveform in the frequency domain for the kth detector is

$$\begin{eqnarray}&&{\tilde{h}}_{k}(f)=[{F}_{+,k}{\tilde{h}}_{+}(f)+{F}_{\times ,k}{\tilde{h}}_{\times }(f)]{e}^{-2\pi {if}{\boldsymbol{N}}\cdot {{\boldsymbol{r}}}_{k}},\end{eqnarray} \tag{ 68 }$$

where ${\tilde{h}}_{+}(f)$ and ${\tilde{h}}_{\times }(f)$ are the Fourier-transformed expressions of h₊ and h_× at Earth's center, and ${F}_{+,k}$ and ${F}_{\times ,k}$ are given by the (practically time-independent) orientation of the detectors shown in Table 1.

Similarly, using the frequency harmonic triplets for eccentric precessing inspiraling binaries in the stationary phase approximation, the measured waveform for the kth detector is

$$\begin{eqnarray}&&{\tilde{h}}_{k}({\boldsymbol{f}})={F}_{+,k}{\tilde{h}}_{+,k}({\boldsymbol{f}})+{F}_{\times ,k}{\tilde{h}}_{\times ,k}({\boldsymbol{f}}),\end{eqnarray} \tag{ 69 }$$

where ${\tilde{h}}_{+,k}({\boldsymbol{f}})$ and ${\tilde{h}}_{\times ,k}({\boldsymbol{f}})$ can be derived from ${\tilde{h}}_{+}({\boldsymbol{f}})$ and ${\tilde{h}}_{\times }({\boldsymbol{f}})$ by multiplying each term of ${\boldsymbol{f}}$ with the phase shift factors ${e}^{-2\pi {{if}}_{n}{\boldsymbol{N}}\cdot {{\boldsymbol{r}}}_{k}}$ and ${e}^{-2\pi {{if}}_{n}^{\pm }{\boldsymbol{N}}\cdot {{\boldsymbol{r}}}_{k}}$ for each harmonic, respectively. More specifically, the measured signal's Fourier phase in Equations (35) and (36) in the kth detector is shifted respectively according to

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{n,k}={{\rm{\Psi }}}_{n}-2\pi {{if}}_{n}{\boldsymbol{N}}\cdot {{\boldsymbol{r}}}_{k},\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{n,k}^{\pm }={{\rm{\Psi }}}_{n}^{\pm }-2\pi {{if}}_{n}^{\pm }{\boldsymbol{N}}\cdot {{\boldsymbol{r}}}_{k}.\end{eqnarray} \tag{ 71 }$$

In this section, we derive the time evolution of different harmonics in the detectors' sensitive frequency band, and for each orbital harmonic, we determine the eccentricity at which the signal enters the detectors' sensitive frequency band. These formulae will be utilized in Appendix D.

The time-dependent GW signal of a precessing eccentric BH binary as measured by the kth detector can be given as

$$\begin{eqnarray}{h}_{k}(t) & = & {h}_{+}(t){F}_{+,k}[{\alpha }_{N}(t),{\beta }_{N},{\alpha }_{L},{\beta }_{L}]\\ & & +{h}_{\times }(t){F}_{\times ,k}[{\alpha }_{N}(t),{\beta }_{N},{\alpha }_{L},{\beta }_{L}],\end{eqnarray} \tag{ 72 }$$

where ${h}_{+}(t)$ and ${h}_{\times }(t)$ are given in Equations (1) and (2), and F_+,k and ${F}_{\times ,k}$ are quantified by Equation (67). We ignore spins in this study, and therefore the angular momentum vector direction $({\alpha }_{L},{\beta }_{L})$ is conserved during the eccentric inspiral (Cutler & Flanagan 1994). The polar angle of the source α_N relative to the detector depends on the rotation phase of the Earth during the day. We ignore the Earth's rotation, since the total duration of an eccentric inspiral from e = 0.9 to merger is of order $[{({\rho }_{{\rm{p}}0}/40)}^{4}{(4\eta )}^{-1}{M}_{\mathrm{tot}}/(20\,{M}_{\odot })\,\min ]$ as shown in Figure 3 in O'Leary et al. (2009).¹⁶

The waveform of an eccentric binary in the stationary phase approximation is a sum over harmonics n for each component of the frequency triplet (${f}_{n},{f}_{n}^{+},{f}_{n}^{-}$) with different reference times $({t}_{n},{t}_{n}^{+},{t}_{n}^{-})$ (see Section 2.2). During the evolution, e(t) (and ${\rho }_{p}(t)$) shrinks strictly monotonically in time (Peters 1964); therefore, its inverse function t(e) is well-defined and determines t_n, ${t}_{n}^{+}$, and ${t}_{n}^{-}$. For inspiraling circular binaries, time t can be expressed using the frequency of the emitted GW signal $f=2\nu$ as

$$\begin{eqnarray}&&t(f)={t}_{c}-{\int }_{f}^{\infty }\displaystyle \frac{{df}^{\prime} }{\dot{f}^{\prime} }={t}_{c}-5{(8\pi f)}^{-8/3}{{ \mathcal M }}^{-5/3}\end{eqnarray} \tag{ 73 }$$

(see Equations (2.13) and (2.19) in Cutler & Flanagan 1994 for details), where the constant of integration, t_c, is defined by the requirement that $t\to {t}_{c}$ as $f\to \infty$. We generalize Equation (73) for eccentric inspirals by changing the integration variable from f to e in Equation (73),

$$\begin{eqnarray}&&t(e)={t}_{c}+{\int }_{0}^{e}\displaystyle \frac{{de}^{\prime} }{\dot{e}(e^{\prime} )}={t}_{c}-\tau {I}_{t}(e),\end{eqnarray} \tag{ 74 }$$

where $\dot{e}$ is given by Equation (15), and we introduced

$$\begin{eqnarray}\tau & = & \displaystyle \frac{15}{304}{{ \mathcal M }}^{-5/3}{(2\pi {c}_{0})}^{-8/3}\\ & & =\displaystyle \frac{15}{304}\displaystyle \frac{{M}_{\mathrm{tot}}^{8/3}}{{{ \mathcal M }}^{5/3}}\displaystyle \frac{{[(1-{e}_{\mathrm{LSO}}){\rho }_{\mathrm{pLSO}}({e}_{\mathrm{LSO}})]}^{4}}{H{({e}_{\mathrm{LSO}})}^{8/3}},\end{eqnarray} \tag{ 75 }$$

and substituted ${c}_{0}$ using Equation (22). Here, ${\rho }_{\mathrm{pLSO}}({e}_{\mathrm{LSO}})$ and $H({e}_{\mathrm{LSO}})$ are given by Equations (21) and (19), and I_t(e) in Equation (74) is of the form

$$\begin{eqnarray}&&{I}_{t}(e)={\int }_{0}^{e}\displaystyle \frac{{x}^{29/19}{\left(1+\tfrac{121}{304}{x}^{2}\right)}^{\tfrac{1181}{2299}}}{{(1-{x}^{2})}^{3/2}}\,{dx}\end{eqnarray} \tag{ 76 }$$

(see Mikóczi et al. 2012 for an analytic result). Using Equation (74), we obtain the total duration of the nth harmonic in the detector's sensitive frequency band

$$\begin{eqnarray}&&{T}_{n}=t({e}_{\min ,n})-t({e}_{\max ,n})=[{I}_{1}({e}_{\max ,n})-{I}_{1}({e}_{\min ,n})]\tau ,\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&{T}_{n}^{+}=t({e}_{\min ,n}^{+})-t({e}_{\max ,n}^{+})=[{I}_{1}({e}_{\max ,n}^{+})-{I}_{1}({e}_{\min ,n}^{+})]\tau ,\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}&&{T}_{n}^{-}=t({e}_{\min ,n}^{-})-t({e}_{\max ,n}^{-})=[{I}_{1}({e}_{\max ,n}^{-})-{I}_{1}({e}_{\min ,n}^{-})]\tau .\end{eqnarray} \tag{ 79 }$$

Here, ${e}_{\min ,n}$ (${e}_{\min ,n}^{+},{e}_{\min ,n}^{-}$) refers to the eccentricity at which the harmonic f_n (${f}_{n}^{+},\,{f}_{n}^{-}$) reaches the LSO or when it exists the detectable highest frequency for the given detector, and ${e}_{\max ,n}$ (${e}_{\max ,n}^{+},{e}_{\max ,n}^{-}$) refers to the eccentricity at which the signal related to f_n (${f}_{n}^{+},\,{f}_{n}^{-}$) first enters the detector's sensitive frequency band or when it forms within the band. Thus,

$$\begin{eqnarray}&&{e}_{\min ,n}=\max ({e}_{\mathrm{LSO}},{e}_{\det ,n}^{\min }),\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&{e}_{\max ,n}=\min ({e}_{0},{e}_{\det ,n}^{\max }),\end{eqnarray} \tag{ 81 }$$

where

$$\begin{eqnarray}&&{e}_{\det ,n}^{\min }={\nu }_{n}^{-1}({f}_{\det ,\max }),\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&{e}_{\det ,n}^{\max }={\nu }_{n}^{-1}({f}_{\det ,\min }),\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&{\nu }_{n}(e)=n\nu (e).\end{eqnarray} \tag{ 84 }$$

Here, $\nu (e)$ is given analytically by Equation (18), ${\nu }_{n}^{-1}(\cdot )$ denotes the inverse function of ${\nu }_{n}(e)$, ${f}_{\det ,\min }$ and ${f}_{\det ,\max }$ are the lower and upper limits of the detector's sensitive frequency band (typically $10\,\mathrm{Hz}$ and 10⁴ Hz, respectively), and e_LSO is determined from Equation (21). Similarly, we define the parameters corresponding to ${f}_{n}^{+}$ and ${f}_{n}^{-}$ as

$$\begin{eqnarray}&&{e}_{\min ,n}^{+}=\max ({e}_{\mathrm{LSO}},{e}_{\det ,n+}^{\min }),\end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray}&&{e}_{\max ,n}^{+}=\min ({e}_{0},{e}_{\det ,n+}^{\max }),\end{eqnarray} \tag{ 86 }$$

$$\begin{eqnarray}&&{e}_{\min ,n}^{-}=\max ({e}_{\mathrm{LSO}},{e}_{\det ,n-}^{\min }),\end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray}&&{e}_{\max ,n}^{-}=\min ({e}_{0},{e}_{\det ,n-}^{\max }),\end{eqnarray} \tag{ 88 }$$

where

$$\begin{eqnarray}&&{e}_{\det ,n+}^{\min }={\nu }_{n+}^{-1}({f}_{\det ,\max }),\end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray}&&{e}_{\det ,n+}^{\max }={\nu }_{n+}^{-1}({f}_{\det ,\min }),\end{eqnarray} \tag{ 90 }$$

$$\begin{eqnarray}&&{e}_{\det ,n-}^{\min }={\nu }_{n-}^{-1}({f}_{\det ,\max }),\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}&&{e}_{\det ,n-}^{\max }={\nu }_{n-}^{-1}({f}_{\det ,\min }),\end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray}&&{\nu }_{n\pm }(e)=n\nu (e)\pm \displaystyle \frac{\dot{\gamma }(e)}{\pi },\end{eqnarray} \tag{ 93 }$$

where $\dot{\gamma }(e)$ is given by Equation (25), and ${\nu }_{n\pm }^{-1}(\cdot )$ is the inverse function of ${\nu }_{n\pm }(e)$ given by Equation (93).

In practice, the second term is negligible in Equations (77)–(79). We find that ${T}_{n}^{+}$ and ${T}_{n}^{-}$ are within ≤20% of T_n for any fixed n. The top panel of Figure 9 shows the total time the GW signal spends in an aLIGO-type detector's sensitive frequency band for different harmonics. Higher harmonics enter the aLIGO band earlier and that depending on ρ_p0, the first 10 orbital harmonics spend between seconds to minutes in the detector's sensitive frequency band for a 30 ${M}_{\odot }$–30 ${M}_{\odot }$ precessing highly eccentric BH binary.

The bottom panel of Figure 9 shows the fraction of the squared S/N that accumulates in different orbital harmonics for various ρ_p0 for aLIGO. For high ρ_p0, the signal effectively circularizes by the time it enters the detector's sensitive frequency band and the n = 2 harmonic dominates. However, the contribution of $n\ne 2$ is significant for ρ_p0 ≲ 20 for a 30 ${M}_{\odot }$–30 ${M}_{\odot }$ precessing highly eccentric BH binary.

Zoom In Zoom Out Reset image size

In this section, we derive numerically efficient formulae to calculate the S/N and the Fisher matrix for individual detectors. We first ignore pericenter precession, then extend the calculations for precessing eccentric sources.

D.1. Signal-to-noise Ratio

D.1.1. Eccentric Inspirals Without Precession

The NoPrec signal measured by a detector at position ${\boldsymbol{r}}$ is given in Fourier space from Equations (31) and (32) as

$$\begin{eqnarray}&&{\tilde{h}}_{\mathrm{NoPrec}}=\displaystyle \sum _{n=1}^{\infty }{L}_{n}(e,{f}_{n}){{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e){e}^{i{{\rm{\Psi }}}_{n}(e,{f}_{n})},\end{eqnarray} \tag{ 94 }$$

where Ψ_n is the Fourier phase at the origin of the coordinate system set to the Earth's center, given by Equation (35), $\dot{\gamma }\equiv 0$, γ ≡ γ_c, ${{\rm{\Theta }}}_{{\rm{H}}}(\cdot )$ denotes the Heaviside function, which is zero and unity for negative and positive arguments,¹⁷ respectively, and

$$\begin{eqnarray}{L}_{n}(e,{f}_{n}) & = & -\left[\displaystyle \frac{{A}_{n}{\sin }^{2}{\rm{\Theta }}}{4}+\displaystyle \frac{(1+{\cos }^{2}{\rm{\Theta }})}{4}\right.\\ & & \times ({B}_{n}^{+}{e}^{2i{\gamma }_{c}}-{B}_{n}^{-}{e}^{-2i{\gamma }_{c}})]\times {h}_{0}{F}_{+}{{\rm{\Lambda }}}_{n}{e}^{i({\rm{\Delta }}{{\rm{\Phi }}}_{n}-\pi /4)}\\ & & -\displaystyle \frac{{{ih}}_{0}{F}_{\times }{{\rm{\Lambda }}}_{n}\cos {\rm{\Theta }}}{2}({B}_{n}^{+}{e}^{2i{\gamma }_{c}}+{B}_{n}^{-}{e}^{-2i{\gamma }_{c}})\\ & & \times {e}^{i({\rm{\Delta }}{{\rm{\Phi }}}_{n}+\pi /4)},\end{eqnarray} \tag{ 97 }$$

where h₀ and Λ_n are given by Equations (33) and (34), A_n and ${B}_{n}^{\pm }$ are given by Equations (5) and (7), γ_c specifies the argument of pericenter, which is assumed to be fixed here, and F₊ and ${F}_{\times }$ are the antenna factors given by Equation (67). The factor ${\rm{\Delta }}{{\rm{\Phi }}}_{n}=-2\pi {f}_{n}{\boldsymbol{N}}\cdot {\boldsymbol{r}}$ gives the phase shift of the measured signal between the position of the detector ${\boldsymbol{r}}$ and the origin of the coordinate system for the nth harmonic (Appendix B). L_n depends on f_n implicitly through ΔΦ_n and h₀.

In Equation (94), ${{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e)$ accounts for the start of the waveform when the binary forms with initial eccentricity¹⁸ e₀. Along the same lines, a similar term ${{\rm{\Theta }}}_{{\rm{H}}}(e-{e}_{\mathrm{LSO}})$ could be incorporated to account for the end of the eccentric inspiral, where the waveform transitions to a plunge and ringdown phase. However, we conservatively do not account for such a term, since the waveform near the end of the inspiral is sensitive to higher-order PN corrections, which are not known and ignored here (Kocsis & Levin 2012; Loutrel & Yunes 2017). Nevertheless, the inspiral rate is sensitive to e_LSO in Equations (18), (19), and (22), which affects L_n and Ψ_n.

For each detector, the square of the S/N for the NoPrec waveform, ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{2}$, can be obtained by substituting ${\tilde{h}}_{\mathrm{NoPrec}}$ into Equation (42). We find that the product of sums in ${\tilde{h}}_{\mathrm{NoPrec}}\,{\tilde{h}}_{\mathrm{NoPrec}}^{* }$ is dominated by the elements such that¹⁹

$$\begin{eqnarray}&&{({\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}})}^{2}\approx 4\displaystyle \sum _{n=1}^{\infty }{\int }_{{f}_{\min ,n}}^{{f}_{\max ,n}}\displaystyle \frac{| {L}_{n}(e({f}_{n}),{f}_{n}){| }^{2}}{{S}_{n}({f}_{n})}{{df}}_{n},\end{eqnarray} \tag{ 98 }$$

where ${f}_{\min ,n}$ is the frequency at which the $n\mathrm{th}$ harmonic first enters the detector's sensitive frequency band or when it forms in the band, and similarly ${f}_{\max ,n}$ is the frequency at which the signal exits the detector's sensitive frequency band or when it reaches the LSO,

$$\begin{eqnarray}&&{f}_{\max ,n}=\min ({\nu }_{n}({e}_{\mathrm{LSO}}),{f}_{\det ,\max }),\end{eqnarray} \tag{ 99 }$$

$$\begin{eqnarray}&&{f}_{\min ,n}=\max ({\nu }_{n}({e}_{0}),{f}_{\det ,\min }).\end{eqnarray} \tag{ 100 }$$

Computationally, it is practical to change the integration variable from f_n to e as

$$\begin{eqnarray}&&{{df}}_{n}=n\left|\displaystyle \frac{d\nu }{{de}}\right|\,{de},\end{eqnarray} \tag{ 101 }$$

thus Equation (98) can be rewritten generally as

$$\begin{eqnarray}&&{({\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}})}^{2}\approx 4\displaystyle \sum _{n=1}^{{n}_{\max }({e}_{0})}{\int }_{{e}_{\min ,n}}^{{e}_{\max ,n}}\displaystyle \frac{n| {L}_{n}(e){| }^{2}}{{S}_{n}(n\nu (e))}\left|\displaystyle \frac{d\nu }{{de}}\right|\,{de},\end{eqnarray} \tag{ 102 }$$

where $\nu (e)$ is given analytically by Equation (18), and L(e) may be obtained from Equation (97) by substituting ${f}_{n}=n\nu (e)$. The integration bounds ${e}_{\min ,n}$ and ${e}_{\max ,n}$ are given by Equations (80) and (81). We truncate the calculation beyond a maximum spectral harmonic ${n}_{\max }(e)$ defined in Equation (37).

${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$ accurately recovers Equation (98) generally for any value of the initial eccentricity in the range $0\lt {e}_{0}\lt 1$; however, the number of considered harmonics ${n}_{\max }({e}_{0})$ increases rapidly for high e₀ (Equation (37)), and ${n}_{\max }({e}_{0})\to \infty$ in the limit ${e}_{0}\to 1$. Therefore, ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$ is computationally efficient for low to moderate initial eccentricities (${e}_{0}\lesssim 0.8$), and it is inefficient for higher e₀. In order to make ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$ computationally efficient for high initial eccentricities, we reverse the order of the sum and the integral in Equation (102) and truncate the sum over harmonics at ${n}_{\max }(e)$ (O'Leary et al. 2009).²⁰ Thus, we get

$$\begin{eqnarray}&&{({\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{High}})}^{2}\approx 4{\int }_{{e}_{\min ,n}}^{{e}_{\max ,n}}\displaystyle \sum _{n=1}^{{n}_{\max }(e)}\displaystyle \frac{n| {L}_{n}(e){| }^{2}}{{S}_{n}(n\nu (e))}\left|\displaystyle \frac{d\nu }{{de}}\right|\,{de}.\end{eqnarray} \tag{ 103 }$$

We use the above introduced trick to derive computationally efficient formulae in the high initial eccentricity limit for Fisher matrix elements in the precession-free case (Appendix D.2.1) and for the S/N and the Fisher matrix elements in the precessing case (Appendices D.1.2 and D.2.2).

D.1.2. Eccentric Inspirals with Precession

We derive the S/N of the precessing model in this section. The Fourier-transformed waveform given by Equation (69) can be rewritten as

$$\begin{eqnarray}{\tilde{h}}_{\mathrm{Prec}} & = & \displaystyle \sum _{n=1}^{\infty }{K}_{n}(e,{f}_{n}){{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e){e}^{i{{\rm{\Psi }}}_{n}(e,{f}_{n})}\\ & & +\displaystyle \sum _{n=1}^{\infty }{K}_{n}^{+}(e,{f}_{n}^{+}){{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e){e}^{i{{\rm{\Psi }}}_{n}^{+}(e,{f}_{n}^{+})}\\ & & +\displaystyle \sum _{n=1}^{\infty }{K}_{n}^{-}(e,{f}_{n}^{-}){{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e){e}^{i{{\rm{\Psi }}}_{n}^{-}(e,{f}_{n}^{-})},\end{eqnarray} \tag{ 104 }$$

where the terms K_n, ${K}_{n}^{+}$, and ${K}_{n}^{-}$ are defined as

$$\begin{eqnarray}&&{K}_{n}(e,{f}_{n})=-\displaystyle \frac{{h}_{0}{\sin }^{2}{\rm{\Theta }}}{4}{A}_{n}{{\rm{\Lambda }}}_{n}{F}_{+}{e}^{i({\rm{\Delta }}{{\rm{\Phi }}}_{n}-\pi /4)},\end{eqnarray} \tag{ 105 }$$

$$\begin{eqnarray}{K}_{n}^{+}(e,{f}_{n}^{+}) & = & -\displaystyle \frac{{h}_{0}\cos {\rm{\Theta }}}{2}{B}_{n}^{-}{{\rm{\Lambda }}}_{n}^{+}{F}_{\times }{e}^{i({\rm{\Delta }}{{\rm{\Phi }}}_{n}^{+}+\pi /4)}\\ & & +\displaystyle \frac{{h}_{0}(1+{\cos }^{2}{\rm{\Theta }})}{4}{B}_{n}^{-}{{\rm{\Lambda }}}_{n}^{+}{F}_{+}{e}^{i({\rm{\Delta }}{{\rm{\Phi }}}_{n}^{+}-\pi /4)},\end{eqnarray} \tag{ 106 }$$

$$\begin{eqnarray}{K}_{n}^{-}(e,{f}_{n}^{-}) & = & -\displaystyle \frac{{h}_{0}\cos {\rm{\Theta }}}{2}{B}_{n}^{+}{{\rm{\Lambda }}}_{n}^{-}{F}_{\times }{e}^{i({\rm{\Delta }}{{\rm{\Phi }}}_{n}^{-}+\pi /4)}\\ & & -\displaystyle \frac{{h}_{0}}{4}{B}_{n}^{+}{{\rm{\Lambda }}}_{n}^{-}{F}_{+}(1+{\cos }^{2}{\rm{\Theta }}){e}^{i({\rm{\Delta }}{{\rm{\Phi }}}_{n}^{-}-\pi /4)}.\end{eqnarray} \tag{ 107 }$$

The terms K_n and ${K}_{n}^{\pm }$ depend on ν through h₀, which are expressed with ${f}_{n},{f}_{n}^{\pm }$ using Equations (29) and (30). Furthermore, these equations depend on f_n and ${f}_{n}^{\pm }$ through ${\rm{\Delta }}{{\rm{\Phi }}}_{n}$ and ${\rm{\Delta }}{{\rm{\Phi }}}_{n}^{\pm }$ (e.g., ${\rm{\Delta }}{{\rm{\Phi }}}_{n}=-2\pi {f}_{n}{\boldsymbol{N}}\cdot {\boldsymbol{r}})$ and ${\rm{\Delta }}{{\rm{\Phi }}}_{n}^{\pm }=-2\pi {f}_{n}^{\pm }{\boldsymbol{N}}\cdot {\boldsymbol{r}}$).

Next, we substitute this waveform ${\tilde{h}}_{\mathrm{Prec}}$ into Equation (42). Similarly to that of the NoPrec signal (Equation (94)), the cross-terms in the product of sums in ${\tilde{h}}_{\mathrm{Prec}}\,{\tilde{h}}_{\mathrm{Prec}}^{* }$ have negligible contributions to ${({\rm{S}}/{{\rm{N}}}_{\mathrm{Prec}})}^{2}$, and so

$$\begin{eqnarray}{({\rm{S}}/{{\rm{N}}}_{\mathrm{Prec}})}^{2} & \approx & 4\displaystyle \sum _{n=1}^{\infty }{\displaystyle \int }_{{f}_{\min ,n}}^{{f}_{\max ,n}}\displaystyle \frac{| {K}_{n}(e({f}_{n}),{f}_{n}){| }^{2}}{{S}_{h}({f}_{n})}\,{{df}}_{n}\\ & & +4\displaystyle \sum _{n=1}^{\infty }{\displaystyle \int }_{{f}_{\min ,n}^{+}}^{{f}_{\max ,n}^{+}}\displaystyle \frac{| {K}_{n}^{+}(e({f}_{n}^{+}),{f}_{n}^{+}){| }^{2}}{{S}_{h}({f}_{n}^{+})}\,{{df}}_{n}^{+}\\ & & +4\displaystyle \sum _{n=1}^{\infty }{\displaystyle \int }_{{f}_{\min ,n}^{-}}^{{f}_{\max ,n}^{-}}\displaystyle \frac{| {K}_{n}^{-}(e({f}_{n}^{-}),{f}_{n}^{-}){| }^{2}}{{S}_{h}({f}_{n}^{-})}\,{{df}}_{n}^{-},\end{eqnarray} \tag{ 108 }$$

where ${f}_{\max ,n}$ and ${f}_{\min ,n}$ are defined in Equations (99) and (100). The integration bounds for the integrals over ${f}_{n}^{\pm }$ are defined similarly to ${f}_{\max ,n}$ and ${f}_{\min ,n}$ in Equations (99) and (100),

$$\begin{eqnarray}&&{f}_{\max ,n}^{\pm }\,=\,\min ({\nu }_{n\pm }({e}_{\mathrm{LSO}}),{f}_{\det ,\max }),\end{eqnarray} \tag{ 109 }$$

$$\begin{eqnarray}&&{f}_{\min ,n}^{\pm }\,=\,\max ({\nu }_{n\pm }({e}_{0}),{f}_{\det ,\min }),\end{eqnarray} \tag{ 110 }$$

where ${\nu }_{n\pm }$ is defined in Equation (93).

Next, we change the integration variables from f_n to e using Equation (101) and similarly from ${f}_{n}^{\pm }$ to e using Equation (30) as

$$\begin{eqnarray}&&{{df}}_{n}^{\pm }=\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|\left|\displaystyle \frac{d\nu }{{de}}\right|\,{de}.\end{eqnarray} \tag{ 111 }$$

Here, $d\dot{\gamma }/d\nu$ is given by Equation (25) as

$$\begin{eqnarray}\displaystyle \frac{d\dot{\gamma }}{d\nu } & = & \displaystyle \frac{\partial \dot{\gamma }(\nu ,e)}{\partial \nu }+\displaystyle \frac{\partial \dot{\gamma }(\nu ,e)}{\partial e}\displaystyle \frac{1}{d\nu /{de}}\\ & = & \left(\displaystyle \frac{5}{3\nu }-\displaystyle \frac{2e}{1-{e}^{2}}\displaystyle \frac{1}{d\nu /{de}}\right)\dot{\gamma }.\end{eqnarray} \tag{ 112 }$$

After truncating the sum over the harmonics to the relevant range as in Equation (102), Equation (108) can be written as

$$\begin{eqnarray}&&{({\rm{S}}/{{\rm{N}}}_{\mathrm{Prec}}^{\mathrm{Gen}})}^{2}\approx 4\displaystyle \sum _{n=1}^{{n}_{\max }({e}_{0})}{\displaystyle \int }_{{e}_{\min ,n}}^{{e}_{\max ,n}}\displaystyle \frac{n| {K}_{n}(e){| }^{2}}{{S}_{n}({\nu }_{n}(e))}\left|\displaystyle \frac{d\nu }{{de}}\right|{de}\\ &&\quad +\,4\displaystyle \sum _{n=1}^{{n}_{\max }({e}_{0})}{\displaystyle \int }_{{e}_{\min ,n}^{+}}^{{e}_{\max ,n}^{+}}\displaystyle \frac{| {K}_{n}^{+}(e){| }^{2}}{{S}_{n}({\nu }_{n+}(e))}\left|n+\displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|\left|\displaystyle \frac{d\nu }{{de}}\right|{de}\\ &&\quad +\,4\displaystyle \sum _{n=1}^{{n}_{\max }({e}_{0})}{\displaystyle \int }_{{e}_{\min ,n}^{-}}^{{e}_{\max ,n}^{-}}\displaystyle \frac{| {K}_{n}^{-}(e){| }^{2}}{{S}_{n}({\nu }_{n-}(e))}\left|n-\displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|\left|\displaystyle \frac{d\nu }{{de}}\right|{de},\end{eqnarray} \tag{ 113 }$$

where ${K}_{n}(e)\equiv {K}_{n}(e,{\nu }_{n}(e))$ and ${K}_{n}^{\pm }(e)\equiv {K}_{n}^{\pm }(e,{\nu }_{n\pm }(e))$. The integration bounds are given by Equations (85) and (88).

Similarly to ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$, ${\rm{S}}/{{\rm{N}}}_{\mathrm{Prec}}^{\mathrm{Gen}}$ is computationally efficient only for low to moderate initial eccentricities (${e}_{0}\lesssim 0.8$). For high initial eccentricities, the computationally efficient form of ${\rm{S}}/{{\rm{N}}}_{\mathrm{Prec}}^{\mathrm{Gen}}$, ${\rm{S}}/{{\rm{N}}}_{\mathrm{Prec}}^{\mathrm{High}}$, can be given by reversing the order of the sum and the integral as

$$\begin{eqnarray}&&{({\rm{S}}/{{\rm{N}}}_{\mathrm{Prec}}^{\mathrm{High}})}^{2}\approx 4{\displaystyle \int }_{{e}_{\min ,n}}^{{e}_{\max ,n}}\displaystyle \sum _{n=1}^{{n}_{\max }(e)}\displaystyle \frac{n| {K}_{n}(e){| }^{2}}{{S}_{n}({\nu }_{n}(e))}\left|\displaystyle \frac{d\nu }{{de}}\right|\,{de}\\ &&\quad +\,4{\displaystyle \int }_{{e}_{\min ,n}^{+}}^{{e}_{\max ,n}^{+}}\displaystyle \sum _{n=1}^{{n}_{\max }(e)}\displaystyle \frac{| {K}_{n}^{+}(e){| }^{2}}{{S}_{n}({\nu }_{n+}(e))}\left|n+\displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|\left|\displaystyle \frac{d\nu }{{de}}\right|{de}\\ &&\quad +\,4{\displaystyle \int }_{{e}_{\min ,n}^{-}}^{{e}_{\max ,n}^{-}}\displaystyle \sum _{n=1}^{{n}_{\max }(e)}\displaystyle \frac{| {K}_{n}^{-}(e){| }^{2}}{{S}_{n}({\nu }_{n-}(e))}\left|n-\displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|\left|\displaystyle \frac{d\nu }{{de}}\right|\,{de}.\end{eqnarray} \tag{ 114 }$$

D.2. Fisher Matrix

Due to the similarity of the equations defining the S/N (see Equation (42)) and the Fisher matrix (see Equation (44)), we may follow the same procedure to derive numerically efficient formulae for the Fisher matrix in the limit of high initial eccentricity. Similarly to Appendix D.1, we start the analysis with the NoPrec model and then generalize the calculation to the Prec model, which accounts for the precessing case.

D.2.1. Eccentric Inspirals Without Precession

Let us substitute Equation (94) into Equation (44). Similarly to the product of the sums of orbital harmonics ${\tilde{h}}_{\mathrm{NoPrec}}{\tilde{h}}_{\mathrm{NoPrec}}^{* }$, we find numerically that the cross-terms in ${\partial }_{j}{\tilde{h}}_{\mathrm{NoPrec}}\,{\partial }_{k}{\tilde{h}}_{\mathrm{NoPrec}}^{* }$ with $j\ne k$ have a negligible contribution to ${{\rm{\Gamma }}}_{{jk}}$. Thus, we find that the stationary phase approximation is applicable if we drop the cross-terms, and thus in Equation (44) we may use

$$\begin{eqnarray}&&{\partial }_{j}{\tilde{h}}_{\mathrm{NoPrec}}\,{\partial }_{k}{\tilde{h}}_{\mathrm{NoPrec}}^{* }=\displaystyle \sum _{n=1}^{\infty }{\tilde{L}}_{n,j}(e)\,{\tilde{L}}_{n,k}^{* }(e){{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e)\\ &&\qquad +\,{\delta }_{{e}_{0},j}\displaystyle \sum _{n=1}^{\infty }{\tilde{L}}_{n}(e)\,{\tilde{L}}_{n,k}^{* }(e)\delta (e-{e}_{0})\\ &&\qquad +\,{\delta }_{{e}_{0},k}\displaystyle \sum _{n=1}^{\infty }{\tilde{L}}_{n,j}(e)\,{\tilde{L}}_{n}^{* }(e)\delta (e-{e}_{0})\\ &&\qquad +\,{\delta }_{{e}_{0},k}{\delta }_{{e}_{0},j}\displaystyle \sum _{n=1}^{\infty }| {L}_{n}(e){| }^{2}{[{\partial }_{{e}_{0}}{{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e)]}^{2}.\end{eqnarray} \tag{ 115 }$$

Here,

$$\begin{eqnarray}&&{\tilde{L}}_{n}(e)={L}_{n}({f}_{n}(e),e)\,{e}^{i{{\rm{\Psi }}}_{n}({f}_{n}(e),e)},\end{eqnarray} \tag{ 116 }$$

and

$$\begin{eqnarray}&&{\tilde{L}}_{n,j}(e)={\partial }_{j}[{L}_{n}({f}_{n}){e}^{i{{\rm{\Psi }}}_{n}({f}_{n})}](e),\end{eqnarray} \tag{ 117 }$$

where

$$\begin{eqnarray}&&{L}_{n}({f}_{n})\equiv {L}_{n}(e({f}_{n}),{f}_{n}),\quad {{\rm{\Psi }}}_{n}({f}_{n})\equiv {{\rm{\Psi }}}_{n}(e({f}_{n}),{f}_{n}),\end{eqnarray} \tag{ 118 }$$

and ${L}_{n}(e,{f}_{n})$ and ${{\rm{\Psi }}}_{n}(e,{f}_{n})$ are given by Equations (97) and (35). We first differentiate the expressions in the bracket $[\,]$ in Equation (117) with respect to λ_j, then change the variable from f_n back to e. Note that for all n and e, ${\tilde{L}}_{n}(e)$ is independent of e₀, and so ${\tilde{L}}_{n,j}(e)=0$ for ${\lambda }_{j}={e}_{0}$. In Equation (115), ${\delta }_{a,b}$ in the second, third, and fourth terms denote the Kronecker δ, defined to be unity if a = b and zero otherwise. In Equation (115), the second and third terms arise due to the Heaviside function in the waveform in Equation (94), which represents the start of the waveform with eccentricity e₀. The e₀ derivative of this function is $\delta ({e}_{0}-e)$, which denotes the Dirac δ function. Note that we use a smoothed version of ${{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e)$ over a scale $\delta {e}_{0}$, which is given in Equation (95), whose derivative is approximately²¹

$$\begin{eqnarray}{[{\partial }_{{e}_{0}}{{\rm{\Theta }}}_{{\rm{H}}}({e}_{0}-e)]}^{2} & \approx & \displaystyle \frac{1}{{(\delta {e}_{0})}^{2}}\,\mathrm{if}\,\,{e}_{0}-\delta {e}_{0}\leqslant e\lt {e}_{0}\\ & \approx & \displaystyle \frac{\delta ({e}_{0}-e)}{\delta {e}_{0}}.\end{eqnarray} \tag{ 119 }$$

To avoid confusion, note that the numerator denotes the Dirac δ function, which has a unit integral over e ≈ e₀, and $\delta {e}_{0}$ in the denominator is the quantity given by Equation (96).

Furthermore, we note that

$$\begin{eqnarray}&&{\mathfrak{R}}({\tilde{L}}_{n}({e}_{0}){\tilde{L}}_{n,k}^{* }({e}_{0}))=\displaystyle \frac{1}{2}{\partial }_{k}{| {L}_{n}({e}_{0})| }^{2}\end{eqnarray} \tag{ 120 }$$

in Equation (115). In these equations, f_n enters when substituting ${f}_{n}/n$ for ν. By substituting Equation (115) into Equation (44), changing the integration variable from f_n to e respectively for each harmonic using²² ${f}_{n}=n\nu (e)$ and Equations (18) and (22), and truncating the sum over the harmonics to the relevant range, the Fisher matrix becomes²³

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{jk}}^{\mathrm{NoPrec},\mathrm{Gen}}\approx 4\displaystyle \sum _{n=1}^{{n}_{\max }({e}_{0})}{\displaystyle \int }_{{e}_{\min ,n}}^{{e}_{\max ,n}}\displaystyle \frac{{\mathfrak{R}}({\tilde{L}}_{n,j}(e){\tilde{L}}_{n,k}^{* }(e))}{{S}_{n}({\nu }_{n}(e))}\left|n\displaystyle \frac{d\nu }{{de}}\right|\,{de}\\ &&\qquad +\,2{\delta }_{k,{e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{j}| {L}_{n}({e}_{0}){| }^{2}}{{S}_{n}[{\nu }_{n}({e}_{0})]}\\ &&\qquad +\,2{\delta }_{j,{e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{k}| {L}_{n}({e}_{0}){| }^{2}}{{S}_{n}[{\nu }_{n}({e}_{0})]}\\ &&\qquad +\,4\displaystyle \frac{{\delta }_{j,{e}_{0}}{\delta }_{k,{e}_{0}}}{\delta {e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{| {L}_{n}({e}_{0}){| }^{2}}{{S}_{n}({\nu }_{n}({e}_{0}))}.\end{eqnarray} \tag{ 121 }$$

The limits of integration in Equation (121) are defined by Equations (80) and (81). Here, the four terms correspond respectively to the four terms in Equation (115). The first term is directly analogous to that appearing in the S/N (see Equation (103)). Note that in particular, the elements corresponding to the $j={e}_{\mathrm{LSO}}$ and $k={e}_{\mathrm{LSO}}$ terms are nonzero. The ${e}_{\mathrm{LSO}}$ dependence enters in $\nu (e)$ as shown in Equations (18) and (22). However, the first term in Equation (121) is zero for the $j={e}_{0}$ and $k={e}_{0}$ elements. If the binary forms in the detector's sensitive frequency band, the second, the third, and the fourth terms in Equation (115) contribute to this element of the Fisher matrix. The eccentricity integral in the Fisher matrix may be carried analytically over the δ function, which yields the second and third terms in Equation (121). There, ${\delta }_{j,{e}_{0}}$ is the Kronecker δ, which is zero unless j corresponds to the parameter e₀, and similarly for ${\delta }_{k,{e}_{0}}$. Note further that only harmonics with ${f}_{\max ,\det }/\nu ({e}_{0})\geqslant n\geqslant {f}_{\min ,\det }/\nu ({e}_{0})$ contribute to these boundary terms, since otherwise ${S}_{n}({\nu }_{n}({e}_{0}))=\infty$.

Similarly to the ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$, ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{NoPrec}}$ is generally valid for any initial eccentricity in the range $0\lt {e}_{0}\lt 1$, but it is computationally efficient only for low to moderate initial eccentricities (${e}_{0}\lesssim 0.8$). In order to make the calculation computationally efficient for high initial eccentricities, we reverse the order of the sum and the integral in the first term in Equation (121), and truncate the sum over harmonics at ${n}_{\max }(e)$ as in Appendix D.1.1. We get

$$\begin{eqnarray}{{\rm{\Gamma }}}_{{jk}}^{\mathrm{NoPrec},\mathrm{High}} & \approx & 4{\displaystyle \int }_{{e}_{\min ,n}}^{{e}_{\max ,n}}\displaystyle \sum _{n=1}^{{n}_{\max }(e)}\displaystyle \frac{{\mathfrak{R}}({\tilde{L}}_{n,j}(e){\tilde{L}}_{n,k}^{* }(e))}{{S}_{n}({\nu }_{n}(e))}\left|n\displaystyle \frac{d\nu }{{de}}\right|\,{de}\\ & & +2{\delta }_{k,{e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{j}| {L}_{n}({e}_{0}){| }^{2}}{{S}_{n}[{\nu }_{n}({e}_{0})]}\\ & & +2{\delta }_{j,{e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{k}| {L}_{n}({e}_{0}){| }^{2}}{{S}_{n}[{\nu }_{n}({e}_{0})]}\\ & & +4\displaystyle \frac{{\delta }_{j,{e}_{0}}{\delta }_{k,{e}_{0}}}{\delta {e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{| {L}_{n}({e}_{0}){| }^{2}}{{S}_{n}({\nu }_{n}({e}_{0}))}.\end{eqnarray} \tag{ 122 }$$

D.2.2. Eccentric Inspirals with Precession

Following the steps of Appendix D.2.1 for the precession-free model, we may generalize the calculation of the Fisher matrix to include precession similar to Appendix D.1.2. The Fisher matrix, which is computationally efficient for low to moderate initial eccentricities (${e}_{0}\lesssim 0.8$), can be given as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{jk}}^{\mathrm{Prec},\mathrm{Gen}}={{\rm{\Gamma }}}_{{jk}}^{\mathrm{Gen},n}+{{\rm{\Gamma }}}_{{jk}}^{\mathrm{Gen},n+}+{{\rm{\Gamma }}}_{{jk}}^{\mathrm{Gen},n-},\end{eqnarray} \tag{ 123 }$$

where

$$\begin{eqnarray}{{\rm{\Gamma }}}_{{jk}}^{\mathrm{Gen},n} & \approx & 4\displaystyle \sum _{n=1}^{{n}_{\max }({e}_{0})}{\displaystyle \int }_{{e}_{\min ,n}}^{{e}_{\max ,n}}\displaystyle \frac{{\mathfrak{R}}({\tilde{K}}_{n,j}(e){\tilde{K}}_{n,k}^{* }(e))}{{S}_{n}({\nu }_{n}(e))}\left|n\displaystyle \frac{d\nu }{{de}}\right|\,{de}\\ & & +2{\delta }_{j,{e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{k}{| {K}_{n}({e}_{0})| }^{2}}{{S}_{n}({\nu }_{n}({e}_{0}))}\\ & & +2{\delta }_{k,{e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{j}{| {K}_{n}({e}_{0})| }^{2}}{{S}_{n}({\nu }_{n}({e}_{0}))}\\ & & +4\displaystyle \frac{{\delta }_{j,{e}_{0}}{\delta }_{k,{e}_{0}}}{\delta {e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{| {K}_{n}({e}_{0}){| }^{2}}{{S}_{n}({\nu }_{n}({e}_{0}))},\end{eqnarray} \tag{ 124 }$$

and

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{jk}}^{\mathrm{Gen},n\pm }\approx 4\,\displaystyle \sum _{n=1}^{{n}_{\max }({e}_{0})}{\displaystyle \int }_{{e}_{\min ,n}^{\pm }}^{{e}_{\max ,n}^{\pm }}{de}\,\displaystyle \frac{{\mathfrak{R}}({\tilde{K}}_{n,j}^{\pm }(e){\tilde{K}}_{n,k}^{\pm ,* }(e))}{{S}_{n}({\nu }_{n}^{\pm }(e))}\\ &&\quad \times \,\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|\left|\displaystyle \frac{d\nu }{{de}}\right|\\ &&\quad +\,2{\delta }_{j,{e}_{0}}\displaystyle \sum _{n}{\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|}_{{e}_{0}}{\left|\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{k}{| {K}_{n}^{\pm }({e}_{0})| }^{2}}{{S}_{n}({\nu }_{n}^{\pm }({e}_{0}))}\\ &&\quad +\,2{\delta }_{k,{e}_{0}}\displaystyle \sum _{n}{\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|}_{{e}_{0}}{\left|\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{j}{| {K}_{n}^{\pm }({e}_{0})| }^{2}}{{S}_{n}({\nu }_{n}^{\pm }({e}_{0}))}\\ &&\quad +\,4\displaystyle \frac{{\delta }_{j,{e}_{0}}{\delta }_{k,{e}_{0}}}{\delta {e}_{0}}\displaystyle \sum _{n}{\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|}_{{e}_{0}}{\left|\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{| {K}_{n}^{\pm }({e}_{0}){| }^{2}}{{S}_{n}({\nu }_{n}^{\pm }({e}_{0}))},\end{eqnarray} \tag{ 125 }$$

where the integration bounds are given by Equations (85)–(88), and

$$\begin{eqnarray}&&{\tilde{K}}_{n,j}(e)={\partial }_{j}[{K}_{n}({f}_{n})\,{e}^{i{{\rm{\Psi }}}_{n}({f}_{n})}](e),\end{eqnarray} \tag{ 126 }$$

$$\begin{eqnarray}&&{\tilde{K}}_{n,j}^{\pm }(e)={\partial }_{j}[{K}_{n}^{\pm }({f}_{n}^{\pm })\,{e}^{i{{\rm{\Psi }}}_{n}^{\pm }({f}_{n}^{\pm })}](e).\end{eqnarray} \tag{ 127 }$$

Similar to the precession-free case, here

$$\begin{eqnarray}&&{K}_{n}({f}_{n})\equiv {K}_{n}(e({f}_{n}),{f}_{n}),\quad {{\rm{\Psi }}}_{n}({f}_{n})\equiv {{\rm{\Psi }}}_{n}(e({f}_{n}),{f}_{n}),\end{eqnarray} \tag{ 128 }$$

$$\begin{eqnarray}&&{K}_{n}^{\pm }({f}_{n})\equiv {K}_{n}^{\pm }(e({f}_{n}^{\pm }),{f}_{n}^{\pm }),\quad {{\rm{\Psi }}}_{n}^{\pm }({f}_{n}^{\pm })\equiv {{\rm{\Psi }}}_{n}^{\pm }(e({f}_{n}^{\pm }),{f}_{n}^{\pm }),\end{eqnarray} \tag{ 129 }$$

where ${K}_{n}(e,{f}_{n})$ and ${K}_{n}^{\pm }(e,{f}_{n}^{\pm })$ are given by Equations (105)–(107), ${{\rm{\Psi }}}_{n}(e,{f}_{n})$ and ${{\rm{\Psi }}}_{n}^{\pm }(e,{f}_{n}^{\pm })$ are expressed by Equations (35) and (36), and we first differentiate the expressions in the bracket [ ] in Equations (126) and (127) with respect to λ_j, then change variables from f_n and ${f}_{n}^{\pm }$ back to e. Similar to the precession-free case, only harmonics with ${f}_{\max ,\det }/\nu ({e}_{0})\geqslant n\geqslant {f}_{\min ,\det }/\nu ({e}_{0})$ and ${f}_{\max ,\det }/{\nu }^{\pm }({e}_{0})\geqslant n\geqslant {f}_{\min ,\det }/{\nu }^{\pm }({e}_{0})$ contribute to the boundary terms in Equations (124) and (125), since otherwise ${S}_{n}({\nu }_{n}({e}_{0}))=\infty$ and ${S}_{n}({\nu }_{n}^{\pm }({e}_{0}))=\infty$.

For high initial eccentricities, the computationally efficient form of ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{Prec},\mathrm{Gen}}$, ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{Prec},\mathrm{High}}$, can be derived by reversing the order of the sum and the integral in the first term in Equations (124) and (125), and truncating the sum over harmonics at ${n}_{\max }(e)$ as in Appendix D.1.1. Thus, ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{Prec},\mathrm{High}}$ can be given as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{jk}}^{\mathrm{Prec},\mathrm{High}}={{\rm{\Gamma }}}_{{jk}}^{\mathrm{High},n}+{{\rm{\Gamma }}}_{{jk}}^{\mathrm{High},n+}+{{\rm{\Gamma }}}_{{jk}}^{\mathrm{High},n-},\end{eqnarray} \tag{ 130 }$$

where

$$\begin{eqnarray}{{\rm{\Gamma }}}_{{jk}}^{\mathrm{High},n} & \approx & 4{\displaystyle \int }_{{e}_{\min ,n}}^{{e}_{\max ,n}}\displaystyle \sum _{n=1}^{{n}_{\max }(e)}\displaystyle \frac{{\mathfrak{R}}({\tilde{K}}_{n,j}(e){\tilde{K}}_{n,k}^{* }(e))}{{S}_{n}({\nu }_{n}(e))}\left|n\displaystyle \frac{d\nu }{{de}}\right|\,{de}\\ & & +2{\delta }_{j,{e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{k}{| {K}_{n}({e}_{0})| }^{2}}{{S}_{n}({\nu }_{n}({e}_{0}))}\\ & & +2{\delta }_{k,{e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{j}{| {K}_{n}({e}_{0})| }^{2}}{{S}_{n}({\nu }_{n}({e}_{0}))}\\ & & +4\displaystyle \frac{{\delta }_{j,{e}_{0}}{\delta }_{k,{e}_{0}}}{\delta {e}_{0}}\displaystyle \sum _{n}{\left|n\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{| {K}_{n}({e}_{0}){| }^{2}}{{S}_{n}({\nu }_{n}({e}_{0}))},\end{eqnarray} \tag{ 131 }$$

and

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{jk}}^{\mathrm{High},n\pm }\approx 4{\displaystyle \int }_{{e}_{\min ,n}^{\pm }}^{{e}_{\max ,n}^{\pm }}{de}\,\displaystyle \sum _{n=1}^{{n}_{\max }(e)}\displaystyle \frac{{\mathfrak{R}}({\tilde{K}}_{n,j}^{\pm }(e){\tilde{K}}_{n,k}^{\pm ,* }(e))}{{S}_{n}({\nu }_{n}^{\pm }(e))}\\ &&\qquad \times \,\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|\left|\displaystyle \frac{d\nu }{{de}}\right|\\ &&\qquad +\,2{\delta }_{j,{e}_{0}}\displaystyle \sum _{n}{\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|}_{{e}_{0}}{\left|\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{k}{| {K}_{n}^{\pm }({e}_{0})| }^{2}}{{S}_{n}({\nu }_{n}^{\pm }({e}_{0}))}\\ &&\qquad +\,2{\delta }_{k,{e}_{0}}\displaystyle \sum _{n}{\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|}_{{e}_{0}}{\left|\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{{\partial }_{j}{| {K}_{n}^{\pm }({e}_{0})| }^{2}}{{S}_{n}({\nu }_{n}^{\pm }({e}_{0}))}\\ &&\qquad +\,4\displaystyle \frac{{\delta }_{j,{e}_{0}}{\delta }_{k,{e}_{0}}}{\delta {e}_{0}}\displaystyle \sum _{n}{\left|n\pm \displaystyle \frac{1}{\pi }\displaystyle \frac{d\dot{\gamma }}{d\nu }\right|}_{{e}_{0}}{\left|\displaystyle \frac{d\nu }{{de}}\right|}_{{e}_{0}}\displaystyle \frac{| {K}_{n}^{\pm }({e}_{0}){| }^{2}}{{S}_{n}({\nu }_{n}^{\pm }({e}_{0}))}.\end{eqnarray} \tag{ 132 }$$

E.1. Analytic Circular Limit Without Precession

First, we study the circular limit of the eccentric waveform, ${\tilde{h}}_{\times }({\boldsymbol{f}})$ and ${\tilde{h}}_{+}({\boldsymbol{f}})$, defined by Equations (31) and (32) for the NoPrec model. For $e\to 0$, ${B}_{n}^{+}={\delta }_{n,2}$ and ${{\rm{\Psi }}}_{n}^{\pm }={{\rm{\Psi }}}_{n}\mp 2{\gamma }_{c},$ implying that γ_c is degenerate with Φ_c, which we henceforth omit. After integrating the term $\nu ^{\prime} /\dot{\nu }(\nu ^{\prime} )$ over $\nu ^{\prime}$ in Equation (26) according to previous considerations, substituting ν with $f/2$, and expanding the expression of ${h}_{0}{{\rm{\Lambda }}}_{2}$, the polarization components become

$$\begin{eqnarray}&&{h}_{\times }(f)=-2i\sqrt{\displaystyle \frac{5}{96}}\displaystyle \frac{{{ \mathcal M }}_{z}^{5/6}{f}^{-7/6}{e}^{i{\rm{\Psi }}}}{{\pi }^{2/3}{D}_{{\rm{L}}}}\cos {\rm{\Theta }},\end{eqnarray} \tag{ 133 }$$

$$\begin{eqnarray}&&{h}_{+}(f)=-\sqrt{\displaystyle \frac{5}{96}}\displaystyle \frac{{{ \mathcal M }}_{z}^{5/6}{f}^{-7/6}{e}^{i{\rm{\Psi }}}}{{\pi }^{2/3}{D}_{{\rm{L}}}}(1+{\cos }^{2}{\rm{\Theta }}),\end{eqnarray} \tag{ 134 }$$

where the phase function Ψ can be given as

$$\begin{eqnarray}&&{\rm{\Psi }}=2\pi {{ft}}_{c}-{{\rm{\Phi }}}_{c}-\pi /4+\displaystyle \frac{3}{4}{(8\pi {{ \mathcal M }}_{z}f)}^{-5/3}.\end{eqnarray} \tag{ 135 }$$

These are indeed the well-known frequency-domain polarization components of circular binaries in leading order (Cutler & Flanagan 1994). The parameter set characterizing this waveform is

$$\begin{eqnarray}&&{{\boldsymbol{\lambda }}}_{\mathrm{circ}}=\{\mathrm{ln}({D}_{L}),\mathrm{ln}({{ \mathcal M }}_{z}),{\theta }_{N},{\phi }_{N},{\theta }_{L},{\phi }_{L},{t}_{c},{{\rm{\Phi }}}_{c}\}.\end{eqnarray} \tag{ 136 }$$

For validation tests, we calculate the S/N of circular binaries, ${\rm{S}}/{{\rm{N}}}_{\mathrm{circ}}$, for a single aLIGO detector. ${\rm{S}}/{{\rm{N}}}_{\mathrm{circ}}$ is calculated by substituting Equations (133) and (134) into Equation (68) and then using Equation (42). The Fisher matrix of the circular binaries for the parameter set ${{\boldsymbol{\lambda }}}_{\mathrm{circ}}$ for each detector ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{circ}}$ is calculated by substituting Equations (133) and (134) into Equation (68) and then using Equation (44).

E.2. Eccentric Inspiral Without Pericenter Precession

Next, we discuss the validation tests performed for the codes using the NoPrec waveform model. In this case, the parameters are

$$\begin{eqnarray}&&{{\boldsymbol{\lambda }}}_{\mathrm{NoPrec}}=\{\mathrm{ln}({D}_{L}),\mathrm{ln}({{ \mathcal M }}_{z}),{\theta }_{N},\\ &&\qquad {\phi }_{N},{\theta }_{L},{\phi }_{L},{t}_{c},{{\rm{\Phi }}}_{c},{c}_{0},{e}_{0},{\gamma }_{c}\}.\end{eqnarray} \tag{ 137 }$$

Compared to ${{\boldsymbol{\lambda }}}_{\mathrm{circ}}$, ${{\boldsymbol{\lambda }}}_{\mathrm{NoPrec}}$ includes c₀ (set by ${M}_{\mathrm{tot},z}$ and e_LSO; see Equation (22)) and e₀.

E.2.1. Signal-to-noise Ratio

First, we generate a set of source parameters for comparison $({m}_{A},{m}_{B},{D}_{L},{\theta }_{N},{\phi }_{N},{\theta }_{L},{\phi }_{L})$, and compare the output of ${\rm{S}}/{{\rm{N}}}_{\mathrm{circ}}$ with the output of ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$ in the circular limit for a single aLIGO detector (Table 1). Here and in further validation tests, the set of fiducial source parameters generated for comparison are $({m}_{A},{m}_{B},{D}_{L},{\theta }_{N},{\phi }_{N},{\theta }_{L},{\phi }_{L})$, where $\cos {\theta }_{N}$ and $\cos {\theta }_{L}$ are drawn from a uniform distribution between [−1, 1], the set of ${\phi }_{N}$ and ϕ_L are drawn from a uniform distribution between [0, 2π], m_A and m_B are drawn from a uniform distribution between [5 ${M}_{\odot }$, 100 ${M}_{\odot }$], and D_L is drawn from a uniform distribution between [100 Mpc, 1000 Mpc]. The generation of other source parameters are described in detail in the corresponding paragraph. We assume the fiducial value ${t}_{c}={{\rm{\Phi }}}_{c}={\gamma }_{c}=0$ for each binary, and in practice we set ${e}_{0}={10}^{-4}$ and ρ_p0 = 1000 when considering the circular limit (Appendix A). We find that the relative discrepancy between ${\rm{S}}/{{\rm{N}}}_{\mathrm{circ}}$ and ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$ is less than 10⁻³ in all cases.

Next, we examine if the output of ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$ agrees with Figure 11 in O'Leary et al. (2009), which shows the source sky position- and orientation-averaged rms S/N for a single aLIGO detector as a function of ${\rho }_{{\rm{p}}0}$ and M_tot. O'Leary et al. (2009) used the Fourier domain orbit-averaged leading-order waveform (Peters & Mathews 1963), which corresponds to our NoPrec model. We set the sensitivity curve to that used in O'Leary et al. (2009), set e = 0.95, and find the results for several ${\rho }_{{\rm{p}}0}$ and M_tot in the range [5,100] and [10 ${M}_{\odot }$, 1000 ${M}_{\odot }$], respectively. For each ρ_p0 and M_tot, we generate random Monte Carlo samples of source sky location (${\theta }_{N},{\phi }_{N}$) and binary orientation (${\theta }_{L},{\phi }_{L}$) as introduced above in this section. We find that the rms of the ${({\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}})}^{2}$ distributions are in agreement with the results of Figure 11 in O'Leary et al. (2009).

Finally, we generate a set of source parameters for comparison, and compare the output of ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$ and ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{High}}$ for several high e₀ and ${\rho }_{{\rm{p}}0}$, and for a single aLIGO detector. In particular, e₀ and ${\rho }_{{\rm{p}}0}$ are drawn from a uniform distribution between ]0.9, 1[ and [5, 1000] in these calculations, respectively. The relative discrepancy between ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{Gen}}$ and ${\rm{S}}/{{\rm{N}}}_{\mathrm{NoPrec}}^{\mathrm{High}}$ is less than 10⁻³ in all cases.

E.2.2. Fisher Matrix

Since the parameter set of the leading-order circular and eccentric binaries, ${{\boldsymbol{\lambda }}}_{\mathrm{circ}}$ and ${{\boldsymbol{\lambda }}}_{\mathrm{NoPrec}}$ differ—see Equations (136) and (137)—we cannot simply compare the output of ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{circ}}$ with the output of ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{NoPrec},\mathrm{Gen}}$ in the circular limit. Therefore, we first restrict ${{\boldsymbol{\lambda }}}_{\mathrm{NoPrec}}$ to the parameter set ${{\boldsymbol{\lambda }}}_{\mathrm{circ}}$. Next, we generate a set of source parameters for comparison, and compare the output of ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{circ}}$ with the output of ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{NoPrec},\mathrm{Gen}}$ in the circular limit for the Fisher matrix elements corresponding to the circular parameters ${{\boldsymbol{\lambda }}}_{\mathrm{circ}}$ and for a single aLIGO detector. The relative discrepancy between ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{circ}}$ and ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{NoPrec},\mathrm{Gen}}$ is less than 10⁻² in all cases.

Finally, we generate a set of source parameters for comparison, and compare the output of ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{NoPrec},\mathrm{Gen}}$ and ${{\rm{\Gamma }}}_{{jk}}^{\mathrm{NoPrec},\mathrm{High}}$ for the Fisher matrix elements corresponding to the parameter set ${{\boldsymbol{\lambda }}}_{\mathrm{NoPrec}}$ in the high-eccentricity limit for several ${\rho }_{{\rm{p}}0}$ and for a single aLIGO detector. Similarly to Appendix E.2.1, e₀ and ${\rho }_{{\rm{p}}0}$ are drawn from a uniform distribution between ]0.9, 1[ and [5, 1000] in these calculations, respectively. The relative discrepancy between Fisher matrix elements is less than 10⁻³ in all cases.

E.3. Eccentric Pericenter-precessing Binary Waveforms