Physics of eccentric binary black hole mergers: A numerical relativity perspective

TL;DR

This paper presents a catalog of 89 numerical relativity waveforms for eccentric binary black hole mergers, providing insights into their physics and aiding gravitational wave detection efforts.

Contribution

It introduces a new catalog of waveforms for non-spinning black hole binaries with high eccentricity, enhancing understanding of their gravitational wave signatures.

Findings

Quantified energy and angular momentum loss during mergers

Analyzed properties of black hole remnants including mass, spin, and recoil velocity

Discussed implications for gravitational wave detection and modeling

Abstract

Gravitational wave observations of eccentric binary black hole mergers will provide unequivocal evidence for the formation of these systems through dynamical assembly in dense stellar environments. The study of these astrophysically motivated sources is timely in view of electromagnetic observations, consistent with the existence of stellar mass black holes in the globular cluster M22 and in the Galactic center, and the proven detection capabilities of ground-based gravitational wave detectors. In order to get insights into the physics of these objects in the dynamical, strong-field gravity regime, we present a catalog of 89 numerical relativity waveforms that describe binary systems of non-spinning black holes with mass-ratios $1\leq q \leq 10$, and initial eccentricities as high as $e_0=0.18$ fifteen cycles before merger. We use this catalog to quantify the loss of energy and angular momentum through gravitational radiation, and the astrophysical properties of the black hole remnant, including its final mass and spin, and recoil velocity. We discuss the implications of these results for gravitational wave source modeling, and the design of algorithms to search for and identify eccentric binary black hole mergers in realistic detection scenarios.

Tables

Table 1. Table 1: The table presents, from left to right, the mass-ratio of binary black hole mergers, q 𝑞 q , the recoil velocity of quasi-circular binary black hole mergers, | v | kick e 0 → 0 subscript superscript 𝑣 → subscript 𝑒 0 0 kick \lvert v\rvert^{e_{0}\rightarrow 0}_{\textrm{kick}} , and the minimum and maximum recoil velocities of our numerical relativity catalog for a given mass-ratio population, [ | v | kick min , | v | kick max ] subscript superscript 𝑣 min kick subscript superscript 𝑣 max kick [\lvert v\rvert^{\texttt{min}}_{\textrm{kick}},\,\lvert v\rvert^{\texttt{max}}_{\textrm{kick}}] .

$q$	${\left|v\right|}_{\text{kick}}^{e_0\to 0}\left[\text{km}/\text{s}\right]$	$[{\left|v\right|}_{\text{kick}}^{\text{min}},{\left|v\right|}_{\text{kick}}^{\text{max}}]\left[\text{km}/\text{s}\right]$
1.0	0.0	[0.0, 0.0] [0.0, 0.0]
1.5	107.4	[94.6, 101.8] [104.0, 109.9]
2.0	156.7	[136.9, 149.0] [115.1, 157.5]
2.5	173.5	[152.8, 165.0] [132.9, 180.7]
3.0	174.1	[161.4, 170.9] [151.8, 179.1]
3.5	167.1	[154.4, 173.7] [144.0, 178.8]
4.0	156.9	[143.1, 166.8] [140.0, 173.8]
4.5	145.6	[137.9, 154.6] [133.6, 157.1]
5.0	134.4	[121.0, 137.1] [117.6, 136.8]
6.0	113.9	[99.5, 121.0] [104.5, 121.3]
8.0	82.7	[88.7, 96.5] [70.6, 90.4]
10.0	61.9	[47.0, 78.6] [53.9, 66.4]

Table 2. Table 2: ( e 0 , ℓ 0 , x 0 ) subscript 𝑒 0 subscript ℓ 0 subscript 𝑥 0 (e_{0},\,\ell_{0},\,x_{0}) represent the measured values of eccentricity, mean anomaly, and dimensionless orbital frequency parameters. These quantities are computed upon removing the first 60 ​ M 60 𝑀 60M of evolution of the numerical relativity waveforms, as described in  Habib and Huerta ( 2019 ) .

