Black Hole Mergers in Galactic Nuclei Induced by the Eccentric Kozai-Lidov Effect

I. Introduction

Recently, the Advanced Laser Interferometer Gravitational-Wave Observatory http://www.ligo.org/ (LIGO) has directly detected gravitational waves (GWs) from at least five inspiraling black hole-black hole (BH-BH) binaries in the Universe (LIGO Scientific and Virgo Collaboration 2016a; LIGO Scientific and Virgo Collaboration 2016b; LIGO Scientific and Virgo Collaboration 2017; The LIGO Scientific Collaboration et al. 2017; Abbott et al. 2017). With ongoing further improvements to LIGO and the commissioning of additional instruments VIRGO http://www.virgo-gw.eu/ and KAGRA http://gwcenter.icrr.u-tokyo.ac.jp/en/, hundreds of BH-BH binary sources may be detected within the decade, opening the era of gravitational wave astronomy (The LIGO Scientific Collaboration et al. 2016).

The primary astrophysical origin of BH-BH mergers is still under debate. Merging BH binaries may be produced dynamically in dense star clusters such as globular clusters or nuclear star clusters at the center of galaxies (Portegies Zwart & McMillan 2000; Wen 2003; O’Leary et al. 2006; Antonini et al. 2014; Rodriguez et al. 2016b; O’Leary et al. 2016; O’Leary et al. 2009; Kocsis & Levin 2012; Bartos et al. 2016; Stone et al. 2016). These merging BH binaries may also be a result of isolated binary evolution in the galactic field due to special modes of stellar evolution and evolution in active galactic nuclei (Mandel & de Mink 2016; de Mink & Mandel 2016; Belczynski et al. 2016a; Marchant et al. 2016). Some studies also suggested that these LIGO detections are sourced in the first stars (Kinugawa et al. 2014; Kinugawa et al. 2016; Hartwig et al. 2016; Inayoshi et al. 2016; Dvorkin et al. 2016), cores of massive stars (Reisswig et al. 2013; Loeb 2016; Woosley 2016), or dark matter halos comprised of primordial black holes (Bird et al. 2016; Clesse & García-Bellido 2016; Sasaki et al. 2016).

Here we focus on sources produced in galactic nuclei around massive black holes (MBHs). The large escape speed in nuclear star clusters (NSCs) creates an ideal environment to accumulate a population of stellar-mass BHs. Thus, even stellar-mass BH binaries can survive in this configuration despite of supernova kicks (Lu & Naoz, in prep.). Antonini & Rasio 2016 investigated the evolution of binaries in NSCs without MBHs in their centers, and found a rate of $1.5\,\rm Gpc^{-3}\rm yr^{-1}$. Furthermore, O’Leary et al. 2009 considered binaries that form due to gravitational wave emission during close encounters between single BHs (Kocsis & Levin 2012; Gondán et al. 2017, see also). In this channel, the merger rate is dominated by massive BHs over $25\rm M_{\odot}$, which are delivered to the center by dynamical friction and relax to form steep density cusps near the MBH. Mergers of low mass BHs, and neutron stars may be rare in this channel (Tsang 2013).

We investigate the secular evolution of stellar-mass BH binaries in galactic nuclei which include an MBH in their centers. These BH-BH binaries undergo large amplitude eccentricity oscillations due to the Eccentric Kozai-Lidov (Naoz 2016, EKL, e.g.,) mechanism in the presence of the MBH (Antonini et al. 2010; Antonini & Perets 2012, e.g.,). If the eccentricity reaches a sufficiently high value, GW emission drives the binary to merge.Antonini & Perets 2012 studied the merger rate of BH-BH binaries in the presence of an MBH in galactic nuclei and estimated the merger rates to be $1.7-4.8\times 10^{-4}$ Myr-1. VanLandingham et al. 2016 studied the merger rate of BH-BH binaries in the presence of intermediate mass black holes $\sim 10^{3-4}\rm M_{\odot}$ and found it to be very high. Here we also focus on BH-BH mergers in the presence of an MBH in galactic nuclei. However, unlike Antonini & Perets 2012, which assumes that the heaviest BH mass in the cluster is $10\,\rm M_{\odot}$, we explore a large range of BH masses because recent LIGO observations have shown that there is a much greater range of BH masses than previously thought. EKL effects are stronger when there is a greater mass difference between the two BHs in the binary. Thus, Antonini & Perets 2012 concluded that EKL was not important and based their rate estimate on the semi-analytical timescale given in Thompson 2011, using only the quadrupole order in a multipole expansion. The range in the rate estimation of Antonini & Perets 2012 represents core to mass-segregated distributions of the binaries, and corresponds to $\sim 0.002-0.48$ Gpc-3 yr-1. Here we use the EKL mechanism which represents the secular approximation up to the octupole level of approximation. Furthermore, we adopt a BH population consistent with a BH mass function that extends to higher masses (O’Leary et al. 2009; Kocsis & Levin 2012).

