Characterization of numerical relativity waveforms of eccentric binary black hole mergers

TL;DR

This paper presents a new method to characterize eccentric binary black hole mergers in numerical relativity, demonstrating circularization tendencies and emphasizing the importance of higher-order waveform modes for accurate gravitational wave modeling.

Contribution

It introduces a robust technique for quantifying eccentricity and other parameters in numerical relativity simulations, and analyzes the significance of higher-order modes in waveform accuracy.

Findings

Moderately eccentric mergers circularize within 15 cycles before merger.

Higher-order modes significantly improve waveform overlap accuracy.

Eccentricity and mode inclusion impact gravitational wave detection and modeling.

Abstract

We introduce a method to quantify the initial eccentricity, gravitational wave frequency, and mean anomaly of numerical relativity simulations that describe non-spinning black holes on moderately eccentric orbits. We demonstrate that this method provides a robust characterization of eccentric binary black hole mergers with mass-ratios $q\leq10$ and eccentricities $e_0\leq0.2$ fifteen cycles before merger. We quantify the circularization rate of a variety of eccentric numerical relativity waveforms introduced in [1] by computing overlaps with their quasi-circular counterparts, finding that $50M$ before merger they attain overlaps ${\cal{O}}\geq0.99$, furnishing evidence for the circularization of moderately eccentric binary black hole mergers with mass-ratios $q\leq10$. We also quantify the importance of including higher-order waveform modes for the characterization of eccentric binary black hole mergers. Using two types of numerical waveforms, one that includes $(\ell, \, |m|)= \{(2,\,2),\, (2,\,1),\, (3,\,3),\, (3,\,2), \, (3,\,1),\, (4,\,4),\, (4,\,3),\, (4,\,2),\,(4,\,1)\}$ and one that only includes the $\ell=|m|=2$ mode, we find that the overlap between these two classes of waveforms is as low as ${\cal{O}}=0.89$ for $q=10$ eccentric binary black hole mergers, underscoring the need to include higher-order waveform modes for the description of these gravitational wave sources. We discuss the implications of these findings for future source modeling and gravitational wave detection efforts.

Tables

Table 1. Table 1: The table presents the mass-ratio, q 𝑞 q , and the measured values of initial eccentricity, mean anomaly and dimensionless orbital frequency, ( e ^ 0 , l ^ 0 , x ^ 0 ) subscript ^ 𝑒 0 subscript ^ 𝑙 0 subscript ^ 𝑥 0 (\hat{e}_{0},\,\hat{l}_{0},\,\hat{x}_{0}) , at a time, t 0 subscript 𝑡 0 t_{0} , when the numerical relativity simulations are free from junk radiation.

Simulation	$q$	${\hat{e}}_0$	${\hat{l}}_0$	${\hat{x}}_0$
E0001	1.0	0.052	3.0	0.0770
K0006	4.0	0.068	3.0	0.0826
L0009	4.5	0.052	3.0	0.0839
L0016	5.0	0.140	2.9	0.0862
P0001	6.0	0.050	3.0	0.0867
P0017	8.0	0.060	3.0	0.0927
P0006	8.0	0.080	2.9	0.0931
P0007	8.0	0.100	2.9	0.0926
P0009	10.0	0.060	2.9	0.0971
P0022	10.0	0.080	2.9	0.0979
P0023	10.0	0.120	2.9	0.0968
P0024	10.0	0.180	3.0	0.0957

Equations

$$h^{(\ell =|m|=2)}(t)$$

$$M\omega$$

$$\left(h(t)|s(t)\right)$$

$$\left(\hat{h}(t)|\hat{s}(t)\right)$$

$$\mathcal{O}\left(h,s\right)$$

$$h(t,\theta ,\varphi )=h_{+}-\text{i}{h}_{\times }=\sum _{\ell \ge 2}\sum _{m\ge -\ell }^{m\le \ell }{h}^{\ell m}{}_{-2}{Y}_{\ell m}\left(\theta ,\varphi \right),$$

Full text

I Introduction