Simulation	$q$	$e_0$	${\ell }_0$	$x_0$
E0001	1	0.052	3.0	0.0770
E0009	2	0.052	3.0	0.0794
E0013	2.5	0.050	3.0	0.0813
E0017	3	0.050	3.0	0.0831
F0002	1	0.066	3.0	0.0780
F0010	2	0.066	3.0	0.0803
F0014	2.5	0.068	3.0	0.0822
F0018	3	0.068	3.0	0.0842
G0003	1	0.094	3.0	0.0788
H0004	1	0.140	3.0	0.0826
H0011	2	0.092	3.0	0.0795
H0015	2.5	0.094	3.0	0.0812
H0019	3	0.094	3.0	0.0832
I0004	1	0.140	3.0	0.0765
I0012	2	0.140	3.0	0.0791
I0016	2.5	0.140	3.0	0.0811
I0020	3	0.140	3.0	0.0824
I0028	4	0.140	2.9	0.0865
J0005	1.5	0.050	3.0	0.0779
J0006	1.5	0.064	3.0	0.0782
J0007	1.5	0.100	3.1	0.0762
J0008	1.5	0.140	3.0	0.0768
J0037	1	0.058	3.0	0.0768
J0038	1	0.076	3.0	0.0762
J0039	1	0.120	3.1	0.0749
J0040	1	0.160	3.0	0.0761
J0041	1.5	0.056	3.0	0.0777
J0042	1.5	0.074	3.0	0.0771
J0043	1.5	0.120	3.1	0.0756
J0044	1.5	0.160	2.9	0.0778
J0045	2	0.056	3.0	0.0793
J0046	2	0.076	3.0	0.0787
J0047	2	0.100	3.0	0.0778
J0048	2	0.160	2.9	0.0794
J0049	2.5	0.058	3.0	0.0811
J0050	2.5	0.078	3.0	0.0806
J0051	2.5	0.120	3.0	0.0795
J0052	2.5	0.160	2.9	0.0817
J0053	3	0.058	3.0	0.0829
J0054	3	0.080	3.0	0.0823
J0055	3	0.120	3.0	0.0816
J0056	3	0.160	2.9	0.0829
J0061	4	0.060	3.0	0.0855
J0062	4	0.080	3.1	0.0847
J0063	4	0.120	3.0	0.0841
J0064	4	0.160	2.9	0.0863
J0065	4.5	0.058	3.0	0.0878
J0066	4.5	0.080	3.0	0.0870
J0067	4.5	0.120	3.0	0.0858
J0068	4.5	0.180	2.9	0.0874
K0001	3.5	0.060	3.0	0.0802
K0002	3.5	0.080	3.0	0.0808
K0003	3.5	0.094	3.1	0.0800
K0004	3.5	0.140	3.0	0.0810
K0005	4	0.054	3.0	0.0817
K0006	4	0.068	3.0	0.0826
K0007	4	0.094	3.0	0.0823
K0008	4	0.140	2.9	0.0833
K0016	5	0.140	2.9	0.0868
K0017	3.5	0.060	3.0	0.0801
K0018	3.5	0.080	3.1	0.0801
K0019	3.5	0.120	3.1	0.0789
K0020	3.5	0.160	2.9	0.0829
K0021	4	0.060	3.0	0.0821
K0022	4	0.080	3.0	0.0823
K0023	4	0.120	3.0	0.0817
K0024	4	0.160	2.9	0.0856
K0032	5	0.160	2.8	0.0888
L0009	4.5	0.052	3.0	0.0839
L0010	4.5	0.070	3.0	0.0841
L0011	4.5	0.100	3.0	0.0837
L0012	4.5	0.140	2.9	0.0849
L0013	5	0.052	3.0	0.0854
L0014	5	0.080	3.0	0.0856
L0015	5	0.100	3.0	0.0853
L0016	5	0.140	2.9	0.0862
L0017	5.5	0.060	3.0	0.0869
L0018	5.5	0.068	3.0	0.0878
L0019	5.5	0.100	3.0	0.0869
L0020	5.5	0.140	2.9	0.0882
L0029	4.5	0.058	3.0	0.0844
L0030	4.5	0.080	3.1	0.0835
L0031	4.5	0.120	3.1	0.0827
L0032	4.5	0.180	3.0	0.0849
L0033	5	0.060	3.0	0.0852
L0034	5	0.080	3.0	0.0852
L0037	5.5	0.060	3.0	0.0870
L0038	5.5	0.080	3.0	0.0870
L0039	5.5	0.120	2.9	0.0867
L0040	5.5	0.180	2.9	0.0894
P0001	6	0.050	3.0	0.0867
P0004	6	0.140	2.9	0.0867
P0006	8	0.080	2.9	0.0931
P0007	8	0.100	2.9	0.0926
P0008	8	0.140	2.9	0.0910
P0009	10	0.060	2.9	0.0971
P0013	6	0.054	3.0	0.0871
P0014	6	0.078	2.9	0.0885
P0016	6	0.160	2.8	0.0900
P0017	8	0.060	3.0	0.0927
P0020	8	0.180	2.9	0.0936
P0022	10	0.080	2.9	0.0979
P0023	10	0.120	2.9	0.0968
P0024	10	0.180	3.0	0.0957

Equations

$$\Delta E$$

$$\Delta J$$

$$N_{\ell m}(t)$$

$$\Delta E^{\prime }$$

$$\Delta J^{\prime }$$

$$M_f$$

$$M_{\text{irr}}$$

$$q_f=\frac{1}{8\pi }{\oint }_S{K}_{ij}{s}^i{\varphi }^{j}dA,$$

