Eccentricity distributions of eccentric binary black holes in galactic nuclei

TL;DR

This paper models the eccentricity distributions of binary black holes formed in galactic nuclei, showing that gravitational capture produces more eccentric binaries detectable by LIGO than the Kozai-Lidov mechanism, and discusses how future observations can distinguish these channels.

Contribution

It provides the first detailed eccentricity distribution predictions for binary black holes from two main formation channels in galactic nuclei, aiding future gravitational wave data interpretation.

Findings

Approximately 10% of binaries from the Kozai-Lidov mechanism have eccentricities > 0.1 at 10 Hz.

About 75% of binaries from gravitational capture have eccentricities > 0.1 at 10 Hz.

Few tens of detections can constrain formation channel proportions within 0.2 confidence interval.

Abstract

Galactic nuclei are expected to be one of the main sites for formations of eccentric binary black holes (EBBHs), with an estimated detection rate of $\mathcal{O}(1-100$ yr$^{-1})$ with Advanced LIGO (aLIGO) detectors operating at design sensitivity. The two main formation channels of these binaries are gravitational capture and the secular Kozai-Lidov mechanism, with expectedly commensurable formation rates. We used Monte Carlo simulations to construct the eccentricity distributions of EBBHs formed through these channels in galactic nuclei, at the time their gravitational-wave signals enter the aLIGO band at $10$ Hz. We have found that the proportion of binary black holes entering the aLIGO band with eccentricities larger than $0.1$ is $\sim 10$ percent for the secular Kozai-Lidov mechanism, and $\sim 75$ percent for gravitational capture. We show that if future EBBH detection rates with aLIGO will be dominated by EBBHs formed in galactic nuclei, then the proportions of EBBHs formed through the two main channels can be constrained to a $\Delta \mathcal{F}= 0.2$ wide one-sigma confidence interval with a few tens of observations, even if parameter estimation errors are taken into account at realistic levels.

Tables

Table 1. Table 1: Fraction of EBBHs having e LIGO > 0.1 subscript 𝑒 LIGO 0.1 e_{\mathrm{LIGO}}>0.1 at f GW = 10 subscript 𝑓 GW 10 f_{\mathrm{GW}}=10 Hz ( F e LIGO > 0.1 subscript 𝐹 subscript 𝑒 LIGO 0.1 F_{e_{\mathrm{LIGO}}>0.1} ), of EBBHs starting their GW-dominated evolution within the aLIGO band ( F LIGO subscript 𝐹 LIGO F_{\mathrm{LIGO}} ), and of EBBHs having eccentricities larger than 0.1 0.1 0.1 at their last stable orbit ( F e LSO > 0.1 subscript 𝐹 subscript 𝑒 LSO 0.1 F_{e_{\mathrm{LSO}}>0.1} ). The table rows correspond to the gravitational capture (GC) and the secular Kozai-Lidov (KL) formation channels, respectively.

Formation	$F_{e_{\mathrm{LIGO}}>0.1}$	$F_{\mathrm{LIGO}}$	$F_{e_{\mathrm{LSO}}>0.1}$
GC	$75\%$	$28\%$	$11\%$
KL	$11\%$

Equations

$$\delta E=-\frac{85\pi {\eta }^2{M}^{9/2}}{12\sqrt{2}{r}_{p0}^{7/2}},$$

$$\delta L=-\frac{6\pi {\eta }^2{M}^4}{r_{p0}^2},$$

$$r_{p0}={\left(\sqrt{\frac{1}{b^2}+\frac{M^2}{b^4{w}^4}}+\frac{M}{b^2{w}^2}\right)}^{-1},$$

$$a_0=-\frac{\eta M^2}{2E_{\mathrm{fin}}},$$

$$e_0=\sqrt{1+2\frac{E_{\mathrm{fin}}{b}^2{w}^2}{M^3\eta }}.$$

$$b_{\min }\equiv \frac{4M}{w}<b<{\left(\frac{340\pi \eta }{3}\right)}^{1/7}\frac{M}{w^{9/7}}\equiv b_{\max }.$$

$$\frac{\mathrm{d}e}{\mathrm{d}t}=-\frac{304\eta e(1-e{)}^{3/2}}{15M{\rho }_{\mathrm{p}}^4(1+e{)}^{5/2}}\left(1+\frac{121}{304}{e}^2\right).$$

$$\frac{\mathrm{d}{\rho }_{\mathrm{p}}}{\mathrm{d}t}=-\frac{64\eta e(1-e{)}^{1/2}}{5M{\rho }_{\mathrm{p}}^3(1+e{)}^{7/2}}\left(1+\frac{73}{24}{e}^2+\frac{37}{96}{e}^4\right)+$$

$$+\frac{304\eta e(1-e{)}^{3/2}}{15M{\rho }_{\mathrm{p}}^3(1+e{)}^{5/2}}\left(1+\frac{121}{304}{e}^2\right).$$

$${\rho }_{\mathrm{p}}=\frac{c_0{e}^{12/19}}{(1+e)M}{\left(1+\frac{121}{304}{e}^2\right)}^{\frac{870}{2299}},$$

$${\rho }_{\mathrm{p}}(e_{\mathrm{LSO}})=\frac{6+2e_{\mathrm{LSO}}}{1+e_{\mathrm{LSO}}},$$

$$n_{\mathrm{peak}}\approx \frac{2(1+e{)}^{1.1954}}{(1-e^2{)}^{3/2}}.$$

$$f_{\mathrm{GW}}=n_{\mathrm{peak}}\frac{M^{1/2}}{2\pi a^{3/2}},$$

$$t_{\mathrm{merge}}(m_1,m_2,r_{p0},e_0)=\frac{15c_0^4}{304M^3\eta }\mathcal{I}(e_{\mathrm{LSO}},e_0),$$

$$\mathcal{I}(e_{\mathrm{LSO}},e_0)=\underset{e_{\mathrm{LSO}}}{\overset{e_0}{\int }}\frac{e^{29/19}(1+\frac{121}{304}{e}^2{)}^{\frac{1181}{2299}}}{(1-e^2{)}^{3/2}}\mathrm{d}e.$$

$$t_{\mathrm{enc}}\approx \frac{w^3}{12\pi M^2{n}_{\mathrm{tot}}},$$

$$\mathcal{H}=\frac{m_1{m}_2}{2a_1}+\frac{(m_1+m_2)m_3}{2a_2}+$$

$$+\frac{1}{a_2}\sum _{j=2}^{\infty }{\alpha }^j{M}_j{\left(\frac{r_1}{a_1}\right)}^j{\left(\frac{a_2}{r_2}\right)}^{j+1}{P}_j(\cos \Phi ),$$

$$M_j=m_1{m}_2{m}_3\frac{m_1^{j-1}+(-1{)}^j(m_2{)}^{j-1}}{(m_1+m_2{)}^j}.$$

$$r_{\max }=\frac{M_{\mathrm{SMBH}}}{{\sigma }_{\ast }^2},$$

$$M_{\mathrm{SMBH}}\simeq 1.3\times 10^8{M}_{\odot }{\left(\frac{{\sigma }_{\ast }}{200km/s}\right)}^4.$$

$$\frac{\mathrm{d}{N}_{\mathrm{SMBH}}}{\mathrm{d}{M}_{\mathrm{SMBH}}}\propto {\left(\frac{M_{\mathrm{SMBH}}}{M_{\bullet }}\right)}^{-1.25}{e}^{-M_{\mathrm{SMBH}}/M_{\bullet }},$$

$$\frac{\mathrm{d}{N}_{\mathrm{BH}}}{\mathrm{d}{m}_{\mathrm{BH}}}=\frac{(1-\beta )m_{\mathrm{BH}}^{-\beta }}{m_{\mathrm{BH},\max }^{1-\beta }-m_{\mathrm{BH},\min }^{1-\beta }},$$

$$N_{\mathrm{BH}}\approx 20,000\times \frac{M_{\mathrm{SMBH}}}{M_{{\mathrm{SgrA}}^{\ast }}},$$

$$f_m(\mathbf{r},\mathbf{v})=C_{m}E(r,v{)}^{p_m},$$

$$E(r,v)=\frac{M_{\mathrm{SMBH}}}{r}-\frac{v^2}{2}$$

$$p_m=p_0\frac{m}{m_{\mathrm{BH},\max }},$$

$$n_m(r)=\int f_m(\mathbf{r},\mathbf{v}){\mathrm{d}}^{3}v=n_m(r_{\max }){\left(\frac{r}{r_{\max }}\right)}^{-{\alpha }_m},$$

$${\alpha }_m=\frac{3}{2}+p_m.$$

$$n_{\mathrm{MS}}(r_{\max })\simeq 1.38\times 10^5{\mathrm{pc}}^{-3}\sqrt{\frac{10^6{M}_{\odot }}{M_{\mathrm{SMBH}}}}.$$

$$\phi (\bar{v},m)=C_{\bar{v}}(1-{\bar{v}}^2{)}^{p_m},$$

Full text

