Precessional dynamics of black hole triples: binary mergers with near-zero effective spin

Introduction

Since the first detection of merging black hole (BH) binaries, there has been a proliferation of astrophysical models for producing such systems. To narrow the range of possibilities, it has been shown that the misalignment between the binary’s orbital angular momentum and the BH spins can serve as an important discriminant between these different formation channels (e.g., Rodriguez et al. 2016a; Farr et al. 2017a; Farr et al. 2017b). Here, we consider binary BH mergers formed through the evolution of stellar triples in the field. In this scenario, the binary is driven to merger by the presence of a third BH companion. This dynamical configuration can induce high eccentricities in the inner BH binary via the Lidov-Kozai (LK) mechanism (Kozai 1962; Lidov 1962), eliminating the need for the common-envelope phase often invoked in standard binary evolution models (Antonini et al. 2014; Silsbee & Tremaine 2017; Antonini et al. 2017). We specifically consider the secular evolution of the effective spin parameter, $\chi_{\rm eff}$ (see Eq. 6), the combination of BH spins best measured by current gravitational-wave (GW) detectors (Abbott et al. 2016b; Abbott et al. 2016a; Abbott et al. 2017b).

To date, the BH binaries detected by LIGO/VIRGO have all exhibited small $\chi_{\rm eff}$, with all but one–GW151226–being consistent with $\chi_{\rm{eff}}=0$ (Abbott et al. 2016a; Abbott et al. 2017a; Abbott et al. 2017b). The individual component spins in isolated field binaries are expected to be nearly perpendicular to the binary orbital plane (Kalogera 2000), meaning that current measurements would imply slow rotation of the BHs (Belczynski et al. 2017, if these are formed in the field;). Low spins, however, are in contrast with current observational constraints and theoretical expectations, both favouring rapid BH rotation at formation (Gammie et al. 2004; Miller et al. 2011). In this letter, we show that a more consistent explanation is possible: the small values of $\chi_{\rm eff}$ are a consequence of the spin-orbit tilt produced by the binary’s long-term interaction with a distant companion.

Orbit-averaged Equations

We consider a BH binary with total mass $M=m_{1}+m_{2}$, semi-major axis $a$ and eccentricity $e$, orbited by a tertiary BH with mass $m_{3}$ on an outer orbit with semi-major axis $a_{\rm out}$ and eccentricity $e_{\rm out}$. We work in terms of the dimensionless inner-orbit angular-momentum vector, $\bm{j}={\sqrt{1-e^{2}}\bm{\hat{j}}}$, and the eccentricity vector, ${\bm{e}}=e\bm{\hat{e}}$, defined in Jacobi coordinates. We define the circular angular momenta for the inner and outer orbits as $L=\mu\sqrt{GMa}$ and $L_{\rm out}=\mu_{\rm out}\sqrt{G(M+m_{3})a_{\rm out}}$ respectively, with $\mu=m_{1}m_{2}/M$, and $\mu_{\rm out}=Mm_{3}/(M+m_{3})$. Finally, the total angular momentum of the triple system is $\bm{J}=L\bm{j}+L_{\rm out}\bm{j}_{\rm out}$, where $\bm{j}_{\rm out}$ is the outer-orbit angular-momentum vector. In the absence of dissipation, $\bm{J}$ is constant.

A tertiary companion highly inclined with respect to the binary orbit by an angle, denoted by $I$, can induce large amplitude, periodic oscillations in the inner binary eccentricity (see Naoz 2016, and references therein). We describe this evolution using the secular equations at the octupole level of approximation (Liu et al. 2015; Petrovich 2015). We also add the 1 and 2.5 post-Newtonian (pN) terms, describing the (Schwarzschild) precession of the argument of periapsis and the orbital decay due to gravitational-wave emission respectively (Peters 1964). The evolution of the inner binary orbit is determined by the set of equations:

where $\nu=\sqrt{GM/a^{3}}$ and the LK terms are given explicitly by Eq. (17)-(20) in Liu et al. 2015. The associated timescale of the LK oscillations is (Antognini 2015, e.g.,)

Finally, we add the spin-orbit interaction terms (Schnittman 2004, e.g.,):

with $\bm{S}_{1}$ and $\bm{S}_{2}$ the spins of $m_{1}$ and $m_{2}$ respectively. The spin vectors can also be written in terms of the dimensionless spin parameter $\bm{\chi}$ as $\bm{S}_{1(2)}=\bm{\chi}_{1(2)}Gm_{1(2)}^{2}/c$, with $|\bm{\chi}_{1(2)}|\leq 1$. We do not account for the back-reaction torque from $\bm{S}$ on $\bm{L}$, as well as the spin-spin precessional terms. Because these terms depend on the first power of the spin-angular momentum and $jL\gg S$ during the LK oscillations, they can be safely neglected.

