Convergence of Fourier-domain templates for inspiraling eccentric compact binaries

TL;DR

This paper analyzes the convergence and accuracy of frequency-domain inspiral templates for eccentric black hole binaries, crucial for gravitational wave detection with LISA, focusing on systematic and statistical errors up to 2PN order and eccentricity $e_0^6$.

Contribution

It evaluates the convergence of eccentric inspiral templates and assesses the impact of neglecting higher-order terms on systematic errors in gravitational wave analysis.

Findings

Eccentric waveforms increase statistical errors due to parameter correlations.

Including terms of order $e_0^2$ or higher keeps systematic errors negligible.

Templates converge well up to 2PN order and $e_0^6$ in eccentricity.

Abstract

The space-based detector LISA may observe gravitational waves from the early inspiral of stellar-mass black hole binaries, some of which could have significant eccentricity. Current gravitational waveform templates are only valid for small orbital velocities (i.e., in a post-Newtonian expansion) and small initial eccentricity $e_0$ ("post-circular" expansion). We conventionally define $e_0$ as the eccentricity corresponding to an orbital frequency of $5 \text{ mHz}$, and we study the convergence properties of frequency-domain inspiral templates that are accurate up to 2PN and order $e_0^6$ in eccentricity. We compute the so-called "unfaithfulness" between the full template and "reduced" templates obtained by dropping some terms in the phasing series; we investigate the conditions under which systematic errors are negligible with respect to statistical errors, and we study the convergence properties of statistical errors. In general, eccentric waveforms lead to larger statistical errors than circular waveforms due to correlations between the parameters, but the error estimates do not change significantly as long as we include terms of order $e_0^2$ or higher in the phasing.

Figures

Tables

Table 1. Table 1: Template naming conventions. According to the alphabetical naming convention, the template obtained by dropping terms corresponding to the letters B, C and D from the phase is called “template B”, and so on (see text). In the curly bracket convention, the { e 0 2 } superscript subscript 𝑒 0 2 \left\{e_{0}^{2}\right\} template is obtained by retaining terms of order up to e 0 2 superscript subscript 𝑒 0 2 e_{0}^{2} (i.e., the two leftmost columns). Additionally, “template H” corresponds to dropping cells C, D and E from the fiducial template.

2PN ($x^2$)	A	B	C	D
1.5PN ($x^{3/2}$)				E
1PN ($x^1$)				F
Newtonian ($x^0$)				G
	$e_0^0$	$e_0^2$	$e_0^4$	$e_0^6$

Table 2. Table 2: Maximum detection SNR – computed using Eq. ( 22 ) and the median faithfulness – such that systematic errors are smaller than statistical errors when we drop the term corresponding to each box. Numbers correspond to e 0 = 0.01 , 0.04 , 0.07 subscript 𝑒 0 0.01 0.04 0.07 e_{0}=0.01,0.04,0.07 and 0.1 0.1 0.1 . ∞ \infty means that the median unfaithfulness is zero within machine roundoff errors.

2PN $x^2$	$3.7,3.7,3.7,3.7$	$(100,6.0,1.7,0.76)\times {10}^2$	$(2300,8.3,0.80,0.17)\times {10}^4$	$(2.1,0.26,0.009,0.001)\times {10}^8$
1.5PN ($x^{3/2}$)	$2.6,2.6,2.6,2.6$	$170,16,7.0,5.0$	$(6000,23,2.2,0.58)\times {10}^2$	$(\infty ,130,4.3,0.50)\times {10}^4$
1PN ($x^1$)	$2.6,2.6,2.6,2.6$	$28,4.0,3.1,2.9$	$(500,1.8,0.21,0.081)\times {10}^2$	$(210000,60,2.0,0.22)\times {10}^3$
Newtonian ($x^0$)	$2.7,2.7,2.7,2.6$	$2.7,2.6,2.6,2.6$	$(270,3.7,2.7,2.6)$	$(17000,4.0,0.18,0.050)\times {10}^2$
	$e_0^0$	$e_0^2$	$e_0^4$	$e_0^6$

Table 3. Table 3: Maximum e 0 subscript 𝑒 0 e_{0} such that detections with ρ < 25 𝜌 25 \rho<25 have systematic errors smaller than statistical errors, when terms corresponding to the respective cells in the table are dropped. A dash means that the mismatch is so low that systematic errors are larger than statistical errors for all e 0 subscript 𝑒 0 e_{0} in the range we consider.

2PN ($x^2$)	$-$	$0.16$	$0.21$	$0.30$
1.5PN ($x^{3/2}$)	$-$	$0.03$	$0.12$	$0.24$
1PN ($x^1$)	$-$	$0.008$	$0.06$	$0.18$
Newtonian ($x^0$)	$-$	$0.0007$	$0.015$	$0.05$
	$e_0^0$	$e_0^2$	$e_0^4$	$e_0^6$

Equations

$$x\equiv {\left(\frac{2\pi G{m}_z{f}_{\mathrm{orb}}}{c^3}\right)}^{2/3},$$

$$\tilde{h}(f)=\frac{\sqrt{3}}{2}\tilde{\mathcal{A}}{\left(\frac{G{m}_z\pi f}{c^3}\right)}^{-7/6}\sum _{j=1}^6\xi (f,j){\left(\frac{j}{2}\right)}^{2/3}{e}^{-i(\pi /4+\Psi (f,j)-{\varphi }_D(f,j))},$$

$$\tilde{\mathcal{A}}$$

$$\xi (f,j)$$

$$\Gamma (f,j)=F_{+}(f,j)C_{+}^j(f)+F_{\times }(f,j)C_{\times }^j(f),$$

$$\Sigma (f,j)=F_{+}(f,j)S_{+}^j(f)+F_{\times }(f,j)S_{\times }^j(f).$$

$$\Psi (f,j)$$

$$e(f,j)$$

$$-\frac{11042329}{1108992}\chi (f,j{)}^{-19/6}+\frac{50259743}{6653952}\chi (f,j{)}^{-95/18})e_0^5,$$

$$\mathcal{C}(f,j)=1-\frac{2355}{1462}{e}_0^2\chi (f,j{)}^{-19/9}+\left(-\frac{2608555}{444448}\chi (f,j{)}^{-19/9}+\frac{5222765}{998944}\chi (f,j{)}^{-38/9}\right)e_0^4$$

