Black Hole Mergers from Hierarchical Triples in Dense Star Clusters

In this study, we consider triple systems that satisfy the stability condition of Mardling & Aarseth (2001):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{a}_{\mathrm{out}}}{{a}_{\mathrm{in}}}\gt \displaystyle \frac{3.3}{1-{e}_{\mathrm{out}}}\left[\displaystyle \frac{2}{3}\left(1+\displaystyle \frac{{m}_{2}}{{m}_{\mathrm{bin}}}\right)\right.\\ \quad \times {\left.\displaystyle \frac{1+{e}_{\mathrm{out}}}{{\left(1-{e}_{\mathrm{out}}\right)}^{1/2}}\right]}^{2/5}(1-0.3I/\pi ).\end{array}\end{eqnarray} \tag{ 1 }$$

If this stability criterion is satisfied, a system is long lived, and thus the Hamiltonian of the system can be described as two separate orbits with a perturbative interaction term (Kozai 1962; Lidov 1962). Assuming a circular outer orbit, one may find the following conserved quantity involving eccentricity, mutual inclination, and argument of pericenter (Antognini 2015):

$$\begin{eqnarray}&&{C}_{\mathrm{KL}}={e}_{\mathrm{in}}^{2}\left(1-\displaystyle \frac{5}{2}{\sin }^{2}I{\sin }^{2}{\omega }_{\mathrm{in}}\right).\end{eqnarray} \tag{ 2 }$$

This conserved quantity has an explicit dependency on the inner eccentricity, inner argument of pericenter, and mutual inclination. As a result, the outer binary induces inclination and eccentricity oscillations as well as pericenter precession. This occurs on a timescale of

$$\begin{eqnarray}&&{T}_{\mathrm{KL}}\approx \displaystyle \frac{8}{15\pi }\displaystyle \frac{{m}_{\mathrm{bin}}}{{m}_{2}}\displaystyle \frac{{P}_{\mathrm{out}}^{2}}{{P}_{\mathrm{in}}}{\left(1-{e}_{\mathrm{out}}^{2}\right)}^{3/2},\end{eqnarray} \tag{ 3 }$$

where ${P}_{\mathrm{in}}$ and ${P}_{\mathrm{out}}$ are the orbital periods of the inner and outer orbits, respectively. This process is known as the KL mechanism. ¹⁰ In this level of approximation (the quadrupole approximation) the KL oscillations can only occur within a mutual inclination window of roughly 40°–140°. Allowing for an eccentric outer orbit necessitates an octupole-level approximation. The strength of the octupole-level interaction with respect to the quadrupole interaction can be quantified by

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{m}_{0}-{m}_{1}}{{m}_{0}+{m}_{1}}\displaystyle \frac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}}\displaystyle \frac{{e}_{\mathrm{out}}}{1-{e}_{\mathrm{out}}^{2}}\end{eqnarray} \tag{ 4 }$$

(e.g., Naoz et al. 2013a). The inclusion of the octupole terms in the equations of motion allows for much more complex dynamical evolution, such as orbit flips, oscillations from nearly coplanar orbits, and chaotic behavior (Katz et al. 2011; Lithwick & Naoz 2011; Naoz et al. 2011, 2013a; Li et al. 2014). As this behavior is qualitatively different from the effects presented by Kozai and Lidov, we refer to this as the eKL mechanism. The full octupole equations of motion can be found in Naoz et al. (2013a).

At 1PN order, relativistic precession will be induced on the inner binary, which happens on a timescale (Eggleton & Kiseleva-Eggleton 2001; Naoz et al. 2013b; Liu et al. 2015; Rodriguez & Antonini 2018)

$$\begin{eqnarray}&&{T}_{\mathrm{GR}}={\left(d{\omega }_{\mathrm{in}}/{dt}\right)}^{-1}={P}_{\mathrm{in}}\displaystyle \frac{{c}^{2}{a}_{\mathrm{in}}(1-{e}_{\mathrm{in}}^{2})}{{{Gm}}_{\mathrm{bin}}}.\end{eqnarray} \tag{ 5 }$$

This precession can quench the eKL oscillations. If ${T}_{\mathrm{GR}}\sim {T}_{\mathrm{KL}}$, then the maximum eccentricity attainable by the inner orbit is limited. If ${T}_{\mathrm{GR}}\ll {T}_{\mathrm{KL}}$, then the KL mechanism is completely suppressed and inner binary eccentricity is unaffected by the tertiary.

At 1.5 PN order, each of the spin vectors ${\boldsymbol{S}}$ ₀ and ${\boldsymbol{S}}$ ₁ of the inner binary precess due to torques from the inner binary. Thus, it is necessary to include the spin–orbit interaction terms (Apostolatos et al. 1994)

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{S}}}_{0}}{{dt}}=\displaystyle \frac{2G\mu \nu }{{c}^{2}{a}_{\mathrm{in}}{j}^{3}}\,\left(1+\displaystyle \frac{3{m}_{1}}{4{m}_{0}}{\boldsymbol{j}}\times {{\boldsymbol{S}}}_{0}\right),\end{eqnarray} \tag{ 6 }$$

where μ is the reduced mass, ν is the Keplerian orbital frequency, and $j=| {\boldsymbol{j}}|$, replacing indices $0\mapsto 1$ for the other BH. For a fixed binary, this describes uniform precession around the binary angular momentum vector. In the presence of the KL oscillation, the binary angular momentum will itself precess, thus allowing for a much more interesting behavior. This can be quantified by the evolution of the binary effective spin parameter

$$\begin{eqnarray}&&{\chi }_{\mathrm{eff}}=\displaystyle \frac{{m}_{0}{\chi }_{0}\cos {\theta }_{0}+{m}_{1}{\chi }_{1}\cos {\theta }_{1}}{{m}_{\mathrm{bin}}},\end{eqnarray} \tag{ 7 }$$

where $\cos {\theta }_{i}=\hat{{{\boldsymbol{S}}}_{i}}\cdot \hat{{\boldsymbol{j}}}$. Though this quantity is conserved for an isolated binary, it evolves over a KL cycle due to the change in the direction of angular momentum. Note that the spin–orbit coupling introduces additional precession on the inner binary, but this depends on the in-plane spin components. Since ${\chi }_{\mathrm{eff}}$ only contains information about the spin components perpendicular to the orbital plane, another quantity must be defined with information about the in-plane spin components, known as the effective precession parameter (Schmidt et al. 2015)

$$\begin{eqnarray}&&{\chi }_{{\rm{p}}}=\max \left\{| {\chi }_{0}\sin {\theta }_{0}| ,\kappa | {\chi }_{1}\sin {\theta }_{1}| \right\},\end{eqnarray} \tag{ 8 }$$