When describing our results, we will often refer to the evolution of the angle between the two spin vectors and the spin-orbit misalignment using $\theta_{\rm ss}=\cos^{-1}\bm{{\hat{S}}}_{1}\cdot\bm{{\hat{S}}}_{2}$ and $\theta_{\rm 1(2)}=\cos^{-1}\bm{\hat{S}}_{\rm{1(2)}}\cdot\bm{\hat{j}}$ respectively. Similarly, the misalignment between the spins and the total angular momentum of the triple is given by $\Theta_{\rm 1(2)}=\cos^{-1}\bm{\hat{S}}_{\rm{1(2)}}\cdot\bm{\hat{J}}$. Finally, we define the binary effective spin parameter as

In the following sections, unless otherwise specified, we work under the assumption that $\bm{S}_{1}$, $\bm{S}_{2}$ and $\bm{j}$ are initially nearly aligned with each other. Spin-orbit alignment is generally observed in solitary (solar-type) binaries with relatively large separation (Hale 1994, $\lesssim 40\rm AU$; ) and could either be primordial in origin (Corsaro et al. 2017) or produced through tidal evolution and mass transfer during stellar evolution prior to BH formation (Kalogera 2000, e.g.,). We note that while a significant spin-orbit tilt has been measured in a number of cases (Triaud et al. 2013; Albrecht et al. 2013; Albrecht et al. 2014), such misalignment is typically attributed to long-term triple evolution (Naoz & Fabrycky 2014, e.g.,).

Suppression of chaos

To the lowest pN order, each of the two BH spin vectors precess in response to torques from the binary at the orbit-averaged rate (Apostolatos et al. 1994, e.g.,):

For fixed $\bm{j}$, the precession described by Eq. (5) has the form of uniform precession of the spin vector about the binary angular momentum vector. At the same time, during the LK oscillations, the binary angular momentum vector precesses around the total angular momentum of the triple system at a rate: $\Omega_{\rm LK}\approx{t_{\rm LK}^{-1}}\ .$

There are three possible regimes for the evolution of each of the two BH spins (Storch et al. 2014; Storch & Lai 2015; Liu & Lai 2017; Lai et al. 2018): (i) for $R_{1(2)}\gg 1$ (adiabatic regime), the spin follows $\bm{\hat{j}}$ adiabatically, maintaining an approximately constant spin-orbit alignment angle $\theta_{\rm 1(2)}$ and consequently a constant $\chi_{\rm eff}$; (ii) for $R_{1(2)}\ll 1$ (non-adiabatic regime), $\bm{\hat{S}}_{1(2)}$ effectively precesses about $\bm{\hat{J}}$, maintaining an approximately constant angle $\Theta_{\rm 1(2)}$; and (iii) for $R_{1(2)}\approx 1$ (trans-adiabatic regime), the spin precession rate matches the orbital precession rate and the evolution of the spin orbit orientation can become more complicated.

We might expect that at $R\approx 1$ the spin-orbit angle will exhibit chaotic evolution due to overlapping resonances, possibly leading to a wide range of final spin-orbit angles (Liu & Lai 2017). As discussed next, we find that such chaotic behaviour is suppressed due to the 1pN precession of the periapsis.

To lowest order, the Schwarzschild contribution is $\omega_{\rm 1pN}=\frac{3GM}{c^{2}aj^{2}}\nu.$ By setting $\omega_{\rm 1pN}/\pi=|dj/dt|/j\approx t_{\rm LK}^{-1}/j$, we find the critical angular momentum below which LK oscillations are strongly quenched by relativistic precession:

At $j\leq j_{\rm GR}$, LK oscillations are damped by the in-plane precession caused by the 1pN terms. An approximate criterion for Schwarzschild precession to fully quench the LK oscillations can be obtained by setting $j_{\rm GR}\geq 1$ in the previous equation, which gives:

See also Blaes et al. 2002 and Silsbee & Tremaine 2017 for similar derivations. We expect that for systems which satisfy this condition the spin-orbit dynamics leading to large misalignment will also be somewhat suppressed.

