Gravitational-wave Merging Events from the Dynamics of Stellar-mass Binary Black Holes around the Massive Black Hole in a Galactic Nucleus

Denote E = −GM_•/a₂ and $J=\sqrt{{{GM}}_{\bullet }{a}_{2}(1-{e}_{2}^{2})}$ as the energy and angular momentum of the outer orbit of a BBH. Both E and J are evolved under the NR. J of the BBH is additionally affected by RR. The NR is due to the weak gravitational interactions between a BBH and field stars, where we treat the BBH as a point-mass particle. The NR (or two-body relaxation) timescale is given by (Binney & Tremaine 1987)

$$\begin{eqnarray}&&{T}_{\mathrm{NR}}=\displaystyle \frac{0.34{\sigma }^{3}}{{\left({{Gm}}_{\star }\right)}^{2}{n}_{\star }\mathrm{ln}{\rm{\Lambda }}},\end{eqnarray} \tag{ 2 }$$

where Λ ≃ ${M}_{\bullet }$/m_⋆ (Bahcall & Wolf 1976). Here σ ∝ r^−1/2, and thus ${T}_{\mathrm{NR}}\propto {r}^{{\alpha }_{\star }-3/2}$. The time evolutions of the diffusion process in the E−J space due to two-body relaxation can be described by the Fokker–Planck equation and simulated by a Monte Carlo scheme according to their orbit-averaged diffusion coefficients, i.e., D_EE, D_E, D_JJ, D_J, and D_EJ (Lightman & Shapiro 1977). Here D_EE (D_JJ) and D_E (D_J) describe the orbit-averaged scatterings of the energy (angular momentum) and its drift, respectively. D_EJ describes the correlations between the scatterings of energy and angular momentum. The details of the formalisms are given in Appendix A.1. We adopt a Monte Carlo scheme similar to Shapiro & Marchant (1978) and Bar-Or & Alexander (2016) to calculate the two-body relaxation in E−J of BBHs due to the field stars. More details will be explained in Section 2.5.

In the inner part of the cluster (but not too close to the MBH, so the general relativistic precession is still not significant), RR becomes important, and the resonant torques between the orbits of the field stars can change quickly the orbital angular momentum of the BBHs. The RR changes both the amplitude (scalar RR) and the direction of the angular momentum (vector RR). The timescale of scalar RR is given by (Hopman & Alexander 2006)

$$\begin{eqnarray}&&{T}_{\mathrm{RR}}^{s}=\displaystyle \frac{{A}_{\mathrm{RR}}^{s}}{{t}_{\omega }}{\left(\displaystyle \frac{{M}_{\bullet }}{{m}_{\star }}\right)}^{2}\displaystyle \frac{{P}^{2}({a}_{2})}{{N}_{\star }(\lt {a}_{2})},\end{eqnarray} \tag{ 3 }$$

where ${A}_{\mathrm{RR}}^{s}$ ≃ 3.56 (Rauch & Tremaine 1996). Here t_ω = 2π/ν_p is the timescale of orbital precession, and ν_p is given by Equation (37). P is the orbital period, and N_⋆(<a₂) is the number of field stars within a distance of a₂. The scalar RR can be suppressed owing to the rapid general relativistic orbital precession below the locus called the "Schwarzschild barrier" (SB; Merritt et al. 2011; Antonini & Merritt 2013), where the orbital precession frequency equals the coherence frequency, i.e., ν_p = 2π/T_c(a₂) (Bar-Or & Alexander 2016). Here T_c(a₂) is the coherence timescale given by Equation (41).

The timescale of vector RR is given by (Hopman & Alexander 2006)

$$\begin{eqnarray}&&{T}_{\mathrm{RR}}^{v}\,=\,2{A}_{\mathrm{RR}}^{v}\left(\displaystyle \frac{{M}_{\bullet }}{{m}_{\star }}\right)\displaystyle \frac{P({a}_{2})}{{N}_{\star }^{1/2}(\lt {a}_{2})},\end{eqnarray} \tag{ 4 }$$

where ${A}_{\mathrm{RR}}^{v}\simeq 0.31$ (Rauch & Tremaine 1996). By a method similar to Bar-Or & Alexander (2016), we can consider the scalar RR by calculating the diffusion coefficients in angular momentum ${D}_{{JJ}}^{\mathrm{RR}}$ and ${D}_{J}^{\mathrm{RR}}$. For vector RR, we simply consider it as a diffusion process. The details of the scalar and vector RR are shown in Appendix A.2.

During the evolution, if the BBH approaches too close to the MBH, they will likely be disrupted by the tidal force of the MBH. The tidal radius of binary stars is given by

$$\begin{eqnarray}&&{r}_{t}={\left(\displaystyle \frac{3{M}_{\bullet }}{{m}_{\mathrm{BBH}}}\right)}^{1/3}{a}_{1},\end{eqnarray} \tag{ 5 }$$

where m_BBH = m_A + m_B is the total mass of the binary. The loss cone of the angular momentum is ${J}_{\mathrm{lc}}^{2}$ ≃ 2 ${M}_{\bullet }$r_t. The empty loss cone region requires that the change of angular momentum per orbital period of the BBH is smaller than the size of the lose cone, i.e., dJ < J_lc, and the full loss cone region requires that dJ > J_lc. In the full lose cone region, the binary can jump in and out of the lose cone multiple times before it encounters the MBH. In the empty loss cone region, the binary takes multiple periods to move into the lose cone and may have encountered the MBH multiple times. As the probability of tidal disruption is a function of distance from the MBH (e.g., in the case of binary stars; Zhang et al. 2010), the BBH may experience multiple encounters with the MBH before the tidal disruption. There are also chances that the BBH moves away from the lose cone after the multiple encounters with the MBH. Zhang et al. (2010) have showed that the multiple encounters can change both the SMA and the eccentricity of the binary, which thus may lead to merging events.

To consider accurately such effects, we use three-body numerical integrations to calculate the outcome of each encounter between the BBH and the MBH. Considering that the MBH can impose strong impacts on the inner orbit of the BBH even before it reaches the loss cone, we set up the three-body integration when a₂(1−e₂) < 3r_t. At the beginning of these explicit three-body simulations, initially the binary is placed at 10⁴a₁ (or at apocenter of the outer orbit if a₂ is smaller). After the encounter, the change of the inner orbit of the binary at distance 10⁴a₁ from the MBH is recorded for the next encounters. Such a process repeats until either the binary is disrupted or merged owing to the GW emission, or moving out of the loss cone.

The BBH–MBH system forms a natural triple system, and the tidal force of the MBH can trigger KL oscillation on the inner and outer orbits of the BBH. The KL effect will increase the eccentricity of the inner orbit of the binary under some preferred orbital configurations (e.g., Naoz et al. 2013; Naoz 2016). The period of KL oscillation is given by (Kiseleva et al. 1998)

$$\begin{eqnarray}\begin{array}{rcl}{T}_{K} & = & \displaystyle \frac{2{P}_{2}^{2}}{3\pi {P}_{1}}{\left(1-{e}_{2}^{2}\right)}^{3/2}\displaystyle \frac{{m}_{A}+{m}_{B}+{M}_{\bullet }}{{M}_{\bullet }}\\ & \simeq & \ \displaystyle \frac{2}{{3}^{1/3}\pi }{\left(\displaystyle \frac{{r}_{p}}{{r}_{t}}\right)}^{3/2}{P}_{2}.\end{array}\end{eqnarray} \tag{ 6 }$$

The GR procession of the inner orbit of the binary is given by

$$\begin{eqnarray}&&{T}_{\mathrm{GR}}=\displaystyle \frac{2\pi {a}_{1}^{5/2}{c}^{2}(1-{e}_{1}^{2})}{3{G}^{3/2}{\left({m}_{A}+{m}_{B}\right)}^{3/2}}.\end{eqnarray} \tag{ 7 }$$

If T_GR < T_K, or

$$\begin{eqnarray}&&{r}_{p}\lesssim 0.58{\left(\displaystyle \frac{{T}_{\mathrm{GR}}}{{P}_{2}}\right)}^{2/3}{r}_{t},\end{eqnarray} \tag{ 8 }$$

the KL effect will be suppressed (e.g., Ford et al. 2000; Naoz et al. 2013). In most cases, the KL effect is only significant when the pericenter of the outer orbit of the BBH is close to the tidal radius of the MBH. In this work, we consider the KL effect only when r_p falls in the range of 3r_t to min[20, 0.58(T_GR/P₂)^2/3]r_t. When r_p < 3r_t we use the explicit three-body integration to calculate the interactions between the BBHs and the MBH (See Section 2.2).

