Double neutron star mergers from hierarchical triple-star systems

TL;DR

This study explores how hierarchical triple-star systems can significantly enhance neutron star merger rates through Lidov-Kozai oscillations, especially if neutron stars are born with low natal kicks, addressing gaps in previous binary evolution models.

Contribution

It introduces detailed simulations of massive triple-star systems incorporating stellar evolution, secular dynamics, and supernova effects to evaluate their impact on neutron star merger rates.

Findings

Triple systems can increase NS-NS merger rates if natal kicks are low.

High natal kicks reduce the influence of triples on merger rates.

Triple-induced mergers could explain observed LIGO event rates if low kicks are common.

Abstract

The isolated binary evolution model for merging neutron stars (NSs) involves processes such as mass transfer, common-envelope evolution, and natal kicks, all of which are poorly understood. Also, the predicted NS-NS merger rates are typically lower than the rates inferred from the LIGO GW170817 event. Here, we investigate merger rates of NS and black hole (BH)-NS binaries in hierarchical triple-star systems. In such systems, the tertiary can induce Lidov-Kozai (LK) oscillations in the inner binary, accelerating its coalescence, and potentially enhancing compact object merger rates. However, since compact objects originate from massive stars, the prior evolution should also be taken into account. Natal kicks, in particular, could significantly reduce the rates by unbinding the tertiary before it can affect the inner binary through LK evolution. We carry out simulations of massive triples taking into account stellar evolution starting from the main sequence, secular and tidal evolution, and the effects of supernovae. For large NS birth kicks (\sigma_k = 265 km/s), we find that the triple NS-NS merger rate (several hundred Gpc^{-3} yr^{-1}) is lower by a factor of ~2-3 than the binary rate, but for no kicks (\sigma_k = 0 km/s), the triple rate (several thousand Gpc^{-3} yr^{-1}) is comparable to the binary rate. Our results indicate that a significant fraction of NS-NS mergers could originate from triples if a substantial portion of the NS population is born with low kick velocities, as indicated by other work. However, uncertainties and open questions remain because of our simplifying assumption of dynamical decoupling after inner binary interaction has been triggered.

Figures

Tables

Table 1. Table 1 Description of important quantities and their initial value(s) and/or distributions assumed in the population synthesis.

Symbol	Description	Initial value(s) and/or distribution in population synthesis
$m_1$	Mass of the primary star.	$8-50M_{\odot }$ with a Salpeter initial mass function (Salpeter 1955, i.e., $dN/d{m}_1\propto m_1^{-2.35}$).
$m_2$	Mass of the secondary star.	$m_1{q}_1$, where $q_1$ has a flat distribution, and with $m_2>4M_{\odot }$.
$m_3$	Mass of the tertiary star.	Either (1) $m_3=q_2\left(m_1+m_2\right)$, or (2) $m_3=q_2^{\prime }{m}_2$, where both $q_2$ and $q_2^{\prime }$ have a flat distribution, subject to $m_3>0.1M_{\odot }$.
$Z_i$	Metallicity of star $i$.	$0.02$
$R_i$	Radius of star $i$.	From stellar evolution code.
$P_{s,i}$	Spin period of star $i$.	$10d$
${\theta }_{s,i}$	Obliquity (spin-orbit angle) of star $i$.	$0^{\circ }$
$t_{V,i}$	Viscous time-scale of star $i$.	Computed from the stellar properties using the prescription of Hurley et al. (2002).
$k_{\mathrm{AM},i}$	Apsidal motion constant of star $i$.	0.014
$r_{g,i}$	Gyration radius of star $i$.	0.08
$V_k$	SNe kick speed.	Maxwellian distribution with ${\sigma }_k=0$, 40, or $265\mathrm{km}{s}^{-1}$.
$P_{\mathrm{orb},i}$	Orbital period of orbit $i$ (inner orbit: $i=1$; outer orbit: $i=2$).	Flat distribution in ${\log }_{10}\left(P_{\mathrm{orb},i}\right)$, with , and subject to dynamical stability constraints.
$a_i$	Semimajor axis of orbit $i$.	Computed from $P_{\mathrm{orb},i}$ and the $m_i$ using Kepler’s law.
$e_i$	Eccentricity of orbit $i$.	Flat distribution between 0 and 0.9.
$i_i$	Inclination of orbit $i$.	$0-{180}^{\circ }$ (flat distribution in $\cos i_i$)
${\omega }_i$	Argument of periapsis of orbit $i$.	$0-{360}^{\circ }$ (flat distribution in ${\omega }_i$)
${\Omega }_i$	Longitude of the ascending node of orbit $i$.	$0-{360}^{\circ }$ (flat distribution in ${\Omega }_i$)

Table 2. Table 2 Outcome fractions of the simulations with the tertiary mass assuming a flat distribution in q 2 ≡ m 3 / ( m 1 + m 2 ) subscript 𝑞 2 subscript 𝑚 3 subscript 𝑚 1 subscript 𝑚 2 q_{2}\equiv m_{3}/(m_{1}+m_{2}) . The first three data columns correspond to triple systems, whereas the last three data columns correspond to the binary case, i.e., taking the same triple systems but in the absence of the tertiary star. For each case, we show results for the three different values of the kick velocity dispersion, σ k subscript 𝜎 k \sigma_{\mathrm{k}} . Fractions are based on N MC = 10 3 subscript 𝑁 MC superscript 10 3 N_{\mathrm{MC}}=10^{3} simulations, and quoted errors are based on Poisson statistics. Outcomes that do not apply (e.g., dynamically unstable systems in the binary case) are marked with ‘—’. Some of the stellar evolutionary states are G – giant star (including red giants and asymptotic giants), and core helium burning star (CHeB). Refer to the text in Section  4.1.1 for a description of the different channels.

	Fraction of all systems
	Triple			Binary
	${\sigma }_k/\mathrm{km}{s}^{-1}$			${\sigma }_k/\mathrm{km}{s}^{-1}$
	0	40	265	0	40	265
No Interaction	$0.024\pm 0.005$	$0.004\pm 0.002$	$0.004\pm 0.002$	$0.253\pm 0.016$	$0.011\pm 0.003$	$0.008\pm 0.003$
RLOF $\star$1	$0.607\pm 0.025$	$0.445\pm 0.021$	$0.442\pm 0.021$	$0.612\pm 0.025$	$0.612\pm 0.025$	$0.612\pm 0.025$
MS	$0.068\pm 0.008$	$0.063\pm 0.008$	$0.061\pm 0.008$	$0.056\pm 0.007$	$0.056\pm 0.007$	$0.056\pm 0.007$
G	$0.444\pm 0.021$	$0.314\pm 0.018$	$0.313\pm 0.018$	$0.460\pm 0.021$	$0.460\pm 0.021$	$0.460\pm 0.021$
CHeB	$0.095\pm 0.010$	$0.068\pm 0.008$	$0.068\pm 0.008$	$0.096\pm 0.010$	$0.096\pm 0.010$	$0.096\pm 0.010$
RLOF $\star$2	$0.024\pm 0.005$	$0.021\pm 0.005$	$0.021\pm 0.005$	$0.021\pm 0.005$	$0.009\pm 0.003$	$0.002\pm 0.001$
MS	$0.021\pm 0.005$	$0.021\pm 0.005$	$0.021\pm 0.005$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
G	$0.002\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.021\pm 0.005$	$0.009\pm 0.003$	$0.002\pm 0.001$
CHeB	$0.001\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
RLOF $\star$3	$0.009\pm 0.003$	$0.009\pm 0.003$	$0.008\pm 0.003$	—	—	—
MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
G	$0.004\pm 0.002$	$0.004\pm 0.002$	$0.004\pm 0.002$	—	—	—
CHeB	$0.005\pm 0.002$	$0.005\pm 0.002$	$0.004\pm 0.002$	—	—	—
Dynamical inst.	$0.028\pm 0.005$	$0.003\pm 0.002$	$0.004\pm 0.002$	—	—	—
MS+MS	$0.002\pm 0.001$	$0.002\pm 0.001$	$0.003\pm 0.002$	—	—	—
G+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
CHeB+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
NS+MS	$0.012\pm 0.003$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
NS+G	$0.003\pm 0.002$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
NS+CHeB	$0.003\pm 0.002$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
NS+NS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
BH+MS	$0.004\pm 0.002$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
BH+G	$0.001\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
BH+CHeB	$0.001\pm 0.001$	$0.001\pm 0.001$	$0.001\pm 0.001$	—	—	—
BH+NS	$0.002\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
BH+BH	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
Semisecular	$0.004\pm 0.002$	$0.004\pm 0.002$	$0.004\pm 0.002$	—	—	—
MS+MS	$0.004\pm 0.002$	$0.004\pm 0.002$	$0.004\pm 0.002$	—	—	—
Secular collision	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
MS+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
Unbound (SNe)	$0.304\pm 0.017$	$0.514\pm 0.023$	$0.517\pm 0.023$	$0.114\pm 0.011$	$0.368\pm 0.019$	$0.378\pm 0.019$
MS+MS	$0.000\pm 0.000$	$0.262\pm 0.016$	$0.265\pm 0.016$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
G+MS	$0.000\pm 0.000$	$0.004\pm 0.002$	$0.004\pm 0.002$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
CHeB+MS	$0.000\pm 0.000$	$0.017\pm 0.004$	$0.017\pm 0.004$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
NS+MS	$0.219\pm 0.015$	$0.136\pm 0.012$	$0.136\pm 0.012$	$0.024\pm 0.005$	$0.263\pm 0.016$	$0.269\pm 0.016$
NS+G	$0.000\pm 0.000$	$0.001\pm 0.001$	$0.001\pm 0.001$	$0.000\pm 0.000$	$0.001\pm 0.001$	$0.001\pm 0.001$
NS+CHeB	$0.022\pm 0.005$	$0.020\pm 0.004$	$0.020\pm 0.004$	$0.000\pm 0.000$	$0.029\pm 0.005$	$0.029\pm 0.005$
NS+NS	$0.007\pm 0.003$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.084\pm 0.009$	$0.001\pm 0.001$	$0.000\pm 0.000$
BH+MS	$0.019\pm 0.004$	$0.056\pm 0.007$	$0.056\pm 0.007$	$0.000\pm 0.000$	$0.046\pm 0.007$	$0.054\pm 0.007$
BH+G	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+CHeB	$0.001\pm 0.001$	$0.005\pm 0.002$	$0.005\pm 0.002$	$0.000\pm 0.000$	$0.004\pm 0.002$	$0.004\pm 0.002$
BH+NS	$0.028\pm 0.005$	$0.005\pm 0.002$	$0.005\pm 0.002$	$0.006\pm 0.002$	$0.014\pm 0.004$	$0.010\pm 0.003$
BH+BH	$0.008\pm 0.003$	$0.008\pm 0.003$	$0.008\pm 0.003$	$0.000\pm 0.000$	$0.010\pm 0.003$	$0.011\pm 0.003$
Inner bound	$0.286\pm 0.017$	$0.304\pm 0.017$	$0.295\pm 0.017$	—	—	—

