Black Hole and Neutron Star Binary Mergers in Triple Systems: Merger Fraction and Spin–Orbit Misalignment

If we introduce the instantaneous separation between the inner bodies as ${\boldsymbol{r}}\equiv r\hat{{\boldsymbol{r}}}$, and the separation between the external perturber and the center of mass of the inner bodies as ${{\boldsymbol{r}}}_{{\rm{out}}}\equiv {r}_{{\rm{out}}}{\hat{{\boldsymbol{r}}}}_{{\rm{out}}}$, then the complete Hamiltonian of the system can be written as (e.g., Harrington 1968)

$$\begin{eqnarray}&&{ \mathcal H }=\displaystyle \frac{1}{2}\mu | \dot{{\boldsymbol{r}}}{| }^{2}+\displaystyle \frac{1}{2}{\mu }_{\mathrm{out}}| {\dot{{\boldsymbol{r}}}}_{\mathrm{out}}{| }^{2}-\displaystyle \frac{{{Gm}}_{1}{m}_{2}}{r}-\displaystyle \frac{{{Gm}}_{12}{m}_{3}}{{r}_{\mathrm{out}}}+{\rm{\Phi }},\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}{\rm{\Phi }} & = & -{{Gm}}_{1}{m}_{2}{m}_{3}\displaystyle \sum _{l\,=\,2}^{\infty }\left[\displaystyle \frac{{m}_{1}^{l-1}+{(-1)}^{l}{m}_{2}^{l-1}}{{m}_{12}^{l}}\right]\\ & & \times \ \displaystyle \frac{{r}^{l}}{{r}_{\mathrm{out}}^{l\,+\,1}}{P}_{l}(\cos \,\theta ).\end{eqnarray} \tag{ 13 }$$

Here P_l(x) is the Legendre polynomial of degree l and θ is the angle between ${\boldsymbol{r}}$ and ${{\boldsymbol{r}}}_{\mathrm{out}}$.

2.1.1. Double-averaged Secular Equations

For sufficiently hierarchical systems, the angular momenta of the inner and outer binaries exchange periodically over a long timescale (longer than the companion's orbital period), while the exchange of energy is negligible. The orbital evolution of the triple system can be studied by expanding the Hamiltonian to the octupole order and averaging over both the inner and outer orbits (double averaging), i.e., Φ = Φ_quad + Φ_oct. The quadrupole (l = 2) piece is given by

$$\begin{eqnarray}\langle \langle {{\rm{\Phi }}}_{\mathrm{quad}}\rangle \rangle & = & \displaystyle \frac{\mu {{\rm{\Phi }}}_{0}}{8}[1-6{e}^{2}-3(1-{e}^{2}){(\hat{{\boldsymbol{L}}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}})}^{2}\\ & & +\ 15{e}^{2}{(\hat{{\boldsymbol{e}}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}})}^{2}],\end{eqnarray} \tag{ 14 }$$

and the octupole (l = 3) potential is

$$\begin{eqnarray}\langle \langle {{\rm{\Phi }}}_{\mathrm{oct}}\rangle \rangle & = & \displaystyle \frac{15\mu {{\rm{\Phi }}}_{0}{\varepsilon }_{\mathrm{oct}}}{64}\{e(\hat{{\boldsymbol{e}}}\cdot {\hat{{\boldsymbol{e}}}}_{\mathrm{out}})[8{e}^{2}-1\\ & & -\ 35{e}^{2}{(\hat{{\boldsymbol{e}}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}})}^{2}+5(1-{e}^{2}){(\hat{{\boldsymbol{L}}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}})}^{2}]\\ & & +\ 10e(1-{e}^{2})(\hat{{\boldsymbol{e}}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}})(\hat{{\boldsymbol{L}}}\cdot {\hat{{\boldsymbol{e}}}}_{\mathrm{out}})(\hat{{\boldsymbol{L}}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}})\},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{0}\equiv \displaystyle \frac{{{Gm}}_{3}{a}^{2}}{{a}_{\mathrm{out}}^{3}{(1-{e}_{\mathrm{out}}^{2})}^{3/2}},\end{eqnarray} \tag{ 16 }$$

and

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{oct}}\equiv \displaystyle \frac{{m}_{1}-{m}_{2}}{{m}_{12}}\left(\displaystyle \frac{a}{{a}_{\mathrm{out}}}\right)\displaystyle \frac{{e}_{\mathrm{out}}}{1-{e}_{\mathrm{out}}^{2}}.\end{eqnarray} \tag{ 17 }$$

The explicit expressions for ${(d{\boldsymbol{L}}/{dt})}_{\mathrm{LK}}$, (de/dt)_LK and for (dL_out/dt)_LK, (de_out/dt)_LK are provided in Liu et al. (2015). In general, ${\dot{{\boldsymbol{L}}}}_{\mathrm{in},\mathrm{out}}$ and ${\dot{{\boldsymbol{e}}}}_{\mathrm{in},\mathrm{out}}$ consist of two pieces: a quadrupole term and an octupole term. The quadrupole term induces the oscillations in the eccentricity and mutual orbital inclination on the timescale of

$$\begin{eqnarray}&&{t}_{\mathrm{LK}}=\displaystyle \frac{1}{n}\displaystyle \frac{{m}_{12}}{{m}_{3}}{\left(\displaystyle \frac{{a}_{\mathrm{out},\mathrm{eff}}}{a}\right)}^{3},\end{eqnarray} \tag{ 18 }$$

where the effective outer binary separation is defined as

$$\begin{eqnarray}&&{a}_{\mathrm{out},\mathrm{eff}}\equiv {a}_{\mathrm{out}}\sqrt{1-{e}_{\mathrm{out}}^{2}}.\end{eqnarray} \tag{ 19 }$$

The octupole piece is quantified by terms proportional to ε_oct, which measures the strength of the octupole potential relative to the quadrupole one.

For systems that can be correctly described by the DA equations, the timescale for eccentricity variation of the inner binary must be longer than the period of the companion's orbit. Otherwise, the secular equations may break down (e.g., Seto 2013; Antonini et al. 2014). Note that when the eccentricity of the inner binary is excited to the maximum value e_max, the eccentricity vector  e evolves on a timescale ${t}_{\mathrm{LK}}\sqrt{1-{e}_{\max }^{2}}$ (e.g., Anderson et al. 2016), much shorter than the quadrupole LK period (∼t_LK). Thus, for the DA secular equations to be valid, we require

$$\begin{eqnarray}&&{t}_{\mathrm{LK}}\sqrt{1-{e}_{\max }^{2}}\gtrsim {P}_{\mathrm{out}},\end{eqnarray} \tag{ 20 }$$

where P_out is the period of the outer binary.

2.1.2. Single-averaged Secular Equations

For moderately hierarchical systems, the change in the angular momentum of the inner binary may be significant within one period of the outer orbit, and the short-term (≲P_out) oscillations of the system cannot be ignored. In this case, the DA secular equations break down, and we can use the SA secular equations (only averaging over the inner orbital period).

Averaging over the inner orbit, the quadrupole term in Equation (13) becomes

$$\begin{eqnarray}&&\langle {{\rm{\Phi }}}_{\mathrm{quad}}\rangle =\displaystyle \frac{\mu {\rm{\Phi }}{{\prime} }_{0}}{4}[-1+6{e}^{2}+3{({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})}^{2}-15{({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})}^{2}],\end{eqnarray} \tag{ 21 }$$

and the octupole term is

$$\begin{eqnarray}\langle {{\rm{\Phi }}}_{\mathrm{oct}}\rangle & = & \displaystyle \frac{5\mu {\rm{\Phi }}{{\prime} }_{0}\varepsilon {{\prime} }_{\mathrm{oct}}}{16}[(3-24{e}^{2})({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})\\ & & -\ 15{({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})}^{2}({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})+35{({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})}^{3}],\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{j}}\equiv j\hat{{\boldsymbol{L}}}=\sqrt{1-{e}^{2}}\hat{{\boldsymbol{L}}}\end{eqnarray} \tag{ 23 }$$

is the dimensionless angular momentum vector, and the coefficients Φ'₀ and ε'_oct are given by

$$\begin{eqnarray}&&{\rm{\Phi }}{{\prime} }_{0}=\displaystyle \frac{{{Gm}}_{3}{a}^{2}}{{r}_{\mathrm{out}}^{3}},\end{eqnarray} \tag{ 24 }$$

and

$$\begin{eqnarray}&&\varepsilon {{\prime} }_{\mathrm{oct}}=\displaystyle \frac{{m}_{1}-{m}_{2}}{{m}_{12}}\displaystyle \frac{a}{{r}_{\mathrm{out}}}.\end{eqnarray} \tag{ 25 }$$

In terms of the averaged potentials, the equations of motion for the inner orbital vectors  j and  e are (e.g., Tremaine et al. 2009)

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{j}}}{{dt}}=-\displaystyle \frac{1}{\mu \sqrt{{{Gm}}_{12}a}}({\boldsymbol{j}}\times {{\rm{\nabla }}}_{{\boldsymbol{j}}}\langle {\rm{\Phi }}\rangle +{\boldsymbol{e}}\times {{\rm{\nabla }}}_{{\boldsymbol{e}}}\langle {\rm{\Phi }}\rangle ),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{e}}}{{dt}}=-\displaystyle \frac{1}{\mu \sqrt{{{Gm}}_{12}a}}({\boldsymbol{j}}\times {{\rm{\nabla }}}_{{\boldsymbol{e}}}\langle {\rm{\Phi }}\rangle +{\boldsymbol{e}}\times {{\rm{\nabla }}}_{{\boldsymbol{j}}}\langle {\rm{\Phi }}\rangle ).\end{eqnarray} \tag{ 27 }$$

Substituting Equation (21) into (26) and (27), the quadrupole-level equations can be obtained:

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d{\boldsymbol{j}}}{{dt}}\right|}_{\mathrm{quad}}=\displaystyle \frac{3}{2{t}_{\mathrm{LK}}^{{\prime} }}[5({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}){\boldsymbol{e}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}-({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}){\boldsymbol{j}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}],\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}{\left.\displaystyle \frac{d{\boldsymbol{e}}}{{dt}}\right|}_{\mathrm{quad}} & = & \displaystyle \frac{3}{2{t}_{\mathrm{LK}}^{{\prime} }}[5({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}){\boldsymbol{j}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}\\ & & -\ ({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}){\boldsymbol{e}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}-2{\boldsymbol{j}}\times {\boldsymbol{e}}].\end{eqnarray} \tag{ 29 }$$

Similarly, the octupole contributions are

$$\begin{eqnarray}{\left.\displaystyle \frac{d{\boldsymbol{j}}}{{dt}}\right|}_{\mathrm{oct}} & = & \displaystyle \frac{15\varepsilon {{\prime} }_{\mathrm{oct}}}{16{t}_{\mathrm{LK}}^{{\prime} }}[10({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}){\boldsymbol{j}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}\\ & & -\ (1-8{e}^{2}){\boldsymbol{e}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}+5{({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})}^{2}{\boldsymbol{e}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}\\ & & -\ 35{({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})}^{2}{\boldsymbol{e}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}],\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}{\left.\displaystyle \frac{d{\boldsymbol{e}}}{{dt}}\right|}_{\mathrm{oct}} & = & \displaystyle \frac{15\varepsilon {{\prime} }_{\mathrm{oct}}}{16{t}_{\mathrm{LK}}^{{\prime} }}[16({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}){\boldsymbol{j}}\times \hat{{\boldsymbol{e}}}-(1-8{e}^{2}){\boldsymbol{j}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}\\ & & +\ 5{({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})}^{2}{\boldsymbol{j}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}-35{({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})}^{2}{\boldsymbol{j}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}\\ & & +\ 10({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}})({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}){\boldsymbol{e}}\times {\hat{{\boldsymbol{r}}}}_{\mathrm{out}}].\end{eqnarray} \tag{ 31 }$$