$$+\left(-\frac{1326481225}{10134144}\chi (f,j{)}^{-19/9}+\frac{173355248095}{455518464}\chi (f,j{)}^{-38/9}-\frac{75356125}{3326976}\chi (f,j{)}^{-19/3}\right)e_0^6.$$

$${\varphi }_D(f,j)=2\pi Rf\sin {\bar{\theta }}_S\cos \left(\bar{\varphi }(f,j)-{\bar{\varphi }}_S\right),$$

$$\bar{\varphi }(f,j)=\frac{2\pi t(f,j)}{{\mathrm{T}}_0}=\frac{1}{{\mathrm{T}}_0}[2\pi t_c+\frac{G{m}_z}{\eta c^3}{\left(\frac{2G{m}_z\pi f}{j{c}^3}\right)}^{-8/3}(-\frac{5}{256}+\frac{785e_0^2}{11008}\chi (f,j{)}^{-19/9}+e_0^4(-\frac{5222765}{14475264}\chi (f,j{)}^{-38/9}$$

$$+\frac{2608555}{10039296}\chi (f,j{)}^{-19/9})+e_0^6\left(\frac{75356125}{35487744}\chi (f,j{)}^{-19/3}+\frac{17355248095}{6600720384}\chi (f,j{)}^{-38/9}+\frac{1326481225}{2288959488}\chi (f,j{)}^{19/9}\right))].$$

$$T_{\text{rr}}=\frac{5G\mathcal{M}}{96c^3}{\left(\frac{G\mathcal{M}\omega }{c^3}\right)}^{-8/3}\left[\frac{{\left(1-e^2\right)}^{7/2}}{1+\frac{73}{24}{e}^2+\frac{37}{96}{e}^4}\right].$$

$$M$$

$$(h_1,h_2)$$

$${\rho }^2(h)$$

$$S_h(f)$$

$$+S_c(f),$$

$$P_{\mathrm{OMS}}$$

$$P_{\mathrm{acc}}$$

$$\left[1+{\left(\frac{f}{8\mathrm{mHz}}\right)}^4\right]\mathrm{H}{\mathrm{z}}^{-1},$$

$$S_c(f)$$

$$\left[1+\tanh (\gamma (f_k-f))\right]{\mathrm{Hz}}^{-1},$$

$${\tau }_{AB}=4\Re {\int }_0^{\infty }\left[{\frac{\partial \tilde{h}(f)}{\partial p_A}}^{\ast }\frac{\partial \tilde{h}(f)}{\partial p_B}\right]\frac{1}{S_h(f)}df.$$

$$4\Re {\int }_0^{\infty }\left[{\frac{\partial {\tilde{h}}_1(f)}{\partial p_A}}^{\ast }\frac{\partial {\tilde{h}}_1(f)}{\partial p_B}\right]\frac{1}{S_h(f)}df,$$

$$4\Re {\int }_0^{\infty }\left[{\frac{\partial {\tilde{h}}_2(f)}{\partial p_A}}^{\ast }\frac{\partial {\tilde{h}}_3(f)}{\partial p_B}\right]\frac{1}{S_h(f)}df,$$

$$1-E(M_{\mathrm{sta}})=\frac{(D-1)}{2{\mathrm{SNR}}^2},$$

$$1-M\ll \frac{(D-1)}{2{\mathrm{SNR}}^2}.$$

$$\tilde{h}(f)$$

$${\tilde{h}}_j(f)$$

$${\Psi }_j(f)$$

$$\frac{\partial {\tilde{h}}_j(f)}{\partial p_A}$$

$$\Re {\int }_{f_{\min }}^{f_{\max }}\frac{1}{S_n(f)}\frac{\partial {\tilde{h}}_{j_1}(f)}{\partial p_A}\frac{\partial {\tilde{h}}_{j_2}(f{)}^{\ast }}{\partial p_B}df$$

$${\Psi }_{j_1}-{\Psi }_{j_2}$$

Full text

I Introduction

The first two observing runs by Advanced LIGO and Virgo LIGO Scientific Collaboration and Virgo Collaboration (2018) detected gravitational waves (GWs) from 11 compact binary coalescence events, and the era of GW astronomy has begun. In the near future, GW astronomy has the potential to answer several important open questions: it will test general relativity in the strong gravity regime, probe the neutron star equation of state, shed light on astrophysical formation scenarios of compact object binaries, and potentially resolve outstanding open problems in cosmology Sathyaprakash and Schutz (2009); Yunes and Siemens (2013); Will (2014); Berti et al. (2015); Barack et al. (2018).

This paper is motivated by the possibility of identifying the formation channels of binary black holes (BBHs) using a combination of Earth- and space-based GW detectors Sesana (2016); Wong et al. (2018); Cutler et al. (2019); Gerosa et al. (2019), and in particular by the prospect of using eccentricity measurements to distinguish between two of the main proposed formation channels Nishizawa et al. (2016); Breivik et al. (2016); Nishizawa et al. (2017): field and dynamical formation. Most BBHs are expected to circularize by the time they enter the most sensitive band of ground based detectors, so we must rely on measurements of masses, spins, redshifts and kicks to distinguish between different formation scenarios Gerosa et al. (2013); Stevenson et al. (2015); Kovetz et al. (2017); Fishbach and Holz (2017); Zevin et al. (2017); Barrett et al. (2018); Gerosa and Berti (2017); Talbot and Thrane (2018); Wysocki et al. (2018); Gerosa et al. (2018). However, typical BBH eccentricities are larger for binaries formed dynamically than for binaries formed in the field (see e.g. Fig. 1 of Nishizawa et al. (2017)), at the typical frequencies $\sim 10^{-2}$ Hz targeted by the planned Laser Interferometer Space Antenna (LISA) Amaro-Seoane et al. (2013). Therefore LISA has the potential to measure BBH eccentricities and to shed light on their formation channel in a way that is complementary to ground based detectors Nishizawa et al. (2016); Breivik et al. (2016); Nishizawa et al. (2017); Samsing (2018); Samsing et al. (2018); Samsing and D’Orazio (2018); D’Orazio and Samsing (2018); Samsing and D’Orazio (2019); Zevin et al. (2019); Rodriguez et al. (2018).