††pubyear: 2018††pagerange: Eccentricity distributions of eccentric binary black holes in galactic nuclei–Eccentricity distributions of eccentric binary black holes in galactic nuclei

1 Introduction

Second generation gravitational-wave (GW) detectors are state-of-the-art, ground-based, L-shaped interferometers with kilometer-scale arms. The network of Advanced LIGO (aLIGO) GW detectors (Aasi et al., 2015b), consisting of aLIGO-Hanford and aLIGO-Livingston, finished its first observing run (O1) in January 2016. During O1, this network achieved the first detections of GWs from coalescences of binary black holes (BHs) (Abbott et al., 2016). The two aLIGO detectors carried out their second observing run between November 2016 and August 2017, with an improved network sensitivity compared to O1, achieving further observations of binary BHs (Abbott et al., 2017b, d) and the first observation of a binary neutron star (Abbott et al., 2017c). Second generation GW detectors also include Advanced Virgo (Acernese et al., 2015), which joined the aLIGO observing run in August 2017. Two additional GW detectors are planned to join aLIGO and Advanced Virgo: the Japanese KAGRA is under construction and is planned to begin its first observing run in 2020, while LIGO-India is expected to become operational in 2024 at the earliest (Abbott et al., 2018b).

The most common sources for GW detections with second generation detectors are binaries of stellar-mass BHs (Abadie et al., 2010). These systems can form through binary stellar evolution (e.g. Postnov & Yungelson, 2014), as well as through dynamical processes (e.g. Park et al., 2017). Dynamical formation involves non-dissipative processes, such as strong gravitational encounters and dynamical three-body formation, but it also includes gravitational capture, as a dissipative process. Binary BHs formed dynamically in dense stellar environments generally exhibit higher eccentricities compared to those forming in isolated environments in galactic fields (see e.g. Breivik et al., 2016). Since GW emission tends to circularize the orbit of binary BHs, most of the binaries are expected to possess negligible eccentricities (usually $e\lesssim 10^{-3}$) at the time their GW frequency reaches $10$ Hz, i.e. the nominal low-frequency limit of aLIGO detectors (Breivik et al., 2016; Rodriguez et al., 2016). Several studies have shown that signal search methods using circular binary templates (see e.g. Capano et al., 2016, and references therein) are unsuitable for finding binary BHs with $e>0.1$ when entering the aLIGO band (Brown & Zimmerman, 2010; Huerta & Brown, 2013). Some processes, including the ones we studied here in details, are predicted to produce such $e>0.1$ binaries. Hence within the scope of our work, we define eccentric binary BHs (EBBHs) as binary BHs having eccentricities $e>0.1$ at the time their GW frequency reaches $10$ Hz.

According to theoretical models, we expect most EBBHs to form in dense stellar clusters, such as galactic nuclei (GNs) and globular clusters (GCs) (Portegies Zwart & McMillan, 2000; Kocsis et al., 2006; O’Leary et al., 2009). Two of the main formation channels of these binaries are the following: (i) Initially unbound BHs can form binaries during close encounters, if they lose enough kinetic energy to become bound. In almost all cases, this energy loss occurs through GW emission, possibly resulting with an EBBH (e.g. Peters & Mathews, 1963; O’Leary et al., 2009). (ii) EBBHs can also form in hierarchical triples through the secular Kozai-Lidov mechanism (e.g. Wen, 2003; Antonini & Perets, 2012; Antonini et al., 2014; Antognini & Thompson, 2016; Hoang et al., 2018). During this process, binary BHs are expected to reach very high eccentricities (with $1-e$ ranging from $10^{-2}$ to $10^{-5}$, similarly to EBBHs formed through gravitational capture, see Gondán et al., 2018b) with a merger rate comparable with that of EBBHs formed through gravitational capture.

In GNs the expected merger rate density of EBBHs is $\sim 0.8$ Gpc*-3yr-1* for EBBHs formed through gravitational capture (O’Leary et al., 2009). Hoang et al. (2018) showed that the total merger rate density of all binary BHs formed through the secular Kozai-Lidov mechanism in GNs is $\sim 1-3$ Gpc*-3yr-1*, which, based on considerations derived from our results described in Section 4, corresponds to a merger rate density of $\mathcal{O}(0.1$ Gpc*-3yr-1*) for EBBHs. Samsing et al. (2018) and Rodriguez et al. (2018) recently showed that the formation rate of EBBHs produced by binary-single interactions in GCs may be comparable to the formerly mentioned formation rates with an expected merger rate density of $0.05-0.5$ Gpc*-3yr-1*. Another alternative for EBBH formation in GNs and GCs is through multi-body interactions (e.g. Antonini & Rasio, 2016; Zevin et al., 2019), although the merger rate from this channel is expected to be at least an order of magnitude below the rates for the previously mentioned channels.

Binary BHs formed through gravitational capture are expected to be detectable up to a few Gpc by aLIGO detectors operating at design sensitivity, with estimated total detection rates of $\mathcal{O}(1-100$ yr*-1*) for GNs and $\mathcal{O}$(1 yr*-1*) for GCs (see e.g. Kocsis et al., 2006; O’Leary et al., 2009). Although the exact ratio of detection rates for the different EBBH formation channels is yet unknown, we also expect a commensurable detection rate for binary BHs formed through the secular Kozai-Lidov mechanism (see e.g. Antonini et al., 2014; Hoang et al., 2018). Note that these detection rates are non-negligible compared to the rates predicted for circular binary BHs. In fact, the reconstructed spin configuration of one of the observed binary BHs (GW170104, see Abbott et al., 2017a) suggests that it may have formed dynamically in a dense stellar cluster. The relatively high masses of the observed binary BHs are also consistent with the assumption that these binaries were formed through dynamical processes in such environments (see e.g. McKernan et al., 2018).

According to the leading order results in Peters & Mathews (1963), compared to circular binaries, EBBHs are detectable from greater distances and in a wider mass range due to the strength and spectral richness of their GW signals (O’Leary et al., 2009; Kocsis & Levin, 2012). By detecting these signals, in principle, many parameters of their sources can be reconstructed, including their orbital parameters (eccentricity and pericenter distance) at the time the frequency of the binary’s GW signal reaches $10$ Hz. Numerical tools for simulating GW signals of EBBHs already exist (see e.g. Levin et al., 2011; Csizmadia et al., 2012; East et al., 2013) and are available for testing search algorithms. There are also search methods available aiming to find GW signals of highly eccentric EBBHs (Tai et al., 2014; Tiwari et al., 2016) and EBBHs with $e<0.6$ (Coughlin et al., 2015). For reconstructing circular binary parameters, several algorithms are used, e.g. BAYESTAR (Singer & Price, 2016), LALInference (Veitch et al., 2015), and GstLAL (Cannon et al., 2012; Privitera et al., 2014). Coherent WaveBurst (Klimenko et al., 2016) and BayesWave (Cornish & Littenberg, 2015; Bécsy et al., 2017; Millhouse et al., 2018) algorithms offer ways for model-independent parameter reconstruction. The development of algorithms aiming to recover EBBH parameters is still underway.

In the work presented in this paper we used Monte Carlo simulations to construct the eccentricity distributions of EBBHs formed through gravitational capture and the secular Kozai-Lidov mechanism in GNs, at the time their GW frequency enters the aLIGO band. We applied a set of numerical simulations introduced in Gondán et al. (2018b) to determine the expected parameter distributions of EBBHs formed through gravitational capture, however we extended them by introducing randomized masses for the central supermassive BHs (SMBHs). For EBBHs formed through the secular Kozai-Lidov mechanism, we solved the octupole-order equations of motion of Naoz et al. (2013a) in which we included gravitational radiation, while making reasonable assumptions on the initial parameter distributions of triple systems. In this paper we also show that if future EBBH detection rates with aLIGO will be dominated by EBBHs formed in GNs, then it is feasible to constrain the proportions of EBBHs formed through the two main channels (i.e. the branching ratios). We performed a statistical test to determine the minimum number of EBBH observations needed to constrain the branching ratio of gravitational capture (from now on, denoted by $\mathcal{F}\mathrm{{}_{GC}}$) to an arbitrarily chosen $\Delta\mathcal{F}\mathrm{{}_{GC}}=0.2$ wide one-sigma confidence interval (note that for symmetrical errors $\Delta\mathcal{F}\mathrm{{}_{GC}}=0.2$ corresponds to a $\pm 10$ percent range of the total $[0,1]$ interval). First, we carried out our tests with the idealized case of having no errors on the estimated parameters. We then repeated our analysis, assuming statistical errors on the estimated parameters. Since EBBH parameter estimation methods are still in development, we only considered a set of plausible errors consistent with the results of previous studies (see Gondán et al., 2018a; Lower et al., 2018; Gondán & Kocsis, 2019).