When the eccentricity of the inner orbit, e₁, is excited to a high value owing to the KL effect (or other effects), the GW radiation becomes important and the BBH is likely to be merged in a very short timescale. The GW orbital decay timescale is given by (Peters 1964)

$$\begin{eqnarray}&&{T}_{\mathrm{GW}}=\displaystyle \frac{5}{64}\displaystyle \frac{{c}^{5}{a}^{4}}{{m}_{A}{m}_{B}({m}_{A}+{m}_{B}){G}^{3}}\displaystyle \frac{{\left(1-{e}_{1}^{2}\right)}^{7/2}}{1+\tfrac{73}{24}{e}_{1}^{2}+\tfrac{37}{96}{e}_{1}^{4}}.\end{eqnarray} \tag{ 9 }$$

The evolution of the SMA and eccentricity of the inner binary by GW radiation is given by (Peters 1964)

$$\begin{eqnarray}\begin{array}{rcl}{\dot{a}}_{1}^{\mathrm{GW}} & = & -\displaystyle \frac{64}{5}\displaystyle \frac{{m}_{A}{m}_{B}({m}_{A}+{m}_{B}){G}^{3}}{{c}^{5}{a}_{1}^{3}}\displaystyle \frac{1+\tfrac{73}{24}{e}_{1}^{2}+\tfrac{37}{96}{e}_{1}^{4}}{{\left(1-{e}_{1}^{2}\right)}^{7/2}},\\ {\dot{e}}_{1}^{\mathrm{GW}} & = & -\displaystyle \frac{304}{15}{e}_{1}\displaystyle \frac{{m}_{A}{m}_{B}({m}_{A}+{m}_{B}){G}^{3}}{{c}^{5}{a}_{1}^{4}}\displaystyle \frac{1+\tfrac{121}{304}{e}_{1}^{2}}{{\left(1-{e}_{1}^{2}\right)}^{5/2}}.\end{array}\end{eqnarray} \tag{ 10 }$$

Here the dot means the derivative with respect to time t. When the eccentricity is high, the argument of pericenter (denoted as g₁) of the BBH's inner orbit can precess significantly. For simplicity, we consider only the post-Newtonian GR effect, i.e.,

$$\begin{eqnarray}&&{\dot{g}}_{1}^{\mathrm{PN}}=\displaystyle \frac{3{\left({m}_{A}+{m}_{B}\right)}^{3/2}{G}^{3/2}}{{c}^{2}{a}_{1}^{5/2}(1-{e}_{1}^{2})}.\end{eqnarray} \tag{ 11 }$$

We consider the evolution of the inner and outer orbits of the BBHs under the above two effects and the KL oscillations. For the KL effect, we adopt the time evolution formalism in Naoz et al. (2013). The equations of motion in a₁, e₁, and g₁ now become

$$\begin{eqnarray}&&{\dot{a}}_{1}={\dot{a}}_{1}^{\mathrm{GW}},\ {\dot{e}}_{1}={\dot{e}}_{1}^{\mathrm{GW}}+{\dot{e}}_{1}^{\mathrm{KL}},\ {\dot{g}}_{1}={\dot{g}}_{1}^{\mathrm{PN}}+{\dot{g}}_{1}^{\mathrm{KL}}.\end{eqnarray} \tag{ 12 }$$

Here ${\dot{e}}_{1}^{\mathrm{KL}}$, ${\dot{g}}_{1}^{\mathrm{KL}}$ can be found in Naoz et al. (2013). For the KL evolution of the other orbital elements of BBHs, see Naoz et al. (2013).⁴ Recent studies show that these analytical formulae capture only the secular terms of the evolution, while the nonsecular terms may be important for a BBH–MBH triple system (Grishin et al. 2018). The nonsecular terms can enhance the maximum eccentricity in oscillations and increase the merging rates by up to an order of magnitude (Fragione et al. 2018). Thus, our merging rates due to the KL effect may have been underestimated.

It would be interesting to investigate the eccentricity and the SMA when the peak of GW frequency enters the observing band of ground-based GW detectors, i.e., f_GW ≳ 1 Hz. We describe the method in Appendix B, which is mainly based on Equations (43) and (44). Figure 2 illustrates the two possible cases of merging a BBH. The first one is shown by the top panels, i.e., when the GW frequency approaches f_GW = 1 Hz, the GW dominates the evolution of the BBH and the eccentricity of the BBH reduces to e₁ ∼ 0.1. Alternatively, as the bottom panels indicate, f_GW can be larger than 1 Hz before the GW dominates, and the eccentricity of the BBH is extremely high, i.e., 1−e₁ ∼ 10⁻³–10⁻⁴. We find that the latter case is very rare, i.e., <0.1% of all the BBH mergers in our simulation. However, if the entering frequency can be as low as about ∼0.1 Hz, then the probability can be much higher; the in-spiraling phase of BBHs, before the GW dominates the evolution, can now be observable (see the top right panel of Figure 2). If such phenomena can be detected, probably in decihertz GW detectors (Mandel et al. 2018), it can be strong evidence that those GW mergers are due to KL oscillations.

Zoom In Zoom Out Reset image size

The BBH will experience encounters with the field stars (or stellar remnants) if the stellar density is high. Such encounters can change both a₁ and e₁, or even lead to exchange of the binary components. There is a critical velocity that defines the hard and soft regions of the binary encounter (Hut & Bahcall 1983),

$$\begin{eqnarray}&&{v}_{c}=\sqrt{\displaystyle \frac{{{Gm}}_{T}}{{a}_{1}}\displaystyle \frac{{\mu }_{12}}{{m}_{3}}},\end{eqnarray} \tag{ 13 }$$

where μ₁₂ = m_Am_B/(m_A + m_B) is the reduced mass of the binary, m_T = m_A + m_B + m₃ is the total mass, and m₃ is the mass of the incoming star. Suppose that the velocity of the incoming star with respect to the mass center of BBHs is ${v}_{\infty }$; if ${v}_{\infty }\gt {v}_{c}$, the binary is soft, and it is hard if ${v}_{\infty }\lt {v}_{c}$. After the encounter, a hard binary always becomes harder, while for a soft binary it could become harder or softer, or even ionized (disrupted).

The event rates of binary–single encounter depend on the density profiles, whether the binary is hard or soft, and other details of the cluster. The detailed calculations can be found in Appendix C. Assuming that α_⋆ = 1.75, for a soft binary, the event rate is given by

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{EC}} & = & 2.23\,\times \,{10}^{-7}\,{\mathrm{yr}}^{-1}\,{\rm{\Theta }}{\left(\displaystyle \frac{p}{0.1\mathrm{au}}\right)}^{2}{\left(\displaystyle \frac{{a}_{2}}{{10}^{3}\mathrm{au}}\right)}^{-9/4}\\ & & \times \ {\left(\displaystyle \frac{{M}_{\bullet }}{4\times {10}^{6}{M}_{\odot }}\right)}^{7/8}\left(\displaystyle \frac{10\,{M}_{\odot }}{{m}_{\star }}\right),\end{array}\end{eqnarray} \tag{ 14 }$$

and for a hard binary, we have

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{EC}} & = & 1.88\,\times \,{10}^{-9}\,{\mathrm{yr}}^{-1}\,{\rm{\Phi }}\left(\displaystyle \frac{p}{0.1\,\mathrm{au}}\right){\left(\displaystyle \frac{{a}_{2}}{{10}^{4}\mathrm{au}}\right)}^{-5/4}\\ & & \times \ {\left(\displaystyle \frac{{M}_{\bullet }}{4\times {10}^{6}{M}_{\odot }}\right)}^{-1/8}\left(\displaystyle \frac{10\,{M}_{\odot }}{{m}_{\star }}\right)\left(\displaystyle \frac{{m}_{T}}{30\,{M}_{\odot }}\right),\end{array}\end{eqnarray} \tag{ 15 }$$

where Θ and Φ are dimensionless constants of order unity, which depends slightly on the orbital eccentricity e₂ (see Appendix C for more details). Here p is the pericenter distance of the incoming star with respect to the BBH, and it relates to the impact parameter b by

$$\begin{eqnarray}&&b\,=\,p\sqrt{1\,+\,2{m}_{T}G/({v}_{\infty }^{2}p)}.\end{eqnarray} \tag{ 16 }$$

We can see that, for a typical galactic nucleus, the binary–single encounter should be considered in the simulation, as they could be frequent, especially in the inner regions of the nucleus cluster (see also Figure 1).