The octupole correction drives chaotic variations in the eccentricity, increases the probability of close encounters, and thus enhances the efficiency of BH-BH mergers. We show that this mechanism results in a merger rate which is higher than the rates found in the case without MBH, and it is coincidentally comparable to estimated globular cluster rates.

We describe our simulations in Section II, and present our results, predictions and merger rate in Section III. Finally, we offer our discussions in Section IV.

II. Numerical setup

We run several large sets of Monte-Carlo simulations to investigate the effects of the SMBH’s gravitational perturbations on binary BHs, with variations between sets of simulations to account for different effects (see Table 1 for details of all simulations). We study the secular dynamical evolution of binary BHs around the MBH in galactic nuclei, starting from the BH binary phase. We include the secular equations up to the octupole-level of approximation (Naoz 2016, e.g.,), general relativity precession of the inner and outer orbits (Naoz et al. 2013, e.g.,), and gravitational wave emission (Peters 1964). We consider this as a simple proof of concept to investigate the effects of the EKL mechanism Note that we do not assume any specific mechanism for the formation or delivery of those binaries. We later suggest that these binaries might be the result of a continuous stellar formation in the center but we do not exclude other delivery scenarios. If the former is in effect, then in some cases some binaries might merge during the main sequence evolution (Prodan et al. 2015; Stephan et al. 2016). In general, we neglect the effects of Newtonian precession, which causes the precession of the outer orbit, and thus does not yield a suppression of the EKL mechanism (Li et al. 2015). We have rerun the merged systems via EKL processes for the GC case (see Table 1) with Newtonian precession, and demonstrated that it indeed has a negligible effect on the EKL mechanism for the physical picture we considered (see the discussion in Section III.3).

The number of BHs, their mass distribution, and number density are poorly known in NSCs. Theoretically, a single-mass distribution of objects forms a power law density cusp around a massive object with $n(r)\propto r^{-1.75}$ (Bahcall & Wolf 1976), where $n(r)$ is the number density and $r$ is the distance from the MBH. For multi-mass distributions, lighter and heavier objects develop shallower ($\propto r^{-1.5}$) and steeper cusps (typically $\propto r^{-2}$ to $r^{-2.2}$, and $r^{-3}$ in extreme cases), respectively (Bahcall & Wolf 1977; Hopman & Alexander 2006b; Freitag et al. 2006; Keshet et al. 2009; Aharon & Perets 2016). Recent observations of the stellar distribution in the Milky Way NSC identify a cusp with $n(r)\propto r^{-1.25}$ (Gallego-Cano et al. 2017; Schödel et al. 2017) consistent with the profile after a Hubble time (Baumgardt et al. 2017). As BHs are heavier than typical stars, they are expected to relax into the steeper cusps. The relaxation time of BH populations is much shorter: $0.1$–$1\,$Gyr (O’Leary et al. 2009), although it can become much longer than that in the case of a shallow stellar density profile (Dosopoulou & Antonini 2017). In this paper, we assume that the BH number density follows a cusp with either $n(r)\propto r^{-2}$ or $r^{-3}$ in our two sets of calculations.The BH mass in the two cases is set arbitrarily to $10^{7}$ and $4\times 10^{6}\,\rm M_{\odot}$, respectively, and we refer to the two models as “Bahcall-Wolf-like” (BW) and “Galactic Center” (GC) examples. Note that here “Galactic Center” refers to the assumed MBH mass (Ghez et al. 2008) and the observed stellar distribution (see below). In reality, the cusp distribution varies with BH mass. Thus, we have also generated initial conditions for the BW case so that $\beta$ in the number density distribution , $n(r)\propto r^{-(3/2+\beta)}$, is calculated by $\beta(m)=m/4M_{0}$, where $m$ is the binary mass and $M_{0}$ is the weight average mass (Keshet et al. 2009; Alexander & Hopman 2009; Aharon & Perets 2016, e.g.). We found that, due to the stability conditions, the initial condition distribution does not change significantly. For the GC case, the initial conditions distribution will change more significantly if we allow $\beta$ to vary with mass. However, we keep $\beta~=~3$ to investigate the effects of a steep number density distribution on the rates. As we will show later in the paper, the choice of number density distribution will have a very limited effect on the merger rate.