This paper is organized as follows. In Section 2 we review the mechanism of EBBH formation through gravitational capture, as well as through the secular Kozai-Lidov mechanism, and describe the evolution of EBBH orbits using the leading order results of Peters (1964). In Section 3 we introduce the model we used for GNs in our work and describe the numerical methods we applied in our Monte Carlo simulations to generate EBBHs. We present our results in Section 4. Finally, we offer our conclusions and summarize the implications of our work in Section 5.

2 Eccentric Binary Black Holes

In this section we summarize the formational mechanism of EBBHs through gravitational capture (see Section 2.1) and through the secular Kozai-Lidov mechanism (see Section 2.3). First, we introduce the equations to calculate the initial orbital parameters for EBBHs formed through gravitational capture, as well as their evolution during the eccentric inspiraling of the two BHs. We then consider the possible disruptions of formed binaries by gravitational interactions with third objects. Finally, we describe the Kozai-Lidov mechanism as an alternative formational process of EBBHs, and show the importance of using equations of motion beyond the quadrupole approximation.

In our calculations, we use the geometric unit system ($G=c=1$). We denote the total mass of the two BHs forming a binary by $M=m_{1}+m_{2}$, and their symmetrical mass ratio by $\eta=m_{1}m_{2}/M^{2}$, where $m_{1}$ and $m_{2}$ are the masses of the BHs, defined as $m_{1}\geq m_{2}$. The reduced mass of the binary is $\mu=M\eta$. We describe the orbital evolution in a quasiperiodic manner. We define the dimensionless pericenter distance as $\rho_{\mathrm{p}}=r_{\mathrm{p}}/M$, where $r_{\mathrm{p}}$ is the pericenter distance. The mean orbital frequency from Kepler’s law is $f_{\mathrm{orb}}=(2\pi M)^{-1}\rho_{\mathrm{p}}^{-3/2}(1-e)^{3/2}$, where $e$ is the eccentricity of the orbit. Note that in this quasiperiodic treatment, the time evolution of the system is incorporated in the time-dependence of the orbital parameters $e$ and $\rho_{\mathrm{p}}$. The following discussion is based on O’Leary et al. (2009) and Gondán et al. (2018b).

2.1 Formation through gravitational capture

During a close encounter, two BHs can form a binary if they emit enough energy in the form of GWs to become bound. Due to the relativistic nature of such events and the fact that the velocity dispersion in GNs is $\sigma\ll 1$, the encounters are always nearly parabolic (Quinlan & Shapiro, 1987; Lee, 1993), and the binaries formed this way are always highly eccentric ($1-e_{0}^{2}\sim 10^{-1}-10^{-5}$, where $e_{0}$ is the initial eccentricity of the formed binary, see e.g. Gondán et al., 2018b). The change in energy and angular momentum due to GW emission during the encounter is

{ceqn}

[TABLE]

{ceqn}

[TABLE]

where

{ceqn}

[TABLE]

we denoted the initial values of orbital parameters with a ’0’ symbol in the lower index, $b$ is the impact parameter of the encounter, and $w$ is the relative velocity of the two BHs (see Peters & Mathews, 1963; Turner, 1977).

The initial orbital parameters of the binary are determined by the final energy ($E_{\mathrm{fin}}$) and final angular momentum ($L_{\mathrm{fin}}$) of the system after the first encounter. O’Leary et al. (2009) showed that $\delta L$ (see Eq. 2) is negligible for nearly all first encounters, therefore, for simplicity, we set $\delta L=0$. In case the BHs emit enough energy to have $E_{\mathrm{fin}}<0$, the system becomes bound after the first encounter, with an initial semi-major axis

{ceqn}

[TABLE]

and initial eccentricity

{ceqn}

[TABLE]

Note that there is an upper and lower limit for the $b$ impact parameter of encounters that produce EBBHs for fixed values of $w$. The upper limit is set by the condition $E_{\mathrm{fin}}<0$, since a distant encounter of BHs does not produce bound binary system. Additionally, there is a lower limit set by the fact that we require BHs to avoid direct collision during their first encounter (Kocsis et al., 2006; O’Leary et al., 2009). Therefore, in leading order, $b$ must satisfy,

{ceqn}

[TABLE]

The leading order approximation in the allowed range of impact parameters is in excellent agreement with $2.5$ and $3.5$ order post-Newtonian simulations (Kocsis & Levin, 2012). Note that in realistic GNs, relativistic corrections only modify the capture cross section by less than $10$ percent (Kocsis et al., 2006; O’Leary et al., 2009).

2.2 Evolution of isolated eccentric binaries

The time derivative of the eccentricity in the leading order approximation is given by the equation (Peters, 1964)

{ceqn}

[TABLE]

By utilizing the formula for d$a/$d$t$ (see Eq. 5.6 in Peters, 1964), the time derivative of $\rho_{\mathrm{p}}$ can be expressed as

{ceqn}

[TABLE]

Dividing Eq. (8) by Eq. (7) yields a separable differential equation for $\rho_{\mathrm{p}}(e)$, which gives the solution (Peters, 1964):

{ceqn}

[TABLE]

where $c_{0}$ is a function of the initial parameters that can be determined by setting $e=e_{0}$ and $\rho_{\mathrm{p}}=\rho_{\mathrm{p0}}$ in Eq. (9).

The orbital equations Eq. (7-9) are only valid for $\rho_{\mathrm{p}}\gg 2$, when the BHs are relatively far from each other. At a certain point the binary reaches its last stable orbit (LSO), after which the evolution is no longer quasi-periodic, and the BHs in the binary quickly merge. In the leading order approximation for zero spins and a finite mass ratio, $e_{\mathrm{LSO}}$ can be given by solving the following equation numerically (Cutler et al., 1994):

{ceqn}

[TABLE]

where $\rho_{\mathrm{p}}(e_{\mathrm{LSO}})$ is obtained by substituting $e_{\mathrm{LSO}}$ to Eq. (9). Unlike GW signals of circular binaries, the ones from EBBHs cover a broad band of frequencies, where the maximum power occurs at a harmonic $n_{m}\gg 1$ for $e\approx 1$ ($f_{m}=n_{m}f_{\mathrm{orb}}$ is the frequency of the GW signal’s $m^{\mathrm{th}}$ harmonic in Fourier space). Wen (2003) showed that an approximate expression can be derived for this "peak harmonic" in the form of

{ceqn}

[TABLE]

In this context, GWs of EBBHs are often considered entering the aLIGO band when the GW frequency,

{ceqn}

[TABLE]

reaches $10\>\mathrm{Hz}$ (see e.g. Wen, 2003; Gondán et al., 2018b, and references therein), where both $n_{\mathrm{peak}}=n_{\mathrm{peak}}(e)$ and $a=a(e)$ are functions of $e$. Hence, we also use this criterion in our study.

By substituting Eq. (9) to Eq. (7), the merger time can be expressed as an integral over eccentricities from $e_{\mathrm{LSO}}$ to $e_{0}$ (Gondán et al., 2018b):

{ceqn}

[TABLE]

where

{ceqn}

[TABLE]

GNs are dense stellar environments, therefore there is a possibility that a third object interacts with the binary during the inspiral. Binary-single interactions can be divided into three major categories: weak perturbations, strong perturbations, and close interactions (Samsing et al., 2014). Weak and strong perturbations ultimately lead to the merger of the binary, while close encounters can lead to several different outcomes. According to O’Leary et al. (2009), the characteristic timescale of a binary-single interaction can be approximated as

{ceqn}

[TABLE]

where $n_{\mathrm{tot}}=n_{\mathrm{tot}}(r)$ is the combined number density of all objects residing in a GN (we assume that GNs are spherically symmetric, and $r$ denotes the radial distance from their center). Similarly to Gondán et al. (2018b), we assume that binaries satisfying the condition $t_{\mathrm{enc}}>t_{\mathrm{merge}}$ avoid disruption by a third object, and those that do not are disrupted. Note that this assumption overestimates the number of those binary-single interactions that lead to the disruption of the binary. However, Gondán et al. (2018b) showed that even this assumption results in a negligible number of disrupted binaries.

2.3 The secular Kozai-Lidov mechanism

Triple systems are also believed to play an important role in the formation of EBBHs (see e.g. Miller & Hamilton, 2002; Thompson, 2011; Wen, 2003) through secular evolution (i.e. coherent perturbations on timescales much longer than the orbital period), and specifically through the so-called secular Kozai-Lidov mechanism (Kozai, 1962; Lidov, 1962). Inside a GN, where the dynamics is determined by a central SMBH, binary BHs automatically form a triple system with the SMBH. Previous studies suggested a considerable aLIGO detection rate of $\mathcal{O}$(1 yr*-1*) for mergers of EBBHs formed through the Kozai-Lidov mechanism (see e.g. Antonini et al., 2014; Hoang et al., 2018). Wen (2003) calculated the eccentricity distribution for such binaries merging inside GCs using orbit-averaged approximations, Antonini et al. (2014) however showed the breakdown of these assumptions for triples with high inclination and with an outer BH at a moderate distance. Antonini et al. (2014) and Fragione et al. (2018) also showed that compared to the secular approach, direct N-body simulations predict higher eccentricities for EBBHs entering the aLIGO band. Consistently with other studies (e.g. Antonini & Perets, 2012; VanLandingham et al., 2016), they found that a fraction of these binaries can have $e>0.1$ when they reach the aLIGO band. However, they did not study the formation of such hierarchical triples, but instead, they made assumptions for the distributions of initial parameters.

