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

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.


INTRODUCTION
The recent breakthrough in the detection of gravitational waves (GWs) from merging black hole (BH) binaries by advanced LIGO (Abbott et al. 2016a(Abbott et al. ,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. 2010Belczynski et al. , 2016Mandel & 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, (where m 1,2 are the BH masses, a 1,2 = cS 1,2 /(Gm 2 1,2 ) are the dimensionless BH spins, andL is the unit orbital angular momentum vector), which can be measured from the phase evolution of GWs (Abbott et al. 2016b(Abbott et al. , 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 prod-uct 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 ( 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 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).

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 µ in ≡ m 1 m 2 /m 12 , with m 12 ≡ m 1 + m 2 . Similarly, the outer binary has µ out ≡ (m 12 m 3 )/m 123 with m 123 ≡ m 12 + m 3 . The orbital semimajor axes and eccentricities are denoted by a in,out and e in,out , respectively. The orbital angular momenta of the inner and outer binaries are whereL in,out are unit vectors. The relative inclination betweenL in andL out is denoted by I. For convenience, we will frequently omit the subscript "in". The merger time due to GW radiation of an isolated binary with initial a 0 and e 0 = 0 is given by 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: Here the "Lidov-Kozai" (LK) terms are given explicitly in (Liu et al. 2015a) (we also evolve L out and e out ), and the associated timescale of LK oscillation is where n = Gm 12 /a 3 is the mean motion of the inner binary and a out,eff ≡ a out 1 − e 2 out . General Relativity (1-PN correction) induces pericenter precession 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 ( 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 • ). We see that the inner binary undergoes cyclic excursions to maximum eccentricity e 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.

Eccentricity Excitation and Merger Time
In the quadrupole approximation, the maximum eccentricity e max attained in the LK oscillations (starting from an initial I 0 and e 0 0) can be calculated analytically using energy and angular momentum conservation, according to the equation (Anderson et al. 2017a) 3 8 -The maximum eccentricity of the inner BH binary versus the initial inclination I 0 of the tertiary companion, calculated using Equation (9). The inner binary has m 1 = m 2 = 30M , a = 0.1AU, and initial e 0 0. The companion has a circular orbit and its mass and semimajor axis are as labeled. The emax(I 0 ) curve depends mainly on m 3 /a 3 out . The horizontal (e lim ) and vertical (I ± ) lines are given by Equations (11) and (12), respectively.
where j min ≡ 1 − e 2 max , η ≡ (L/L out ) e=0 and ε GR is given by Note that in the limit of η → 0 and ε GR → 0, Equation (9) yields the well-known relation e max = 1 − (5/3) cos 2 I 0 . The maximum possible e max for all values of I 0 , called e lim , is given by and is reached at cos I 0 = (η/10)(4j 2 lim − 5). Eccentricity excitation (e max ≥ 0) occurs within a window of inclinations (cos This window vanishes when ε GR ≥ 9 4 + 3 80 η 2 (no eccentricity excitation). (13) Figure 2 shows some examples of the e max (I 0 ) curves. For η 1, these curves depend mainly on m 3 /a 3 out,eff (for given inner binary parameters). As ε GR increases (with decreasing m 3 /a 3 out,eff ), the LK window shrinks and e lim decreases. Eccentricity excitation is suppressed (for all I 0 's) when Equation (13)  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 lim given by Equation (11) remains valid even for ε oct = 0 (Liu et al. 2015a;Anderson & Lai 2017b). Therefore Equation (13) still provides a good criterion for "negligible eccentricity excitation" (note that when ε oct = 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 m compared to the circular merger time T m,0 (Figure 3). We compute T 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 max , which can be calculated using Equation (9). We see that in general T m /T m,0 can be approximated by (1 − e 2 max ) α , with α 1.5, 2 and 2.5 for e max = (0, 0.6), (0.6, 0.8) and (0.8, 0.95), respectively.

Spin-Orbit Coupling
We now study how the BH spin evolves during the binary merger, considering only S 1 = S 1Ŝ1 (whereŜ 1 is the unit vector). The de Sitter precession ofŜ 1 around L is (Barker & O'Connell 1975) Note that there is a back-reaction torque from S 1 on L; this can be safely neglected since L S 1 . Before presenting our numerical results, it is useful to note the different regimes for the evolution of the spinorbit misalignment angle θ SL (the angle between S 1 and L). In general, the inner binary axisL precesses (and nutates when e = 0) around the total angular momentum J = L + L out (recall that S 1 , S 2 L, and J is constant in the absence of GW dissipation). The related precession rate Ω L is of order t −1 LK at e ∼ 0 (Equation 17), but increases with e. Depending on the ratio of Ω dS and Ω L , we expect three possible spin behaviors: (i) For Ω L Ω dS ("nonadiabatic"), the spin axisŜ 1 cannot "keep up" with the rapidly changingL, and thus effectively precesses aroundĴ, keeping θ SJ ≡ cos −1 (Ŝ 1 ·Ĵ) constant. (ii) For Ω dS Ω L ("adiabatic"),Ŝ 1 closely "follows"L, maintaining an approximately constant θ SL . (iii) For Ω dS ∼ Ω L ("trans-adiabatic"), the spin evolution can be chaotic due to overlapping resonances. Since both Ω dS and Ω 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 ofL is governed by the equation whereĴ is the unit vector along J = L + L out , and We can define an "adiabaticity parameter" As the binary orbit decays, the system may transition from "non-adiabatic" (A 1) at large a's to "adiabatic" (A 1) at small a's, where the final spin-orbit misalignment angle θ f SL is "frozen". Note that A is directly related to ε GR by Thus, when the initial value of ε GR (at a = a 0 ) satisfies ε GR,0 9/4 (a necessary condition for LK eccentricity excitation; see Equation 13), we also have A 0 (3m 2 + µ)/(2m 12 ) ∼ 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, θ SL evolves in a chaotic way, with the final value θ f SL depending sensitively on the initial conditions (the precise initial θ SL and I 0 ). Also note that retrograde spin (θ f SL > 90 • ) can be produced even though the binary always remains prograde with respect to the outer companion (I < 90 • ). 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 θ f SL for various triple parameters. In our "population synthesis" study, we consider initial conditions such that S 1 is parallel toL, the binary inclinations are isotropically distributed (uniform distribution in cos I 0 ), and the orientations of e and e out (for systems with e out = 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 out = 0 (so that the octupole effect vanishes), and several values of a out . We see that when the eccentricity of the inner binary is excited (I 0 lies in the LK window), a wide range of θ f SL is generated, including appreciable fraction of retrograde (θ f SL > 90 • ) systems (see the a out = 3 AU case, for which e lim = 0.84). The "memory" of chaotic spin evolution is evident, as slightly different initial inclinations lead to vastly different θ f SL . The regular behavior of θ f SL around I 0 = 90 • (again for the a out = 3 AU case) is intriguing, but may be understood using the theory developed in Ref. . For systems with no eccentricity excitation, θ f 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 = m 2 and e out = 0, for which octupole terms may affect the orbital evolution. We see that eccentricity excitation and the corresponding reduction in T m occur outside the analytic (quadupole) LK window (see the e out = 0.8 case, for which e lim = 0.66). As in the equal-mass case (Figure 4), a wide range of θ f SL values are produced whenever eccentricity excitation occurs. A larger fraction (23%) of systems attain retrograde spin (θ f SL > 90 • ). Again, for systems with negligible eccentricity excitation, θ f SL 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 θ f 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 thatL rotates around theĴ axis at the rate (−Ω L ). In this rotating frame, the spin evolution equation (15) transforms to Note that in the absence of GW dissipation,L andL out are constants (in the rotating frame), and thusŜ 1 precesses with a constant Ω eff . The relative inclination between Ω eff andL is given by The inner binary has fixed initial a 0 = 0.1 AU and e 0 = 0.001, and eout = 0 for the outer binary. In the top panel, the dots are the result of numerical integration for the aout = 3 AU system (a total of 400 runs on a uniform cos I 0 grid), and the dashed curves are the analytical results for circular orbits, as given by Equations (23) and (24). The initial value of the adiabaticity parameter A 0 (Equation 18 with a = a 0 ) is also given. The vertical lines (I ± ) shown in the middle panel correspond to the LK window of eccentricity excitation (Equation 12). The bottom panel shows the distribution of the final spin-orbit misalignment angle for the system with aout = 3AU. Now if we include GW dissipation,L ·L out = cos I is exactly conserved, and Ω eff becomes a slowly changing vector. When the rate of change of Ω eff is much smaller than |Ω eff |, the angle between Ω eff andŜ 1 is adiabatic invariant, i.e. θ eff,S1 constant (adiabatic invariant).
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, A 0 (3m 2 + µ)/(2m 12 ) ∼ 1 (see Equations 13 and 19), and thus θ f SL θ 0 eff,L 20 • , i.e., only modest spin-orbit misalignment can be generated. For systems with A 0 1 (e.g., very distant companion), we have θ f SL 1. The above analytical result can be easily generalized to the situation of non-zero initial spin-orbit misalignment. It shows that θ f SL θ 0 SL for A 0 1.

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 θ f SL depending sensitively on the initial conditions. A wide range of θ f SL (including θ f SL > 90 • ) can be produced from an initially aligned (θ 0 SL = 0) configuration (see . 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 ( 20 • ) spin-orbit misalignment can be produced -we have derived an analytic expression for θ f SL 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 A 0 (see Eq.18).
The BH binaries detected by aLIGO so far (Abbott et al. 2016b(Abbott et al. , 2017 have relatively small χ 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 spinorbit 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.