where $\kappa =q(4q+3)/(4+3q)$ and $0\lt q\leqslant 1$ is the mass ratio of the binary. We note that we ignore the backreaction terms, since the binary angular momentum $L\gg S$ during a KL oscillation (e.g., Rodriguez & Antonini 2018).

In this study we neglect the 2 PN order effects. Including or discarding these terms in the secular equations of motion does not greatly change the dynamics of the system (Naoz et al. 2013b).

Finally, at the 2.5 PN order, the inner binary experiences GW radiation reaction, which is described by the following equations (Peters 1964):

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d{\boldsymbol{e}}}{{dt}}\right|}_{\mathrm{GW}}=-\displaystyle \frac{304}{15}\displaystyle \frac{{G}^{3}{m}_{0}{m}_{1}{m}_{\mathrm{bin}}}{{c}^{5}{a}_{\mathrm{in}}^{4}{j}^{5}}\,\left(1+\displaystyle \frac{121}{304}{e}_{\mathrm{in}}^{2}\right){\boldsymbol{e}}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{{da}}_{\mathrm{in}}}{{dt}}\right|}_{\mathrm{GW}}=-\displaystyle \frac{64}{5}\displaystyle \frac{{G}^{3}{m}_{0}{m}_{1}{m}_{\mathrm{bin}}}{{c}^{5}{a}_{\mathrm{in}}^{3}{j}^{7}}\,\left(1+\displaystyle \frac{73}{24}{e}_{\mathrm{in}}^{2}+\displaystyle \frac{37}{96}{e}_{\mathrm{in}}^{4}\right).\end{eqnarray} \tag{ 10 }$$

For an isolated binary, these describe a gradual inspiral through the emission of gravitational radiation. The lifetime of such systems until coalescence can be converted to the following integral:

$$\begin{eqnarray}&&{T}_{\mathrm{GW}}=\displaystyle \frac{12}{19}\displaystyle \frac{{c}_{0}^{4}}{\beta }{\int }_{0}^{{e}_{0}}\displaystyle \frac{{e}^{29/19}{[1+(121/304){e}^{2}]}^{1181/2299}}{{\left(1-{e}^{2}\right)}^{3/2}}{de},\end{eqnarray} \tag{ 11 }$$

where e₀ is the initial eccentricity of the binary and the constants in the prefactor are defined as

$$\begin{eqnarray*}&&{c}_{0}\equiv {a}_{0}{e}_{0}^{-12/19}(1-{e}_{0}^{2}){\left[1+\displaystyle \frac{121}{304}{e}_{0}^{2}\right]}^{-870/2299}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&\beta \equiv \displaystyle \frac{64}{5}\displaystyle \frac{{G}^{3}{m}_{0}{m}_{1}{m}_{\mathrm{bin}}}{{c}^{5}}.\end{eqnarray*}$$

For a typical quasi-circular binary composed of two $30\,\,{M}_{\odot }$ BHs to merge within ${t}_{\mathrm{Hubble}}$, the separation must be within $\sim 0.22\,\mathrm{au}$. However, at extremely high eccentricities, the efficiency of GW radiation is greatly increased, and the merger time greatly reduced. In the case of the example system, for ${a}_{0}=0.22\,\mathrm{au}$ and e₀ = 0.99, the time to merger decreases from ${t}_{\mathrm{Hubble}}$ to $\sim 2.3\times {10}^{4}\,$ yr. For an even more extreme eccentricity of ${e}_{0}=0.999$, the merger time decreases to just 8 yr. While eccentricities this high are usually inaccessible for an isolated compact binary due to circularization during common-envelope evolution, the KL mechanism may be able to drive binaries to such eccentricities, depending on the initial conditions.

We make use of the GC models produced using the Cluster Monte Carlo code. The CMC code is a Hénon-type Monte Carlo code (Hénon 1971, 1975) to treat the long-term evolution of GCs (Joshi et al. 2000, 2001; Fregeau et al. 2003; Fregeau & Rasio 2007; Chatterjee et al. 2010; Umbreit et al. 2012; Pattabiraman et al. 2013; Morscher et al. 2015; Rodriguez et al. 2016, 2018; Kremer et al. 2020b). This includes detailed treatments of stellar evolution via SSE and BSE (Hurley et al. 2000, 2002; Chatterjee et al. 2010) with updated prescriptions for compact object formation, two-body relaxation (Joshi et al. 2000), single–single capture (Samsing et al. 2020), three-body binary formation (Morscher et al. 2013), direct stellar collisions (Fregeau & Rasio 2007), galactic tides (Chatterjee et al. 2010; Pattabiraman et al. 2013), and the direct integration of strong three- and four-body encounters (Fregeau & Rasio 2007). Note that in these cluster models, all BHs are assumed to form with no natal spin (Fuller & Ma 2019). Full details for all the prescriptions used for these CMC models can be found in Kremer et al. (2020b).

In Kremer et al. (2020b), 148 cluster simulations ¹¹ were produced, varying the total number of particles (single stars plus binaries; $N=2\times {10}^{5}$, 4 × 10⁵, 8 × 10⁵, $1.6\times {10}^{6}$, and $3.2\times {10}^{6}$), initial cluster virial radius (${r}_{v}/\mathrm{pc}=0.5,1,2,4$), metallicity ($Z/{Z}_{\odot }=0.01,0.1,1$), and galactocentric distance (${R}_{\mathrm{gc}}/\mathrm{kpc}=2,8,20$). All cluster models initially assume a King potential with King parameter ${W}_{0}=5$ (King 1962) and a primordial binary fraction of 5%. Three of the clusters in the catalog do not form any triples composed entirely of BHs due to collisional runaway before the end of the main-sequence lifetime of the most massive stars.

Primordial triples are not included in our models. However, stable hierarchical triples can be formed as the result of strong binary–binary encounters (Rasio et al. 1995). Current limitations of CMC require that triples are broken up into a binary and a single at the end of each integration time step. Nevertheless, CMC outputs information on the triple, including masses, stellar types, radii, inner, and outer semimajor axes and eccentricity, as well as the formation time and properties of the cluster core at formation time. Since we lack information about the mutual orientation of the two orbits, we compensate by creating 10 different realizations with different mutual orientations (Antonini et al. 2016). In principle, these triple components can be a part of many different triples due to this limitation over the lifetime of the cluster during the simulation.