There is a large body of work extending the pioneering study of GWs from eccentric compact binaries by Peters and Mathews Peters and Mathews (1963), where the binary dynamics was treated at Newtonian order, to higher post-Newtonian (PN) orders. The first analytic Fourier-domain templates were calculated within the stationary-phase approximation for arbitrary initial eccentricity, including the effect of periastron advance, in Mikoczi et al. (2012). Analytical expressions for the decay of the orbital parameters under radiation reaction using hypergeometric functions were derived in Pierro and Pinto (1996); Pierro et al. (2002); Mikoczi et al. (2015) at Newtonian order, and in Mikoczi et al. (2015) at 1PN order. A Kepler-like parametrization was introduced and extended up to 3PN order in Damour and Deruelle (1985); Damour and Schafer (1988); Memmesheimer et al. (2004), and the decay of the orbital parameters under radiation reaction was computed in Schäfer (1985); Iyer and Will (1995); Königsdörffer et al. (2003); Nissanke and Blanchet (2005). Now the evolution of the orbital parameters under radiation reaction is known up to 3.5 PN order Damour et al. (2004); Königsdörffer and Gopakumar (2006). A time-domain template for compact binaries where the orbital elements were evolved via numerical integration was introduced in Ref. Tanay et al. (2016), and it can be regarded as the eccentric extension of the popular TaylorT4 approximant for quasicircular binaries.

In the analysis of quasicircular binary inspirals it is common to use analytic Fourier-domain templates computed within the stationary-phase approximation (SPA), such as the TaylorF2 template. Reference Yunes et al. (2009) generalized these Fourier-domain templates to eccentric binaries, computing templates which are valid at Newtonian order and up to order $e_{0}^{8}$ in a small-eccentricity expansion of the phase (here $e_{0}$ is defined to be the eccentricity at some reference orbital frequency $f_{\rm orb}=f_{\rm orb,\,0}$).

This work was extended to 2PN-$e_{0}^{6}$ in Ref. Tanay et al. (2016), which will be the starting point of our study. Moore et al. Moore et al. (2016) extended Ref. Tanay et al. (2016) to 3PN order, but only at leading order in $e_{0}$. All of the templates above are valid at Newtonian order in amplitude and do not include the effect of periastron advance in the phase. Reference Moore et al. (2018) constructed analytic templates which are valid for arbitrary initial eccentricity $e_{0}$ within the SPA using a truncated sum of harmonics and hypergeometric functions. This work was recently extended to 3PN accuracy Moore and Yunes (2019), and efforts are underway to extend Ref. Yunes et al. (2009) up to 3PN-$e_{0}^{6}$ order in phase and 1PN order in amplitude, incorporating periastron advance effects Tiwari et al. (2019).

With so many parallel efforts on analytic Fourier-domain eccentric templates underway, it is crucial to investigate the convergence of these proposed GW templates. This is the main goal of our work. We study the convergence properties of the 2PN-$e_{0}^{6}$ accurate template proposed in Ref. Tanay et al. (2016). This “fiducial template” is a sum over harmonics (labeled by $j$), where each harmonic has a phase which itself is a bivariate series in the initial eccentricity $e_{0}$ and in the PN parameter

[TABLE]

where $f_{\rm orb}$ is the orbital frequency for a circular binary, and $m_{z}=(1+z)m$ is the redshifted total mass of the binary. A preliminary, more limited investigation of the convergence of this bivariate series in the context of parameter estimation can be found in Nishizawa et al. (2016). Our work should be helpful in guiding future efforts to extend the above templates to higher orders, and it can readily be generalized as soon as more accurate templates become available.

We focus on the convergence of the expansion of the phasing (rather than the amplitude) because the phasing is known to have greater impact on detectability and parameter estimation. We first drop some terms from the fiducial template to get (presumably) less accurate “reduced” templates, then we perform calculations of the so-called “unfaithfulness” between these reduced templates and the fiducial 2PN-$e_{0}^{6}$ accurate template to assess the importance of the dropped term(s). We also investigate the conditions under which systematic errors due to dropping high-order terms from the fiducial template exceed the statistical errors. This is useful because whenever systematic errors are negligible with respect to statistical errors, one can choose a truncated template to improve computational efficiency in parameter estimation. Finally we study the convergence of statistical errors.

The paper is organized as follows. Section II describes the waveform model and gives details on our calculation of matches and Fisher matrices. Section III is an overview of data analysis concepts that are relevant for our study. Our results on unfaithfulness are presented in Sec. IV.1, the comparison of systematic and statistical errors is shown in Sec. IV.2, and the convergence of statistical errors is studied in Sec. IV.3. In Sec. V we summarize our main results and outline directions for future work. To improve readability, some technical details are relegated to the Appendices. Appendix A illustrates why certain cross terms in the Fisher matrix integrands can be neglected due to their oscillatory nature, and Appendix B defines beam pattern functions and other quantities appearing in the calculation of the GW strain. The code used in our analysis is publicly available online Ecc .

II Waveforms, match and Fisher matrix calculations

In this section we describe our waveform model, and then we give details of our unfaithfulness and Fisher matrix calculations.

Our analytic frequency-domain GW template for compact binaries inspiraling in eccentric orbits uses Newtonian amplitudes, the first six harmonics $(j=1,\dots,6)$, and a 2PN-$e_{0}^{6}$ accurate phase. Here $e_{0}$ is the eccentricity at which the orbital frequency of the binary is $f_{0}/2$, and we (somewhat arbitrarily) set $f_{0}=10$ mHz. In other words, when the binary has eccentricity $e_{0}$, the second harmonic (which dominates the signal for small eccentricities) has frequency $f_{0}$. From Appendix B and Eq. (3.11) of Yunes et al. (2009) it follows that if we include six harmonics, the expression of $e$ used in the amplitude should be accurate up to ${\cal}{O}(e_{0}^{3})$. Below we list the relevant expressions only at Newtonian order for illustration, but in the actual calculations we retained all terms up to 2PN-$e_{0}^{6}$ order; these can be found in Appendix A of Tanay et al. (2016). The waveform in the stationary phase approximation has the form

[TABLE]

where the symmetric mass ratio $\eta=m_{1}m_{2}/m^{2}$,

[TABLE]

the quantity $D_{L}$ is the luminosity distance to the source, and