As the outcome of a binary–single encounter is complex and we find that the change of the inner orbit of the BBH may have significant impacts on their GW merging rates, such a process should be calculated as accurately as possible. Here we use explicit three-body numerical integration to calculate the outcome of the binary–single encounters. The incoming star is assumed to be coming from infinity with ${v}_{\infty }=\sqrt{{{GM}}_{\bullet }/(2{a}_{2})}$ and impact parameter 0 < b < b_max, which follows a distribution f(b) ∝ b. Here b_max relates to p_max according to Equation (16). We set ${p}_{\max }=\,\max [6{({m}_{\star }/{m}_{\mathrm{BBH}})}^{1/2},2]{a}_{1}(1+{e}_{1}/2)$ such that the simulation results can be converged. As m_⋆/m_BBH ∼ 1–0.1, p_max ∼ 2 a₁–9 a₁. The initial orbital orientation of the incoming star is set to be random. Initially the incoming star is put at a distance of 50 a₁ away from the BBH. After each encounter, we determine the outcomes according to the relative energies and separations between each particle (black holes or incoming stars).

Figure 3 illustrates the four most likely outcomes of binary–single encounters in our simulation: (1) Exchange of one of the binary components with the incoming object. This event is common if m_BBH ≲ m_⋆, i.e., when the incoming object is massive. The incoming object is most likely a star considering that the number fraction of black holes is very small (≲10⁻²; see Hopman & Alexander 2006); thus, the exchange event usually leads to a non-BBH object. In the current work, we simply remove such binaries in the simulation. However, we notice that in the future our work can be easily expanded to further consider the evolution of binaries of which one component is a star or a neutron star and the other is a black hole. (2) Ionization, which means that the binary is destroyed owing to the encounter. This event is only possible if ${v}_{\infty }\gt {v}_{c}$, i.e., when the binary is soft. (3) Flyby, where the binary maintains integrity after the encounter. However, the energy and angular momentum of the inner binary, or equivalently a₁ and e₁, are both changed. (4) Three bodies experience chaotic evolutions, and finally a binary is formed and the third star is ejected from the system. These events are only possible if v₃ < v_c, i.e., when the binary is hard. There are other possible outcomes, e.g., hierarchical resonants (Heggie & Hut 1993), however only for hard binaries. As most of the binaries in our work are soft binaries, we do not discuss them in more detail here.

Zoom In Zoom Out Reset image size

We use the Monte Carlo schemes that are similar to Shapiro & Marchant (1978) to simulate the evolution of the inner and outer orbits of a BBH around the MBH. Our method has combined all of the effects described above, i.e., the NR and RR of the BBHs (Section 2.1), the tidal disruption of the BBHs (Section 2.2), the KL oscillations and the GW radiation of the inner orbit of the BBHs (Section 2.3), and the binary–single encounters (Section 2.4). The details of the method are described as follows.

Define x = E/E₀ as the dimensionless energy, where E₀ =−GM_•/r_h is the characteristic energy. We assume that the energy distribution of the field stars is given by (Lightman & Shapiro 1977; Magorrian & Tremaine 1999) f(E) = ${(2\pi {\sigma }_{0}^{2})}^{-3/2}{n}_{0}{g}_{\star }(x)$, where

$$\begin{eqnarray}&&{g}_{\star }(x)=\displaystyle \frac{{\rm{\Gamma }}({\alpha }_{\star }+1)}{{\rm{\Gamma }}({\alpha }_{\star }-1/2)}{x}^{{\alpha }_{\star }-3/2},\,\,3\gt {\alpha }_{\star }\gt 1/2,\end{eqnarray} \tag{ 17 }$$

and x > 0, where α_⋆ is the power-law index of the density profile of the field stars. Here ${\sigma }_{h}^{2}={{GM}}_{\bullet }/{r}_{h}$, and n₀ is the number density of the stars at position r = r_h. If x < 0, we simply set g_⋆(x) = e^x (Lightman & Shapiro 1977).

We initially put all the BBHs near the edge of the influence radius of the MBH, i.e., r = r_i = r_h (or x = 0.5). The initialization of the inner orbits of the binaries depends on the problems concerned and will be introduced in more detail at the beginning of Sections 3 and 4. The dimensionless distribution function of the energy of the BBHs' outer orbit is defined by

$$\begin{eqnarray}&&{g}_{\mathrm{BBH}}(x)=\displaystyle \frac{n(x)}{n(1)}\displaystyle \frac{{\rm{\Gamma }}({\alpha }_{\star }+1)}{{\rm{\Gamma }}({\alpha }_{\star }-1/2)}{x}^{5/2},\end{eqnarray} \tag{ 18 }$$

where n(x) is the number distribution of BBHs. This is obtained according to n(E) ∝ f(E)E^−5/2 in an isotropic cluster, where n(E) is the number density of BBHs as a function of the outer orbits' energy (Hopman 2009). In the above equation, g_BBH is normalized for the convenience of model comparison, as g_BBH(1) = g_⋆(1) for all models with the same α_⋆.

If t is the current time of the simulation, we take the following steps for each BBH:

• 1.  
  A time step size δt is determined according to Section 2.6.
• 2.  
  We trace the NR of the energy (E = −GM_•/(2a₂)) and NR and RR of the angular momentum ($J=\sqrt{{{GM}}_{\bullet }{a}_{2}(1-{e}_{2}^{2})}$) of the outer binary by the method described in Section 2.1. The change of energy and angular momentum is given by

  $$\begin{eqnarray}\begin{array}{rcl}{dE} & = & {D}_{E}^{\mathrm{NR}}\delta t+{y}_{1}\sqrt{{D}_{{EE}}^{\mathrm{NR}}\delta t},\\ {dJ} & = & \left({D}_{J}^{\mathrm{NR}}+{D}_{J}^{\mathrm{RR}}\right)\delta t,\\ & & +\ {y}_{2}\sqrt{{D}_{{JJ}}^{\mathrm{NR}}\delta t}+{y}_{3}\sqrt{{D}_{{JJ}}^{\mathrm{RR}}\delta t}.\end{array}\end{eqnarray} \tag{ 19 }$$

  Here y₁, y₂, and y₃ are unity normal random variables. Note that y₁ and y₂ have correlations $\rho ={D}_{{EJ}}/\sqrt{| {D}_{{EE}}{D}_{{JJ}}| }$. They can be obtained by generating a bi-normal distribution and then taking the transformation, ${y}_{2}\to {y}_{1}\rho +{y}_{2}\sqrt{1-{\rho }^{2}}$.
• 3.  
  If the pericenter of the outer orbit of a BBH is within r_p < 3r_t, and when int(t/P) − int[(t − δt)/P] > 0, i.e., this BBH is near the loss cone for an orbital period, we start a three-body simulation, which traces the interactions between this BBH and the MBH (see Section 2.2 for more details).
• 4.  
  The expected number of collisions between each BBH with the field stars is calculated according to k = R_ECδt, where R_EC is given by Equation (52). The realized number of collisions, n_EC, for each BBH is generated according to the Poisson distribution P_EC(k) (given by Equation (53)). We then take a number of n_EC successive three-body integrations to consider these encounter events (the details can be found in Section 2.4). Usually, the time step is set such that n_EC ∼ 1 (see Section 2.6). Note that any change of the inner orbit of the BBH is saved for the next steps. After each encounter, the BBH is considered to be ionized if $\epsilon \lt {10}^{-3}{m}_{\star }| E|$; here $\epsilon \,={m}_{\mathrm{BBH}}^{2}/(2{a}_{1})$.
• 5.  
  If the pericenter of the outer orbit of a BBH is $3{r}_{t}\lt {r}_{p}\lt \min [20,0.58{({T}_{\mathrm{GR}}/{P}_{2})}^{2/3}]{r}_{t}$, we consider the KL oscillation of the inner binary that lasts for a time period of δt (see Section 2.3 for more details). We calculate the GW radiation timescale T_GW for this BBH at any moment, and the BBH is considered to be merged owing to GW radiation if δt−t' < 50 T_GW, where 0 < t' < δt is the current time of simulation in the KL subroutines. We calculate the eccentricity and SMA of the BBH when it reaches f_GW = 10 Hz according to Appendix B. If $| {\dot{e}}_{1}^{\mathrm{KL}}| \gt | {\dot{e}}_{1}^{\mathrm{GW}}|$ and f_GW ≥ 10 Hz, the BBH is now entering the LIGO frequency band, and we record the current value of eccentricity and SMA of the BBH.
