Black hole and neutron star mergers in Galactic Nuclei: the role of triples

TL;DR

This study explores how triple systems of black holes and neutron stars near supermassive black holes in galactic nuclei can significantly enhance merger rates, potentially explaining observed gravitational wave events.

Contribution

It introduces a new scenario involving triple compact objects in galactic nuclei, showing they can increase merger rates compared to binary systems.

Findings

Merger fraction up to 5-8 times higher in triples than binaries.

Triples can merge on shorter timescales due to dynamical interactions.

Triples could contribute significantly to LIGO/VIRGO gravitational wave detections.

Abstract

Nuclear star clusters that surround supermassive black holes (SMBHs) in galactic nuclei are thought to contain large numbers of black holes (BHs) and neutron stars (NSs), a fraction of which form binaries and could merge by Kozai-Lidov oscillations (KL). Triple compact objects are likely to be present, given what is known about the multiplicity of massive stars, whose life ends either as a NS or a BH. In this paper, we present a new possible scenario for merging BHs and NSs in galactic nuclei. We study the evolution of a triple black hole (BH) or neutron star (NS) system orbiting an SMBH in a galactic nucleus by means of direct high-precision $N$-body simulations, including Post-Newtonian terms. We find that the four-body dynamical interactions can increase the KL angle window for mergers compared to the binary case and make BH and NS binaries merge on shorter timescales. We show that the merger fraction can be up to $\sim 5$--$8$ times higher for triples than for binaries. Therefore, even if the triple fraction is only $\sim 10\%$--$20\%$ of the binary fraction, they could contribute to the merger events observed by LIGO/VIRGO in comparable numbers.

Figures

Tables

Table 1. Table 1: Description of important quantities used in the text.

Symbol	Description
$M_{\mathrm{SMBH}}$	Mass of the SMBH
$m_{11}$	Mass of the first component of the inner binary in the BH/NS triple
$m_{12}$	Mass of the second component of the inner binary in the BH/NS triple
$m_3$	Mass of the third companion in the BH/NS triple
$a_{\mathrm{out}}$	Semi-major axis of the orbit of the triple around the SMBH
$e_{\mathrm{out}}$	Eccentricity of the orbit of the triple around the SMBH
$a_{12}$	Semi-major axis of the inner orbit in the BH/NS triple
$e_{12}$	Eccentricity of the inner orbit in the BH/NS triple
$a_3$	Semi-major axis of the outer orbit in the BH/NS triple
$e_3$	Eccentricity of the outer orbit in the BH/NS triple

Table 2. Table 2: Models: triple type, slope of the BH mass function ( β 𝛽 \beta ), slope of the distribution of a out subscript 𝑎 out {a_{\rm out}} ( α 𝛼 \alpha ), distribution of e 12 subscript 𝑒 12 {e_{12}} and e 3 subscript 𝑒 3 {e_{3}} ( f ​ ( e ) 𝑓 𝑒 f(e) ), maximum outer semi-major axis of the triple ( a 3 , max subscript 𝑎 3 max {a_{\rm 3,max}} ), merger fraction from the N 𝑁 N -body simulations ( f merge subscript 𝑓 merge f_{\rm merge} ). For comparison, the second part of the table summarises the results of Fragione et al. ( 2018 ) .

Triple Type	$\beta$	$\alpha$	$f\left(e\right)$	$a_{3,\max }$ (AU)	$f_{\mathrm{merge}}$
BH	$1$	$2$	uniform	$50$	$0.287$
BH	$1$	$2$	uniform	$200$	$0.342$
BH	$1$	$2$	uniform	$1000$	$0.346$
BH	$1$	$2$	thermal	$50$	$0.358$
BH	$2$	$2$	uniform	$50$	$0.256$
BH	$3$	$2$	uniform	$50$	$0.233$
BH	$4$	$2$	uniform	$50$	$0.276$
BH	$1$	$3$	uniform	$50$	$0.304$
NS	-	$2$	uniform	$50$	$0.140$
NS	-	$2$	uniform	$200$	$0.158$
NS	-	$2$	uniform	$1000$	$0.152$
NS	-	$2$	thermal	$50$	$0.188$
NS	-	$1.5$	uniform	$50$	$0.141$
Binary Type	$\beta$	$\alpha$	$f\left(e\right)$	$a_{3,\max }$ (AU)	$f_{\mathrm{merge}}$
BH	$1$	$2$	uniform	-	$0.045$
BH	$2$	$2$	uniform	-	$0.050$
BH	$3$	$2$	uniform	-	$0.036$
BH	$4$	$2$	uniform	-	$0.041$
BH	$1$	$3$	uniform	-	$0.051$
NS	-	$2$	uniform	-	$0.028$
NS	-	$1.5$	uniform	-	$0.032$

Equations

$$T_{\mathrm{KL}}=\frac{8}{15\pi }\frac{m_{\mathrm{tot}}}{m_{3b}}\frac{P_{3b}^2}{P_{\mathrm{bin}}}{\left(1-e_{3b}^2\right)}^{3/2},$$

$$e_{in,max}=\sqrt{1-\frac{5}{3}{\cos }^2{i}_0}.$$

$$T_{\mathrm{oct}}=\frac{1}{ϵ}{T}_{\mathrm{LK}},$$

$$ϵ=\frac{m_{\mathrm{bin},1}-m_{\mathrm{bin},2}}{m_{\mathrm{bin},1}+m_{\mathrm{bin},2}}\frac{a_{\mathrm{bin}}}{a_{3b}}\frac{e_{3b}}{1-e_{3b}^2},$$

$$T_{GR}=\frac{a_{\mathrm{bin}}^{5/2}{c}^2(1-e_{\mathrm{bin}}^2)}{3G^{3/2}(m_{\mathrm{bin},1}+m_{\mathrm{bin},2}{)}^{3/2}}.$$

$$T_{\mathrm{KL}}^{12M}$$

$$T_{\mathrm{KL}}^{123M}$$