[TABLE]

The definitions of $F_{+},F_{\times},C_{+}^{j},C_{\times}^{j},S_{+}^{j},S_{\times}^{j}$ are given in Appendix B. The Fourier phase $\Psi(f,j)$ and $e(f,j)$, up to the leading PN order and sixth order in $e_{0}$, are given by

[TABLE]

with $\chi(f,j)=2f/(jf_{0})$ and

[TABLE]

We also introduced the Doppler phase

[TABLE]

where $R=1AU$ and the orbital phase $\bar{\phi}(f,j)$ of LISA’s barycenter around the Sun is

[TABLE]

Here $T_{0}=1$ yr, the angles $(\bar{\theta}_{S},\bar{\phi}_{S})$ define the direction of the source in the Solar barycenter frame, and $t_{c}$ and $\phi_{c}$ denote the time and phase at coalescence, respectively Berti et al. (2005); Yagi and Tanaka (2010). The amplitudes $C_{+}^{j},C_{\times}^{j},S_{+}^{j},S_{\times}^{j}$ are computed by keeping the first two terms in the expansion of Eq. (7): cf. Appendix B. This is consistent with retaining a Newtonian amplitude and six harmonics. The beam pattern functions $F_{+}$ and $F_{\times}$ depend on $f$ and $j$ through $\bar{\phi}(f,j)$. The equations listed above are of order $x^{0}$ in a PN expansion, hence they do not depend on $x$.

The above template is slightly modified with respect to Ref. Tanay et al. (2016): the amplitude has an extra factor of $\sqrt{3}/2$ to account for the $60^{\circ}$ opening angle of the LISA arms Berti et al. (2005), and the Fourier phase has been changed from $\Psi$ to $(\Psi-\phi_{D})$ to account for the Doppler phase due to the motion of the detector around the Sun Nishizawa et al. (2016); Yagi and Tanaka (2010).

II.1 Fiducial template and truncated templates

We will now introduce the structure of the templates used in our calculations. We refer to the 2PN-$e_{0}^{6}$ order accurate template of Tanay et al. (2016) as the “fiducial template”. To assess convergence, we also consider various “truncated templates,” i.e. templates derived from the fiducial one by dropping certain terms in the phase. The amplitude of all templates is accurate at Newtonian order and ${\cal{O}}(e^{4})$, because the phase plays a more important role than the amplitude for detection and parameter estimation (see e.g. Damour et al. (2004); Königsdörffer and Gopakumar (2006); Yunes et al. (2009); Tanay et al. (2016)).

Table 1 illustrates the difference between the various templates. Each cell in the table represents a term of a certain order in the PN frequency parameter $x$ (rows) and in the initial eccentricity $e_{0}$ (columns). We will use two different notations to distinguish between templates.

A “letter-based” template means that we drop all terms of order greater than or equal to the corresponding cell in the table. For example, “waveform A” is obtained by neglecting cells A, B, C and D, i.e. all of the 2PN corrections to the waveform; “waveform B” is obtained by neglecting cells B, C and D, i.e. all 2PN corrections of order $e_{0}^{2}$ and higher in the initial eccentricity; and “waveform C” is obtained by neglecting cells C and D. Similarly, “waveform G” is obtained by neglecting cells G, F, E and D; “waveform F” is obtained by neglecting cells F, E and D; and “waveform E” is obtained by neglecting cells E and D. “Waveform D” corresponds to neglecting only the 2PN, $e_{0}^{6}$ term.

We will also use a curly bracket notation $\left\{y^{n}\right\}$, meaning that the phase is $y^{n}$ accurate in the parameter $y$, where $y$ stands either for $e_{0}$ or for $x$ in the bivariate series for the Fourier phase of Eq. (6). For example, the $\left\{e_{0}^{2}\right\}$ template is accurate up to order $e_{0}^{2}$ (and 2PN) in phase, i.e. we retain the first two columns from the left in Table 1. Likewise, for the $\left\{x^{1}\right\}$ template we retain terms up to 1PN (and order $e_{0}^{6}$) in phase, i.e. the two bottom rows in Table 1.

The fiducial template of Eq. (2) is a series of harmonics labeled by the integer $j$, where the phase $\Psi(f,j)$ of each harmonic is itself a bivariate series in $x$ and $e_{0}$: cf. Eqs. (6) and (8). Here we focus on the convergence of $\Psi(f,j)$ as a bivariate series because, as already mentioned, the phase of a GW template is more important than the amplitude (as long as $e_{0}$ is small, so that a small-$e_{0}$ expansion is valid). 111At Newtonian order, the radiation reaction timescale $T_{\text{rr}}=\omega/\dot{\omega}$ (where $\omega=2\pi f_{\rm orb}$ is the angular frequency) is

$\displaystyle T_{\text{rr}}=\frac{5G{\cal{M}}}{96c^{3}}\left({\frac{G\cal{M}\omega}{c^{3}}}\right)^{-8/3}\left[\frac{\left(1-e^{2}\right)^{7/2}}{1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}}\right].$

The quantity in square brackets – say, $Z(e)$ – must be expanded for small $e$ to finally arrive at the expression of the Fourier phase $\Psi$ which occurs in Eq. (2) Yunes et al. (2009). Any Taylor series has a radius of convergence equal at most to the distance from the expansion point (here $e=0$) and the nearest singularity in the complex plane Arfken and Weber (2005), which here is located at $e\sim\pm~{}0.58~{}i$. Therefore none of our templates should be trusted beyond $e_{0}\sim 0.58$ (although they may become unfaithful for much smaller values of $e_{0}$). In this work we have retained only the first six harmonics $(j=1,\dots,6)$ in all of our templates. The convergence of the harmonic expansion is an interesting topic for future work.

II.2 Match calculations

Our convergence analysis of PN, small-eccentricity waveforms is based on some data analysis concepts that we introduce below. First of all, we define the “faithfulness” $M$ between two GW signals $h_{1}(t)$ and $h_{2}(t)$ as the following integral, maximized over the time and phase of coalescence $t_{c}$ and $\phi_{c}$:

[TABLE]

and the “unfaithfulness” as $(1-M)$. The inner product between two waveforms $(h_{1},h_{2})$ is defined as

[TABLE]