Note that $dN=4\pi r^{2}n(r)dr=4\pi r^{3}n(r)d(\ln r)$, where $N$ is the number of objects, and thus we choose to have the initial outer binary semi-major axis follow a uniform distribution in $a_{2}$ and $\ln a_{2}$ in the BW and GC model, respectively. We set the minimum $a_{2}$ to be one at which the relaxation timescale equals the outer binary gravitational wave merger timescale. The maximum $a_{2}$ is chosen to be $0.1$ pc, which corresponds to the value at which the eccentric Kozai-Lidov timescale is equal to the timescale on which accumulated fly-bys from single stars tend to unbind the binary (see Eq. 3 and 5). However, the BH-binary semi-major axis distribution changes after applying the stability criteria (see below).

Motivated by the recent LIGO detections, in both examples the mass of each of the BHs is chosen from a distribution uniform in logspace between $6-100$ M⊙ (i.e. $dN/dm\propto m^{-1}$, Miller 2002; Will 2004). For comparison, we have also run additional Monte-Carlo simulations in the case of the BW distribution, keeping all other parameters the same but using a mass distribution uniform in logspace between $5-15$ M⊙ instead of $6-100$ ⊙ (Belczynski et al. 2004). This had a negligible effect on the EKL induced mergers (see table 1). The binary separation $a_{1}$ is drawn from a uniform in log distribution between $0.1-50$ AU. This is consistent with Sana et al. 2012 distribution which favors short period binaries, although those BHs have already gone through stellar evolution and thus their exact distribution is unknown. We note that our stability criteria largely modifies the chosen distribution and yields a steeper distribution with a cut-off for systems beyond $50$ AU, as shown in Figure 1 (see below). The lower limit is such as to avoid mass-transfer before the supernova took place. We note that the vast majority of our binaries are soft, with only $2\%$ ($5\%$) of the GC (BW) binaries being hard binaries (Quinlan 1996, e.g.,).

The eccentricity of the BH binaries is chosen from a uniform distribution (Raghavan et al. 2010) and taking the outer eccentricity distribution to be thermal (Jeans 1919). Furthermore, the inner and outer argument of periapsis $\omega_{1}$ and $\omega_{2}$ are chosen from a uniform distribution between $0$ to $2\pi$. In addition, the mutual inclination $i$ is drawn from an isotropic distribution (uniform in $\cos i$).

After drawing these initial conditions we require that the systems satisfy dynamical stability, such that the hierarchical secular treatment is justified. We use two stability criteria. First we set

which is a measure of the relative strengths of the octupole and quadrupole level of approximations (Naoz 2016). Secondly, we require that the inner binary does not cross the Roche limit of the central MBH:

(Naoz & Silk 2014, e.g.,) See Antonini & Perets 2012 for quasi-secular evolution.. These stability criteria may significantly alter the distribution of BH binary systems that can survive to long timescales around the MBH (Stephan et al. 2016, similar to ). We show the before and after-stability distributions in Figure 1 for both examples. We use these distributions to initialize our runs. We note that the distribution of the angles ($\omega_{1},\Omega_{2}$ and $i$) remained the same after applying the stability criteria.

The main difference between the BW and GC models is the stable outer binary semi-major axis ($a_{2}$) distribution. As depicted in Figure 1, after the stability criteria $\sim 80\%$ of the BW systems remained stable, and their orbital parameter distribution did not significantly change. On the other hand, the choice of uniform in $\rm ln(a_{2})$ initial distribution for the GC case results in a steeper distribution than the nominal Bahcall & Wolf 1976 density profile, with $\sim 57\%$ of the systems remaining stable. Coincidentally, this type of steeper distribution is consistent with the de-projected density profile of the disk of massive stars (Bartko et al. 2009; Lu et al. 2009). Moreover, this distribution represents a strongly mass segregated cluster (Keshet et al. 2009).

Gravitational perturbations from the MBH can lead to eccentricity excitation of the binary BH which may lead to mergers, via the eccentric Kozai-Lidov (EKL) mechanism (Naoz 2016, see for review). We integrate the EKL equations (Naoz 2016, e.g.,) including GR precession and GW emission for 1000 (1500) systems in the GC (BW) case either until they merge or until they become unbound, whichever happens first. The latter takes place on the order of the typical timescale at which close encounters with other stars in the cluster cause the binary to unbind, (Binney & Tremaine 1987, e.g.,). This evaporation timescale has the form:

(see, Genzel et al. 2010; Tremaine et al. 2002, for the two cases) is the density of the surrounding stars, $m_{1}$ and $m_{2}$ are the masses of the two stellar mass BHs, $m_{3}$ is the average mass of the background stars, $G$ is the universal gravitational constant, $\ln\Lambda$ is the Coulomb logarithm, and $\sigma$ is the velocity dispersion. We adopted $\ln\Lambda=15$, $\sigma=280$ km s$\rm{}^{-1}\sqrt{0.1~pc/a_{2}}$, and $m_{3}=1~\rm M_{\odot}$ (Kocsis & Tremaine 2011). For the BW case, $\sigma_{0}=200$ km s-1 and $M_{0}=3\times 10^{8}~\rm M_{\odot}$ are constants. Note that the density distributions of the background stars are flatter than the density distribution of the BHs.

In the GC case, the average evaporation timescale is about $200$ Myr, but there is a very broad distribution from $\sim 1$ Myr to $\sim 3$ Gyr. In the BW case, the average evaporation timescale is about $120$ Myr, but again there is a broad distribution from $\sim 1$ Myr to $\sim 2$ Gyr. Binaries that do not merge are all evaporated by 10 Gyr, consistent with Stephan et al. 2016.

The above evaporation time assumes that the unbinding is taking place through interaction with the surrounding background stars. However, interaction with background stellar mass black holes may result in different evaporation timescales. If the stellar mass BHs follow the stellar density profile in Equation 4, then the evaporation timescale is similar to the evaporation timescale due to interactions with stars, up to a numerical factor that comes from the BH mass. However, if the density profile of the stellar BHs is steeper, like the one that might be expected from mass segregation over long timescales (Keshet et al. 2009, e.g.), the evaporation timescale might be shorter. Thus, we have also calculated two representative evaporation timescales resulting from steeper profile for $10$ $M_{\odot}$ and $30$ $M_{\odot}$ black holes, using densities from Aharon & Perets 2016 for these background black holes (see table 1). These evaporation timescales are much shorter on average than the evaporation timescales resulting from the blackground stars, and thus the absolute number of mergers is reduced. However, the merger rate is not reduced, as explained in Section III.3.

III. Results: EKL-induced mergers and GW-only mergers

The EKL mechanism has been shown to play an important role in producing short period binaries and merged systems (Thompson 2011; Naoz et al. 2011; Antonini & Perets 2012; Prodan et al. 2013; Naoz & Fabrycky 2014; Prodan et al. 2015; Antonini & Rasio 2016; Stephan et al. 2016; Naoz et al. 2016; Naoz 2016, e.g.,). The high eccentricity values achieved during the binary evolution leads to a shorter GW emission timescale, which may cause a merger before the binary become unbound as shown in Figure 2. If a merger takes place due to high eccentricity excitation we denote this merger as an EKL-induced merger.

The eccentricity excitation due to EKL can be inhibited by GR precession if the timescale for the latter is much shorter than the former (Naoz et al. 2013, e.g.,). The EKL timescale at the quadrupole level of approximation is estimated as:

(Antognini 2015, e.g.,), where $P$ denotes the orbital period and the GR precession of the inner orbit timescale is:

(Naoz et al. 2013, e.g.,) where $c$ is the speed of light. If $t_{\rm GR,inner}<t_{\rm quad}$, then EKL effects are negligible, the binary evolves due to close encounters with other stars in the cluster, and the inspiral is caused only by GW emission. The binary may approach merger if the GW timescale, $t_{\rm GW}$, (Peters 1964) is shorter than the timescale it takes the binary to become unbound (see Equation (3)). We identify those as GW-only mergers.

In other words we identify two channels for mergers. In the GW-only merger we have $t_{\rm GR}<t_{\rm quad}$ and $t_{\rm ev}>t_{\rm GW}$, and in the EKL-induced merger $t_{\rm quad}<t_{\rm GR}$, in which the eccentricity can be excited to near unity.

The EKL-induced high eccentricity excitations usually appear in a distinctive regime in the $\epsilon-i$ parameter space. This can be seen in Figure 3, where the EKL-induced mergers inhabit a specific regime of higher $\epsilon$ (between $10^{-4}$ and $10^{-2}$) and relatively large inclinations Although we note, that some EKL-induced mergers take place beyond the nominal Kozai angles, (i.e., $i<40^{\circ}$ and $i>140^{\circ}$), which is consistent with the near co-planar behavior (Li et al. 2014).. The EKL-induced mergers represent $\sim 7\%$ ($1.7\%$) of all GC (BW) Monte-Carlo systems. The systems that merge via the GW-only process (where no significant eccentricity excitations took place) occupy this parameter space uniformly. They represent about $9\%$ ($9.1\%$) of all the GC (BW) systems. The two sub-panels in Figure 3 show that the EKL yields a systematically shorter merger timescale. In other words, out of all systems that merged, $44\%$ are EKL-induced mergers in the GC case and about $16\%$ are EKL-induced mergers for the BW density profile.