Table 3. Table 3 Similar to Table  2 , here for simulations with the tertiary mass sampled assuming a flat distribution in q 2 ′ ≡ m 3 / m 2 superscript subscript 𝑞 2 ′ subscript 𝑚 3 subscript 𝑚 2 q_{2}^{\prime}\equiv m_{3}/m_{2} .

	Fraction of all systems
	Triple			Binary
	${\sigma }_k/\mathrm{km}{s}^{-1}$			${\sigma }_k/\mathrm{km}{s}^{-1}$
	0	40	265	0	40	265
No Interaction	$0.026\pm 0.005$	$0.006\pm 0.002$	$0.006\pm 0.002$	$0.221\pm 0.015$	$0.011\pm 0.003$	$0.008\pm 0.003$
RLOF $\star$1	$0.618\pm 0.025$	$0.622\pm 0.025$	$0.621\pm 0.025$	$0.595\pm 0.024$	$0.612\pm 0.025$	$0.612\pm 0.025$
MS	$0.060\pm 0.008$	$0.060\pm 0.008$	$0.061\pm 0.008$	$0.056\pm 0.007$	$0.056\pm 0.007$	$0.056\pm 0.007$
G	$0.457\pm 0.021$	$0.462\pm 0.021$	$0.460\pm 0.021$	$0.448\pm 0.021$	$0.460\pm 0.021$	$0.460\pm 0.021$
CHeB	$0.101\pm 0.010$	$0.100\pm 0.010$	$0.100\pm 0.010$	$0.091\pm 0.010$	$0.096\pm 0.010$	$0.096\pm 0.010$
RLOF $\star$2	$0.017\pm 0.004$	$0.016\pm 0.004$	$0.016\pm 0.004$	$0.018\pm 0.004$	$0.004\pm 0.002$	$0.002\pm 0.001$
MS	$0.014\pm 0.004$	$0.014\pm 0.004$	$0.014\pm 0.004$	$0.000\pm 0.000$	$0.001\pm 0.001$	$0.000\pm 0.000$
G	$0.002\pm 0.001$	$0.002\pm 0.001$	$0.002\pm 0.001$	$0.018\pm 0.004$	$0.003\pm 0.002$	$0.002\pm 0.001$
CHeB	$0.001\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
RLOF $\star$3	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
G	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
CHeB	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
Dynamical inst.	$0.035\pm 0.006$	$0.006\pm 0.002$	$0.005\pm 0.002$	—	—	—
MS+MS	$0.004\pm 0.002$	$0.004\pm 0.002$	$0.004\pm 0.002$	—	—	—
G+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
CHeB+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
NS+MS	$0.013\pm 0.004$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
NS+G	$0.005\pm 0.002$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
NS+CHeB	$0.002\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
NS+NS	$0.001\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
BH+MS	$0.005\pm 0.002$	$0.001\pm 0.001$	$0.000\pm 0.000$	—	—	—
BH+G	$0.001\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
BH+CHeB	$0.002\pm 0.001$	$0.001\pm 0.001$	$0.001\pm 0.001$	—	—	—
BH+NS	$0.002\pm 0.001$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
BH+BH	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	—	—	—
Semisecular	$0.003\pm 0.002$	$0.002\pm 0.001$	$0.002\pm 0.001$	—	—	—
MS+MS	$0.002\pm 0.001$	$0.002\pm 0.001$	$0.002\pm 0.001$	—	—	—
Secular collision	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
MS+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
Unbound (SNe)	$0.301\pm 0.017$	$0.348\pm 0.019$	$0.350\pm 0.019$	$0.160\pm 0.013$	$0.373\pm 0.019$	$0.378\pm 0.019$
MS+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
G+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
CHeB+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
NS+MS	$0.218\pm 0.015$	$0.245\pm 0.016$	$0.245\pm 0.016$	$0.048\pm 0.007$	$0.266\pm 0.016$	$0.268\pm 0.016$
NS+G	$0.000\pm 0.000$	$0.001\pm 0.001$	$0.001\pm 0.001$	$0.000\pm 0.000$	$0.001\pm 0.001$	$0.001\pm 0.001$
NS+CHeB	$0.024\pm 0.005$	$0.028\pm 0.005$	$0.028\pm 0.005$	$0.000\pm 0.000$	$0.028\pm 0.005$	$0.029\pm 0.005$
NS+NS	$0.009\pm 0.003$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.099\pm 0.010$	$0.001\pm 0.001$	$0.000\pm 0.000$
BH+MS	$0.020\pm 0.004$	$0.050\pm 0.007$	$0.053\pm 0.007$	$0.003\pm 0.002$	$0.048\pm 0.007$	$0.055\pm 0.007$
BH+G	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+CHeB	$0.000\pm 0.000$	$0.004\pm 0.002$	$0.004\pm 0.002$	$0.000\pm 0.000$	$0.004\pm 0.002$	$0.004\pm 0.002$
BH+NS	$0.026\pm 0.005$	$0.008\pm 0.003$	$0.008\pm 0.003$	$0.009\pm 0.003$	$0.012\pm 0.003$	$0.010\pm 0.003$
BH+BH	$0.004\pm 0.002$	$0.012\pm 0.003$	$0.011\pm 0.003$	$0.001\pm 0.001$	$0.013\pm 0.004$	$0.011\pm 0.003$
Inner bound	$0.269\pm 0.016$	$0.012\pm 0.001$	$0.003\pm 0.002$	—	—	—

Table 4. Table 4 Fractions for merger outcomes in the ‘isolated binary’ simulations, which are based on the systems in which RLOF was triggered in the inner binary system. Data in this table are based on the high-mass tertiary case, i.e,. m 3 = q 2 ​ ( m 1 + m 2 ) subscript 𝑚 3 subscript 𝑞 2 subscript 𝑚 1 subscript 𝑚 2 m_{3}=q_{2}(m_{1}+m_{2}) ; see Table  5 for the low-mass tertiary case. The fractions in this table are given with respect to the systems undergoing RLOF triggered by the primary or secondary star; the fractions of the latter cases corresponding to all systems were given in Table  2 . Here, ‘He’ refers to a helium-burning star.

	Fraction of RLOF systems ($\star 1$ or $\star$ 2)
	Triple			Binary
	${\sigma }_k/\mathrm{km}{s}^{-1}$			${\sigma }_k/\mathrm{km}{s}^{-1}$
	0	40	265	0	40	265
RLOF $\to$ Merger	$0.743\pm 0.034$	$0.835\pm 0.042$	$0.905\pm 0.044$	$0.708\pm 0.033$	$0.842\pm 0.037$	$0.914\pm 0.039$
MS+MS	$0.079\pm 0.011$	$0.101\pm 0.015$	$0.102\pm 0.015$	$0.039\pm 0.008$	$0.040\pm 0.008$	$0.041\pm 0.008$
CHeB+MS	$0.035\pm 0.007$	$0.030\pm 0.008$	$0.032\pm 0.008$	$0.021\pm 0.006$	$0.021\pm 0.006$	$0.020\pm 0.006$
CHeB+G	$0.005\pm 0.003$	$0.004\pm 0.003$	$0.004\pm 0.003$	$0.003\pm 0.002$	$0.003\pm 0.002$	$0.003\pm 0.002$
G+MS	$0.052\pm 0.009$	$0.052\pm 0.011$	$0.052\pm 0.011$	$0.063\pm 0.010$	$0.064\pm 0.010$	$0.065\pm 0.010$
G+CHeB	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
G+G	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
He+MS	$0.016\pm 0.005$	$0.002\pm 0.002$	$0.004\pm 0.003$	$0.011\pm 0.004$	$0.011\pm 0.004$	$0.016\pm 0.005$
He+G	$0.005\pm 0.003$	$0.004\pm 0.003$	$0.004\pm 0.003$	$0.003\pm 0.002$	$0.003\pm 0.002$	$0.003\pm 0.002$
WD+MS	$0.073\pm 0.011$	$0.049\pm 0.010$	$0.050\pm 0.010$	$0.070\pm 0.010$	$0.072\pm 0.011$	$0.073\pm 0.011$
WD+G	$0.013\pm 0.004$	$0.017\pm 0.006$	$0.017\pm 0.006$	$0.022\pm 0.006$	$0.021\pm 0.006$	$0.021\pm 0.006$
WD+He	$0.002\pm 0.002$	$0.000\pm 0.000$	$0.002\pm 0.002$	$0.006\pm 0.003$	$0.006\pm 0.003$	$0.008\pm 0.004$
WD+WD	$0.002\pm 0.002$	$0.004\pm 0.003$	$0.004\pm 0.003$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
NS+MS	$0.209\pm 0.018$	$0.423\pm 0.030$	$0.576\pm 0.035$	$0.207\pm 0.018$	$0.417\pm 0.026$	$0.578\pm 0.031$
NS+CHeB	$0.011\pm 0.004$	$0.006\pm 0.004$	$0.004\pm 0.003$	$0.032\pm 0.007$	$0.011\pm 0.004$	$0.008\pm 0.004$
NS+G	$0.082\pm 0.011$	$0.047\pm 0.010$	$0.000\pm 0.000$	$0.073\pm 0.011$	$0.047\pm 0.009$	$0.008\pm 0.004$
NS+He	$0.016\pm 0.005$	$0.006\pm 0.004$	$0.004\pm 0.003$	$0.017\pm 0.005$	$0.014\pm 0.005$	$0.005\pm 0.003$
NS+WD	$0.013\pm 0.004$	$0.009\pm 0.004$	$0.002\pm 0.002$	$0.008\pm 0.004$	$0.005\pm 0.003$	$0.003\pm 0.002$
NS+NS	$0.105\pm 0.013$	$0.045\pm 0.010$	$0.009\pm 0.004$	$0.104\pm 0.013$	$0.058\pm 0.010$	$0.015\pm 0.005$
BH+MS	$0.008\pm 0.004$	$0.011\pm 0.005$	$0.013\pm 0.005$	$0.011\pm 0.004$	$0.011\pm 0.004$	$0.011\pm 0.004$
BH+CHeB	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+G	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.002\pm 0.002$	$0.002\pm 0.002$	$0.002\pm 0.002$
BH+He	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+WD	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+NS	$0.010\pm 0.004$	$0.021\pm 0.007$	$0.022\pm 0.007$	$0.005\pm 0.003$	$0.023\pm 0.006$	$0.021\pm 0.006$
BH+BH	$0.002\pm 0.002$	$0.002\pm 0.002$	$0.002\pm 0.002$	$0.002\pm 0.002$	$0.002\pm 0.002$	$0.002\pm 0.002$
Other	$0.008\pm 0.004$	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.009\pm 0.004$	$0.010\pm 0.004$	$0.010\pm 0.004$