Numerical relativity (NR) has played a key role in the discovery and interpretation of gravitational wave (GW) observations Abbott et al. (2016a, b, 2017a, 2017b, 2017c); The LIGO Scientific Collaboration et al. (2018); Abbott et al. (2017d). As the LIGO Abbott et al. (2016c); The LIGO Scientific Collaboration et al. (2015) and Virgo Acernese et al. (2015) observatories continue to probe the GW spectrum, NR will provide key insights to infer the properties of binary black hole (BBH) mergers whose GWs exhibit strong spin-precession, or high-order waveform modes Chu et al. (2016); Mroué et al. (2013); Kumar et al. (2016); Abbott et al. (2016d); The LIGO Scientific Collaboration et al. (2017); Lange et al. (2017); Huerta et al. (2018); Rebei et al. (2018); Hinder et al. (2010); Rosofsky et al. (2018); Huerta et al. (2017); Sperhake et al. (2008); Hinder et al. (2018); Ott (2009); Kotake (2013); Hinderer et al. (2018); Foucart et al. (2018); Radice et al. (2018). These improved studies will be key to infer the formation channels of these objects, and to ascertain whether they are accurately described by general relativity.

BBH systems that form through massive stellar evolution in the field of galaxies are expected to enter the frequency band of LIGO-type detectors with nearly circular orbits Postnov and Yungelson (2014). Under this assumption, NR groups have produced thousands of NR waveforms to get insights into the physics of these GW sources Mroué et al. (2013); Boyle et al. (2019); Healy et al. (2019, 2017a); Jani et al. (2016). These NR waveforms have been used to calibrate semi-analytical models Bohé et al. (2017); Taracchini et al. (2014); Buonanno et al. (2009); Pan et al. (2010); Cotesta et al. (2018); London et al. (2018); Husa et al. (2016); Khan et al. (2016); Huerta et al. (2017), produce fast interpolators using gaussian process regression Huerta et al. (2018); Varma et al. (2019), and surrogate models Blackman et al. (2015, 2017) to inform the development of signal-processing tools for GW searches Usman et al. (2016); Sachdev et al. (2019); George and Huerta (2018, 2018); George et al. (2017); Shen et al. (2017); Allen et al. (2019); Shen et al. (2019); George and Huerta (2017); Wei and Huerta (2019); Rebei et al. (2018), and more recently to directly infer the astrophysical properties of BBH mergers through GW observations Abbott et al. (2016d); Kumar et al. (2018); Lange et al. (2017).

BBHs may also form in dense stellar environments, such as globular clusters and galactic nuclei Samsing and Ramirez-Ruiz (2017); Samsing et al. (2014); Samsing (2018); Leigh et al. (2018); Samsing et al. (2018a, b, c); Randall and Xianyu (2018a); Huerta and Gair (2009); Samsing et al. (2017); Samsing and Ilan (2018); Huerta et al. (2015); Samsing and Ilan (2019); Samsing (2018); Samsing et al. (2018d); Huerta et al. (2014); Antonini et al. (2014); Samsing and D’Orazio (2018a); D’Orazio and Samsing (2018); Samsing and D’Orazio (2018b); Zevin et al. (2018); Rodriguez et al. (2018a); Kremer et al. (2018); Lopez et al. (2018a); Hoang et al. (2018); Gondán et al. (2017); Hoang et al. (2017); Randall and Xianyu (2018b); Mikóczi et al. (2012); Naoz et al. (2013); Gondán and Kocsis (2018); Antonini and Rasio (2016); Huerta and Brown (2013); Arca-Sedda et al. (2018); Takács and Kocsis (2018); Gondán et al. (2018); Antonini et al. (2018, 2016); Randall and Xianyu (2018a, b); Samsing et al. (2019); Lopez et al. (2018b); Rodriguez et al. (2018b); Zevin et al. (2019). Electromagnetic observations provide evidence for their existence in galactic clusters, and in the center of the Milky Way Hailey et al. (2018); Sippel and Hurley (2013); Strader et al. (2012); Samsing et al. (2018b). This increasing body of observational evidence has sparked the interest of the community to better understand these sources. This program includes GW source modeling, formation channels and merger rates, and astrophysically motivated scenarios where GW observations may be used to confirm or rule out the existence of BBHs in dense stellar environments. On this latter point, it is now widely accepted that the measurement of orbital eccentricity through GW observations would be the cleanest signature to furnish evidence for the existence of compact binary populations in dense stellar environments. It is for this reason that the GW source modeling community is sharpening its analytical and numerical tools to infer the imprints of orbital eccentricity in GW searches. East et al. (2013); Cao and Han (2017); Hinderer and Babak (2017); Huerta et al. (2017); Huerta and Gair (2009); Huerta et al. (2015, 2018, 2014); Moore et al. (2018); Loutrel and Yunes (2017a, b); Yunes et al. (2009).