where $\tilde{h}_{i}(f)$ stands for the Fourier transform of $h_{i}(t)$ and $S_{h}(f)$ is the noise power spectral density of the LISA detector Robson et al. (2018). The signal-to-noise ratio (SNR) $\rho$ of template $h$ in the detector can be estimated by

[TABLE]

Unlike Robson et al. (2018), we do not use a sky-averaged response, and thus we should not include the corresponding factor of 5 in the noise curve. We also treat the two channels separately (so our noise curve differs by an extra factor of 2) and we include the geometrical factor $\sqrt{3}/2$ in the waveform definition (2), yielding an extra factor of 3/4. For these reasons, we use the noise curve of Robson et al. (2018) without the overall factor 10/3, i.e.

[TABLE]

where $L=2.5~{}\mathrm{Gm}$ and $f_{*}=19.09~{}\mathrm{mHz}$. Furthermore $P_{\mathrm{OMS}}(f)$, $P_{\mathrm{acc}}$ and the confusion noise $S_{c}(f)$ are given by

[TABLE]

where $A=9\times 10^{-45}$ and all parameters have been chosen corresponding to an observation time of 2 years: $\alpha=0.165~{}\mathrm{Hz}^{-1},~{}\beta=299~{}\mathrm{Hz}^{-1},~{}\kappa=611~{}\mathrm{Hz}^{-1},~{}\gamma=1340~{}\mathrm{Hz}^{-1}$ and $f_{k}=0.00173~{}\mathrm{Hz}$.

II.3 Fisher matrix and cross terms

Our eccentric GW templates depend on 11 parameters: $\ln{\cal{M}}_{z},\ln\eta,t_{c},\phi_{c},\ln{D_{L}},e_{0},\bar{\theta}_{S},\bar{\phi}_{S},\bar{\theta}_{L},\bar{\phi}_{L}$ and $\beta$. Here ${\cal{M}}_{z}=\eta^{3/5}m_{z}$ is the redshifted chirp mass, $\eta=m_{1}m_{2}/m^{2}$ is the symmetric mass ratio, $t_{c}$ and $\phi_{c}$ are the coalescence time and orbital phase, $D_{L}$ is the luminosity distance and $e_{0}$ is the initial eccentricity, defined as the eccentricity corresponding to an orbital frequency $f_{\rm orb}=5$ mHz (so that the second harmonic of the radiation is at $10$ mHz). The angles $(\bar{\theta}_{S},\bar{\phi}_{S})$ and the angles $(\bar{\theta}_{L},\bar{\phi}_{L})$ define the direction of the source and of the binary’s orbital angular momentum in the Solar barycenter frame. Finally, $\beta$ represents the position of the pericenter in the binary’s orbital frame Yunes et al. (2009); Martel and Poisson (1999).

Let ${\bf p}=\{p_{A}\}$ be a vector whose components are any of these eleven parameters. The Fisher matrix is defined as the matrix with elements

[TABLE]

In the high SNR regime, the statistical errors associated with estimating the parameters can be approximated by the square root of the diagonal elements of the inverse Fisher matrix Vallisneri (2008).

Each of our templates $\tilde{h}(f)$ is obtained by summing over the first six harmonics, i.e. $\tilde{h}(f)=\sum_{j=1}^{6}\tilde{h}_{j}$. Therefore the integrand above contains 36 terms: 6 “diagonal terms” where both harmonics are the same, and 30 “cross terms” involving different harmonics. Examples of a diagonal term and of a cross term are

[TABLE]

and

[TABLE]

respectively.

Numerical calculations show that the cross terms oscillate rapidly and that they do not significantly contribute to the integral, so they can be dropped. An analytical justification for this approximation can be found in Appendix A.

II.4 Binary catalog and cosmology

To investigate the statistical properties of our waveforms we use a catalog of $1000$ systems. The individual source-frame masses of the binary components are uniformly distributed between 5$$\textup{M}_{\odot} and 45$$\textup{M}_{\odot}. The angles $\bar{\theta}_{S},\bar{\phi}_{S},\bar{\theta}_{L},\bar{\phi}_{L}$ are uniformly distributed over the sphere, and the angle $\beta$ is uniformly distributed between $[0,2\pi]$. We set $t_{c}$ and $\phi_{c}$ equal to [math]. We truncate the Fisher matrix integrals of Eq. (18) at $f_{\rm max}=1$ Hz, and we choose $f_{\rm min}$ so that the observation time is 2 years. We assume a $\Lambda$CDM cosmology and a spatially flat universe with $H_{0}=67.36$ km/s/Mpc, $\Omega_{M}=0.3153$ and $\Omega_{\Lambda}=0.6847$ and we fix $z=0.1$ (corresponding to $D_{L}=447.8$ Mpc) for all binaries.

III Some data analysis background and motivation

There are two kinds of parameter estimation errors: systematic and statistical. Systematic errors are due (e.g.) to mismodeling of GW signals, and statistical errors are due to the noise in “ideal” detectors (real detectors usually contribute to systematic errors as well) Creighton and Anderson (2011). If the systematic errors associated with neglecting some terms in the waveform template are smaller than statistical errors, we can safely neglect those terms and increase the computational efficiency of parameter recovery without compromising its accuracy. Here we follow Appendix G of Chatziioannou et al. (2017) and we introduce a criterion to decide whether systematic errors are smaller than statistical errors.

Consider a detection scenario where systematic errors are negligible, so all parameter estimation errors are statistical and due to noise. Assume also that the SNR for this detection is large enough that the posterior probability distribution is sharply peaked close to the true parameter values $\bm{p}_{0}$ Creighton and Anderson (2011). For a parameter vector $\bm{p}$ close to $\bm{p}_{0}$, the unfaithfulness is $1-M_{\mathrm{sta}}$, and the mean value of the unfaithfulness over the posterior probability is Chatziioannou et al. (2017)

[TABLE]

where $E$ denotes the expectation value, $D=11$ is the dimension of our parameter space, and the subscript “sta” stands for “statistical.”

Let us now include systematic errors. If we demand systematic errors to be negligible with respect to statistical errors, then the unfaithfulness $(1-M)$ due to GW mismodeling should be negligible with respect to the expected value of the unfaithfulness in Eq. (21), i.e.

[TABLE]