In the above, we have defined the SA (quadrupole) LK timescale as

$$\begin{eqnarray}&&{t}_{\mathrm{LK}}^{{\prime} }=\displaystyle \frac{1}{n}\displaystyle \frac{{m}_{12}}{{m}_{3}}{\left(\displaystyle \frac{{r}_{\mathrm{out}}}{a}\right)}^{3}.\end{eqnarray} \tag{ 32 }$$

The evolution equations of  L and  e are

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d{\boldsymbol{L}}}{{dt}}\right|}_{\mathrm{LK}}=\mu \sqrt{{{Gm}}_{12}a}\left({\left.\displaystyle \frac{d{\boldsymbol{j}}}{{dt}}\right|}_{\mathrm{quad}}+{\left.\displaystyle \frac{d{\boldsymbol{j}}}{{dt}}\right|}_{\mathrm{oct}}\right),\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d{\boldsymbol{e}}}{{dt}}\right|}_{\mathrm{LK}}={\left.\displaystyle \frac{d{\boldsymbol{e}}}{{dt}}\right|}_{\mathrm{quad}}+{\left.\displaystyle \frac{d{\boldsymbol{e}}}{{dt}}\right|}_{\mathrm{oct}}.\end{eqnarray} \tag{ 34 }$$

On the other hand, for the external companion, the dynamics is governed by

$$\begin{eqnarray}&&{\mu }_{\mathrm{out}}\displaystyle \frac{{d}^{2}{{\boldsymbol{r}}}_{\mathrm{out}}}{{{dt}}^{2}}={{\rm{\nabla }}}_{{{\boldsymbol{r}}}_{\mathrm{out}}}\left(\displaystyle \frac{{{Gm}}_{12}{m}_{3}}{{r}_{\mathrm{out}}}\right)-{{\rm{\nabla }}}_{{{\boldsymbol{r}}}_{\mathrm{out}}}(\langle {{\rm{\Phi }}}_{\mathrm{quad}}\rangle +\langle {{\rm{\Phi }}}_{\mathrm{oct}}\rangle ).\end{eqnarray} \tag{ 35 }$$

The explicit form can be obtained by substituting Equations (21) and (22) into (35). Equations (28)–(31) and (33)–(35), together with Equations (6)–(9), completely determine the evolution of the triple system in the single averaging approximation.

The SA equations are applicable to a wider range of system parameters than the DA equations. Nevertheless, their validity still requires that the timescale for eccentricity evolution at e ∼ e_max be longer than the orbital period of the inner binary, i.e.,

$$\begin{eqnarray}&&{t}_{\mathrm{LK}}\sqrt{1-{e}_{\max }^{2}}\gtrsim {P}_{\mathrm{in}}.\end{eqnarray} \tag{ 36 }$$

Before exploring the LK-induced mergers systematically (Section 3), we summarize some key analytical results for LK eccentricity excitations. It is well known that the effects of short-range forces (such as GR-induced apsidal precession; see Equation (6)) play an important role in determining the maximum eccentricity e_max in LK oscillations (e.g., Holman et al. 1997; Fabrycky & Tremaine 2007). Analytical expression for e_max for general hierarchical triples (arbitrary masses) can be obtained in the DA secular approximation when the disturbing potential is truncated to the quadrupole order (Liu et al. 2015; Anderson et al. 2016, 2017).

In the absence of energy dissipation, the evolution of the triple is governed by two conservation laws. The first is the total orbital angular momentum of the system,  L_tot =  L +  L_out. In the quadrupole approximation, e_out is constant, and the conservation of $| {{\boldsymbol{L}}}_{\mathrm{tot}}|$ implies

$$\begin{eqnarray}&&K\equiv j\cos I-\displaystyle \frac{\eta }{2}{e}^{2}=\mathrm{constant},\end{eqnarray} \tag{ 37 }$$

where $j=\sqrt{1-{e}^{2}}$, I is the angle between $\hat{{\boldsymbol{L}}}$ and ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$, and we have defined

$$\begin{eqnarray}&&\eta \equiv {\left(\displaystyle \frac{L}{{L}_{\mathrm{out}}}\right)}_{e=0}=\displaystyle \frac{\mu }{{\mu }_{\mathrm{out}}}{\left[\displaystyle \frac{{m}_{12}a}{{m}_{123}{a}_{\mathrm{out}}(1-{e}_{\mathrm{out}}^{2})}\right]}^{1/2}.\end{eqnarray} \tag{ 38 }$$

In the limit L ≪ L_out, Equation (37) reduces to the usual "Kozai constant," $\sqrt{1-{e}^{2}}\cos I\,=$ constant.

The second conserved quantity is the total energy. In the double averaging approximation, it is given by (to the quadrupole order)

$$\begin{eqnarray}&&{\rm{\Phi }}=\langle \langle {{\rm{\Phi }}}_{\mathrm{quad}}\rangle \rangle +\langle {{\rm{\Phi }}}_{\mathrm{GR}}\rangle ,\end{eqnarray} \tag{ 39 }$$

where $\langle \langle {{\rm{\Phi }}}_{\mathrm{quad}}\rangle \rangle$ is given by Equation (14), and $\langle {{\rm{\Phi }}}_{\mathrm{GR}}\rangle$ is given by

$$\begin{eqnarray}&&\langle {{\rm{\Phi }}}_{\mathrm{GR}}\rangle =-\displaystyle \frac{3{G}^{2}\mu {m}_{12}^{2}}{{a}^{2}{c}^{2}j}=-{\varepsilon }_{\mathrm{GR}}\displaystyle \frac{\mu {{\rm{\Phi }}}_{0}}{j},\end{eqnarray} \tag{ 40 }$$

with

$$\begin{eqnarray}{\varepsilon }_{\mathrm{GR}} & = & \displaystyle \frac{3{{Gm}}_{12}^{2}{a}_{\mathrm{out},\mathrm{eff}}^{3}}{{c}^{2}{a}^{4}{m}_{3}}\\ & & \simeq \ 3.6\,\times \,{10}^{-5}{\left(\displaystyle \frac{{m}_{12}}{60{M}_{\odot }}\right)}^{2}{\left(\displaystyle \frac{{m}_{3}}{30{M}_{\odot }}\right)}^{-1}\\ & & \times \ {\left(\displaystyle \frac{{a}_{\mathrm{out},\mathrm{eff}}}{{10}^{3}\mathrm{au}}\right)}^{3}{\left(\displaystyle \frac{a}{{10}^{2}\mathrm{au}}\right)}^{-4}.\end{eqnarray} \tag{ 41 }$$

Using Equations (37) and (39), the maximum eccentricity e_max attained in the LK oscillation (starting from an initial I₀ and e₀ ≃ 0) can be calculated analytically (Liu et al. 2015; Anderson et al. 2017):

$$\begin{eqnarray}&&\displaystyle \frac{3}{8}\displaystyle \frac{{j}_{\min }^{2}-1}{{j}_{\min }^{2}}\left[5{\left(\cos {I}_{0}+\displaystyle \frac{\eta }{2}\right)}^{2}\right.\\ &&\,-\ \left.\left(3+4\eta \cos {I}_{0}+\displaystyle \frac{9}{4}{\eta }^{2}\right){j}_{\min }^{2}+{\eta }^{2}{j}_{\min }^{4}\right]\\ &&\,+\ {\varepsilon }_{\mathrm{GR}}(1-{j}_{\min }^{-1})=0,\end{eqnarray} \tag{ 42 }$$

where ${j}_{\min }\,\equiv \,\sqrt{1-{e}_{\max }^{2}}$. Note that in the limit $\eta \to 0$ and ${\varepsilon }_{\mathrm{GR}}\to 0$, Equation (42) yields the well-known relation ${e}_{\max }=\sqrt{1-(5/3){\cos }^{2}{I}_{0}}$. For general η, the maximum possible e_max for all values of I₀, called e_lim, is achieved at I_0,lim that satisfies de_max/dI₀ = 0, i.e.,

$$\begin{eqnarray}&&\cos {I}_{0,\mathrm{lim}}=\displaystyle \frac{\eta }{2}\left(\displaystyle \frac{4}{5}{j}_{\mathrm{lim}}^{2}-1\right).\end{eqnarray} \tag{ 43 }$$

Substituting Equation (43) into (42), we find that the limiting eccentricity e_lim, the maximum of the e_max(I₀) curve, is determined by

$$\begin{eqnarray}&&\displaystyle \frac{3}{8}({j}_{\mathrm{lim}}^{2}-1)\left[-3+\displaystyle \frac{{\eta }^{2}}{4}\left(\displaystyle \frac{4}{5}{j}_{\mathrm{lim}}^{2}-1\right)\right]+{\varepsilon }_{\mathrm{GR}}(1-{j}_{\mathrm{lim}}^{-1})=0.\end{eqnarray} \tag{ 44 }$$

On the other hand, eccentricity excitation (e_max ≥ 0) occurs only within a window of inclinations (cos I₀)₋ ≤ cos I₀ ≤ (cos I₀)₊, where (Anderson et al. 2017)

$$\begin{eqnarray}&&{(\cos {I}_{0})}_{\pm }=\displaystyle \frac{1}{10}\left(-\eta \pm \sqrt{{\eta }^{2}+60-\displaystyle \frac{80}{3}{\varepsilon }_{\mathrm{GR}}}\right).\end{eqnarray} \tag{ 45 }$$

This window exists only when

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{GR}}\leqslant \displaystyle \frac{9}{4}+\displaystyle \frac{3}{80}{\eta }^{2}.\end{eqnarray} \tag{ 46 }$$

In other words, no eccentricity excitation is possible when relation (46) is not satisfied.

Figure 1 shows some examples of the e_max(I₀) curves. For η ≲ 1, these curves depend mainly on ${m}_{3}/{a}_{\mathrm{out},\mathrm{eff}}^{3}$ (for given inner binary parameters). We see that the excitation of eccentricity can only happen within a finite range of I₀, and the achieved maximum e cannot exceed e_lim for any values of η.

Zoom In Zoom Out Reset image size

For systems with ${m}_{1}\ne {m}_{2}$ and ${e}_{\mathrm{out}}\ne 0$, ε_oct is non-negligible, the octupole effect may become important (e.g., Ford et al. 2000; Blaes et al. 2002; Katz et al. 2011; Naoz et al. 2011, 2013; Naoz 2016). This tends to widen the inclination window for large eccentricity excitation. However, the analytic expression for e_lim given by Equation (44) remains valid even for ${\varepsilon }_{\mathrm{oct}}\ne 0$ (Liu et al. 2015; Anderson & Lai 2017). In other words, because of the effect of short-range forces due to GR, the maximum eccentricity cannot exceed e_lim even when the octupole potential is significant. Higher eccentricity may be achieved when the double averaging approximation breaks down (see Section 3.2).