Full text

I Introduction

The gravitational wave (GW) detection of several binary black hole (BBH) mergers Abbott et al. (2016a, b, 2017a, 2017b, 2017c); The LIGO Scientific Collaboration et al. (2018), and the first multi-messenger observation of two colliding neutron stars (NSs) in gravitational and electromagnetic waves Abbott et al. (2017d) have shed light into the nature of gravity in the most extreme astrophysical settings, and has unveiled the identity of the central engines that power the most energetic electromagnetic explosions in the Universe Abbott et al. (2017e); Coulter et al. (2017); The LIGO Scientific Collaboration et al. (2017a); Abbott et al. (2017f), while also providing the means to put at work visionary methods to use GWs to quantify the rate of expansion of the Universe Abbott et al. (2017g); Fishbach et al. (2018); Schutz (1986); Holz and Hughes (2005).

Along this trail of discovery, it has also become evident that numerical relativity (NR) plays a central role to understand the physics of GW sources, and to inform the development of signal-processing algorithms to detect and characterize these astrophysical events Chu et al. (2016); Mroué et al. (2013); Kumar et al. (2016); Abbott et al. (2016c); The LIGO Scientific Collaboration et al. (2017b); Lange et al. (2017), and astrophysical sources that still await discovery George and Huerta (2018, 2018); Shen et al. (2017, 2019); George et al. (2018); Huerta et al. (2018); Rebei et al. (2019); Hinder et al. (2010); 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).

In preparation for the characterization of BBH mergers whose astrophysical properties span a parameter space that has not yet been probed by existing GW detections, several NR groups are working in earnest to construct large-scale NR waveform catalogs Mroué et al. (2013); Healy et al. (2017a); Jani et al. (2016). Since these activities have thus far focused on the study of quasi-circular BBH mergers, in this article we fill in a critical void in the literature by presenting a comprehensive study of the physics of moderately eccentric BBH mergers.

The rationale for this study is multifold. From the perspective of electromagnetic observations, recent findings are consistent with the existence of stellar-mass BHs in the vicinity of the Galactic center, and in the Galactic Cluster M22 Hailey et al. (2018); Sippel and Hurley (2013); Strader et al. (2012); Samsing et al. (2018a). These observations have triggered the development of numerical models that provide a realistic description of the formation and retention of BBHs in dense stellar environments, correcting previous calculations based on $N$-body simulations that did not include post-Newtonian corrections Blanchet (2014) to model the orbital dynamics of these systems, thereby underestimating the merger rates of these systems by orders of magnitude Samsing and Ramirez-Ruiz (2017); Samsing et al. (2014); Samsing (2018); Leigh et al. (2018); Samsing et al. (2018b, a, 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. (2018); Kremer et al. (2018); Lopez et al. (2018); Hoang et al. (2018a); Gondán et al. (2018a); Hoang et al. (2018b); 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. (2018b); Antonini et al. (2018, 2016). In summary, we have evidence for the existence of stellar-mass BHs that may form eccentric compact binary systems in dense stellar environments, and consequently be detected through GW emission. Through this study, we provide new insights into the physics of these astrophysically motivated sources.

Furthermore, as discussed in Klimenko et al. (2016), no matched-filtering algorithm has been presented in the literature that is tailored for the detection of $\ell=\lvert m\rvert=2$ eccentric waveforms Klimenko et al. (2016). However, signal processing algorithms based on deep neural networks have been used to demonstrate that moderately eccentric BBH mergers can be detected and characterized from real LIGO noise, considering both NR waveforms that only include the leading order quadrupole mode $\ell=\lvert m\rvert=2$ George and Huerta (2018, 2018); Shen et al. (2017); George et al. (2018), and higher-order waveform multipoles Rebei et al. (2019). We expect that this NR waveform catalog may be used to quantify the sensitivity of burst searches, and of next-generation neural network models that are tailored to detect and characterize eccentric BBH mergers. In summary, activities around modeling, detection and characterization of eccentric BBHs are reaching the required level of maturity to establish or rule out the existence of compact binary populations in dense stellar environments.

To advance our understanding of the physics of compact binary populations in dense stellar environments, in this article we introduce a NR waveform catalog that describes eccentric BBH mergers, and utilize it to get insights into the dynamics of these GW sources, e.g., the energy and angular momentum loss through GW emission, and the astrophysical properties of the BH remnant, i.e., its final mass and spin as a function of initial eccentricity and mass-ratio, as well as the recoil velocity of the BH remnant. These studies will inform ongoing GW modeling efforts, and the development of signal-processing algorithms to search for and identify these sources. This article is organized as follows. Section II describes the properties of our NR catalog. In Section III we compute the energy and angular momentum radiated away through GW emission, making pair-wise comparisons between NR waveforms that include the $\ell=\lvert m\rvert=2$ mode or higher-order waveform multipoles. In Section IV we compute the astrophysical properties of the BH remnants, and compare these results with those obtained for quasi-circular BBH mergers. We describe the relevance of these analyses in terms of GW modeling efforts for eccentric BBH mergers in Section V. We summarize our findings and outline future directions of work in Section VI.

II Numerical relativity catalog

We have produced a catalog of 89 simulations with the open source, NR software, the Einstein Toolkit etw ; Loffler et al. (2012); Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1998); Baker et al. (2006); Campanelli et al. (2006); Pollney et al. (2011); Wardell et al. (2016); Pollney et al. (2011); Thomas and Schnetter (2010); Löffler et al. (2012); Ansorg et al. (2004); Diener et al. (2007); Dreyer et al. (2003); Schnetter et al. (2004); Thornburg (2004); Brown et al. (2009); Husa et al. (2006); Kranc . This catalog describes non-spinning BBHs with mass-ratios $1\leq q\leq 10$ and eccentricities as high as $e_{0}=0.18$ fifteen cycles before merger. A visualization of this catalog may be found at NCSA Gravity Group (2019a, b). We have post-processed the data products of these simulations using the open source software stack POWER Johnson et al. (2018), and extracted the modes $(\ell,\,\lvert m\rvert)=\{(2,\,2),\,(2,\,1),\,(3,\,3),\,(3,\,2),\,(3,\,1),\,(4,\,4),\,(4,\,3),\,(4,\,2)$, $(4,\,1)\}$. As described in Appendix B, each of these simulations was produced with several levels of resolution to quantify convergence. The real part of the $\ell=\lvert m\rvert=2$ mode, extracted at future null infinity, for each NR waveform is presented in Figure 1. The properties of these NR waveforms are listed in Table LABEL:properties.