Table 5. Table 6 Fractions for merger outcomes in the ‘isolated binary’ simulations, which are based on the systems in which the tertiary becomes unbound from the inner binary due to a SNe event. Data in this table are based on the high-mass tertiary case, i.e,. m 3 = q 2 ​ ( m 1 + m 2 ) subscript 𝑚 3 subscript 𝑞 2 subscript 𝑚 1 subscript 𝑚 2 m_{3}=q_{2}(m_{1}+m_{2}) ; see Table  7 for the low-mass tertiary case. The fractions in this table are given with respect to the systems in which the tertiary becomes unbound, but with the inner binary still bound; the latter fraction is given in the last row of Table  2 .

	Fraction of unbound triples (stable inner binary)
	Triple
	${\sigma }_k/\mathrm{km}{s}^{-1}$
	0	40	265
Tertiary unbound
$\to$ inner binary merger	$0.476\pm 0.041$	$0.931\pm 0.055$	$0.969\pm 0.057$
MS+MS	$0.000\pm 0.000$	$0.076\pm 0.016$	$0.078\pm 0.016$
CHeB+MS	$0.000\pm 0.000$	$0.003\pm 0.003$	$0.003\pm 0.003$
CHeB+G	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
G+MS	$0.000\pm 0.000$	$0.023\pm 0.009$	$0.024\pm 0.009$
G+CHeB	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
G+G	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
He+MS	$0.000\pm 0.000$	$0.043\pm 0.012$	$0.041\pm 0.012$
He+G	$0.000\pm 0.000$	$0.003\pm 0.003$	$0.003\pm 0.003$
WD+MS	$0.000\pm 0.000$	$0.109\pm 0.019$	$0.112\pm 0.019$
WD+G	$0.000\pm 0.000$	$0.016\pm 0.007$	$0.017\pm 0.008$
WD+He	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
WD+WD	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
NS+MS	$0.164\pm 0.024$	$0.513\pm 0.041$	$0.600\pm 0.045$
NS+CHeB	$0.003\pm 0.003$	$0.036\pm 0.011$	$0.037\pm 0.011$
NS+G	$0.000\pm 0.000$	$0.026\pm 0.009$	$0.003\pm 0.003$
NS+He	$0.000\pm 0.000$	$0.007\pm 0.005$	$0.003\pm 0.003$
NS+WD	$0.000\pm 0.000$	$0.003\pm 0.003$	$0.000\pm 0.000$
NS+NS	$0.252\pm 0.030$	$0.013\pm 0.007$	$0.007\pm 0.005$
BH+MS	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+CHeB	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+G	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+He	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+WD	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
BH+NS	$0.056\pm 0.014$	$0.053\pm 0.013$	$0.034\pm 0.011$
BH+BH	$0.000\pm 0.000$	$0.000\pm 0.000$	$0.000\pm 0.000$
Other	$0.000\pm 0.000$	$0.007\pm 0.005$	$0.007\pm 0.005$

Table 6. Table 8 Compact object merger rates according to our simulations. Top part: NS-NS mergers; middle part: BH-NS mergers. We also include the rate of NS-MS mergers (Thorne- Z ˙ ˙ Z \mathrm{\dot{Z}} ytkow objects) in the bottom part. We include mergers from triples and binaries, and for different kick dispersions σ k subscript 𝜎 k \sigma_{\mathrm{k}} . The assumed channels are mergers following RLOF, and after the tertiary becomes unbound from the inner binary (only applies to triples). The ‘Total’ row gives the sum of the two channels (note: quoted numbers have been rounded to integers). Results are shown for sets of simulations with the tertiary mass sampled according to m 3 = q 2 ​ ( m 1 + m 2 ) subscript 𝑚 3 subscript 𝑞 2 subscript 𝑚 1 subscript 𝑚 2 m_{3}=q_{2}(m_{1}+m_{2}) (‘high-mass tertiary’ case), and according to m 3 = q 2 ′ ​ m 2 subscript 𝑚 3 superscript subscript 𝑞 2 ′ subscript 𝑚 2 m_{3}=q_{2}^{\prime}m_{2} (‘low-mass tertiary’ case), where both q 2 subscript 𝑞 2 q_{2} and q 2 ′ superscript subscript 𝑞 2 ′ q_{2}^{\prime} have flat distributions.

	Rate (${\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}$)
	Triple			Binary
	${\sigma }_k/\mathrm{km}{s}^{-1}$			${\sigma }_k/\mathrm{km}{s}^{-1}$
	0	40	265	0	40	265
NS-NS
High-mass tertiary
RLOF	1741	551	109	3874	2119	542
Unbound	1893	104	54	—	—	—
Total	3634	655	164	3874	2119	542
Low-mass tertiary
RLOF	1937	909	309	4149	1731	627
Unbound	1856	0	0	—	—	—
Total	3793	909	309	4149	1731	627
BH-NS
High-mass tertiary
RLOF	166	257	268	186	840	759
Unbound	421	423	263	—	—	—
Total	586	680	531	186	840	759
Low-mass tertiary
RLOF	109	346	345	274	747	823
Unbound	514	21	0	—	—	—
Total	622	367	345	274	747	823
NS-MS
High-mass tertiary
RLOF	3464	5178	7006	7710	15238	20883
Unbound	1232	4097	4650	—	—	—
Total	4697	9275	11656	7710	15238	20883
Low-mass tertiary
RLOF	3928	7475	10351	8299	16088	23368
Unbound	1564	28	0	—	—	—
Total	5493	7504	10351	8299	16088	23368

Equations

$$\frac{\mathrm{d}N}{\mathrm{d}{V}_{\mathrm{k}}}=\sqrt{\frac{2}{\pi }}\frac{V_{\mathrm{k}}^2}{{\sigma }_{\mathrm{k}}^3}\exp \left(-\frac{V_{\mathrm{k}}^2}{2{\sigma }_{\mathrm{k}}^2}\right).$$