$$T_{\mathrm{KL}}^{123}$$

$${\tau }_{1+3}$$

$$\frac{E_t}{m_{\mathrm{t}}{\sigma }^2}\lesssim 1,$$

$$T_{\mathrm{EV}}$$

$$T_{\mathrm{VRR}}=\frac{\pi a_{\mathrm{out}}^{3/2}}{G^{1/2}}\frac{M_{\mathrm{SMBH}}^{1/2}}{m_{\ast }\sqrt{N_{\ast }(a_{\mathrm{out}})}},$$

$$T_{\mathrm{SRR}}=\frac{\pi a_{\mathrm{out}}^{3/2}}{2G^{1/2}}\frac{M_{\mathrm{SMBH}}^{1/2}}{m_{\ast }},$$

$$\frac{dN}{dm}\propto M^{-\beta }$$

$$\frac{R_{\mathrm{p},\mathrm{per}}}{a_{\mathrm{in}}}\ge 2.8{\left[\left(1+\frac{m_{\mathrm{per}}}{m_{\mathrm{in},\mathrm{tot}}}\right)\frac{1+e_{\mathrm{per}}}{\sqrt{1-e_{\mathrm{per}}}}\right]}^{2/5}\left(1.0-0.3\frac{I}{\pi }\right),$$

$${\ddot{\text{r}}}_i=-G\sum _{j\ne i}\frac{m_j({\text{r}}_i-{\text{r}}_j)}{{\left|{\text{r}}_i-{\text{r}}_j\right|}^3},$$

$$T=\min (T_{\mathrm{EV}},1\mathrm{Myr}).$$

$${\mathcal{R}}_{12M}=\frac{T_{KL}^{123}}{T_{KL}^{12M}}$$

$${\mathcal{R}}_{123M}=\frac{T_{KL}^{123}}{T_{KL}^{123M}}$$

$${\mathcal{R}}_M=\frac{T_{KL}^{123M}}{T_{KL}^{12M}},$$

$$f_{\mathrm{GW}}=\frac{\sqrt{G(m_{11}+m_{12})}}{\pi }\frac{(1+e_{12}{)}^{1.1954}}{[a_{12}(1-e_{12}^2){]}^{1.5}}.$$

Full text

1 Introduction

The LIGO-Virgo collaboration has recently released a catalogue of compact object (CO) mergers due to gravitational wave (GW) emission, comprised of ten merging black hole (BH-BH) binaries and one merging neutron star (NS-NS) binary (The LIGO Scientific Collaboration & the Virgo Collaboration, 2018). Studying the possible mechanisms that lead to BH and NS mergers is currently an active area of research. Several scenarios for their origins have been proposed. The possibilities include isolated binary evolution through a common envelope phase (Belczynski et al., 2016), chemically homogeneous evolution in short-period stellar binaries (Mandel & de Mink, 2016; Marchant et al., 2016), stellar triples (Antonini et al., 2017; Silsbee & Tremaine, 2017), quadruple systems (Fragione & Kocsis, 2019; Liu & Lai, 2019) and Kozai-Lidov (KL) mergers of binaries in galactic nuclei (Antonini & Perets, 2012; Hamers et al., 2018; Hoang et al., 2018; Fragione et al., 2018; Grishin et al., 2018). These mergers could either occur in active galactic nuclei (AGN) disks (e.g. McKernan et al., 2018; Secunda et al., 2018) or in star clusters (Askar et al., 2017; Banerjee, 2018; Fragione & Kocsis, 2018; Rodriguez et al., 2018). While these models typically predict roughly the same rates ($\sim\ \mathrm{few}$ Gpc*-3* yr*-1*), their respective contributions could potentially be disentagled using the observed eccentricity, spin, mass and redshift distributions (see e.g. O’Leary et al., 2016; Samsing et al., 2018; Zevin et al., 2018; Perna et al., 2019).

The wide literature on dynamically-induced mergers has focused on GW signals generated from merging BH-BH binaries in open and globular clusters. However, galactic nuclei typically have total stellar masses comparable to the combined stellar masses of open cluster systems, and a significant fraction of globular clusters in a given galaxy (Böker et al., 2004; Côté et al., 2006). Thus, galactic nuclei represent a promising environment where CO binaries can form dynamically, then harden and finally merge due to the emission of GW radiation (O’Leary, Kocsis & Loeb, 2009; Antonini & Rasio, 2016; Petrovich & Antonini, 2017; Fragione, 2019; Rasskazov & Kocsis, 2019).

Compared to open and globular clusters, studies focused on galactic nuclei have generally either concentrated on how large-scale nuclear star cluster dynamics ($\gtrsim 1$ pc) can drive CO binaries to merge within a Hubble time (Antonini & Rasio, 2016; Petrovich & Antonini, 2017; Gondán et al., 2018), or on how these binaries merge due to eccentric KL oscillations imparted from the enormous SMBH potential ($\lesssim 1$ pc; Antonini & Perets, 2012; Hoang et al., 2018; Fragione et al., 2018; Hamers et al., 2018).

Triple stars are likely not rare in open clusters and possibly even globular clusters, in particular for massive stars whose lives end either as a NS or a BH (Duquennoy & Mayor, 1991; Raghavan et al., 2010; Sana et al., 2013; Tokovinin, 2014a, b; Dunstall et al., 2015; Sana, 2017). Little attention has been devoted to the role of merging binary COs in triple systems in galactic nuclei, where they can orbit a central SMBH. Previous studies have primarily focused on determining the merger rates from isolated triple BHs (Antonini et al., 2017; Silsbee & Tremaine, 2017) or triple BHs in globular clusters (Antonini et al., 2016). Compared to binary systems, triples have an additional degree of freedom and KL oscillations can be imparted to the inner binary either by the SMBH or by the third companion (Lidov, 1962; Kozai, 1962). The interplay between these three KL mechanisms working together can give rise to a rich and complicated dynamical evolution for the inner binary (Grishin et al., 2018). This could potentially lead to a larger number of mergers, thtat can be further enhanced by possible resonances (Hamers & Lai, 2017).