An accurate description of the physics of eccentric BBHs throughout the late-inspiral, merger and ringdown requires NR Hinder et al. (2008, 2010); East et al. (2012); East and Pretorius (2012); Gold et al. (2012); Gold and Brügmann (2013); Paschalidis et al. (2015); East et al. (2016); Radice et al. (2016); The LIGO Scientific Collaboration et al. (2016); Lewis et al. (2017); Huerta et al. (2019). Once the data products of NR simulations are post-processed, and NR waveforms are extracted Johnson et al. (2018), it is necessary to characterize them, i.e., we need to quantify the eccentricity and other orbital parameters that uniquely identify them. One can address this task using a variety of methods. If BBHs are on nearly quasi-circular orbits, then one could use the approach introduced in Healy et al. (2017b), which combines information about the orbital separation of the BHs, and waveform phase and amplitude of the Weyl scalar $\psi_{4}$. This methodology only includes ${\cal{O}}(e)$ corrections to measure orbital eccentricity, which limits its applicability to characterize moderately eccentric BBH mergers. Some other methods try to infer orbital eccentricity based on the trajectories of the BHs in the NR simulation, which is not a sound approach given that these trajectories are gauge-dependent.

We circumvent the aforementioned limitations by introducing a gauge-invariant method that characterizes a NR waveform by comparing to a large array of ENIGMA waveforms Huerta et al. (2018). The inspiral-merger-ringdown ENIGMA model consists of an inspiral evolution that encodes higher-order post-Newtonian (PN) corrections to the motion of compact sources on eccentric orbits, combined with self-force and BH perturbation corrections. This approach ensures that the dynamics of quasi-circular and moderately eccentric BBHs are accurately described. Assuming that moderately eccentric systems circularize prior to merger, we attach a stand-alone quasi-circular merger waveform to the inspiral evolution. The merger waveforms are produced using a Gaussian Process Emulator Moore et al. (2016) that is trained with NR waveforms describing quasi-circular BH mergers.

In summary, our goal is to constrain the PN parameters that produce the optimal overlap between ENIGMA waveforms and their NR counterparts. This study is timely and relevant if we are to use NR waveforms to properly characterize future observations of eccentric BBH mergers. Once we showcase the application of this method, we also quantify the circularization of eccentric BBHs near merger, and quantify the impact of higher-order waveform modes in the morphology of eccentric NR waveforms.

This paper is organized as follows. In Section II we describe how we adapted the ENIGMA model to characterize eccentric NR waveforms, and present results of this method in Section III. In Section IV we study the circularization of moderately eccentric BBH mergers that retain eccentricity a few cycles before merger. We use overlap calculations to quantify the impact of higher-order waveform modes on the morphology of eccentric BBH mergers in Section V. We summarize our findings and future directions of work in Section VI.

II Methods

We measure the orbital eccentricity of NR waveforms using the $\ell=\lvert m\rvert=2$ mode. The rationale for this choice is that we are characterizing NR waveforms using the ENIGMA waveform model, which only includes the $\ell=\lvert m\rvert=2$ mode.

We explored a variety of gauge-invariant objects to directly compare NR and ENIGMA waveforms, and found that the dimensionless object $M\omega$, where $\omega$ is the mean orbital frequency, and $M$ stands for the total mass of the BBH, provides a robust approach to capture the signatures of eccentricity. To compute $M\omega$ using our NR waveforms, we first extract the waveform modes $h_{\ell,\,m}$ at future null infinity using the open source software POWER Johnson et al. (2018). As mentioned above, since we will be comparing $M\omega$ with a waveform model that only includes the $\ell=\lvert m\rvert=2$ mode, we take only the waveform mode $h^{(\ell=\lvert m\rvert=2)}(t)$ to compute $M\omega$ using the following relations:

[TABLE]

where $\dot{\phi}^{(\ell=\lvert m\rvert=2)}$ is the unwrapped GW phase of $h^{(\ell=\lvert m\rvert=2)}(t)$. On the other hand, the quantity $M\omega$ is one of the building blocks of the ENIGMA model, which includes tail, tails-of-tails and non-linear memory corrections at 1.5PN, 2.5PN and 3PN order Huerta et al. (2017, 2018). To be self-contained, it is worth mentioning that the mean orbital frequency, $\omega$, is related to the mean motion, $n$, through the relation $\omega={\cal{K}}n$, where ${\cal{K}}$ is the periastron precession. Furthermore, the mean anomaly, $l$, is related to the mean motion, $n$, through the relation $M\dot{l}=Mn$ Blanchet (2006).

II.1 Optimization algorithm

In order to characterize an NR waveform, we compare it to ENIGMA waveforms generated algorithmically. Initial conditions of the ENIGMA simulation are varied until the ENIGMA time evolution of $M\omega$ agrees with the NR evolution to a specified degree of accuracy. The ENIGMA parameters of concern are initial orbital eccentricity $e_{0}$, initial GW frequency $f_{0}$, and mean anomaly $l_{0}$. This method is schematically described in Figure 1. Following Königsdörffer and Gopakumar (2006), we parameterize the orbital eccentricity, $e$, in the ENIGMA model using the PN time eccentricity parameter, i.e., $e^{\textrm{PN}}_{t}\rightarrow e$.

The NR waveform is first preprocessed by cutting off initial noise (junk radiation). We then apply a Savitsky-Golay filter to remove high frequency noise from the $M\omega$ time-series data. We have used the Python implementation of the Savitsky-Golay filter, which is extensively used to remove high frequency noise from data while preserving the original features of the signal better than other types of filtering approaches, such as moving averages techniques SciPy Cookbook . Applying this filter also allows for more accurate detection of local extrema in the signal. Figure 2 presents the input and output time-series data when we use this filter to remove high-frequency noise from the $M\omega$ time-series data.

From a modeling perspective, we know that both eccentricity and gravitational wave frequency are coupled (this is the reason why many semi-analytical approaches have been explored in recent years using the stationary phase approximation Moore and Yunes (2019); Moore et al. (2018); Huerta et al. (2014); Yunes et al. (2009)). In view of this observation, we optimize $(e_{0},f_{0})$ simultaneously. On the other hand, we know that the net effect of varying $l_{0}$ is to introduce what amounts to a minor shift in the position of the peaks of the orbital frequency $M\omega$. To reduce computational time, we perform a separate, coarser search for $l_{0}$ instead of optimizing it simultaneously with the other parameters. Since we know experimentally that values close to $l_{0}\sim\pi$ will produce ENIGMA waveforms that tend to be aligned with their NR counterparts, we start the optimization with this seeded value for $l_{0}$ and then do minor refinement to this value later on.

Since input parameters are required to produce ENIGMA waveforms, seed values are initialized for $f_{0}$ and $l_{0}$. We provide an informed guess of GW frequency using the relation $\omega_{0}/\left(M\pi\right)$, where $\omega(t_{0})=\omega_{0}$ and $t_{0}$ is the time at which the NR waveform is free from junk radiation. Mean anomaly is initialized to $\pi$, a value manually determined to be optimal through verification of a few individually sampled NR waveforms. Orbital eccentricity does not require a seed value since the range of possible values is consistent for all catalogued waveforms, so the following grid search directly samples eccentricity from a predefined range.

The algorithm starts with a grid search in the 2D parameter space of $(f_{0},e_{0})$, and iteratively refines the resulting parameter guess. To generate the grid, we densely sample the frequency range $f\in[f_{0}-5\textrm{Hz},\,f_{0}+5\textrm{Hz}]$ and the eccentricity range $e_{0}\in[0.1,\,0.3]$. For each coordinate pair, an ENIGMA $M\omega$ is produced using the specified $(f_{0},\,e_{0})$ values and the seeded $l_{0}$. The resulting $M\omega$ time evolution is then compared to that of the original NR waveform. Parameters are chosen that minimize two cost functions. Specifically, the optimal coordinate pair is that which minimizes the two cost functions $|t^{FC}_{\texttt{{ENIGMA}}}-t^{FC}_{NR}|$ and $|A^{*}_{\texttt{ENIGMA}}-A^{*}_{NR}|$, where $t^{FC}$ is the time duration of the first orbital cycle (time between two consecutive maxima, as shown in Fig. 3) and $A^{*}$ is the amplitude (difference in $M\omega$ from first maximum to first minimum, as shown in Fig. 3) respectively. Because the eccentricity is most easily measurable early on in the waveform time evolution, only the first orbital cycle is considered in these cost functions. Throughout this grid search, $l_{0}$ is held constant.