Figure 2 summarizes the parameter regimes of BH triples in terms of the mass (m₃) and semimajor axis (a_out) of the tertiary companion. For concreteness, we consider a fixed set of inner binary parameters (m₁ = 30 M_⊙, m₂ = 20 M_⊙, and a₀ = 100 au), with e_out = 0 (upper panel) and e_out = 0.6 (lower panel). The stability of the triple requires (e.g., Mardling & Aarseth 2001)

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{\mathrm{out}}}{a}\gt 2.8{\left(1+\displaystyle \frac{{m}_{3}}{{m}_{12}}\right)}^{2/5}\displaystyle \frac{{(1+{e}_{\mathrm{out}})}^{2/5}}{{(1-{e}_{\mathrm{out}})}^{6/5}}\left(1-\displaystyle \frac{0.3{I}_{0}}{180^\circ }\right).\end{eqnarray} \tag{ 47 }$$

In Figure 2, several regions have been identified (color-coded) and the boundaries are given by the various criteria (relations (20), (36), (46), and (47); see the dotted curves). We see that, in the rightmost region, the perturber is so far away that no LK oscillations occur (no e-excitation). In the "DA region," the dynamics of the system can be well described by the DA secular equations. The numbers shown on the dashed curves indicate the values of log₁₀(1 − e_lim), suggesting the extent of the excitation of eccentricity (Equation (44)). In the "SA region," the outer averaging fails, but the SA secular equations are valid.

Zoom In Zoom Out Reset image size

We first consider the cases when the octupole effect is negligible (ε_oct ≃ 0; see Equation (17)). These apply when the tertiary companion has a circular orbit (e_out = 0) or when the inner BHs have equal masses (m₁ = m₂). Figure 3 summarizes our results for a given set of binary and companion parameters (m₁ = m₃ = 30 M_⊙, m₂ = 20 M_⊙, a₀ = 100 au, and a_out = 4500 au) as a function of the initial mutual inclination angle I₀. All initial systems satisfy the criterion of double averaging for triples (relation (20)).

Zoom In Zoom Out Reset image size

The top panel of Figure 3 shows e_max for a grid of inclinations (uniformly distributed in cos I₀) in the absence of GW emission (this panel is similar to Figure 1, but is an enlarged version). We see that the eccentricity can be driven to be as large as ${e}_{\max }\simeq 1\mbox{--}{10}^{-6}$, at I₀ = I_0,lim ≃ 9216 (see Equation (43)). In the second panel, we plot ${I}_{{e}_{\max }}$, which is the instantaneous inclination at e = e_max, as a function of I₀. When e_max(I₀) achieves the maximum, ${I}_{{e}_{\max }}$ becomes very close to I_0,lim, implying that the range of oscillation in I is small (i.e., $\hat{{\boldsymbol{L}}}$ exhibits negligible nutation).

In the third panel of Figure 3, we turn on orbital decay due to gravitational radiation. The merger of the inner binary is achieved within the Hubble time (${T}_{{\rm{m}}}\lesssim {10}^{10}$ yr) for a range of inclinations around I_0,lim. The eccentricity excitation leads to a shorter binary merger time T_m compared to the "circular" merger time T_m,0 (see Equation (11)). In Liu & Lai (2017), we found that the merger timescale in LK-induced mergers can be described by the fitting formula ${T}_{{\rm{m}}}={T}_{{\rm{m}},0}{(1-{e}_{\max }^{2})}^{\alpha };$ the coefficient α depends on e_max (from Equation (42)), with α ≃ 1.5, 2, and 2.5 for e_max = (0, 0.6), (0.6, 0.8), and (0.8, 0.95), respectively. Here we consider the regime where e_max is much closer to unity, and we find that T_m can be best fitted by α = 3, i.e.,

$$\begin{eqnarray}&&{T}_{{\rm{m}}}\simeq \,{T}_{{\rm{m}},0}{(1-{e}_{\max }^{2})}^{3}.\end{eqnarray} \tag{ 48 }$$

This scaling can be understood as follows: the intrinsic GW-induced orbital decay rate $| \dot{a}/a{| }_{\mathrm{GW}}$ is proportional to (1 − e²)^−7/2 (Equation (10)). In an LK-induced merger, the orbital decay mainly occurs at e ≃ e_max. During the LK cycle, the binary only spreads a fraction ($\sim \sqrt{1-{e}_{\max }^{2}}$) of the time near e ≃ e_max. Thus, the LK-averaged orbital decay rate is of order ${T}_{{\rm{m}},0}^{-1}{(1-{e}_{\max }^{2})}^{-3}$, as indicated by Equation (48). Using Equation (48), we can define the "merger eccentricity" e_m via

$$\begin{eqnarray}&&{T}_{{\rm{m}},0}{(1-{e}_{{\rm{m}}}^{2})}^{3}={T}_{\mathrm{crit}}.\end{eqnarray} \tag{ 49 }$$

Thus, only systems with e_max ≳ e_m can have the merger time T_m less than T_crit—throughout this paper, our numerical results refer to T_crit = 10¹⁰ yr (see Section 3.3). For the systems shown in Figure 3, we find $1-{e}_{{\rm{m}}}\simeq {10}^{-4}$, and the merger window of initial inclinations is ${I}_{0,\mathrm{merger}}^{-}\leqslant {I}_{0}\leqslant {I}_{0,\mathrm{merger}}^{+}$, with ${I}_{0,\mathrm{merger}}^{-}=91\buildrel{\circ}\over{.} 56$ and ${I}_{0,\mathrm{merger}}^{+}=92\buildrel{\circ}\over{.} 76$ (Equation (42)), in agreement with the direct numerical results. As expected, the width of the merger window (${I}_{0,\mathrm{merger}}^{+}-{I}_{0,\mathrm{merger}}^{-}\simeq 1\buildrel{\circ}\over{.} 2$) is rather small. Also note that T_m shows a constant distribution around I₀ ∼ I_0,lim. This is the result of a "one-shot" merger, where the system only undergoes the first LK cycle, then "suddenly" merges during the high-e phase.

Figures 4–6 show a few examples of the orbital evolution for systems inside the merger window, for which the initial inclination is I₀ = 9252, 9233, and 9218, respectively. The evolution of BH spin is also shown, and this will be discussed in Section 4. In the three upper panels of Figure 4, we see that the inner binary undergoes cyclic excursions to the maximum eccentricity ${e}_{\max }$, with accompanying oscillations in the inclination I. As the binary decays, the range of eccentricity oscillations becomes smaller, and the eccentricity "freezes" at a large value. In the final phase, GW dissipation causes the orbit to circularize and its semimajor axis to shrink.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figures 5 and 6, I₀ is closer to I_0,lim, so that e_max achieved during LK oscillations is closer to e_lim. The GW-induced orbital decay is more efficient (Equations (8) and (9)), so the binary only experiences a few or even less than one LK cycles before merging. In Figure 5, the orbit undergoes the usual freezing of eccentricity oscillations as in Figure 4. In Figure 6, a decays abruptly, and the binary merges in the first high-eccentricity episode ("one-shot merger").

Figure 7 shows another example for a system with a more distant companion (a_out = 6700 au). Even though I₀ ≃ I_0,lim for this example, the inner BH binary does not attain sufficiently large e_max to enable a "one-shot" merger.

Zoom In Zoom Out Reset image size

For all the examples considered in Figures 3–7, we find that the merging BH binaries have a negligible eccentricity (e ≲ 0.01) when entering the aLIGO band.

The lower panel of Figure 8 shows the merger window (in terms of cos I₀) as a function of the effective semimajor axis of the tertiary companion. From Section 2.2 (see Figure 1), we have found that in the quadrupole approximation, the eccentricity excitation depends on m₃, a_out, and e_out through the ratio ${m}_{3}/{a}_{\mathrm{out},\mathrm{eff}}^{3}$ (where a_out,eff is given by Equation (19)). We therefore introduce the dimensionless scaled semimajor axis

$$\begin{eqnarray}{\bar{a}}_{\mathrm{out},\mathrm{eff}} & \equiv & \left(\displaystyle \frac{{a}_{\mathrm{out},\mathrm{eff}}}{1000\,\mathrm{au}}\right){\left(\displaystyle \frac{{m}_{3}}{30{M}_{\odot }}\right)}^{-1/3}\\ & = & \left(\displaystyle \frac{{a}_{\mathrm{out}}\sqrt{1-{e}_{\mathrm{out}}^{2}}}{1000\,\mathrm{au}}\right){\left(\displaystyle \frac{{m}_{3}}{30{M}_{\odot }}\right)}^{-1/3}\end{eqnarray} \tag{ 50 }$$

to characterize the "strength" of the outer perturber (note that Figure 8 neglects the octupole effect, which can complicate the single dependence of the merger window on ${\bar{a}}_{\mathrm{out},\mathrm{eff}};$ see Section 3.2). For a given BH binary (m₁ = 30 M_⊙, m₂ = 20 M_⊙, a₀ = 100 au), we fix e_out = 0 and m₃ = 30 M_⊙ or 10 M_⊙, but vary a_out. For each ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$, we consider 3000 values of I₀ spaced equally in cos I₀ ∈ (−1, 1), evolve the systems numerically, and record every merger event. The results obtained from the double- and single-averaged secular equations are marked by dark and light colors, respectively.

Zoom In Zoom Out Reset image size

The upper panel of Figure 8 shows the merger fraction from the mergers shown in the lower panel, which can also be characterized by the analyzed expression

$$\begin{eqnarray}&&{f}_{\mathrm{merger}}({a}_{0},{a}_{\mathrm{out}},{e}_{\mathrm{out}})=\displaystyle \frac{1}{2}| \cos {I}_{0,\mathrm{merger}}^{+}-\cos {I}_{0,\mathrm{merger}}^{-}| .\end{eqnarray} \tag{ 51 }$$

In the upper panel, the merger fraction is around ∼1%, and gradually decreases as ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ increases. The merger window is closed when ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ exceeds a certain critical value. Figure 8 also shows results for a different value of m₃ (10 M_⊙). This gives different values of η and I_0,lim (Equation (43)), but the merger window is qualitatively similar to the m₃ = 30 M_⊙ case, except shifted to relatively large values of I₀. Both merger windows (for the two values of m₃) are well described by Equations (42) and (49). The merger fractions, f_merger, are essentially identical (see the dashed curve in the upper panel). These indicate that the fitting formula (49), together with Equation (42), can be used to predict what types of systems will undergo a merger in less than 10¹⁰ yr, at least in the quadrupole order. For "pure" quadrupole systems, simple scaling relations for f_merger (as a function of a₀, m₁, m₂, and T_crit) can be obtained (see Section 3.3).

For ${m}_{1}\ne {m}_{2}$ and eccentric companions (${e}_{\mathrm{out}}\ne 0$), the octupole effect becomes important when ε_oct (Equation (17)) is appreciable, and some of the analytical expressions given in Section 2.2 break down. Previous works (Liu et al. 2015; Anderson et al. 2017) have shown that the main effect of the octupole potential is to broaden the range of the initial I₀ for extreme eccentricity excitations (e_max = e_lim), while the quadrupole expression for limiting eccentricity e_lim (Equation (44)) remains valid.