In this paper, we propose a new scenario for merging BHs and NSs in galactic nuclei. For the first time, we study the evolution of a triple composed of COs orbiting an SMBH in a galactic nucleus and compare it to the CO binary case (Fragione et al., 2018). We consider a four-body system consisting of a triple BH-BH-BH/NS-NS-NS and an SMBH (see Figure 1 and Table 1). We make a systematic and statistical study of these systems by means of direct high-precision $N$-body simulations, including Post-Newtonian (PN) terms up to order PN2.5. Compared to secular evolution, our approach allows us to correctly follow the large eccentricity values attained by CO binaries via KL cycles (Antognini et al., 2014).

The paper is organized as follows. In Section 2, we discuss what is to date known about the multiplicity properties of BHs and NSs in galactic nuclei and star clusters. In Section 3, we discuss the relevant timescales for the systems considered in this paper. In Section 4, we present our numerical methods to determine the rate of BH and NS mergers in triples orbiting massive black holes, for which we discuss the results. Finally, in Section 5, we discuss the implications of our findings for future observations of mergers involving BHs and NSs, and draw our conclusions.

2 Multiplicity

In this section we briefly review what is known from the literature about multiplicity (i.e., binaries and triples) in the Galactic field, open clusters, globular clusters and nuclear star clusters.

2.1 The Galactic field

Raghavan et al. (2010) updated the seminal work of Duquennoy & Mayor (1991) by selecting from the Hipparcos catalogue a volume-limited sample of field solar-type primary stars in the solar neighbourhood. The authors find that the observed fractions of objects that are single, double, triple and higher-order systems are, respectively, 56 $\pm$ 2, 33 $\pm$ 2, 8 $\pm$ 1 and 3 $\pm$ 1 per cent. Subsequent work showed that most triples composed of low-mass stars tend to be approximately co-planar, whereas the inner and outer orbits of triples composed of massive stars tend to be more inclined relative to each other (e.g. Moe & Kratter, 2018).

2.2 Open clusters

The young star-forming region Taurus-Auriga has an exceedingly low density, such that dynamical interactions involving multiples should be very rare, if they occur at all. In fact, Taurus-Auriga is the only known open cluster with a triple fraction higher than that observed in the Galactic field (Raghavan et al., 2010; Leigh & Geller, 2013). Kraus et al. (2011) performed a high-resolution imaging study to characterize the multiple-star populations in Taurus-Auriga. They found that between $\sim$ $2/3$–$3/4$ of all targets are multiples composed of at least two stars. Thus, only $\sim$ $1/4$–$1/3$ of their objects are single stars.

More nearby open clusters (OCs) are also known to have comparably high multiplicity fractions (see Leigh & Geller (2013) for a more detailed review of the following studies). The Hyades (Patience et al., 1998), Pleiades (Mermilliod et al., 1992; Bouvier et al., 1997) and Praesepe (Mermilliod & Mayor, 1999; Bouvier et al., 2001) have binary fractions of, respectively, 35%, 34% and 40%, and triple fractions of, respectively, 6%, 3% and 6%.

Finally, old open clusters, with masses and densities typically higher than the above OCs, have also been observed to harbour high binary fractions. For example, for the old ($\sim$ 4 Gyr) OC M67, Fan et al. (1996) observed a binary fraction of 45%. For comparison, the old ($\sim$ 7 Gyr) OC NGC 188 has an observed multiplicity fraction of 27% (Geller et al., 2009). Following Leigh & Geller (2013), if we adopt the observed value for the ratio between the fractions of binaries and triples from Latham & Milone (1996), or $f_{\rm t}$/$f_{\rm b}$ $\sim$ 0.1, we infer triple fractions for M67 and NGC 188 of, respectively, 5% and 2%.

2.3 Globular clusters

In globular clusters (GCs), the situation is quite different. Due to their much higher central densities compared to OCs, the multiple star populations in massive GCs can most efficiently be studied using photometry. The pioneering study of Milone et al. (2012) analyzed the main-sequence (MS) binary populations in a sample of 59 GCs. The authors find photometric binary fractions ranging from less than a percent to a few tens of percent, and a previously reported (Sollima et al., 2007) anticorrelation between the binary fraction and the total cluster mass. Sollima (2008) argued that this trend can arise assuming an universal initial binary fraction combined with the dynamical disruption of binaries by (mostly) interloping single stars. Here, the disruption of soft binaries in the cluster core generates the observed anticorrelation, combined with the evaporation of single stars from the cluster outskirts which serves to increase the binary fraction throughout the cluster in preferentially low-mass clusters evaporating at the fastest rates (e.g. Fregeau et al., 2009).

Other than identifying binaries photometrically above and/or to the red of the main sequence in the cluster colour-–magnitude diagram (e.g. Milone et al., 2012), binaries can also be found at higher energies as exotic objects like low-mass X-ray binaries (LMXBs; e.g. Hut et al., 1991), millisecond pulsars (MSPs; e.g. Verbunt et al., 1987) and cataclysmic variables (CVs; e.g. Pooley & Hut, 2006; Cohn et al., 2010). Blue straggler (BS) formation is also thought to involve binary or triple stars (e.g. Perets & Fabrycky, 2009; Knigge et al., 2009; Leigh & Sills, 2011). This can occur in several different ways, such as mass transfer within a binary, collisions during encounters involving binaries or even some triple-based mechanism (e.g. Leigh & Sills, 2011; Geller et al., 2013) such as KL-induced mergers (e.g. Perets & Fabrycky, 2009; Naoz & Fabrycky, 2014).

To date, only one triple star system is known to exist in the dense environments characteristic of GCs (e.g. Prodan & Murray, 2012). The system in question, called 4U 1820-30, lives near the centre of the GC NGC 6624. It consists of a low-mass X-ray binary with a NS primary and a WD secondary, in orbit with a period $\sim$ $685$ s. There is also a large luminosity variation for this system (a factor $\sim$ 2) with a period of $\sim$ 171 days. This longer period is thought to be due to the presence of a tertiary companion (Grindlay et al., 1988).

