Spin-Orbit Misalignment of Merging Black-Hole Binaries with Tertiary Companions

TL;DR

This paper investigates how tertiary companions influence the spin-orbit misalignment in merging black-hole binaries, revealing chaotic spin evolution and potential observational signatures in gravitational wave signals.

Contribution

It provides a detailed analysis of the effects of external companions on black-hole binary evolution, including an analytic expression for final misalignment in non-eccentric cases.

Findings

Close inclined companions induce Lidov-Kozai oscillations, shortening merger times.

Chaotic spin evolution leads to wide misalignment angles from initially aligned systems.

Analytic formula for final misalignment in systems without eccentricity excitation.

Abstract

We study the effect of external companion on the orbital and spin evolution of merging black-hole (BH) binaries. An sufficiently close by and inclined companion can excite Lidov-Kozai (LK) eccentricity oscillations in the binary, thereby shortening its merger time. During such LK-enhanced orbital decay, the spin axis of the BH generally exhibits chaotic evolution, leading to a wide range ($0^\circ$-$180^\circ$) of final spin-orbit misalignment angle from an initially aligned configuration. For systems that do not experience eccentricity excitation, only modest ($\lesssim 20^\circ$) spin-orbit misalignment can be produced, and we derive an analytic expression for the final misalignment using the principle of adiabatic invariance. The spin-orbit misalignment directly impacts the gravitational waveform, and can be used to constrain the formation scenarios of BH binaries and dynamical influences of external companions.

Equations

$${\chi }_{\mathrm{eff}}=\frac{m_1{\mathbf{a}}_1+m_2{\mathbf{a}}_2}{m_1+m_2}\cdot \hat{\text{L}},$$

$${\text{L}}_{\mathrm{in}}={\mathrm{L}}_{\mathrm{in}}{\hat{\text{L}}}_{\mathrm{in}}={\mu }_{\mathrm{in}}\sqrt{G{m}_{12}{a}_{\mathrm{in}}(1-e_{\mathrm{in}}^2)}{\hat{\text{L}}}_{\mathrm{in}},$$

$${\text{L}}_{\mathrm{out}}={\mathrm{L}}_{\mathrm{out}}{\hat{\text{L}}}_{\mathrm{out}}={\mu }_{\mathrm{out}}\sqrt{G{m}_{123}{a}_{\mathrm{out}}(1-e_{\mathrm{out}}^2)}{\hat{\text{L}}}_{\mathrm{out}},$$

$$\begin{split}T_{\mathrm{m,0}}&=\frac{5c^{5}a_{0}^{4}}{256G^{3}m_{12}^{2}\mu}\\ &\simeq 10^{10}\bigg{(}\frac{60M_{\odot}}{m_{12}}\bigg{)}^{2}\bigg{(}\frac{15M_{\odot}}{\mu}\bigg{)}\bigg{(}\frac{a_{0}}{0.202\mathrm{AU}}\bigg{)}^{4}{\rm yrs}.\end{split}$$

$$\displaystyle\frac{d\textbf{L}}{dt}=\frac{d\textbf{L}}{dt}\bigg{|}_{\mathrm{LK}}+\frac{d\textbf{L}}{dt}\bigg{|}_{\mathrm{GW}}~{},$$

$$\displaystyle\frac{d\textbf{e}}{dt}=\frac{d\textbf{e}}{dt}\bigg{|}_{\mathrm{LK}}+\frac{d\textbf{e}}{dt}\bigg{|}_{\mathrm{GR}}+\frac{d\textbf{e}}{dt}\bigg{|}_{\mathrm{GW}}~{}.$$

$$t_{\mathrm{LK}}=\frac{1}{n}\frac{m_{12}}{m_{3}}\bigg{(}\frac{a_{\mathrm{out},\mathrm{eff}}}{a}\bigg{)}^{3},$$

$$\frac{d\textbf{e}}{dt}\bigg{|}_{\mathrm{GR}}=\Omega_{\mathrm{GR}}\hat{\textbf{L}}\times\textbf{e},~{}~{}~{}\Omega_{\mathrm{GR}}=\frac{3Gnm_{12}}{c^{2}a(1-e^{2})}.$$