In our models, we only consider triples produced from binary–binary encounters in the cluster core. We find that the distribution of inner and outer eccentricity for the triples produced by CMC are approximately thermal. Binary–binary encounters can be resonant or nonresonant. In nonresonant encounters, the tighter binary will exchange into the wider binary, with the replaced object kicked out. We find that nonresonant exchange is typical for the BH triples we form in our models. When this happens, according to energy conservation, the newly formed system will have an outer semimajor axis

$$\begin{eqnarray}&&{a}_{\mathrm{out}}\simeq {a}_{\mathrm{wide}}\displaystyle \frac{{m}_{\mathrm{bin}}}{{m}_{\mathrm{esc}}},\end{eqnarray} \tag{ 12 }$$

where ${a}_{\mathrm{wide}}$ is the semimajor axis of the wider binary and ${m}_{\mathrm{esc}}$ is the mass of the ejected single (Sigurdsson & Phinney 1993). For more details, and to see the accuracy of this relation for a wider variety of triple archetypes, see the bottom panel of Figure 1 and the associated text of Paper I. Since the new semimajor axis depends on this mass ratio, the new orbit can become much wider if the single is, for example, a low-mass star. As a result, even if binaries in the cluster core are very compact, the triple orbits can be very wide. In some cases, the binary–binary interaction can be much more complicated, involving resonant interactions (Zevin et al. 2019). In those cases, it is not straightforward to make such predictions for the end state of individual systems. For details, see Paper I.

Zoom In Zoom Out Reset image size

When the single is ejected, the triple will also be imparted a recoil velocity ${v}_{\mathrm{rec}}$. If ${v}_{\mathrm{rec}}\gt {v}_{\mathrm{esc}}$, the escape speed from the cluster, the triple will be completely ejected from the cluster. If on the other hand ${v}_{\mathrm{rec}}\lt {v}_{\mathrm{esc}}$, the triple's cluster-centric radial orbit will have a new apocenter

$$\begin{eqnarray}&&{r}_{\mathrm{apo}}\approx {r}_{{\rm{c}}}\sqrt{\displaystyle \frac{{v}_{\mathrm{esc}}^{4}}{{\left({v}_{\mathrm{esc}}^{2}-{v}_{\mathrm{rec}}^{2}\right)}^{2}}-1},\end{eqnarray} \tag{ 13 }$$

where ${r}_{{\rm{c}}}$ is the core radius of the cluster. Here, we assume that the encounter occurs at the center of the cluster. This is a reasonable assumption because the potential is approximately constant in the cluster core where these BH triples are formed (Sigurdsson & Hernquist 1993). While this formula assumes a Plummer potential for the cluster, it is still valid across a wide range of values of W₀ for the King models that our clusters assume (Antonini et al. 2019).

In the dense environment of a GC, triples may be perturbed through stellar encounters. To account for this, we calculate for each triple the typical encounter time within the cluster core (Ivanova et al. 2008)

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{enc}} & \approx & 8.5\times {10}^{12}\,\mathrm{yr}\,{P}_{\mathrm{out}}^{-4/3}{M}_{\mathrm{tri}}^{-2/3}{\sigma }_{10}^{-1}{n}_{5}^{-1}\\ & & \times \left[1+913\displaystyle \frac{{M}_{\mathrm{tri}}+\langle M\rangle }{2{P}_{\mathrm{out}}^{2/3}{M}_{\mathrm{tri}}^{1/3}{\sigma }_{10}^{2}}\right],\end{array}\end{eqnarray} \tag{ 14 }$$

where ${P}_{\mathrm{out}}$ is the outer binary period in days, ${m}_{\mathrm{trip}}$ is the total mass of the triple in ${M}_{\odot }$, ${\sigma }_{10}$ is the central velocity dispersion of the cluster in units of 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$, $\langle M\rangle$ is the average mass of an object in the cluster in ${M}_{\odot }$, and n₅ is the number density of objects in units of ${10}^{5}\,{\mathrm{pc}}^{-3}$. In order to account for the recoil of triples into orbits that may extend outside the core, we can replace the quantities in Equation (14) describing the cluster with their local versions, which are monotonically declining functions of r (King 1962). Then we can average as follows in order to account for the triple's total orbit:

$$\begin{eqnarray}&&\langle {T}_{\mathrm{enc}}\rangle ={\int }_{0}^{{r}_{\mathrm{apo}}}{r}^{2}{T}_{\mathrm{enc}}(r){dr}.\end{eqnarray} \tag{ 15 }$$

In reality, this approach does not take into account further evolution of the orbit due to dynamical friction, though it is sufficient since mergers due to eKL will typically happen before significant evolution can take place. For a triple that is completely ejected from the cluster, we treat the encounter time as infinite.

We create 10 different realizations of the orbital orientations, sampling the mutual inclination from an isotropic distribution, i.e., uniform in $\cos I$ between 0 and 1 and the other two Euler angles uniformly between 0 and 2π. Finally, we sample the recoil velocities of each triple 10 times. The procedure for this is described in detail in Paper I. Thus, we create for each triple a total of 100 different realizations.

Of the triples in the cluster catalog, 63,508 are composed of three BHs. A small number of the triples are only stable at certain eccentricities, so after resampling each triple 100 times, we once again evaluate their stability using Equation (1), finding a total sample size of 6,090,030. Those that remain stable are integrated forward in time numerically using the publicly available secular code Kozai ¹² (Antonini et al. 2018; Rodriguez & Antonini 2018) until either the integration time reaches $\min ({t}_{\mathrm{Hubble}}-{t}_{\mathrm{formation}},1000{T}_{\mathrm{KL}},\langle {T}_{\mathrm{enc}}\rangle )$ or the triple reached the LVK frequency band ${f}_{\mathrm{GW}}=10\,\mathrm{Hz}$, producing a merger at time ${T}_{\mathrm{merger}}$. We also keep track of the eccentricity of the inner binary when the peak frequency passes through the values 0.01 Hz, 1 Hz, and 10 Hz, characteristic frequencies for the LISA, DECIGO, and LVK detectors, respectively. We calculate the frequency as the highest harmonic produced by the inner binary as defined by Wen (2003):

$$\begin{eqnarray}&&{f}_{\mathrm{GW}}=\displaystyle \frac{\sqrt{{{Gm}}_{\mathrm{bin}}}}{\pi }\displaystyle \frac{{\left(1+{e}_{\mathrm{in}}\right)}^{1.1954}}{{[{a}_{\mathrm{in}}(1-{e}_{\mathrm{in}}^{2})]}^{1.5}}.\end{eqnarray} \tag{ 16 }$$