Figure 9 shows some examples of the merger windows for different values of ε_oct. To illustrate the effect of octupole perturbation, we consider four cases with the same ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ (Equation (50)), but different e_out (= 0, 0.3, 0.6, 0.9) and thus different a_out. All parameters in these examples satisfy the double averaging approximation (relation (20)). The initial longitude of the periapse ω_out is randomly chosen in the range (0, 2π).⁵ We find that, when the octupole effect gets stronger, the eccentricity excitation becomes increasingly erratic as a function of I₀, and more systems have the maximum eccentricity driven to be $1-{e}_{\max }\lesssim {10}^{-4}$. Consequently, more mergers over the Hubble timescale can be generated, and the merger window becomes noticeably broader. Because of the erratic variation of e_max, the merger events are not uniformly spaced in cos I₀. In this situation, the merger window cannot be described by the fitting formula (Equation (48); see the orange dashed curves in Figure 9).

Zoom In Zoom Out Reset image size

Figure 10 depicts an example of the time evolution of binary merger for ${\varepsilon }_{\mathrm{oct}}\ne 0$. Because of the octupole effect, the maximum eccentricities reached in successive (quadrupole) LK cycles increase. Eventually e_max becomes sufficiently large and the binary merges quickly.

Zoom In Zoom Out Reset image size

When ε_oct is sufficiently large, the orbital evolution of the inner binary becomes chaotic, and the evolution shows a strong dependence on the initial conditions (e.g., Lithwick & Naoz 2011; Li et al. 2014). Figure 11 illustrates this chaotic behavior. We see that the octupole-induced extreme eccentricity excitation occurs in an irregular manner, shortening or extending the time for mergers. As a result, T_m exhibits an irregular dependence on I₀, as seen in Figure 9.

Zoom In Zoom Out Reset image size

The merger windows shown in Figure 9 are based on the DA secular equations. For close and more eccentric companions, these DA equations break down, and we can use SA equations (see Section 2.1.2). Figure 12 shows sample numerical results for the merger windows computed using SA equations and DA equations. Here, ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ is chosen to be relatively small (≃5.6), where the system lies near the boundary of parameter regions between single and double averaging. The upper panel shows the eccentricity excitation for a grid of I₀ values when the system is non-dissipative (i.e., GW emission is turned off). We find that a portion of the systems computed from SA equations can reach higher e_max, even beyond e_lim (which is derived from the double averaging approximation). In particular, when gravitational radiation is included, a larger number of mergers with ${T}_{{\rm{m}}}\lesssim ({\rm{a}}\,{\rm{f}}{\rm{e}}{\rm{w}})\times {10}^{8}\,{\rm{y}}{\rm{r}}$ occur due to the extreme e_max, as depicted in the lower panel; such rapid mergers are relatively rare in the calculations based on the DA equations.

Zoom In Zoom Out Reset image size

Figure 13 shows the merger windows and merger fractions as a function of ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ (see Equation (50)) for different values of e_out. In our calculations, the orientation of the initial  e_out (for a given I₀) is random (i.e., ω_out,0 is uniformly distributed in 0–2π). We see that, for a given e_out, the merger window shows a general trend of widening as ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ decreases, resulting in an increase in f_merger. Moreover, for the same value of ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ (thus the same quadrupole effect), the merger window and merger fraction can be very different for different e_out. In general, the larger the eccentricity e_out, the stronger the octupole effect, and therefore the wider the window. Compared to the analytical expressions based on the quadrupole approximation (see Section 3.1), f_merger can be enhanced by a factor of a few. Note that for some values of ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$, an irregular distribution of merger events inside the merger window appears; this results from the chaotic behaviors of the octupole-order LK oscillations (see also the examples in Figure 9, particularly the e_out = 0.6 case).

Zoom In Zoom Out Reset image size

Although in this paper we focus on BH binaries, a similar analysis can be done for NS binary mergers induced by tertiary companions. A double NS merger event (GW170817) has recently been detected through GWs and electromagnetic radiation (e.g., Abbott et al. 2017d). We can expect more such detections in the future.

NS binaries differ from BH binaries in that the NS mass is much smaller than the BH mass, and thus for the same initial a_in = a₀ (≳1 au), a larger eccentricity excitation is required to induce NS binary merger. Moreover, since the masses of the two members of NS binaries are typically quite similar, the octupole LK effect is negligible (ε_oct ≃ 0). Therefore, the mergers of NS binaries in the presence of distant companions can be well described in the quadrupole approximation (see Section 3.1). Thus, the maximum eccentricity required for mergers (within time T_crit) can be obtained from Equation (49), and the required initially mutual inclination can be calculated using Equation (42) by replacing e_max with e_m. In other words, the merger window and merger fraction for NS binaries can be calculated analytically (Equation (51)), without the need for numerical integrations of the single- or double-averaged equations.

Figure 14 presents the results for the merger window and merger fraction for equal-mass NS binaries (m₁ = m₂ = 1.4 M_⊙). All systems shown satisfy the stability criterion. We choose three different initial semimajor axes (a₀ = 1, 10, and 100 au).⁶ For each NS binary, we consider a variety of tertiary bodies (different m₃ and e_out, as labeled). We find that, for a given a₀, different m₃ and e_out (with the same ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$) affect the position of the merger window (i.e., the range of cos I₀) but not the value of f_merger (see Figure 8). On the other hand, the merger windows and fractions have a strong dependence on the initial semimajor axis (e.g., f_merger for a₀ = 1 au is about 100 times larger than that for a₀ = 100 au). This is because, for the small a₀, the induced eccentricity in the LK oscillations does not have to be too large to produce mergers within 10¹⁰ yr (e.g., $1-{e}_{{\rm{m}}}\simeq {10}^{-3}$ for a₀ = 1 au, 10⁻⁴ for a₀ = 10 au, and 10⁻⁶ for a₀ = 100 au). In addition, the range of ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ producing a merger is different for different a₀.

Zoom In Zoom Out Reset image size

The result of Figure 14 (upper panel) for the merger fraction can be summarized by the fitting formula for ${f}_{\mathrm{merger}}^{\max }$, the maximum value of f_merger (for a given a₀), and ${\bar{a}}_{\mathrm{out},\mathrm{eff}}^{\max }$, the maximum value of ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ for a merger to be possible. Figure 15 shows that for the parameters of Figure 14 (with m₁ = m₂ = 1.4 M_⊙, ${T}_{\mathrm{crit}}={10}^{10}$ yr), we have

$$\begin{eqnarray}&&{f}_{\mathrm{merger}}^{\max }\,\simeq \,5.8 \% {\left(\displaystyle \frac{{a}_{0}}{\mathrm{au}}\right)}^{-0.67},\,{\bar{a}}_{\mathrm{out},\mathrm{eff}}^{\max }\simeq 0.23{\left(\displaystyle \frac{{a}_{0}}{\mathrm{au}}\right)}^{1.1}.\end{eqnarray} \tag{ 52 }$$

Zoom In Zoom Out Reset image size

Scaling relations for general quadrupole systems. Equation (52) can be generalized to other types of systems (with different m₁, m₂) and different value of merger time T_crit. From Equations (11) and (49), we see that the critical eccentricity e_m required for merger within time T_crit depends on ${\mu {T}_{{\rm{crit}}}({m}_{12}/{a}_{0}^{2})}^{2}$. From Equation (42) we see that for η ≪ 1 (a good approximation), the critical inclinations (${I}_{0,\mathrm{merger}}^{\pm }$; see Equation (51)) for a given e_max = e_m depend only on ε_GR, or the combination ${({m}_{12}/{a}_{0}^{2})}^{2}({a}_{\mathrm{out},\mathrm{eff}}^{3}/{m}_{3})$. Thus the merger fraction f_merger depends on m₁, m₂, a₀, and T_crit only through ${m}_{12}/{a}_{0}^{2}$ and μT_crit. We therefore expect from Equation (52) that ${f}_{\mathrm{merger}}^{\max }\,\propto \,{({a}_{0}/{m}_{12}^{0.5})}^{-0.67}{(\mu {T}_{\mathrm{crit}})}^{\alpha }$ and ${\bar{a}}_{\mathrm{out},\mathrm{eff}}^{\max }\,\propto \,{({a}_{0}/{m}_{12}^{0.5})}^{1.1}{(\mu {T}_{\mathrm{crit}})}^{\beta }$, where α, β are fitting parameters. Figure 16 shows the fitting. We find

$$\begin{eqnarray}{f}_{{\rm{m}}{\rm{e}}{\rm{r}}{\rm{g}}{\rm{e}}{\rm{r}}}^{{\max }} & \simeq & 5.8{\rm{ \% }}{\left[\left(\displaystyle \frac{{a}_{0}}{{\rm{a}}{\rm{u}}}\right){\left(\displaystyle \frac{{m}_{12}}{2.8{M}_{\odot }}\right)}^{-0.5}\right]}^{-0.67}\\ & & \times \,{\left(\displaystyle \frac{\mu }{0.7{M}_{\odot }}\displaystyle \frac{{T}_{{\rm{c}}{\rm{r}}{\rm{i}}{\rm{t}}}}{{10}^{10}{\rm{y}}{\rm{r}}}\right)}^{0.16},\end{eqnarray} \tag{ 53 }$$

and

$$\begin{eqnarray}{\bar{a}}_{{\rm{o}}{\rm{u}}{\rm{t}},{\rm{e}}{\rm{f}}{\rm{f}}}^{max} & \simeq & 0.23{\left[\left(\displaystyle \frac{{a}_{0}}{{\rm{a}}{\rm{u}}}\right){\left(\displaystyle \frac{{m}_{12}}{2.8{M}_{\odot }}\right)}^{-0.5}\right]}^{1.1}\\ & & \times \,{\left(\displaystyle \frac{\mu }{0.7{M}_{\odot }}\displaystyle \frac{{T}_{{\rm{c}}{\rm{r}}{\rm{i}}{\rm{t}}}}{{10}^{10}{\rm{y}}{\rm{r}}}\right)}^{0.06}.\end{eqnarray} \tag{ 54 }$$

These fitting formulae are valid for any type of LK-induced BH/NS merger in the quadrupole order.

Zoom In Zoom Out Reset image size

We now study how the BH spin evolves during LK-induced binary mergers. We present the evolution equation for ${{\boldsymbol{S}}}_{1}={S}_{1}{\hat{{\boldsymbol{S}}}}_{1}$ (where S₁ is the magnitude of the spin angular momentum of m₁ and ${\hat{{\boldsymbol{S}}}}_{1}$ is the unit vector). The de Sitter precession of ${\hat{{\boldsymbol{S}}}}_{1}$ around $\hat{{\boldsymbol{L}}}$ (1.5 PN effect) is governed by (e.g., Barker & O'Connell 1975)

$$\begin{eqnarray}&&\displaystyle \frac{d{\hat{{\boldsymbol{S}}}}_{1}}{{dt}}={{\rm{\Omega }}}_{\mathrm{SL}}\hat{{\boldsymbol{L}}}\times {\hat{{\boldsymbol{S}}}}_{1},\end{eqnarray} \tag{ 55 }$$

with the orbital-averaged spin precession rate

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{SL}}=\displaystyle \frac{3{Gn}({m}_{2}+\mu /3)}{2{c}^{2}a(1-{e}^{2})}.\end{eqnarray} \tag{ 56 }$$