2.4 Nuclear clusters

Much less is known of the multiple star populations in nuclear star clusters (NSCs). Binaries can only be detected at home, in the Milky Way’s Galactic nucleus. But here, significant reddening intervenes along the relevant line of sight, rendering observations of the Galactic Centre challenging at most wavelengths. If present, however, they should experience direct dynamical interactions with single stars in such dense environments on short timescales, and most likely be disrupted in the process (e.g. Leigh et al., 2016, 2018).

There are at least 3 known stellar binaries within the central $0.2$ pc of the Galactic Centre (e.g. Naoz et al., 2018; Ott et al., 1999; Rafelski et al., 2007). What’s more, it is still possible that many of the observed S-stars are unresolved binaries (e.g. Naoz et al., 2018). Even S0-2, which has the closest pericentre passage to Sgr A* and has received relatively careful scrutiny in recent years, could still be an unresolved binary (Chu et al., 2018).

More recently, Hailey et al. (2018) discovered $13$ X-ray sources within $\sim 1$ pc from Sgr A*, at least in projection. The authors argue that the accretors must be BHs instead of NSs for the majority of their sample. This argument is based on attempting to relate the properties of the observed spectra to those of other NSs known to be present in the Galactic Centre. These authors argue that the new X-ray sources show, in this regard, unique spectral features relative to the other known NSs in the Galactic Centre.

As for triples, to the best of our knowledge, there are none yet known to exist in our Galactic Centre. If they do exist, they likely formed from a relatively recent episode of star formation, which there is plenty of observational evidence to support (Bartko et al., 2009). This is because in the high-density, high-velocity dispersion environment of the Galactic Centre, triples should be rapidly destroyed due to dynamical interactions with single stars (e.g. Leigh et al., 2018).

3 Timescales

A triple system made up of an inner binary that is orbited by an outer companion undergoes KL oscillations in eccentricity whenever the initial mutual orbital inclination of the inner and outer orbit is in the window $i_{0}\sim 40^{\circ}$–$140^{\circ}$. The KL oscillations occur on a secular quadrupole timescale (Antognini, 2015; Naoz, 2016)

[TABLE]

where $m_{\rm 3b}$ is the mass of the outer body and $e_{\rm 3b}$ its eccentricity, $m_{\rm tot}$ is the total mass of the triple system, and $P_{{\rm bin}}\propto a_{{\rm bin}}^{3/2}$ and $P_{{\rm 3b}}\propto a_{{\rm 3b}}^{3/2}$ are the orbital periods of the inner and outer binary, respectively. At the quadruple order of approximation (inner test particle and outer circular orbit), the maximal eccentricity is a function of the initial mutual inclination

[TABLE]

We note that there is also a dependence on the relative orbital orientation (Kinoshita & Nakai, 1999). At this order of approximation, the inner binary eccentricity approaches unity when $i_{0}$ approaches $90^{\circ}$ (if not quenched by other sources of precession). In the case the outer orbit is eccentric, the inner eccentricity can reach almost unity even if the initial inclination is outside of the $i_{0}\sim 40^{\circ}$-$140^{\circ}$ KL range (octupole order of approximation; Naoz et al., 2013). This happens over the octupole timescale

[TABLE]

where the octupole parameter is defined as

[TABLE]

with $m_{\rm bin,1}$ and $m_{\rm bin,2}$ being the masses of the objects in the inner binary111According to Antognini (2015), the octupole timescale is rather $T_{\rm oct}=T_{\rm LK}/\epsilon^{1/2}$.. In the case of a CO binary, the large values reached by the eccentricity of the inner binary during the KL excursions make its merger time shorter since it efficiently dissipates energy when $e\sim e_{\rm max}$ (e.g. Blaes et al., 2002; Thompson, 2011; Antonini & Perets, 2012; Antognini, 2015). Tidal bulges and relativistic precession can suppress KL cycles (e.g. Naoz, 2016). In particular, for COs the relevant process is general relativistic precession that operates on a timescale

[TABLE]

Thus, the KL oscillations of the orbital elements are damped by relativistic effects in the region of parameter space where $T_{KL}>T_{GR}$. On the other hand, Naoz et al. (2013) showed that there can also be resonant behaviour when the KL and 1PN (Post-Newtonian) timescales are comparable, which could lead to eccentricity excitations.

At the quadrupole order, the Hamiltonian representing the evolution of our systems can be described as composed of three three-body Hamiltonians (Hamers et al., 2015): (i) SMBH and the innermost binary of the CO triple; (ii) SMBH and the binary made up of the centre-of-mass of the inner binary and the third companion of the CO triple; (iii) third companion and the inner binary of the CO triple. Therefore, for the systems studied in this paper, there are three relevant (quadrupole) KL timescales. The first can be calculated by considering the SMBH as a perturber of the inner binary of the CO triple. In this case $m_{\rm tot}\approx{M_{\rm SMBH}}$ and

[TABLE]

The second KL timescale is obtained when considering the effect of the SMBH on the binary made up of the centre-of-mass of the inner binary and the third companion of the CO triple. Also in this case, $m_{\rm tot}\approx{M_{\rm SMBH}}$ and

[TABLE]

where $m_{\mathrm{t}}=m_{11}+m_{12}+m_{3}$. Of course, $T_{\rm KL}^{123M}<T_{\rm KL}^{12M}$. The last KL timescale refers to tidal perturbations of the third BH/NS on the inner binary of the CO triple (”internal” KL timescale of the triple)

[TABLE]