As implied from Figure 3 these two merger channels can yield different predictions for the statistics of BH binary mergers. The EKL-induced mergers take place preferentially in systems that are closer to the MBH (see left panel of Figure 4), whereas the GW-only systems have no apparent trend. We note that the cut-off in GW-only mergers for $a_{2}<200$ AU takes place because systems below this value either have GR precession timescale much longer than the EKL, or unbind before the binary can merge.

The EKL-induced mergers systematically happen on shorter timescales compared to the GW-only mergers (see inset on Figure 5). On average EKL-induced mergers in the GC case take place within about $25$ Myrs while GW-only mergers merge on average after $296$ Myrs (see Figure 5 for the merger time histogram of the two channels). On the other hand, the BW case yields a shorter unbinding timescale and thus we get that the average EKL-induced mergers in the BW is about 1 Myr while the GW-only merger is about $166$ Myr.

As expected, in both of our examples, none of the EKL-induced mergers were initially hard. On the other hand, out of the GW-only mergers, $14\%$ in the GC example and $43\%$ of the BW example, were, in fact, hard binaries initially. Hardening channels (Quinlan 1996, e.g.,) may result in a shorter merger time for those binaries.

III.2. Rate Estimate

We also estimate the merger rate of BH-BH systems. The total merger rate is dominated by the merger rate of the soft binaries because the hard binaries represent a very small percentage of the total number of systems in our simulation ($2\%$ for the GC case and $5\%$ for the BW case). Thus, we only perform the following calculations for the soft binaries.The merger rate per unit volume is defined as:

where $n_{g}$ is the density of galaxies, $\Gamma$ is the merger rate per galaxy, and $f_{\rm MBH}$ is the fraction of galaxies containing a MBH. We adopt $n_{g}=0.02\,\rm Mpc^{-3}$ (Conselice et al. 2005) and $f_{\rm MBH}=0.5$ following (Antonini et al. 2015b; Antonini et al. 2015a). Note that this is a rather conservative number and we expect this fraction to be higher.

The merger time for our systems span a very broad range of timescales with a slight systematic trend for shorter merger time closer to the MBH (i.e., smaller $a_{2}$) as can be seen in Figure 6. However, these short-lived systems represent a small fraction of the total number of merged systems, and the apparent trend may be a result of statistical bias. To ascertain whether the merger rate has a dependency on $a_{2}$, we perform a bootstrap resampling of the merger time versus the $a_{2}$ distribution. We find that a bootstrapped analysis (both binned and not binned in $a_{2}$) yields a dependency in $a_{2}$ (see Figure 6). However, the percentage of systems close to the SMBH is small. Thus, we proceed by calculating the cumulative fraction of mergers as a function of time to calculate an average merger rate.

Figure 7 shows the cumulative fraction of merged systems as a function of time for various cases. We find that the cumulative number of merged systems fits a power law of the form:

where $a$ and $b$ are constants. We perform a bootstrap resampling on the cumulative number of mergers versus merger time distribution and generated 500 resampled data sets for the GC and BW cases. We fit the above power law to all of these data sets, which gives us a range of possible fits. The average best fit values for $a$ and $b$ are 2.8 and 0.19 for the GC case, and 2.3 and 0.19 for the BW case. We note that we also fit a combination of different exponential functions to the cumulative number of mergers (not shown), which perhaps follows a exponential decay-like process. These fits were systematically less good, but nonetheless yield consistent results with those shown below. We can define a ”half-life”, $t_{1/2}$, for each sample as the time it takes for half of the systems that will merge to merge. Thus, the merger frequency at $t=t_{1/2}$ is:

We use this frequency as the characteristic merger frequency for each data set.

Assuming a steady state number of BH binaries in the galactic nucleus, $N_{\rm steady}$, we can then estimate the merger rate per galaxy $\Gamma$ as:

where $f_{\rm merge}$ is the merger fraction from our simulation, equal to 0.15 for the GC case and 0.07 for the BW case. $N_{\rm steady}$ is highly uncertain, so we set it to be a free parameter between 1 and 500. The right panel of Figure 7 shows the total merger rate as a function of $N_{\rm steady}$. We use a nominal value of $N_{\rm steady}~=~200$ to calculate the merger rates for the GC and BW cases. We find a nominal average merger rate of $\sim$ 2 $\,\rm Gpc^{-3}\rm yr^{-1}$ for both the GC and BW cases. The average merger rate values are represented by the red and black solid lines in the right panel of Figure 7; the merger rate range obtained from the bootstrap is depicted by the shaded areas.

Our estimated merger rates for the GC and BW cases are on the same order as the merger rates estimated for globular clusters, $5$ Gpc-3 yr-1 (Rodriguez et al. 2016a), isolated triples in galaxies 0.14–$6\,\rm Gpc^{-3}\rm yr^{-1}$ (Silsbee & Tremaine 2016), and mergers following close encounters of initially unbound BHs 0.04–$3\,\rm Gpc^{-3}\rm yr^{-1}$ (O’Leary et al. 2009) Here we use the rates per single galaxy in O’Leary et al. 2009 Table 1 and multiply by $\xi=10$ and $n_{g}=0.02\,\rm Mpc^{-3}$. However, that mechanism is weakly sensitive to the NSC and MBH mass. Dwarf galaxies dominate the rates with a much higher $n_{g}$.. The current LIGO/VIRGO detection constrain the total merger rate of circular BH binaries to within 12–$240\,\rm Gpc^{-3}\rm yr^{-1}$ (LIGO Scientific and Virgo Collaboration 2017).

Note that Antonini & Perets 2012 rate estimation is based on the quadrupole based semi-analytical timescale which does not capture the system’s full dynamical behavior. The octupole level of approximation, used here, shortens the merger timescale for $16\%-40\%$ of the EKL-induced mergers, for the GC and BW distributions, respectively (see Figure 2). Furthermore, Antonini & Perets 2012 used the semi-analytical timescale given in Thompson 2011 to estimate the merger timescale. However, this does not capture the correct merger timescale for each binary BH as can be deduced from the inset in Figure 2.

III.3. Additional Tests

As stated in Section II, we have also considered the effects of using a narrower black hole mass distribution ranging from $5-15$ $M_{\odot}$ instead of $6-100$ $M_{\odot}$ for the BW case, while keeping all other parameters the same; and the effects of using the evaporation timescale resulting from interactions with $10$ $M_{\odot}$ and $30$ $M_{\odot}$ background black holes instead of background stars. For the $5-15$ $M_{\odot}$ black hole mass distribution, the average evaporation time is shorter than the nominal GC and BW cases because of the reduced BH binary mass (see Equation 3 and Table 2). We performed 500 Monte Carlo simulations and found that the number of GW-only mergers are greatly reduced, due to the fact that they generally take longer to merge. However, the short-timescale EKL mergers are unaffected, and we find that using a $5-15$ $M_{\odot}$ black hole mass distribution results in a merger population that predominantly consists of short-timescale EKL mergers. We also run a set quadrupole-only simulation for the EKL mergers in this case, and found that about $\sim 25\%$ take longer to merge than with the octupole, or do not merge at all (the total percentage of non-mergers with only quadrupole are $\sim 13\%$). Fitting Equation (8) to the number of mergers as a function of time results in a shorter merger half-life (this can be seen from the left panel of Figure 7), and therefore a higher merger frequency as defined in Equation (9). Thus, performing the calculation for $\Gamma_{\rm tot}$ results in a merger rate that is higher than before, as can be seen in the right panel of Figure 7. This makes sense, considering that a majority of the mergers in this case are EKL mergers that take place very quickly. However, a shorter average evaporation timescale and average merger timescale for a fixed steady-state number of binaries also means that we require a quicker BH binary replenishment mechanism (on a timescale roughly equal to the average evaporation and average merger timescale) to maintain steady state. Otherwise the high rate of merger seen in Figure 7 cannot be maintained. For this particular case, a replenishment timescale of about $\rm 10~Myr$ is required to maintain steady state.