A triple system consists of an inner binary with masses $m_{1}$ and $m_{2}$ (defined as $m_{1}\geq m_{2}$), and a third object with mass $m_{3}$ (in our case the third object is always the central SMBH). It is convenient to describe the orbit using Jacobi coordinates (see e.g. Murray & Dermott, 2000, p. 441-443). In this treatment, the dominant motion of the triple system can be divided into two separate Keplerian orbits: one is the motion of objects 1 and 2 relative to each other, and the other one is the motion of object 3 relative to the center of mass of the inner binary. We denote the relative position vectors, semi-major axes, and eccentricities of these systems by $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, $a_{1}$ and $a_{2}$, and $e_{1}$ and $e_{2}$, respectively. If the third object is sufficiently distant from the inner binary, a perturbative approach can be used to desribe the evolution of the system. In the usual secular approximation (e.g. Marchal, 1990), the two orbits exchange angular momentum, but not energy. Hence, the eccentricities and orientations of the orbits can change, but the semi-major axes cannot. The coupling term in the Hamiltonian can be written in a power series expansion of parameter $\alpha=a_{1}/a_{2}$ (see e.g. Harrington, 1968), where the Hamiltonian is

[TABLE]

{ceqn}

[TABLE]

$P_{j}$ is the $j$th-degree Legendre polynomial, $\Phi$ is the angle between $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, and

{ceqn}

[TABLE]

When analyzing these hierarchical triple systems, it is particularly convenient to adopt the canonical coordinates known as Delaunay’s elements (e.g. Valtonen & Karttunen, 2006). The coordinates are the mean anomalies, $l_{1}$ and $l_{2}$, the longitudes of ascending nodes, $h_{1}$ and $h_{2}$, and the arguments of periastron, $g_{1}$ and $g_{2}$. One can then eliminate the short-period terms by a canonical transformation to get equations that describe the long-term evolution of the triple system. This technique is known as the Von Zeipel transformation (see e.g. Brouwer, 1959). Naoz et al. (2013a) showed that going beyond the quadrupole-order approximation is necessary in order to acquire highly eccentric Kozai-excitations. They also pointed out a common mistake in the Hamiltonian treatment of these hierarchical systems in previous studies, which could lead to erroneous conclusions in the test particle limit. The full octupole-order equations of motion can be found in Naoz et al. (2013a).

We utilized a fourth-order Runge-Kutta method with adaptive stepsize control to solve the octupole-order equations of motion of Naoz et al. (2013a), including terms from gravitational radiation effects. Fig. 1 shows the eccentricity evolution of the inner binary in a typical hierarchical triple system, where the general relativistic terms become dominant after a few Kozai-cycles. Using the quadrupole-order approximation results in closed, periodic orbits in the phase space of parameters. However, by expanding our equations to the octupole-order, the behaviour of triple systems becomes chaotic. This can lead to extremely high eccentricities for the inner binaries, in which case the gravitational radiation term becomes dominant. As a result, the inner binary ultimately merges. Note that going beyond the octupole-order expansion is unnecessary in most cases, since this approximation already exhibits the system’s chaotic nature, and all higher-order terms take effect on a longer timescale. In addition, yet another drastic change in the qualitative behaviour of the system is not expected, since all the remaining integrals of motion beyond the octupole-expansion are the total energy and the components of the total angular momentum, which remain being integrals of motion even in an exact solution of the equations of motion. Naoz et al. (2013a) also showed a great agreement between the secular approximation and direct N-body integrations in most cases.

3 The Monte Carlo Simulation

In this section we introduce our method of simulating EBBHs in GNs. First, we give an overview of stellar populations residing in GNs in Section 3.1. Then, we introduce the GN model of Bahcall & Wolf (1977), as well as we calculate the probability distributions needed for our Monte Carlo simulations (see Section 3.2). We also describe our Monte Carlo code for simulating EBBHs formed through gravitational capture (see Section 3.3), and introduce our approach for simulating EBBHs formed in hierarchical triple systems in Section 3.4.

3.1 Galactic nuclei and their stellar populations

The main contributors to EBBH formation are expected to be GNs, which are assumed to be spherical, relaxed stellar systems bound by gravity to a SMBH in their center. The stellar populations of GNs consist of main sequence stars (MSs), white dwarves (WDs), neutron stars (NSs), and BHs. An exhaustive investigation of relaxation around SMBHs has been carried out by various studies (e.g. Bahcall & Wolf, 1977; Freitag et al., 2006; Alexander & Hopman, 2009; O’Leary et al., 2009), who found that within a certain radius (called the radius of influence of the SMBH), stellar objects undergo dynamical mass segregation, and form a power-law density profile ($n(r)\propto r^{-\alpha}$, where $\alpha$=$\alpha(m)$, having higher values for larger $m$ masses). The radius of influence is given by the following formula:

{ceqn}

[TABLE]

where $M_{\mathrm{SMBH}}$ is the mass of the SMBH, and $\sigma_{*}$ is the velocity dispersion of the stellar system. The velocity dispersion can be obtained using the $M_{\mathrm{SMBH}}-\sigma_{*}$ relationship (Tremaine et al., 2002; Norris et al., 2014):

{ceqn}

[TABLE]

We used a multi-mass population of SMBHs, where the number distribution of SMBHs can be approximated as

{ceqn}

[TABLE]

where $M_{\bullet}=1.3\times 10^{8}$ $M_{\odot}$, as was shown by the measurements of Aller & Richstone (2002). For the mass range of SMBHs, we chose a lower limit of $10^{5}$ $M_{\odot}$, and we set a conservative upper limit of $2\times 10^{7}$ $M_{\odot}$ (see Barth et al., 2005; Greene & Ho, 2006; Gondán et al., 2018b).

Amongst all the stellar objects in GNs, our main interests are BH populations. There is some uncertainty in the present-day mass function of stellar mass BHs in GNs, since star formation inside GNs can be quite different from the formation of field stars. Hence, there is a significant variation in the mass function suggested by models above $\sim 10$ $M_{\odot}$ (Alexander & Hopman, 2009). Similarly to previous studies (Alexander & Hopman, 2009; Gondán et al., 2018b), we assumed single mass MS, WD, and NS populations with component masses of $1$ $M_{\odot}$, $0.6$ $M_{\odot}$, and $1.4$ $M_{\odot}$, respectively. However, as it was shown by Gondán et al. (2018b), the distributions of EBBH parameters are not sensitive to this choice. We chose the relative fractions of MSs, WDs, and NSs to be $1:0.1:0.01$ (Alexander, 2005). For stellar-mass BHs, we implemented a multi-mass population, assuming a power-law mass function $m_{\mathrm{BH}}^{-\beta}$. The mass distribution of BHs is then

{ceqn}

[TABLE]

where $m_{\mathrm{BH,min}}$ is the lower and $m_{\mathrm{BH,max}}$ is the upper mass limit of BHs, and we chose $\beta=2.35$ similarly to Gondán et al. (2018b). Note that this $\beta$ exponent need not be related to the Salpeter initial mass function (Salpeter, 1955). We set the minimum mass of BHs to $m_{\mathrm{BH,min}}=5$ $M_{\odot}$ based on observations of X-ray binaries (see Belczynski et al., 2012, and references therein), and the maximum mass to $m_{\mathrm{BH,max}}=50$ $M_{\odot}$ (based on Belczynski et al., 2016), consistent with GW observations of binary BHs (Abbott et al., 2018a). The total number of BHs in GNs is proportional to the mass of the central SMBH (Miralda-Escudé & Gould, 2000), therefore, it can be approximated as

{ceqn}

[TABLE]

where $M_{\mathrm{SgrA}^{*}}=4.3\times 10^{6}$ $M_{\odot}$ is the mass of the Sgr A∗ SMBH (Gillessen et al., 2009).

We followed the method of Gondán et al. (2018b) to describe EBBHs formed through gravitational capture. We used the one-body phase space distribution of Bahcall & Wolf (1977), generalized for multi-mass systems (O’Leary et al., 2009) to acquire the distribution of stellar objects within relaxed GNs. We utilized these analytical distributions in our Monte Carlo simulations (see Section 3.3).

3.2 The Bahcall & Wolf model of galactic nuclei

In our simulations we assumed relaxed GNs with spherical symmetry. Hence, we applied the one-body phase space distribution $f(\mathbf{r},\mathbf{v})$ of Bahcall & Wolf (1977) generalized for multi-mass systems (O’Leary et al., 2009), which takes the following form for objects with mass $m$:

{ceqn}

[TABLE]

where $r$ and $v$ denote the magnitudes of the position vector $\mathbf{r}$ and the velocity vector $\mathbf{v}$, respectively,

{ceqn}

[TABLE]