$$\displaystyle\!\!\!\frac{3}{8}\frac{j^{2}_{\mathrm{min}}-1}{j^{2}_{\mathrm{min}}}\bigg{[}5\left(\cos I_{0}+\frac{\eta}{2}\right)^{2}-\Bigl{(}3+4\eta\cos I_{0}+\frac{9}{4}\eta^{2}\Bigr{)}j^{2}_{\mathrm{min}}$$

$$\displaystyle\quad+\eta^{2}j^{4}_{\mathrm{min}}\bigg{]}+\varepsilon_{\mathrm{GR}}\left(1-j_{\mathrm{min}}^{-1}\right)=0,$$

$$\displaystyle\varepsilon_{\mathrm{GR}}=t_{\rm LK}\Omega_{\rm GR}\Bigr{|}_{e=0}=\frac{3Gm_{12}^{2}a_{\mathrm{out},\mathrm{eff}}^{3}}{c^{2}a^{4}m_{3}}$$

$$\displaystyle\simeq 0.96\bigg{(}\!\frac{m_{12}}{60M_{\odot}}\!\bigg{)}^{\!\!2}\!\bigg{(}\!\frac{m_{3}}{30M_{\odot}}\!\bigg{)}^{\!\!-1}\!\bigg{(}\!\frac{a_{\mathrm{out},\mathrm{eff}}}{3\mathrm{AU}}\!\bigg{)}^{\!\!3}\!\bigg{(}\!\frac{a}{0.1\mathrm{AU}}\!\bigg{)}^{\!\!-4}.$$

$$\frac{3}{8}(j_{\lim }^2-1)\left[-3+\frac{{\eta }^2}{4}\left(\frac{4}{5}{j}_{\lim }^2-1\right)\right]+{\epsilon }_{\mathrm{GR}}\left(1-j_{\lim }^{-1}\right)=0,$$

$$(\cos I_{0})_{\pm}=\frac{1}{10}\Big{(}-\eta\pm\sqrt{\eta^{2}+60-\frac{80}{3}\varepsilon_{\mathrm{GR}}}\Big{)}.$$

$${\epsilon }_{\mathrm{GR}}\ge \frac{9}{4}+\frac{3}{80}{\eta }^2(\mathrm{no}\mathrm{eccentricity}\mathrm{excitation}).$$

$${\epsilon }_{\mathrm{oct}}\equiv \frac{m_1-m_2}{m_{12}}\left(\frac{a}{a_{\mathrm{out}}}\right)\frac{e_{\mathrm{out}}}{1-e_{\mathrm{out}}^2}$$

$$\frac{d{\hat{\text{S}}}_1}{dt}={\Omega }_{\mathrm{dS}}\hat{\text{L}}\times {\hat{\text{S}}}_1,{\Omega }_{\mathrm{dS}}=\frac{3Gn(m_2+\mu /3)}{2c^{2}a(1-e^2)}.$$

$$\frac{d\hat{\bf L}}{dt}\bigg{|}_{\mathrm{LK},e=0}=-\Omega_{\mathrm{L}}\hat{\bf L}_{\mathrm{out}}\times\hat{\bf L}=-\Omega_{\mathrm{L}}^{\prime}\hat{\bf J}\times\hat{\bf L},$$

$$\Omega_{\mathrm{L}}=\frac{3}{4t_{\mathrm{LK}}}\big{(}\hat{\textbf{L}}\cdot\hat{\textbf{L}}_{\mathrm{out}}\big{)},\quad\Omega_{\mathrm{L}}^{\prime}=\Omega_{\mathrm{L}}{J\over L_{\mathrm{out}}}.$$

$$\begin{split}\mathcal{A}&\equiv\bigg{(}\frac{\Omega_{\mathrm{dS}}}{\Omega_{\mathrm{L}}}\bigg{)}_{e,I=0}\simeq 0.37\bigg{[}\frac{(m_{2}+\mu/3)}{35M_{\odot}}\bigg{]}\bigg{(}\frac{m_{12}}{60M_{\odot}}\bigg{)}\\ &\times\bigg{(}\frac{m_{3}}{30M_{\odot}}\bigg{)}^{-1}\bigg{(}\frac{a_{\mathrm{out},\mathrm{eff}}}{3\mathrm{AU}}\bigg{)}^{3}\bigg{(}\frac{a}{0.1\mathrm{AU}}\bigg{)}^{-4}.\end{split}$$

$$\frac{\mathcal{A}}{{\epsilon }_{\mathrm{GR}}}=\frac{2}{3}\frac{m_2+\mu /3}{m_{12}}.$$