After completing the grid search, the chosen $(f_{0},\,e_{0})$ parameters are further refined iteratively using a hill-climbing approach. In this stage, the initial GW frequency and orbital eccentricity are independently varied stepwise. In order to avoid the hill-climbing problem of finding a solution at a local minimum of the cost function, the initial guesses of $(f_{0},\,e_{0})$ are randomly sampled within the range confined by the grid search, and the greedy search is repeated until parameters are found that vary by $\leq 1\%$. This “greedy” search is run for a preset number of iterations. At each iteration, an ENIGMA waveform is generated using the current $(f_{0},\,e_{0})$ guesses and the seeded $l_{0}$. Properties of the $M\omega$ evolution are then compared between the ENIGMA and NR models to determine how to increment $e_{0}$ and $f_{0}$. Initial eccentricity is evaluated depending on the amplitude of the first orbital cycle $A^{*}$, calculated again as the difference in $M\omega$ from first maximum to first minimum. If $A^{*}_{\texttt{ENIGMA}}>A^{*}_{\textrm{NR}}$, then the current $e_{0}$ guess is too large, so $e_{0}$ is decreased; vice versa is true for $A^{*}_{\texttt{ENIGMA}}<A^{*}_{\textrm{NR}}$. Initial frequency is varied based on the time $t^{*}$ each simulation reaches the threshold frequency, which corresponds to the time at which a quasi-circular waveform would be attached in the ENIGMA model. If $t^{*}_{\texttt{ENIGMA}}-t^{*}_{\textrm{NR}}>0$, the current $f_{0}$ guess is too low, so $f_{0}$ is increased; conversely, if $t^{*}_{\texttt{ENIGMA}}-t^{*}_{\textrm{NR}}<0$, the current $f_{0}$ guess is decreased. To increase parameter precision, increment step size is one degree of precision higher than that of the grid search resolution. Using as input seeds the values found for $(f_{0},\,e_{0})$, we constrain $l_{0}$ using an additional grid search. This second grid search finds the $l_{0}$ (sampled from the range $l_{0}\in[0,\,2\pi)$) minimizing the cost function $|t^{M}_{\textrm{ENIGMA}}-t^{M}_{\textrm{NR}}|$, where $t^{M}$ is the difference in time of first orbital cycle occurrence between the ENIGMA and NR waveforms.

Upon finding the optimal values $(\hat{f}_{0},\,\hat{e}_{0},\,\hat{\ell}_{0})$, we recast the initial parameter into one that is more commonly used to describe NR waveforms, i.e., $\hat{x}_{0}=\left(M\omega_{0}\right)^{2/3}$. In Figure 3 we present results of this optimization procedure for a sample of NR waveforms. Note that the $M\omega$ time-series datasets presented therein have already been cleaned from high-frequency noise using a Savitsky-Golay filter, as shown in Figure 2. We used this method to characterize the 89 NR waveforms presented in Huerta et al. (2019). The properties of the BBHs considered in this study are summarized in Table 1.

III Characterization of eccentric numerical relativity waveforms

To demonstrate that our optimization algorithm provides an accurate characterization of eccentric NR waveforms, in Figure 4 we present overlap calculations between ENIGMA waveforms using the corresponding optimal triplet $(\hat{e}_{0},\,\hat{l}_{0},\,\hat{x}_{0})$ and their NR counterparts. We note that the overlap between the two classes of waveforms is ${\cal{O}}\geq 0.95$. An important finding of these results is that ENIGMA is correctly characterizing NR waveforms of moderately eccentric BBH mergers that have rather asymmetric mass-ratios.