is the Keplerian binding energy per unit mass of an object in the field of the SMBH ($E(r,v)>0$ for all objects inside a relaxed GN), and $C_{m}$ is a normalization constant. The $p_{m}$ exponent is mass-dependent, with which we are able to describe mass segregation. For all stellar objects,

{ceqn}

[TABLE]

where $p_{0}=0.5$ in our simulations, based on the Fokker-Planck simulations of O’Leary et al. (2009). Hence for MSs, WDs, and NSs $p_{m}\approx 0$.

The 3D number density of objects with mass $m$ and at radius $r$ can be obtained from the phase space distribution function:

{ceqn}

[TABLE]

where

{ceqn}

[TABLE]

For MSs of $m=1$ $M_{\odot}$, the number density at the radius of influence emerges from the $M_{\mathrm{SMBH}}-\sigma_{*}$ relationship of Eq. (19), assuming that the stellar mass within the radius of influence of the SMBH is $2M_{\mathrm{SMBH}}$ (see e.g. O’Leary et al., 2009):

{ceqn}

[TABLE]

The number density normalization of WD, NS and BH populations is achieved as it is described in Section 2.2 of Gondán et al. (2018b).

The steady-state description of a stellar population within a GN is only valid inside a certain region, where $r_{\mathrm{min}}<r<r_{\mathrm{max}}$. The maximum radius is the radius of influence of the SMBH defined by Eq. (18). Inside a certain region, called the “loss cone” (Shapiro & Lightman, 1976; Syer & Ulmer, 1999), stellar objects get disrupted by the central SMBH and eventually merge into it. Hence $r_{\mathrm{min}}$ is set by the outer boundary of this loss cone, where the density cusp exhibits a cutoff. We calculated $r_{\mathrm{min}}$ according to the assumptions of Gondán et al. (2018b), requiring that (i) the number distribution of stellar populations reaches a relaxed profile within the age of the galaxy; and (ii) the relaxation time is shorter than the timescale of inspiral into the SMBH (see Appendix A of Gondán et al., 2018b).

The velocity distribution of stellar objects can also be obtained from the phase space distribution function:

{ceqn}

[TABLE]

where $\bar{v}=v/v_{\mathrm{max}}(r)$ is the dimensionless velocity, with $v_{\mathrm{max}}$ being the escape velocity at radius $r$, and $C_{\bar{v}}$ is a normalization factor. In our Monte Carlo simulations, one of the steps was to generate randomized values of relative velocities, hence constructing the distribution of relative velocities was essential. Using the velocity distributions $f_{1}(\mathbf{v}_{1})$ and $f_{2}(\mathbf{v}_{2})$ of objects $1$ and $2$, the distribution of the magnitude of relative velocities $w$ is

[TABLE]

{ceqn}

[TABLE]

where $\delta(\cdot)$ is the Dirac-delta function. Using Eq. (29), the distribution of the magnitude of relative velocities in the Bahcall & Wolf model can be expressed with a single integral (see Appendix B of Gondán et al., 2018b):

[TABLE]

{ceqn}

[TABLE]

where $\bar{w}=w/v_{\mathrm{max}}(r)$ ($0\leq\bar{w}\leq 2$), and $p_{1}$ and $p_{2}$ are the $p_{m}$ exponents corresponding to $m_{1}$ and $m_{2}$, respectively. The remaining integral can be evaluated analytically for any integer or half-integer $p_{1}$ and $p_{2}$, although in a general case it can only be calculated numerically.

3.3 EBBH formation rates and the Monte Carlo method

We simulated the formations of EBBHs through gravitational capture using a Monte Carlo code. In our simulation, we used analytical expressions for EBBH formation rates. The infinitesimal formation rate (measured in s*-1* units) for a fixed $M_{\mathrm{SMBH}}$ is

{ceqn}

[TABLE]

where $1$ and $2$ index the two BHs in the EBBH, and $\sigma=\sigma(w)=\pi[b^{2}_{\mathrm{max}}(w)-b^{2}_{\mathrm{min}}(w)]$ is the cross section for BHs forming an EBBH through gravitational capture (see Eq. (6) for $b_{\mathrm{min}}$ and $b_{\mathrm{max}}$). We used the results of Gondán et al. (2018b), who carried out the analytical calculations to determine the distribution of EBBH formation rate as a function of different parameters. Compared to the simulations of Gondán et al. (2018b), we expanded our parameter space by picking the $M_{\mathrm{SMBH}}$ values randomly too. The formation rate as a function of $M_{\mathrm{SMBH}}$ can be calculated as

{ceqn}

[TABLE]

where $\left\langle\Gamma\right\rangle_{M_{\mathrm{SMBH}}}$ is the formation rate of EBBHs in a single GN with $M_{\mathrm{SMBH}}$ in its centre. We show the distribution of the EBBH formation rate (red dashed line) and the SMBH number distribution (blue solid line) as functions of $M_{\mathrm{SMBH}}$ in Fig. 2. Note that the formation rate distribution is mainly determined by the mass distribution of SMBHs and is only weakly affected by the $M_{\mathrm{SMBH}}$ dependency of the formation rate of EBBHs in a single GN (see the black dash-dotted line in Fig. 2). A thorough investigation on the $M_{\mathrm{SMBH}}$ dependence of the formation rate of EBBHs in a single GN gives $\left\langle\Gamma\right\rangle_{M_{\mathrm{SMBH}}}\propto M_{\mathrm{SMBH}}^{\gamma}$ with $\gamma\approx 0.1-0.25$ for single-mass BH populations with different BH masses (see Appendix B of Gondán et al., 2018b).

In our Monte Carlo simulation, we applied the following steps. First, we used the EBBH formation rate distribution shown in Fig. 2 to generate a number of random $M_{\mathrm{SMBH}}$ values, thus simulating GNs with different SMBH masses. Then for each $M_{\mathrm{SMBH}}$, we drew an $m_{1}$-$m_{2}$ pair from their corresponding two-dimensional distribution, where both masses range from $5$ $M_{\odot}$ to $50$ $M_{\odot}$ (see Appendix B.3 in Gondán et al., 2018b, for details). For each $m_{1}$-$m_{2}$ pairs, we randomly drew $r$ radii values from the range [$r^{1,2}_{\mathrm{min}},r_{\mathrm{max}}$] (see Eq. (18) and Appendix A of Gondán et al., 2018b, for details). Lastly, for each $r$ we randomly drew pairs of $w$ and $b$ values from ranges [$0,2v_{\mathrm{max}}$] and [$b_{\mathrm{min}},b_{\mathrm{max}}$], respectively.

The analytical estimates we used in our simulations can be violated by EBBHs that interact with a third object between their formation and merger. Gondán et al. (2018b) showed that this affects only a small fraction ($<1$ percent) of binaries over the considered ranges of parameters. EBBHs can form with gravitational capture if the impact parameter of the encounter between the BHs is smaller than $b_{\mathrm{max}}$ (so that enough energy is emitted during the encounter for the binary to form), and larger than $b_{\mathrm{min}}$ (to avoid the direct collision of the BHs). These parameters scale as $b_{\mathrm{min}}\propto w^{-1}$ and $b_{\mathrm{max}}\propto w^{-9/7}$ with the relative velocity of the BHs. As a consequence, in the inner regions of GNs, near $r_{\mathrm{min}}$, the condition $b_{\mathrm{max}}>b_{\mathrm{min}}$ can be violated. We discarded all such $(w,b)$ pairs, however, as it was shown by Gondán et al. (2018b), this corresponds only to a small fraction ($<1$ percent) of binaries over the ranges of parameters we considered in our simulations.

3.4 Numerical approach for hierarchical triples

The formation of hierarchical triple systems of BHs is hardly manageable with analytical calculations, since they can form through various different processes, most of which include multiple interactions between binary systems that can easily become untraceable from an analytical point of view. Therefore, we applied a Monte Carlo approach different from the one employed for gravitational captures. We adopted a set of distributions to set the initial parameters of triple systems. We then used these initial conditions as inputs for our differential equation integrator, and solved the octupole-order equations of motion of Naoz et al. (2013a), in which we included gravitational radiation. For solving the differential equations we utilized a fourth-order Runge-Kutta differential equation integrator with adaptive stepsize control to increase the speed of the integration while not losing precision.

We picked the masses of BHs making up the inner binary independently from the distribution defined by Eq. (21) using an exponent of $\beta=2.35$. The third object in the triple system was always chosen as the central SMBH. The mass of the SMBH was picked from the distribution defined by Eq. (20). The lower and upper limits for masses were chosen to be the same as for gravitational captures. The initial semi-major axis of the inner binary, $a_{1}$, was drawn uniformly in the logarithmic domain (which corresponds to d$N/$d$a\propto 1/a$). The minimum and maximum values of $a_{1}$ were chosen as $0.1$ AU and $50$ AU, respectively, similarly to the setup of Hoang et al. (2018). The initial semi-major axis of the outer binary was drawn from a uniform distribution (which corresponds to a 3D number density distribution of $n(r)\propto r^{-2}$) between $100$ AU and $0.1$ pc, based on the simulational setup of Hoang et al. (2018). Initial eccentricities of the inner and outer binaries were drawn between [math] and $1$ from uniform (Raghavan et al., 2010), and thermal distributions (d$N/$d$e\propto e$, see Jeans, 1919), respectively. $g_{1}$ and $g_{2}$ were picked randomly from uniform distributions, while the initial inclination $I$ was drawn from a distribution d$N/$d$I\propto\cos(I)$ over the range of $[65^{\circ},125^{\circ}]$. The range of initial inclination was based on the results of Silsbee & Tremaine (2017), who found that mergers typically occur with $-0.5<\cos(I)<0.5$.