• 6.  
  If r_p > min[20, 0.58 (T_GR/P₂)^2/3]r_t, we calculate the GW radiation of the inner binary according to Equation (10). Similarly, we consider the BBH to be merged if 10 T_GW < δt, and we calculate the eccentricity and SMA of the BBH when it reaches f_GW = 10 Hz according to Appendix B.
• 7.  
  We remove any BBH that has a dimensionless energy of x < 0.1 or x > 10⁴. The inner boundary corresponds to a distance of about 10 au if ${M}_{\bullet }$ = 10⁶–10⁷ ${M}_{\odot }$.
• 8.  
  If the BBH is not destroyed (due to either merger of GW radiation or ionization, or tidal disruption, etc.), repeat steps 1 to 7 until the density profile of BBHs, i.e., g_BBH, reaches an equilibrium state.⁵ For any BBH that is destroyed, its simulation ends and the information is saved for statistics. The information includes the inner and outer orbits of the BBH, the current simulation time, the number of different dynamical events it has experienced, and so on.

We trace the evolution of the ascending node (Ω₁, Ω₂), inclination (I₁, I₂), and the argument of periapsis (g₁, g₂) for both the inner and outer orbits of the BBHs. Note that all these elements are defined with respect to an arbitrary selected reference plane and a direction within it, which are the same for all BBHs in a simulation. We assume that the vector RR process affects only the orientation of the outer orbit (Ω₂ and I₂), but not the inner orbit. During the binary–single encounter, the orbital elements of the inner orbit are changed and recorded for successive encounters. Currently, we do not consider the change of the outer orbit due to the binary–single encounter, but defer it to future studies. Such simplification should not significantly affect our results, as the velocity dispersion in the cluster is usually much larger than the velocity of the inner orbits of BBHs. All angles are considered and traced in the KL oscillations and the BBH–MBH encounters. Note that the additional precession of g₂ due to distributed mass has not been included in the current simulation (see Section 5 for how it affects our results).

As sometimes the density of the binaries drops rapidly near the MBH, to increase the number statistics, we use the "clone scheme" similar to Shapiro & Marchant (1978). We generate a number of Π − 1 clones of a BBH (either an original or a clone one) when its energy crosses the boundary x = 10ⁱ from lower x, where i = 1, 2, 3. Here the number Π is the amplification factor in the scheme and is usually a number selected between 2 and 10 in each run, such that the number of BBHs in both outer and inner regions of the cluster is sufficiently large for statistics. If a clone particle crosses such a boundary from higher x, it will be removed from the simulation (its information is not saved for statistics). When a BBH is destroyed because of any step in (3)–(7), if it is the original particle, a new particle (regarded also as a new original particle) at r = r_i is generated. But if it is a clone particle, it is removed from the simulation. Note that any BBH located at position 10ⁱ < x < 10ⁱ⁺¹ has a statistical weight of Π⁻ⁱ.

We use the code DORPI5 based on the explicit fifth (fourth)-order Runge–Kutta method (Dormand & Prince 1980; Hairer et al. 1987) to calculate the three-body dynamics of the binary–single encounters and binary–MBH encounters and the equations of the motion of KL oscillation. The level of integration accuracy is always below 10⁻¹², which should be small enough for the convergence of the simulation results.

The time step in our simulation should be small enough that the results can converge. We use a step control similar to that of Shapiro & Marchant (1978). The time step δt satisfies

$$\begin{eqnarray}\begin{array}{rcl}\min \left(\delta t\left|{D}_{E}^{\mathrm{NR}}\right|,\sqrt{\delta {{tD}}_{{EE}}^{\mathrm{NR}}}\right) & \leqslant & 0.15| E| ,\\ \max \left(\sqrt{\delta {{tD}}_{{JJ}}^{\mathrm{NR}}},\sqrt{\delta {{tD}}_{{JJ}}^{\mathrm{RR}}}\right) & \leqslant & \min \\ & & \ [0.1{J}_{{\rm{c}}},0.4\left(1.0075{J}_{c}-J\right)],\\ \max \left(\sqrt{\delta {{tD}}_{{JJ}}^{\mathrm{NR}}},\sqrt{\delta {{tD}}_{{JJ}}^{\mathrm{RR}}}\right) & \leqslant & \max \\ & & \ \left(0.25| J-{J}_{\mathrm{lc}}| ,0.1{J}_{\mathrm{lc}}\right).\end{array}\end{eqnarray} \tag{ 20 }$$

Also, if r_p = a₂(1−e₂) < 3r_t, we set δt ≤ P, where P is the orbital period, such that the three-body encounter between the BBH and the MBH is considered for each period. Under these conditions, the step size is small enough that the change of angular momenta is always smaller than the size of the lose cone.

Also, if the binary–single encounters are considered, we require that the time step size is small enough that in each step the number of binary–single collisions is of order unity,

$$\begin{eqnarray}&&\delta t\leqslant {R}_{\mathrm{CO}}^{-1}.\end{eqnarray} \tag{ 21 }$$

Here R_CO is given by Equation (52). If the vector relaxation is considered, to avoid a large change of the orientation of the angular momentum, we set the time step such that the change of angles is less than 30°, i.e.,

$$\begin{eqnarray}&&\delta t\leqslant 0.25{T}_{\mathrm{RR}}^{v},\end{eqnarray} \tag{ 22 }$$

where ${T}_{\mathrm{RR}}^{v}$ is given by Equation (4).

As the number density is dominated by the field stars, the expected density profile of the BBHs should be the same as that of the field stars if m_BBH = m_⋆ and show mass segregation effects if m_BBH > m_⋆ (Bahcall & Wolf 1977). The expected index of the density profile is given by (Alexander & Hopman 2009)