We also need a measure of the overall statistical uncertainty. Following e.g. Lyons Lyons (1989), we can associate an $n$-dimensional error ellipsoid (about the maximum of the posterior distribution) with a Gaussian posterior in an $n$-dimensional parameter space. This ellipsoid is the region of $1~{}\sigma$ confidence interval, within which the parameter vector can be found with probability $\sim 0.68$. The square roots of the diagonal elements of the inverse Fisher matrix are the projections of this error ellipse on the parameter axes. Define $\epsilon$ as the product of Fisher errors on all parameters, and $\epsilon_{0}$ as the volume of the error ellipsoid. It can be shown that $\epsilon_{0}$ is given by the square root of the determinant of the inverse Fisher matrix and that $\epsilon\geq\epsilon_{0}$, where the equality is realized when the parameters are uncorrelated with each other, while having $\epsilon/\epsilon_{0}>1$ means that there is some correlation between the parameters. The parameter $\epsilon_{0}$ can be considered an overall measure of statistical errors and it is unaffected by a linear transformation of the parameters (or equivalently, by a rotation of the parameter axes): in fact, $\epsilon_{0}$ measures the volume of the error ellipsoid, which should be independent of the orientation of the parameter axes. The ratio $\epsilon/\epsilon_{0}$ (“correlation factor”) quantifies the degree of correlation among the parameters, or the misalignment of the error ellipsoid with the parameter axes.

IV Results

In this section we present our results on the convergence properties of the 2PN-${e_{0}}^{6}$ accurate bivariate template (the fiducial template) proposed in Ref. Tanay et al. (2016), as measured in terms of the unfaithfulness and Fisher matrix errors. In subsection IV.1 we compute the unfaithfulness due to truncating the fiducial template in various ways (described below), to get an idea of the relative importance of various terms. Subsection IV.2 identifies criteria under which systematic errors due to neglecting certain terms in the phase are smaller than statistical errors. Finally, subsection IV.3 discusses the convergence properties of the statistical errors.

IV.1 Unfaithfulness of truncated templates

We compute unfaithfulness distributions for the 1000 binaries in our catalog. We compare the eight truncated templates A-H defined in Sec. II.1 (cf. Table 1) against the fiducial template for three selected initial eccentricities $(e_{0}=0.04,\,0.07,\,0.1)$. A large value of the unfaithfulness indicates that the dropped term(s) are significant.

Our results are shown in Fig. 1. The histograms in the top panel show the effect of dropping terms of various orders in the initial eccentricity $e_{0}$ at 2PN order (i.e., we move “horizontally” along the top row Table 1). As expected, the unfaithfulness decreases as we move from template A to template D: template D is the closest to the fiducial template. Similarly, the middle panel shows the effect of dropping terms of various PN orders at order $e_{0}^{6}$ in the initial eccentricity (i.e., we move vertically along the right column of Table 1). Again, as we move from template G to template D (thereby dropping all terms higher than the Newtonian, 1PN, 1.5PN and 2PN order terms, respectively) the unfaithfulness decreases. Finally, the bottom panel corresponds to moving diagonally inwards along Table 1, starting from the top-right corner. In each panel, the unfaithfulness gets larger as we increase $e_{0}$: this is expected, since all of these waveforms are small-eccentricity expansions.

An interesting exception is template A. In this case we are dropping a circular term of order 2PN and $e_{0}^{0}$, so the unfaithfulness is largely independent of $e_{0}$ (as it should be). The small faithfulness ($M\sim 0.65$) of template A means that the 2PN-$e_{0}^{0}$ term is very important, and it is suggestive of the necessity to include higher PN orders for circular ($e_{0}^{0}$ order) templates. This is well known in the GW data analysis community, and it is indeed implemented in the $3.5$PN accurate circular template TaylorF2 Creighton and Anderson (2011).

IV.2 Systematic errors vs. statistical errors

In Fig. 2 we study the conditions under which systematic errors are smaller than statistical errors for selected truncated templates. Each panel shows the unfaithfulness (left y-axis) and the SNR obtained when we replace the inequality in (22) by an equality (right y-axis) of selected truncated templates as a function of $e_{0}$. Solid lines correspond to the median unfaithfulness (or SNR) over our sample of 1000 binaries. Shaded areas correspond to the 25th and 75th percentiles, so they give an idea of the spread in the data.

Consider, for example, a detection with SNR $\rho=20$ and initial eccentricity $e_{0}=0.02$. If the corresponding point in one of these plots lies below the unfaithfulness-SNR curve of the corresponding template, then systematic errors are negligible with respect to statistical errors. The faithfulness by itself is not sufficient to decide whether systematic errors are negligible: we also need Eq. (22) to determine the maximum SNR beyond which systematic errors dominate. Note also that it would be incorrect to use these plots for low SNRs, since large SNRs were assumed to derive Eq. (22).

Another way to read these plots is as follows. Suppose that we want to compute the posterior distribution for a detection with maximum likelihood corresponding to $e_{0}=2\times 10^{-2}$, $\rho\sim 20$. If we want systematic errors to be negligible with respect to statistical errors when we construct the posterior distribution, we can choose a template whose unfaithfulness-SNR curve at $e_{0}\sim 2\times 10^{-2}$ (as shown in Fig. 2) gives $\rho>20$. For example, templates B, C and D satisfy this criterion, whereas template A does not. Let us remark once again that Fig. 2 should not be trusted for low SNRs, therefore (for example) the unfaithfulness-SNR curve for template A cannot be trusted in this example.

Curves corresponding to low unfaithfulness (or high $\rho$) mean that the GW template will have negligible systematic errors (recall that a point must lie below the unfaithfulness-SNR curve for systematic errors to be negligible). Figure 2 implies that systematic errors become negligible as we move from templates A to D, G to D and H to D, i.e. as we move towards templates which are closer to the fiducial template, and our mismodeling errors become smaller. Furthermore, the ratio of systematic to statistical errors becomes smaller as $e_{0}$ decreases ($\rho$ gets higher as $e_{0}$ decreases). An exception is template A, for which the dominant dropped term is independent of $e_{0}$.