$$\bigg{(}\frac{d\hat{\textbf{S}}_{1}}{dt}\bigg{)}_{\mathrm{rot}}=\bm{\Omega}_{\mathrm{eff}}\times\hat{\textbf{S}}_{1},~{}~{}~{}\bm{\Omega}_{\mathrm{eff}}\equiv\Omega_{\mathrm{dS}}\hat{\textbf{L}}+\Omega^{\prime}_{\mathrm{L}}\hat{\textbf{J}}.$$

$$\tan\theta_{\mathrm{eff},\mathrm{L}}=\frac{\Omega_{\mathrm{L}}\sin I}{\Omega_{\mathrm{dS}}+\big{(}\eta+\cos I\big{)}\Omega_{\mathrm{L}}}.$$

$${\theta }_{\mathrm{eff},{\mathrm{S}}_1}\simeq \mathrm{constant}(\mathrm{adiabatic}\mathrm{invariant}).$$

$${\theta }_{\mathrm{SL}}^{\mathrm{f}}\simeq {\theta }_{\mathrm{eff},\mathrm{L}}^0.$$

$$\tan {\theta }_{\mathrm{eff},\mathrm{L}}^0=\frac{\sin I_0}{({\mathcal{A}}_0/\cos I_0)+{\eta }_0+\cos I_0}.$$

Full text

1. Introduction

The recent breakthrough in the detection of gravitational waves (GWs) from merging black hole (BH) binaries by advanced LIGO (Abbott et al., 2016a, b, 2017) has generated renewed interest in understanding the formation mechanisms of compact BH binaries, from the evolution of massive stellar binaries (Lipunov et al., 1997; Belczynski et al., 2010, 2016; Mandel & de Mink, 2016; Lipunov et al., 2017) and triples (Silsbee & Tremaine, 2017; Antonini et al., 2017) in the galactic fields, to dynamical interactions in galactic nuclei (Antonini & Perets, 2012; Petrovich & Antonini, 2017) and in the dense core of globular clusters (Miller & Hamilton, 2002; Rodriguez et al., 2015; Chatterjee et al., 2017).

Because of the uncertainties associated with various formation channels (e.g., common envelope evolution in the standard binary channel), it is difficult to distinguish the different formation mechanisms based on BH mass measurement alone. The detection of eccentric systems would obviously indicates some dynamical processes at work (e.g., Antonini & Perets (2012); Silsbee & Tremaine (2017)). However, because of the efficient eccentricity damping by gravitational radiation, the vast majority of compact binaries will likely be circular when entering the LIGO sensitivity band regardless of the formation channels. It has been suggested that the BH spin and the spin-orbit misalignment angle may be an important discriminant. In particular, the spin-orbit misalignment directly impacts the projected spin parameter of the merging binaries,

[TABLE]

(where $m_{1,2}$ are the BH masses, ${\bf a}_{1,2}=c{\bf S}_{1,2}/(Gm_{1,2}^{2})$ are the dimensionless BH spins, and $\hat{\textbf{L}}$ is the unit orbital angular momentum vector), which can be measured from the phase evolution of GWs (Abbott et al., 2016b, 2017).

In this paper, we study the merger and spin-orbit misalignment of BH binaries in the presence of tertiary companion. Such triple BH systems could be a direct product of massive triple stars in the galactic field (Silsbee & Tremaine, 2017; Antonini et al., 2017), or could be produced dynamically in a dense cluster (Miller & Hamilton, 2002; Rodriguez et al., 2015; Antonini & Rasio, 2016). For binaries formed near the center of a galaxy, the third body could be a supermassive BH (Antonini & Perets, 2012; Petrovich & Antonini, 2017).

It is well known that a tertiary body on an inclined orbit can accelerate the decay of an inner binary by inducing Lidov-Kozai (LK) eccentrcity/inclination oscillations (Lidov, 1962; Kozai, 1962). This has been studied before in the contexts of supermassive BH binary merger (Blaes et al., 2002) and stellar mass BH binaries (e.g., Miller & Hamilton (2002); Thompson (2011); Antonini et al. (2014); Silsbee & Tremaine (2017)). We focus on the latter in this paper. We show that as the BH binary undergoes LK-enhanced decay from a wide orbit and eventually enters the LIGO band, the spin axis of individual BHs can experience chaotic evolution, so that a significant spin-orbit misalignment can be produced prior to merger even for binaries formed with zero initial misalignment. We derive relevant analytic relations and quantify how the final spin-orbit misalignment angle depends on various parameters of the system (binary and external companion).