Fig. 1 shows bifurcation diagrams giving for each value of $R$ (or $a_{\rm out}$) the corresponding value of $\chi_{\rm eff}$ at every eccentricity maximum ($\chi_{\rm eff;max}$) during 100 LK oscillations. When $R$ is near unity, the 1pN apsidal precession terms affect the evolution of $\chi_{\rm eff;max}$ in important ways. If these terms are not included in our calculations (bottom panels), $\chi_{\rm eff;max}$ can attain large values and its evolution is often chaotic as illustrated by the high degree of scatter for a single value of $R$. Such chaotic behaviour is well known (Storch et al. 2014; Storch & Lai 2015; Anderson et al. 2016), and corresponds to the trans-adiabatic regime discussed above. However, when all 1pN terms are added, $\chi_{\rm eff;max}$ remains close to unity at all $R\gtrsim 0.1$. We conclude that the chaotic spin-orbit dynamics seen in the bottom panels, is effectively suppressed by the 1pN relativistic precession of the inner BH binary orbit. This is a consequence of the fact that $\Omega_{\rm S}\approx\omega_{\rm 1pN}$ for any $a$ and $j$, so that when $R\approx 1$ the inequality Eq. (10) is also satisfied. In this situation, the LK oscillations are quenched and only modest misalignment can be induced.

The results displayed in Fig. 1 indicate (and the numerical simulations below confirm) that only when $R\ll 1$ the spin-orbit orientation does change significantly while the binary simultaneously experiences extreme eccentricity excitation that can lead to its coalescence through energy loss by GW radiation.

Population synthesis model

In this section we derive the spin-orbit misalignment of binary BHs driven to a merger by a distant BH companion using a population synthesis approach. We evolved massive stellar triples to BH triples using a modified form of the Binary Stellar Evolution (BSE) package (Hurley et al. 2002). The tertiary star was evolved simultaneously using the single stellar evolution subset of BSE.

We employed 11 different stellar metallicities, logarithmically-spaced from $1.5Z_{\odot}$ to $0.01Z_{\odot}$. We sampled the primary star mass for the inner binary from a Kroupa initial mass function (Kroupa & Weidner 2003) in the range $22\leq m\leq 150\ M_{\odot}$. The masses of the secondary and tertiary stars were assigned by assuming flat mass ratio ($m_{2}/m_{1}$ and $m_{3}/(m_{1}+m_{2})$) distributions between $0$ and $1$. We take $N\propto\log\left(P_{1}/{\rm days}\right)^{-0.55}$ for the inner orbital period, and take the outer orbital periods to be flat in log-space with $a_{\rm out}\leq 10^{5}\rm AU$. The eccentricities of the inner binary was drawn from a $N\propto e^{-0.42}$ distribution with $e$ from 0 to 0.9, while the eccentricity of the outer orbit was drawn from a thermal distribution, $N\propto 2e$. Our choice of initial conditions is consistent with observations of nearby young clusters and associations (Sana et al. 2012, e.g.,). All the angles defining the triple (arguments of periapsis, longitudes of the ascending node, and the inclination) were drawn from isotropic distributions.

During the main-sequence evolution we followed the changes to the initial orbital properties due to the mass loss and super-nova kicks (Fryer et al. 2012, based on) experienced by each component of the triple during the formation of a BH (Toonen et al. 2016, e.g.,). Our prescriptions for mass-loss, stellar winds and natal kicks are identical to those in Rodriguez et al. 2016b, and include the latest prescriptions for the pulsational pair-instability in massive stars (Belczynski et al. 2016). We reject any systems for which either the inner or outer binaries collide or in which the triple becomes secularly unstable at any point (Mardling & Aarseth 2001).

The distributions of initial conditions for the stellar progenitors and for the BH triples produced by our models are displayed in Fig. 2. For two metallicities, we show in these figures the distribution of masses, eccentricities and semi-major axes for the bound BH triple systems that are hierarchically stable, and are not dominated by Schwarzschild precession (i.e., triples with $j<j_{\rm GR}$).

Our methodology does not consider the possibilities of mass accretion between the inner binary and the tertiary, or any dynamical interaction between the inner and outer binaries. Such physics, while interesting, is significantly beyond the scope of this letter (though see Antonini et al. 2017, for an analysis of such triples using a self-consistent method)

Finally, we assume that the spin vectors of the BHs at the time of their formation are aligned with the stellar progenitor spins, and that these progenitor spins are initially aligned with $\bm{{\hat{j}}}$. Although our models include the effect of natal kicks on the orientation of the orbits, we find that 97% (99%) of our initial population of BH triples have post-kick spin-orbit misalignments less than $0.1^{\circ}$ ($6^{\circ}$).