$$\begin{eqnarray}{\alpha }_{\mathrm{BBH}}=\left\{\begin{array}{cc}3/2+{m}_{\mathrm{BBH}}/(4{m}_{\star }) & \mathrm{if}\,{m}_{\mathrm{BBH}}/{m}_{\star }\lesssim 4\\ 9/2-{\alpha }_{\star } & \mathrm{if}\,{m}_{\mathrm{BBH}}/{m}_{\star }\gtrsim 4.\end{array}\right.\end{eqnarray} \tag{ 23 }$$

The yellow solid lines in Figure 4 show the distribution function g_BBH(x) in models M1, M2, M3a, and M4 when only the NR is considered, and without the loss cone. For model M1, in which the mass of the binary m_BBH = 10 ${M}_{\odot }$ is equal to the field stars, the density profile of BBHs is the same as that of the field stars, i.e., α_⋆ = 7/4, consistent with Equation 23. For models M2, M3a, and M4, the density profile of the BBHs also follows Equation (23), consistent with the theoretical expectations.

Zoom In Zoom Out Reset image size

When additionally the loss cone is considered, the results of density profiles in different models are shown with the blue solid lines in Figure 4. We can see that the loss cone can reduce the number of BBHs in the inner regions. The RR process can cause an additional decrease of BBHs in the inner regions (see the red solid lines in Figure 4), as also suggested in other studies (e.g., Hopman & Alexander 2006). This is because the RR can excite the eccentricity of the outer orbit of BBHs and move them more efficiently to the loss cone regions. The decrease of BBHs due to the NR or RR effects is only significant when m_BBH ≃ m_⋆. If ${m}_{\mathrm{BBH}}\gtrsim 4{m}_{\star }$, the BBHs destroyed in the loss cone are quickly replenished by the BBHs moving from other regions under the strong mass segregation effect.

The binary–single encounters can reduce further the number of BBHs in the inner region (see the green lines in Figure 4), mainly because the BBHs in the inner regions become soft binaries in all models. We define a dimensionless parameter ξ (Hopman 2009),

$$\begin{eqnarray}&&\xi =\displaystyle \frac{\epsilon }{{m}_{\star }{\sigma }^{2}}\simeq \displaystyle \frac{{m}_{A}{m}_{B}{a}_{2}}{2{m}_{\star }{M}_{\bullet }{a}_{1}}.\end{eqnarray} \tag{ 24 }$$

Then, the binary is soft if ξ ≪ 1 and hard if ξ ≫ 1. Table 2 shows the value of ξ at a₂ = r_h. We can see that for models M1 and M3a-c, within the cluster the BBHs are mostly soft binaries. For models M2 and M4, the binary is hard near the edge of the cluster. However, they will be soft in the inner regions of the cluster.

The soft binaries are likely to be ionized after encounter. Figure 5 shows the distribution of the SMA of the inner orbit of the BBHs for different models (note that KL and GW orbital decay are still ignored in these figures). We can see that in the hard binary region, i.e., ξ > 1, the binary always becomes harder after encounter. However, in the soft binary region, the binary–single encounters can increase the SMA (a₁) of the inner binary. For model M2, as m_BBH ≫ m_⋆, a significant number of BBHs sink into the inner regions of the cluster under the strong mass segregation effect, where the rates of binary–single encounters are dramatically increased. In all models, in the soft binary region, such encounters can still decrease a₁ down to 10⁻² au, which will help increase the event rates of GW mergers of the BBHs.

Zoom In Zoom Out Reset image size

Similarly, Figure 6 shows the distribution of the eccentricity of the inner orbit of the BBHs. We can see that in all models the eccentricity of the binaries will be dramatically changed owing to the binary–single encounters. The maximum eccentricity can be up to ≳0.99, which will further increase the event rate of GW mergers of the BBHs.

Zoom In Zoom Out Reset image size

The magenta lines in Figure 4 show the density distribution of BBHs when all of the dynamical effects (including the vector RR relaxation, KL oscillations, and the orbital decay due to GW radiations) are considered. Comparing the green lines with the magenta lines in Figure 4, we find that considering the GW orbital decay and KL effects can reduce further the number of BBHs in the inner regions, which means that many BBHs are merged. We find that the vector RR relaxation process has insignificant effects on the number density of BBHs. In all models, the BBHs can survive if they are 10²–10⁴ au away from the MBH. Only in model M2, where the mass segregation effects are the most significant, can the BBHs sink into the more inner region with the distance of ∼100 au from the MBH.

We find that the density profile of BBHs in models with α_⋆ = 1 (e.g., M3b) is similar to those of α_⋆ = 7/4, and only in the case that m_BBH ≫ m_⋆ do we find that the BBHs are more concentrated in the central regions for the case α_⋆ = 1. These simulation results are consistent with what is expected in Equation (23).

There are a number of possible dynamical channels to merge BBHs in the cluster containing an MBH. Here we exemplify three of them, which are presented in Figures 7–9. The details are described as follows.

In the first case, the BBHs experience hundreds to thousands of binary–single encounters and then merge under the combined effect of KL oscillations and the binary–single encounters. This channel is more common for those models with strong mass segregation effects, i.e., m_BBH ≫ m_⋆. Figure 7 shows an example of the evolutions of the BBH's inner and outer orbits in model M2. The outer orbit of the BBH first gradually migrates from the outer parts of the cluster (a₂ ∼ 10⁵ au) into the inner regions (a₂ ∼ 10³ au). Along the way the binary–single encounters become more and more frequent. Each binary–single encounter changes the inner orbit slightly, as the incoming field star is much lighter than the components of the BBHs, i.e., m_A, m_B ≫ m_⋆. From the top right panel of Figure 7, we can see that the encounters gradually increase the e₁ of the BBH up to ∼0.8–0.9, and the SMA of the BBH also increases up to about 2 au. When the pericenter of the BBH's outer orbit approaches to about tens of r_t, the KL effects become important, which can increase e₁ rapidly. Note that simultaneously the binary–single encounters also change the inner orbits. Under these two effects, the inner orbit of the BBH finally reaches to a very high eccentricity, and eventually the BBH merges owing to GW radiation.

Zoom In Zoom Out Reset image size

The bottom panels of Figure 7 show the evolutions of the orbital orientations of both the inner and outer orbits of the BBH. The orientation of the outer orbit evolves under vector RR and the KL effect (bottom left panel of Figure 7). The orientation of the inner orbit is changed cumulatively during multiple binary–single encounters and becomes quite rapid during the KL precesses (bottom right panel of Figure 7). The rapid oscillations of the inner orbital orientations can help to trigger the merging of BBHs by the KL effect. For example, before merging, I₂ ∼ π and I₁ ∼ π/2, the relative inclination between the inner and outer orbits is close to π/2, which can help to enhance the eccentricity oscillations and merging of the BBHs (Wen 2003).

Figure 8 shows another example of the BBH's evolutions in model M2. In Figure 8, the SMA of the BBH always tends to decrease owing to multiple binary–single encounters. This is because the BBH is mostly remaining in the hard binary regions, and the binary–single encounters always decrease the SMA for a hard binary. A hard binary makes the frequency of binary–single encounters drop and prevents another encounter even if it reaches the inner regions.

Zoom In Zoom Out Reset image size

Different from model M2, model M3a assumes more massive background field stars. Such initial conditions have two consequences: (1) the change of the inner orbit of the BBHs can be dramatic, and (2) the number of encounters is much smaller in M3a than in M2, as the number of field stars is much smaller (due to their massiveness) according to Equation (14) or Equation (15). In these cases, the BBHs can experience only a few binary–single encounters, which will dramatically increase the eccentricity of the binary and soon lead to merger owing to the GW radiation.

In some rare cases, the GW merging event can be triggered owing to the multiple encounters between the BBH and the MBH, as shown in Figure 9. In the final revolution, from the bottom panels of Figure 9 we find that the relative inclination between the plane of the inner and outer orbits is ∼π/2. As the tidal force of the MBH is nearly perpendicular to the orbital plane of the BBH, it tends to suppress the BBH (Zhang et al. 2010) and help to merge the BBH.

Zoom In Zoom Out Reset image size

This usually happens when the BBH is at the outer parts of the cluster, such that it can approach the loss cone region without being ionized by the binary–single encounters and that the effect of KL oscillations is not so significant. However, we notice that if the binary–single encounters are switched off, then this case is very frequent. The eccentricities of the BBHs can be excited to very high values owing to the multiple encounters between the BBH and MBH, and the GW is so strong that they can merge before the next encounters with the MBH.

In all models, we find that the SB does not affect much the steady-state flux of binaries into the loss cone, which may be important for those single stars (Bar-Or & Alexander 2016). As we can see from Figures 8 and 9, in most cases the SB is below the tidal radius of BBHs with a₁ ≳ 0.1 au, and the SB is only effective if a₂ ≲ 300 au, where most BBHs cannot penetrate into. As RR is effective below a₂ ∼ 10⁴ au (See Figure 1), the tidal rates of BBHs are not affected much by the SB, and the RR process can still drive most of the BBHs effectively into loss cone regions.

At the end of the simulation, suppose that the number of integrated binaries that survived in the simulation is N and the rate of merging events is ${\dot{N}}_{e}$ = dN_e/dt, where dN_e is the total number of merging events during time span dt. We can define a normalized merging event rate in the nucleus cluster, ${\mathscr{R}}={\dot{N}}_{e}/N$. It means that if we observed only one BBH in the real galactic nucleus clusters at the current moment, the merging event rate in such a cluster is then given by ${\mathscr{R}}$. Note here that N and ${\dot{N}}_{e}$ are calculated after considering the weight of each event in the clone scheme (see Section 2.5).

Based on ${\mathscr{R}}$, the merging rate of the BBHs in a real galactic nucleus can then be estimated after some scaling with the realistic number of BBHs in each galactic nucleus. If assuming a constant number fraction, i.e., f_nb, of BBHs in the cluster with respect to the field stars, the event rate is given by

$$\begin{eqnarray}&&R=\displaystyle \frac{{f}_{{nb}}{M}_{\bullet }}{{m}_{\star }}{\mathscr{R}}.\end{eqnarray} \tag{ 25 }$$

If alternatively assuming a constant mass fraction, i.e., f_mb, of BBHs with respect to the field stars, it becomes

$$\begin{eqnarray}&&R=\displaystyle \frac{{f}_{{mb}}{M}_{\bullet }}{\langle {m}_{\mathrm{BBH}}\rangle }{\mathscr{R}},\end{eqnarray} \tag{ 26 }$$

where $\langle {m}_{\mathrm{BBH}}\rangle$ is the mean mass of the BBHs.

In this work, we assume that the mass and the number fraction of all stellar black holes are much smaller than the field stars. If assuming a constant mass fraction, we set f_mb = 10⁻³. If assuming a constant number fraction of BBHs, we set ${f}_{{nb}}=\min ({10}^{-2}{m}_{\star }/\langle {m}_{\mathrm{BBH}}\rangle ,{10}^{-3})$, such that the total mass of the BBHs does not exceed ∼0.01 ${M}_{\bullet }$, which is likely the mass fraction of stellar-mass black holes in a nucleus cluster under standard initial mass function of stars (Hopman & Alexander 2006).

Table 3 shows the estimated GW event rates according to Equations (25) and (26) in different models. The first row for each model shows the result if assuming a constant number fraction of BBHs and ${f}_{{nb}}=\min ({10}^{-3},{10}^{-2}{m}_{\star }/\langle {m}_{\mathrm{BBH}}\rangle )$, while the second row assumes a constant mass fraction of BBHs and f_mb = 10⁻³. We can see that, if all the dynamical effects are considered, the event rates range from 1 to 10³ Gyr⁻¹ in the six models investigated here. Merging rates assuming a constant number fraction are usually larger than those assuming a constant mass fraction, as usually the field stars are considered lighter than the BBHs. The change of α_⋆ parameter does not affect significantly the merging rates, as it changes only the density profiles in the cases that m_BBH ≫ m_⋆. If alternatively the initial eccentricity is not zero, but follows the thermal distribution, the merging rates will be about twice larger. This can be straightforwardly obtained as the larger the eccentricity, the easier the BBHs can be merged.

Table 3.  Dependence of the Merging Rate of BBHs on Different Dynamical Effects

Model	ALL	Vec RR Off	KL Off	BySn Off
M1	3.3 ± 0.2a	2.5 ± 0.2	2.6 ± 0.2	1.5 ± 0.2
	3.3 ± 0.2b	2.5 ± 0.2	2.6 ± 0.2	1.5 ± 0.2
M2	1046 ± 76	939 ± 72	899 ± 70	362 ± 48
	104.6 ± 7.6	93.9 ± 7.2	89.9 ± 7.0	36.2 ± 4.8
M3a	11.1 ± 0.3	11.0 ± 0.3	10.5 ± 0.3	3.3 ± 0.2
	5.5 ± 0.1	5.5 ± 0.1	5.3 ± 0.1	1.7 ± 0.1
M3b	9.0 ± 0.4	7.8 ± 0.4	8.2 ± 0.4	3.3 ± 0.2
	4.5 ± 0.2	3.9 ± 0.2	4.1 ± 0.2	1.7 ± 0.1
M3c	22.1 ± 0.9	19.6 ± 0.8	20.3 ± 0.8	6.9 ± 0.5
	11.0 ± 0.4	9.8 ± 0.4	10.1 ± 0.4	3.4 ± 0.2
M4	169 ± 10.6	175 ± 11	163 ± 10	48 ± 5.7
	28.1 ± 1.8	29.3 ± 1.8	27.3 ± 1.7	8.0 ± 1.0

Notes. Note that the rates in this table are only for the purpose of demonstrating the effects of different dynamical processes. For each model, we have assumed Milky Way–like galaxies with MBH ${M}_{\bullet }$ = 4 × 10⁶ ${M}_{\odot }$ and have oversimplified initial conditions of the BBHs (see Table 2). For more realistic event rate estimations in the local universe, see Section 4 or Table 5.

^aMerging rates in units of Gyr⁻¹, assuming a constant number fraction of BBHs, i.e., ${f}_{{nb}}=\min ({10}^{-3},{10}^{-2}\,{m}_{\star }/\langle {m}_{\mathrm{BBH}}\rangle )$. ^bMerging rates in units of Gyr⁻¹, assuming a constant mass fraction of BBHs, i.e., f_mb = 10⁻³.

Download table as:  ASCIITypeset image

We find that the exchange of the BBH with the incoming stars is common in some models, e.g., M1 and M3a, as usually the exchange event is frequent if m_⋆ ∼ m_BBH. These exchange events decrease the total GW event rate significantly. By removing the exchanged BBHs from the simulation, we find that the GW event rate is about twice smaller.

As we have shown in Figures 7–9, the mergers of BBHs are results of the combinations of various dynamical effects. To explore their importance in the merger of BBHs, we also run similar simulations for each model, but excluding some dynamical processes. The results are shown in Table 3. We do not see a clear difference when the vector RR effect is being considered or not. If the vector RR effects are switched off, the event rates usually drop only slightly for all models. As we consider only ${M}_{\bullet }$ = 4 × 10⁶ ${M}_{\odot }$ in this section, this is consistent with Hamers et al. (2018), that the vector RR effect is more important for MBHs with ${M}_{\bullet }$ = 10⁴–10⁵ ${M}_{\odot }$ but quickly becomes negligible with increasing M_•.

Many of the BBHs are ionized or tidally disrupted along the way from the outer to the inner regions. Thus, the KL effect does not contribute much to the total merging rate, as few of the BBHs can survive to the vicinity of MBHs where it is most efficient. If the KL effects are switched off, we do not see a significant difference on the GW merging rates (see Table 3). On the other hand, many of the BBHs may coalesce owing to multiple binary–single encounters. Our results suggest that the binary–single encounters can be an important channel that leads to the mergers of the BBHs.

If we switch off the binary–single effects, then the inner orbits of BBHs are affected only by KL oscillations or the strong encounters with the MBH. We find that in these cases most of the BBHs will merge after the excitement of eccentricity owing to both the KL oscillations and the BBH–MBH encounters. Note that the contribution from the latter one to the merging rates is only significant when the binary–single encounters are switched off. However, as these two effects are less efficient in merging the BBHs than those of the binary–single encounters, there is a significant drop in the merging rates for all models (see results in Table 3).

We find that the BBHs in our simulations are all affected by the dynamical effects mentioned above before they merged. We do not find any BBHs ending up with a merger by evolving isolatedly⁶ for all models except M3c. In model M3c, the merging rate of BBHs evolved isolatedly is ∼0.5 Gyr⁻¹. Thus, the merging rates in Table 3 are unique results for BBHs in Galactic nuclei and are not affected by BBHs evolved in isolation.

The BBH maintains some eccentricity when the frequency of its gravitational wave enters into the LIGO band, i.e., f_GW = 10 Hz. We denote the eccentricity at this moment as e_{10 Hz}. For BBHs in the galactic field, their orbital eccentricities are very close to zero when entering the LIGO band (Peters 1964). If a BBH has a significant residual eccentricity, say, e_{10 Hz} ≳ 0.1, the modulation in the waveform will be easily recognized by the matched-filtering techniques (Hinderer & Babak 2017). In contrast, as we will see below, the BBHs from galactic nuclei have a non-negligible possibility to have a significant e_{10 Hz}. This will be a smoking-gun effect to distinguish BBH mergers from galactic fields and nucleus clusters when enough events are detected to perform a statistical study.

As we mentioned in earlier sections, the excitement of the eccentricity of BBHs is mainly through binary–single encounters, BBH–MBH multiple encounters, or KL oscillations. Different channels excite the eccentricity of the BBHs in different ways and to different degrees, and thus the distribution of e_{10 Hz} depends on each or combinations of the above three effects: (1) Binary–single encounters excite e₁ dramatically if m_⋆ ∼ m_BBH, and slowly if m_⋆ ≪ m_BBH. When this effect is switched on, the excitement of e₁ due to BBH–MBH encounters mentioned below will be significantly suppressed. (2) BBH–MBH encounters excite e₁ dramatically if the encounter is close, e.g., r_p ≲ r_t, and smoothly if ${r}_{p}\gtrsim {r}_{t}$ (Zhang et al. 2010). Models with massive field stars have large perturbations in the outer orbits of BBHs, and the BBHs can move fast in the lose cone region. Thus, for these models their excitement of e₁ can be dramatic. (3) The KL effect can excite e₁ to very high values, however, only if some preferred orbital configurations are satisfied (especially for the inclination angle). In most cases, the KL effect changes e₁ very smoothly.

Figure 10 shows their probability distribution function (pdf) of e_{10 Hz} when all dynamical effects are considered or without some of the specific ones. Table 4 shows the fraction of BBHs merged with e_{10 Hz} > 0.01 (or ${e}_{10\mathrm{Hz}}\gt 0.05)$ in different models. When all the dynamical effects are included, we find that ∼2%–15% (0.3%–5%) of the BBHs have e_{10 Hz} > 0.01 (e_{10 Hz} > 0.05). M1 has the highest probability, as in this model the binary–single encounters excite e₁ very dramatically. Note that KL effects can affect the e₁ of BBHs; however, it is not the most possible final merging channel in all models. Thus, switching it off does not affect much of the pdf of e_{10 Hz}. We also do not find significant difference by switching vector RR off in the simulation.

Zoom In Zoom Out Reset image size

If the effect of binary–single encounter is switched off, we can see that for all models (except M2), the e_{10 Hz} is higher than the case when the binary–single encounter is switched on. For these models (except M2), the BBH–MBH multiple encounters can excite the eccentricity of BBHs very dramatically, such that e₁ can be very high. For model M2, as m_⋆ ≪ m_BBH, the outer orbits of BBHs are perturbed very smoothly, such that they move in and out of the KL region very frequently (3r_t < r_p < min[20, 0.58(T_GR/P2)^2/3]r_t). In some rare cases the KL excites the e₁ of BBHs to very high values, and thus they merge quickly. However, in most cases the BBHs move out of the KL region, and their eccentricities are changed only to median values. They then evolve isolatedly and merge before it enters the KL region again; thus, the overall e_{10 Hz} remains low.

The results of normalized merging rates for single galaxies are shown in Table 5. We can see that ${\mathscr{R}}$ is a decreasing function of ${M}_{\bullet }$ and m_⋆ (see also the top panel of Figure 11). This is mainly because the larger the mass of the central MBH and the mass of the field stars, the softer the BBHs becomes, and thus the easier for them to be disrupted owing to binary–single encounters. On the other hand, the tidal radius increases with the MBH mass; thus, the BBHs will be easier to disrupt around high-mass MBHs. Figure 12 shows the distribution of the position where the BBHs merged. We can see that for an MBH, the BBH only merges at the outskirt of the nuclei, but for an MBH with a smaller mass, the BBH can merge in the inner parts of the cluster. Given the same MBH, if the mass of field stars is smaller, the BBHs can reach to deeper regions of the cluster.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To estimate the total event rate in the nearby universe, we need to integrate the merge events over all galaxies. The number density of MBHs in the universe is given by (Aller & Richstone 2002)

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{{dM}}_{\bullet }}={c}_{0}{\left(\displaystyle \frac{{M}_{\bullet }}{{m}_{0}}\right)}^{-1.25}\exp \left(-\displaystyle \frac{{M}_{\bullet }}{{m}_{0}}\right),\end{eqnarray} \tag{ 27 }$$