Figure 2 shows these three relevant KL timescales as a function of distance from the SMBH, for different parameters of the triple. For reference, we fix $m_{11}=m_{12}=m_{3}=10\,\rm M_{\odot}$. The internal KL timescale for the CO triple is mainly set by the (outer) semi-major axis $a_{3}$. For the case ${a_{12}}=10$ AU and ${a_{3}}=50$ AU, $T_{\rm KL}^{123}$ is smaller than $T_{\rm KL}^{12M}$ and $T_{\rm KL}^{123M}$ when ${a_{\rm out}}\gtrsim 10^{3}$ AU for ${M_{\rm SMBH}}=4\times 10^{6}\,\rm M_{\odot}$ (i.e. a Milky Way-like galaxy). For the same ${a_{12}}$ and ${a_{3}}=1000$ AU, $T_{\rm KL}^{123}$ becomes the smallest only for ${a_{\rm out}}\gtrsim 10^{4}$ AU–$10^{5}$ AU. On the other hand, if we reduce the inner semi-major axis of the triple to ${a_{12}}=1$ AU, the typical distance at which the internal triple KL mechanism operates on the shortest timescale increases by a factor of $\sim 2$.

In Figure 2, we also plot two additional dynamical timescales relevant for triples in galactic nuclei. The first one is the typical encounter timescale for triples with other (single) stars (Leigh & Sills, 2011),

[TABLE]

where $f_{b}$ and $f_{t}$ are the fractions of binaries and triples in the galactic nucleus, $\rho$ is the density in stars and $\sigma$ is the velocity dispersion, and $m_{*}$ the mass of a typical object in the system. For timescales $\gtrsim\tau_{1+3}$, the triple survives around the SMBH, but with orbital elements different from the initial ones. The second timescale is the evaporation timescale $T_{\rm EV}$ due to the subsequent dynamical interactions with field stars in the dense environment of a galactic nucleus. A triple can evaporate when

[TABLE]

where $E_{t}$ is the internal orbital energy of the triple and $m_{\mathrm{t}}=m_{11}+m_{12}+m_{3}$. This process happens on an evaporation timescale222To compute the relevant quantities in Eq. 9 and Eq. 11, we have assumed the stars are distributed according to a Bahcall & Wolf (1976) cusp and normalization as given by Eq. 4 in Hoang et al. (2018). (Binney & Tremaine, 1987)

[TABLE]

where $\ln\Lambda$ is the Coulomb logarithm. Both the encounter timescale and the evaporation timescale are of the order of $\sim$ Myr for a Milky Way-like galaxy. For timescales $\sim T_{\rm EV}$, the triple is expected to be broken up by the many dynamical interactions with other stars.

Finally, there are two other relevant timescales, that we show in Figure 2, which are the timescales associated with resonant relaxation (RR). This effect is the result of torques originating from correlated encounters with stars orbiting the SMBH. For a given star or CO, RR affects both the eccentricity (scalar RR; SRR) and the inclination (vector RR; VRR) of the orbital plane with respect to the SMBH on timescales (e.g. Hamers et al., 2018)

[TABLE]

and

[TABLE]

respectively, where $N_{*}$ is the number of stars orbiting the SMBH within ${a_{\rm out}}$.

In the case of a CO binary orbiting an SMBH, we expect most of the mergers due to the KL mechanism to happen when $i_{0}\sim 90$ (Fragione et al., 2018). When a CO triple is present, the relative interplay among the three different KL mechanisms can bring the CO orbits into the relevant KL window, and hence induce mergers via GW emission because of a faster energy dissipation rate supplied by the increasing eccentricity. This is also favoured by the possible resonance between nodal precession and KL oscillations, whenever these two mechanisms operate on comparable timescales (Hamers & Lai, 2017; Grishin et al., 2018).

4 N-Body Simulations

In this section we quantify the physical properties of the merging triples with direct N-body simulations in order to compute their merger rates and assess whether they are enhanced with respect to the binary case.

4.1 Initial conditions

In our study, we consider an SMBH orbited by a triple CO, made up of either three BHs or three NSs. We assume an SMBH mass of $4\times 10^{6}\,\rm M_{\odot}$ for a Milky-Way like nucleus.

For the assumed distribution of BH masses, the initial mass function in triples is highly uncertain (and also likely to vary depending on the formation mechanism). Therefore we use a generic, parameterized form corresponding to a negative power-law distribution

[TABLE]

in the mass range $5\,\rm M_{\odot}$–$100\,\rm M_{\odot}$333Note that pulsational pair instabilities may limit the maximum mass to $\sim 50\,\rm M_{\odot}$ (Belczynski et al., 2016)., and sample the three BH masses independently. To study how the results depend on our assumptions, we vary the slope of the BH mass function in the range $\beta=1$, $2$, $3$, $4$ (O’Leary et al., 2016). For NSs, we fix the mass to $1.3\,\rm M_{\odot}$ (e.g. Lattimer & Prakash, 2005; Fragione, Pavlík & Banerjee, 2018).

The numbers and spatial profiles of BHs and NSs surrounding SMBHs are unconstrained from an observational point of view. Stars and COs tend to form a power-law density cusp ($n(r)\propto r^{-\alpha}$) around an SMBH, where lighter and heavier objects develop shallower and steeper cusps, respectively (Bahcall & Wolf, 1976). Typically, stars tend to have $\alpha\sim 1.5$–$1.75$, whereas BHs tend to have $\alpha\sim 2$–$3$ due to mass segregation (Alexander, 2017). We assume that the BH and NS number densities follow a cusp with $\alpha=2$, and we study the effects of the cusp slope, by considering a steeper cusp ($\alpha=3$) for BHs, and a shallower cusp ($\alpha=1.5$) for NSs. We take the maximum outer semi-major axis to be $0.1$ pc444This value corresponds to the distance at which the octupole timescale is equal to the timescale on which accumulated fly-bys from stars tend to unbind the binary (Hoang et al., 2018). (Hoang et al., 2018), and sample the outer orbital eccentricity from a thermal distribution (Jeans, 1919).