$$\displaystyle f_{\mathrm{fb}}=\left\{\begin{array}[]{cc}0,&M_{\mathrm{CO}}<5\,\mathrm{M}_{\odot};\\ 0.378M_{\mathrm{CO}}-1.889,&5.0\,\mathrm{M}_{\odot}\leq M_{\mathrm{CO}}<7.6\,\mathrm{M}_{\odot};\\ 1,&M_{\mathrm{CO}}\geq 7.6\,\mathrm{M}_{\odot},\end{array}\right.$$

$$\sqrt{1-e_1}<5\pi \frac{m_3}{m_1+m_2}{\left[\frac{a_1}{a_2\left(1-e_2\right)}\right]}^3.$$

$$M_{\mathrm{s}}={\int }_{m_{\mathrm{low}}}^{m_{\mathrm{up}}}m\frac{\mathrm{d}N}{\mathrm{d}m}\mathrm{d}m=N_{\mathrm{s}}\tilde{M},$$

$$M_{\mathrm{tot}}\simeq \left({\alpha }_{\mathrm{s}}+\frac{3}{2}{\alpha }_{\mathrm{bin}}+\frac{9}{4}{\alpha }_{\mathrm{tr}}\right)N_{\mathrm{sys}}\tilde{M}=\left(1+\frac{1}{2}{\alpha }_{\mathrm{bin}}+\frac{5}{4}{\alpha }_{\mathrm{tr}}\right)N_{\mathrm{sys}}\tilde{M}.$$

$$N$$

$$R_{\mathrm{merge}}\simeq \frac{\alpha R_{\mathrm{SFR}}}{\left(1+\frac{1}{2}{\alpha }_{\mathrm{bin}}+\frac{5}{4}{\alpha }_{\mathrm{tr}}\right)\tilde{M}}{f}_{\mathrm{calc}}{f}_{\mathrm{merge}}.$$

Full text

1. Introduction

The detection of a double neutron star (NS) merger by LIGO’s O2 observing run (Abbott et al., 2017b) has definitively shown that double neutron stars (NSs) can merge in the Universe, and the accompanying electromagnetic signals in gamma rays (e.g., Abbott et al., 2017a; Goldstein et al., 2017; Savchenko et al., 2017), X-rays (e.g., Margutti et al., 2017), optical (e.g., Coulter et al., 2017; Shappee et al., 2017; Soares-Santos et al., 2017; Nicholl et al., 2017), near-infrared (e.g., Chornock et al., 2017; Cowperthwaite et al., 2017; Pian et al., 2017), and radio wavelengths (e.g., Alexander et al., 2017), have revealed a wealth of information on kilonova transients (e.g., Lattimer & Schramm, 1974; Freiburghaus et al., 1999; Metzger et al., 2010; Rosswog, 2015; Just et al., 2015), and short gamma-ray bursts (e.g., Eichler et al., 1989; Piran, 1999; Berger, 2014). The origin of compact object mergers is unclear, since in standard stellar evolution theory, the progenitor stars would merge before evolving to compact objects in an orbit that would cause them to merge due to gravitational wave (GW) emission within a Hubble time. Two main scenarios have been proposed for compact object mergers: isolated binary evolution (e.g., Tutukov & Yungelson, 1973; Tutukov et al., 1992; Tutukov & Yungelson, 1993; Lipunov et al., 1997; Belczynski et al., 2002; Voss & Tauris, 2003; Kalogera et al., 2007; Dominik et al., 2012, 2013; Belczynski et al., 2014, 2016a, 2017; Stevenson et al., 2017; Chruslinska et al., 2018; Giacobbo & Mapelli, 2019) and dynamical interactions, such as those in triple-star systems (e.g., Thompson, 2011; Hamers et al., 2013; Antonini et al., 2017; Silsbee & Tremaine, 2017; Liu & Lai, 2017, 2018; Toonen et al., 2018; Liu et al., 2019), globular clusters (e.g., Sigurdsson & Hernquist, 1993; Portegies Zwart & McMillan, 2000; O’Leary et al., 2006; Ziosi et al., 2014; Rodriguez et al., 2015, 2016; Kimpson et al., 2016; Mapelli, 2016; Samsing & Ramirez-Ruiz, 2017; Samsing et al., 2018a, b; Samsing, 2018; Rodriguez et al., 2018), and galactic nuclei (e.g., Antonini & Perets, 2012; Antonini et al., 2014; Prodan et al., 2015; Stephan et al., 2016; VanLandingham et al., 2016; Antonini & Rasio, 2016; Petrovich & Antonini, 2017; Arca-Sedda & Gualandris, 2018; Randall & Xianyu, 2018a, b; Hamers et al., 2018; Gondán et al., 2018; Hoang et al., 2018; Fragione et al., 2018; Arca-Sedda & Capuzzo-Dolcetta, 2019).

The isolated binary channel involves binary interactions such as mass transfer, common-envelope (CE) evolution, and supernovae (SNe) kicks associated with the birth of NSs or black holes (BHs). The details of these processes are highly uncertain, yet, taking into account these uncertainties, it is challenging in the isolated binary model to predict rates (e.g., Kalogera et al. 2004b, a, see also Abadie et al. 2010 for a review) that are on the same order of magnitude as derived from the LIGO observation of GW170817, $1540_{-1220}^{+3200}\,\mathrm{Gpc^{-3}\,yr^{-1}}$ (Abbott et al., 2017b). For example, Giacobbo & Mapelli (2018) find rates of up to $1\times 10^{3}\,\mathrm{Gpc^{-3}\,yr^{-1}}$, but only assuming low SNe kicks, and a high CE efficiency.

Triple systems have been invoked in a large variety of contexts to explain an enhanced rate of interaction in binary systems (see, e.g., Naoz 2016 for a review). This occurs through oscillations of the inner binary eccentricity as a result of the torque of the tertiary, known as Lidov-Kozai (LK) oscillations (Lidov, 1962; Kozai, 1962). In particular, in the context of compact objects, it has been suggested that LK-driven eccentricity excitation can accelerate the mergers of compact objects (e.g., Blaes et al., 2002; Miller & Hamilton, 2002; Wen, 2003; Thompson, 2011; Antonini et al., 2017; Hoang et al., 2018; Bonetti et al., 2018; Stephan et al., 2019; Fragione & Loeb, 2019). However, since compact objects originate from massive stars111Excluding the possibility of primordial compact objects., the prior stellar evolution should also be taken into account, similarly to studies of binary evolution (e.g., Hamers et al., 2013; Stephan et al., 2016; Naoz et al., 2016; Toonen et al., 2018; Hamers, 2018a; Stephan et al., 2019). Natal kicks, in particular, could significantly reduce the rates by unbinding the tertiary before it can affect the inner binary through LK evolution (e.g., Lu & Naoz, 2019; Hamers & Thompson, 2019).

In this paper, we investigate the interplay between these processes in massive triple-star systems, taking into account stellar evolution starting from the main sequence (MS), secular and tidal evolution, and the effects SNe. In contrast to previous studies which focused on wide inner binaries that do not interact in the absence of a tertiary star (e.g., Hamers et al., 2013; Antonini et al., 2017; Hamers & Thompson, 2019), we here consider triples with no restrictions on the inner binary orbital separation. In particular, this implies that in the ‘binary case’ of our triples, i.e., when the effect of the tertiary on the inner binary is ignored, the system can in fact interact, and possibly produce a double NS merger via the standard binary evolution channel.

This approach requires to model binary processes such as mass transfer and CE evolution, as well as the effects associated with the gravitational perturbations from the tertiary (i.e., LK oscillations). Combining these processes is challenging and complicated (see, e.g., Hamers & Dosopoulou 2019 for an exploratory study for the case of mass transfer taking into account orbital effects due to mass transfer and LK evolution). Here, we take a simpler approach in which we model the evolution of the tertiary initially using a secular code taking into account dynamical, stellar and tidal evolution, but not mass transfer nor CE evolution. We track the onset of mass transfer (i.e., Roche-Lobe overflow, RLOF), and in this case, we continue the evolution of the inner binary system using a dedicated binary population synthesis code which includes all the required binary physics, but not any effects associated with the tertiary star. Here, we make the assumption that, after the onset of RLOF in the inner binary, the latter is dynamically decoupled from the tertiary. This assumption is typically justified in the case of CE evolution, in which case the inner binary shrinks significantly (see Fig. 8 of Hamers et al. 2013). However, in the case of mass transfer, the inner orbit can also expand (if the donor has become less massive than the accretor). We here ignore this complication, but instead take the simpler, decoupled approach.

Several types of compact object mergers have been studied in triple systems, including WD-WD mergers (Thompson, 2011; Katz & Dong, 2012; Hamers et al., 2013; Toonen et al., 2018; Stephan et al., 2019), BH-BH mergers (Silsbee & Tremaine, 2017; Antonini et al., 2017; Liu & Lai, 2017, 2018; Hoang et al., 2018; Liu et al., 2019), and BH-NS mergers (Stephan et al., 2019; Fragione & Loeb, 2019). To our knowledge, NS-NS mergers in triples with stellar-mass tertiaries (taking into account stellar evolution and dynamical evolution; see Stephan et al. 2019 for a comparable study, but with supermassive BH tertiaries) have not been studied. Our focus is therefore on mergers of double NSs.

The structure of this paper is as follows. We describe our methodology in Section 2, and the initial conditions of our simulations in Section 3. We present our results, most notably the merger rates, in Section 4. We discuss our findings in Section 5, and conclude in Section 6.

2. Numerical method

We use a hybrid method in which we use SecularMultiple (Hamers & Portegies Zwart, 2016a; Hamers, 2018b) to model the secular dynamical, tidal, and stellar evolution of a binary or hierarchical triple star system, and the binary stellar evolution code BSE (Hurley et al., 2000, 2002) to model the evolution of systems in which we consider the tertiary to be unimportant. The latter case, which we will refer to as ‘isolated’ binary evolution, includes interacting systems that undergo mass transfer in the inner binary, and systems in which the tertiary star becomes unbound from the inner binary due to a SNe event but with the inner binary remaining bound, and possibly merging at a later time. Both codes, SecularMultiple and BSE, are implemented within the AMUSE framework (Portegies Zwart et al., 2013; Pelupessy et al., 2013).

2.1. SecularMultiple

In SecularMultiple, we model the evolution of a binary or triple system, starting from MS stars, and taking into account stellar evolution, tidal evolution, and secular dynamical evolution (in the case of triples, and up to and including fifth order in the expansion of the separation ratio of the inner to the outer binary, see Hamers & Portegies Zwart 2016a; Hamers 2018b222The reference frame in SecularMultiple is an arbitrary frame, rather than the invariable plane.). The modeling in SecularMultiple is similar to that in previous works in which we coupled secular dynamical evolution to stellar and tidal evolution (Hamers & Portegies Zwart, 2016b; Hamers, 2018a; Hamers & Thompson, 2019).

Stellar evolution is taken into account by using SSE (Hurley et al., 2000) (as implemented in AMUSE), which is based on analytical fits to detailed stellar evolution models, and which uses the same stellar tracks as (i.e., is consistent with) BSE. We set the metallicity to Solar, $Z=0.02$. Quantities that are used from SSE include the stellar type and (convective envelope) mass and radius as a function of age. These quantities are used to take into account mass loss from the system, assumed to occur either adiabatically, or because of an impulsive change due to SNe. In the former case, we assume isotropic and adiabatic mass loss, i.e., $a_{i}M_{i}$ and $e_{i}$ are constant (Huang, 1956, 1963) for an orbit $i$ in the system ($i=1$ and $i=2$, for the inner and outer orbits, respectively, and if applicable). Here, $M_{1}=m_{1}+m_{2}$ for the inner orbit, and $M_{2}=M_{1}+m_{3}$ for the outer orbit.

In the case of SNe, we use the routines incorporated into SecularMultiple and described in Hamers (2018b) to compute the effects on the inner and outer orbits of the (assumed to be instantaneous) mass loss and (possible) kick to the newly-formed NS or BH, assuming a random orbital phase of both orbits at the moment of SNe. We assume that the kick distribution is a Maxwellian distribution, i.e., the probability density function of the kick speed $V_{\mathrm{k}}$ is given by

[TABLE]

Here, $\sigma_{\mathrm{k}}$ characterizes the typical kick speed, and is assumed to be given (see Section 3 below). For BHs, we use the prescription of Fryer et al. (2012, Section 4.1) to rescale the sampled speed $V_{\mathrm{k}}$ in equation (1) according to $V_{\mathrm{k},\,\mathrm{BH}}=V_{\mathrm{k}}(1-f_{\mathrm{fb}})$, where the fallback factor $f_{\mathrm{fb}}$ is given by

[TABLE]

where $M_{\mathrm{CO}}$ is the mass of the carbon-oxygen core of the proto-BH, which is extracted from SSE.

Tidal evolution is taken into account with the assumption of the equilibrium tide model (Hut, 1981; Eggleton et al., 1998). Specifically, we use equations (81) and (82) of Eggleton et al. (1998), with the non-dissipative terms $X$, $Y$ and $Z$ given explicitly by equation (10)-(12) of Eggleton & Kiseleva-Eggleton (2001), and the dissipative terms $V$ given explicitly by equations (A7)-(A11) of Barker & Ogilvie (2009). We apply these terms to both the inner and outer orbits (if applicable), in the latter case treating the inner binary as a point mass. The stellar spins are included in the tidal calculations; the initial spin-orbit angles are assumed to be zero (i.e., zero obliquities). The tidal dissipation efficiency is computed as a function of the stellar parameters as part of the set of ordinary differential equations using the prescription of Hurley et al. (2002). For all stars, we assume a fixed apsidal motion constant of $k_{\mathrm{AM},i}=0.014$, and a gyration radius of 0.08 (these parameters are not provided by SSE in AMUSE).

Many uncertainties still remain in tidal evolution (see, e.g., Ogilvie 2014 for a review). Especially for high eccentricities, which could be excited by LK oscillations, the equilibrium tide model could break down. A more sophisticated treatment of tides is beyond the scope of this work.

During the integration of the system in SecularMultiple, we check for a number of stopping conditions. These are listed and explained below.

For both binaries and triples, we check for RLOF in the inner binary using analytic fits (Sepinsky et al., 2007, eqs. 47-52). Specifically, we check for RLOF at periapsis as a function of the current radii, spins, and inner binary semimajor axis and eccentricity. If RLOF occurs, we stop the simulation within SecularMultiple, and continue the evolution of the inner binary using BSE; see Section 2.2 below for further details. We also check for RLOF by the tertiary star on the inner binary, treating the binary as a point mass and using the same fits of Sepinsky et al. (2007), although this scenario occurs much less commonly compared to RLOF in the inner binary (see Tables 2 and 3). 2. 2.

For both binaries and triples, after each SNe event, we check if the system is still bound (i.e,. if $a_{i}>0$ for $i=1$ and $i=2$ if applicable). If not, we continue the evolution of the inner binary with BSE if the outer orbit became unbound, but the inner binary is still bound (and could interact and merge at a later time); see Section 2.2 below. 3. 3.

For triples, we check for dynamical stability of the system using the criterion of Mardling & Aarseth (2001), and stop the simulation in SecularMultiple if the system is unstable according to this criterion. We do not consider the subsequent evolution of the system. 4. 4.

For triples, we check if at any point in the evolution if the secular approximation made in the equations of motion breaks down. Such a breakdown can occur if the timescale for the inner binary angular momentum to change significantly becomes comparable to the inner or outer orbital periods (e.g., Ćuk & Burns 2004; Ivanov et al. 2005; Katz & Dong 2012; Antonini & Perets 2012; Seto 2013; Antonini et al. 2014; Antognini et al. 2014; Bode & Wegg 2014; Luo et al. 2016; Lei et al. 2018; Grishin et al. 2018). To evaluate whether or not this has occurred, we use the criterion of Antonini et al. (2014), i.e.,

[TABLE]

Although we check for this regime (also known as the semisecular regime), it occurs very rarely in our simulations since extremely high eccentricities are required, and in most cases this implies that the stars in the inner binary would undergo mass transfer well before. 5. 5.

The inner binary components collide directly, i.e., $a_{1}(1-e_{1})<R_{1}+R_{2}$, where $R_{1}$ and $R_{2}$ are the primary and secondary radius, respectively. This typically can only occur for compact objects, since the larger sizes of pre-compact objects generally imply that, when the eccentricity is high, usually RLOF occurs or possibly the semisecular regime is triggered. 6. 6.

The system exceeds an age of 10 Gyr.

In contrast to previous work related to white dwarfs (Hamers, 2018a; Hamers & Thompson, 2019), we here do not include the effects of flybys on the binary or triple system. We find that most NS-NS mergers occur relatively early, i.e., within $\sim 100\,\mathrm{Myr}$ (see Section 4 below). For systems in the field, flybys are typically unimportant on such time-scales. Mergers occurring as a result of perturbations from passing stars are beyond the scope of this paper (see, e.g., Kaib & Raymond 2014; Michaely & Perets 2019).

2.2. BSE

As described in Section 2, we consider two situations in the simulations with SecularMultiple in which we assume for triples that the inner binary can subsequently be decoupled from the tertiary star. We then continue the evolution with the binary population synthesis code BSE. The latter code includes prescriptions for binary-star evolution, most notably mass transfer and CE evolution, which are both not modeled in SecularMultiple. In BSE, we assume a CE binding energy factor of $\lambda=0.5$, and a CE efficiency of $\alpha=1$. Furthermore, in BSE we assume the same NS kick speed distribution as in SecularMultiple (see equation 1). Due to limitations of BSE, we do not include the fallback correction to kicks applied to BHs, as described by equation (5). Instead, in BSE, we assume zero kick speeds for BHs. This is strictly inconsistent with the prior simulations in SecularMultiple, but our main focus in this work is on merging NSs, so we do not expect this inconsistency to strongly affect our conclusions.

The two situations in which we hand off the evolution in SecularMultiple to BSE correspond to stopping conditions (1) and (2) described in Section 2.1. In situation (1), the inner binary undergoes RLOF, which can be triggered by several factors: expansion of the stars in the inner binary due to stellar evolution, possibly combined with tidal evolution, and/or high eccentricity due to secular evolution. Note that, contrary to previous studies, we do not restrict to wide systems. Therefore, many systems will undergo RLOF even in the ‘binary’ case, i.e., in absence of a tertiary star.

The crucial assumption made in situation (1) for triples is that the subsequent evolution of the inner binary is completely decoupled from the tertiary star. This is a strong assumption, and is certainly not correct for all systems. The motivation for this approach is simplicity, since we currently do not have the tools to self-consistently simulate the evolution of the triple system taking into account the secular dynamical evolution with mass transfer, in particular if the orbit is eccentric and driven to varying eccentricity due to LK oscillations with similar RLOF and LK time-scales. Preliminary work (Hamers & Dosopoulou, 2019) indicates that the evolution in this regime is complicated, and merits future study. The decoupling assumption likely gives rise to a significant systematic uncertainty in our results. We further discuss this issue in Section 5.1.

3. Initial conditions

Here, we describe the initial conditions assumed in our simulations. We consider several sets of initial conditions; in each set, we generate $N_{\mathrm{MC}}=10^{3}$ systems through Monte Carlo sampling. Each set is characterized by the type of system – ‘Triple’ or ‘Binary’, where ‘Binary’ refers to the same system, but in absence of the tertiary star –, the dispersion $\sigma_{\mathrm{k}}$ assumed in the kick speed distribution (see equation 1), and the type of assumption on the tertiary mass distribution (see Section 3.1 below).

Below, we describe in more detail the assumptions made in the Monte Carlo sampling. An overview of our notation and assumptions is given in Table 1.

3.1. Masses

We sample the mass of the primary star, $m_{1}$, from a Salpeter distribution (Salpeter, 1955), i.e., $\mathrm{d}N/\mathrm{d}m_{1}\propto m_{1}^{-2.35}$. Since our interest is in NS mergers, we sample with the range $8<m_{1}/\mathrm{M}_{\odot}<50$. The mass of the secondary star, $m_{2}$ is sampled assuming a flat distribution in $q_{1}\equiv m_{2}/m_{1}$, consistent with observations of massive stars (e.g., Sana et al., 2012; Duchêne & Kraus, 2013; Kobulnicky et al., 2014). Here, we set the lower limit on $m_{2}$ to be $4\,\mathrm{M}_{\odot}$. The tertiary mass $m_{3}$ is sampled according to two methods: sampling from (1) a flat distribution in the mass ratio $q_{2}\equiv m_{3}/(m_{1}+m_{2})$, and (2) a flat distribution in the mass ratio $q_{2}^{\prime}\equiv m_{3}/m_{2}$. In either case, we set the lower limit on $m_{3}$ to be $0.1\,\mathrm{M}_{\odot}$. These two choices give different results in our simulations, since in case (1) the tertiary star can be more massive than the primary star, and therefore can evolve first. Consequently, the tertiary star can become unbound from the inner binary as it explodes in a SNe event. In contrast, in method (2), the tertiary star is always less massive than the primary star, such that the inner binary always evolves first.

In Fig. 1, we illustrate the difference between the two assumptions on the tertiary mass ratio distribution. With method (1), the median tertiary mass is $\approx 5\,\mathrm{M}_{\odot}$, whereas with method (2) it is $\approx 12\,\mathrm{M}_{\odot}$.

3.2. Orbits

For both inner and outer orbits, we assume a flat distribution in the logarithm of the orbital period, i.e., flat in $\log_{10}(P_{\mathrm{orb},i})$ (also known as an Öpic distribution, Öpik 1924). The limits for both inner and outer orbits are set to $1<\log_{10}(P_{\mathrm{orb},i}/\mathrm{d})<10$. The eccentricities $e_{i}$ of both inner and outer orbits are sampled from flat distributions in $e_{i}$, with $0<e_{j}<0.9$. The orbital period and eccentricity distributions are roughly consistent with observations of massive binary stars (Kobulnicky et al., 2014). Given the uncertainties in the observed orbital distributions of massive triples, we do not take into account a more sophisticated initial distribution, even though the observed distribution for massive MS triples reflects the initial distribution (Rose et al., 2019).

We reject a sampled system if it is unstable according to the criterion of Mardling & Aarseth (2001). The orientations of the orbits are taken to be random, i.e., the inclinations $i_{i}$ are sampled from flat distributions in $\cos i_{i}$, and the arguments of periapsis $\omega_{i}$ and longitudes of the ascending node $\Omega_{i}$ are sampled from flat distributions. There are indications that lower-mass triples have more aligned inner and outer orbits if the system is compact ($a_{2}\lesssim 100\,\,\textsc{au}$), whereas the wider systems are more isotropically distributed (Tokovinin, 2017). However, it is unclear if this trend also persists for higher-mass triples (primaries with masses $\gtrsim 8\,\mathrm{M}_{\odot}$). For simplicity, we here restrict to randomly-oriented triples.

3.3. Kick distributions

As discussed in Section 2, we sample the kick distribution for SNe from a Maxwellian distribution (e.g., Hansen & Phinney, 1997). We take the dispersion $\sigma_{\mathrm{k}}$ to be fixed for each set of Monte Carlo simulations, and adopt three values: 0, 40, and 265 $\mathrm{km\,s^{-1}}$. The value $\sigma_{\mathrm{k}}=0\,\mathrm{km\,s^{-1}}$ is to evaluate the importance of mass loss associated with SNe only (i.e., only the Blaauw kick, Blaauw 1961; Boersma 1961). The value $\sigma_{\mathrm{k}}=265\,\mathrm{km\,s^{-1}}$ is a commonly-adopted value inferred from the proper motions of pulsars (Hobbs et al., 2005). However, NS kicks are uncertain and their magnitude is still highly debated. For example, Arzoumanian et al. (2002) found a two-component distribution of kick speed distributions with characteristic velocities of $90$ and $500\,\mathrm{km\,s^{-1}}$ based on the velocities of isolated radio pulsars, and Beniamini & Piran (2016) also found evidence for a bimodal distribution based on observed binary NSs, with a low-kick population with $V_{\mathrm{k}}<30\,\mathrm{km\,s^{-1}}$, and a high-kick population with kicks up to $400\,\mathrm{km\,s^{-1}}$. Rather than adopting a bimodal distribution, we here choose to carry out another set with $\sigma_{\mathrm{k}}=40\,\mathrm{km\,s^{-1}}$, to evaluate the importance of low but non-zero kicks.

4. Results

4.1. Outcome fractions

4.1.1 Main channels

In Tables 2 and 3, we show the fractions for the main channels in our simulations assuming a flat distribution in $q_{2}\equiv m_{3}/(m_{1}+m_{2})$ and $q_{2}^{\prime}\equiv m_{3}/m_{2}$, respectively (henceforth, we refer to the latter two assumptions on the tertiary mass ratio as the ‘high mass tertiary’ and ‘low mass tertiary’ cases, respectively). These data are based on simulations with SecularMultiple only. The first three data columns correspond to triple systems, whereas the last three data columns correspond to the binary case, i.e., taking the same triple systems but in the absence of the tertiary star. For each case, we show results for the three different values of the kick velocity dispersion, $\sigma_{\mathrm{k}}$.

The channels shown in these tables include ‘No interaction’, i.e,. the triple or binary survived for 10 Gyr without triggering interaction such as mass transfer, or dynamical instability or instability as a result of SNe. In most cases, this ‘inert’ channel is unlikely; instead, much more common are RLOF, or the unbinding of the system due to SNe. The only notable exception is for the binary case with $\sigma_{\mathrm{k}}=0\,\mathrm{km\,s^{-1}}$, in which case the non-interacting fraction is $\sim 0.2$. In the equivalent triple case, the non-interacting fraction is only $\sim 0.02$, i.e., a factor 10 times smaller. This can be attributed largely due to the Blaauw kick in the inner binary which can keep the inner binary bound, but make the tertiary unbound (compare the ‘NS+MS’ Unbound fractions between the triple and binary cases with $\sigma_{\mathrm{k}}=0\,\mathrm{km\,s^{-1}}$).

For RLOF, we distinguish between RLOF of the primary, secondary, or tertiary star (in the latter case, the inner binary is treated as a point mass in the fits of Sepinsky et al. 2007). RLOF is predominantly triggered by the primary star, with a fraction of $\sim 0.6$ assuming $\sigma_{\mathrm{k}}=0\,\mathrm{km\,s^{-1}}$, and $\sim 0.4$ assuming $\sigma_{\mathrm{k}}>0\,\mathrm{km\,s^{-1}}$ for the triple systems with higher-mass tertiaries. The decrease in the RLOF fraction with increasing $\sigma_{\mathrm{k}}$ in the latter case can be ascribed to the higher fraction of unbound systems, which in turn is mostly due to kicks imparted on the tertiary star when it evolves first (see the unbound MS+MS fractions for triples in Table 2). This trend is much less pronounced in the case of lower-mass tertiaries (Table 3) – in this case, the tertiary, which is always as massive or less massive than the secondary star, does not evolve first, and therefore is less likely to unbind the triple system due to its SNe kick.

In contrast to triples, for binaries, the RLOF star 1 fraction is independent of $\sigma_{\mathrm{k}}$. This can be explained by noting that for binaries, the primary star is always the most massive and evolves the fastest; whether or not star fills its Roche lobe is completely determined by the initial masses ($m_{1}$ and $m_{2}$), $a_{1}$, and $e_{1}$, and independent of $\sigma_{\mathrm{k}}$. This is no longer the case for the secondary star, e.g., RLOF of the secondary star can be triggered by a SNe event of the primary star, whose properties in turn are set by $\sigma_{\mathrm{k}}$.

It may be surprising that the RLOF fraction is typically $\sim 0.6$, for both triple and binary cases. It might be expected that the RLOF fraction would be much higher for triple systems, since RLOF can be triggered by LK evolution. To explain this, we note that in most previous studies of the onset of RLOF in triples (e.g., Hamers et al., 2013; Antonini et al., 2017; Hamers & Thompson, 2019), the inner binaries were assumed to be wide enough to avoid interaction in the absence of a tertiary star. In contrast, we here include tight systems as well. Therefore, a significant fraction of systems already interact in the ‘binary case’. In addition, as shown by the higher unbound fractions for the triple cases compared to the binary cases, orbital changes due to SNe associated with the tertiary star play an important role in triples with NSs.

The fraction of dynamically unstable systems for triples is relatively small, with the fraction being at most $\sim 0.03$ if $\sigma_{\mathrm{k}}=0\,\mathrm{km\,s^{-1}}$. Dynamical instability is typically triggered by mass loss from the inner binary, i.e., the traditional triple dynamical instability scenario (Perets & Kratter, 2012). Instability can also be triggered by SNe events, but this is a rare event. The dynamical instability fraction in fact decreases with increasing $\sigma_{\mathrm{k}}$, which can be understood from the larger fraction of systems becoming unbound due to SNe before a dynamical instability can be triggered.

As discussed in Section 2, we also check for the semisecular regime in our simulations. This regime is triggered only very rarely, and predominantly with the inner binary consisting of two MS stars. The reason for the rarity of this channel is that very high eccentricities are required to trigger it, which usually instead lead to RLOF in eccentric orbits.

We also check for direct collisions (indicated with ‘Secular collision’ in the tables). These do not occur for pre-compact objects, since RLOF is expected to ensue before direct collision. However, direct collisions could occur during later stages, when compact objects have formed, and when the inner orbit is excited in eccentricity due to secular evolution (e.g., Thompson, 2011; Katz & Dong, 2012). Nevertheless, we find no such direct collisions of compact objects. This can be attributed to the high fractions of systems that undergo RLOF or become unbound. In other words, the probability that the system survives without interacting or becoming unstable due to SNe and a collision is triggered at later stages, is small. We emphasize that our simulations are limited in terms of the number of systems. We therefore cannot exclude that direct collisions would occur if $N_{\mathrm{MC}}$ were increased. However, we here focus on the largest contribution to NS-NS mergers, which we find originate from interacting or unbound systems (i.e., RLOF-induced and ‘unbound’ mergers). The latter are discussed in further detail below.

4.1.2 Mergers from RLOF systems

As described in detail in Section 2.2, we continue the evolution of the inner binary after RLOF occurs in the inner binary using the binary population synthesis code BSE. In Tables 5 and 5 for the high- and low-mass tertiary simulations, respectively, we show the fractions for merger outcomes of the subsequent ‘isolated binary’ simulations. These fractions are with respect to the systems undergoing RLOF triggered by the primary or secondary star; the fractions of the latter cases corresponding to all systems were given in Tables 2 and 3 for the high- and low-mass tertiary cases, respectively.

The evolution of the inner binary after RLOF in BSE is driven by several processes, including mass transfer, CE evolution, and SNe. There is a high likelihood, of $\sim 0.7$ to $\sim 0.9$ of the RLOF systems depending on $\sigma_{\mathrm{k}}$, that the inner binary eventually merges. It is evident from Tables 5 and 5 that a large number of merger outcomes are possible. Generally, there are little differences in the fractions between the ‘triple’ and ‘binary’ cases. This can be attributed to the generally only small differences in the properties of the inner binaries when RLOF is triggered (see Section 4.2 below).

A dominant merger channel is a NS-MS merger, with a fraction of up to $\sim 0.6$ of RLOF systems if $\sigma_{\mathrm{k}}=265\,\mathrm{km\,s^{-1}}$. The resulting stars, Thorne-$\mathrm{\dot{Z}}$ytkow objects (Thorne & Zytkow, 1977), are also found in isolated binary evolution studies (e.g., Brandt & Podsiadlowski, 1995); our results show that this channel is also possible (and relatively likely) in triples. These objects could also evolve to become X-ray binaries, which are also believed to form in triples (without taking into account isolated binary evolution) through the mass-loss induced eccentric LK mechanism (Shappee & Thompson, 2013; Naoz et al., 2016).

The channel of most interest here is the merger of two NSs, which occurs for a relatively large fraction of $\sim 0.1$ of RLOF systems if $\sigma_{\mathrm{k}}=0\,\mathrm{km\,s^{-1}}$, but drops quickly to a tenth of this fraction, to $\sim 0.01$, if $\sigma_{\mathrm{k}}=265\,\mathrm{km\,s^{-1}}$. As expected, the kick dispersion has a large impact on the fraction of NS-NS mergers. Also, these fraction are mostly sensitive to $\sigma_{\mathrm{k}}$, and do not depend very strongly on the assumption of the tertiary mass ratio distribution (it should be taken into account, however, that the latter does affect the overall fraction of RLOF systems).

4.1.3 Mergers from systems with unbound tertiaries

In addition to considering the subsequent ‘isolated binary’ evolution of systems that undergo RLOF, we also consider the ‘isolated binary’ evolution of systems in which the tertiary becomes unbound from the system due to a SNe event, but the inner binary remains bound. In the latter case, the inner binary can subsequently merge due to ‘isolated binary’ evolution, which we take into account by evolving these systems with BSE (see Section 2.2). Similarly to Tables 5 and 5, we show in Tables 7 and 7 the outcome fractions of several merger channels for these ‘unbound tertiary’ systems. Evidently, in this case there are no equivalent ‘binary’ systems (since this channel originates exclusively from triples). The fractions in these tables are given with respect to the systems in which the tertiary becomes unbound, but with a bound inner binary; the latter fractions with respect to all systems are given in the bottom rows in Tables 2 and 3 for the high-mass and low-mass tertiaries, respectively. Note that the latter fraction is very small for the low-mass tertiary simulations and non-zero $\sigma_{\mathrm{k}}$.

Similarly to RLOF-induced mergers, mergers originating from unbound tertiary systems are dominated by NS-MS mergers, in particular for non-zero $\sigma_{\mathrm{k}}$. The fraction of NS-NS mergers is relatively high at $\sim 0.25$ for $\sigma_{\mathrm{k}}=0\,\mathrm{km\,s^{-1}}$, but rapidly decreases with increasing $\sigma_{\mathrm{k}}$, dropping to $\sim 0.01$ for $\sigma_{\mathrm{k}}=265\,\mathrm{km\,s^{-1}}$ in the high-mass tertiary case, and to zero (within our statistical certainty) for $\sigma_{\mathrm{k}}=265\,\mathrm{km\,s^{-1}}$ in the low-mass tertiary case.

4.2. Orbital properties

We further discuss the RLOF-induced (Section 4.1.2) and tertiary unbound (Section 4.1.3) channels by showing in Figs 2 and 3 the distributions of the inner binary semimajor axis $a_{1}$ and eccentricity $e_{1}$ at the moment of the onset of RLOF, and the unbinding of the tertiary star. Here, Figs 2 and 3 correspond to the high- and low-mass tertiary cases, respectively.

As expected, RLOF-induced systems tend to have significantly smaller inner binary semimajor axes compared to all systems (compare the red and black lines in the figures), with RLOF systems having a median of $a_{1}\sim 1\,\,\textsc{au}$, compared to $a_{1}\sim 10\,\,\textsc{au}$ for all systems overall. Also, SNe kicks have very little impact on the orbital properties RLOF systems, as expected (since in our simulations RLOF typically occurs before stars evolve to compact objects). In addition, for the RLOF systems there are little differences in terms of $a_{1}$ between the triple and binary cases. There is some difference in terms of the eccentricity – in the triple case, RLOF-induced systems tend to have slightly higher eccentricities compared to the binary case. This can be explained by eccentricity excitation by LK oscillations. There are no noticeable differences in the orbital properties of RLOF systems between the high- and low-mass tertiary simulations.

The ‘tertiary unbound’ systems show typically larger semimajor axes compared to all systems (compare the green and black lines in the figures). With higher kicks, the semimajor axes tend to be larger compared to without ($\sigma_{\mathrm{k}}=0\,\mathrm{km\,s^{-1}}$); this can be understood by noting that kicks tend to unbind the inner binary if it is wide, so, for the inner binary to remain bound, the inner binary semimajor axis should be smaller. Note that the number of unbound systems in the low-mass tertiary simulations for $\sigma_{\mathrm{k}}>0\,\mathrm{km\,s^{-1}}$ are very small (see the bottom row in Table 3), causing the jagged behavior in Fig. 3.

4.3. Delay-time distributions

In Figs 4 and 5, we show delay-time distributions (DTDs) for various merger channels for the triple and binary cases, respectively. We include the RLOF-induced (red lines), and tertiary unbound (dashed green lines) cases. Each panel corresponds to a different channel, labeled in the top left, and we show results only for channels with a significant number of systems ($\geq 5$). Data apply to the simulations with $\sigma_{\mathrm{k}}=40\,\mathrm{km\,s^{-1}}$, and the high-mass tertiary assumption.

Generally, RLOF-induced mergers occur earlier than tertiary unbound mergers. This can be attributed to the generally tighter orbits (see Section 4.2). There are generally no large differences between the triple (Fig. 4) and binary (Fig. 5) RLOF-induced DTDs, which can be explained by the similarities in the initial orbits.

Of most interest here are the NS-NS mergers. Most of these mergers occur within $\sim 100\,\mathrm{Myr}$, although there is a small tail with delay times of several Gyr.

The DTDs for simulations with different parameters (low-mass tertiary, and different $\sigma_{\mathrm{k}}$; not shown here) are qualitatively similar to those in Figs 4 and 5.

4.4. Merger rates

Here, we compute the NS-NS merger rates based on the merger fractions from our simulations. To convert merger fractions to absolute rates, we start with the local star formation rate (SFR) density, which we take to be $R_{\mathrm{SFR}}=0.025\times 10^{9}\,\mathrm{M}_{\odot}\,\mathrm{Gpc^{-3}\,yr^{-1}}$ (Bothwell et al., 2011). This is the mass of all stars formed per unit volume and time (we are agnostic about the type of galaxy in which the stars are formed and consider the local Universe only; this is justified by the typically short delay times, see Section 4.3). We assume that all stellar systems consist of either single, binary or triple stars (and ignore high-order systems). Consider a population of stars with $N_{\mathrm{sys}}$ stellar systems. The number of single, binary, and triple systems is then $\alpha_{\mathrm{s}}N_{\mathrm{sys}}$, $\alpha_{\mathrm{bin}}N_{\mathrm{sys}}$, and $\alpha_{\mathrm{tr}}N_{\mathrm{sys}}$, respectively, and where $\alpha_{\mathrm{s}}+\alpha_{\mathrm{bin}}+\alpha_{\mathrm{tr}}=1$. We assume $\alpha_{\mathrm{s}}=0.19$, $\alpha_{\mathrm{bin}}=0.56$, and $\alpha_{\mathrm{tr}}=0.25$ (Sana et al., 2014), independent of mass. We remark that other recent observational studies (e.g., Moe & Di Stefano, 2017) indicate that these fractions could be more heavily biased towards more high-multiplicity systems for higher primary masses.

We assume a Kroupa (2001) initial mass function (IMF), i.e., $\mathrm{d}N/\mathrm{d}m\propto m^{-\alpha}$, where $\alpha=0.3$ for $0.01<m/\mathrm{M}_{\odot}<0.08$, $\alpha=1.3$ for $0.08<m/\mathrm{M}_{\odot}<0.5$, and $\alpha=2.3$ for $m>0.5\,\mathrm{M}_{\odot}$333The slope of $-2.3$ here is inconsistent with the assumed slope for the massive stars in our simulations (primary masses $8<m_{1}/\mathrm{M}_{\odot}<50$; slope $-2.35$). However, the difference in slope is very small. . A single-star population with $N_{\mathrm{s}}$ stars with such an IMF has a total mass of

[TABLE]

where $N_{\mathrm{s}}$ is the number of (single) stars, $m_{\mathrm{low}}=0.01\,\mathrm{M}_{\odot}$ and $m_{\mathrm{up}}=150\,\mathrm{M}_{\odot}$, and we calculate $\tilde{M}$ to be $\tilde{M}\simeq 0.38\,\mathrm{M}_{\odot}$ for a Kroupa (2001) IMF. In our mixed population with single, binary, and triple stars, the mass contribution from single stars is therefore $\alpha_{\mathrm{s}}N_{\mathrm{sys}}\tilde{M}$. Assuming a flat mass ratio distribution (i.e., flat in $q_{1}\equiv m_{2}/m_{1}$), the binaries have a mass contribution which is approximately $\left(1+\frac{1}{2}\right)\tilde{M}\alpha_{\mathrm{bin}}N_{\mathrm{sys}}=\frac{3}{2}\alpha_{\mathrm{bin}}N_{\mathrm{sys}}\tilde{M}$. For triples, the mass contribution depends on our assumption on the tertiary mass ratio; we assumed distributions that are flat in either $q_{2}\equiv m_{3}/(m_{1}+m_{2})$, or in $q_{2}^{\prime}=m_{3}/m_{2}$. In the former case, the mass contribution from triples is $\left(\frac{3}{2}+\frac{1}{2}\frac{3}{2}\right)\tilde{M}\alpha_{\mathrm{tr}}N_{\mathrm{sys}}=\frac{9}{4}\alpha_{\mathrm{tr}}\tilde{M}N_{\mathrm{sys}}$. In the latter case, it is $\left(\frac{3}{2}+\frac{1}{2}\frac{1}{2}\right)\tilde{M}\alpha_{\mathrm{tr}}N_{\mathrm{sys}}=\frac{7}{4}\alpha_{\mathrm{tr}}\tilde{M}N_{\mathrm{sys}}$. Adding up the contributions from all hierarchies, the total mass of our population is (in the high-mass tertiary case)

[TABLE]

The number of binaries/triples is therefore given by

[TABLE]

where $\alpha=\alpha_{\mathrm{bin}}$ for binaries, and $\alpha=\alpha_{\mathrm{tr}}$ for triples. The corresponding occurrence rates (per unit volume and time) are obtained by replacing $M_{\mathrm{tot}}$ in equation (9) by $R_{\mathrm{SFR}}$.

To convert the total occurrence rate of binaries/triples to (NS-NS) merger rates, we need to take into account the fraction of systems that were taken into account in the simulations compared to all astrophysically-occurring systems (which we denote as $f_{\mathrm{calc}}$), and the actual merger fractions in the simulations (which we denote as $f_{\mathrm{merge}}$). Since we did not restrict the orbital separation distributions in the Monte Carlo simulations (see Section 3), the calculated fraction $f_{\mathrm{calc}}$ is determined solely by the mass cutoffs. Specifically, we restricted the primary mass to the range $8<m_{1}/\mathrm{M}_{\odot}<50$ and the secondary mass to $4<m_{2}/\mathrm{M}_{\odot}<50$, whereas there was no restriction on $m_{3}$ (within the mass ranges of the IMF of Kroupa 2001). For the adopted Kroupa (2001) IMF, this implies that $f_{\mathrm{calc}}\simeq 2.6\times 10^{-3}$, which applies to both triples and binaries. The merger fractions $f_{\mathrm{merge}}$ can be inferred from Tables 2, 3, 5, 5, 7, and 7.

The merger rate for binaries/triples (number per unit volume and time) is then given by

[TABLE]

Note that this equation applies to the high-mass tertiary case; in the low-mass tertiary case, the factor $\frac{5}{4}$ in equation (10) should be replaced by $\frac{3}{4}$. The resulting NS-NS merger rates for all our simulations are given in Table 8. We also give the BH-NS merger rates in the same table. We do not include BH-BH rates, since the associated fractions $f_{\mathrm{merge}}$ in the simulations are at the noise level. The NS-NS DTD distributions normalized to the rates are shown in Fig. 6, and the total rates are plotted as a function of $\sigma_{\mathrm{k}}$ in Fig. 7. We discuss our rates further in Section 5.2.

5. Discussion

5.1. The decoupling assumption

As mentioned in Section 2.2, a (strong) assumption that we made in our modeling is that, after RLOF is triggered in the inner binary system, the subsequent evolution of the inner binary is completely decoupled from the tertiary star. This approach was taken for simplicity, and current lack of the available framework to model this phase self-consistently. Here, we briefly investigate this assumption, by considering the systems that, after RLOF is triggered in the inner binary, survive the binary evolution, i.e., do not merge within 10 Gyr.

To evaluate whether or not the inner system can be decoupled from the tertiary, we compute for these systems the ratio of the timescale $t_{\mathrm{1PN}}$ for the lowest-order post-Newtonian (PN) precession (e.g., Weinberg, 1972), to the LK timescale $t_{\mathrm{LK}}$ (e.g., Kinoshita & Nakai, 1999; Antognini, 2015). If $t_{\mathrm{1PN}}\ll t_{\mathrm{LK}}$, we expect 1PN precession to dominate and the decoupling assumption to be justified, whereas if $t_{\mathrm{1PN}}\gg t_{\mathrm{LK}}$, this is no longer the case (e.g., Blaes et al., 2002; Wu & Murray, 2003; Fabrycky & Tremaine, 2007; Thompson, 2011; Naoz et al., 2013; Liu et al., 2015).

In our simulations, we do not model the subsequent evolution of the tertiary star and orbit after the onset of RLOF. Here, we assume, for simplicity, that the tertiary mass does not change between the onset of RLOF and 10 Gyr, and make two assumptions on the outer orbit semimajor axis $a_{2}$ after 10 Gyr: we either take it to be the value at the onset of RLOF, or the value if all mass in the inner binary between the time of the onset of RLOF and 10 Gyr were lost adiabatically, i.e., with $a_{2}(m_{1}+m_{2}+m_{3})$ constant.

In Fig. 8, we show the resulting distributions of the ratio $t_{\mathrm{1PN}}/t_{\mathrm{LK}}$, for two values of $\sigma_{\mathrm{k}}$, and the high-mass tertiary simulations. As shown, there is a non-negligible fraction of systems with $t_{\mathrm{1PN}}/t_{\mathrm{LK}}\gg 1$, in which case the decoupling assumption is not well justified. The fraction of these systems depends on $\sigma_{\mathrm{k}}$; the distribution of $t_{\mathrm{1PN}}/t_{\mathrm{LK}}$ shifts to larger values with increasing $\sigma_{\mathrm{k}}$.

This indicates that not all systems are well decoupled from the tertiary. This implies there is a systematic error in our results, in particular of the NS-NS merger rates. We also note that we here only considered the surviving systems; the systems undergoing RLOF may also not be truly decoupled. At this point, it is difficult to estimate whether fully self-consistent modeling would lead to lower or higher merger rates. Such an endeavor is left for future work.

5.2. Rates

Here, we comment on our NS-NS merger rates (see Section 4.4). First, as also discussed above, it should be noted that there are likely significant systematic errors in our rates as a consequence of simplifications made in the modeling. Therefore, our results should be interpreted as estimates.

Interestingly, and most importantly, the triple rates, although always lower, are typically comparable to the binary rates (see Table 8). Both the triple and binary rates are highly sensitive to the SNe kick dispersion $\sigma_{\mathrm{k}}$, as expected. However, the triple rates are even more sensitive to $\sigma_{\mathrm{k}}$ than the binary rates. This can be explained by the fact that a third star adds more possibilities for the system to become unbound, and this effect increases in importance with increasing $\sigma_{\mathrm{k}}$. It should also be mentioned that we assumed a conservative triple fraction ($\alpha_{\mathrm{tr}}=0.25$) compared to the binary fraction ($\alpha_{\mathrm{bin}}=0.56$). If a higher triple fraction were assumed relative to the binary fraction, for which there is observational evidence for massive multiple systems (Moe & Di Stefano, 2017), then the importance of triple systems to the NS-NS merger rate is even larger. In particular, the ratio of the triple to the binary merger rate is $\propto\alpha_{\mathrm{tr}}/\alpha_{\mathrm{bin}}$ (see equation 10); therefore, if we assumed, e.g., $\alpha_{\mathrm{tr}}=0.5$ and $\alpha_{\mathrm{bin}}=0.4$ instead of $\alpha_{\mathrm{tr}}=0.25$ and $\alpha_{\mathrm{bin}}=0.56$, this ratio would increase by a factor of $\simeq 3$.

Our rates vary from typically several hundred to several thousand $\mathrm{Gpc^{-3}\,yr^{-1}}$, strongly depending on $\sigma_{\mathrm{k}}$. The rate inferred by GW170817 is $1540_{-1220}^{+3200}\,\mathrm{Gpc^{-3}\,yr^{-1}}$ (Abbott et al., 2017b), which falls well within our ranges (for both triples and binaries). Based on observations of short gamma-ray bursts, Coward et al. (2012) find a rate of $\sim 8-1800\,\mathrm{Gpc^{-3}\,yr^{-1}}$, Petrillo et al. (2013) of $\sim 500-1500\,\mathrm{Gpc^{-3}\,yr^{-1}}$, Siellez et al. (2014) of $\sim 92-1154\,\mathrm{Gpc^{-3}\,yr^{-1}}$, and Fong et al. (2015) of $\sim 90-1850\,\mathrm{Gpc^{-3}\,yr^{-1}}$, also consistent with our rates.

Theoretical studies based on isolated binary evolution typically find NS-NS merger rates of several tens to hundreds $\mathrm{Gpc^{-3}\,yr^{-1}}$, up to order thousand in extreme cases. For example, Belczynski et al. (2016b) found rates of $\sim 50-150\,\mathrm{Gpc^{-3}\,yr^{-1}}$, and Kruckow et al. (2018) found rates of $\sim 10-400\,\mathrm{Gpc^{-3}\,yr^{-1}}$ (the upper limit being ‘rather optimistic’). The rates are highly sensitive to the CE $\alpha$ parameter, as shown by Giacobbo & Mapelli (2018), who found rates up to several hundred $\mathrm{Gpc^{-3}\,yr^{-1}}$ and up to $\sim 10^{3}\,\mathrm{Gpc^{-3}\,yr^{-1}}$ if a high CE $\alpha$ value is assumed, and Chruslinska et al. (2018), who found similar results. Another crucial ingredient is the assumption on natal kicks (e.g., Giacobbo & Mapelli, 2019).

Generally, isolated binary evolution studies find rates that are lower than our ‘binary’ rates, as listed in Table 8. This can be understood by noting that our ‘binary’ population consists of the inner binaries of a triple population. The requirement of dynamical stability implies that the inner binary separation distribution becomes skewed towards smaller values. As can be seen in Fig. 2 (see the black dotted lines), the inner binary semimajor axis distribution of the binaries in our simulations has a median value of $\approx 10\,\,\textsc{au}$. In contrast, for a distribution flat in $\log_{10}(a_{1})$ without the restriction of dynamical stability imposed by the tertiary, the median semimajor axis is $a_{\mathrm{low}}(a_{\mathrm{up}}/a_{\mathrm{low}})^{1/2}\simeq 32\,\,\textsc{au}$ (setting $a_{\mathrm{low}}=0.02\,\,\textsc{au}$, and $a_{\mathrm{up}}=5\times 10^{4}\,\,\textsc{au}$, see Fig. 2). A more compact distribution of semimajor axes implies a higher formation rate of double NS binaries, explaining the higher base binary NS-NS merger rate in our simulations. Another way of understanding this trend is by noting that a more compact initial semimajor axis distribution is somewhat analogous to assuming a higher CE efficiency. The latter indeed leads to higher merger rates, of several hundred and up to $\sim 10^{3}$ $\mathrm{Gpc^{-3}\,yr^{-1}}$ (Giacobbo & Mapelli, 2018).

6. Conclusions

In this paper, we estimated the rates of mergers of double neutron stars (NSs) in triple systems. We took into account secular, stellar, tidal and binary evolution, and the effects of supernovae (SNe) on the orbits, starting with main-sequence (MS) stars, until the merger of two NSs. We made different assumptions on the properties of the massive triple progenitors, including different kick distributions, and different assumptions on the tertiary mass distribution. Our main conclusions are given below.

1. Contrary to previous studies of the secular and stellar evolution of triples which focused on wide inner binaries that do not interact in the absence of a tertiary star (e.g., Hamers et al., 2013; Antonini et al., 2017; Hamers & Thompson, 2019), we found that the tertiary does not significantly affect the probability of Roche-lobe overflow (RLOF) in the inner binary system. This can be attributed to the fact that we did not restrict to wide, non-interacting systems, but instead considered the entire range of orbital separations. Consequently, the inner binary interacts also in the absence of the tertiary star. Typically, about $60\%$ of systems undergo RLOF in our simulations, for both triple and binary cases (in the latter case, the same inner binary parameters are adopted, but without taking into account the effects of the tertiary star). We modeled the subsequent evolution of RLOF systems using a binary population synthesis code, neglecting the effect of the tertiary. Subsequently, the stars in the inner binary can merge as two NSs after periods of mass transfer and/or common-envelope evolution. Future work should not make the simplifying decoupling assumption, but model the system self-consistently (see Section 5.1).

2. For $\sim 10$ up to $\sim 50\%$ of systems, the inner and/or outer binary becomes unbound due to the effects of SNe (instantaneous mass loss and/or the effects of SNe kicks). Kicks are typically more important for triples, since a third star adds more possibilities for the system to become unbound, and this effect increases in importance with increasing kick dispersion. In up to $\sim 30\%$ of simulated triples, the inner binary remains bound after a SNe event whereas the tertiary star becomes unbound. The inner binaries of these systems can still potentially merge at a later time due to ‘isolated binary’ evolution. We continued the evolution of these isolated binaries using a binary population synthesis code (and considered these systems to be part of the original ‘triple’ population). The fraction of triple systems without strong interactions is small, typically a few per cent. This can be attributed to the large RLOF fraction (due to the inclusion of tighter systems), and the importance of SNe kicks, especially for NSs.

3. We considered two pathways for NS-NS mergers in triples: following binary interactions after RLOF (possibly induced by the tertiary star through Lidov-Kozai oscillations), or following binary interactions after the inner binary became unbound from the tertiary due to a SNe event. For the equivalent binaries (with the same properties as the inner binaries of the triples that we modeled), we considered only the channel of mergers following binary interactions after RLOF. Our rates (see Fig. 7 for a visual summary) vary from typically several hundred to several thousand $\mathrm{Gpc^{-3}\,yr^{-1}}$, and are within the rate estimates of LIGO based on GW170817, $1540_{-1220}^{+3200}\,\mathrm{Gpc^{-3}\,yr^{-1}}$ (Abbott et al., 2017b). We find that the rates decrease strongly with increasing SNe kick speed, $\sigma_{\mathrm{k}}$. Also, the ratio of the triple to binary NS-NS merger rate decreases with increasing $\sigma_{\mathrm{k}}$. Our ‘binary’ rates are higher compared to dedicated isolated binary evolution studies. This can be understood from the more compact inner binaries in our simulations, which are the result of the requirement of dynamical stability of the corresponding triple system.

4. Most of the NS-NS mergers in our models occur relatively early, with a delay-time distribution (DTD) peaked around several tens of Myr. Some mergers can occur at late times, of several Gyr.

5. We also considered mergers of other types of stars (see Figs. 4 and 5). In particular, we found a large fraction of NS-MS mergers in both triple and binary cases. Such mergers result in Thorne-$\mathrm{\dot{Z}}$ytkow objects (Thorne & Zytkow, 1977), and are also found in dedicated isolated binary evolution studies (e.g., Brandt & Podsiadlowski, 1995). Our results show that this channel is also possible (and relatively likely) in triples; in fact, the fractions of Thorne-$\mathrm{\dot{Z}}$ytkow objects in our simulations are very similar in the triple and binary cases, giving formation rates on the order of several thousand $\mathrm{Gpc^{-3}\,yr^{-1}}$ (see Table 8). Taking into account the occurrence rates of triples and binaries, the relative formation rate of these objects formed in our simulations ranges between $R_{\mathrm{T\dot{Z},triple}}/R_{\mathrm{T\dot{Z},binary}}\sim 0.4$ to $\sim 0.6$ (depending on the assumed kick speed and tertiary mass ratio distributions). In addition, we found BH-NS mergers at a rate of typically several hundred $\mathrm{Gpc^{-3}\,yr^{-1}}$ (see Table 8).

Acknowledgements

We thank the anonymous referee for a helpful report. A.S.H. gratefully acknowledges support from the Institute for Advanced Study, and the Martin A. and Helen Chooljian Membership. T.A.T. is supported in part by a Simons Foundation Fellowship, an IBM Einstein Fellowship from the Institute for Advanced Study, NSF grant 1313252, and Scialog Scholar grant 24216 from the Research Corporation.