Note that in this paper we focus on triple systems with relatively small binary separations ($\lesssim 0.2$ AU for the inner binaries), so that the inner binary can merger within $10^{10}$ years either by itself or through modest ($e\lesssim 0.99$) LK eccentricity excitation. Such compact triple systems likely have gone through a complex (and highly uncertain) sequence of common envelop or mass transfer evolution – we do not study such evolution in this paper and therefore do not address issues related to the occurrence rate of compact triples. Our goal is to use such triple systems to illustrate the complex spin dynamics of the individual BHs. We expect that similar spin dynamics may take place in other types of triple systems, e.g., those with much larger initial separations, but experience extreme eccentricity excitation due to non-secular forcing from tertiary companions (Silsbee & Tremaine, 2017; Antonini et al., 2017).

2. Lidov-Kozai cycles in BH triples with gravitational radiation

2.1. Setup and Orbital Evolution

Consider a hierarchical triple system, consisting of an inner BH binary with masses $m_{1}$, $m_{2}$ and a relatively distant companion of mass $m_{3}$. The reduced mass for the inner binary is $\mu_{\mathrm{in}}\equiv m_{1}m_{2}/m_{12}$, with $m_{12}\equiv m_{1}+m_{2}$. Similarly, the outer binary has $\mu_{\mathrm{out}}\equiv(m_{12}m_{3})/m_{123}$ with $m_{123}\equiv m_{12}+m_{3}$. The orbital semimajor axes and eccentricities are denoted by $a_{\mathrm{in},\mathrm{out}}$ and $e_{\mathrm{in},\mathrm{out}}$, respectively. The orbital angular momenta of the inner and outer binaries are

[TABLE]

where $\hat{\textbf{L}}_{\mathrm{in},\mathrm{out}}$ are unit vectors. The relative inclination between $\hat{\bf L}_{\mathrm{in}}$ and $\hat{\bf L}_{\mathrm{out}}$ is denoted by $I$. For convenience, we will frequently omit the subscript “$\mathrm{in}$”.

The merger time due to GW radiation of an isolated binary with initial $a_{0}$ and $e_{0}=0$ is given by

[TABLE]

A sufficiently inclined external companion can raise the binary eccentricity through Lidov-Kozai oscillations, thereby reducing the merger time or making an otherwise non-merging binary merge within $10^{10}$ years. To study the evolution of merging BH binary under the influence of a companion, we use the secular equations to the octupole level in terms of the angular momentum L and eccentricity e vectors:

[TABLE]

Here the “Lidov-Kozai” (LK) terms are given explicitly in (Liu et al., 2015a) (we also evolve $\textbf{L}_{\mathrm{out}}$ and $\textbf{e}_{\mathrm{out}}$), and the associated timescale of LK oscillation is

[TABLE]

where $n=\sqrt{Gm_{12}/a^{3}}$ is the mean motion of the inner binary and $a_{\mathrm{out},\mathrm{eff}}\equiv a_{\mathrm{out}}\sqrt{1-e^{2}_{\mathrm{out}}}$. General Relativity (1-PN correction) induces pericenter precession

[TABLE]

We include GW emission (2.5-PN effect) that causes orbital decay and circularization (Shapiro & Teukolsky, 1983), but not the extreme eccentricity excitation due to non-secular effects (Antonini et al., 2014; Silsbee & Tremaine, 2017).

The top three panels of Figure 1 show an example of the orbital evolution of a BH binary with an inclined companion (initial $I_{0}=80^{\circ}$). We see that the inner binary undergoes cyclic excursions to maximum eccentricity $e_{\mathrm{max}}$, with accompanying oscillations in the inclination $I$. As the binary decays, the range of eccentricity oscillations shrinks. Eventually the oscillations freeze and the binary experiences “pure” orbital decay/circularization governed by GW dissipation.

2.2. Eccentricity Excitation and Merger Time

In the quadrupole approximation, the maximum eccentricity $e_{\mathrm{max}}$ attained in the LK oscillations (starting from an initial $I_{0}$ and $e_{0}\simeq 0$) can be calculated analytically using energy and angular momentum conservation, according to the equation (Anderson et al., 2017a)

[TABLE]

where $j_{\mathrm{min}}\equiv\sqrt{1-e_{\mathrm{max}}^{2}}$, $\eta\equiv(\mathrm{L}/\mathrm{L}_{\mathrm{out}})_{e=0}$ and $\varepsilon_{\mathrm{GR}}$ is given by

[TABLE]