The triple (BH-BH-BH or NS-NS-NS) inner and outer semi-major axis and eccentricity distributions are not well known observationally. Apart from their initial distributions, it should be taken into account that the dense environments characteristic of galactic nuclei or star clusters could cause both distributions to change or diffuse over time (Hopman, 2009). We sample the triple inner and outer semi-major axis ${a_{12}}$ and ${a_{3}}$, respectively, by adopting a log-uniform distribution. Given the uncertainties in the orbital parameters of the triples, we explore a wide range of values for the outer triple semi-major axis; in particular, we study models with $a_{\rm 3,max}=50$ AU–$200$ AU–$1000$ AU. For the eccentricity, we consider models where we sample both $e_{12}$ and $e_{3}$ from uniform and thermal distributions. This is done to roughly reflect the different formation pathways of CO triples. Though the eccentricity mostly depends on the natal-kick and/or the common-envelope phase, it also depends on scattering and perturbations of the CO binaries by local COs and stars. As a consequence, we expect that in-situ formation should naively favour circular systems, while CO triples that migrate from farther out should most likely prefer a thermal distribution, as a consequence of energy and angular momentum exchange through many dynamical encounters.

The relative inclinations of the inner and outer orbits are also poorly constrained observationally. In the field, observations have shown that inner and outer orbital planes of triples composed of massive stars tend to be inclined relative to each other (Moe & Kratter, 2018), which may hint at moderately high relative inclinations between the inner and outer orbital planes of the triple COs. For the sake of generality, we draw the mutual inclination angles between the inner binary orbit with respect to the SMBH orbit ($i$) and the third companion in the triple’s orbit ($i_{3}$) from an isotropic distribution.

After we sample from the relevant distributions, we check that the Mardling & Aarseth (2001) criterion, namely

[TABLE]

is satisfied for each hierarchy in our four-body system in order to ensure orbital stability. In the previous equation, $m_{\rm in,tot}$ and $m_{\rm per}$ are the total inner mass and the perturber mass, respectively, $R_{\rm p,per}$ is the pericenter of the perturber, $a_{\rm in}$ is the inner semi-major axis, and $I$ (either $i$ or $i_{3}$ in our system) is the relative orbital inclination (in radians).

Given the above sets of initial parameters, we integrate the system of differential equations of motion of the four bodies

[TABLE]

with $i=1$,$2$,$3$,$4$. The integrations are performed using the archain code (Mikkola & Merritt, 2006, 2008), a fully regularized code able to model the evolution of systems of arbitrary mass ratios and large eccentricities with extreme accuracy. archain includes PN corrections up to order PN2.5. As a criterion for merger, we require that the relative distance between the COs is smaller than three times the sum of their relative radii.

For each of the $13$ different models (see Tab. 2), we run $\sim 800$–$1000$ realizations, for a total of $\sim 11000$ simulations. We fix the total integration time to $T=1$ Myr. Over this timescale, VRR may efficiently operate and reorient the binary centre-of-mass orbital plane with respect to the SMBH, thus affecting the relative inclination of the inner and outer orbital planes (we do not take into account VRR in our simulations). In this regard, VRR can affect the KL evolution (of the outer binary of the CO triple) even if the KL timescale is significantly shorter than $T_{\rm VRR}$. As shown in Hamers et al. (2018), the distribution of maximum eccentricities is already significantly different (higher eccentricities) if $T_{\rm LK}/T_{\rm VRR}\sim 0.01$. This also affects the relative LK dynamics and renders the $4$-body approximation insufficient. The in-plane precession induced by the nuclear cluster potential and departures from spherical symmetry in the galactic nucleus would make the CO center of mass orbit precess even faster than is quantified by vector resonant relaxation alone in a MW-like nucleus (Petrovich & Antonini, 2017). However, a relevant constraint for the triples lifetime in galactic nuclei comes from their evaporation as a result of the encounters with other stars, which takes place on timescales $\lesssim$ Myr for a Milky Way-like galaxy. Therefore, our final criterion for the total integration time is

[TABLE]

In our simulations, the CO triple can undergo several fates: (i) the inner binary of the CO triple merges producing a GW merger event; (ii) the CO triple survives with orbital elements perturbed with respect to the initial ones; (iii) the CO triple is tidally broken apart by differential forces exerted by the SMBH, and its components will either be captured by the SMBH or ejected from the galactic nucleus. In particular, the latter can be caused by KL cycles of the outer orbit of the BH/NS triple, which becomes highly eccentric and destabilizes eventually disrupting the triple. This type of evolution was recently found to be relevant in population synthesis studies of 3+1 quadruples (Hamers, 2018, 2019).

4.2 Inclination distribution

We show in Fig. 3 the distribution of the initial orbital inclinations $i$ of the inner binaries that merge in our simulations with respect to the SMBH (top panel), for $\beta=1$, $\alpha=2$ and different ${a_{\rm 3,max}}$. For the case of a binary CO orbiting an SMBH, most of the binaries that merge have initial inclinations $i\sim 90^{\circ}$ with respect to the SMBH, where the enhancement of the maximum eccentricity is expected to be the largest due to KL oscillations whenever not suppressed by GR precession (Fragione et al., 2018). The distribution thus results in a sharp peak at $\sim 90^{\circ}$ and only a few mergers have larger or smaller initial inclinations. For a triple CO orbiting an SMBH, this distribution is no longer sharply peaked, but instead is almost isotropic.

On the other hand, the distribution of the orbital inclination of the third companion with respect to the inner binary in the triple $i_{3}$ (bottom panel) is found to peak at $\sim 100^{\circ}$, but with non-negligible tails. For comparison, isolated triples that merge due to the KL mechanism show a very pronounced peak at $\sim 90^{\circ}$, with only a few mergers happening in low-inclination systems (Antonini et al., 2017). There are two possible caveats to this. First, possible resonances between nodal precession and KL oscillations can arise. As shown in Hamers & Lai (2017), this could make even low-inclination systems merge. Second, there are three different KL mechanisms competing, thus the eccentricity of the inner binary of the CO triple can be pumped up by the torque either of the SMBH or the third companion in the CO triple, for which the KL timescale is shorter (e.g. Hamers et al., 2015). We define the following ratios