Since we use secular approximations when describing triple systems, we require our systems to satisfy dynamical stability in order to justify the hierarchical treatment. Hence, we use two stability criteria (see e.g. Stephan et al., 2016; Hoang et al., 2018). First we require that the relative strengths of the octupole and quadrupole terms fulfil the following criterion (Naoz, 2016):

{ceqn}

[TABLE]

The second condition requires that the inner binary does not cross the SMBH’s Roche limit (e.g. Naoz & Silk, 2014):

{ceqn}

[TABLE]

where $M=m_{1}+m_{2}$ is the total mass of the inner binary. We discarded all triple systems that failed to satisfy these conditions.

The highly eccentric excitations can be inhibited by the general relativistic precession of the inner binary, as was shown by Naoz et al. (2013b). This happens when the timescale of Kozai-Lidov oscillations at the quadrupole level of approximation ($t_{\mathrm{quad}}$, see e.g. Antognini, 2015) is much longer than the timescale of general relativistic precession ($t_{\mathrm{GR,inner}}$, see e.g. Naoz et al., 2013b). We consider this phenomenon in our simulations by discarding all triple systems for which $t_{\mathrm{GR,inner}}<t_{\mathrm{quad}}$, as the inner binaries of these systems may not experience highly eccentric excitations.

We terminated our simulation if at any time $a_{1}$ became smaller than $10^{-4}$ AU, which we assumed to guarantee a merger within a short period of time. Close encounters with other stellar objects inside the GN can unbind the inner binary during its evolution. This evaporation timescale is

{ceqn}

[TABLE]

where in the inner parts ($\leq 0.1$ pc)

{ceqn}

[TABLE]

is the density of stars, $m_{\mathrm{ave}}=1$ $M_{\odot}$ is the average mass of background stars, $\ln\Lambda=15$ is the Coulomb logarithm, $\sigma=\sqrt{0.1\>\mathrm{pc}/a_{2}}$ $280$ km/s is the velocity dispersion, $G$ is the universal gravitational constant, while $\sigma_{0}=200$ km/s and $M_{0}=3\times 10^{8}$ $M_{\odot}$ are constants (see Hoang et al., 2018). We simply terminated our simulation and discarded the triple system if the inner binary has not produced a merger after $t_{\mathrm{ev}}$. The inner binary can merge without experiencing high-eccentricity excitations. However, we are only interested in those mergers that take place due to Kozai-excitations. Therefore we discarded those triple systems, in which the merger did not happen due to high-eccentricity excitations. In addition to these conditions, we kept only those EBBHs, for which effects due to GW emission dominated over secular effects at $f_{\mathrm{GW}}=10$ Hz, and this condition held throughout their evolution from that point onward (see Section 4.1 for details).

4 Results

In this section we present our results. First, we discuss how we acquired the distributions of EBBH eccentricities at $f_{\mathrm{GW}}=10$ Hz for all formed EBBHs by evolving EBBH orbits in our simulations described in Section 3.3 and 3.4. We then construct the distributions of eccentricities for future EBBH observations with design aLIGO detectors (see Section 4.1). In Section 4.2 and 4.3 we also show that if future EBBH detection rates with aLIGO will be dominated by EBBHs formed in GNs, then it is feasible to constrain the branching ratios between the two main channels. We performed a statistical test to determine the minimum number of EBBH observations needed to constrain the branching ratio of gravitational capture to an arbitrarily chosen $\Delta\mathcal{F}\mathrm{{}_{GC}}=0.2$ wide one-sigma confidence interval (note that due to the underlying symmetry, our results would be the same if the branching ratio for the secular Kozai-Lidov mechanism was constrained instead). We present our results without accounting for errors of the reconstruction of EBBH orbital parameters in Section 4.2. We then discuss the implications of non-zero parameter reconstruction errors in Section 4.3.

4.1 Eccentricity distributions at $f_{\mathrm{GW}}=10$ Hz

We produced the distributions of eccentricities of EBBHs at the time their peak GW frequency, $f_{\mathrm{GW}}$, reaches $10$ Hz, from now on denoted as $e_{\mathrm{LIGO}}$. For EBBHs formed through gravitational capture, we produced the $e_{\mathrm{LIGO}}$ distribution for all formed EBBHs by solving Eqs. (9) and (12) for $f_{\mathrm{GW}}=10$ Hz, using initial orbital parameters obtained from our simulation described in Section 3.3. We also acquired the $e_{\mathrm{LIGO}}$ distribution for EBBHs formed through the Kozai-Lidov mechanism from our simulations described in Section 3.4. Up until the last few Kozai-cycles before the merger, the evolution of the inner binary is dominated by perturbing effects of the outer object, after which the domination is taken over by the energy and angular momentum loss due to gravitational radiation (see Fig. 1 as an example). According to our results, the proportion of EBBHs formed through the Kozai-Lidov mechanism, for which secular effects still cannot be ignored at the time their GW signals enter the aLIGO band is $<1$ percent. For such EBBHs, both finding the GW signal and reconstructing EBBH parameters pose a complex problem. Thus, for simplicity, we leave out these EBBHs from our analysis.

The distribution of physical parameters is different for detected EBBHs than for formed EBBHs, hence we needed to produce the $e_{\mathrm{LIGO}}$ distribution for all detected EBBHs in order to carry out our analysis. This is not equivalent with the $e_{\mathrm{LIGO}}$ distribution of all formed EBBHs, which must be weighted by the volume from which EBBHs can be observed, and the maximum distance to which EBBHs are observable (called the horizon distance, which we denote by $D_{\mathrm{H}}$) differs for EBBHs with different component masses and eccentricities. O’Leary et al. (2009) showed that BH binaries with higher initial eccentricities are detectable only to slightly larger distances than circular binaries for low-mass ($M<60$ $M_{\odot}$) binaries, and the difference in the horizon distance between circular and eccentric binaries becomes substantial only at higher masses ($M\approx 100$ $M_{\odot}$). For EBBHs formed through the secular Kozai-Lidov mechanism the proportion of EBBHs with $M>60$ $M_{\odot}$ is only $\sim 1$ percent. Therefore for simplicity, for EBBHs formed through this channel, we defined $D_{\mathrm{H}}$ the same way as it is defined for circular binaries (see e.g. Abadie et al., 2010):

{ceqn}

[TABLE]

where

{ceqn}

[TABLE]

is the chirp mass, $f_{\mathrm{low}}$ is set by the detector noise cutoff (which in our case is $f_{\mathrm{low}}=10$ Hz), $f_{\mathrm{high}}=(6\sqrt{6}\pi M)^{-1}$ is the frequency at the innermost stable circular orbit, and $S_{n}(f)$ is the power spectral density of the detector noise, for which we used the design aLIGO projection (see Aasi et al., 2015a). For EBBHs formed through gravitational capture, we used the horizon distance calculated for eccentric binaries in Eq. (58) of O’Leary et al. (2009), since in this channel the proportion of binaries with total masses of $M>60$ $M_{\odot}$ is $\sim 10$ percent. The number of detectable EBBHs with certain parameter values should be proportional to the volume enclosed by the corresponding horizon distance. Hence, to convert the $e_{\mathrm{LIGO}}$ distribution of all formed EBBHs to the distribution of all detected EBBHs, we weighted each EBBH parameter set by $V_{\mathrm{H}}=4\pi D_{\mathrm{H}}^{3}/3$. The horizon distance $D_{\mathrm{H}}$ as a function of BH masses for low-mass binaries ($m_{2}\leq m_{1}\leq 30$ $M_{\odot}$, defined by Eq. (38)) is shown in Fig. 3. The binaries with the smallest masses are detectable to $\sim 1.3$ Gpc, while binaries with $M=60$ $M_{\odot}$ have a corresponding horizon distance of $\sim 4.4$ Gpc.