A similar equation applies to the spinning body 2. Note that Ω_SL is of the same order as Ω_GR (Equation (7)) for m₁ ∼ m₂. There are also back-reaction torques from ${\hat{{\boldsymbol{S}}}}_{1}$ on $\hat{{\boldsymbol{L}}}$ and $\hat{{\boldsymbol{e}}}$:

$$\begin{eqnarray}{\left.\displaystyle \frac{d{\boldsymbol{L}}}{{dt}}\right|}_{\mathrm{LS}} & = & {{\rm{\Omega }}}_{\mathrm{LS}}{\hat{{\boldsymbol{S}}}}_{1}\times {\boldsymbol{L}},\\ {\left.\displaystyle \frac{d{\boldsymbol{e}}}{{dt}}\right|}_{\mathrm{LS}} & = & {{\rm{\Omega }}}_{\mathrm{LS}}{\hat{{\boldsymbol{S}}}}_{1}\times {\boldsymbol{e}}-3{{\rm{\Omega }}}_{\mathrm{LS}}(\hat{{\boldsymbol{L}}}\cdot {\hat{{\boldsymbol{S}}}}_{1})\hat{{\boldsymbol{L}}}\times {\boldsymbol{e}},\end{eqnarray} \tag{ 57 }$$

where

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{LS}}={{\rm{\Omega }}}_{\mathrm{SL}}\displaystyle \frac{{S}_{1}}{L}=\displaystyle \frac{{{GS}}_{1}(4+3{m}_{2}/{m}_{1})}{2{c}^{2}{a}^{3}{(1-{e}^{2})}^{3/2}}.\end{eqnarray} \tag{ 58 }$$

We include this effect in our calculations, although it (Equation (57) is usually negligible since S₁ ≪ L.⁷ The spin–spin coupling (2 PN correction) is always negligible until the final phase of the merger, and will be ignored in our calculations. In addition, the de Sitter precession of ${\hat{{\boldsymbol{S}}}}_{1}$ induced by the tertiary companion is neglected as well (since we consider m₁ ∼ m₂, and m₃ is not much larger, but a_out ≫ a). In all our calculations, we use χ₁ = χ₂ = 0.1 for concreteness. But note that the values of χ₁ and χ₂ do not affect the results of our paper (except Figure 20, which assumes χ₁ = χ₂). This is because (i) the de Sitter precession frequency Ω_SL (Equation (56)) is independent of χ₁, χ₂, and (ii) the inequality S₁, S₂ ≪ L is well satisfied (see footnote 7).

In terms of the inner BH binary axis $\hat{{\boldsymbol{L}}}$, the effect of the companion is to induce precession of $\hat{{\boldsymbol{L}}}$ around ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$ with nutation (when ${\text{}}e\ne 0$). In the quadrupole order, we have

$$\begin{eqnarray}{\left.\displaystyle \frac{d{\boldsymbol{L}}}{{dt}}\right|}_{\mathrm{LK},\mathrm{quad}} & = & -\displaystyle \frac{3L}{4{t}_{\mathrm{LK}}\sqrt{1-{e}^{2}}}\\ & & \times \ [({\boldsymbol{j}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}}){\hat{{\boldsymbol{L}}}}_{\mathrm{out}}\times {\boldsymbol{j}}+5({\boldsymbol{e}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}}){\boldsymbol{e}}\times {\hat{{\boldsymbol{L}}}}_{\mathrm{out}}].\end{eqnarray} \tag{ 59 }$$

An approximate expression for the rate of change of $\hat{{\boldsymbol{L}}}$ is given by (Anderson et al. 2016)

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{L}}}={\left|\displaystyle \frac{d\hat{{\boldsymbol{L}}}}{{dt}}\right|}_{\mathrm{LK},\mathrm{quad}}\simeq \displaystyle \frac{3(1+4{e}^{2})}{8{t}_{\mathrm{LK}}\sqrt{1-{e}^{2}}}| \sin 2I| .\end{eqnarray} \tag{ 60 }$$

The spin evolution is determined by two competing processes: ${\hat{{\boldsymbol{S}}}}_{1}$ precesses around $\hat{{\boldsymbol{L}}}$ at the rate Ω_SL, and $\hat{{\boldsymbol{L}}}$ varies at the rate Ω_L. There are three possible spin behaviors depending on the ratio Ω_SL/Ω_L:

• (i)  
  For Ω_L ≫ Ω_SL ("nonadiabatic" regime), the spin axis ${\hat{{\boldsymbol{S}}}}_{1}$ cannot "keep up" with the rapidly changing $\hat{{\boldsymbol{L}}}$, which precesses around a fixed ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$ (for L_out ≫ L). Thus ${\hat{{\boldsymbol{S}}}}_{1}$ effectively precesses around ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$, keeping the misalignment angle between ${\hat{{\boldsymbol{S}}}}_{1}$ and ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$, ${\theta }_{\mathrm{sb}}\equiv {\cos }^{-1}({\hat{{\boldsymbol{S}}}}_{1}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}})$, approximately constant.
• (ii)  
  For Ω_SL ≫ Ω_L ("adiabatic" regime), ${\hat{{\boldsymbol{S}}}}_{1}$ is strongly coupled to $\hat{{\boldsymbol{L}}}$. The spin axis ${\hat{{\boldsymbol{S}}}}_{1}$ closely "follows" $\hat{{\boldsymbol{L}}}$, maintaining an approximately constant spin–orbit misalignment angle ${\theta }_{\mathrm{sl}}\equiv {\cos }^{-1}({\hat{{\boldsymbol{S}}}}_{1}\cdot \hat{{\boldsymbol{L}}})$.
• (iii)  
  For Ω_SL ∼ Ω_L ("trans-adiabatic" regime), the spin evolution can be complex, potentially generating large spin–orbit misalignment θ_sl. Since both Ω_SL and Ω_L depend on e during the LK cycles, the precise transitions between these regimes can be fuzzy.

To help characterize the spin dynamics, we introduce an "adiabaticity parameter" as

$$\begin{eqnarray}&&{ \mathcal A }\equiv \left|\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{SL}}}{{{\rm{\Omega }}}_{{\rm{L}}}}\right|={{ \mathcal A }}_{0}\displaystyle \frac{1}{(1+4{e}^{2})\sqrt{1-{e}^{2}}| \sin 2I| },\end{eqnarray} \tag{ 61 }$$

where

$$\begin{eqnarray}{{ \mathcal A }}_{0} & \equiv & {\left|\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{SL}}}{{{\rm{\Omega }}}_{{\rm{L}}}}\sin 2I\right|}_{e=0}=\displaystyle \frac{4G({m}_{2}+\mu /3){m}_{12}{a}_{\mathrm{out},\mathrm{eff}}^{3}}{{c}^{2}{m}_{3}{a}^{4}}\\ & \simeq & \ 2.76\times {10}^{-5}\left[\displaystyle \frac{({m}_{2}+\mu /3)}{35\,{M}_{\odot }}\right]\left(\displaystyle \frac{{m}_{12}}{60\,{M}_{\odot }}\right)\\ & & \times \ {\left(\displaystyle \frac{{m}_{3}}{30{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{{a}_{\mathrm{out},\mathrm{eff}}}{{10}^{3}\mathrm{au}}\right)}^{3}{\left(\displaystyle \frac{a}{{10}^{2}\mathrm{au}}\right)}^{-4}.\end{eqnarray} \tag{ 62 }$$

Note that ${ \mathcal A }$ has a steep dependence on the eccentricity e and inclination I, and it is time-varying, while ${{ \mathcal A }}_{0}$ is an intrinsic indicator for identifying which system may undergo potentially complicated spin evolution. Since ${{ \mathcal A }}_{0}$ depends sensitively on a, during the orbital decay a system may transit from "nonadiabatic" at large a to "adiabatic" at small a, where the final spin–orbit misalignment angle ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ is "frozen." Note that ${{ \mathcal A }}_{0}$ is directly related to ε_GR (see Equation (41)) by

$$\begin{eqnarray}&&\displaystyle \frac{{{ \mathcal A }}_{0}}{{\varepsilon }_{\mathrm{GR}}}=\displaystyle \frac{4}{3}\displaystyle \frac{{m}_{2}+\mu /3}{{m}_{12}}.\end{eqnarray} \tag{ 63 }$$

Thus, when the initial value of ε_GR (at a = a₀) satisfies ε_GR,0 ≲ 9/4 (a necessary condition for LK eccentricity excitation; see relation (46)), we also have the initial ${{ \mathcal A }}_{0}\lesssim (3{m}_{2}+\mu )/{m}_{12}\sim 1$. This implies that any system that experiences enhanced orbital decay due to LK oscillations must go through the "trans-adiabatic" regime and therefore possibly complicated spin evolution (Liu & Lai 2017).

In our previous study (Liu & Lai 2017), we considered initially compact BH binaries (with a₀ ∼ 0.2 au), which can merge by themselves without the aid of a tertiary companion. We focused on systems with initial ${{ \mathcal A }}_{0}$ not much less than unity, and showed that such systems can experience complex/chaotic spin evolution during the LK-enhanced mergers. In this paper, we consider the inner BH binaries with large initial semimajor axis (a₀ = 20 or 100 au) and initial ${{ \mathcal A }}_{0}\ll 1$. As we shall see, such systems exhibit a variety of different spin evolutionary behaviors during the LK-induced mergers.

The bottom panel of Figure 4 shows a representative example of the spin evolution during the LK-induced orbital decay. The system begins with ${{ \mathcal A }}_{0}\sim {10}^{-3}\ll 1$ (Equation (62)). At the early stage of the evolution, Ω_SL ≪ Ω_L, leading ${\theta }_{\mathrm{sb}}$ to be nearly constant. Because of the large variation of $\hat{{\boldsymbol{L}}}$, the spin–orbit angle θ_sl oscillates with a large amplitude. As the orbit decays and circularizes, $\hat{{\boldsymbol{L}}}$ becomes frozen relative to ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$ (with final inclination I ≃ 125°), while $\hat{{\boldsymbol{S}}}$ precesses rapidly around $\hat{{\boldsymbol{L}}}$, with θ_sl settling down to the final value (≃90°). A non-zero final spin–orbit misalignment has been produced from the originally aligned configuration. This is only one example of the complex BH spin evolutionary paths during LK-induced mergers (Section 4.2).

The problem we study here is similar to the problem of the dynamics of stellar spin driven by a giant planet undergoing LK oscillations and migration (Storch et al. 2014, 2017; Storch & Lai 2015; Anderson et al. 2016). However, there is an important difference: the de Sitter precession of the BH spin is always prograde with respect to the orbit (the precession rate vector is ${{\rm{\Omega }}}_{\mathrm{SL}}\hat{{\boldsymbol{L}}}$), while the Newtonian precession of the stellar spin driven by the planet arises from the rotation-induced stellar oblateness and depends on $\cos {\theta }_{\mathrm{sl}}$ (the precession rate vector is along the direction of $-\cos {\theta }_{\mathrm{sl}}\hat{{\boldsymbol{L}}}$). This difference implies that the (Newtonian) stellar spin axis is prone to resonant (and potentially chaotic) excitation of spin–orbit misalignment, even for a circular orbit (e.g., Lai 2014; Lai et al. 2018), while the BH spin evolution is more regular: the nodal precession of the inner orbit driven by the external companion (i.e., the precession of $\hat{{\boldsymbol{L}}}$ and ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$) is retrograde (see Equation (59)), whereas the precession of $\hat{{\boldsymbol{S}}}$ around $\hat{{\boldsymbol{L}}}$ is prograde, so secular resonance does not usually happen when the orbital evolution is regular.

In the case of NS binaries, a Newtonian effect due to the oblateness of the NS (m₁) also contributes to the spin precession. Equation (55) is changed to

$$\begin{eqnarray}&&\displaystyle \frac{d{\hat{{\boldsymbol{S}}}}_{1}}{{dt}}={{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{dS})}\hat{{\boldsymbol{L}}}\times {\hat{{\boldsymbol{S}}}}_{1}+{{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{Newtonian})}\hat{{\boldsymbol{L}}}\times {\hat{{\boldsymbol{S}}}}_{1},\end{eqnarray} \tag{ 64 }$$

where

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{Newtonian})}=-\displaystyle \frac{3{{Gm}}_{2}({I}_{3}-{I}_{1})}{2{a}^{3}{(1-{e}^{2})}^{3/2}}\displaystyle \frac{\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}}{{S}_{1}}.\end{eqnarray} \tag{ 65 }$$