[TABLE]

[TABLE]

[TABLE]

which we plot in Fig. 4 for the BH-BH binaries that merge in our simulations in Models MW, $\alpha=2$, ${a_{\rm 3,max}}=50$ AU and different values of $\beta$. Clearly, the shape of the three distributions does not depend on the assumed value of $\beta$, and we find that $\mathcal{R}_{12M}$, $\mathcal{R}_{123M}$ and $\mathcal{R}_{M}$ peak at $\sim 10^{-4}$, $\sim 5\times 10^{-3}$ and $\sim 10^{-1}$, respectively. We note that, however, KL cycles in a given CO triple may be inactive if the relative inclinations are not in the KL window. As a consequence, the initial dynamical evolution of the CO triple can ultimately be dictated by a KL cycle that takes place on a long timescale. The latter can in turn excite the relative inclination of one of the other orbits, which could activate KL cycles that take place on a shorter timescale. Thus, the CO triple may experience rather different and rich dynamical paths, which tend to lead to chaotic behaviour (Grishin, Lai & Perets, 2018).

4.3 Mass distribution

Figure 5 (top panel) illustrates the distribution of the total mass of merging BH-BH-BH triples orbiting a $4\times 10^{6}\,\rm M_{\odot}$ SMBH in a Milky Way-like nucleus, for different values of $\alpha$, $\beta$ and ${a_{\rm 3,max}}$. For binaries orbiting an SMBH (Fragione et al., 2018), the resulting total mass distribution is not significantly affected by the slope of the density profile $\alpha$ around the SMBH, with a roughly constant shape in the range $\sim 25\,\rm M_{\odot}$–$125\,\rm M_{\odot}$ and a tail extending up to $\sim 180\,\rm M_{\odot}$. Also, the extent of the triple does not influence the total mass distribution. The only parameter that sets the distribution properties is the slope $\beta$ of the BH mass function. As expected, we find that the shallower the slope of the BH mass function, the larger the typical total mass of merging BH-BH binaries. For $\beta=1$, the distribution is peaked at a total mass of $\sim 30\,\rm M_{\odot}$ with a tail that extends up to $\sim 125\,\rm M_{\odot}$. For larger values of $\beta$, the tail gradually disappears and the distribution is peaked at $\sim 15\,\rm M_{\odot}$. We find that $\sim 90\%$ of the mergers have total masses $\lesssim 20\,\rm M_{\odot}$, $\lesssim 30\,\rm M_{\odot}$, $\lesssim 70\,\rm M_{\odot}$, and $\lesssim 100\,\rm M_{\odot}$ for $\beta=1$, $2$, $3$, and $4$, respectively.

In the bottom panel of Fig. 5, we show a comparison between the total mass of merging BHs in binaries (Fragione et al., 2018) and triples orbiting a $4\times 10^{6}\,\rm M_{\odot}$ SMBH, assuming $\alpha=2$ and different values for the parameter $\beta$ (for triples we fix ${a_{\rm 3,max}}=50$ AU). We find that the total mass of merging binaries is nearly independent of the multiplicity of the system. Thus, the main parameter that governs the shape of the final distribution for the total mass of merged BHs is the slope of the BH mass function.

4.4 Merger times

Due to KL cycles, the eccentricity of the inner binary in the triple can approach almost unity when $i_{0}\sim 90^{\circ}$ and the inner binary merger time becomes shorter due to efficient energy dissipation at pericentre via KL oscillations. In the top panel of Fig. 6, we show the cumulative distribution functions (CDFs) of merger times for BH-BH binaries in triples orbiting an SMBH of mass ${M_{\rm SMBH}}=4\times 10^{6}\,\rm M_{\odot}$. Different values for the parameters $\alpha$, $\beta$ and ${a_{\rm 3,max}}$ do not affect the distribution of merger times.

In the bottom panel of Fig. 6, we compare the CDFs of BHs merging in binaries (Fragione et al., 2018) to those in triples orbiting a Milky Way-like SMBH, for different values of $\beta$. BHs merge in triples on shorter timescales compared to BHs that merge in binaries, independent of the slope of the BH mass function (${a_{\rm 3,max}}=50$ AU, $\alpha=2$). We find that $\sim 50\%$ of the mergers have typical merger timescales of $\lesssim 3\times 10^{3}$ yr and $\lesssim 3\times 10^{4}$ yr, respectively, for the binary and triple cases. The most likely reason for these deviations comes from a richer dynamical evolution, with three KL cycles competing, combined with possible resonances that can arise in the case of triples orbiting an SMBH.

4.5 Eccentricity in the LIGO band

Hierarchical configurations are expected to have eccentricities when entering the LIGO band ($10$ Hz) that are larger than for binaries that merge in isolation (see e.g. Fragione & Loeb, 2019). For the CO binaries that merge in our simulations, we compute a proxy for the GW frequency, i.e. the frequency corresponding to the harmonic that gives the maximum GW emission (Wen, 2003)

[TABLE]

In Figure 7, we illustrate the distribution of eccentricities at the moment the CO binaries enter the LIGO frequency band. We also plot the minimum $e_{\rm 10Hz}=0.081$ where the LIGO/VIRGO/KAGRA network may distinguish eccentric sources from circular sources (Gondán & Kocsis, 2019). A large fraction of systems that merge have a significant eccentricity in the LIGO band, compared to binaries that merge in isolation. We note that a similar signature could be found in CO binaries that merge near intermediate-mass black holes (Fragione & Bromberg, 2019), in the GW capture scenario in star clusters (Samsing, 2018), and in isolated hierarchical triples (Antonini et al., 2017) and quadruples (Fragione & Kocsis, 2019).

5 Discussion and conclusions