The $e_{\mathrm{LIGO}}$ distributions of detected EBBHs for the two different formation channels are shown in Fig. 4. For EBBHs formed through the secular Kozai-Lidov mechanism, we show the $e_{\mathrm{LIGO}}$ probability distribution functions resulting from both the quadrupole- and octupole-order equations of motion. Using the octupole-order equations generally produces EBBHs with higher eccentricities. This is understandable, since we expect EBBHs to begin their inspiral from higher initial eccentricities due to extreme Kozai-excitations caused by the octupole term. Note that the $e_{\mathrm{LIGO}}$ distributions for detected EBBHs only slightly differ from the corresponding distributions for all formed EBBHs (see the right panel of Fig. 4). We also find, supplementing the results of Hoang et al. (2018), that about 10 percent of binary BHs formed through the Kozai-Lidov mechanism in GNs have eccentricities $e_{\mathrm{LIGO}}>0.1$ (note that the range and distribution of BH masses in our simulations were different from that of Hoang et al., 2018, and we did not use explicitly the 1PN terms describing orbital precession, these, however, should not change the percentage value significantly).

For EBBHs formed through gravitational capture two main groups can be distinguished in the cumulative distribution function (see the right panel of Fig. 4). More than half of the EBBHs form outside the sensitive band of aLIGO detectors, with the initial frequency of the GW signal being $f_{\mathrm{GW,0}}<10$ Hz. The eccentricities of these binaries at $f_{\mathrm{GW}}=10$ Hz fall between $10^{-2}<e_{\mathrm{LIGO}}<1$. The other group consists of EBBHs that form inside the sensitive band of aLIGO detectors, with $f_{\mathrm{GW,0}}\geq 10$ Hz. These binaries form a sharp peak in the probability distribution function at eccentricities $e_{\mathrm{LIGO}}>0.95$. For binary BHs formed through gravitational capture in GNs, we find that about 75 percent have eccentricities $e_{\mathrm{LIGO}}>0.1$. Table 1 gives further information about the fraction of EBBHs under various conditions (similarly to Table 1 in Gondán et al., 2018b).

EBBHs tend to have higher eccentricities at lower signal frequencies. Hence we may also be able to constrain these formation channels using future GW observations with detectors operating at lower frequencies (i.e. with LISA or DECIGO). However, as it was recently shown by several studies, binary black holes that enter the aLIGO band with $e>5\times 10^{-3}$ may completely elude the LISA band in the epoch their evolution is dominated by GW emission (see e.g. Antonini et al., 2017; Samsing et al., 2018; D’Orazio & Samsing, 2018, and references therein). EBBHs formed in hierarchical triples may still enter the LISA band multiple times during the time secular evolution dominates over gravitational radiation. The parameter reconstruction of these systems, however, pose a more complex problem than in the case of an EBBH in its merger phase. The proportion of EBBHs formed through the secular Kozai-Lidov mechanism, for which gravitational radiation does not dominate over secular effects of the outer object in the DECIGO band, may also be considerably higher.

The $e_{\mathrm{LIGO}}$ distributions of detected EBBHs determined in this section are idealized forms of the distributions of observed $e_{\mathrm{LIGO}}$ values, since they do not take into account the potential errors of the parameter estimation. We considered a set of plausible errors consistent with the results of previous studies (see Gondán et al., 2018a; Lower et al., 2018; Gondán & Kocsis, 2019) to account for the magnitudes of errors on the parameter estimations in our statistical tests (see Section 4.3).

4.2 Constraining EBBH formation channels with zero parameter estimation errors

While this may not be true (see Section 1 for further explanations), in this section and in Section 4.3 we assume that future EBBH detection rates with aLIGO will be dominated by EBBHs formed in GNs. We perform a statistical test to demonstrate the feasibility of constraining the branching ratios between gravitational capture and the secular Kozai-Lidov mechanism in GN host environments, when only these formation channels are taken into account. To carry out our analysis, we needed to produce the distributions of physically measurable EBBH parameters. In principle, reconstructible parameters of non-spinning EBBHs include sky location, BH masses, spins, and orbital parameters ($e$, $\rho_{\mathrm{p}}$) at a given point of EBBH evolution. In realistic scenarios, however, the reconstruction of these parameters is liable to a number of limitations. First of all, the reconstruction of EBBH parameters may involve degeneracies. For example, for EBBHs formed through gravitational capture, initial orbital parameters may only be reconstructed if the peak GW frequency of the initial orbit is large enough to be in the sensitive frequency band of GW detectors (see Gondán et al., 2018a). Parameter estimation errors are also inherently involved in the reconstruction due to the finite sensitivity of GW detectors (see e.g. Bécsy et al., 2017). These errors may also be affected by theoretical uncertainties of the waveform model (Cutler & Vallisneri, 2007). Kocsis & Levin (2012) showed that post-Newtonian calculations of different order lead to substantially different orbital evolutions at small separations ($\rho_{p}<20$), although they do not deviate considerably from the leading order results for higher separations. On the other hand, as already mentioned in Section 4.1, the reconstruction of parameters of EBBHs formed through the secular Kozai-Lidov mechanism may only be accurate if the inner binary has achieved a low enough separation to have its evolution governed by gravitational radiation instead of secular effects.

Considering these factors, we chose the parameter most suitable for our analysis to be $e_{\mathrm{LIGO}}$, since at $f_{\mathrm{GW}}=10$ Hz the two BHs have a high enough separation that the post-Newtonian calculations are still accurate, while they have already shrunk to an orbit, where secular effects caused by a third object can be neglected. We chose the distribution of an EBBH orbital parameter over distributions of BH masses, since distributions of orbital parameters only weakly depend on the choice of the distribution of BH masses (see Gondán et al., 2018b). The results of Kocsis & Levin (2012) suggest that the choice of $f_{\mathrm{GW}}=10$ Hz corresponds to a point of EBBH evolution where the leading-order approximations are still highly accurate. On the other hand, we kept only those EBBHs in hierarchical triples, for which effects due to GW emission dominated over secular effects at $f_{\mathrm{GW}}=10$ Hz, and this condition held throughout their evolution from that point onward.

We carried out our statistical test in the following manner. We constructed a mixed distribution of $e_{\mathrm{LIGO}}$ values by weighting the $e_{\mathrm{LIGO}}$ PDF for gravitational captures (see Fig. 4) by a factor of $\mathcal{F}\mathrm{{}_{GC}}$, and the $e_{\mathrm{LIGO}}$ PDF for the secular Kozai-Lidov mechanism by a factor of $1-\mathcal{F}\mathrm{{}_{GC}}$, and then summing them up together. We then picked $N$ number of $e_{\mathrm{LIGO}}$ values randomly from this mixed distribution, and maximized the likelihood of this $N$-element sample over all possible values of $\mathcal{F}\mathrm{{}_{GC}}\in[0,1]$ using the Anderson-Darling test (Anderson & Darling, 1952). We denote the resulting $\mathcal{F}\mathrm{{}_{GC}}$ value as $\mathcal{F}\mathrm{{}_{GC}^{obs}}$. We repeated this process $10^{5}$ times for a given $N$, to map the PDF of $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ for a fixed value of $\mathcal{F}\mathrm{{}_{GC}}$. For each $N$ value, we performed this test for different values of $\mathcal{F}\mathrm{{}_{GC}}$ ranging from $\mathcal{F}\mathrm{{}_{GC}}=0$ to $\mathcal{F}\mathrm{{}_{GC}}=1$, obtaining the corresponding PDFs of $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ values. We then determined the PDF of $\mathcal{F}\mathrm{{}_{GC}}$ for a fixed value of $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ using conditional probability, and determined the 68 percent and 90 percent confidence intervals of these distributions for different values of $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ ranging from $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.1$ to $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.9$. We also repeated the whole process for different values of $N$ and determined the minimum number of EBBH observations needed to constrain $\mathcal{F}\mathrm{{}_{GC}}$ to an arbitrarily chosen $\Delta\mathcal{F}\mathrm{{}_{GC}}=0.2$ wide $68$ percent (one-sigma) confidence interval (note that for symmetrical errors $\Delta\mathcal{F}\mathrm{{}_{GC}}=0.2$ corresponds to a $\pm 10$ percent range of the total $[0,1]$ interval).

Fig. 5 summarizes our results. We find that the PDFs of $\mathcal{F}\mathrm{{}_{GC}}$ for fixed values of $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ peak at $\mathcal{F}\mathrm{{}_{GC}^{obs}}$, which means that our statistic is unbiased. The widths of the 68 percent and 90 percent intervals decrease with increasing values of $|\mathcal{F}\mathrm{{}_{GC}^{obs}}-0.5|$, hence the widest intervals correspond to $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.5$. Also note that the plots for different values of $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ are symmetric to $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.5$, which means that they would practically be the same if instead of gravitational capture, the secular Kozai-Lidov mechanism was considered. As expected, the 68 percent and 90 percent confidence intervals shrink with increasing numbers of simulated observations. The widths of the 68 percent confidence intervals reach $\Delta\mathcal{F}\mathrm{{}_{GC}}=0.2$ at $N=21$ for $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.1$ (and for $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.9$), at $N=43$ for $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.3$ (and for $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.7$), and at $N=46$ for $\mathcal{F}\mathrm{{}_{GC}^{obs}}=0.5$. Therefore, we conclude that in the idealized case of having no parameter reconstruction errors, a few tens of EBBH observations are sufficient to constrain $\mathcal{F}\mathrm{{}_{GC}}$ to a $\Delta\mathcal{F}\mathrm{{}_{GC}}=0.2$ wide one-sigma confidence interval.