In Figure 4, we show the eccentricity spectrum produced by different populations of merging binaries at a range of frequencies representing the peak frequencies of LVK (10 Hz), DECIGO (1 Hz), and LISA (0.01 Hz). In addition to the KL mergers and non-KL mergers from this study, we also show the binary mergers from the same cluster models, subdivided into four different categories (e.g., Samsing & D'Orazio 2018; Zevin et al. 2019; Kremer et al. 2020b). The ejected mergers are mergers of binaries ejected from the cluster (Portegies Zwart & McMillan 2000). In-cluster mergers are defined as mergers of binaries formed by dynamical encounters that occur within the cluster (Rodriguez et al. 2018). Few-body captures are the mergers that occur during a resonant three- or four-body encounter (Samsing et al. 2014; Zevin et al. 2019). Finally, single–single captures occur when two BHs on initially hyperbolic orbits come sufficiently close for GW radiation to create a binary (Samsing et al. 2020). We also include a population of field binaries from Kremer et al. (2019), which in turn were computed using COSMIC (Breivik et al. 2020).

Zoom In Zoom Out Reset image size

Note that the capture merger channels form binaries with initially high frequency. In order to compute the formation frequency for each binary, we follow the procedure outlined in detail in Zevin et al. (2019). In short, we use Equation (16) and

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{{{da}}_{\mathrm{in}}}{{{de}}_{\mathrm{in}}}\right\rangle =\displaystyle \frac{12}{19}\displaystyle \frac{{a}_{\mathrm{in}}}{{e}_{\mathrm{in}}}\displaystyle \frac{[1+(73/24){e}_{\mathrm{in}}^{2}+(37/96){e}_{\mathrm{in}}^{4}]}{(1-{e}_{\mathrm{in}}^{2})[1+(121/304){e}_{\mathrm{in}}^{2}]},\end{eqnarray} \tag{ 17 }$$

which is found from the coupled Equations (9)–(10) in Peters (1964). These equations are not differentiable at a given frequency, ${f}_{\mathrm{GW}}$, when they form above that frequency, as $e({f}_{\mathrm{GW}})\gt 1$. Therefore, we compute the pericenter distance R_p at a reference formation eccentricity $1-e={10}^{-3}$ and compare this to the semimajor axis ${a}_{\mathrm{circ}}$ of the binary if it were on a circular orbit with frequency ${f}_{\mathrm{orb}}={f}_{\mathrm{GW}}/2$. If ${R}_{p}[{f}_{\mathrm{GW}}]\lt {a}_{\mathrm{circ}}[{f}_{\mathrm{orb}}]$, then the binary formed above f_GW. As a result, the LISA panel of this figure does not show the lines representing the latter two of these four categories as these mergers form with initial frequencies above 0.01 Hz. Similarly, we exclude mergers which form above 1 Hz in the DECIGO panel. However, in the LVK panel, we assign mergers that form above 10 Hz ${e}_{\mathrm{in}}=0.999$.

The top-left panel shows that the KL mechanism produces much more eccentric mergers in the LVK frequency band compared to in-cluster and ejected mergers, but also shows that the KL mergers are in general less eccentric than the few-body and single–single captures. All of the non-capture binary-mediated mergers have eccentricity less than ∼10⁻³. The non-KL triple merger eccentricity distribution reflects the fact that it is a subset of the standard ejected and in-cluster binary merger channels. On the other hand, $\sim 30 \%$ of the KL mergers have eccentricity above 10⁻² and a few percent are in excess of ∼0.1. This is consistent with previous results from Antonini et al. (2016). However, when compared to the few-body and GW capture scenarios, it is clear that the contribution of eccentric sources by KL mergers will be subdominant. In particular, almost $\sim 20 \%$ of the few-body mergers produced during resonant interactions were formed over ∼10 Hz.

When using a secular code to evaluate KL-driven evolution, there are many cases for which the secular equations of motion cannot accurately reproduce the evolution of the system during the peaks of the KL oscillation. In addition to producing more mergers with shorter merger times, direct N-body integration produces a much higher peak in the eccentricity spectrum than when using a secular code (Antonini et al. 2016; Grishin et al. 2018; Fragione et al. 2019a). This is because the inner binary angular momentum can become arbitrarily small during the peak of the KL oscillation, such that GW emission only begins to dominate the evolution once the inner binary frequency is near or within the LVK frequency band. Note that this would only affect the shape of the solid black curve showing the KL merging systems because the other populations do not experience the high eccentricity oscillations that would allow this to happen. We discuss this further in Section 4.6.

From the bottom panel, it is clear that all the KL mergers have high eccentricities at ${f}_{\mathrm{GW}}\sim 0.01\,\mathrm{Hz}$. More than $\sim 80 \%$ of these binaries have eccentricities in excess of ∼0.9. Remarkably, however, KL mergers will potentially enter and leave the LISA sensitivity range many times before merger, whereas all the other merger channels will evolve with monotonically increasing frequency. Randall & Xianyu (2019), Hoang et al. (2019), and Emami & Loeb (2020) have shown that it is possible for LISA to directly detect the eccentricity oscillations of hierarchical triples, whether the tertiary is a third stellar-mass body or an SMBH. ¹³ This is only possible when they are close to the peak of the eccentricity oscillation, when the rapid change in eccentricity results in a rapid change in the emitted peak harmonic frequency. Both works find that the typical timescale for changes in the characteristic strain will happen over the timescale of ${ \mathcal O }({10}^{2})$ days. Emami & Loeb (2020) showed that these methods can detect BH triples with stellar-mass components as far away as M87. While there could in principle be many Galactic triples undergoing KL oscillation, the population of triples detectable in this way will be dominated by triples in the Galactic field. We note that while the bottom panels show that $\sim 20 \%$ of non-KL triples have ${e}_{\mathrm{in}}\gt 0.99$, these sources will not show the same increase and subsequent decrease in eccentricity since the evolution is dominated by GW emission.

In order to have a better estimate for the number of eccentric sources and track their evolution prior to merger, it is necessary to probe the sub-hertz frequencies between the frequency ranges of LVK and LISA with future decihertz range detectors such as DECIGO (Kawamura et al. 2011). We find that $\sim 20 \%$ of the KL-driven mergers will have ${e}_{\mathrm{in}}\gt 0.1$ in the DECIGO band, comparable to the fraction of few-body captures and single–single captures that have ${e}_{\mathrm{in}}\gt 0.1$. Moreover, since some sources will form within the DECIGO frequency range, it will be possible to further distinguish the nature of different sources. However, as the absolute number of few-body capture and single–single capture mergers from a given cluster is higher than the number of KL mergers, we expect that the contribution of the KL mergers to the eccentric merger rate will be subdominant. We once again reiterate that the secular code underestimates the peak of the eccentricity spectrum, and thus we expect that a significant fraction of the KL-driven mergers will have high eccentricity in the DECIGO band.

These results are of particular interest in light of the recent detection of GW190521, which has component masses consistent with a dynamical origin (LIGO Scientific Collaboration & Virgo Collaboration 2020b; LIGO Scientific Collaboration & Virgo Collaboration 2020d). Romero-Shaw et al. (2020) have argued that the waveform could be consistent with a moderately eccentric binary at 10 Hz, depending on the priors used. This detection emphasizes the necessity of quantifying the contributions of the various dynamical channels to the population of eccentric merger events.

In Figure 5, we compare the total masses, chirp masses, and mass ratios of the merging binaries through both triple and binary channels, where the chirp mass is defined as

$$\begin{eqnarray}&&{m}_{\mathrm{chirp}}\equiv \displaystyle \frac{{\left({m}_{0}{m}_{1}\right)}^{3/5}}{{\left({m}_{0}+{m}_{1}\right)}^{1/5}}.\end{eqnarray} \tag{ 18 }$$

In the top panels, we compare the combined (KL and non-KL) triple sample to the binary mergers from Kremer et al. (2020b). We find a peak chirp mass of $\sim 17\,{M}_{\odot }$ and a peak total mass of $\sim 40\,{M}_{\odot }$ for the triple mergers. To first order, the mass and mass-ratio spectrum between binaries and triples are similar. Only one difference stands out: compared to the binary mergers, mergers in triples show a diminution of low-mass mergers with ${m}_{\mathrm{bin}}\lt 25\,{M}_{\odot }$ by approximately 40%. This could happen because low-mass binaries are less likely to take part in the binary–binary interactions necessary to produce the triples in our models. In the bottom panels, we compare the KL and non-KL systems to each other and the combined triple sample. Once again, we find minimal differences between these two distributions. This shows the influence of a tertiary companion has a minimal effect on the properties of merging binaries, whereas these properties are determined primarily by the far more frequent binary-mediated dynamical interactions within the cluster core.

Zoom In Zoom Out Reset image size

We find that our models produce heavier BH mergers when compared to Antonini et al. (2016), who did not produce any BH mergers with total masses above $\sim 50\,{M}_{\odot }$, though they did find similar peak total and chirp masses. This is most likely due to the wider parameter space for the initial conditions of our GC models. The difference is even more stark when compared to the study of field triples from Antonini et al. (2017), who were only able to produce binaries in the mass range of $13-20\,{M}_{\odot }$, emphasizing the key role of mass segregation in producing low-mass mergers. However, Rodriguez & Antonini (2018) were able to produce binaries in field triples with a total mass spectrum similar to the one presented here. The difference between this study and earlier studies arises due to a higher limit on the maximum stellar mass in their models.

In Figure 6, we present the component masses and the total mass versus mass ratio for KL mergers, non-KL mergers, and binaries from Kremer et al. (2020b) colored by metallicity. For comparison, we also present the properties of the confirmed BBH detections from GWTC-1 and the recent announcement of the detection of GW190412, the first binary BH system with asymmetric masses at high confidence (LIGO Scientific Collaboration & Virgo Collaboration 2019a, 2019b, 2020a). ¹⁴ We find that the three different populations present minimal differences, as expected from the previous discussion. Moreover, we find that the BBH mergers that have been announced so far are consistent with production in a GC, with all the GWTC-1 events that lie within dense regions of both diagrams. On the other hand, GCs are not as efficient at producing events in the mass/mass-ratio range as low as that of GW190412 due to the efficiency of mass segregation. More importantly, we find that merger via a hierarchical triple does not increase the percentage of similar events relative to the binary channel.

Zoom In Zoom Out Reset image size

However, recent works by Rodriguez et al. (2020) and Gerosa et al. (2020) have shown that the relative fraction of low-mass-ratio merger events can be enhanced by second- and third-generation mergers in dynamical environments, producing GW190412-like events as a result, though analysis by Kimball et al. (2020) casts doubt on the hierarchical merger scenario for GW190412. Due to the efficiency of mass segregation, the remnant of a previous merger is expected to be twice as massive as the 1G BHs in the cluster (see also Samsing & Hotokezaka 2020). In Section 4.5, we show how triples may contribute to increase the number of second-generation mergers.

As the initial distribution of BH spins is unknown, we adopt four different models of spin magnitude to bracket uncertainties. The different models are summarized in Table 2. For each model, we sample the spin of both inner binary components 10 times each, sampling the orientations isotropically. Models A, B, and C set the natal χ from a uniform distribution between 0 and 1, 0.5, and 0.2 respectively. Model D sets the birth spins using the following function which reflects work by Belczynski et al. (2020) set by the BH mass:

$$\begin{eqnarray}&&\chi =\displaystyle \frac{{p}_{1}-{p}_{2}}{2}\tanh \left({p}_{3}-\displaystyle \frac{m}{\,{M}_{\odot }}\right)+\displaystyle \frac{{p}_{1}+{p}_{2}}{2}.\end{eqnarray} \tag{ 19 }$$

The quantities in the previous formula have values ${p}_{1}\,=0.86\pm 0.06$, ${p}_{2}=0.13\pm 0.13$, and ${p}_{3}=29.5\pm 8.5$. Following Gerosa et al. (2018), we sample values in between the lower and upper limits of the parameters for a BH of a given mass m. This last model has the advantage of assigning lower spins to more massive BHs and vice versa.

Table 2.  χ Distributions of Different Models

Model	χ
A	Uniform in $[0,1)$
B	Uniform in $[0,0.5)$
C	Uniform in $[0,0.2)$
D	Equation (19)

Note. All models sample spin vector orientations isotropically.

Download table as:  ASCIITypeset image

We compare the ${\chi }_{\mathrm{eff}}$ distributions from different models in the top panel of Figure 7. All models peak around ${\chi }_{\mathrm{eff}}\approx 0$, which is to be expected from an isotropic distribution of spin–orbit misalignment angles. This is consistent with the announced detections from O1/O2, which were all consistent with having near zero effective spin. In the absence of GW190412, a model for natal spin such as B or C would be favored if all the events were generated by dynamically assembled binaries. However, it is extremely unlikely for models B and C to produce a ${\chi }_{\mathrm{eff}}$ consistent with that of GW190412. Additionally, as the detection of GW190412 allowed for a well-constrained measurement of ${\chi }_{{\rm{p}}}$, we compare these distributions for different spin models in the bottom panel. Models A and B are strongly favored as they peak within the 90% probability region, while model C is once again strongly disfavored. As such, if we assume that both components of GW190412 are 1G BHs, then a model that allows for higher natal spin is necessary.

Zoom In Zoom Out Reset image size

We find no significant difference between the initial and final ${\chi }_{\mathrm{eff}}$ spin distributions across any of our models, as shown in Figure 8. This is in agreement with previous findings from Rodriguez & Antonini (2018) and Fragione et al. (2019b). The spin vectors in the inner binary evolve in response to two different torques. On the one hand, they are precessing around the inner binary angular momentum according to the relativistic spin–orbit coupling (Equation (6)). On the other hand, the KL mechanism causes the precession of the inner orbit angular momentum around the outer orbit angular momentum. If the rate of precession from the relativistic spin–orbit coupling is faster than the rate of precession due to KL, then ${\chi }_{\mathrm{eff}}$ evolution is completely suppressed (Storch et al. 2014; Storch & Lai 2015; Liu & Lai 2017, 2018; Antonini et al. 2018). The individual triple systems span the full space of allowed values of initial versus final ${\chi }_{\mathrm{eff}}$, even if the initial and final distributions of ${\chi }_{\mathrm{eff}}$ do not appear to change (see Figure 8). As a result, we can safely conclude that spin relativistic precession does not suppress the evolution of ${\chi }_{\mathrm{eff}}$ in most cases. In summary, while spin precession does take place in these triples, since the triples are dynamically assembled, the final ${\chi }_{\mathrm{eff}}$ distribution will remain isotropic.

Zoom In Zoom Out Reset image size

Recent work by Belczynski et al. (2020) comparing different stellar evolution models in light of the detections of GWTC-1 has shown that the results from the LIGO–Virgo collaboration strongly favor low natal spins. However, more recent work by Olejak et al. (2020) has shown that GW190412 is consistent with an isolated binary formation scenario, though this claim is disputed by Safarzadeh & Hotokezaka (2020). It is not clear if isolated stellar evolutionary models alone can consistently produce the ${\chi }_{\mathrm{eff}}$ distribution in the detected events so far. Recent works by Rodriguez et al. (2020) and Gerosa et al. (2020) have also shown that it would be possible to reproduce these detections with repeated mergers of previous merger products in extremely large clusters. This would require that BHs are born with low natal spins, so as to retain more of the merger remnants. Further discussion of the importance of spin with respect to retaining merger remnants and producing mergers with 2G black holes will follow in Section 4.5.

Other formation channels lessen the sensitivity to the assumed natal spin model. In particular, previous work on field triples has shown that the spin evolution of these systems can take a wide range of evolutionary paths depending on initial assumptions. While the natal spin assumptions are important, initial spin–orbit misalignment and the properties of the outer orbit play as large a role in determining the final state of the system (Fragione et al. 2019b; Liu et al. 2019). Nevertheless, the field triple channel tends to produce distributions peaked at ${\chi }_{\mathrm{eff}}\simeq 0$ (Antonini et al. 2018; Rodriguez & Antonini 2018). Another proposed scenario involves a hierarchical 3 + 1 quadruple system, wherein the innermost orbit is driven to merger by secular chaotic evolution and then the merger remnant merges with another compact object due to the eKL mechanism (Safarzadeh et al. 2020). While the rates for these events are extremely uncertain, this is another potential formation channel for events like GW190412 (Hamers & Safarzadeh 2020). Similarly, 2 + 2 quadruples could also produce GW190412-like events (Fragione et al. 2020a). Another possibility is through the merger of BHs in AGN disks (Tagawa et al. 2018, 2020a, 2020b). In particular, recent work by Tagawa et al. (2020a) has shown that the spin distribution of the O1/O2 events is consistent with mergers in AGN disks, while a hierarchical merger within AGN disks could potentially reproduce the properties of GW190412.

Following the method of Rodriguez et al. (2015), the cumulative merger rate is given by

$$\begin{eqnarray}&&R(z)={\int }_{0}^{z}{ \mathcal R }(z^{\prime} )\displaystyle \frac{{{dV}}_{c}}{{dz}^{\prime} }{\left(1+z^{\prime} \right)}^{-1}{dz}^{\prime} .\end{eqnarray} \tag{ 20 }$$

Here, ${{dV}}_{c}/{dz}^{\prime}$ is the comoving volume at redshift z and ${ \mathcal R }(z^{\prime} )$ is the comoving source merger rate, where the comoving rate is given by

$$\begin{eqnarray}&&{ \mathcal R }=f\times {\rho }_{\mathrm{GC}}\times \displaystyle \frac{{dN}(z)}{{dt}}.\end{eqnarray} \tag{ 21 }$$

We assume that the volumetric number density of clusters has a constant value of ${\rho }_{\mathrm{GC}}=2.3\,{\mathrm{Mpc}}^{-3}$ (Rodriguez et al. 2015; Rodriguez & Loeb 2018). f is an r_v -dependent scaling factor to incorporate the high-end tail of the cluster mass function not covered in the catalog. Finally, ${dN}(z)/{dt}$ is the number of mergers per unit time at z.

In order to compute ${dN}(z)/{dt}$, we compile a complete list of ${T}_{\mathrm{mergers}}$ for each realization of the merging triples in our sample. For each merger, we resample their merger times by sampling 10 random host cluster ages for each merger and define the effective merger time ${t}_{\mathrm{effective}}={t}_{\mathrm{Hubble}}-{t}_{\mathrm{age}}+{T}_{\mathrm{merger}}$ using the host cluster ages from the metallicity-dependent age distributions of El-Badry et al. (2019). We then bin this list of merger times into redshift bins. We must also scale down these rates to correct for oversampling. We divide the rates by a factor of 1000 to account for the resampling with respect to cluster age, triple recoil velocity, and orbital orientation. We must also divide by the total number of cluster models in order to correct for drawing from a large set of cluster models, weighting all models equally for simplicity.

The scaling factor f is included in order to account for the low-mass tail of the cluster mass function not covered by our models (consistent with Kremer et al. 2020b). In practice, high-mass clusters ($M\gtrsim 5\times {10}^{5}\,\,{M}_{\odot }$) contribute roughly four times the number of mergers compared to low-mass clusters ($M\lesssim 5\times {10}^{5}\,\,{M}_{\odot }$), so $f\approx 4$ across different r_v values. Full details can be found in Section 9.2 and Table 5 of Kremer et al. (2020b).

In Figure 9, we compare the cumulative and comoving merger rates, finding that the triple merger rate remains roughly two orders of magnitude less than the binary merger rate from GCs across all redshifts. Thus, we find a local-universe eKL-assisted triple merger rate in GCs $\approx 0.35$ Gpc⁻³ yr⁻¹. This is consistent with the previous merger rate found by Antonini et al. (2016). The non-KL merger rate on the other hand is $\approx 0.62$ Gpc⁻³ yr⁻¹, with a possible range of ≈0.2–2 Gpc⁻³ yr⁻¹. We also show rates for events that may have a detectable eccentricity in the LVK band at design sensitivity. ¹⁵ In agreement with the lack of eccentric events thus far (LIGO Scientific Collaboration & Virgo Collaboration 2019c), we find low rates $\approx 0.025\mbox{--}0.035$ Gpc⁻³ yr⁻¹ for the triple channel and $\approx 1$ Gpc⁻³ yr⁻¹ for the few-body and single–single capture channels. This should be compared to the following merger rate estimates from the LIGO–Virgo collaboration and from studies looking at hierarchical triples in the field:

• 1.  
  LVC rate: 9.7–101 Gpc⁻³ yr⁻¹ (LIGO Scientific Collaboration & Virgo Collaboration 2019a, 2019b).
• 2.  
  Field triples: 0.14–6 Gpc⁻³ yr⁻¹ (Silsbee & Tremaine 2017); 0.3–1.3 Gpc⁻³ yr⁻¹ (Antonini et al. 2017); 2–25 Gpc⁻³ yr⁻¹ (Rodriguez & Antonini 2018); 0.02–24 Gpc⁻³ yr⁻¹ (Fragione et al. 2019b).
• 3.  
  Nuclear star cluster triples: ≈1–10² detections yr⁻¹ (O'Leary et al. 2009); 0.6–15 Gpc⁻³ yr⁻¹ (Petrovich & Antonini 2017); 1–3 Gpc⁻³ yr⁻¹ (Hoang et al. 2018; Fragione et al. 2019a; Stephan et al. 2019).

This rate should be taken as a conservative estimate, as this study only focuses on triples composed of three BHs. We do not analyze the other classes of triples, which may well become BBHs with other companions, as we are unable to self-consistently track the stellar evolution of three separate bodies undergoing the eKL mechanism (see, e.g., Di Stefano 2020a, 2020b).

Zoom In Zoom Out Reset image size

While it is clear that KL mergers are an addition to the previous merger rate found by Kremer et al. (2020b), it is not clear if all non-KL mergers were included in this previous estimate. In principle, since all of these binaries merged solely due to GW radiation, regardless of whether or not the tertiary companion were present, these events have already been accounted for in the previous binary rate estimate. In Section 4.5, we will explain why it is still important to keep track of triple systems whether or not the eKL mechanism dominates the triple evolution.

When two BHs merge, the merger remnant will experience a kick due to the asymmetric emission of GW radiation, depending on the mass asymmetry and the alignment of the spin vectors with the binary orbital plane (Campanelli et al. 2007; Lousto & Zlochower 2008, 2013; Lousto et al. 2010, 2012). The contribution due to the mass asymmetry takes the form

$$\begin{eqnarray}&&{v}_{\mathrm{kick}}=A{\eta }^{2}\displaystyle \frac{1-q}{1+q}(1+B\eta ),\end{eqnarray} \tag{ 22 }$$

where $\eta ={m}_{0}{m}_{1}/{\left({m}_{0}+m1\right)}^{2}$, q is the mass ratio, and A and B are numerically calibrated constants, with $A=1.2\,\times {10}^{4}$ km s⁻¹ and $B=-0.93$. The numerical relativity simulations from those same studies have shown that merging BHs with even moderate spin greatly increase the magnitude of the kick velocity, depending on their mutual misalignment and the misalignment with the orbital plane, reaching magnitudes of up to thousands of km s⁻¹. Because all the models of the CMC Catalog assume that all first -generation (1G) BHs are born with no spin (Fuller & Ma 2019), we elect to only consider the mass asymmetry term in the recoil velocity.

If the kick velocity is greater than the local escape speed, then the remnant will be ejected from the cluster. However, if the merger happens within a triple, then the remnant must also overcome the potential of the tertiary body. As a result, the merger remnant must overcome the local escape speed, which we calculate as follows:

$$\begin{eqnarray}&&{v}_{\mathrm{esc},\mathrm{eff}}={v}_{\mathrm{outer}}+{v}_{\mathrm{esc}}-\alpha (r){v}_{\mathrm{rec}}.\end{eqnarray} \tag{ 23 }$$

Here, ${v}_{\mathrm{esc}}$ is the escape speed from the cluster core at the time of merger, ${v}_{\mathrm{outer}}$ is the escape speed from the tertiary companion, and ${v}_{\mathrm{rec}}$ is the recoil kick that the triple experiences at formation. Due to the uncertainty regarding the triple's location within the cluster at the time of merger, we include a factor $\alpha (r)\in [0,1]$ to account for the phase of the orbit and the effect of dynamical friction in reducing the orbit's apocenter (Binney & Tremaine 2008). We define this parameter such that $\alpha =1$ if the triple is at the apocenter of its orbit and has experienced no dynamical friction and $\alpha =0$ if the triple is located within the cluster core.

In Figure 10, we compare the retention fraction for both KL and non-KL systems including and not including the potential of the tertiary. Since we do not know the location of the triple within the cluster at the time of the merger, we show the range of possibilities by showing curves where $\alpha =0$ and $\alpha =1$. In reality, the merger retention fraction will be somewhere in between these values. We find that including the tertiary's potential has an appreciable effect on merger remnant retention. We find that for KL mergers, the fraction of retained systems increases by $\sim 10 \%$, whereas the retention fraction increases by $\sim 6 \%$ for non-KL mergers. The effect of the tertiary is much more pronounced for KL-induced mergers since the orbits are more compact. We note that the solid green line in the bottom panel is equivalent to the merger remnant retention fraction for the binary merger channel, as it corresponds to mergers in the cluster core without a tertiary companion. From this, it is clear that triple systems enhance the number of 2G BHs retained in the host cluster.

Zoom In Zoom Out Reset image size

These results are extremely interesting in light of the recent detection of GW190521 (LIGO Scientific Collaboration & Virgo Collaboration 2020d), a merger event with a total mass of $150\,{M}_{\odot }$, whose components could be 2G BHs. Many studies have discussed the frequency and importance of successive mergers with 2G and higher BHs in massive clusters such as GCs, super star clusters, and nuclear star clusters (Antonini et al. 2019; Rodriguez et al. 2019, 2020; Fragione & Silk 2020; Gerosa et al. 2020; Samsing & Hotokezaka 2020). Mergers involving 2G BHs are expected to be more massive, potentially falling within the mass range of $\sim 46\mbox{--}133\,\,{M}_{\odot }$, the so-called upper mass gap, where no BHs are expected to be formed due to pair instability supernovae (Woosley & Heger 2015; Woosley 2017; Wen 2019). As a result, if any BH mergers are detected with a component in the upper mass gap, the only way they could have been formed is through prior stellar mergers or successive mergers of BHs (Di Carlo et al. 2020; Fragione & Silk 2020; Kremer et al. 2020a). This possibility has been recently confirmed with the detection of GW190521, which is consistent with a dynamical origin (LIGO Scientific Collaboration & Virgo Collaboration 2020c). Second, it is expected that a 2G+1G merger will have very asymmetric masses. Mass segregation in globular and higher-mass clusters causes binaries to form primarily between two low-mass members. Since it is unlikely that there will be more than one 2G BH in a cluster at any time, the 2G BH will necessarily form with a 1G BH of much lower mass (Rodriguez et al. 2019). Rodriguez et al. (2019) have shown that when assuming BHs have no natal spin 2G+1G mergers typically contribute $\sim 20 \%$ of detectable mergers from GCs, while 2G+2G mergers only contribute a few percent. Thus, keeping track of successive mergers is important to understand the true mass and mass-ratio spectrum from GCs. It is also worth noting that since the LVK detectors are more sensitive to more massive BHs, mergers involving 2G or higher BHs will make up a larger fraction of the detected BHs.

In addition to the mass gap, Baibhav et al. (2020) described the existence of the "spin gap" in the context of successive mergers. If two BHs are non-spinning at birth, the only way to form a BH where one or both components have high spin is through successive mergers. Even if BHs are born with some natal spin, studies by Berti & Volonteri (2008), Gerosa & Berti (2017), and Fishbach et al. (2017) have shown that, regardless of the model of the natal spin of 1G BHs, the spin distribution of mergers involving 2G BHs peaks at $\chi \simeq 0.7$. Successive mergers would be the only way to form very high ${\chi }_{\mathrm{eff}}$ mergers, requiring the dynamical assembly of binaries. Constraining the rate of higher-generation mergers using observations populating the upper mass gap and the spin gap will be important to disentangle the isolated field formation channel from the dynamical formation channels. Since 2G+1G and 2G+2G events are not expected to make up a large fraction of detections, it is crucial to accurately constrain the rates of these events.

At the most extreme eccentricity maxima of the eKL oscillations, the system may enter a dynamical regime wherein the secular approximation breaks down (Antonini & Perets 2012; Antonini et al. 2016; Fragione et al. 2019a). Furthermore, some systems close to the stability limit may also have angular momentum changes faster than the orbital period of the inner binary (Antonini et al. 2014). In this regime, the secular equations of motion could no longer describe the system correctly, thus necessitating a direct n-body integration of the orbits. Antonini et al. (2016) and Fragione et al. (2019a) have demonstrated that using an n-body code instead of a secular code for the same sample increased the total number of mergers and increased the fraction of merging systems entering the LVK frequency band with a finite eccentricity. We confirm these results by simulating a small number (∼100) of the merging systems using a direct integration with ARCHAIN ¹⁶ (Mikkola & Merritt 2006, 2008; Chassonnery et al. 2019).

In Figure 11, we show an example of the dynamical evolution of one such system until merger. The system shown on the left is semi-hierarchical with ${a}_{\mathrm{out}}(1-{e}_{\mathrm{out}})/{a}_{\mathrm{in}}\approx 12.7$ and the value of the octupole parameter is $\epsilon =0.0034$. Since the system is semi-hierarchical, during certain points in the eKL oscillation, the change in angular momentum of the inner binary happens on a timescale faster than the outer orbital period. This translates to the following numerical condition on the separation between the inner and outer binary (Equation (9) of Antonini et al. 2016):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{a}_{\mathrm{out}}(1-{e}_{\mathrm{out}})}{{a}_{\mathrm{in}}}\lesssim \displaystyle \frac{3}{1+{e}_{\mathrm{out}}}{\left(\displaystyle \frac{{m}_{2}}{{m}_{\mathrm{bin}}}\right)}^{2/3}\\ \quad \times \,{\left(\displaystyle \frac{{m}_{\mathrm{bin}}}{{m}_{\mathrm{trip}}}\right)}^{1/3}{\left(\displaystyle \frac{{a}_{\mathrm{in}}}{{10}^{9}\mathrm{cm}}\right)}^{1/3}.\end{array}\end{eqnarray} \tag{ 24 }$$

In this case, the condition is satisfied. As a result, the secular integration will not be able to reproduce the most extreme peaks of the eccentricity oscillations. In both the secular and the n-body integration, with each peak in the eccentricity oscillation, the semimajor axis decreases by a small but appreciable amount. However, while the two different numerical methods show agreement during the first two eKL oscillations, the center panel shows that the semimajor axis begins decreasing more rapidly in the n-body integration, implying a larger peak in eccentricity compared to the secular integration. Each successive eccentricity oscillation accelerates the divergence of the results between these two methods. While the secular integration shows a gradual decline in ${a}_{\mathrm{in}}$ as GW radiation begins to dominate, the direct n-body integration shows the inner binary merges in about half the time compared to the secular integration. Furthermore, we find direct integration indicates an eccentricity of $\sim 2\times {10}^{-3}$ at $10\,\mathrm{Hz}$ versus $\sim 7\times {10}^{-4}$ from the secular integration.

Zoom In Zoom Out Reset image size

On the other hand, the system shown on the right is more hierarchical, with ${a}_{\mathrm{out}}(1-{e}_{\mathrm{out}})/{a}_{\mathrm{in}}\approx 30.5$ and Equation (24) is not satisfied. As such, there is much better agreement between the two numerical schemes.

To summarize, using an n-body integration method instead of a secular integration method increases the total number of mergers and the fraction of merging systems entering the LVK frequency band with a high eccentricity (e.g., Antonini & Perets 2012; Antonini et al. 2016; Fragione et al. 2019a), which we have confirmed by simulating a small sample of our triple systems. In light of this and also due to the large computational cost from the number of systems requiring integration, we choose to primarily restrict ourselves to using a secular code. Finally, we note that there some issues may arise with incorporating the conservative post-Newtonian effects into n-body simulations, since the Newtonian energy is no longer explicitly conserved in this approximation (see, e.g., Figure 2 of Rodriguez et al. 2018). As a result, the incorporation of the post-Newtonian effects could introduce an increase in energy during the pericenter passage. In hierarchical triples, this could cause the triple to move further into the chaotic regime and increase the amount of GW emission, possibly changing the inspiral waveform.