Characterizing the properties of the NR waveforms presented in Table LABEL:properties requires the construction of a method to quantify the orbital eccentricity of these simulations. Using the orbital evolution of these simulations to obtain an estimate of the orbital eccentricity is inadequate due to the gauge-dependent nature of the binary’s orbit. On the other hand, methods to construct initial data for spinning BHs on quasi-circular orbits have also introduced definitions of orbital eccentricity, based on orbital separations and waveform phase and amplitude of the Weyl scalar $\psi_{4}$ Healy et al. (2017b). However, while the scope of the method introduced in Healy et al. (2017b) is to construct high-quality initial data for quasi-circular mergers, and therefore, using ${\cal{O}}(e)$ approximations to model the effect of eccentricity may suffice, we aim to measure larger values of orbital eccentricity.

To address this matter, we have used the inspiral-merger-ringdown ENIGMA waveform model introduced in Huerta et al. (2018) to determine the eccentricity, mean anomaly and gauge-invariant frequency parameters, ($e_{0},\,\ell_{0},\,x_{0}$), that optimally describe each NR waveform in our catalog. We do this by finding the ($e_{0},\,\ell_{0},\,x_{0}$) triplet that maximizes the overlap between each NR waveform and its ENIGMA counterpart. In Habib and Huerta (2019) we quantified the optimal time window to remove junk radiation while keeping intact the signatures of eccentricity at early times in the NR waveforms. Such time range is given by $t\leq 60M$. A detailed description of this method, including the corresponding open source software stack for its use to characterize NR waveforms catalogs at scale, is presented in an accompanying article Habib and Huerta (2019).

In brief, we construct our method using the inspiral evolution of the ENIGMA waveform model, which contains state-of-the-art post-Newtonian corrections for eccentric binaries, which include eccentricity corrections in the conservative and radiative pieces up ${\cal{O}}(e^{12})$, including instantaneous, tails and tails-of-tails contributions, and a contribution due to nonlinear memory; and quasi-circular corrections both from post-Newtonian, self-force and perturbative calculations up to ${\cal{O}}(x^{6})$ Huerta et al. (2017, 2018). Furthermore, in Huerta et al. (2018), we have demonstrated that $e_{0}=0$ ENIGMA waveforms capture the dynamics of quasi-circular BBH mergers with excellent accuracy. We showed this by computing overlaps between quasi-circular ENIGMA waveforms and their quasi-circular Effective One Body (EOB) counterparts Bohé et al. (2017). Assuming advanced LIGO’s Zero Detuned High Power sensitivity Barsotti et al. (2018), and using an initial GW frequency of 15 Hz to compute the overlaps, Figure 2 in Huerta et al. (2018) shows that the overlap between quasi-circular ENIGMA and EOB waveforms is ${\cal{O}}\geq 0.99$. Since the waveforms we are characterizing in this study are much shorter than those used to assess the accuracy of the ENIGMA model in the quasi-circular limit, it follows that ENIGMA will capture the dynamics of moderately eccentric systems with excellent accuracy.