4.3 Results with non-zero parameter estimation errors

In a realistic situation, reconstructed $e_{\mathrm{LIGO}}$ values are liable to estimation errors. We investigated their effect in our statistical test (described in Section 4.2) by assuming two different magnitudes for estimation errors on $e_{\mathrm{LIGO}}$ (from now on, denoted by $\Delta e_{\mathrm{LIGO}}$), which we take into account by exchanging each no-error $e_{\mathrm{LIGO}}$ value with a randomized value taken from a Gaussian distribution centered on $e_{\mathrm{LIGO}}$ and with $\sigma=\Delta e_{\mathrm{LIGO}}$ (note that these Gaussian distributions were truncated outside the allowed $e_{\mathrm{LIGO}}\in[0,1]$ range). The two $\Delta e_{\mathrm{LIGO}}$ values we used in our investigations were $\Delta e_{\mathrm{LIGO}}=10^{-2}$ and $\Delta e_{\mathrm{LIGO}}=10^{-1}$. Choosing $\Delta e_{\mathrm{LIGO}}=10^{-2}$ is supported by predictions on plausible errors for observations of EBBHs at $\sim 1$ Gpc given in Gondán et al. (2018a) and Gondán & Kocsis (2019). Thus $\Delta e_{\mathrm{LIGO}}=10^{-1}$ is an overestimation of the $e_{\mathrm{LIGO}}$ estimation errors realistically achievable in the future.

By repeating the analysis described in Section 4.2 with $N$-element samples of such non-zero error $e_{\mathrm{LIGO}}$ values, we have found that $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ becomes a biased estimator of $\mathcal{F}\mathrm{{}_{GC}}$, even when $N$ goes to infinity. This is due to the fact that the $e_{\mathrm{LIGO}}$ values estimated with errors do not follow the zero-error $e_{\mathrm{LIGO}}$ distributions shown in Fig. 4 that we use as bases in our statistical test. For an example case of choosing $N=50$, the peak of the $\mathcal{F}\mathrm{{}_{GC}}-\mathcal{F}\mathrm{{}_{GC}^{obs}}$ distribution ($\Delta\mathcal{F}\mathrm{{}_{GC}^{p}}$) shifts around $(\mathcal{F}\mathrm{{}_{GC}}-\mathcal{F}\mathrm{{}_{GC}^{obs}})=0$ by a varying amount that is a function of $\mathcal{F}\mathrm{{}_{GC}^{obs}}$. Fig. 6 shows the amount of this shift as a function of $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ for both $\Delta e_{\mathrm{LIGO}}=10^{-2}$ and $\Delta e_{\mathrm{LIGO}}=10^{-1}$.

We have also found that the widths of the 68 percent confidence intervals of the $\mathcal{F}\mathrm{{}_{GC}}$ distributions do not change significantly when estimation errors are introduced, but remain essentially the same as in the zero-error case described in Section 4.2. Since the dependence of $\Delta\mathcal{F}\mathrm{{}_{GC}^{p}}$ on $\mathcal{F}\mathrm{{}_{GC}^{obs}}$ is heavily dependent on $\Delta e_{\mathrm{LIGO}}$, we can only draw conclusions on the magnitude of these shifts for the different $\Delta e_{\mathrm{LIGO}}$ values. For example, we can compare the magnitude of the shifts to the one-sigma confidence intervals of the $\mathcal{F}\mathrm{{}_{GC}}$ distributions in the zero-error case (shown by a shaded green area in Fig. 6), which closely resembles the corresponding one-sigma confidence intervals for the two non-zero error cases. We have found that with $\Delta e_{\mathrm{LIGO}}=10^{-1}$ the shift of the peak can reach $\Delta\mathcal{F}\mathrm{{}_{GC}^{p}}>0.2$. However with $\Delta e_{\mathrm{LIGO}}=10^{-2}$, this shift remains within the one-sigma confidence interval of statistical errors, with a maximum shift of $|\Delta\mathcal{F}\mathrm{{}_{GC}^{p}}|\simeq 0.06$ when $\mathcal{F}\mathrm{{}_{GC}^{obs}}\in[0.1,0.85]$. Hence, we conclude that in practice $e_{\mathrm{LIGO}}$ estimation errors should be kept at the $\Delta e_{\mathrm{LIGO}}\leq 10^{-2}$ level to be able to constrain the branching ratios between different EBBH formation channels with this method to a $\Delta\mathcal{F}=0.2$ wide one-sigma confidence interval.

5 Conclusions

We calculated the probability density functions of $e_{\mathrm{LIGO}}$ ($e$ at $f_{\mathrm{GW}}=10$ Hz) of EBBHs formed in GNs through gravitational capture and through the secular Kozai-Lidov mechanism (see Section 4.1), using Monte Carlo simulations described in Section 3.3 and 3.4. We have found, supplementing the results of Hoang et al. (2018), that $\sim 10$ percent of binary BHs formed through the Kozai-Lidov mechanism in GNs have eccentricities $e_{\mathrm{LIGO}}>0.1$ (see Section 4.1), and that $\sim 75$ percent of binary BHs formed through gravitational capture in GNs have eccentricities $e_{\mathrm{LIGO}}>0.1$. Our simulations of EBBHs formed through gravitational capture extend the simulations of Gondán et al. (2018b) by introducing randomized masses for the central SMBHs. It is important to note that the $e_{\mathrm{LIGO}}$ distributions presented here may be different if we assume a different underlying model for GNs, for example if we assume that BHs form disks inside GNs as Szölgyén & Kocsis (2018) has suggested, or if we assume a different radial dependence of the number density of BHs compared to the one we considered in Section 3.2.

Beyond our main results described in the previous paragraph, we have also carried out a statistical test that can be used in the future to constrain the branching ratios (denoted by $\mathcal{F}\in[0,1]$ for a certain formation channel) between the two main EBBH formation channels in GNs with future GW observations of EBBHs (see Section 4.2 and 4.3). Using these tests we have found that with design aLIGO, a few tens of observations are sufficient for every observed value of $\mathcal{F}\mathrm{{}^{obs}}$ to constrain $\mathcal{F}$ to a $\Delta\mathcal{F}=0.2$ wide one-sigma confidence interval (see Section 4.2 for details). We also supplemented our study with additional tests in which we considered non-zero parameter estimation errors. We have found that $\mathcal{F}\mathrm{{}^{obs}}$ becomes a biased estimator of $\mathcal{F}$ even if the number of EBBH observations approaches infinity. We conclude that $e_{\mathrm{LIGO}}$ estimation errors should be kept at the feasible $\Delta e_{\mathrm{LIGO}}\leq 10^{-2}$ level (Gondán et al., 2018a; Gondán & Kocsis, 2019) to be able to constrain the branching ratios between the two main channels to a $\Delta\mathcal{F}=0.2$ wide one-sigma confidence interval (see Section 4.3 for details).

The expected detection rate of EBBHs formed through gravitational capture is $\mathcal{O}(1-100$ yr*-1*) with aLIGO detectors operating at design sensitivity (O’Leary et al., 2009). This aLIGO detection rate may also be higher than $\sim$100 yr*-1* if the BH mass function extends to masses over $25$ $M_{\odot}$ (O’Leary et al., 2009). Several LIGO and Virgo observations of such high mass BHs have been made recently (e.g. Abbott et al., 2016, 2017a, 2017b). The detection rate of EBBHs formed through the secular Kozai-Lidov mechanism is expected to be of the same order as the detection rate of EBBHs formed through gravitational capture. Based on these detection rates, the minimum observation time needed to constrain the branching ratios between these two channels to a $\Delta\mathcal{F}=0.2$ wide one-sigma confidence interval is $\mathcal{O}(0.1-10$ yr), if $e_{\mathrm{LIGO}}$ estimation errors can be kept at the feasible $\Delta e_{\mathrm{LIGO}}\leq 10^{-2}$ level. Note that our tests were carried out while assuming that future observations of EBBHs will be dominated by EBBHs formed in GNs, and thus our results should be treated with due caution.

Note that during our analysis we assumed the search pipeline applied in the future for finding GW signals of EBBHs and reconstructing their parameters is equally sensitive for different eccentricities. This assumption may be idealistic, and must be tested in the future when such pipelines become available. Nevertheless, the selection effect that the limited sensitivity of a specific pipeline introduces can be taken into account as corrections in our predictions for the detected $e_{\mathrm{LIGO}}$ distributions, and thus in our statistical tests as well.

Acknowledgements

We would like to thank László Gondán and Bence Kocsis for their valuable comments on the manuscript. We would also like to thank László Gondán for providing us a set of numerical codes we used in our simulations. János Takátsy, Bence Bécsy, and Péter Raffai were supported by the ÚNKP-16-1 (J.T.), ÚNKP-17-2 (B.B.), and ÚNKP-17-4 (P.R.) New National Excellence Programs of the Ministry of Human Capacities.