Note that in the limit of $\eta\rightarrow 0$ and $\varepsilon_{\mathrm{GR}}\rightarrow 0$, Equation (9) yields the well-known relation $e_{\mathrm{max}}=\sqrt{1-(5/3)\cos^{2}I_{0}}$. The maximum possible $e_{\mathrm{max}}$ for all values of $I_{0}$, called $e_{\mathrm{lim}}$, is given by

[TABLE]

and is reached at $\cos I_{0}=(\eta/10)(4j_{\mathrm{lim}}^{2}-5)$. Eccentricity excitation ($e_{\mathrm{max}}\geq 0$) occurs within a window of inclinations $(\cos I_{0})_{-}\leqslant\cos I_{0}\leqslant(\cos I_{0})_{+}$, where (Anderson et al., 2017a)

[TABLE]

This window vanishes when

[TABLE]

Figure 2 shows some examples of the $e_{\mathrm{max}}(I_{0})$ curves. For $\eta\lesssim 1$, these curves depend mainly on $m_{3}/a_{\mathrm{out},\mathrm{eff}}^{3}$ (for given inner binary parameters). As $\varepsilon_{\mathrm{GR}}$ increases (with decreasing $m_{3}/a_{\mathrm{out},\mathrm{eff}}^{3}$), the LK window shrinks and $e_{\mathrm{lim}}$ decreases. Eccentricity excitation is suppressed (for all $I_{0}$’s) when Equation (13) is satisfied.

For systems with $m_{1}\neq m_{2}$ and $e_{\mathrm{out}}\neq 0$, so that

[TABLE]

is non-negligible, the octupole effect may become important (e.g. Ford et al. (2000); Naoz (2016)). This tends to widen the inclination window for large eccentricity excitation. However, the analytic expression for $e_{\mathrm{lim}}$ given by Equation (11) remains valid even for $\varepsilon_{\rm oct}\neq 0$ (Liu et al., 2015a; Anderson & Lai, 2017b). Therefore Equation (13) still provides a good criterion for “negligible eccentricity excitation” (note that when $\varepsilon_{\rm oct}\neq 0$, the inner binary cannot be exactly circular; e.g., Anderson & Lai (2017b); Liu et al. (2015b)).

Eccentricity excitation leads to a shorter binary merger time $T_{\mathrm{m}}$ compared to the circular merger time $T_{\mathrm{m},0}$ (Figure 3). We compute $T_{\mathrm{m}}$ by integrating the secular evolution equations of the BH triples with a range of $I_{0}$. Each $I_{0}$ run has a corresponding $e_{\mathrm{max}}$, which can be calculated using Equation (9). We see that in general $T_{\mathrm{m}}/T_{\mathrm{m},0}$ can be approximated by $(1-e_{\mathrm{max}}^{2})^{\alpha}$, with $\alpha\simeq 1.5,~{}2$ and 2.5 for $e_{\mathrm{max}}=(0,0.6),~{}(0.6,0.8)$ and $(0.8,0.95)$, respectively.

3. Spin-Orbit dynamics in merging BH binary with an external companion

3.1. Spin-Orbit Coupling

We now study how the BH spin evolves during the binary merger, considering only $\textbf{S}_{1}=\mathrm{S_{1}}\hat{\textbf{S}}_{1}$ (where $\hat{\textbf{S}}_{1}$ is the unit vector). The de Sitter precession of $\hat{\textbf{S}}_{1}$ around $\hat{\textbf{L}}$ is (Barker & O’Connell, 1975)

[TABLE]

Note that there is a back-reaction torque from ${\bf S}_{1}$ on ${\bf L}$; this can be safely neglected since $\mathrm{L}\gg\mathrm{S}_{1}$.