We have a similar consequence resulting from using the evaporation timescales from interactions with background black holes of masses $10$ $M_{\odot}$ and $30$ $M_{\odot}$ rather than background stars. To do this, we first calculate the new evaporation timescales using the number densities for $10$ $M_{\odot}$ and $30$ $M_{\odot}$ found in Aharon & Perets 2016, then we compare the merger time of each of our soft binaries with the new evaporation timescales, and if the merger time is longer, we count it as a non-merger. We then fit the new merger distributions as a function of time to Equation (8) as done previously, and then calculate $\Gamma_{\rm tot}$ as a function of $N_{\rm steady}$, shown in the right panel of Figure 7. We find that using the evaporation timescale dominated by the $10$ $M_{\odot}$ black holes results in a merger rate very similar to the nominal BW and GC cases, but using the evaporation timescale dominated by the $30$ $M_{\odot}$ black holes results in a much higher merger rate than before (see right panel of Figure 7). In the $30$ $M_{\odot}$ background black holes case, the average evaporation timescale is very short compared to the nominal GC and BW cases (see Table 2). Similar to what happened with the $5-15$ $M_{\odot}$ black hole mass distribution, only BH binaries with short merger timescales can merge before they are evaporated, so the average merger timescale is very short, leading to a higher $\Gamma_{\rm tot}$. Furthermore, like in the $5-15$ $M_{\odot}$ black hole mass distribution case, this higher $\Gamma_{tot}$ will require a shorter binary replenishment timescale (roughly one Myr) to continue in steady state with a fixed number of BHs. For the $10$ $M_{\odot}$ background black holes case, the average evaporation timescale is very similar to the the average evaporation timescale in the nominal GC and BW cases (see Table 2). Thus, the average merger timescale is very similar to before, leading to a similar $\Gamma_{\rm tot}$. Note that since the evaporation timescale has a weak dependence on perturber mass, the significantly shorter evaporation timescale resulting from the $30$ $M_{\odot}$ is primarily consequence of their steeper density profile as compared to the $10$ $M_{\odot}$ background black holes.

Finally, we have also re-run our EKL mergers run in the GC case with Newtonian precession added. First, one should differentiate between Newtonian precession that occurs due to a spherical mass distribution (which was what we assumed in the paper), versus precession due to deviations from spherical symmetry as explored in Petrovich & Antonini 2017. Thus, in the following analysis we constrain ourselves to the spherically symmetric Newtonian precession. The Newtonian precession in our case causes the outer orbit to precess. While it may take place on similar timescale as the octupole timescale, it will not affect the inner orbit precession, which has a much greater effect on the EKL mechanism. However, to quantify how much Newtonian precession will affect our systems, we re-run the 63 EKL mergers for the GC case with Newtonian precession included, following Equation (44) in Tremaine (2005). We found that over 90% of these runs remained mergers even with Newtonian precession included. Furthermore, the merger time distribution remains the same.

IV. Discussion

We investigated the secular evolution of stellar-mass BH binaries in the neighborhood of the MBH in galactic nuclei using Monte Carlo simulations. Our equations included the hierarchical secular effects up to the octupole-level of approximation (Naoz 2016, the so-called EKL mechanism,), and general relativity precession of the inner and outer orbits, and gravitational wave emission between the two stellar BHs. During their evolution, binary stellar mass BHs may undergo large eccentricity excitations which can drive them to merge (see for example Figure 2). While BHs are expected to segregate toward the center of the galactic nuclei, the power law index of the number density of the cusp is still uncertain. As a proof-of-concept, we explored two cases: $r^{-2}$, denoted as BW; and $r^{-3}$, denoted as GC. We find a consistent merger rate in both cases.

We identified two channels for BH-BH mergers. In the first type of merger, the general relativity precession timescale is shorter than the quadrupole timescale, suppressing eccentricity excitations. However, the initial GW merger timescale is shorter than the evaporation timescale, leading to a merger. We call these GW-only mergers. These systems will merge in the absence of an MBH in the center of a galaxy. For the second type of merger, the EKL mechanism produces large eccentricity oscillations in the inner binary orbit, driving the binary BHs to merge. We call these EKL-induced mergers.

GW-only and EKL-induced mergers occupy different parts of the parameter space for both number density profiles we examined. Specifically, considering the $\epsilon-i$ parameter space, where $\epsilon$ is given by Eq. (1). EKL-induced mergers will preferentially occupy a regime of high $\epsilon$ and close to $90^{\circ}$ mutual inclination. On the other hand, GW-only mergers are uniformly spread in $i$ and $\rm ln(\epsilon)$ (see Figure 3). This yields a prediction for the merger distribution as a function of distance from the MBH: the EKL-induced mergers will be systematically closer to the MBH. This is depicted in Figure 4.