where c₀ = 3.2 × 10⁻¹¹ ${M}_{\odot }$⁻¹ Mpc⁻³ and m₀ = 1.3 × 10⁸ ${M}_{\odot }$. The cosmological event rate is then

$$\begin{eqnarray}&&{R}_{\mathrm{tot}}={\int }_{{10}^{5}\,{M}_{\odot }}^{{10}^{9}\,{M}_{\odot }}R({M}_{\bullet })\displaystyle \frac{{dn}}{{{dM}}_{\bullet }}{{dM}}_{\bullet }.\end{eqnarray} \tag{ 28 }$$

Here R is given by Equation (25) or Equation (26), depending on the model assumption, and we have assumed negligible merging events from a nucleus cluster with ${M}_{\bullet }$ > 10⁹ ${M}_{\odot }$ (see Figure 11).

To obtain the total event rates by Equation (28), we need to expand the results of ${\mathscr{R}}$ shown in Table 5 to arbitrary MBH masses from 10⁵ to 10⁹ ${M}_{\odot }$. For MPW10 and MUN10 model series, we have

$$\begin{eqnarray}&&\mathrm{log}\left({\mathscr{R}}/{\mathrm{Gyr}}^{-1}\right)\simeq -1.14-0.51\mathrm{log}{M}_{7}+0.04{\left(\mathrm{log}{M}_{7}\right)}^{2}\end{eqnarray} \tag{ 29 }$$