Before presenting our numerical results, it is useful to note the different regimes for the evolution of the spin-orbit misalignment angle $\theta_{\mathrm{SL}}$ (the angle between ${\bf S}_{1}$ and ${\bf L}$). In general, the inner binary axis $\hat{\bf L}$ precesses (and nutates when $e\neq 0$) around the total angular momentum ${\bf J}={\bf L}+{\bf L}_{\mathrm{out}}$ (recall that $\mathrm{S}_{1},\mathrm{S}_{2}\ll\mathrm{L}$, and ${\bf J}$ is constant in the absence of GW dissipation). The related precession rate $\Omega_{\mathrm{L}}$ is of order $t_{\rm LK}^{-1}$ at $e\sim 0$ (Equation 17), but increases with $e$. Depending on the ratio of $\Omega_{\rm dS}$ and $\Omega_{\mathrm{L}}$, we expect three possible spin behaviors: (i) For $\Omega_{\mathrm{L}}\gg\Omega_{\rm dS}$ (“nonadiabatic”), the spin axis $\hat{\bf S}_{1}$ cannot “keep up” with the rapidly changing $\hat{\textbf{L}}$, and thus effectively precesses around $\hat{\textbf{J}}$, keeping $\theta_{\mathrm{SJ}}\equiv\cos^{-1}(\hat{\textbf{S}}_{1}\cdot\hat{\textbf{J}})\simeq$ constant. (ii) For $\Omega_{\mathrm{dS}}\gg\Omega_{\mathrm{L}}$ (“adiabatic”), $\hat{\textbf{S}}_{1}$ closely “follows” $\hat{\textbf{L}}$, maintaining an approximately constant $\theta_{\mathrm{SL}}$. (iii) For $\Omega_{\mathrm{dS}}\sim\Omega_{\mathrm{L}}$ (“trans-adiabatic”), the spin evolution can be chaotic due to overlapping resonances. Since both $\Omega_{\mathrm{dS}}$ and $\Omega_{\mathrm{L}}$ depend on $e$ during the LK cycles, the precise transitions between these regimes can be fuzzy (Storch et al., 2014; Storch & Lai, 2015; Anderson et al., 2016; Storch et al., 2017).

For circular orbits ($e=0$), the precession of $\hat{\bf L}$ is governed by the equation

[TABLE]

where $\hat{\bf J}$ is the unit vector along ${\bf J}={\bf L}+{\bf L}_{\mathrm{out}}$, and

[TABLE]

We can define an “adiabaticity parameter”

[TABLE]

As the binary orbit decays, the system may transition from “non-adiabatic” (${\cal A}\ll 1$) at large $a$’s to “adiabatic” (${\cal A}\gg 1$) at small $a$’s, where the final spin-orbit misalignment angle $\theta_{\rm SL}^{\rm f}$ is “frozen”. Note that ${\cal A}$ is directly related to $\varepsilon_{\mathrm{GR}}$ by

[TABLE]