Finally we estimate the total merger rate to be at the order of unity $\,\rm Gpc^{-3}\rm yr^{-1}$ and perhaps even higher. This rate is coincidentally comparable to the estimated merger rate in globular clusters, which suggests that galactic nuclei may host a significant fraction of the BH-BH mergers. Note that Antonini & Rasio 2016 showed that if the natal kick of BH binaries is higher than $50$ km sec-1, then the efficiency of forming BH binaries in globular clusters is suppressed significantly (Chatterjee et al. 2017) And also above this natal kick the mergers of isolated binaries is significantly suppressed (Belczynski et al. 2016b) , and BH mergers in nuclear star clusters dominate over mergers in globular clusters. We note that, if natal kicks are smaller, so that BH binaries form efficiently in globular clusters, BH mergers in nuclear star cluster with an MBH is still comparable to that of globular clusters.

We note that we have neglected the effects of resonant relaxation, as they should be added self-consistently (Rauch & Tremaine 1996; Sridhar & Touma 2016). These effects operate on longer timescales than the quadrupole timescale. However, they can refill the high inclination EKL regime in some cases, causing further eccentricity excitation. This is beyond the scope of this paper. Recently Petrovich & Antonini 2017 showed that variations from spherical symmetry in the potential may induce eccentricity excitation at large $a_{2}$. Thus, their proposed mechanism may even increase the merger rate discussed here for systems at larger distances. However, their study also neglected vector resonant relaxation which operates on the same timescales and distances. Therefore, further investigation is needed for a more accurate treatment of these effects.

Our results suggest that the key parameter controlling the merger rate is the steady-state number of BH-BH binaries. We find that, for typical BH–BH binary populations often assumed in the literature, the predicted merger rate is $\sim 1-3$$\,\rm Gpc^{-3}\rm yr^{-1}$; but this can become much higher, depending on the exact assumptions made (see Figure 7). Our steady-state assumption requires a binary replenishment mechanism, without which the EKL-induced BH-BH merger rate would rapidly decay (as implied from Figure 5). The BH-BH binary population in real systems may be replenished from the outside, for example from globular clusters that spiral in, or from nearby star formation (Hopman & Alexander 2006a; Hopman & Alexander 2006b; Gnedin et al. 2014; Antonini et al. 2015a; Aharon & Perets 2016, e.g.,). A different channel for recent in-situ formation of BH–BH binaries may exist near the centers of E+A, or post-starburst, galaxies (Dressler & Gunn 1983, e.g.,). These galaxies are special because it appears that they underwent a star formation episode that terminated abruptly $\sim 1$ Gyr ago. Thus, this population may hold recent BH-BH binaries near their nuclei. While these galaxies are a relatively rare subtype of elliptical galaxies, they may increase the stellar mass of the galaxy by $\sim 10\%$ (Swinbank et al. 2012, e.g.,). Thus, due to their recent star formation and overdensity in the nuclei, these galaxies were linked to an enhancement of tidal disruption events (Arcavi et al. 2014; Stone & van Velzen 2016, e.g.,). Therefore, these galaxies may have an overabundance of BH binaries, which can cause enhancement of BH-BH mergers. Recently, Bartos et al. 2017 showed that with a big enough sample of observed BH-BH mergers, a certain BH-BH merger channel can be statistically correlated with rare galaxy types, including E+A galaxies. This may prove a valuable test for our model in the future.

The BH-BH merger scenario presented here, in our proof-of-concept examples, suggests that MBH gravitational perturbations can enhance the merger rate. Since these mergers take place close to the SMBH, the acceleration of the binary around the SMBH may cause a detectable Doppler shift of the waveform (Meiron et al. 2017; Inayoshi et al. 2017). Furthermore, a GW echo of the binary BH merger caused by the SMBH may be detectable (Kocsis 2013). We have provided the parameter space at which these mergers take place, and demonstrated that they occupy a different part of the parameter space. This suggests that this channel may be statistically distinguished with a sufficiently large sample of mergers (Hoang et al. in prep.).

Finally, motivated by the calculation of the rate of BH-BH mergers for $>0.1$ pc by Petrovich & Antonini 2017, we follow their set of assumptions of a compact object formation rate between $2\times 10^{-5}-10^{-4}$ $\rm yr^{-1}$, with a fraction of surviving BH-BH binaries of $2.5-4.5\%$. The merger fraction from Petrovich & Antonini 2017 along with these assumptions yields a merger rate of $0.6-15$ $\,\rm Gpc^{-3}\rm yr^{-1}$. Similarly, with our merger fraction and these assumptions we find a merger rate of $0.7-14$ $\,\rm Gpc^{-3}\rm yr^{-1}$ for our $<0.1$ pc regime. Therefore, assuming that these assumptions are correct, we conclude that the BH-BH merger rate in the proximity of galactic nuclei can be as significant as $\sim 30$ $\,\rm Gpc^{-3}\rm yr^{-1}$.

References