We have presented a new scenario for merging BHs and NSs in galactic nuclei. We have studied the mergers of binary BHs and NSs in triples orbiting an SMBH in a Milky Way-like galactic nucleus. By conducting a series of high-precision numerical integrations, we have quantified the relative differences with the case where the BHs and NSs are in binaries (Fragione et al., 2018). We find that, unlike the binary case in which the initial orbital inclination of the merging BHs/NSs is sharply peaked at $\sim 90^{\circ}$, in the case of triples the distribution is significantly broader, due to the more complex dynamics of the triple-SMBH system, with three KL cycles competing, and possible resonances between nodal precession and KL oscillations.

To accurately determine the global BH merger rate from this channel, we would need to quantify the population of triple BHs and NSs in galactic nuclei and star clusters, which is highly uncertain. For the binary case, this rate is thought to be of the order of $\sim 0.1$–$1\rm{\ Gpc}^{-3}\rm{\ yr}^{-1}$ (Antonini & Perets, 2012; Fragione et al., 2018; Hamers et al., 2018; Hoang et al., 2018). Nevertheless, the multiplicity fraction of high-mass main-sequence stars is large, with each star having two or more companions on average (Sana, 2017). We find that the merger fraction in triple BH and NS systems can be larger than the one for binaries, by up to a factor of $\sim 8$ (see the last column of Tab. 2). Our results thus suggest that dynamically-driven BH and NS mergers in this scenario may be important and contribute to the merger events observed by LIGO/VIRGO. Moreover, as for other hierarchical systems, we have found that binaries merging in CO triples are expected to have eccentricities in the LIGO band ($10$ Hz) larger than binaries that merge in isolation, with a peak at $\sim 1$.

We note that in the case the inner binary merges, we stop our simulations. In reality, however, the merger remnant should be imparted a recoil velocity kick that could unbind the system (e.g. Lousto et al., 2010), now a binary comprised of the merger remnant itself and the third companion in the former triple. If the kick does not unbind the system, a second-generation binary is born that can undergo KL oscillations as a result of the SMBH tidal field, which can eventually result in a second merger. Therefore, the presence of the massive SMBH could further enhance the merger rate from this scenario, by a factor of up to $\sim 2$.

We caution, however, that there are a lot of uncertainties regarding the exact population of stars and CO triples, as well as their mass and orbital parameters distributions. Moreover, we have not taken into account the pre-CO evolution, which could be relevant in determining the properties of triple COs (Hamers et al., 2013; Shappee & Thompson, 2013; Antonini et al., 2017; Toonen et al., 2018). In our simulations, we only consider triples made up of all BHs or all NSs, which implies that their progenitors had a mass ratio near unity. While this is the case for a subset of systems, additional simulations incorporating all possible compositions of the CO triples will be needed to precisely compute how the merger fraction depends on the specific composition of the triple. Importantly, the maximum outer semi-major axis of triples plays an important role in determining the relevant timescales. In particular, both the the 3+1 interaction timescale and the triple evaporation timescale depend crucially on it. As we show, these timescales can be as short as a few Myrs, such that that more compact triples will probably dominate the rate. This would in turn imply that CO triples in galactic nuclei need to be born more compact than in the field, if they were to contribute to the observed GW merger events.

What makes these events especially interesting for the GW community is the fact that they have enhanced chance of being detected also in the electromagnetic (EM) spectrum. For a binary NS merger, the intensity of the gamma-ray radiation, which primarily depends on the amount of mass available for accretion, is primarily determined by the properties (namely mass and equation of state) of the merging stars (Shibata et al., 2006; Rezzolla et al., 2010; Giacomazzo et al., 2013; Hotokezaka et al., 2013). However, if the merger event drives a relativistic jet as evinced by the observations of GW170817/GRB170817A (Abbott et al., 2017) and detailed broadband modeling (Lazzati et al., 2018), then the density $n$ of the medium in which the event occurs also plays an important role. In fact, the relativistic jet will produce bright radiation (known as an afterglow), spanning the full EM spectrum, from hard X-rays to the radio, and the peak of this emission scales as $\propto n^{1/2}$ (Sari et al., 1998). The binary NS mergers from triples orbiting an SMBH studied here typically have short merger times, even shorter than those of isolated binaries alone, as studied by Fragione et al. (2018). As such, they are expected to merge in the innermost regions of a galaxy, where densities are at their highest. For gas densities of a few cm*-3*, the EM brightness would allow for their detection in multiple bands well beyond the LIGO horizon, as computed by Fragione et al. (2018). EM observations are a crucial element allowing for event localization, and hence to test the nuclear origins of these events.

Mergers of binaries made up of two BHs are more likely to be dark, lacking the natural reservoir of accretion matter provided by the tidal disruption of a NS immediately preceding the final coalescence. While some reservoir of material could still be found around one of the BHs as a dead disk remnant from the SN explosion (Perna et al., 2016), this may not be a common scenario. If the SMBH is surrounded by a dense disk as in active galactic nuclei then, if the binary BH merges within the disk, a temporary accretion of matter at high rates can ensue (Bartos et al., 2017), resulting in a possible EM counterpart. Nevertheless, for the case of a binary BH in a triple there is the additional possibility of producing electromagnetic radiation if the third object, instead of being a BH, is a main sequence star that fills its Roche lobe and generates a circumbinary disk around the inner binary BH (Chang & Murray, 2018). As the binary BH merges, the loss of energy via GWs along with a possible merger kick can result in orbital crossings of the disk particles, resulting in disk shocking. Alternatively (or additionally) the circumbinary disk will eventually accrete on a viscous timescale upon merger of the binary. All together, binary mergers from triples near SMBHs may constitute an especially important channel for combined GW/EM transients.

Acknowledgements

We thank the anonymous referee for useful comments. GF is supported by the Foreign Postdoctoral Fellowship Program of the Israel Academy of Sciences and Humanities. GF also acknowledges support from an Arskin postdoctoral fellowship. NL and RP acknowledge support by NSF award AST-1616157. GF thanks Seppo Mikkola for helpful discussions on the use of the code archain. Simulations were run on the Astric cluster at the Hebrew University of Jerusalem. The Center for Computational Astro-physics at the Flatiron Institute is supported by the Simons Foundation.