When are systematic errors negligible with respect to statistical errors? Focus, for example, on template H in the bottom-left panel of Fig. 2 and on the two templates in the bottom-right panel. For template H, systematic errors are negligible in the whole range $0<e_{0}<0.1$ when $\rho<10^{3}$. For the $\left\{e_{0}^{4}\right\}$ template, the bottom-right panels shows that systematic errors become negligible in the range $e_{0}<6\times 10^{-2}$ for SNRs below the blue curve. The corresponding range is smaller for the $\left\{e_{0}^{2}\right\}$ template. In these regions the truncated template can be used for parameter estimation to save computational time and the dominant errors are statistical. The convergence of statistical errors will be the topic of the next subsection.

Table 2 – built through the inequality (22) – shows the maximum detection SNR for which systematic errors are smaller than statistical errors for templates obtained by dropping the term corresponding to a given cell, and for selected values of $e_{0}=0.01,\,0.04,\,0.07,\,0.1$. Similarly, Table 3 shows the maximum $e_{0}$ such that detections with $\rho<25$ have systematic errors smaller than statistical errors when terms corresponding to the respective cells in the table are dropped.

IV.3 Convergence of statistical errors

Unfaithfulness-SNR plots can be used to identify templates for which systematic errors are smaller than statistical errors. Given such a template, statistical errors can be computed (in the high-SNR limit) as the square root of the diagonal elements of the inverse Fisher matrix. We computed statistical errors for the 1000 binaries in our catalog using the $\left\{e_{0}^{0}\right\}$, $\left\{e_{0}^{2}\right\}$, $\left\{e_{0}^{4}\right\}$ and $\left\{e_{0}^{6}\right\}$ templates (where the last one is the fiducial template).222Similar calculations were performed in Nishizawa et al. (2016) for two templates: a circular template and a template at leading order in $e_{0}$ with different amplitudes. In Fig. 3 we plot median statistical errors for these four templates.

Statistical errors change going from template $\left\{e_{0}^{0}\right\}$ to template $\left\{e_{0}^{2}\right\}$, but then they plateau. This convergence of statistical errors was not observed in Nishizawa et al. (2016), where errors were computed only for the circular and ${\cal{O}}(e_{0}^{2})$-accurate templates.

The error on the eccentricity $\Delta e_{0}$ decreases when we go from a circular template $\{e_{0}^{0}\}$ to the $\{e_{0}^{2}\}$ template, but it is roughly constant as we increase the order of the $e_{0}$ expansion. This is because the phase of the circular template $\{e_{0}^{0}\}$ is independent of $e_{0}$, so all information comes from the amplitude alone, leading to large errors. The errors $\Delta t_{c}$, $\Delta\Omega_{S}$, and $\Delta\ln D_{L}$ decrease mildly (within a factor of two) as $e_{0}$ increases for all four templates: these are all extrinsic parameters for which measurement information comes largely from the motion of the detector, which is not significantly affected by the template we use. The mass errors $\Delta\mathcal{M}/\mathcal{M}$ and $\Delta\eta/\eta$ are underestimated by a factor of 5-10 when we use the circular template $\{e_{0}^{0}\}$, and they are largely the same for all eccentric templates $\{e_{0}^{n}\}$ with $n=2,\,4,\,6$; in other words, the simplest eccentric template $\{e_{0}^{2}\}$ already contains enough information to estimate mass measurement errors. Note also that most errors (with the exception of $\Delta e_{0}$) vary by at most factor of 2 as functions of $e_{0}$.

The bottom row of Fig. 3 addresses the question: how do statistical errors change as we increase the order of the $e_{0}$ expansion in the phase? The bottom central panel shows that $|\tau|^{-1/2}$ (the volume of the error ellipsoid, as given by the square root of the determinant of the inverse Fisher matrix) decreases going from the $\left\{e_{0}^{0}\right\}$ to the $\left\{e_{0}^{2}\right\}$ template: the 11-dimensional error ellipsoid shrinks (i.e., statistical errors decrease) with more accurate templates. This is in apparent contradiction with previous plots, showing that many errors on individual parameters increase. The solution to this apparent paradox (as shown in the bottom-right panel) has to do with correlations between parameters, as measured by $\epsilon/\epsilon_{0}$ – see the discussion below Eq. (22): this quantity333We slightly modify the definition of the correlation factor to stand for the product of $\Delta{\cal{M}}/{\cal{M}}$, $\Delta{\cal{\eta}}/{\cal{\eta}}$, $\Delta t_{c},\Delta\phi_{c}$, $\Delta{\cal{D}}_{L}/{\cal{D}}_{L}$, $\Delta\Omega_{S}$, $\Delta\Omega_{L}$, $\Delta\beta$ and $|\tau|^{1/2}$. Note that $|\tau|$ is equal to the reciprocal of the parameter $\epsilon_{0}$ introduced in Sec. III, thus quantifying an overall measure of the statistical errors. increases by a factor of $10^{2}-10^{3}$ going from the $\left\{e_{0}^{0}\right\}$ template to the $\left\{e_{0}^{2}\right\}$ template, and then remains roughly constant. Recall that all plots refer to a fixed redshift $z=0.1$ ($D_{L}=447.8$ Mpc), but errors scale linearly with $D_{L}$ in the large-SNR limit.

In Fig. 4 we plot histograms of the statistical errors for all 1000 binaries using our fiducial template. These histograms essentially confirm the conclusions of Nishizawa et al. (2016). We can measure the initial eccentricity as long as $\Delta e_{0}<e_{0}$: this is true for most binaries when $e_{0}>0.1$. If $e_{0}=0.01$, $\Delta e_{0}<e_{0}$ for about $90\%$ of the binaries in our sample.

Before closing this section, we would like to mention that setting $\beta=0$ (as was done in Nishizawa et al. (2016), thereby reducing the number of parameters from 11 to 10) results in a decrease of $\Delta\Omega_{S}$ by a factor of $\sim 10$. This suggests that the parameter $\beta$ is correlated with sky location parameters, and that setting it to zero can lead to an underestimation of those errors.

V Conclusions

We studied the convergence of the frequency domain 2PN-$e_{0}^{6}$ order accurate “fiducial” GW templates for compact eccentric binaries of Tanay et al. (2016). We built truncated templates by dropping certain terms, and assessed the importance of those terms by computing the unfaithfulness (Fig. 1). Dropping most terms leads to unfaithfulness $<0.02$ for $e_{0}<0.1$. The terms that produce the largest unfaithfulness when dropped are the 2PN$-e_{0}^{0}$ and 0PN$-e_{0}^{6}$ terms; extensions at 0PN$-e_{0}^{n}$ with $n>6$ Yunes et al. (2009) and $m$PN$-e_{0}^{0}$ with $m>2$ Blanchet (2014) (e.g., the TaylorF2 approximant) are already available in the literature.