It is worth highlighting that while the ENIGMA waveform model was originally validated with eccentric NR waveforms that describe BBH mergers with mass-ratios $q\leq 5.5$ and eccentricities $e_{0}\leq 0.18$ twenty cycles before merger Huerta et al. (2018), it is through this analysis, and with the availability of new NR waveforms, that we can now report that the ENIGMA model can accurately describe BBH mergers with mass-ratios up to $q=10$ with $e_{0}\leq 0.18$ fifteen cycles before merger.

III Energy and angular momentum emission of eccentric black hole mergers

For each NR waveform in our catalog, we have quantified the energy, $E$, and angular momentum, $J$, radiated away through GW emission using the relations Damour et al. (2012)

[TABLE]

where $N^{\ell m}(t)$ represents the complex news function at infinity. The integration is done from the time the NR waveform is free from junk radiation, $t_{0}=60M$, to the final sample time of the NR waveform, $t$. For these calculations we have considered the $(\ell,\,\lvert m\rvert)=\{(2,\,2),\,(2,\,1),\,(3,\,3),\,(3,\,2),\,(3,\,1),\,(4,\,4),\,(4,\,3),\,(4,\,2)$, $(4,\,1)\}$ modes. It is worth pointing out that the choice $t_{0}\rightarrow 60M$ is informed by the study presented in Habib and Huerta (2019), which demonstrated that this choice removes high-frequency noise while keeping intact the imprints of eccentricity in the NR waveforms once they are free from junk radiation. Using Eqs. (1)-(3), in Figure 2 we quantify the importance of including higher-order waveform modes to compute the energy and angular momentum carried away by GWs. We do this through pair-wise comparisons between NR waveforms that include either all the modes listed above, or just the $\ell=\lvert m\rvert=2$ mode, using the relations

[TABLE]

Using the two highest resolution runs for each simulation in our catalog, we computed $(\Delta E^{\prime},\,\Delta J^{\prime})$, and found that the largest difference between these two independent measurements is $\leq 5\%$. The values we present in Figure 2 were extracted from the highest resolution runs. We notice that for each mass-ratio BBH population, i.e., if we consider a given set of markers in the panels of Figure 2, $(\Delta E^{\prime},\,\Delta J^{\prime})$ are nearly constant across the eccentricity range that we have explored in this study. In different words, $(\Delta E^{\prime},\,\Delta J^{\prime})$ are constant polynomials in eccentricity for $e_{0}\leq 0.2$. We observe a minor deviation from this pattern at the high-end of the eccentricity range for the most asymmetric mass-ratio BBH systems. To be precise, if we fit a constant polynomial using the two lowest eccentricity samples for each mass-ratio population, we find that the largest deviation occurs for the most eccentric sample of the $q=10$ BBHs, with a fractional error $\leq 8\%$ for the measurement of $\Delta J^{\prime}$.

These results also show that that for systems with $q\geq 5$ it is essential to include higher-order waveform modes to accurately describe the dynamics of eccentric BBH mergers. This result is consistent with recent studies Rebei et al. (2019), which indicate that the inclusion of higher-order modes for non-spinning, eccentric BBH mergers has a more significant impact for GW detection, in the context of signal-to-noise ratio calculations, than for their non-spinning, quasi-circular BBH counterparts. It is worth highlighting that the eccentric NR waveform we have produced for this analysis for $q\geq 5$ are the first of their kind in the literature, so these results shed new light on the importance of including higher-order waveforms modes for the modeling of radiated energy and angular momentum of eccentric BBH mergers.

IV Final mass, spin, and recoil velocity of post-merger black holes

We have computed the final mass, $M_{f}$, and final spin, $q_{f}$, of the BH remnant using the QuasiLocalMeasures thorn of the Einstein Toolkit. The final mass is given by

[TABLE]

$M_{\text{irr}}$ is the irreducible mass, given in terms of the BHs’ event horizon area, $A$. $q_{f}$ is computed as the Komar angular momentum Misner et al. (1973); Wald (1984); Poisson (2009); Jaramillo and Gourgoulhon (2011)

[TABLE]

where the integral is over the surface, $S$, of the apparent horizon, $K_{ij}$ is the extrinsic curvature, $s^{i}$ is a spacelike, outward normal to the horizon, and $\phi^{i}$ is a Killing vector associated with the rotational symmetry around the spin axis.

In Figure 3 we present results for $(M_{f},\,q_{f})$. As before, we used our two highest resolution runs for each NR simulation to compute these observables, and found that these two independent measurements differ by $\leq 3\%$. The results presented for $(M_{f},\,q_{f})$ in Figure 3 were extracted from the highest resolution runs in our catalog. We notice that for each mass-ratio population, i.e., for a given set of markers, the final mass and spin of the BH remnant are nearly independent of eccentricity in the range $e_{0}\leq 0.2$. We conclude this since both $(M_{f},\,q_{f})$ can be described as constant polynomials in eccentricity within the range we have considered in this study. We can directly compare these results using formulae derived for quasi-circular BBH mergers in Baker et al. (2008); Hofmann et al. (2016). Notice that we have included horizontal gray lines in both panels that provide the predictions for $(M_{f},\,q_{f})$ in the $e_{0}\rightarrow 0$ limit.