Here, I₃ and I₁ are principal moments of inertia of the NS. For ${I}_{3}-{I}_{1}\equiv {k}_{q\ast }{m}_{1}{R}_{1}^{2}{\hat{{\rm{\Omega }}}}_{1}^{2}$ and ${S}_{1}={I}_{3}{{\rm{\Omega }}}_{1}={k}_{* }{m}_{1}{R}_{1}^{2}{{\rm{\Omega }}}_{1}$, where R₁ is the NS radius and ${\hat{{\rm{\Omega }}}}_{1}$ is the rotation rate of the NS in units of ${({{Gm}}_{1}/{R}_{1}^{3})}^{1/2}$, we have

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{Newtonian})}=-\displaystyle \frac{3{k}_{q* }}{2{k}_{* }}\left(\displaystyle \frac{{m}_{2}}{{m}_{1}}\right){\left(\displaystyle \frac{{R}_{1}}{a}\right)}^{3}\displaystyle \frac{{{\rm{\Omega }}}_{1}}{{(1-{e}^{2})}^{3/2}}\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}.\end{eqnarray} \tag{ 66 }$$

Thus the ratio ${{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{Newtonian})}/{{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{dS})}$ is

$$\begin{eqnarray}\left|\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{Newtonian})}}{{{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{dS})}}\right| & = & \left(\displaystyle \frac{{k}_{q* }}{{k}_{* }}\right)\displaystyle \frac{{m}_{2}}{\sqrt{{m}_{1}{m}_{12}}}{\left[\displaystyle \frac{{R}_{1}}{a(1-{e}^{2})}\right]}^{1/2}\\ & & \times \ \displaystyle \frac{{R}_{1}{c}^{2}}{G({m}_{2}+\mu /3)}{\hat{{\rm{\Omega }}}}_{1}| \cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}| .\end{eqnarray} \tag{ 67 }$$

For a typical NS, m₁ ≃ 1.4 M_⊙, R₁ ≃ 10 km, k_q* ≃ 0.17, k_* ≃ 0.26,⁸ and ${\hat{{\rm{\Omega }}}}_{1}\simeq 0.023{({P}_{1}/20\mathrm{ms})}^{-1}$ (where P₁ is the rotation period of the NS). Since a(1 − e²) > R₁, we see that $| {{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{Newtonian})}/{{\rm{\Omega }}}_{\mathrm{SL}}^{(\mathrm{dS})}|$ is always ≪1.

We have seem from Section 3 that LK-induced BH binaries can have a variety of orbital evolution paths toward the final merger. Correspondingly, the evolution of BH spin in these binaries also exhibits a rich set of evolutionary behaviors. They can be roughly divided into four cases (see Figure 17).

Zoom In Zoom Out Reset image size

Case I (see Figures 4 and 5). This usually occurs when the initial inclination I₀ is sufficiently different from I_0,lim (see Equation (43)) that (1 − e_max) is much larger than (1 − e_lim). In this case, the inner binary experiences multiple LK oscillations; the amplitude of the eccentricity oscillations shrinks gradually as the orbit decays; eventually the binary circularizes and merges quickly. As shown in the lower panels of Figures 4 and 5, during the early stage, the angle ${\theta }_{\mathrm{sb}}$ is approximately constant (since ${{ \mathcal A }}_{0}\ll 1$), while θ_sl exhibits a larger variation due to the rapid precession of $\hat{{\boldsymbol{L}}}$ around ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}};$ at the later stage, as the orbit decays, $\hat{{\boldsymbol{L}}}$ becomes fixed relative to ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$, while $\hat{{\boldsymbol{S}}}$ precesses rapidly around $\hat{{\boldsymbol{L}}}$ with a fixed final θ_sl close to 90°.

Case II (Figures 10 and 11). This occurs when I₀ is not close to I_0,lim, but e_max is driven to a value close to e_lim due to the octupole effect. As seen from Figures 10 and 11, the inner binary experiences multiple LK cycles, each with increasing e_max driven by the octupole potential; eventually e_max becomes sufficiently large and the orbit decays rapidly. Unlike case I, the spin evolution transitions from the "nonadiabatic" regime to the "adiabatic" regime quickly. Because of the extremely rapid orbital decay, the oscillation of θ_sl continues to the end (by contrast, in Figures 4 and 5, the θ_sl oscillation freezes as the orbit decays), and the final θ_sl lies in the range ${\theta }_{\mathrm{sl}}^{{\rm{f}}}\in (0,\pi )$.

Case III (Figure 7). This occurs when I₀ is close to I_0,lim. Similarly to Case I, the orbit goes through eccentricity oscillations, suppression of the oscillations, and circularization. However, since I₀ ≈ I_0,lim, the orbital inclination oscillates with a small amplitude and passes through 90°. This implies that ${\hat{{\boldsymbol{S}}}}_{1}$ stays fairly close to $\hat{{\boldsymbol{L}}}$ at the early stage (see Figure 17) and θ_sl does not experience large-amplitude (0–π) oscillations. Eventually, θ_sl settles down to a value below 90°.

Case IV (see Figure 6). This also occurs when I₀ is close to I_0,lim. Similarly to Case II, the inner binary experiences extreme eccentricity excitation and merges within one LK cycle (one-shot merger). Because $\hat{{\boldsymbol{L}}}$ basically does not evolve in time (see Figure 17), a small (<90°) ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ is produced.

It is clear that the spin evolution is complicated and depends on various parameters and timescales. Our understanding of the spin behaviors is based largely on the numerical integrations. The four cases discussed above are representative, and do not capture the complete sets of spin evolutionary behaviors.

Figures 3 and 9 (bottom panels) show the final distribution of ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ as a function of cos I₀ in the merger window for several different systems. When e_out = 0, ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ has a regular distribution, and most of the values are found around ≲90°; the spin evolution follows the examples in Case I (${\theta }_{\mathrm{sl}}^{{\rm{f}}}\simeq 90^\circ$), Case III, and Case IV (one-shot merger) discussed above. When ${e}_{\mathrm{out}}\ne 0$, ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ shows a much wider range of values from 0° to 180°, due to the octupole effect (as in Case II discussed above).

Note that, for small e_out, the final spin–orbit misalignment angles ${\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and ${\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$ are strongly correlated; this correlation is particularly strong for the e_out = 0 case (see Figure 9). This arises because the orbital evolution is regular for small e_out. The spin vectors ${\hat{{\boldsymbol{S}}}}_{1}$ and ${\hat{{\boldsymbol{S}}}}_{2}$ evolve independently during the orbital decay (since spin–spin coupling is negligible). Although the de Sitter precession rates of ${\hat{{\boldsymbol{S}}}}_{1}$ and ${\hat{{\boldsymbol{S}}}}_{2}$ are different (since ${m}_{1}\ne {m}_{2}$), the spin evolution is regular as long as ${{ \mathcal A }}_{0}\ll 1$ (corresponding to Cases I, III, and IV discussed above in this section). In particular, the "90° attractor" is a generic feature independent of the precise value of Ω_SL (see Section 4.3). In contrast, for high e_out (see the e_out = 0.9 case in Figure 9), the octupole effect makes the orbital evolution chaotic, which also induces chaotic evolution of the spin–orbit misalignment (see Figure 11). Therefore, ${\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and ${\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$ become largely uncorrelated.

We see from the previous subsections that when the octupole effect is negligible (ε_oct ≪ 1), the spin–orbit misalignment angle (starting from initial ${\theta }_{\mathrm{sl}}^{0}=0^\circ$) often evolves toward ${\theta }_{\mathrm{sl}}^{{\rm{f}}}\simeq 90^\circ$ as the binary orbit decays. What is the origin of this 90° "attractor"?

In Liu & Lai (2017), we used the principle of adiabatic invariance to derive an analytical expression for ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ for the case where the inner BH binary remains circular in the presence of an inclined tertiary companion (i.e., the inner binary merges by itself without eccentricity excitation, although the binary orbital angular momentum axis $\hat{{\boldsymbol{L}}}$ does vary and precess around ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$). We can use a similar idea to understand qualitatively the origin of the 90° attractor in quadrupole LK-induced mergers. In the quadrupole order, the angular momentum axis $\hat{{\boldsymbol{L}}}$ varies at the rate given by Equation (60). This variation involves precession around ${\hat{{\boldsymbol{L}}}}_{\mathrm{out}}$ and nutation (change in I). If we neglect nutation, we have⁹

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d\hat{{\boldsymbol{L}}}}{{dt}}\right|}_{\mathrm{LK},\mathrm{quad}}\simeq -{{\rm{\Omega }}}_{\mathrm{pl}}{\hat{{\boldsymbol{L}}}}_{\mathrm{out}}\times \hat{{\boldsymbol{L}}}=-{{\rm{\Omega }}}_{\mathrm{pl}}\displaystyle \frac{{{\boldsymbol{L}}}_{\mathrm{tot}}}{{L}_{\mathrm{out}}}\times \hat{{\boldsymbol{L}}},\end{eqnarray} \tag{ 68 }$$

where  L_tot =  L +  L_out, and

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{pl}}=\displaystyle \frac{3\hat{{\boldsymbol{L}}}\cdot {\hat{{\boldsymbol{L}}}}_{\mathrm{out}}}{4{t}_{\mathrm{LK}}\sqrt{1-{e}^{2}}}(1+4{e}^{2}).\end{eqnarray} \tag{ 69 }$$

Equation (68) shows that $\hat{{\boldsymbol{L}}}$ rotates around the  L_tot axis. In this rotating frame, the spin evolution Equation (55) transforms to

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d{\hat{{\boldsymbol{S}}}}_{1}}{{dt}}\right)}_{\mathrm{rot}}={{\boldsymbol{\Omega }}}_{\mathrm{eff}}\times {\hat{{\boldsymbol{S}}}}_{1},\end{eqnarray} \tag{ 70 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{\Omega }}}_{\mathrm{eff}}={{\rm{\Omega }}}_{\mathrm{SL}}\hat{{\boldsymbol{L}}}+{{\rm{\Omega }}}_{\mathrm{pl}}\displaystyle \frac{{{\boldsymbol{L}}}_{\mathrm{tot}}}{{L}_{\mathrm{out}}}.\end{eqnarray} \tag{ 71 }$$

Note that the ratio between Ω_SL and Ω_pl is