After generating initial conditions for the BH triples, we integrated them forward in time up to a maximum time of $13.8\rm Gyr$, or until the binary peak GW frequency became larger than $10\rm Hz$. Fig. 3 shows the distribution of $\chi_{\rm{eff}}$ for BH binary mergers produced by the LK mechanism in our simulations (which we define to be those that would not have merged in less than $13.8\rm Gyr$ if evolved as isolated binaries). In order to obtain the $\chi_{\rm{eff}}$ distribution in Fig. 3 we conservatively assumed that both BHs were maximally spinning, but note that this distribution can be trivially rescaled to any $\chi_{1}=\chi_{2}$; the resulting $\chi_{\rm{eff}}$ distribution is peaked around zero with $\left|\chi_{\rm eff}\right|<\chi_{1}\times 0.5(0.25)$ for $\sim 70\%(40\%)$ of the sources. The value of $\chi_{\rm eff}$ for the BH binaries detected by Advanced LIGO/VIRGO are also displayed in the figure, showing that they are clustered in the range of values $\left|\chi_{\rm eff}\right|\lesssim 0.5$, near the peak of our synthetic distribution. The bottom panel of Fig. 3 shows that $\chi_{\rm eff}$ is clustered around zero for initial conditions which are moderately adiabatic, $R\gtrsim 10^{-3}$, demonstrating that the LK process is able to drive $\chi_{\rm eff}$ towards zero for a large portion of parameter space for triples. In our models, most mergers are produced at metallicities $\lesssim 0.25Z_{\odot}$. In fact, the number of mergers is increased by a factor $\sim 100$ at these lower metallicities primarily due to the reduced BH natal kicks and stellar mass loss which increase the chance for a triple to remain bound prior to BH formation. Finally, we find that all our systems become strongly adiabatic (i.e., $R\gg 1$), after which point $\chi_{\rm eff}=const.$, before the $10$Hz frequency band is reached. A consequence of this is that all binaries suffer substantial circularisation and have $e<0.1$ by the time they enter the LIGO/VIRGO frequency window.

What is the origin of the near-zero peak of the $\chi_{\rm eff}$ distribution in Fig. 3? This peak could be due to two processes: (i) the differential precession of the two spin vectors nearly randomise their relative orientation with respect to each other and with respect to $\bm{j}$ (upper panel in Fig. 4); (ii) both spin-orbit angles evolve individually towards $\pi/2$ as the orbit slowly decays by GW emission (lower panel in Fig. 4). In the former case, $\chi_{\rm eff}$ will peak around zero because $\cos\theta_{1}$ and $\cos\theta_{2}$ have uniform and independent distributions. In case (ii), $\chi_{\rm eff}$ will peak around zero because $\cos\theta_{1}\approx\cos\theta_{2}\approx 0$. Case (ii) turns out to be more important. This is shown in Fig. 3, where we see that the final distributions of $\cos\theta_{1}$ and $\cos\theta_{2}$ (not shown) follow closely the distribution of $\chi_{\rm eff}$ and they are also peaked around $0$ The reason for this “attractor” towards $\pi/2$ has been recently identified in Liu & Lai 2018 while this letter was under review..

Conclusions

For BH binary mergers produced from isolated field binaries, it is expected that the individual BH spins should be nearly aligned with the binary angular momentum (Kalogera 2000). In these standard population synthesis models, the relatively small $\chi_{\rm eff}$ of the BH binaries detected so far by LIGO/VIRGO (Abbott et al. 2016b; Abbott et al. 2016a; Abbott et al. 2017b) is more easily explained if the spin magnitudes were nearly zero. This, however, appears to be disfavored by current observational constraints and theoretical models which suggest finite spins for BHs at birth (Gammie et al. 2004; Miller et al. 2011). A triple origin instead, could provide a more consistent explanation for the low $\chi_{\rm eff}$ values (as shown here), as well as for their merger rates (Rodriguez & Antonini 2018). We note, however, that to be a viable explanation for all the LIGO/VIRGO detections, our mechanism would require strong suppression of other channels.

Software: the secular code used in this paper is available at https://github.com/carlrodriguez/kozai.

FA acknowledges support from an STFC E. Rutherford fellowship (ST/P00492X/1), CR from a Pappalardo fellowship at MIT, CP from the J. L. Bishop Fellowship and from the Gruber Foundation Fellowship, CF from a Sir Edward Youde Memorial Fund scholarship and the Undergraduate Research Opportunities Program at MIT. FA and CR acknowledge the hospitality and support of the Aspen Center for Physics (NSF Grant PHY-1607611).

References