For the $M_{f}$ results (top panel in Figure 3) the gray lines present the quasi-circular prediction for the final mass of the BH remnant for the mass-ratios $q=\{1,\,1.5,\,2,\,2.5,\,3,\,3.5,\,4,\,4.5,\,5,\,6,\,8,\,10\}$ from bottom to top, respectively. We notice that the equal mass-eccentric BBH population presents the largest deviation from the quasi-circular prediction. However, this discrepancy is $\leq 1\%$. In the case of the final spin of the BH remnant, the bottom panel of Figure 3 also presents the quasi-circular predictions for this observable. Notice, however, that in this case, the grey lines describe the mass-ratios listed above but now from top to bottom. As in the case of $M_{f}$, our results for the final spin of moderately eccentric BBH mergers are fairly consistent with results obtained from quasi-circular BBH mergers. This can only the case if the eccentric NR waveforms we have produced in this catalog circularize prior to merger. We have explored this scenario in detail, and have found that this is indeed the case. For a sample case, Figure 4 presents two waveform signals produced by BHs that have the same separation, but different initial eccentricity. The eccentric waveform contains all the telltale signatures of eccentricity, i.e., significant modulations in the amplitude and phase at early times, which correspond to periapse (local maxima) and apoapse (local minima) passages. We also observe that the waveform circularizes very rapidly, from $e_{0}=0.18$ fifteen cycles before merger, turning into a quasi-circular waveform signal near the merger event. This is the reason why the results presented in Figure 3 are consistent with their quasi-circular counterparts.

Earlier work on this front includes Hinder et al. (2008), which presented calculations for the final spin and circularization of equal-mass eccentric BBH mergers, and showed that for BBH mergers with larger initial eccentricities than those considered in this work, the final spin of the BH remnant is greater than its quasi-circular counterpart. Additionally, reference Hinder et al. (2018) discussed the circularization of moderately eccentric BBH mergers with $q\leq 3$. In this article, we provide a systematic study of the observables $(q_{f},\,M_{f})$ to furnish evidence for the circularization of eccentric BBH mergers with $q\leq 10$ and e_{0}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}0.18 fifteen cycles before merger.

We have also computed the recoil velocity of eccentric BBH mergers, $\lvert v\rvert_{\textrm{kick}}$, using the PunctureTracker thorn in the Einstein Toolkit. To do this, we have considered the last 100M of evolution of our NR simulations, and a simple first order finite difference formula for the velocity in terms on the measured locations. In Table 1 we present the minima and maxima of the recoil velocity, $\lvert v\rvert_{\textrm{kick}}$, for the range of eccentricities we consider for each mass-ratio. We have also obtained gauge-invariant perturbations following Pollney et al. (2007) (see Eqs. (33)–(39) therein). These gauge-invariant results are presented in the last column of Table 1.

Important observations to be drawn from Table 1 include: (i) the kick velocity of BH remnants produced by quasi-circular BH mergers, $\lvert v\rvert^{e_{0}\rightarrow 0}_{\textrm{kick}}$, was obtained using the formulae presented in Lousto and Zlochower (2013); (ii) the recoil velocity for all the $q=1$ eccentric BBH mergers in our catalog is $\lvert v\rvert_{\textrm{kick}}=0$, which is consistent with results obtained for non-spinning, quasi-circular BBH mergers Koppitz et al. (2007); (iii) the kick velocities for our $1\leq q\leq 10$ population of eccentric BBH mergers are fairly consistent with the expected values of their quasi-circular counterparts, even though the formulae used to estimate $\lvert v\rvert^{e_{0}\rightarrow 0}_{\textrm{kick}}$ was calibrated with quasi-circular BBH mergers with mass-ratios $q\leq 8$.

These results for $(M_{f},\,q_{f},\,\lvert v\rvert_{\textrm{kick}})$ cover an entirely new region of parameter space in the modeling of eccentric BBH mergers, providing new insights into the physics of these GW sources. We discuss the implications of these findings in the following section.

V Implications for the modeling and detection of eccentric mergers

To date, there are only a handful of inspiral-merger-ringdown waveform models that describe the GW emission of eccentric BBH mergers Huerta et al. (2018, 2017); Hinder et al. (2018); Hinderer and Babak (2017); Cao and Han (2017). These models assume that moderately eccentric BBHs circularize prior to the merger event. This assumption is sound, in light of the results presented in the previous section, for BBHs with $q\leq 10$ and whose residual eccentricity is as high as $e_{0}\leq 0.18$ just fifteen cycles before merger.