$$\begin{eqnarray}&&\left|\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{SL}}}{{{\rm{\Omega }}}_{\mathrm{pl}}}\right|={ \mathcal A }\sin I=\displaystyle \frac{{{ \mathcal A }}_{0}}{2(1+4{e}^{2})\sqrt{1-{e}^{2}}| \cos I| },\end{eqnarray} \tag{ 72 }$$

where ${ \mathcal A }$, ${{ \mathcal A }}_{0}$ are given by Equations (61) and (62).

If we assume that ${\hat{{\boldsymbol{\Omega }}}}_{\mathrm{eff}}\equiv {{\boldsymbol{\Omega }}}_{\mathrm{eff}}/| {{\boldsymbol{\Omega }}}_{\mathrm{eff}}|$ varies slowly (much more slowly than $| {{\boldsymbol{\Omega }}}_{\mathrm{eff}}| ;$ see below), then θ_eff,S1, the angle between ${\hat{{\boldsymbol{S}}}}_{1}$ and ${{\boldsymbol{\Omega }}}_{\mathrm{eff}}$, is an adiabatic invariant. Suppose ${\hat{{\boldsymbol{S}}}}_{1}$ and $\hat{{\boldsymbol{L}}}$ are aligned initially (${\theta }_{\mathrm{sl}}^{0}=0^\circ$), then the initial ${\theta }_{\mathrm{eff},{S}_{1}}^{0}$ equals the initial ${\theta }_{\mathrm{eff},{\rm{L}}}^{0}$ (the angle between ${{\boldsymbol{\Omega }}}_{\mathrm{eff}}$ and $\hat{{\boldsymbol{L}}}$), which is given by

$$\begin{eqnarray}&&\tan {\theta }_{\mathrm{eff},{\rm{L}}}^{0}=\displaystyle \frac{\sin {I}_{0}}{({{ \mathcal A }}_{0}/2\cos {I}_{0})+{\eta }_{0}+\cos {I}_{0}},\end{eqnarray} \tag{ 73 }$$

where η₀ is the initial value of η = L/L_out (see Equation (38)).¹⁰ On the other hand, after the binary has decayed, $\eta \to 0$, $| {{\rm{\Omega }}}_{\mathrm{SL}}| \gg | {{\rm{\Omega }}}_{\mathrm{pl}}|$, and thus ${{\boldsymbol{\Omega }}}_{\mathrm{eff}}={{\rm{\Omega }}}_{\mathrm{SL}}\hat{{\boldsymbol{L}}}$, which implies ${\theta }_{\mathrm{sl}}^{{\rm{f}}}\simeq {\theta }_{\mathrm{eff},{S}_{1}}^{{\rm{f}}}$. Therefore, under adiabatic evolution, we have

$$\begin{eqnarray}&&{\theta }_{\mathrm{sl}}^{{\rm{f}}}\simeq {\theta }_{\mathrm{eff},{S}_{1}}\simeq {\theta }_{\mathrm{eff},{\rm{L}}}^{0}.\end{eqnarray} \tag{ 74 }$$

For systems with η₀ ≪ 1, $| \cos {I}_{0}| \ll 1$, and ${{ \mathcal A }}_{0}/| \cos {I}_{0}| \ll 1$, Equation (73) gives ${\theta }_{\mathrm{eff},{\rm{L}}}^{0}\approx 90^\circ$, and thus adiabatic evolution predicts ${\theta }_{\mathrm{sl}}^{{\rm{f}}}\approx 90^\circ$.

Figure 18 shows the evolution of θ_eff,S1 for the four cases discussed in Section 4.2. We see that for Case I (Figures 4 and 5), θ_eff,S1 is approximately constant throughout the evolution of the inner binary, and the adiabatic evolution correctly predicts the 90° attractor in the spin–orbit misalignment. For the other cases (Cases II–IV), θ_eff,S1 undergoes significant change during the inner binary's evolution, especially near the final orbital decay phase; in these cases, Equation (73) does not predict the correct ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$.

Zoom In Zoom Out Reset image size

The validity of adiabatic evolution requires that the rate of change of ${\hat{{\boldsymbol{\Omega }}}}_{\mathrm{eff}}={{\boldsymbol{\Omega }}}_{\mathrm{eff}}/| {{\boldsymbol{\Omega }}}_{\mathrm{eff}}|$ be much slower than $| {{\boldsymbol{\Omega }}}_{\mathrm{eff}}|$, i.e., $| d{\hat{{\boldsymbol{\Omega }}}}_{\mathrm{eff}}/{dt}| \ll | {{\boldsymbol{\Omega }}}_{\mathrm{eff}}|$. In order of magnitude, we have $| d{\hat{{\boldsymbol{\Omega }}}}_{\mathrm{eff}}/{dt}| \sim \,{T}_{\mathrm{GW}}^{-1}$ (see Equation (10)).¹¹ In Case I, e_max induced by the tertiary companion is not too extreme. So the orbital decay is "gentle" and the adiabatic condition is satisfied. In Cases II, III, and IV, the rapid orbital decay at high eccentricity implies ${T}_{\mathrm{GW}}^{-1}\gtrsim | {{\boldsymbol{\Omega }}}_{\mathrm{eff}}|$, so the adiabatic evolution breaks down.

We reiterate that the above analysis cannot be considered rigorous, since the precession rate in Equation (68) is approximate and the nutation of $\hat{{\boldsymbol{L}}}$ has been neglected. Nevertheless, this analysis (especially Figure 18) provides a qualitative understanding as to why θ_sl evolves toward 90° under some conditions.

Having studied the various spin evolutionary paths in the previous subsections, we now calculate the distribution of final spin–orbit misalignment angle for the merging systems studied in Figure 13. We consider the spins of both BHs, and assume that both ${{\boldsymbol{S}}}_{1}$ and ${{\boldsymbol{S}}}_{2}$ are initially aligned with respect to the binary orbital angular momentum axis.

Figure 19 summarizes our results for a₀ = 100 au (the results for a₀ = 20 au are very similar). The range of ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ values considered in this figure all lie in the regime where the double averaging approximation is valid (see Section 2.3). As in Figure 13, four different values of e_out are considered. When e_out = 0 (the top left panels), the spin evolution is regular, following the examples of Case I, Case III, and Case IV (see Section 4.2). For ${\bar{a}}_{\mathrm{out},\mathrm{eff}}\in (4.2,6.5)$, the tertiary companion is relatively close, and many BH binaries inside the merger window pass through successive stages of LK oscillations, LK suppression, and orbital circularization (Case I), producing a large number of systems with ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ around 90°. For ${\bar{a}}_{\mathrm{out},\mathrm{eff}}\in (6.6,8.8)$, the tertiary companion is relatively distant; the eccentricity cannot grow to be as in the case of small ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$, even when I₀ ≃ I_0,lim. Thus, the spin mainly evolves as described in Case III, and ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ lies in the range 0°–90°.

Zoom In Zoom Out Reset image size