We then investigated the conditions under which truncated templates produce systematic errors which are smaller than statistical errors (Fig. 2). This helps us to identify “fast templates” that can be used for parameter estimation considering only statistical errors.

In Fig. 3 we studied the convergence of statistical errors. Statistical errors converge very quickly, and they do not change much as long as we include terms of order $e_{0}^{2}$ in the phasing. More accurate templates yield larger statistical errors than the $\left\{e_{0}^{0}\right\}$ template for most parameters, with the exception of the error $\Delta e_{0}$ on the initial eccentricity. However the error ellipsoid shrinks as we increase the order of the $e_{0}$ expansion: indeed, statistical errors for most of the individual parameters increase because of the larger correlations between parameters. Figure 4 shows statistical errors for the fiducial template, and it confirms the main conclusions of Ref. Nishizawa et al. (2016) (which used slightly different templates).

Several extensions of this work are possible and necessary. An important limitation of our study is that we kept only the three leading-order harmonics in our templates; future work should further explore the convergence of $\tilde{h}(f)$ [Eq. (2)] as the number of harmonics changes. Our analysis is specific to stellar-origin BH binaries observed with LISA, but similar work should be done for second- and third-generation Earth-based detectors and using other templates (see e.g. Damour et al. (2004); Königsdörffer and Gopakumar (2006)). There are ongoing efforts to extend our “fiducial templates” Tanay et al. (2016) to 1PN order in amplitude and 3PN order in phase, including the effects of periastron advance Tiwari et al. (2019). As soon as these templates are available, an extension of our analysis can be used to assess the relative significance of PN amplitude corrections with respect to phase corrections. It will also be important and useful to extend our study to Fourier-domain templates accurate at 3PN and valid for large eccentricities, which are currently under development Moore et al. (2018); Moore and Yunes (2019).

Acknowledgments

We would like to thank Leo Stein and Kaze W.K. Wong for discussions and suggestions. E.B. is supported by NSF Grant No. PHY-1841464, NSF Grant No. AST-1841358, NSF-XSEDE Grant No. PHY-090003, and NASA ATP Grant No. 17-ATP17-0225. This work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 690904. Computational work was performed at the Mississippi Center for Supercomputing Research (MCSR) and at the Maryland Advanced Research Computing Center (MARCC). A.N. is supported by JSPS KAKENHI Grants No. JP17H06358 and No. JP18H04581.

Appendix A Oscillatory cross-terms in the Fisher matrix

As discussed in Sec. II, the cross terms in the integrand of Eq. (18) are highly oscillatory, and thus can be neglected. This can be understood analytically as follows. Let us first truncate the templates of Tanay et al. (2016) at leading order in both the PN parameter $x$ and the initial eccentricity $e_{0}$ (i.e., we consider circular templates). We first decompose the template into its first six harmonics:

[TABLE]

and then we decompose each harmonic into an amplitude and a phase to get

[TABLE]

By Eq. (18), the Fisher matrix elements involve derivatives of the template with respect to the parameters:

[TABLE]

and they can be broken down into a sum of integrals of the form

[TABLE]

as a result of the coupling of different harmonics $(j_{1},\,j_{2})$ and different parameters $(p_{A},\,p_{B})$. The phase term in the exponential reads

[TABLE]

Therefore the integrand has a rapidly oscillatory phase $\Delta\Psi\propto f^{-5/3}$ whenever $j_{1}\neq j_{2}$. For diagonal elements of the Fisher matrix ($p_{A}=p_{B}$) terms with $j_{1}=j_{2}$ always exist, and those integrals dominate. For certain off-diagonal terms (e.g. the $\ln D_{L}$-$\phi_{c}$ term), the integrals with $j_{1}=j_{2}$ exactly vanish. However, the natural scale for these terms is set by the corresponding diagonal elements of the Fisher matrix, and the integrals with $j_{1}\neq j_{2}$ can still be neglected.

Appendix B Beam pattern functions and other quantities appearing in the templates

In this appendix, for completeness, we define certain quantities appearing in the GW strain via [Eqs. (2), (3), (4) and (5)]. It is well known that certain combinations of trigonometric functions of the eccentric anomaly $u$ of a binary with eccentricity $e$ can be written as Fourier-Bessel series Yunes et al. (2009); Maggiore (2007):

[TABLE]

where $J_{k}$ denotes Bessel functions of the first kind

[TABLE]

and $\Gamma$ is the Gamma function. When combined with the well-known relations involving the orbital phase $\phi$ of a Keplerian orbit

[TABLE]

the equations above yield

[TABLE]

Following Yunes et al. (2009), the plus and cross polarizations can be written as

[TABLE]

where $\mu=m_{1}m_{2}/(m_{1}+m_{2})$ and $p$ (the semilatus rectum of the orbit) is related to the orbital angular frequency $\omega$ via

[TABLE]

where the inclination angle ${\iota}$ is defined by $\cos\iota=\hat{L}\cdot\hat{N}$ Martel and Poisson (1999); Yagi and Tanaka (2010), where $\hat{L}$ and $\hat{N}$ are unit vectors in the direction of the orbital angular momentum and in the direction of the source, respectively. By plugging Eqs. (34) and (35) into the expressions for $h_{+}$ and $h_{\times}$ above, we get

[TABLE]

The quantities $C_{+,\times}^{j}$ and $S_{+,\times}^{j}$ read Yunes et al. (2009)

[TABLE]

where $c_{2\beta}=\cos{2\beta},~{}s_{2\beta}=\sin{2\beta},~{}c_{i}=\cos{\iota}$ and $s_{i}=\sin{\iota}$. The GW strain at the detector is the linear combination $h(t)=F_{+}h_{+}+F_{\times}h_{\times}$, and its Fourier transform is given by Eq. (2). The beam pattern functions $F_{+}$ and $F_{\times}$ can be found, e.g., in Berti et al. (2005); Yagi and Tanaka (2010).