and

$$\begin{eqnarray}&&\mathrm{log}\left({\mathscr{R}}/{\mathrm{Gyr}}^{-1}\right)\simeq -0.83-0.61\mathrm{log}{M}_{7}+0.038{\left(\mathrm{log}{M}_{7}\right)}^{2},\end{eqnarray} \tag{ 30 }$$

respectively, where M₇ = ${M}_{\bullet }$/10⁷ ${M}_{\odot }$. For MUR10 and MUN1 model series, we have

$$\begin{eqnarray}&&\mathrm{log}\left({\mathscr{R}}/{\mathrm{Gyr}}^{-1}\right)\simeq -0.045-0.66\mathrm{log}{M}_{7}+0.030{\left(\mathrm{log}{M}_{7}\right)}^{2}\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&\mathrm{log}\left({\mathscr{R}}/{\mathrm{Gyr}}^{-1}\right)\simeq -0.34-0.85\mathrm{log}{M}_{7}+0.0024{\left(\mathrm{log}{M}_{7}\right)}^{2},\end{eqnarray} \tag{ 32 }$$

respectively.

The estimated total rates of BBH merging events are shown in the last two columns of Table 5. If assuming a constant number fraction of BBHs, given r_i = r_h and m_⋆ = 10 ${M}_{\odot }$, we find that the total merging rate is 1.6 and 3.2 Gyr⁻¹ for models assuming PW and UN mass functions, respectively. However, the rate can be up to ∼31 Gyr⁻¹ for a UN model and given m_⋆ = 1 ${M}_{\odot }$. Assuming a small m_⋆ can increase the merging rates as the BBHs become harder in the cluster and can survive in inner regions of the cluster after multiple binary–single encounters. The merging event rates can also be increased to ∼20 Gyr⁻¹ for a UN model but given r_i = 0.1r_h and m_⋆ = 10 ${M}_{\odot }$. As the BBHs are initially located at a distance much closer to the MBH, the KL effects become much more effective in merging these BBHs, resulting in a significant increase of the merging rates. In these models, we find that the contribution of merging rates from KL oscillations is much larger than those from binary–single encounters.

If assuming a constant mass fraction of BBHs, the rates will be dramatically reduced. This is because $\langle {m}_{\mathrm{BBH}}\rangle$ is usually about 20–40 ${M}_{\odot }$, depending on the assumed mass function (see in Table 5). Such a mass is usually larger (or much larger) than the mass of the field stars and reduces significantly the number of BBHs in each galaxy. Nevertheless, the event rate is on order of 1–10 Gpc⁻³ yr⁻¹, which is still not negligible for LIGO detections. These results suggest that the BBH mergers in the center of galaxies can contribute partially to the LIGO observations.

We notice that for MBHs with 10⁸–10⁹ ${M}_{\odot }$, the two-body relaxation timescale in the cluster is much longer than the Hubble timescale. Thus, in reality these clusters may never reach the equilibrium state. Our estimated merging rates for these clusters could be problematic. Nevertheless, the merging rates contributed from these galaxies are small (See Figure 11), and our estimations of the total merging rates should not be significantly affected.

Figure 13 shows the mass distributions of the merged BBHs in different models. We can see that, compared with the initial mass distribution (the dashed magenta lines in each panel), the merged BBHs are likely more massive. For galaxies containing smaller MBHs with ${M}_{\bullet }$ ≲ 10⁷ ${M}_{\odot }$, the merged BBHs are also likely more massive than those galaxies containing larger MBHs with ${M}_{\bullet }$ ≳ 10⁸ ${M}_{\odot }$. More massive BBHs can survive longer in the galaxies with smaller MBHs, and thus they have higher event rates for LIGO. Such preference for massive BBHs is less significant if the BBHs can penetrate into the inner regions of the cluster, as both the massive and less massive BBHs can be quickly merged in the inner regions, as shown in the bottom panels of Figure 13 (for models MUN1 and MUR10).

Zoom In Zoom Out Reset image size

Figure 14 shows the distribution of the eccentricity of the BBHs when the GW frequency approaches f_GW = 10 Hz. We can see that the eccentricities for low-mass MBHs are relatively higher than those around high-mass MBHs. For MBH with mass 10⁵ ${M}_{\odot }$, we find that ∼10%–20% of the merged BBHs have e_{10 Hz} > 0.01. Many BBHs are hard in these cluster, and they are likely more compact after each binary–single encounter. Thus, the eccentricity of them in the merging phase can be high. On the contrary, around higher-mass MBHs the merged BBHs are commonly less eccentric, as many of them are soft binaries, and a significant number of them are merging by evolving in isolation. If the MBH is more massive, i.e., 10⁹ ${M}_{\odot }$, the fraction will be about ≲1%. We also find that if the BBHs are initially located at inner regions, e.g., r_i = 0.1r_h in model MUR10 series, the merging eccentricities can be significantly higher, as the KL oscillations can be much more effective.

Zoom In Zoom Out Reset image size

The percentage of BBH mergers with e_{10 Hz} > 0.01 will be around 3%–10% for all galaxies in the local universe (see the last column of Table 4). As a comparison, the expected percentage of e_{10 Hz} > 10⁻³ in globular clusters is ∼1% (Rodriguez et al. 2016a). For BBH mergers from galactic fields, the final merging eccentricity is practically zero (Peters 1964; Belczynski et al. 2016). Thus, in principle they can be distinguished from BBH mergers from the galactic fields and globular clusters using eccentric waveforms when enough events are accumulated (Hinderer & Babak 2017).