Given that ENIGMA has only been validated, as opposed to calibrated, with NR waveforms, these results indicate that it is possible to capture the physics of eccentric compact binary systems by combining in a consistent manner results from several analytical relativity approaches, such as PN and BH perturbation theory and the self-force program. Much work still remains to be done to provide a robust framework to extend these formalisms to accurately describe the inspiral-merger-ringdown of highly eccentric systems.

As we mentioned above, ENIGMA was constructed under the assumption that moderately eccentric BBH systems circularize prior to merger. When we constructed ENIGMA, we determined the transition point between inspiral and merger by constraining the time window before merger within which state-of-the-art, inspiral-merger-ringdown quasi-circular waveforms Bohé et al. (2017) and quasi-circular ENIGMA waveforms have overlaps ${\cal{O}}\geq 0.99$. From this time window, we selected the attachment time closest to merger. At the time of that study we had produced eccentric NR waveforms with mass-ratios $q\leq 5.5$, which we used to validate this approach Huerta et al. (2018). Having constructed an NR waveform catalog that now covers a deeper parameter space Huerta et al. (2019), we can actually quantify the circularization rate of more asymmetric mass-ratio BBHs that retain non-negligible eccentricities just a few cycles before merger. This is discussed in the following section.

IV Circularization of eccentric binary black hole systems

To quantify the circularization of moderately eccentric BBH mergers, we compare them directly to quasi-circular BBH mergers that have identical mass-ratios. In practice, and assuming a flat power spectral density, we compute the inner product between eccentric NR waveforms, $h(t)=h_{+}-ih_{\times}$, and their quasi-circular counterparts, $s(t)=s_{+}-is_{\times}$, using the relation

[TABLE]

where $t_{0}$ is a fiducial time from which we compute the overlap, and $T$ corresponds to the merger time. We then compute the normalized overlap

[TABLE]

Finally, the quantity we quote for our results below is the maximized overlap

[TABLE]

which is obtained by maximizing the normalized overlap over time and phase of coalescence, $(t_{c}\,,\phi_{c})$, respectively.

We have selected a variety of scenarios to illustrate how mass-ratio and initial eccentricity, $e_{0}$, drive the circularization of eccentric BBHs. Figure 5 indicates that for the more eccentric systems (see Table 1), circularization only happens about $50M$ before merger. We also notice that the increase in overlap, as $t_{0}\rightarrow 0$ in Eq. (4), is not monotonic. Rather, it has an oscillatory behavior that tracks the eccentric trajectory of the BBH system, and which is clearly captured by the waveform amplitude. We have included the waveform amplitude of the eccentric and quasi-circular signals in the panels of Figure 5 to clearly show this finding. Notice that as soon as the waveform amplitude of the eccentric signal becomes increasingly monotonic, so does the overlap.

In view of this analysis, we conclude that modeling eccentric BBH mergers under the assumption that they circularize prior to merger requires two key components: (i) an effective scheme that describes the early inspiral evolution, and which remains accurate at least $50M$ before merger; (ii) a stand-alone quasi-circular merger that can be smoothly attached to the late-inspiral evolution. To accomplish this, it is key that future models go beyond low-order PN approximation to describe the radiative and conservative pieces of the waveform dynamics Huerta et al. (2017, 2018); Moore and Yunes (2019); Loutrel et al. (2019).

V Higher-order waveform modes

The impact of higher-order waveform modes for GW detection in terms of signal-to-noise (SNR) calculations has been explored in Rebei et al. (2018). In this section we now quantify the impact of including the modes $(\ell,\,\lvert m\rvert)=\{(2,\,2),\,(2,\,1),\,(3,\,3),\,(3,\,2),\,(3,\,1),\,(4,\,4),\,(4,\,3),\,(4,\,2)$, $(4,\,1)\}$ in terms of overlap calculations. In practice, what we do is to construct two types of signals, one that includes the $\ell=\lvert m\rvert=2$ mode only, and one that includes all the $(\ell,\,\lvert m\rvert)$ modes listed above using the relation

[TABLE]