When the companions are eccentric (${e}_{\mathrm{out}}\ne 0$), the octupole effect comes into play and the spin may follow the dynamics of Case II. For a given e_out, when the companion is relatively close (${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ is small), the octupole effect becomes more prominent. In the case of e_out = 0.9, the orbital evolution is dominated by the octupole effect, and the distribution of ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ is close to being isotropic (i.e., uniform distribution in $\cos {\theta }_{\mathrm{sl}}^{{\rm{f}}}$), as shown in the bottom-right panel of Figure 19.

Since the two components of the BH binary have comparable masses, the de Sitter precession rates are similar. Thus it is not surprising that the distributions of ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ for both spins are similar. Note that ${\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and ${\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$ are strongly correlated for e_out = 0, and this correlation becomes much weaker as the octupole effect becomes stronger (see Figure 9).

Having obtained the distributions of $\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and $\cos {\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$ for a range of systems with different parameters, we can compute the distribution of the effective spin parameter for the merging binaries (see Equation (1)):

$$\begin{eqnarray}&&{\chi }_{\mathrm{eff}}=\displaystyle \frac{{m}_{1}{\chi }_{1}\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}+{m}_{2}{\chi }_{2}\cos {\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}}{{m}_{12}},\end{eqnarray} \tag{ 75 }$$

where χ_1,2 are the dimensionless BH spins (we set χ₁ = χ₂ = 0.1 in our calculations, although our results for ${\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and ${\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$ are not affected by this choice since S₁, S₂ ≪ L for all the systems considered in this paper). Figure 20 shows two examples (for a₀ = 100 au and 20 au; see Figures 13 and 19), assuming χ₁ = χ₂.¹² To obtain the χ_eff distribution, we consider systems with ${\bar{a}}_{\mathrm{out},\mathrm{eff}}\in (5.6,8.8)$ for the a₀ = 100 au case and ${\bar{a}}_{\mathrm{out},\mathrm{eff}}\in (0.8,1.4)$ for the a₀ = 20 au case, and assume that the eccentricity of the tertiary companion has a uniform distribution in e_out (i.e., e_out = 0, 0.3, 0.6, 0.9 are equally probable), and that the initial mutual inclination is randomly distributed (uniform in cos I₀). We see that although the systems with the most eccentric companion (e_out = 0.9) contribute a substantial fraction of mergers, the overall distribution of χ_eff has a peak around 0. Importantly, our result indicates that LK-induced BH binary mergers can easily have χ_eff < 0 (see also Liu & Lai 2017). This is quite different from the standard isolated binary evolution channel, where we typically expect spin–orbit alignment and χ_eff > 0.

Zoom In Zoom Out Reset image size

If the distributions of $\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and $\cos {\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$ are uncorrelated, as we may expect to be the case for ${m}_{1}\ne {m}_{2}$, when the octupole effect is significant (see Figure 9 and the discussion in the last paragraph of Section 4.2), the distribution of χ_eff can be derived directly from ${P}_{1}(\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}})$ and ${P}_{2}(\cos {\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}})$, the distribution functions of $\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and $\cos {\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$. Define

$$\begin{eqnarray}&&{\bar{\chi }}_{1}\equiv \displaystyle \frac{{m}_{1}{\chi }_{1}}{{m}_{1}{\chi }_{1}+{m}_{2}{\chi }_{2}},\end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray}&&{\bar{\chi }}_{2}\equiv \displaystyle \frac{{m}_{2}{\chi }_{2}}{{m}_{1}{\chi }_{1}+{m}_{2}{\chi }_{2}},\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&{\bar{\chi }}_{\mathrm{eff}}\equiv \displaystyle \frac{{\chi }_{\mathrm{eff}}}{{\chi }_{\mathrm{eff}}^{\max }}\equiv \displaystyle \frac{{m}_{12}{\chi }_{\mathrm{eff}}}{{m}_{1}{\chi }_{1}+{m}_{2}{\chi }_{2}},\end{eqnarray} \tag{ 78 }$$

where ${\chi }_{\mathrm{eff}}^{\max }=({m}_{1}{\chi }_{1}+{m}_{2}{\chi }_{2})/{m}_{12}$ is the maximum possible value of χ_eff for given m₁χ₁ and m₂χ₂ (this maximum is achieved at $\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}=\cos {\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}=1$). Then Equation (75) becomes

$$\begin{eqnarray}&&{\bar{\chi }}_{\mathrm{eff}}={\bar{\chi }}_{1}{\mu }_{1}+{\bar{\chi }}_{2}{\mu }_{2},\end{eqnarray} \tag{ 79 }$$

where ${\mu }_{1}\equiv \cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and ${\mu }_{2}\equiv \cos {\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$. Note that ${\bar{\chi }}_{1}+{\bar{\chi }}_{2}=1$ and ${\bar{\chi }}_{\mathrm{eff}}\in [-1,1]$. Given P₁(μ₁) and P₂(μ₂), the distribution function of ${\bar{\chi }}_{\mathrm{eff}}$ is

$$\begin{eqnarray}P({\bar{\chi }}_{\mathrm{eff}}) & = & {\displaystyle \int }_{-1}^{1}d{\mu }_{1}{P}_{1}({\mu }_{1}){\displaystyle \int }_{-1}^{1}d{\mu }_{2}{P}_{2}({\mu }_{2})\\ & & \times \ \delta ({\bar{\chi }}_{\mathrm{eff}}-{\bar{\chi }}_{1}{\mu }_{1}-{\bar{\chi }}_{2}{\mu }_{2}).\end{eqnarray} \tag{ 80 }$$

In the special case when μ₁ and μ₂ are uniformly distributed, we have P₁ = P₂ = 1/2, and Equation (80) gives

$$\begin{eqnarray}P({\bar{\chi }}_{\mathrm{eff}})=\left\{\begin{array}{rcl}(1-{\bar{\chi }}_{\mathrm{eff}})/(4{\bar{\chi }}_{1}{\bar{\chi }}_{2}), & & {\bar{\chi }}_{\mathrm{eff}}\geqslant {\bar{\chi }}_{1}-{\bar{\chi }}_{2}\\ 1/(2{\bar{\chi }}_{1}), & {\bar{\chi }}_{2}-{\bar{\chi }}_{1} & \leqslant {\bar{\chi }}_{\mathrm{eff}}\leqslant {\bar{\chi }}_{1}-{\bar{\chi }}_{2}\\ (1+{\bar{\chi }}_{\mathrm{eff}})/(4{\bar{\chi }}_{1}{\bar{\chi }}_{2}), & & {\bar{\chi }}_{\mathrm{eff}}\leqslant {\bar{\chi }}_{2}-{\bar{\chi }}_{1},\end{array}\right.\end{eqnarray} \tag{ 81 }$$

where we have assumed ${\bar{\chi }}_{1}\geqslant {\bar{\chi }}_{2}$ without loss of generality. Thus, even for uniform distributions of $\cos {\theta }_{{{\rm{s}}}_{1}{\rm{l}}}^{{\rm{f}}}$ and $\cos {\theta }_{{{\rm{s}}}_{2}{\rm{l}}}^{{\rm{f}}}$ (see the case of e_out = 0.9 in Figure 19), the effective spin parameter χ_eff is preferentially distributed around χ_eff = 0 (see Figure 20).

Note that the spin–orbit misalignment distribution and χ_eff distribution obtained above refer to relatively wide BH binary systems (a₀ ≳ 10–100 au) that experience a merger due to large LK eccentricity excitation. Such systems necessarily have ${{ \mathcal A }}_{0}\ll 1$. For BH binaries with smaller separations (a₀ ≲ 1 au) and ${{ \mathcal A }}_{0}$ not much less than unity, the spin–orbit misalignment distribution can be quite different (see Liu & Lai 2017).

Antonini et al. (2018; see also Rodriguez & Antonini 2018) have carried population studies of BH mergers in triple systems (based on DA secular equations) and have found a similar peak around χ_eff = 0 in the χ_eff distribution. Rodriguez & Antonini (2018) also showed an example of the final spin–orbit misalignment distribution with a peak around 90°, in qualitative agreement with our result. They did not distinguish the difference in the spin–orbit misalignment distributions between largely quadrupole systems (small ε_oct) and strong octupole systems (large ε_oct). We do not agree with the reason(s) they gave for the peak in the χ_eff distribution. In particular, it is important to recognize that, during the orbital decay, the adiabaticity parameter ${ \mathcal A }$ (Equation (61)) transitions from ≪1 to ≫1, and this transition determines ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$.

• 1.  
  For BH binaries with a given initial separation (a₀ ≳ 10 au), the merger window (i.e., the range of initial inclination angles I₀ between the inner binary and the outer companion that induces binary merger within ∼10¹⁰ yr) and merger fraction depend on the effective semimajor axis ${\bar{a}}_{\mathrm{out},\mathrm{eff}}\propto {a}_{\mathrm{out}}\sqrt{1-{e}_{\mathrm{out}}^{2}}/{m}_{3}^{1/3}$ (Equation (50)) and eccentricity e_out of the companion. The results are summarized in Figure 13. Assuming that the inclination of the companion is randomly distributed, we find that the merger fraction (for typical BH masses m₁ = 30 M_⊙, m₂ = 20 M_⊙) increases rapidly with increasing e_out, from ∼1% at e_out = 0 to 10%–20% at e_out = 0.9. This is because, as the octupole potential ($\propto {\varepsilon }_{\mathrm{oct}}\propto {e}_{\mathrm{out}};$ see Equation (17)) of the tertiary companion increases, extreme eccentricity excitation of the inner binaries becomes possible for a wide range of I₀ (see Figure 9). Regardless of the importance of the octupole effect, the maximum ${\bar{a}}_{\mathrm{out},\mathrm{eff}}$ value for which the inner binary has a chance to merge within 10¹⁰ yr (or any other values) can be determined analytically (using Equations (44) and (49), setting e_m to e_lim; see also Equation (54)).
• 2.  
  For systems where the octupole effect is negligible (such as those with m₁ = m₂ or e_out = 0), the merger window and merger fraction can be determined analytically (see Figure 8). In particular, these analytical results can be applied to NS–NS binaries with external companions (see Section 3.3, Figure 14). We have also obtained fitting formulae relevant to the merger fractions of various systems (Equations (53) and (54)).
• 3.  
  On the technical side, we have developed new dynamical equations for the evolution of triples (Section 2.1.2) in the single averaging approximation (i.e., the equations of motion are only averaged over the inner orbit). These SA equations have a wider regime of validity in the parameter space than the usual DA secular equations (see Section 2.3 and Figure 2). For systems where the octupole effect is negligible, we find that the DA equations accurately predict the merger window and merger fractions even in the regime where the equations formally break down (see Figure 8). However, when the octupole effect is strong (large ε_oct), using the SA equations leads to wider merger windows and larger merger fractions (see Figures 12 and 13).
• 4.  
  During the tertiary-induced binary decay, the spin axes of the BHs exhibit a variety of evolutionary behaviors due to the combined effects of spin–orbit coupling (de Sitter precession), LK orbital precession/nutation, and GW emission. These spin behaviors are correlated with the orbital evolution of the BH binary (Section 4.2). Starting from aligned spin axes (relative to the orbital angular momentum axis), a wide range of spin–orbit misalignments can be generated when the binary enters the LIGO/VIRGO band:

  • a.  
    For systems where the octupole effect is negligible (such as those with m₁ ≃ m₂ or e_out ∼ 0), the BH spin axis evolves regularly, with the final spin–orbit misalignment angle ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ depending on the initial companion inclination angle I₀ in a well-defined manner (see Figure 3 and the top left panels of Figure 9).¹³ We find that when I₀ is not too close to I_lim (the initial inclination angle for maximum/limiting eccentricity excitation; see Equations (43) and (44)), the spin–orbit misalignment evolves into a 90° "attractor" (Figures 4 and 5), a feature that can be qualitatively understood using adiabatic invariance (see Section 4.3). When I₀ is close to I_lim, a qualitatively different spin evolution leads to smaller ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ (Figures 6 and 7).
  • b.  
    For systems with a stronger octupole effect (larger ε_oct), the BH spin evolution becomes increasingly chaotic, with the final spin–orbit misalignment angle depending sensitively on the initial conditions (see Figures 10 and 11). As a result, a wide range of ${\theta }_{\mathrm{sl}}^{{\rm{f}}}$ values are produced, including retrograde configurations (see Figure 9). The final spin–orbit misalignment distribution typically peaks around 90°, but becomes isotropic (uniform in $\cos {\theta }_{\mathrm{sl}}^{{\rm{f}}}$) for systems with sufficiently large ε_oct (Figure 19).
• 5.  
  We have computed the distribution of the mass-weighted spin parameter χ_eff (Equation (75)) of merging BH binaries in triples (Figure 20). While details of this distribution depend on various parameters (e.g., distribution of the companion eccentricities), it has a characteristic shape with peak around χ_eff ≃ 0, extending to the maximum possible positive and negative values (see Equation (78)).

The merger fraction f_merger computed in this paper (and particularly the dependence of f_merger on various parameters) can be used to obtain an estimate of the rate of LK-induced BH binary mergers in the galactic field, provided that one makes certain assumptions about the BH populations in triples and their properties. We do not present such an estimate here since such a calculation necessarily contains large uncertainties (see Antonini et al. 2017; Silsbee & Tremaine 2017), like all other scenarios of producing merging BH binaries. Suffice it to say that with our computed f_merger of a few to 10 per cent, it is possible to produce (with large error bars) the observed BH binary merger rate (10–200 Gpc⁻³ yr⁻¹).

As noted in Section 1, the mass-weighted spin parameter χ_eff may serve as a useful indicator of the binary BH formation mechanism. The five discovered BH binaries all have low values of χ_eff, which could be the result of either slowly spinning BHs (e.g., Zaldarriaga et al. 2017) or large spin–orbit misalignments. The event GW170104 has ${\chi }_{\mathrm{eff}}=-{0.12}_{-0.3}^{+0.21}$, which may require retrograde spinning BHs, especially if low individual spins (χ_1,2 ≲ 0.2) can be ruled out. Such a retrograde spin–orbit misalignment would challenge the isolated binary BH formation channel, and point to the importance of some flavors of dynamical formation mechanisms. We note that the LK-induced BH mergers lead to a unique shape of χ_eff distribution (Figure 20) that may be used to distinguish it from other types of dynamical interactions. For example, a completely random distribution of spin–orbit misalignments, as expected from the mechanisms involving multiple close encounters and exchange interactions in dense clusters (e.g., Rodriguez et al. 2015; Chatterjee et al. 2017), would lead to a specific distribution given by Equation (81). As the number of detected BH merger events increases in the coming years, the distribution of χ_eff will be measured experimentally, therefore providing valuable constraints on the binary BH formation mechanisms.

Although we have focused on isolated BH triples in this paper, many aspects of our results (with proper rescalings) can be applied to triples that dynamically form in globular clusters or BH binaries moving around a supermassive BH (e.g., Miller & Hamilton 2002; Wen 2003; Thompson 2011; Antonini & Perets 2012; Antonini et al. 2014; Petrovich & Antonini 2017; Hoang et al. 2018). In a dense cluster, the orbits of a triple system can be perturbed or even disrupted by close fly-bys of other objects. Therefore the survival timescale of the triple may not be as long as 10¹⁰ yr, depending on the mean density of the surroundings. In this case, the merger window and merger fraction may be reduced (see Equations (53) and (54)), and the remaining systems can lead to extremely large eccentricities and shorter merger times. Also, our conclusion on the distribution of near-merger spin–orbit misalignments depends on the initial BH spin orientations. We have assumed initial spin–orbit alignment throughout this paper, but this may not be valid for dynamically formed binaries and triples in dense clusters. We plan to address some of these issues in a future paper.