Furthermore, we have found that for the most extreme sample of our NR catalog, e.g., P0024, which represents BBHs with $q=10$ and $e_{0}=0.18$ fifteen cycles before merger, circularization is only attained right before merger, as shown in Figure 4. In different words, while assuming circularization of moderately eccentric BBH mergers is a reasonable ansatz, this also means that the modeling of these GW sources demands the development of an inspiral evolution scheme that provides an accurate description of the dynamical evolution of these objects throughout the inspiral evolution, and which remains accurate one or two cycles before merger. To accomplish this level of accuracy so late in the inspiral evolution, we showed in Huerta et al. (2018, 2017) that the inspiral evolution should include, at the very least, higher-order eccentric post-Newtonian corrections for the instantaneous and tails and tails-of-tails pieces, as well as contributions due to non-linear memory, and higher-order self-force and BH perturbation theory corrections.

Future source modeling efforts to describe the inspiral evolution of spinning BBHs on eccentric orbits should include new developments from post-Newtonian, self-force and perturbation theory formalisms Bini et al. (2016a, b); Bini and Geralico (2018); Bini et al. (2018); Kavanagh et al. (2017); Bini et al. (2016c); Hinderer and Babak (2017); Akcay and van de Meent (2016); Akcay et al. (2017); Fujita et al. (2017); Le Tiec (2015); Osburn et al. (2014); Warburton et al. (2017); van de Meent (2016, 2017, 2018); van de Meent and Warburton (2018); Ireland et al. (2019). These schemes may be complemented with stand-alone merger models designed with machine learning, or by directly attaching merger waveforms from NR surrogate waveform families Blackman et al. (2015); Varma et al. (2018); Doctor et al. (2017); Moore et al. (2016); Huerta et al. (2018); Moore et al. (2016); Varma et al. (2019). The validation of these models with eccentric NR simulations will be essential to assess their accuracy and reliability for the detection and characterization of compact binary populations in dense stellar environments.

This waveform catalog may also be used to assess the sensitivity of burst searches to detect eccentric BBH mergers Klimenko et al. (2004); Klimenko and Mitselmakher (2004); Klimenko et al. (2008, 2016); Tiwari et al. (2016), and to train neural network models to detect and characterize these GW sources Rebei et al. (2019); George and Huerta (2018, 2018). These studies will be pursued using this NR waveform catalog.

VI Conclusion

We have studied the physics of eccentric BBH mergers using a NR waveform catalog that describes BBH systems with mass-ratios $q\leq 10$, and initial eccentricities $e\leq 0.18$ up to fifteen cycles before merger.

We quantified the importance of including higher-order waveform modes to compute the energy and angular momentum carried away by GWs in eccentric BBH mergers. We have also demonstrated that the properties of BH remnants described by our NR catalog are consistent with their quasi-circular counterparts, which provides evidence for the circularization of moderately eccentric BBH mergers. We have also computed recoil velocities of BH remnants produced by eccentric BBH mergers, and found that these are fairly similar to those computed in the literature for non-spinning, quasi-circular BBH mergers.

Based on these analyses, we have provided evidence that existing source modeling efforts that assume the circularization of moderately eccentric BBH mergers is sound. We have also shown that since circularization takes place close to the merger event, any semi-analytical model that is used to describe these GW sources should include higher-order corrections both to the conservative and radiative pieces of the source’s dynamics, and to the waveform strain.

Recent studies in the literature have identified parameter space degeneracies between orbital eccentricity and spin corrections Huerta et al. (2017). In order to get better insights into this finding, it is essential to understand the dynamics of spinning BHs on eccentric orbits, and then use the NR waveform catalog we have introduced in this study to carefully assess in which regions of parameter space such a degeneracy may be broken to distinguish these two compact binary populations.

The construction of a NR waveform catalog for spinning BHs on eccentric orbits is already underway to shed light on this timely, and astrophysically motivated study. Specific aspects to address in such study will encompass: (i) orbital configurations that significantly shorten the length of waveform signals, i.e., eccentricity and spin anti-aligned configurations; (ii) competing effects to determine the length of waveform signals, i.e., rapidly spinning BHs on spin-aligned configurations (which increase the length of waveforms as compared to non-spinning BBHs) vs moderate values of initial eccentricity (which decrease the length of waveforms as compared to quasi-circular BBHs); and (iii) identify telltale signatures of GW sources that can be used to infer the existence of eccentric compact binary populations through GW observations, e.g., astrophysical properties of the BH remnant, and the coupling of eccentricity and spin-spin and spin-orbit effects at peri-apse passages during the inspiral evolution of these systems.