where ${}_{-2}Y_{\ell m}\left(\theta,\phi\right)$ are the spin-weight–2 spherical harmonics Blanchet (2006), and the reference frame described by $(\theta\,,\phi)$, fixed at the center of mass of the BBH, determines the location of the GW detector. We construct higher-order mode NR waveforms using Eq. (7), and follow  Rebei et al. (2018) to determine the $(\hat{\theta}\,,\hat{\phi})$ combinations that maximize the contribution of $(\ell,\,\lvert m\rvert)$ modes for GW detection. In practice, we densely sampled the $(\theta\,,\phi)$ parameter space, and constrained the regions where the integrated amplitude of the $(\ell,\,\lvert m\rvert)$ waveforms is larger than that of their $\ell=\lvert m\rvert=2$ counterparts.

Upon constructing NR waveforms with the optimal $(\hat{\theta}\,,\hat{\phi})$ combinations, we compute their overlap using Eq. (6), but now setting $t_{0}$ to the time at which the NR waveform is free from junk radiation, and $T$ to the final time sample of the NR waveform (see Eq. (3)). We present results for these calculations for a variety of astrophysically motivated scenarios in Figure 6. These results indicate that the inclusion of higher-order modes does not quantitatively modify the morphology of $\ell=\lvert m\rvert=2$ NR waveforms that describe equal-mass, eccentric BBH mergers. However, NR waveforms that describe asymmetric mass-ratios, eccentric BBH mergers have a much richer topology that requires the inclusion of $(\ell,\,\lvert m\rvert)$.

Note that the overlap results shown in the four panels in Figure 6 have been produced using NR waveforms of different lengths. The key point to extract from these analyses is that searches for eccentric BBH mergers that have asymmetric mass-ratios will require signal-processing tools that include higher-order waveform modes. The construction of such algorithms must be pursued in the near future.

VI Conclusions

In this article we introduced a method to characterize numerical relativity waveforms that describe non-spinning black holes on moderately eccentric orbits. To do this, we construct a catalog of ENIGMA waveforms, and then we sift through them until we find the optimal combination of parameters that produce an ENIGMA waveform that best replicates a given NR waveform. Through this procedure, we optimize three orbital parameters, namely, initial orbital eccentricity, initial mean anomaly, and initial dimensionless orbital frequency. We have demonstrated that when we compute the overlap between our optimized ENIGMA waveforms and their NR counterparts, we obtain overlaps ${\cal{O}}\geq 0.95$.

We also quantified the circularization rate of $\ell=\lvert m\rvert=2$ eccentric NR waveforms by computing the overlap between these signals and their quasi-circular counterparts. By choosing a variety of representative systems from our NR catalog Huerta et al. (2019), we found that all NR waveforms, with eccentricities $e_{0}\leq 0.2$ fifteen cycles before merger, circularize at least $50M$ before merger. These findings have a variety of implications for ongoing source modeling efforts. For instance, modeling eccentric BBH mergers under the assumption of circularization prior to merger would require a scheme that accurately describes the effects of eccentricity during the early inspiral evolution, and which also remains accurate deep into the strong-field regime up to just a few tens of $M$ before merger. We expect that ongoing developments in PN theory and in the self-force program will provide the required elements to further enhance the accuracy of existing waveform models to accomplish this goal.

Finally, we explored the need to include higher-order waveform modes to accurately describe the waveform morphology of eccentric BBH mergers. In previous studies, we quantified the importance of $(\ell,\,\lvert m\rvert)$ modes for GW detection in terms of SNR calculations. In this study we have broadened that initial approach, showing that the inclusion of $(\ell,\,\lvert m\rvert)$ modes is essential for an accurate description of asymmetric mass-ratio, eccentric BBH mergers.

Having completed these studies, it is now in order to start developing in earnest NR waveform catalogs that describe spinning BBHs on eccentric orbits, and assess the interplay of spin and eccentricity in the dynamical evolution of these GW sources. Extracting observable signatures from their NR waveforms will inform future GW searches that may confirm or rule out the existence of these type of compact binary populations in dense stellar environments.

VII Acknowledgements

This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the State of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. We acknowledge support from the NCSA and the SPIN Program at NCSA. We thank the NCSA Gravity Group for useful feedback. NSF-1550514, NSF-1659702, NSF-OAC1659702, and TG-PHY160053 grants are gratefully acknowledged.