Thus, when the initial value of $\varepsilon_{\mathrm{GR}}$ (at $a=a_{0}$) satisfies $\varepsilon_{\mathrm{GR},0}\lesssim 9/4$ (a necessary condition for LK eccentricity excitation; see Equation 13), we also have ${\cal A}_{0}\lesssim(3m_{2}+\mu)/(2m_{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 chaotic spin evolution.

The bottom panel of Figure 1 shows a representative example of the evolution of the misalignment angle as the BH binary undergoes LK-enhanced orbital decay. We see that the BH spin axis can exhibit complex evolution even though the orbital evolution is “regular”. In particular, $\theta_{\mathrm{SL}}$ evolves in a chaotic way, with the final value $\theta_{\mathrm{SL}}^{\rm f}$ depending sensitively on the initial conditions (the precise initial $\theta_{\mathrm{SL}}$ and $I_{0}$). Also note that retrograde spin ($\theta_{\mathrm{SL}}^{\rm f}>90^{\circ}$) can be produced even though the binary always remains prograde with respect to the outer companion ($I<90^{\circ}$). These behaviors are qualitatively similar to the chaotic evolution of stellar spin driven by Newtonian spin-orbit coupling with a giant planet undergoing high-eccentricity migration (Storch et al., 2014; Storch & Lai, 2015; Anderson et al., 2016; Storch et al., 2017).

We carry out a series of numerical integrations, evolving the orbit of the merging BH binary with a tertiary companion, along with spin-orbit coupling, to determine $\theta_{\mathrm{SL}}^{\mathrm{f}}$ for various triple parameters. In our “population synthesis” study, we consider initial conditions such that $\hat{\textbf{S}}_{1}$ is parallel to $\hat{\textbf{L}}$, the binary inclinations are isotropically distributed (uniform distribution in $\cos I_{0}$), and the orientations of e and $\textbf{e}_{\mathrm{out}}$ (for systems with $e_{\mathrm{out}}\neq 0$) are random. All initial systems satisfy the criterion of dynamical stability for triples (Mardling & Aarseth, 2001).

Figure 4 shows our results for systems with equal masses, $e_{\mathrm{out}}=0$ (so that the octupole effect vanishes), and several values of $a_{\mathrm{out}}$. We see that when the eccentricity of the inner binary is excited ($I_{0}$ lies in the LK window), a wide range of $\theta^{\mathrm{f}}_{\mathrm{SL}}$ is generated, including appreciable fraction of retrograde ($\theta_{\mathrm{SL}}^{\mathrm{f}}>90^{\circ}$) systems (see the $a_{\mathrm{out}}=3$ AU case, for which $e_{\mathrm{lim}}=0.84$). The “memory” of chaotic spin evolution is evident, as slightly different initial inclinations lead to vastly different $\theta_{\mathrm{SL}}^{\mathrm{f}}$. The regular behavior of $\theta_{\mathrm{SL}}^{\mathrm{f}}$ around $I_{0}=90^{\circ}$ (again for the $a_{\mathrm{out}}=3$ AU case) is intriguing, but may be understood using the theory developed in Ref. (Storch et al., 2017). For systems with no eccentricity excitation, $\theta^{\mathrm{f}}_{\mathrm{SL}}$ varies regularly as a function of $I_{0}$ – this can be calculated analytically (see below).

Figure 5 shows our results for systems with $m_{1}\neq m_{2}$ and $e_{\mathrm{out}}\neq 0$, for which octupole terms may affect the orbital evolution. We see that eccentricity excitation and the corresponding reduction in $T_{\mathrm{m}}$ occur outside the analytic (quadupole) LK window (see the $e_{\mathrm{out}}=0.8$ case, for which $e_{\mathrm{lim}}=0.66$). As in the equal-mass case (Figure 4), a wide range of $\theta_{\mathrm{SL}}^{\mathrm{f}}$ values are produced whenever eccentricity excitation occurs. A larger fraction ($23\%$) of systems attain retrograde spin ($\theta_{\mathrm{SL}}^{\mathrm{f}}>90^{\circ}$). Again, for systems with negligible eccentricity excitation, $\theta_{\mathrm{SL}}^{\mathrm{f}}$ behaves regularly as a function of $I_{0}$ and agrees with the analytic result (the “fuzziness” of the numerical result in this regime is likely due to the very small eccentricity of the inner binary; see Liu et al. (2015b)).

3.2. Analytical Calculation of $\theta^{\mathrm{f}}_{\mathrm{SL}}$ for Circular Binaries

If the inner binary experiences no eccentricity excitation and remains circular throughout the orbital decay, the final spin-orbit misalignment can be calculated analytically using the principle of adiabatic invariance.

Equation (16) shows that $\hat{\bf L}$ rotates around the $\hat{\bf J}$ axis at the rate $(-\Omega_{\mathrm{L}}^{\prime})$. In this rotating frame, the spin evolution equation (15) transforms to

[TABLE]

Note that in the absence of GW dissipation, $\hat{\textbf{L}}$ and $\hat{\textbf{L}}_{\mathrm{out}}$ are constants (in the rotating frame), and thus $\hat{\bf S}_{1}$ precesses with a constant $\bm{\Omega}_{\mathrm{eff}}$. The relative inclination between $\bm{\Omega}_{\mathrm{eff}}$ and $\hat{\bf L}$ is given by

[TABLE]

Now if we include GW dissipation, $\hat{\textbf{L}}\cdot\hat{\textbf{L}}_{\mathrm{out}}=\cos I$ is exactly conserved, and $\bm{\Omega}_{\mathrm{eff}}$ becomes a slowly changing vector. When the rate of change of $\bm{\Omega}_{\mathrm{eff}}$ is much smaller than $|\bm{\Omega}_{\mathrm{eff}}|$, the angle between $\bm{\Omega}_{\mathrm{eff}}$ and $\hat{\textbf{S}}_{1}$ is adiabatic invariant, i.e.

[TABLE]

This adiabatic invariance requires $|d\bm{\Omega}_{\mathrm{eff}}/dt|/|\bm{\Omega}_{\mathrm{eff}}|\sim T_{m,0}^{-1}\ll|\bm{\Omega}_{\mathrm{eff}}|$, or $|\bm{\Omega}_{\mathrm{eff}}|T_{\mathrm{m,0}}\gg 1$, which is easily satisfied.

Suppose $\hat{\textbf{S}}_{1}$ and $\hat{\textbf{L}}$ are aligned initially, we have $\theta_{\mathrm{eff},\mathrm{S}_{1}}^{0}=\theta_{\mathrm{eff},\mathrm{L}}^{0}$ (the superscript [math] denotes initial value). Equation (22) then implies $\theta_{\mathrm{eff},\mathrm{S}_{1}}\simeq\theta_{\mathrm{eff},\mathrm{L}}^{0}$ at all times. After the binary has decayed, $\eta\rightarrow 0$, $|\Omega_{\mathrm{dS}}|\gg|\Omega_{\mathrm{L}}|$, and therefore $\bm{\Omega}_{\mathrm{eff}}\simeq\Omega_{\mathrm{dS}}\hat{\textbf{L}}$, which implies $\theta_{\mathrm{SL}}^{\mathrm{f}}\simeq\theta_{\mathrm{eff},\mathrm{S}_{1}}$. Thus we find

[TABLE]

That is, the final spin-orbit misalignment angle is equal to the initial inclination angle between $\hat{\bf L}$ and $\bm{\Omega}_{\mathrm{eff}}$, obtained by evaluating Equation (21) at $a=a_{0}$:

[TABLE]

This analytic expression agrees with the numerical results shown in Figures 4-5 in the appropriate regime. Note that for systems that experience no eccentricity excitation for all $I_{0}$’s, ${\cal A}_{0}\gtrsim(3m_{2}+\mu)/(2m_{12})\sim 1$ (see Equations 13 and 19), and thus $\theta_{\mathrm{SL}}^{\mathrm{f}}\simeq\theta_{\mathrm{eff},\mathrm{L}}^{0}\lesssim 20^{\circ}$, i.e., only modest spin-orbit misalignment can be generated. For systems with ${\cal A}_{0}\gg 1$ (e.g., very distant companion), we have $\theta_{\mathrm{SL}}^{\mathrm{f}}\ll 1$.

The above analytical result can be easily generalized to the situation of non-zero initial spin-orbit misalignment. It shows that $\theta_{\mathrm{SL}}^{\mathrm{f}}\simeq\theta_{\mathrm{SL}}^{0}$ for ${\cal A}_{0}\gg 1$.

4. Summary and Discussion

We have studied the effect of external companion on the orbital and spin evolution of merging BH binaries due to gravitational radiation. A sufficiently close by and inclined companion can excite Lidov-Kozai eccentricity oscillation in the binary, shortening its merger time compared to circular orbits [see Fig.3]. We find that during the LK-enhanced orbital decay, the spin axis of the BH generally experiences complex, chaotic evolution, with the final spin-orbit misalignment angle $\theta_{\mathrm{SL}}^{\mathrm{f}}$ depending sensitively on the initial conditions. A wide range of $\theta_{\mathrm{SL}}^{\mathrm{f}}$ (including $\theta_{\mathrm{SL}}^{\mathrm{f}}>90^{\circ}$) can be produced from an initially aligned ($\theta_{\mathrm{SL}}^{0}=0$) configuration (see Figs.4-5). For systems that do not experience eccentricity excitation (because of relatively low orbital inclinations of the companion or/and suppression by GR-induced precession), modest ($\lesssim 20^{\circ}$) spin-orbit misalignment can be produced – we have derived an analytic expression for $\theta_{\mathrm{SL}}^{\mathrm{f}}$ for such systems (Eqs.23-24). Note that while our numerical results refer to stellar-mass companions, our analysis is not restricted to any specific binary formation scenarios, and can be easily adapted to other types of systems (e.g. when the tertiary is a supermassive BH) by applying appropriate scaling relations. The key dimensionless parameter that determines the spin-orbit evolution is ${\cal A}_{0}$ (see Eq.18).

The BH binaries detected by aLIGO so far (Abbott et al., 2016b, 2017) have relatively small $\chi_{\rm eff}$ ($0.06^{+0.14}_{-0.14}$ for GW150914, $0.21^{+0.2}_{-0.1}$ for GW151226, and $-0.12^{+0.21}_{-0.30}$ for GW170104). These small values could be due to the slow rotation of the BHs (Zaldarriaga et al., 2017) or spin-orbit misalignments. The latter possibility would imply a dynamical formation channel of the BH binaries (such as exchange interaction in globular clusters (Rodriguez et al., 2015; Chatterjee et al., 2017)) or, as our calculations indicate, dynamical influences of external companions.

BL thanks Natalia Storch, Feng Yuan and Xing-Hao Liao for discussions. This work is supported in part by grants from the National Postdoctoral Program and NSFC (No.BX201600179, No.2016M601673 and No.11703068). DL is supported by NASA grants NNX14AG94G and NNX14AP31G, and a Simons Fellowship in theoretical physics. This work made use of the Resource in the Core Facility for Advanced Research Computing at SHAO.