The diffusion coefficients D_EE and D_JJ describe the orbit-averaged scatterings of the energy and angular momentum of the orbit for point-like particles around the MBH. Similarly, D_E and D_J describe the orbit-averaged drifts of the energy and angular momentum, respectively, and D_EJ describes the correlations between the two. The calculations of the diffusion coefficients have been done and discussed in many studies (e.g., Cohn & Kulsrud 1978; Shapiro & Marchant 1978; Bar-Or & Alexander 2016). Here we use the formalism in Bar-Or & Alexander (2016), where the mass of the binary can be different from the field stars. Denote the SMA and the eccentricity of the particle cycling around the MBH as a₂ and e₂, and suppose that the dimensionless distribution function of the field stars is given by f(E). Denote m_BBH and m_⋆ as the masses of the binary and the field star, respectively. Then, first the function Γ_ijk and Γ₀ are calculated,

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{{ijk}} & = & {2}^{1+k-i}\displaystyle \frac{\kappa }{\pi }{\displaystyle \int }_{-1}^{1}{dy}{\displaystyle \int }_{1}^{2/(1+y)}{{dsf}}_{a}({sE})\\ & & \times \ \displaystyle \frac{1}{\sqrt{1-{y}^{2}}}\displaystyle \frac{{\left(1+{{ye}}_{2}\right)}^{i}}{{\left({v}^{2}/E\right)}^{k}}{\left({v}_{a}/v\right)}^{j}\\ {{\rm{\Gamma }}}_{0} & = & \kappa {\displaystyle \int }_{-\infty }^{1}{{dsf}}_{a}({sE}).\end{array}\end{eqnarray} \tag{ 33 }$$

Here κ = (4πGm_⋆)² ln Λ, y = (r/a₂ − 1)/e, ${v}_{a}^{2}$ = 2E(2a₂/r−s), v² = 2E(2a₂/r − 1), and Λ = ${M}_{\bullet }$/m_⋆. Then, the diffusion coefficients are given by

$$\begin{eqnarray}\begin{array}{rcl}{D}_{E}/E & = & \displaystyle \frac{{m}_{\mathrm{BBH}}}{{m}_{\star }}{{\rm{\Gamma }}}_{110}-{{\rm{\Gamma }}}_{0},\\ {D}_{{EE}}/{E}^{2} & = & \displaystyle \frac{4}{3}{{\rm{\Gamma }}}_{13-1}+{{\rm{\Gamma }}}_{0},\\ {D}_{J}/({{jJ}}_{c}) & = & \displaystyle \frac{5-3{j}^{2}}{12}{{\rm{\Gamma }}}_{0}-{j}^{2}\displaystyle \frac{{m}_{\mathrm{BBH}}+{m}_{\star }}{2{m}_{\star }}{{\rm{\Gamma }}}_{111}\\ & & +\ {{\rm{\Gamma }}}_{310}-\displaystyle \frac{1}{3}{{\rm{\Gamma }}}_{330},\\ {D}_{{JJ}}/{J}_{c}^{2} & = & \displaystyle \frac{5-3{j}^{2}}{6}{{\rm{\Gamma }}}_{0}+\displaystyle \frac{{j}^{2}}{2}{{\rm{\Gamma }}}_{131}-\displaystyle \frac{{j}^{2}}{2}{{\rm{\Gamma }}}_{111}\\ & & +\ 2{{\rm{\Gamma }}}_{310}-\displaystyle \frac{2}{3}{{\rm{\Gamma }}}_{330},\\ {D}_{{EJ}}/({{EJ}}_{c}) & = & -\displaystyle \frac{2}{3}j({{\rm{\Gamma }}}_{0}+{{\rm{\Gamma }}}_{130}).\end{array}\end{eqnarray} \tag{ 34 }$$

Here j = J/J_c is the dimensionless angular momentum, and ${J}_{c}=\sqrt{{{GM}}_{\bullet }{a}_{2}}$ is the maximum angular momentum. Note that j could be neither zero nor one; otherwise, the integration will be divergent.

The scalar and vector RRs change only the angular momentum of the outer orbit of the binary. For the scalar RR, we take formalisms similar to Bar-Or & Alexander (2016). Defining Q = ${M}_{\bullet }$/m_⋆ as the number of field stars and ν_r = 2π/P(a₂) = ${a}_{2}^{-3/2}{M}_{\bullet }^{1/2}{G}^{1/2}$ as the Keplerian orbital frequency, the RR torque τ_N is given by

$$\begin{eqnarray}&&{\tau }_{N}=0.28\sqrt{1-j}\sqrt{N(2{a}_{2})}{{Gm}}_{\star }/{a}_{2}.\end{eqnarray} \tag{ 35 }$$

Here N(a₂) = ${M}_{\bullet }$/m_⋆(a₂/r_h)^3−α. Denoting ν_J = τ_N/J_c, the diffusion coefficients are given by

$$\begin{eqnarray}&&{D}_{{JJ}}^{\mathrm{RR}}=\displaystyle \frac{2{J}_{c}^{2}{\nu }_{J}^{2}(j){T}_{c}({a}_{2})}{1+{\left[{T}_{c}({a}_{2}){\nu }_{p}({a}_{2},j)\right]}^{2}}.\end{eqnarray} \tag{ 36 }$$

Here

$$\begin{eqnarray}&&{\nu }_{p}=| {\nu }_{M}+{\nu }_{\mathrm{GR}}| \end{eqnarray} \tag{ 37 }$$

is the precession frequency, where

$$\begin{eqnarray}&&{\nu }_{\mathrm{GR}}=\displaystyle \frac{3{r}_{g}}{{a}_{2}{j}^{2}}{\nu }_{r}({a}_{2})\end{eqnarray} \tag{ 38 }$$

is the frequency of the GR-induced precession. ν_M has an exact form when α = 2,

$$\begin{eqnarray}&&{\nu }_{{\rm{M}}}(\alpha =2,j)=-\displaystyle \frac{N({a}_{2})}{Q}\displaystyle \frac{j}{j\,+\,1}{\nu }_{r}.\end{eqnarray} \tag{ 39 }$$

For other values, it can be obtained numerically. We find that

$$\begin{eqnarray}&&{\nu }_{{\rm{M}}}(\alpha ,j)=-\displaystyle \frac{N({a}_{2})}{Q}\displaystyle \frac{{j}^{2-\gamma }(1-{j}^{\gamma })}{1-{j}^{2}}{\nu }_{r}.\end{eqnarray} \tag{ 40 }$$

When α = 7/4, γ = 1.48; when α = 1/2, we find that γ = 4.0; when α = 2, γ = 1, which reduces to Equation (39); when α = 1, we find that γ = 3.37.

T_c(a₂) is given by (note that the coherence time T_c here has not included the GR precession; Bar-Or & Alexander 2016)

$$\begin{eqnarray}&&{T}_{c}({a}_{2})=\sqrt{\pi /2}{\nu }_{M}^{-1}(2{a}_{2},\sqrt{1/2}).\end{eqnarray} \tag{ 41 }$$

The drift term is given by

$$\begin{eqnarray}&&{D}_{J}^{\mathrm{RR}}=\displaystyle \frac{1}{2J}\displaystyle \frac{\partial ({{JD}}_{{JJ}}^{\mathrm{RR}})}{\partial J}.\end{eqnarray} \tag{ 42 }$$

For the vector RR relaxation, we take a very simple model, where the direction of the angular momentum is randomly walking on the surface of the sphere, with the unity change happening on a timescale given by ${T}_{\mathrm{RR}}^{v}$ in Equation (4). We generate a random vector $\hat{{dj}}$ with magnitude $l={({dt}/{T}_{\mathrm{RR}}^{v})}^{1/2}$, where l is the standard deviation of a lognormal distribution. The unit angular momentum of the CM is changed to $\hat{j}^{\prime} =\hat{j}+\hat{{dj}};$ the direction of $\hat{{dj}}$ is set such that $\hat{j}^{\prime}$ is always a unit vector. Here $\hat{j}$ is the unit angular momentum. First, we solve a reference vector $\hat{{{dj}}_{0}}=(0,l\cos s,l\sin s)$, where $| \hat{j}+{\hat{{dj}}}_{0}| =1$ (if ${j}_{y}^{2}\,+{j}_{z}^{2}\lt {l}^{2}/4$, we can set ${\hat{{dj}}}_{0}=(l\cos s,0,l\sin s)$). Then, the random vector $\hat{j}^{\prime}$ is generated after rotating $\hat{j}+\hat{{{dj}}_{0}}$ around the vector $\hat{j}$ by a random angle ϕ ∈ (0, 2π). If l > 1, we simply set the angular momentum randomly distributed.