The modes we have extracted in this study with the open source POWER Johnson et al. (2018) package do not include $m=0$ memory modes. Their extraction requires the use of the Cauchy Characteristic Extraction method Reisswig et al. (2010). Given the importance of these modes for the characterization of eccentric BBH mergers, we will present these modes in a forthcoming study, accompanied by a systematic analysis on the observability of these modes with second and third generation GW detectors.

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 and TG-PHY160053 grants are gratefully acknowledged.

Appendix A Properties of Numerical Relativity Catalog

Table LABEL:properties lists the properties of our numerical relativity catalog.

Appendix B Convergence of the numerical waveforms

We use a grid setup based on the setup used in Wardell et al. (2016). There is a central, mesh refined cubical region of the grid in which Cartesian coordinates are used, surround by 6 regions that make up a cubed sphere grid with constant angular resolution.

We use 8th order finite differencing operators to compute spatial derivatives of the spacetime quantities in the Einstein field equations. This requires the use of $5$ ghost zones, and together with using a classical 44th order Runge-Kutta timestepper implies that each refined region is surrounded by $20$ points that are filled in via prolongation from the next coarser region. We use vertex centered 5th order prolongation operators rather than full 8th order prolongation operators.

The cubical region employs mesh refinement with the resolution on the coarsest grid being $h_{\text{coarse}}=1.92\,M$. Each of the black holes is surrounded by a set of nested moving boxes such that the resolution in the finest box containing the black hole $i$ is $1.2M_{i}/(N_{l}-1)$ where $M_{i}$ is the initial mass parameter of black hole $i$ and $N_{l}$ is the number of points used for the resolution level $l$ simulation. In our simulations we used $N_{l}=32,36,40,44$, where $N_{l}=44$ was only used for simulations with a mass ratio $q>5$. The finest box surrounding each black hole has a radius of $1.2\,M_{i}$ and each coarser box has twice the radius of the next finer one. During the simulation we track the location of each black hole and keep the set of nested refined boxes approximately centered on the black hole. Finally the outer edge of the cubical region is chosen large enough to contain all refined regions including their prolongation regions.

In the spherical region we choose an angular resolution of $h_{\text{angular}}=\pi/\left(4\,N_{l}\right)$ and a radial resolution of $1.92\,M$ which matches the coarsest resolution in the Cartesian grid. The outer boundary is chosen such that it is causally disconnected from the outermost detector at which we extract gravitational waves from.

We use a time step $\Delta t=0.864\,M$ on the coarsest level, corresponding to a Courant-Friedrichs-Lewy condition of $\Delta t/h_{\text{coarse}}=0.45$ which is held constant on the finer levels by decreasing their time step size.

We extract gravitational waves using modes of the Weyl scalar $\psi_{4}$ extracted on coordinate spheres of radius $r_{\text{det},i}=100\,M,115\,M,136\,M,167\,M,214\,M,300\,M,$ $500\,M$.

Using 8th order finite differencing operators our simulations would, under ideal circumstances, converge towards the correct solution with an error term which scales like $h^{8}$, where $h$ is the spatial resolution of the simulation. However due to lower order schemes present in the simulation, for example the interpolation at mesh refinement boundaries which is only 5th order accurate, as well as artifacts caused by the adaptive mesh refinement logic making independent decisions where to refine for each simulation, the observed convergence order typically differs from 8.

To estimate the convergence of each waveform we simulated each set of physical parameters using at least 3 (4) simulations using increasing resolution for waveforms of mass ratio $q\leq 5$ ($q>5$). We then compute the gravitational wave phase $\phi^{(2,2)}$ from the complex $\ell=m=2$ mode of the spherical harmonic decomposition of the outgoing component of the Weyl scalar $\psi_{4}$ and studied its convergence properties.

Figure 5 shows the rescaled phase differences $\phi^{(2,2)}(h_{n})-\phi^{(2,2)}(h_{\text{high}})$ between the gravitational wave phase obtained from the simulation with resolution $h_{n}$ and the highest resolved simulation. Phase differences have been rescaled such that for a convergent simulation the plotted curves overlap. For case J0046 we observe an approximate convergence order of $N\approx 4.5$ which is within the range of expected values.

Not all simulated cases show clean convergence behaviour, with the convergence order for some of them being larger than 8, which may indicate that our lowest resolution simulation does not adequately resolve the features present in the simulation domain, and others swapping the ordering of phases between the low, medium and high resolution simulations, making an estimate of the convergence order impossible.

Given that large number of simulations, this is to be expected and does not necessarily indicate that the obtained results are incorrect but instead demonstrates the difficulty in controlling the various effects that influence the numerically obtained waveform. Since there are multiple sources of numerical error, and we have chosen parameters such that none is dominant so as to make best use of available computing resources without over-resolving a particular feature, different sources of numerical error potentially cancel each other out, giving rise to unrealistically large (or small) convergence orders.