# Stellar triples with chemically homogeneously evolving inner binaries

Andris Dorozsmai1, Silvia Toonen2, Alejandro Vigna-Gómez3, Selma E. de Mink3,2
and Floris Kummer2
1Institute of Gravitational Wave Astronomy and School of Physics and
Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, United
Kingdom
2 Astronomical Institute Anton Pannekoek, University of Amsterdam, Science
Park 904, 1098 XH Amsterdam, The Netherlands
3Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748
Garching, Germany E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Observations suggest that massive stellar triples are common. However, their
evolution is not yet fully understood. We investigate the evolution of
hierarchical triples in which the stars of the inner binary experience
chemically homogeneous evolution (CHE), particularly to understand the role of
the tertiary star in the formation of gravitational-wave (GW) sources. We use
the triple-star rapid population synthesis code TRES to determine the
evolution of these systems at two representative metallicities: $Z=0.005$ and
$Z=0.0005$. About half of all triples harbouring a CHE inner binary (CHE
triples) experience tertiary mass transfer (TMT) episodes, an event which is
rare for classically evolving stars. In the majority of TMT episodes, the
inner binary consists of two main-sequence stars (58-60 per cent) or two black
holes (BHs, 24-31 per cent). Additionally, we explore the role of von Zeipel-
Lidov-Kozai (ZLK) oscillations for CHE triples. ZLK oscillations can result in
eccentric stellar mergers or lead to the formation of eccentric compact
binaries in systems with initial outer pericenters smaller than $\sim$ 1200
$R_{\odot}$. Approximately $24$-$30$ per cent of CHE triples form GW sources,
and in 31 per cent of these, the tertiary star plays a significant role and
leads to configurations that are not predicted for isolated binaries. We
conclude that the evolution of CHE binaries can be affected by a close
tertiary companion, resulting in astronomical transients such as BH-BH
binaries that merge via GW emission orders of magnitude faster than their
isolated binary counterparts and tertiary-driven massive stellar mergers.

###### keywords:

gravitational waves, stars: evolution, stars: massive, stars:black holes,
binaries:close

††pagerange: Stellar triples with chemically homogeneously evolving inner
binaries–B

##  1 Introduction

An accurate and detailed understanding of the evolution of massive stars is
essential for various important open questions in astrophysics, such as
nucleosynthesis of heavy elements, the origin of supernova events, gamma-ray
bursts, and GW sources (e.g. Langer 2012). Observational evidence shows that
the fraction of stars in hierarchical triples or in higher-order multiple-
stellar systems increases with the mass of the primary star (Evans 2011; Sana
et al. 2014). In particular, Moe & Di Stefano (2017) showed that the majority
of O-type stars reside either in triple or quadruple stellar systems. This
implies that in order to understand the evolution of massive stars, and to
correctly interpret the various astrophysical phenomena related to them, we
need to consider stellar interactions in hierarchical triples.

The evolution of hierarchical triples involves a complex interplay between
three-body dynamics, stellar evolution, and stellar interactions (e.g. Toonen
et al., 2016). Three-body interactions can result in e.g. ZLK oscillations
(von Zeipel, 1910; Lidov, 1962; Kozai, 1962; Naoz, 2016), a secular effect
where the eccentricity of the inner binary can be significantly enhanced as a
result of dynamics. ZLK oscillations coupled with various dissipative
processes (e.g., tides, GWs) can shrink the orbit (e.g. Mazeh & Shaham, 1979;
Fabrycky & Tremaine, 2007; Thompson, 2011) and prompt the merger of the inner
binary (e.g. Perets & Fabrycky, 2009; Vigna-Gómez et al., 2022). These type of
mergers can result in astronomical transient events such as Type Ia supernova
(e.g. Katz & Dong, 2012; Hamers et al., 2013; Toonen et al., 2018; Swaruba
Rajamuthukumar et al., 2022) or double compact object mergers (e.g. Antognini
et al., 2014; Antonini et al., 2017; Rodriguez & Antonini, 2018; Hamers &
Thompson, 2019; Fragione & Loeb, 2019; Stegmann et al., 2022a). Furthermore,
stellar evolution can affect the orbital dynamics of the triple. For example,
radial expansion and mass loss can prompt ZLK oscillations or dynamical
instabilities (Perets & Kratter, 2012; Shappee & Thompson, 2013; Michaely &
Perets, 2014; Toonen et al., 2022; Hamers et al., 2022b).

Population synthesis studies of stellar triples show that the inner binaries
in hierarchical triples have increased stellar interactions compared to
isolated binaries (e.g. Toonen et al., 2020; Stegmann et al., 2022b; Hamers et
al., 2022a). Similarly, tertiary-driven dynamics could play an essential role
in double compact object mergers. While GW sources detected by the LIGO/Virgo
collaboration (LVC, e.g. Abbott et al., 2019b, a, 2021; The LIGO Scientific
Collaboration et al., 2021) have been studied in the context of stellar
triples, this has been done so far only in a limited parameter space. For
example, for systems in which the inner binary is wide enough such that
interaction between the two stars in the form of mass exchange can be
neglected (e.g. Silsbee & Tremaine, 2017; Antonini et al., 2017; Rodriguez &
Antonini, 2018; Fragione & Loeb, 2019; Vigna-Gómez et al., 2021; Martinez et
al., 2022), or in which the stars of the inner binary merge during the main
sequence (Stegmann et al., 2022a).There are still major uncertainties and a
need to explore and to understand the population of merging binary BHs from
hierarchical triples.

In this paper, we focus on the evolution of hierarchical triples in which the
stars of the inner binaries are chemically homogeneously evolving. CHE stars
have been discussed in the context of rapidly-rotating stars (Maeder, 1987;
Yoon & Langer, 2005; Yoon et al., 2006), which can experience enhanced mixing
during the MS stage. This mixing allows hydrogen-rich matter in the radiative
envelope to be deposited into the convective core, where it is fused to
helium. At the same time, helium is mixed throughout the star. This prevents
the build-up of a chemical gradient inside the star and the classical core-
envelope structure. As a result, the stars remain very compact over their
lifetime. CHE has been proposed to occur in very close binaries where the
tidal deformation of both stars is strong and they are forced to rotate
rapidly (de Mink et al., 2009; Song et al., 2016). More recently, CHE binaries
received renewed interest as they have been proposed as a new pathway to form
BH binaries that can merge within the age of the universe (de Mink & Mandel,
2016; Mandel & de Mink, 2016; Marchant et al., 2016; du Buisson et al., 2020;
Riley et al., 2021). Recently, Vigna-Gómez et al. (2021) studied triples with
CHE inner binaries in the context of sequential merging BH-BHs with masses
that fall in the pair-instability mass gap. Specifically, they considered
sequential mergers of hierarchical co-planar triples, a simplified approach
which neglected three-body dynamics. In this paper, we remove the constraints
of co-planarity and explore, for the first time, the evolution of massive
stellar triples with CHE inner binaries in the entire parameter space. As
isolated CHE binaries are known to be promising GW progenitors, we will mostly
focus on the role of the tertiary star in the evolution of the inner binary in
the context of GW astronomy.

This paper is structured as follows. In section 2, we introduce TRES, the
triple evolutionary code we use in this study, and the adaptations we have
made to model CHE and contact binaries. In section 3, we discuss the results
of our population synthesis in TRES and identify the most important
evolutionary channels. In section 4, we show that the initial parameters of
the tertiary star are sufficient to predict the evolutionary channel of each
system. In section 5, we use analytical and numerical methods to explore our
synthetic population of stellar triples in the context of GW sources. Finally,
we discuss the main difference between the evolution of triples with and
without CHE stars in their inner binary.

##  2 Methodology

We use TRES to simulate the evolution of our hierarchical triples (see Toonen
et al., 2016, for a detailed description of the code). TRES couples secular
dynamics of stellar triples with stellar evolution, and takes into account
additional physical processes such as stellar interactions and dissipative
processes.

TRES determines the evolution of each star by using the fitting formulae of
Hurley et al. (2000) to the stellar tracks of Pols et al. (1998) from the
rapid binary synthesis code SeBa (Portegies Zwart & Verbunt, 1996; Toonen et
al., 2012), while interactions between the stars are determined by TRES. TRES
treats three-body dynamics in the following way. For secular evolution, we
include secular three body dynamics (subscript ‘3b’) including quadrupole
(Harrington 1968) and octupole terms (Ford et al. 2004 with corrections of
Naoz et al. 2013). Regarding the additional physical processes, we take into
account: i) general relativistic effects (GR) and GW emission (subscript ‘GR’
Peters, 1964; Blaes et al., 2002), ii) tidal friction (subscript ‘TF’ Hurley
et al., 2002), iii) the effects of stellar winds under the assumptions of
fast, adiabatic wind at the mass loss rate provided by SeBa (subscript
‘wind’), iv) precession due to ZLK, GR, tides (subscript ‘tides’ Smeyers &
Willems, 2001) and intrinsic stellar rotation (subscript ‘rotate’ Fabrycky &
Tremaine, 2007), and v) the change in the stellar rotation due to stellar
evolution based on spin angular momentum conservation (subscript ‘I’). This
gives rise to a set of first-order ordinary differential equations, that are
solved numerically. These equations are:

$\displaystyle\left\\{\begin{array}[]{l c l}\dot{a}_{\rm in}&=&\dot{a}_{\rm
in,GR}+\dot{a}_{\rm in,TF}+\dot{a}_{\rm in,wind}\\\ \dot{a}_{\rm
out}&=&\dot{a}_{\rm out,GR}+\dot{a}_{\rm out,TF}+\dot{a}_{\rm out,wind}\\\
\dot{e}_{\rm in}&=&\dot{e}_{\rm in,3b}+\dot{e}_{\rm in,GR}+\dot{e}_{\rm
in,TF}\\\ \dot{e}_{\rm out}&=&\dot{e}_{\rm out,3b}+\dot{e}_{\rm
out,GR}+\dot{e}_{\rm out,TF}\\\ \dot{g}_{\rm in}&=&\dot{g}_{\rm
in,3b}+\dot{g}_{\rm in,GR}+\dot{g}_{\rm in,tides}+\dot{g}_{\rm in,rotate}\\\
\dot{g}_{\rm out}&=&\dot{g}_{\rm out,3b}+\dot{g}_{\rm out,GR}+\dot{g}_{\rm
out,tides}+\\\ &&\dot{g}_{\rm out,rotate}\\\ \dot{h}_{\rm in}&=&\dot{h}_{\rm
in,3b}\\\ \dot{\theta}&=&\frac{-1}{J_{\rm b,in}J_{\rm b,out}}[\dot{J}_{\rm
b,in}(J_{\rm b,in}+J_{\rm b,out}\theta)+\\\ &&\dot{J}_{\rm b,out}(J_{\rm
b,out}+J_{\rm b,in}\theta)]\\\ \dot{\Omega}_{1}&=&\dot{\Omega}_{\rm
1,TF}+\dot{\Omega}_{\rm 1,I}+\dot{\Omega}_{\rm 1,wind}\\\
\dot{\Omega}_{2}&=&\dot{\Omega}_{\rm 2,TF}+\dot{\Omega}_{\rm
2,I}+\dot{\Omega}_{\rm 2,wind}\\\ \dot{\Omega}_{3}&=&\dot{\Omega}_{\rm
3,TF}+\dot{\Omega}_{\rm 3,I}+\dot{\Omega}_{\rm 3,wind}\end{array}\right.$ (14)

where $a$, $e$, $g$, $h$ and $J_{b}$ represent the semimajor axis,
eccentricity, argument of pericenter, line of ascending nodes, and the orbital
angular momentum for the inner (subscript ‘in’) and outer (subscript ‘out’)
orbit. The dot represents the time derivatives. Lastly $\theta\equiv\cos(i)$,
where $i$ is the mutual inclination between the inner and outer orbit, and
$\Omega_{1},\Omega_{2},\Omega_{3}$ the spin frequency of the primary,
secondary and tertiary star respectively. Per definition the primary and
secondary stars are the stars in the inner binary, with the primary star
initially more massive than the secondary star, and the tertiary star orbits
the inner binary.

We highlight three aspects of the orbital evolution of hierarchical triples
that is particularly relevant for the systems we study in this paper. Firstly,
if the apsidal precession of the inner binary due to short range forces, such
as tides ($\dot{g}_{\rm in,tides}$) and GR effects ($\dot{g}_{\rm in,GR}$)
occurs on a much shorter timescale than the precession due to three-body
dynamics ($\dot{g}_{\rm in,3b}$), ZLK oscillations will be quenched (see e.g.
Holman et al., 1997; Eggleton & Kiseleva-Eggleton, 2001; Blaes et al., 2002;
Fabrycky & Tremaine, 2007; Thompson, 2011; Dong et al., 2014; Liu et al.,
2015; Petrovich, 2015; Anderson et al., 2017). The timescale of ZLK
oscilations can be approximated as (e.g. Innanen et al., 1997; Holman et al.,
1997; Kinoshita & Nakai, 1999):

$t_{\rm ZLK}=\left(\frac{M_{1}+M_{2}}{GM_{\rm
out}^{2}}\right)^{1/2}\left(\frac{a_{\rm out}}{a_{\rm
in}^{1/2}}\right)^{3}(1-e_{\rm out}^{2})^{3/2}.$ (15)

The timescale related to the apsidal precession due to tides are (e.g. Smeyers
& Willems, 2001; Liu et al., 2015):

$t_{\rm tides}=\left(\frac{M_{1}}{15k_{\rm am}\mu_{\rm
in}^{1/2}M_{2}}\right)\left(\frac{a_{\rm
in}^{11/2}}{R_{1}^{5}}\right)\left(\frac{(1-e_{\rm
in}^{2})^{5}}{1+\frac{3}{2}e_{\rm in}^{2}+\frac{1}{8}e_{\rm in}^{4}}\right),$
(16)

where $k_{\rm am}$ the apsidal motion constant, which we assume to be 0.0144
for MS and helium stars, $\mu_{\rm in}=G(M_{1}+M_{2})$, i.e. the standard
gravitational parameter for the inner binary and $R_{1}$ is the radius of the
inner star. The timescale related to precession due to general relativistic
effects is (e.g. Misner et al., 1973; Blaes et al., 2002; Miller & Hamilton,
2002):

$t_{\rm GR}=\frac{c^{2}}{3\mu_{\rm in}^{3/2}}a_{\rm in}^{5/2}(1-e_{\rm
in}^{2}).$ (17)

If $t_{\rm ZLK}\gg\rm{min}(t_{\rm GR},t_{\rm tides})$, then three-body
dynamics are suppressed. If the timescales are comparable, then the maximum
eccentricity induced by the ZLK oscillations is diminished. In principle,
rotation-induced oblateness in the inner binary also induces apsidal
precession ($g_{\rm in,rot}$, see e.g. Fabrycky & Tremaine, 2007). However, as
long as the rotational period of the inner stars is not shorter than the
orbital period (which is true for all systems considered here), $\dot{g}_{\rm
tides}\gg\dot{g}_{\rm rot}$ and therefore precession due to stellar rotation
does not play a role in suppressing three-body dynamics (Liu et al., 2015).

Secondly, octupole terms in the three-body dynamics are typically negligible
for CHE triples, as the mass ratio of the inner binary is always very close to
one. Finally, we estimate the time it takes for the inner binary to merge due
to GWs following Peters (1964), if the tertiary is dynamically decoupled from
the inner binary. If ZLK oscillations are still relevant during the inspiral
phase, we follow the approximation of Miller & Hamilton (2002):

$t_{\rm GW}\approx t_{\rm GW,Peters}(a_{\rm in},e_{\rm in,max})(1-e_{\rm
in,max})^{-1/2},$ (18)

where $t_{\rm GW}$ is the time required for the merger, $t_{\rm GW,Peters}$ is
the time to merger based on the relation of Peters (1964), $e_{\rm in,max}$ is
the maximum eccentricity reached during ZLK oscillations and $a_{\rm in}$ is
the initial inner semimajor axis. The approximation in equation 18 is based on
Wen (2003) and it neglects the effects of precession due GR. When the latter
is taken into account, Thompson (2011) finds that equation 18 underestimates
the actual merger timescale typically by a factor of 2-3.

###  2.1 Modelling of chemically homogeneous evolution

We follow Riley et al. (2021) in order to incorporate CHE stars in `TRES`.
That means that we assume a star evolves chemically homogeneously, if the
angular frequency of the spin of the star is above a certain critical value,
i.e. $\omega_{\rm star}>\omega_{\rm CHE,crit}$. Riley et al. (2021) provides a
fit to this critical value based on `MESA` (Paxton et al., 2011) models at
different masses and metallicities. In order to determine whether a star
evolves chemically homogeneously, we check whether our simulated star is
spinning above $\omega_{\rm CHE,crit}$ at every timestep. If a star meets this
criteria we do not evolve its radius during that timestep. We assume that the
star by the end of core hydrogen burning forms a helium star with a mass
$M_{\rm He,ZAMS}=M_{\rm TAMS}$, where $M_{\rm He,ZAMS}$ is the initial mass of
the helium star and $M_{\rm TAMS}$ is the terminal age main sequence mass of
the star. With these assumptions, CHE stars experience an instantaneous drop
in radii at the end of their MS phase (compare main sequence stellar evolution
with helium star evolution in Hurley et al., 2000). This is a simplification
of the results of detailed simulations of CHE stars, where the latter suggests
a gradual contraction of the radius during the MS (e.g. Maeder, 1987). If a
CHE star loses angular momentum (e.g. due to stellar winds), its rotational
frequency decreases. If the frequency reduces to below the critical value, we
assume the evolution of the star transitions back to the classical non-CHE
case.

For simplicity, we only consider systems in which the stars of the inner
binary are CHE from zero-age main sequence (ZAMS). Stars that do not evolve
chemically homogeneously from ZAMS could, in theory, become CHE if they
attained a sufficiently high-spin frequency before a significant chemical
gradient is built up in their interior. This can be achieved for example, if a
star is spun up by accretion during a mass transfer event (e.g., Cantiello et
al., 2007; Ghodla et al., 2022). We neglect such systems in this study.

###  2.2 Contact binaries

We follow the implementation of Riley et al. (2021) for modelling contact
binaries (which is based on the models of Marchant et al., 2016). We assume
that contact binaries, i.e. binaries in which both stars fill their Roche-
lobes, can maintain co-rotation and consequently survive the contact phase
without merging as long as neither of the stars fill the outer Lagrangian
points (L2 and L3). For contact binaries, Marchant et al. (2016) finds that
mass is transferred between the two stars back and forth until they reach an
equal mass ratio. We follow Marchant et al. (2016) and approximate the L2
point as

$\frac{R_{\rm{L2,2}}-R_{\rm{RL,2}}}{R_{\rm{RL,2}}}=0.299\,\rm{tan}^{-1}(1.84q^{0.397}),$
(19)

where $R_{\rm{RL,2}}$ is the Roche-lobe radius of the secondary star, which we
approximate following Eggleton (1983).

If the stars in the inner binary are in contact but without filling their L2
points, we assume that the masses of the binary equalise via a fully
conservative mass transfer phase. We follow Riley et al. (2021) and assume
this mass equalisation occurs instantaneously and readjust the orbit of the
inner binary as (see, e.g. Soberman et al., 1997):

$\frac{a_{\rm fin}}{a_{\rm init}}=\left(\frac{M_{\rm 1,init}M_{\rm
2,init}}{M_{\rm 1,fin}M_{\rm 2,fin}}\right)^{2},$ (20)

where $a_{\rm init}$, $a_{\rm fin}$ are the initial and the final orbital
separation and $M_{\rm 1,init}$, $M_{\rm 2,init}$ are the initial masses of
the primary and the secondary, respectively. The final masses are $M_{\rm
1,fin}$ = $M_{\rm 2,fin}$ = $1/2\cdot(M_{\rm 1,init}+M_{\rm 2,init})$ by
definition. The assumption of mass equalisation for contact binaries results
in the prediction of the CHE channel leading to mostly equal-mass binary BH
mergers (e.g. Marchant et al., 2016).

###  2.3 Stellar winds

The mass loss rates of stellar winds and their effects on the evolution of the
star are determined by `SeBa` (Hurley et al., 2000; Toonen et al., 2012),
while the effects on the orbit of the triple are determined by TRES (equation
14). In this study, we use the same implementation of stellar winds for
massive stars as in Dorozsmai & Toonen (2022) with one difference; the mass
loss rates of helium stars and giants are calculated according to the
empirical formula of Hamann et al. (1995) instead of Sander & Vink (2020).

For reference, we summarise the mass loss rates prescriptions used in this
study. For MS stars, we follow Vink et al. (2001), if $T_{\rm{eff}}\leq 50$ kK
and Nieuwenhuijzen & de Jager (1990), if $T_{\rm{eff}}>50$ kK. For evolved
stars crossing the Hertzsprung gap or core helium burning (CHeB) stars, we
follow Vink et al. (2001), if $T_{\rm{eff}}\geq 8$ kK or the maximum between
Nieuwenhuijzen & de Jager (1990) and Reimers (1975), if $T_{\rm{eff}}<8$ kK.
For evolved stars beyond the Humphreys-Davidson limit, we assume
$\dot{M}_{\rm{LBV}}=1.5\cdot 10^{-4}\ M_{\odot}\rm{yr^{-1}}$ (Belczynski et
al., 2010). For Asymptotic Giant Brach stars and double shell burning
supergiants, we calculate the maximum between Nieuwenhuijzen & de Jager
(1990), Reimers (1975) and Vassiliadis & Wood (1993). Finally, for helium
stars we follow the empirical form from Hamann et al. (1995) in the form
$\dot{M}_{\rm{WR}}=0.5\cdot
10^{-13}\cdot\left(\frac{L}{L_{\odot}}\right)^{1.5}\left(\frac{Z}{Z_{\odot}}\right)^{0.86}$
with a clumping factor of $\eta=0.5$ from Hamann & Koesterke (1998) and a
metallicity scaling of $\dot{M}_{WR}\sim Z^{0.86}$ (Vink & de Koter 2005).

In order to compute the change in the orbit due to stellar winds, we assume
stellar winds are spherically symmetric and fast compared to the orbital
velocity; additionally, we neglect wind accretion by the companions. In that
case the inner and the outer orbit of the triple widens as

$\dot{a}_{\rm in,wind}=\left(\frac{a_{\rm final}}{a_{\rm init}}\right)_{\rm
in}=\frac{M_{\rm 1,init}+M_{\rm 2,init}}{M_{\rm 1,final}+M_{\rm 2,final}},$
(21)

and

$\dot{a}_{\rm out,wind}=\left(\frac{a_{\rm final}}{a_{\rm init}}\right)_{\rm
out}=\frac{M_{\rm 1,init}+M_{\rm 2,init}+M_{\rm 3,init}}{M_{\rm
1,final}+M_{\rm 2,final}+M_{\rm 3,final}},$ (22)

where subscripts ‘init’ and ‘final’ refer to properties before and after the
stellar winds carried mass away from the stars in a given timestep. We assume
that the eccentricity remains unchanged by stellar winds (Huang, 1956, 1963).

We neglect stellar wind accretion by the other stars in the triple system (see
e.g. Bondi & Hoyle, 1944). Neglecting accretion is justified for line-driven
winds due to their large terminal velocities (see e.g. Vink et al., 2001). The
assumptions of a fast and spherically symmetric wind might not always be valid
(e.g. Brookshaw & Tavani, 1993), and rapidly rotating stars might not have
fully symmetric outflows (Georgy et al., 2011). Particularly, stellar winds in
certain binary-configurations might even lead to orbital shrinking (Schrøder
et al., 2021).

###  2.4 Remnant formation

The mass of the compact object remnant is computed based on the delayed
supernova model from Fryer et al. (2012). This prescription gives the mass of
the stellar remnant as a function of CO core mass, where the latter is
determined in SeBa based on the fits of Hurley et al. (2000). The natal kick
velocity for BHs is calculated as

$v_{BH}=(1-f_{\rm b})\left(\frac{M_{NS}}{M_{BH}}\right)v_{\rm kick},$ (23)

where $f_{\rm b}$ is the fallback fraction (Fryer et al., 2012), $M_{\rm NS}$
is the canonical neutron star mass ($M_{\rm NS}=1.4M_{\odot}$) and $v_{\rm
kick}$ is a random kick velocity drawn from the distribution inferred by
Verbunt et al. (2017) from proper motion measurements of pulsars. We determine
the change in the inner and outer orbit due to the core collapse of any of the
stars in the triple system based on the formalism developed in Pijloo et al.
(2012).

Models of Fryer et al. (2012) predict that the most massive stars collapse
directly (typically $M_{\rm ZAMS}\gtrsim 40\,M_{\odot}$), without any ejecta,
and the only mass loss during the remnant formation is due to neutrino losses,
which is assumed to be 10 per cent of the pre-core-collapse mass of the star.
Additionally, we assume that the neutrino emission is spherically symmetrical
and do not impart natal kick onto the BH. In this case, the orbit is only
changed due to the instantaneous mass loss (e.g. via Blaauw kick, see Blaauw,
1961). We note that, if the pre-core-collapse orbit is circular, a Blauuw kick
due to neutrino losses does not lead to a significant change in the inner
orbital elements. However, this is no longer the case for eccentric pre-core-
collapse orbits. In particular, if the core collapse occurs near the
pericenter, the orbit can become significantly wider (e.g. Hills, 1983).

By the onset of core-oxygen burning, the core temperatures of the most massive
stars can reach above $T_{\rm core}\sim 3\times 10^{9}\,K$. Under these
conditions, the emitted gamma-ray photons in the core are energetic enough to
form electron-positron pairs. This leads pair-instability (see e.g. Fowler &
Hoyle 1964, Rakavy & Shaviv 1967, Barkat et al. 1967, Fraley 1968). Depending
on the mass of the star, this instability can result in a pulsation pair
instability supernova, in which the star experiences a series of pulsations
leading to severe mass loss (i.e. or PPISN, see e.g. Yoshida et al., 2016;
Marchant et al., 2016; Woosley, 2017; Renzo et al., 2020), or pair instability
supernova, in which the star is completely disrupted and no remnant is formed
(PISN, see e.g. Yoshida et al., 2016; Marchant et al., 2016; Woosley, 2017;
Renzo et al., 2020). For the treatment of pair-instability in massive stars,
we follow Stevenson et al. (2019). If the mass of the helium star pre-core-
collapse is $M_{\rm HE,pre-SN}\geq 35\,M_{\odot}$, the star is assumed to
undergo PPISN, and its remnant mass is determined by the fitting formula of
Stevenson et al. (2019), based on the detailed stellar simulations of Marchant
et al. (2019). If $60\leq M_{\rm HE,pre-SN}\leq 130\,M_{\odot}$, we assume the
star undergoes PISN, and leaves no remnant behind. In principle, if $M_{\rm
HE,pre-SN}\geq 130\,M_{\odot}$, photo-disintegration prevents the pair
instability supernova to occur and the star collapses directly into a BH (Bond
et al. 1982, Woosley & Weaver 1982, Heger & Woosley 2002, du Buisson et al.
2020), however this does not occur for any of our simulated systems.

###  2.5 Tertiary mass transfer (TMT) episodes

If the tertiary star fills its Roche-lobe, it will transfer mass to the inner
binary. There have been some efforts to study and model this process (de Vries
et al., 2014; Leigh et al., 2020; Comerford & Izzard, 2020; Glanz & Perets,
2021a; Soker & Bear, 2021; Moreno Méndez et al., 2022), but this complex
scenario remains to be fully understood.

In order to calculate the Roche-lobe of the tertiary star, we assume the inner
binary can be approximated as a point mass and estimate the Roche radius with
the fitting formula of Eggleton (1983). This assumption is valid in the regime
where the orbital separation of the outer star is much larger than that of the
inner binary (e.g. $a_{\rm{out}}\gg a_{\rm{in}}$). TRES determines the
stability of TMT based on extrapolating typical methods from binary star
evolution, i.e. by using critical mass ratios (see e.g. Toonen et al., 2016).
This parameter is defined as $q_{\rm crit}=M_{\rm donor}/M_{\rm accetor}$,
i.e. the ratio of the mass of the donor and the mass of the accretor star at
the onset of the mass transfer episode. The mass transfer phase is assumed to
be dynamically unstable, if the mass ratio of the system is above the critical
mass ratio, i.e. $q>q_{\rm crit}$. We obtain $q_{\rm crit}$ for each stellar
evolutionary stage from Hurley et al. (2002) and Claeys et al. (2014). We
quote these values for the two most common donor types in our simulations (a
complete description of our assumptions about $q_{\rm crit}$ can be found in
Toonen et al., 2016). These are $q_{\rm crit}=3$ and $q_{\rm
crit}=(1.37+2[M_{\rm donor,core}/M_{\rm donor}]^{5})/2.13$ for Hertzsprung gap
stars (i.e. hydrogen shell burning stars which have not regained thermal
equilibrium yet) and core helium burning (CHeB) stars, respectively. The term
in the squared bracket is the core mass to total mass ratio of the donor. If
this equals to $\sim 0.45$ \- $0.65$, which is fairly typical for massive CHeB
stars (Dorozsmai & Toonen, 2022), then $q_{\rm crit}\approx 0.7$-$0.75$. This
reflects the assumption made by Hurley et al. (2000), CHeB stars tend to have
deep convective envelopes (cf. Klencki et al., 2020), and are therefore more
likely to experience unstable mass transfer episodes (e.g. Hjellming &
Webbink, 1987).

Stable TMT could be accompanied with the formation of a circumbinary disc or
it could occur in a ballistic accretion fashion. These two types could lead to
significantly different evolution of the inner orbit (de Vries et al., 2014).
We assume that TMT occurs via ballistic accretion, if $a_{\rm in}(1+e_{\rm
in})\geq R_{\rm cd}$ at the onset of the TMT phase, where $R_{\rm cd}$ is
(i.e. adapting the fitting formulas for mass transferring binaries of Lubow &
Shu 1975 and Ulrich & Burger 1976 to triples):

$R_{\rm cd}=0.0425\,a_{\rm out}(1-e_{\rm out})\left[\frac{1}{q_{\rm
out}}\left(1+\frac{1}{q_{\rm out}}\right)\right]^{1/4}.$ (24)

####  2.5.1 TMT: Evolution of the inner orbit

If the tertiary star fills its Roche-lobe, TRES stops the simulation of the
system. However, when discussing potential GW progenitors (Section 5), we
determine the orbital evolution due to TMT by applying simplified assumptions,
if the mass transfer episode is dynamically stable. In this subsection we
describe our assumptions about the evolution of the inner orbit during a
stable phase of TMT, while in subsection 2.5.2 we discuss the evolution of the
outer orbit.

We distinguish three particular TMT configurations cases, based on the
evolutionary stage of the inner binary and on whether or not the transferred
mass forms a circumbinary disc around the inner binary:

  1. 1.

an inner binary with compact objects and with ballistic accretion,

  2. 2.

an inner binary with compact objects and with a circumbinary disc,

  3. 3.

a non-compact inner binary.

(i) An inner binary with compact objects and with ballistic accretion.
Hydrodynamical simulations of de Vries et al. (2014) showed that in case of a
TMT episode with ballistic accretion, the transferred mass eventually engulfs
the inner binary and exerts friction on it. This leads to a scenario that
could be considered similar to the common-envelope evolution of binaries (e.g.
Paczynski, 1976; Ivanova et al., 2013), since in both cases drag forces
exerted by a gaseous medium supplied from the donor star lead to the orbital
shrinking of the binary. Inspired by this similarity, de Vries et al. (2014)
applied a modified version of $\alpha$-formalism (originally developed for
common-envelope evolution, see e.g. Tutukov & Yungelson, 1979; de Kool et al.,
1987; Dewi & Tauris, 2000) to model the inner binary evolution of triples
experiencing TMT (see also Hamers et al., 2021). For the configuration case
(i), we take the same approach.

Below we explain how the post-mass-transfer inner orbit is determined based on
this formalism in detail. $\Delta M_{\rm trnsf}$ is the mass that is
transferred from the tertiary in a timestep $\Delta t$. When $\Delta M_{\rm
trnsf}$ ends up encompassing the inner binary, it has binding energy of
$E_{\rm bind}$. As the inner orbit is shrinking due to the friction during the
TMT episode, the orbital energy of the inner binary changes by $\Delta E_{\rm
orb}$. We assume that a fraction ($\alpha_{\rm TMT}$) of $\Delta E_{\rm orb}$
is used to unbind $\Delta M_{\rm trnsf}$. We can write an equation expressing
the energy balance as:

$\alpha_{\rm TMT}\Delta E_{\rm orb}=E_{\rm bind},$ (25)

with

$\Delta E_{\rm orb}=\frac{GM_{1}M_{2}}{2a_{\rm
in,fin}}-\frac{G\left(M_{1}+\Delta M_{\rm trnsf}/2\right)\left(M_{2}+\Delta
M_{\rm trnsf}/2\right)}{2a_{\rm in,init}},$ (26)

and

$E_{\rm bind}=\frac{-G(M_{1}+M_{2})\Delta M_{\rm trnsf}}{\lambda_{\rm
TMT}a_{\rm init}},$ (27)

where $\lambda_{\rm TMT}$ is a parameter related to the structure of $\Delta
M_{\rm trnsf}$, parameterising its binding energy, $a_{\rm in,init}$ is the
initial orbital separation before $\Delta M_{\rm trnsf}$ is transferred to the
inner binary and $a_{\rm in,fin}$ is the final orbital separation after
$\Delta M_{\rm trnsf}$ is expelled from the inner binary. We assume that the
total mass transferred to the inner binary throughout the entire TMT episode
equals to the mass of the hydrogen envelope of the tertiary $M_{\rm out,env}$
(but see Laplace et al. 2020). Then assuming a constant $\alpha_{\rm TMT}$ and
$\lambda_{\rm TMT}$, the orbit changes due to the entire TMT episode as:

$\frac{a_{\rm in,fin}}{a_{\rm
in,init}}=\frac{M_{1}M_{2}}{\frac{2(M_{1}+M_{2})M_{\rm out,env}}{\alpha_{\rm
TMT}\lambda_{\rm TMT}}+\left(M_{1}+\frac{M_{\rm
out,env}}{2}\right)\left(M_{2}+\frac{M_{\rm out,env}}{2}\right)}.$ (28)

As both $\alpha_{\rm TMT}$ and $\lambda_{\rm TMT}$ are unknown, we combine
them and try three different values: $\alpha_{\rm TMT}\lambda_{\rm
TMT}=0.05,\,0.5,\,5$. Here $\alpha_{\rm TMT}\lambda_{\rm TMT}=5$ is the
fiducial value used in Hamers et al. (2021), which is in a good agreement with
the hydrodynamical simulations of de Vries et al. (2014), in which the inner
stars are on the MS during the TMT episode. We note that we neglect the
possibility of TMT episode with ballistic accretion transitioning to a TMT
episode with a circumbinary disc.

Additionally, for configuration type (i), we assume that the inner binaries
circularise as a result of the mass transfer phase (as $a_{\rm in,new}=a_{\rm
in}(1-e_{\rm in})$). We note that this assumption might not be correct for
highly eccentric inner binaries. For example, Glanz & Perets (2021b) showed
that binaries at the onset of common-envelope events with $e\gtrsim 0.95$
might retain eccentricities as high as $e\sim 0.2$.

(ii) An inner binary with compact objects with circumbinary disc. If a
circumbinary disc is formed during a mass transfer phase towards an inner BH-
BH binary we assume that the orbit of the inner binary remains unchanged. The
actual physics underlying such a process are very complex (see Lai & Muñoz,
2022, for a review on circumbinary accretion from gaseous medium). The
circumbinary disc may exert a torque on the inner binary and extract angular
momentum from it, while the accreted matter can transfer angular momentum onto
the inner binary. Furthermore, the circumbinary disc and the inner binary
could be tidally distorted by the tertiary star. It is commonly assumed that
circumbinary accretion of a BH-BH binary from a gaseous medium leads to the
shrinking of its orbit due to the torques exerted by the circumbinary disc and
due to dynamical friction of the gas (e.g. Bartos et al. 2017; Stone et al.
2017; Antoni et al. 2019; Tiede et al. 2020; Duffell et al. 2020; McKernan et
al. 2020; Rozner & Perets 2022). However, a consesus regarding this physical
process is still missing, with some hydrodynamical simulations suggesting that
accretion from circumbinary disc could even lead to to orbital widening
instead of orbital decay (e.g. Muñoz et al., 2019; Moody et al., 2019).

(iii) A non-compact inner binary. If the mass transfer occurs with a MS-MS
accretor, we assume that this results in the merger of the inner binary. We
make this assumption because these binaries have very short periods and a
sizeable fraction of them are in contact and most likely they would expand due
to TMT, overfilling their L2 point, which would lead to merger (see later
subsection 5). As we discuss in in subsection 5, we do not consider GW sources
from those triple systems, in which the TMT occurs towards a binary with
evolved (i.e. non-MS), non-compact stars.

We do not model unstable phases of TMT (as we will show later, they are very
rare among the systems we discuss in this paper) . We note, however, that
during this type of mass transfer episode, the outer orbital separation is
predicted to rapidly decrease due to the common-envelope-like evolution in
triple system; this could result in a regime where the secular approximation
from the triple is no longer valid (e.g. Glanz & Perets, 2021a; Comerford &
Izzard, 2020; Soker & Bear, 2021).

####  2.5.2 TMT: Evolution of the outer orbit

When determining the evolution of the outer orbit due to a stable phase of
TMT, we apply the same method for all accretor types, irrespective of whether
a circumbinary disc is formed. We calculate the evolution of the outer orbit
during the TMT phase, based on the following relation:

$\frac{\dot{a}_{\rm out}}{a_{\rm
out}}=-2\frac{\dot{M}_{3}}{M_{3}}\left[1-\beta\frac{M_{3}}{M_{1}+M_{2}}-(1-\beta)\left(\gamma+\frac{1}{2}\right)\frac{M_{3}}{M_{\rm
tot}}\right],$ (29)

where $\beta$ is the fraction of mass accreted by the inner binary, $\gamma$
is the specific angular momentum lost from the system as a fraction of the
specific angular momentum of the triple and $\dot{M}_{3}$ is the mass transfer
rate from the tertiary star. Equation 29 can be derived from angular momentum
arguments. It is an adaptation of the relation describing the orbital
evolution of a circular, mass transferring binary comprised of point particles
(see e.g. Soberman et al., 1997), applied to a triple experiencing a TMT
episode. This adaptation is valid, if the tertiary star is sufficiently far
away from the inner binary, such that the inner binary can be treated as a
point particle with a mass of $M_{1}+M_{2}$. We assume that eventually all the
transferred mass is isotropically expelled from the triple ($\beta=0$), from
near the inner binary. This expelled matter thus carries away a specific
angular momentum that is equal to that of the inner binary
($\gamma=M_{3}/(M_{1}+M_{2})$, see also Hamers et al. 2021, for a similar
approach). In this case equation 29 can be expressed as

$\frac{a_{\rm out,fin}}{a_{\rm out,init}}=\frac{M_{\rm tot,init}}{M_{\rm
tot,fin}}\left(\frac{M_{\rm 3,init}}{M_{\rm
3,fin}}\right)^{2}\exp\left(2\frac{M_{\rm 3,fin}-M_{\rm
3,init}}{M_{1}+M_{2}}\right).$ (30)

In case of BH-BH inner accretors, these assumptions might be valid, as the
accretion rate of BHs might be capped by the Eddington-limit, and most of the
mass could indeed be expelled from the system, for example in the form of a
jet (e.g. King et al. 2000; van den Heuvel et al. 2017) . On the other hand,
MS stars are likely to accrete more efficiently, and therefore $\beta=0$ might
no longer be a good approximation.

###  2.6 Initial conditions

We sample $10^{5}$ triples at two representative (moderate and low)
metallicities: Z = 0.005 and Z = 0.0005. We simulate each hierarchical triples
from ZAMS. After drawing the parameters for a given triple system, we further
check, if it is dynamically stable (based on the criteria of Mardling &
Aarseth, 2001) or if the stars in the inner binary are CHE at ZAMS. If any of
the two criteria are not met, we do not evolve the triple system and only take
it into account for the normalisation of event rate calculations. We terminate
the simulation of a triple system when either a Hubble time (assumed to be
13.5 Gyrs) has passed, or when the tertiary star fills its Roche lobe, a
merger occurs, a dynamical instability occurs or if any of the stars becomes
unbound from the triple. We also stop the simulation, if any of the stars in
the inner binary transitions back from CHE to classical evolution. That is, we
only consider triples in which the stars of the inner binary chemically
homogeneously evolve throughout their entire MS lifetimes. We refer to this
population as CHE triple population.

In this study, we motivate the choice of the initial distributions of the
parameters of the inner binaries based on recent surveys of massive binaries
(e.g. Sana et al., 2012; Kobulnicky et al., 2014). In such surveys, a possible
tertiary companion is not always unequivocally identified and therefore it is
not clear whether the inferred distributions also hold for triples or only for
isolated binaries.

We assume the ZAMS mass of the primary star ($M_{\rm 1,ZAMS}$) follows the
power-law mass distribution of Kroupa (2001), i.e. $N$$\sim$$M_{ZAMS}^{-2.3}$
for $M_{\rm ZAMS}\geq 0.5\,M_{\odot}$ and $N$$\sim$$M_{ZAMS}^{-1.3}$ for
$M_{\rm ZAMS}<0.5\,M_{\odot}$. We sample $M_{\rm 1,ZAMS}$ from a mass range of
20-100$\,M_{\odot}$. The lower limit approximately coincides with the lowest
initial mass at which CHE is still possible in a tidally locked binary (e.g.
Riley et al., 2021), while the upper limit is roughly the maximum mass at
which the stellar tracks used in TRES are still reasonably accurate. We assume
a flat inner mass-ratio (i.e. $q_{\rm in,ZAMS}=M_{\rm 2,ZAMS}/M_{\rm 1,ZAMS}$)
distribution, which is in broad agreement with Sana et al. (2012). We restrict
the range of $q_{\rm in,ZAMS}$ to 0.7-1 given that inner binaries in which
both of the stars are chemically homogeneously evolving and have $q_{\rm
in}\leq 0.7$ would merge early during the MS (where we found the lower limit
of 0.7 from our simulations). We sample the inner semimajor axis from a log-
uniform distribution (Öpik 1924; and in broad agreement with Sana et al. 2012)
in the range of 16 to 40 $R_{\odot}$. We assume that the inner binaries are
tidally locked at ZAMS. This has three implications: i) the inner binaries
have circular orbits, ii) their rotational angular frequency is synchronised
with the orbital angular frequency, and iii) the spins of the stars are
aligned with the orbital angular momentum vector.

We draw the properties of the outer binary from the same distributions that we
assume for the inner binaries, with the exception of outer eccentricities.
Observations of hierarchical multiple systems of galactic solar-type stars
support the assumption that the distributions of the initial parameters of the
inner and the outer binaries are the same (Tokovinin, 2014; Tokovinin et al.,
2006). We sample the outer semimajor axis from a loguniform distribution in
the range of 100 to $10^{5}$ $R_{\odot}$. We assume that the distribution of
the outer mass ratio (i.e. $q_{\rm out,ZAMS}=M_{\rm out,ZAMS}/(M_{\rm
1,ZAMS}+M_{\rm 2,ZAMS}$)) is flat on a range of 0.1 to 1, furthermore the mass
of the tertiary is restricted to a range of 5-100$\,M_{\odot}$. We assume non-
spinning tertiary stars. The eccentricities of the outer orbit are drawn from
a thermal distribution (e.g Heggie, 1975). The mutual inclination between the
inner and outer orbit is assumed to be uniform in $\rm{cos}(i_{\rm{ZAMS}})$,
where $\rm{i}_{\rm ZAMS}$ is the initial inclination. The initial argument of
the pericenter is assumed to be uniformly distributed between $-\pi$ and
$\pi$.

In Section 5, we compare our CHE triple population to a CHE isolated binary
population. To this end, we also perform population synthesis of isolated
binaries with CHE stars. We sample $10^{5}$ isolated binaries at Z = 0.005 and
Z = 0.0005 and evolve them with TRES. We sample from the same initial
distributions that we assumed for the inner binaries of our triple population.
Similarly to the triple population, we discard systems that are not CHE at
ZAMS and stop the simulation, if a Hubble time has passed, or if any of the
stars in the binary transitions from CHE to classical evolution. We only
analyse binaries, in which the stars remain CHE throughout their entire MS
lifetime (hereafter CHE binaries).

Throughout the paper, we estimate birth rate and merger rate densities of
different evolutionary channels (discussed in detail in appendix, section B).
In order to determine each of these quantities, one must know how common
single and multiple stellar systems are. We assume two different stellar
populations, with different binary and triple fractions. In the first, we
assume that about 73 per cent of massive stars are found in triples (with
multiplicity fractions of $f_{\rm single}=0.06$, $f_{\rm binary}=0.21$,
$f_{\rm triple}=0.73$, see e.g. Moe & Di Stefano, 2017)111Moe & Di Stefano
(2017) finds that 73 per cent of O stars are either in triples or quadrupoles.
Therefore $f_{\rm triple}=0.73$ should be considered as a rough upper limit.
We also note that Tokovinin et al. (2006) finds that there is a strong
correlation between the inner period and the triple multiplicity; among solar
type stellar systems, 96 per cent of the spectroscopic binaries with periods
less than 3 days has a tertiary companion. Therefore CHE triples, which have
also inner binaries with periods of few days, could have exceptionally high
triple fractions too., whereas in the second test population, we assume there
are no triples and about 70 per cent of massive stars are in binaries ($f_{\rm
single}=0.3$, $f_{\rm binary}=0.7$, $f_{\rm triple}=0$, see e.g. Sana et al.,
2012) 222Strictly speaking, Sana et al. (2012) did not make any statements
about triple fractions, but they found that 70 per cent of massive stars have
companions that are sufficiently close such that mass exchange will occur some
time in their evolution..

##  3 Results of population synthesis simulations

Table 1: An overview of our sampled triples based on the evolutionary type of
the inner binary. For the definitions of the different categories, see text in
section 3.1.

| Z = 0.005 | Z = 0.0005 | Combined
---|---|---|---
| % of simulation | % of CHE at ZAMS | % of simulation | % of CHE at ZAMS | Birth rate [$\rm{Gpc}^{-3}\rm{yr}^{-1}$]
CHE at ZAMS | 10.3 | 100 | 9.9 | 100 | 13.7
\- CHE triple | 7.6 | 73.5 | 7.5 | 75.7 | 9.6
\- Transition to classical evolution | 1.6 | 15.5 | 0 | 0 | 1.9
\- Merges during contact phase | 0.9 | 9 | 1.7 | 17.5 | 1.4
\- Simulation error | 0.2 | 2 | 0.6 | 6.7 | 0.2

In Table 1, we provide an overview of our sampled systems based on the
evolutionary type of the inner binary. Out of our sampled population of
triples, only about 10 per cent of the triples have an inner binary where both
stars evolve chemically homogeneously from ZAMS (CHE at ZAMS triples, see
Table 1), and we follow the further evolution only for these triples. About 75
per cent of CHE at ZAMS triples qualify as CHE triples and we focus on these
systems for the majority of the paper. For the remaining 25 per cent, we
distinguish three scenarios:

  * •

The inner stars transition to classical evolution. As the orbit of the inner
binary widens due to stellar winds, the rotational frequencies of the inner
stars decrease, because the stellar tides enforce synchronization between the
stellar spins and the (new longer) orbital period. If the inner orbit widens
sufficiently, the angular rotational frequencies of the inner stars drop below
$\omega_{\rm CHE}$ and therefore these stars transition to classical
evolution. This occurs only in our moderate metallicity model (15.5 per cent
of all CHE at ZAMS triples at Z = 0.005 and 0 per cent at Z = 0.0005).

  * •

The inner binary does not survive the contact phase during the MS phase of the
inner stars. We assume a merger takes place when both stars overflow their
outer Lagrangian point during the contact phase. This occurs during mass
equalization in the contact phase or due to GW emission, which lead to
shrinkage of the inner orbit. As orbital widening due to stellar winds prevent
mergers, the process occurs more efficiently at low metallicities (about 9 per
cent of all CHE at ZAMS triples at Z = 0.005 metallicity and 17.5 per cent at
Z=0.0005).

  * •

Computational issue. Finally, we note that the simulation of about 2 (6.7) per
cent of CHE at ZAMS triples fails at Z = 0.005 (Z = 0.0005). This can occur
because either no solution is found for the secular orbital evolution of the
system, or the computation time exceeds the allowed CPU time (which is 5000
seconds per system).

###  3.1 Main evolutionary outcomes

Table 2: A summary of the different channels (in bold font) and their sub-channels (in normal font) identified of CHE triples. The rows with bold fonts in the second and third column express the number of systems in each channel as a percentage of all CHE triple systems. The rows with normal fonts in the second and third column express the number of systems in each sub-channel as a percentage of all systems in their respective main channel.  Channel | % at Z = 0.005 | % Z = 0.0005 | Birth rate [$\rm{Gpc}^{-3}\rm{yr}^{-1}$]
---|---|---|---
No post-MS MT | 27.2 | 10.8 | 2.3
\- Inner binary evolution decoupled from tertiary | 97.4 | 98.6 | 2.2
\- Driven by three-body dynamics | 2.6 | 1.4 | 0.1
Stellar merger of the inner binary due to KL | 3.3 | 2.5 | 0.4
\- double helium star accretor | 74.7 | 58.5 | 0.3
\- helium star-MS accretor | 7.9 | 35 | 0.04
\- helium star-BH accretor | 17.4 | 6.5 | 0.06
Tertiary mass transfer | 55.4 | 52.1 | 5.2
\- MS-MS accretor | 58 | 60.8 | 3.1
\- BH-BH accretor | 30.9 | 24.3 | 1.5
\- Other types of accretors | 11.1 | 14.9 | 0.6
Unbound systems | 10.7 | 34 | 1.4
\- Core collapse SN | 100 | 16 | 0.9
\- (P)PISN | 0 | 84 | 0.5
Dynamical instability | 3.5 | 0.6 | 0.3

In Table 2, we show the most common evolutionary outcomes for CHE triples. We
distinguish 5 different evolutionary channels.:

  * •

No post-MS mass transfer phase: During the MS, it may be in a contact, but the
system does not experience any other form of mass transfer events. The inner
binary eventually forms a BH-BH binary in all these triples.

  * •

Stellar merger of the inner binary due to ZLK: Stellar merger occurs in the
inner binary due to ZLK oscillations.

  * •

Tertiary mass transfer (TMT): The tertiary star fills its Roche lobe.

  * •

Unbound systems: This evolutionary outcome takes place, if any of the stars
becomes unbound from the system. This occurs when a stellar remnant is formed
in the system, with three major subtypes: (i) natal kick imparted onto the
remnant object during the SN explosion, (ii) instantaneous mass loss during
pulsational PISN, or (iii) complete disruption of the star due to PISN.

  * •

Dynamical instability: These systems eventually become dynamically unstable,
where the secular approximation is no longer valid.

We discuss these channels in detail in sections 3.5 \- 3.3.

###  3.2 Examples for the evolution of a few selected systems

In the following, we present the evolution of a few selected systems from some
of the channels introduced in section 3.1. In all of these example systems,
the initial parameters of the inner binary are the same: $M_{\rm
1,ZAMS}=M_{\rm 2,ZAMS}=70\,M_{\odot}$, $a_{\rm in,ZAMS}=22.4\,R_{\odot}$.
These have been specifically chosen such that this system would form a GW
source via the binary CHE channel within the Hubble time, if it was an
isolated binary (i.e. in about 8.9 $\rm Gyrs$). The inner binary is tidally
locked and therefore $e_{\rm in,ZAMS}=0$. The stars of the inner binary are in
contact from ZAMS and equalise in mass soon after ZAMS. The initial mutual
inclination is $\rm{i}_{\rm ZAMS}=90^{\circ}$ in all systems discussed below,
which allows for ZLK oscillations to develop, unless they are suppressed by
short range forces (see e.g. Liu et al., 2015).

In order to understand the evolutionary paths of CHE triples introduced below,
we first show which configurations of CHE triples lead to efficient ZLK
oscillations (see Fig. 1). We evolve the previously introduced CHE inner
binary as an isolated system, and take four snapshots during different
evolutionary stages (ZAMS, end of MS, at the onset of core collapse, and at
the formation of an inner BH-BH binary). For each snapshot, we show a range of
possible tertiary companions to this inner binary with different tertiary
masses ($M_{\rm out}$) and outer semi major axes ($a_{\rm out}$) and identify
those regions, where three-body dynamics are relevant.

As shown in the leftmost panel, precession due to tides completely suppresses
three-body dynamics when the inner stars are still on MS for almost the entire
parameter space of CHE triples. The limited number of triples for which this
is not true typically become dynamically unstable later in the evolution (e.g.
compare panel 1 with panel 4). By the time of hydrogen depletion in the inner
stars, the stellar radii of CHE stars shrinks typically by a factor of 3-5
with respect to their ZAMS value. Therefore, at this stage tides become less
efficient (since $t_{\rm tides}\sim R^{-5}$, see equation 16) and precession
due to GR becomes the major limitation to three-body dynamics. For the systems
shown in Fig. 1, ZLK oscillations occur only, if $a_{\rm out}\lesssim
500\,R_{\odot}$ . During the CHeB phase of the inner stars, the the typical
timescale of precession due to GR further increases, as a result of the strong
Wolf-Rayet winds that significantly widen the inner orbit. As long as the
inner orbit widens faster than the outer orbit (which is always true if the
tertiary star is the initially least massive star in the system), the
timescale related to ZLK oscillations will not significantly increase.
Therefore during this stage, the parameter space where three-body dynamics are
relevant increases. This is also shown in the rightmost panel of Fig. 1; by
the time the inner binary forms BHs, triples with $a_{\rm out}\lesssim
2000\,R_{\odot}$ will develop ZLK oscillations.

####  3.2.1 Example for stellar merger of the inner binary due to ZLK
oscillations

First, we discuss the evolution of a CHE triple, in which the inner binary
merges as a double helium star due to strong ZLK oscillations (shown in Fig.
2). This triple has a tertiary with an initial mass of $M_{\rm
out,ZAMS}=32.1\,M_{\odot}$ and a circular outer orbit with $a_{\rm
out,ZAMS}=200\,R_{\odot}$. As indicated by Fig. 1, when the stars of the close
inner binary are still on the MS, precession associated with strong tides
suppresses the effects of the three-body dynamics (see also e.g. Vigna-Gómez
et al., 2022). At 3.9 Myrs, the stars of the inner binary evolve off the MS.
By this time, these stars had lost a small amount of mass due to stellar winds
and the inner orbit had widened by only 2 per cent as a result. Similarly, the
outer orbit also widens only by a negligible amount. Consequently, the
timescale of the ZLK oscillations does not change significantly. On the other
hand, the tidal effects become much weaker, as the radii of the stars had
decreased by a factor of 5 with respect to their ZAMS value. As a result, the
ZLK oscillations are no longer suppressed (see also second panel of Fig 1). At
this stage, there are two competing mechanisms that drive the evolution of the
pericenter: ZLK oscillations and the strong Wolf-Rayet-like winds, which
decrease and increase the pericenter, respectively. For this triple, the ZLK
timescale is extremely short (few years) and a large inner eccentricity of
$e_{\rm in}\approx 0.65$ is reached shortly after the onset of CHeB, during
which the orbital widening due to stellar winds is negligible. At this stage,
the pericenter becomes sufficiently small such that the helium stars fill
their Roche-lobes at the point of closest approach. We assume this results in
the merger of the inner binary.

####  3.2.2 Example for TMT towards an eccentric BH-BH binary

The next triple we discuss experiences a TMT episode towards an eccentric BH-
BH inner binary (shown in Fig 3). This system has the same parameters as the
previously discussed triple, but with a slightly larger initial outer
semimajor axes: $a_{\rm out,ZAMS}=421\,R_{\odot}$. When the inner stars evolve
off MS, ZLK oscillations are quenched by precession due to GR (compare the
second panels of Fig. 1 and Fig. 3). Three-body dynamics become later
effective, as the orbit of the inner binary widens significantly and faster
than the outer orbit due to strong WR winds (compare the third panels of Fig.
1 and Fig. 3, although by this stage the parameters of the inner binary differ
slightly). As a result $t_{\rm GR}$ increases by a factor of 5, while $t_{\rm
ZLK}$ barely changes. As ZLK cycles become only effective once the inner orbit
has sufficiently widened, the inner binary does not come into contact despite
reaching similarly high inner eccentricities as in the previous system.

As the stars of the inner binary have the same mass, they co-evolve, and they
become stellar remnants at the same time. This occurs around 4.2 Myr, when the
inner eccentricity is $e_{\rm in,max}=0.75$. Since core-collapse occurs in an
eccentric orbit, large range of possbile post-supernova orbits are possible
($a_{\rm in}=$ 42-186 $R_{\odot}$) depending on where exactly the stars are in
their orbit. In the particular example shown in Fig 3, the core collapse
occurs while both stars are near the pericenter (which is less likely as they
spend more time near the apocenter). This leads to an inner semi-major axis of
$a_{\rm in}=171\,R_{\odot}$ after BH-BH formation. As the outer orbit is
circular at the onset of the core-collapse, it only widens by a moderate
amount. As the inner period to outer period ratio has increased by a factor of
7, the timescale of the ZLK oscillations also further decrease, making the
three-body dynamics even more relevant for the further evolution of the
system. The evolution of this triple therefore demonstrates, that if the ZLK
oscillations are strong enough to induce eccentricities before the formation
of an inner BH-BH binary, the importance of three-body dynamics can be
significantly increased during the last stages of the evolution of the triple,
depending on (i) where the inner stars are in their orbit when the formation
of the compact objects occur and (ii) on the eccentricity of the outer orbit.

After the formation of the inner BH-BH binary, the tertiary star evolves off
MS, and at 6.1 Myr fills its Roche-lobe and transfers mass to the highly
eccentric ($e_{\rm in}=0.94$) BH-BH binary at a highly inclined orbit
($\rm{i}=71.5^{\circ}$). At this stage, we stop the simulation (but see later
section 5, where we predict the further evolution of some of these systems).
We note, however, that even if the TMT episode does not affect the inner
binary, it still merges due to GWs about a factor of 8 faster than its
isolated binary counterpart, just alone due to the high eccentricities induced
by the ZLK oscillations.

####  3.2.3 Example for TMT towards a circular BH-BH binary

Next, we show the evolution of a CHE triple, which also experiences a TMT
episode towards a BH-BH binary, but in which three-body dynamics remain
suppressed throughout the entire evolution. The initial outer semimajor axis
is $a_{\rm out,ZAMS}=1069\,R_{\odot}$. For this system the timescales of the
ZLK oscillations remain too long with respect to the timescale associated with
precession due to GR effects throughout the entire post-MS phase. At the onset
of the core-collapse, at which the parameter space for ZLK oscillations is the
typically the largest for CHE triples with inner binaries composed of non-
compact objects, the outer semimajor axis is $a_{\rm
out,ZAMS}=1720\,R_{\odot}$ and the tertiary mass is $M_{\rm
out}=31.9\,M_{\odot}$. Third panel in Fig. 1 implies that the three-body
dynamics is just quenched by the relativistic precession at this stage.
Therefore, the inner orbit remains circular when the BHs are formed, and the
inner orbit only widens moderately due to BH formation. The inner and the
outer orbit after the formation of a BH-BH binary are $a_{\rm
in}=46.6\,R_{\odot}$ and $a_{\rm out}=1860\,R_{\odot}$ and therefore the ZLK
oscillations remain quenched. At 6 Myr, the tertiary reaches a radius of 547
$R_{\odot}$ and fills its Roche-lobe while crossing the Hertzsprung gap. The
last two examples suggest (and we will show in section 4 that this is
generally true for the vast majority of CHE triples) that three-body dynamics
are only relevant for the evolution of CHE triples, if the tertiary star is on
a sufficiently short orbit, such that it will eventually fill its Roche-lobe
and initiate a TMT episode. Conversely, if the tertiary star remains detached
throughout the evolution of the triple, the inner binary evolves effectively
as an isolated binary for the vast majority of CHE triples.

Figure 1: We illustrate the parameter space where ZLK oscillations develop for
typical CHE triples at different evolutionary stages at Z = 0.005. Each panel
represents triples with a specific inner binary (with masses $M_{1}$, $M_{2}$
and inner separation $a_{\rm in}$ as indicated in the top left of each panel)
but with varying tertiary masses, $M_{\rm out}$ (x-axis) and outer semimajor
axes, $a_{\rm out}$ (y-axis, logscale). The parameters of the inner binary in
each panel are the same as that of an isolated CHE binary with $M_{\rm
1,ZAMS}$ = $M_{\rm 2,ZAMS}$ = $70\,M_{\odot}$, $a_{\rm ZAMS}=22.4\,R_{\odot}$
at different evolutionary stages (see text in section 3.2). First panel: two
CHE main sequence stars at zero-age main sequence, second panel: helium stars
at the onset of core-helium burning (HeMS), third panel: two helium stars at
the end of core-helium burning (HeRGB), fourth panel: at the formation of a
BH-BH inner binary. We assume circular inner and outer orbits and a mutual
inclination of $\rm{i}=90^{\circ}$. The light blue and the dark blue lines
show regions, where the ZLK timescale equals to the timescale of apsidal
precession due to tides and GR effects, respectively (see equation 15, 17, and
16). Consequently, triples above any of these two lines do not exhibit ZLK
oscillations, as they are quenched by these short range forces. The green line
represents the boundary of dynamical instability. The shaded grey region
represents triples where ZLK oscillations are effective. Figure 2: A schematic
drawing showing the evolution of a triple system in which the stars in the
inner binary experience stellar merger due to ZLK oscillations. The first line
below each drawing shows the evolutionary stage of each star. The first is for
the primary, the second is for the secondary, and the third is for the
tertiary star. CHE MS is chemically homogeneously evolving MS star, HeMS is
core-helium burning helium star, HeRGB is a helium star that has finished
core-helium burning. The parameters of the triple at ZAMS are: $M_{\rm
1,ZAMS}=M_{\rm 2,ZAMS}=70\,M_{\odot}$, $a_{\rm in,ZAMS}=22.4\,R_{\odot}$,
$\rm{i}_{\rm ZAMS}=90^{\circ}$, $M_{\rm out,ZAMS}=32.1\,M_{\odot}$, $a_{\rm
out,ZAMS}=200\,R_{\odot}$, $e_{\rm out}=0$. The outer eccentricity remains
$e_{\rm out}\lesssim 0.01$ throughout the evolution.  Figure 3: A schematic
drawing showing the evolution of a triple system with that eventually
experiences a TMT episode (i.e. the most common evolutionary outcome, see
Table 2). The triple system with the following initial parameters: $M_{\rm
1,ZAMS}=M_{\rm 2,ZAMS}=70\,M_{\odot}$, $a_{\rm in,ZAMS}=22.4\,R_{\odot}$,
$\rm{i}_{\rm ZAMS}=90^{\circ}$, $M_{\rm out,ZAMS}=32.1\,M_{\odot}$, $a_{\rm
out,ZAMS}=200\,R_{\odot}$, $e_{\rm out}=0$. The outer eccentricity remains
$e_{\rm out}\lesssim 0.01$ throughout the evolution.

###  3.3 No post-MS mass transfer

In these triples, the tertiary star remains bound and detached, while the
stars of the inner binary form a BH-BH binary. The inner stars are in contact
in the majority of the cases (e.g. around 90 per cent at Z = 0.005). There are
no any other mass transfer phases during the evolution of these systems (by
definition). About 27 per cent of CHE triples evolve this way in our moderate
metallicity model (see Table 2). This decreases to 11 per cent at Z = 0.0005.
The main reason for this difference is the larger number of PISN that occurs
at lower metallicities, which prevent the formation of BHs.

After the formation of the BH-BH binary, the system may merge due to GW
emission within a Hubble time. This occurs for all systems of this type at
$Z=0.0005$. However at $Z=0.005$, the stellar winds are strong enough such
that 32 per cent of the inner binaries of these triples end up with orbits
that are too wide to merge within a Hubble time due to GW emission. We note
that these are not necessarily all of the GW sources from our simulations, as
triples in other channels discussed here can also potentially form merging
binary BHs (see discussion in section 5).

For the majority of these triples ($>97$ per cent), the inner binary evolves
essentially unaffected by the tertiary star (see also section 4). Therefore,
the properties of the inner binaries of this channel are nearly
indistinguishable from those of isolated CHE binaries. The initial outer
pericenters of the triples of this channel are large enough such that the
outer star remains detached (i.e. $a_{\rm p,out,ZAMS}\gtrsim
2000$-$3000\,R_{\odot}$ at Z = 0.005 , see also section 4). At such large
tertiary separations, the three-body dynamics remain suppressed during the
entire evolution of the triple.

The properties of the subgroup in which three-body dynamics drive the
evolution of the inner binary are very different. Firstly, they have very
short initial outer pericenters (i.e, $a_{\rm p,out,ZAMS}\approx
100$-$700\,R_{\odot}$), and secondly, the tertiary has a relatively low mass
(typically $M_{\rm out,ZAMS}$ = $10$-$30$ $M_{\odot}$). In these systems, the
ZLK oscillations drive the eccentricity of the inner BH-BH binary up to large
values (e.g. $e_{\rm in}\gtrsim 0.7$-$0.9$). Above a given eccentricity, the
GW emission becomes so efficient that the inner binary decouples from the
tertiary and it plunges due to GWs (see e.g. Silsbee & Tremaine, 2017;
Rodriguez & Antonini, 2018). These systems typically have a relatively low-
mass tertiary star compared to the stars in the inner binary, such that the
inner binary merges as a BH-BH binary due to GW emission before the tertiary
star would evolve off the MS and fill its Roche-lobe. Overall, the parameter
space for this subgroup is very small, and therefore we predict a negligible
GW merger rate (see later discussion in section 5).

###  3.4 Stellar merger of the inner binary due to ZLK

In this scenario, the inner binary merges due to three-body dynamics, before
it would form a BH-BH binary. At Z = 0.005, about 3.3 per cent of the CHE
triple population evolves this way. In our low metallicity model, this
fraction decreases slightly, to 2.2 per cent. This is because at lower
metallicities, the inner period to outer period ratio increases less due to
the weaker stellar winds, and therefore ZLK oscillations remain less efficient
(see equation 15).

Mergers in this channel occur in inner binaries, in which one or both of the
stars have already evolved off MS, otherwise the strong tidal effects
typically quench the ZLK oscillations (see section 3.2). As shown in Table 2,
most of the merger occurs between two helium stars (59-75 per cent). The rest
occurs between a helium star - MS star or helium star - BH binaries. The
majority of the double helium star mergers (>90 per cent) originate from
triples, in which the stars in the inner binary were in contact during MS and
co-evolved. This also implies that the majority of them have equal masses at
the time of the merger. The masses of these helium inner stars typically range
from 29 to 94 $M_{\odot}$ at Z = 0.005.

The outer orbital period of the triples from this channel has to be
sufficiently short, such that the ZLK oscillations are strong enough such that
they prompt the inner binary to merge. The outer pericenter at the moment of
the merger typically ranges from 100 to 200 $R_{\odot}$ and it does not
exceeds 700 $R_{\odot}$. The eccentricities of the inner binary at the moment
of the merger typically have values of $e_{\rm in}\approx$ 0.5-0.9. For all of
these triples, the tertiary is a MS star at the time of the merger and less
massive than the stars of the inner binary, otherwise it would evolve faster
than the stars in the inner binary and would fill its Roche-lobe, while the
inner stars are still on MS. If the outer orbit does not significantly change
after the merger, the tertiary star is expected to transfer mass to the merger
product, once it evolved off the MS.

###  3.5 Systems with tertiary mass transfer (TMT)

Among CHE triples, this is the most common evolutionary pathway. In these
systems, the outer star eventually initiates a mass transfer phase while the
inner binary is detached or in contact. Approximately 55 (52) per cent at
$Z=0.005$ ($Z=0.0005$) of all CHE triples experience this type of evolution
(see Table 2). This means that a TMT episode would eventually occur in about
40 per cent of all stellar systems containing a binary with CHE stars (with
$f_{\rm binary}=0.21$, $f_{\rm triple}=0.73$). While systems containing
binaries with CHE stars are rare (see e.g. typical birth rates in Table 1),
they form GW sources very efficiently (e.g. Mandel & de Mink, 2016; de Mink &
Mandel, 2016; Marchant et al., 2016). Therefore, our predictions suggest that
the evolution of a non-negligible fraction of potential GW progenitors could
experience a TMT episode. This is an interesting result, as TMT is thought to
be very uncommon for classically evolving hierarchical triples, which would
have implied that they play a limited role in important astrophysical
phenomena (see e.g. de Vries et al., 2014; Toonen et al., 2020; Hamers et al.,
2022a; Kummer et al., 2023). In particular Toonen et al. (2020) found that
about 1 per cent of triples with primaries in the intermediate mass range
belong to this evolutionary channel. Similarly, de Vries et al. (2014)
predicts that about only 1 per cent of the observed 725 triples in the
catalogue of Tokovinin (2010) would eventually initiate TMT.

In the following sections (3.5.1-3.5.6), we discuss the properties of the
triples of this channel at the onset of TMT. While predicting the outcome of a
TMT episode is currently extremely challenging, highlighting several important
aspects of these systems (e.g. dynamical stability of TMT, timescales of TMT
epsiodes, the amount of transferred mass, the type of accretors, etc.) helps
to better understand the nature of these systems and the role they potentially
play in the evolution of GW progenitors.

####  3.5.1 Donors of TMT episodes

Here, we discuss the stellar evolutionary stage of the donor star at the onset
of the mass transfer episode. as it is highly relevant for determining if the
mass transfer episode occurs in a dynamically stable or unstable way (which
has dramatic effect on the outcome of the TMT, compare e.g. de Vries et al.,
2014; Glanz & Perets, 2021b).

In particular, convective envelopes can be developed by core-helium-burning or
asymptotic giant branch stars. Mass transfer episodes initiated by such cool-
giant donors with deep convective envelopes are more likely to occur in a
dynamically unstable way than mass transfer phases initiated by giant donors
with mostly radiative envelopes (e.g., Hjellming & Webbink, 1987; Soberman et
al., 1997; Klencki et al., 2020).

Figure 4: We show the properties of the donor star at the onset of tertiary
mass transfer episode. The left panels shows the amount of mass transferred
from the tertiary to the inner binary for systems undergoing TMT at Z = 0.005,
and the right panels show the outer mass ratio, $q_{\rm out}=M_{\rm
out}/(M_{1}+M_{2})$, at the onset of the tertiary mass transfer phase.
Different colours correspond to different evolutionary stages of the donor as
shown by the legend. We exclude cases where the donor is a MS star.
Furthermore, we distinguish BH-BH inner binaries (upper panels) and MS-MS
inner binaries (lower panels), which are the main types of accreting systems.
The unfilled histogram in the panels on the left shows the transferred mass as
a fraction of the total mass of the inner binary for all donor types. All
histograms shown have been normalised with respect to the population of CHE
evolving triples.

At $Z=0.005$, around 80 per cent of the donors of TMT systems are stars
crossing the Hertzsprung gap. At this metallicity, the largest expansion in
the radius of the star occurs during this evolutionary phase, which makes
binary interaction during this stage the most probable. The second most common
donor type is CHeB star with 11.3 per cent, while the rest are either stars on
the first giant branch (when the tertiary $M_{\rm out,ZAMS}\lesssim
8\,M_{\odot}$) or stars on the asymptotic giant branch.

At lower metallicities, CHeB donors are more prevalent. At $Z=0.0005$, only 58
per cent of the tertiary donors are HG stars while 40 per cent are CHeB stars;
this is because the onset of CHeB occurs at a higher effective temperature
with respect to systems at Z = 0.005. Consequently, at lower metallicities,
the onset of CHeB is followed by a larger increase in radius with respect to
their higher metallicity counterparts. This in turn implies that stars are
more likely to fill their Roche-lobes at this evolutionary stage.

####  3.5.2 Stability of TMT episodes

The vast majority of mass transfer episodes in this channel occur in a
dynamically stable way (99.9 per cent at Z = 0.005 and 98.8 per cent at Z =
0.0005). This is due to the relatively low mass ratios at the onset of the
mass transfer phase (i.e. typically $q_{\rm out}<q_{\rm crit}$, see right
panel of Fig. 4 for our moderate metallicity model, and Fig. 13 for our low
metallicity model). Typical mass ratios for systems with HG donors are $q_{\rm
out}$ = 0.4-0.8, while for CHeB donors, they are $q_{\rm out}$ = 0.3-0.5 . The
values for CHeB donors are smaller because of the strong LBV winds that CHeB
star experience decrease the mass ratios over time. Unstable mass transfer
phases exclusively occur with CHeB donors in our simulations.

These low mass ratios also imply that the expansion due to stellar evolution
drives the TMT episodes (e.g. Soberman et al., 1997). Consequently, we expect
TMT episodes with HG donors to last $10^{4}$ yrs, while TMT epiosdes with CHeB
donor could last much longer up to $10^{5}$-$10^{4}$ years.

####  3.5.3 Accretors of TMT episodes

In this subsection, we discuss the type of accretors of TMT episodes. The
evolutionary stage of the inner binaries has a crucial role in the outcome of
TMT episodes. If the inner binary comprises CHE MS stars, a TMT episode
probably leads to the merger, as CHE binaries have very short periods and the
majority of them are in contact at the onset of the TMT (see also Braudo et
al., 2022). On the other hand, if the inner binary consists of BHs, TMT
epsiode is unlikely to lead to merger by itself, however, in principle, could
be a source of (observable) X-ray emission (e.g. Kaaret et al., 2017).

As shown in Table 2, the two most common types of accretors are MS-MS and BH-
BH binaries. In only 11-15 per cent of CHE triples experience TMT with
different accretors, such as an inner binary consisting of two helium stars or
a helium star with a MS or BH companion.

We highlight the relatively large fraction of BH-BH accretors (24-31 per cent
of CHE triples experiencing TMT). For classically evolving triples, mass
transfer towards a BH-BH binary is highly unlikely. Firstly, in systems in
which a TMT episode were to occur towards a BH-BH inner binary, the stars of
the inner binary need to be more massive than the tertiary, such that they
form BHs before the outer star fills its Roche-lobe. Secondly, the outer star
has to be sufficiently close, otherwise it would remain detached throughout
its evolution. This, in turn, puts a limit on the largest possible inner
orbit, if the system is to remain dynamically stable. The maximum inner orbit
for such systems is so small that classically evolving inner stars (which
eventually expand) would initiate mass transfer and would most likely merge,
which would reduce the triple to a binary and a tertiary mass transfer would
never occur (see e.g. Fig. 14 in Toonen et al., 2020). On the other hand, if
the triple has CHE inner stars, the stars will not expand and not merge with
one another, instead the system will evolve to contain a BH-BH binary by the
time the tertiary fills its Roche-lobe.

####  3.5.4 Mass transferred towards the inner binary

We discuss the amount of mass that is transferred during the TMT episode. This
is an important aspect, as the relative transferred mass determines angular
momentum reservoir available to change the orbit of the inner binary.

Assuming that the entire envelope of the donor star is transferred towards the
inner binary, the amount of transferred mass ranges between 1-40$\,M_{\odot}$
for BH-BH accretors and between 10-50$\,M_{\odot}$ for MS-MS accretors (see
left panel of Fig. 4 for Z = 0.005 and 13 in section of A of Appendix for Z =
0.0005). Systems with MS-MS accretors typically receive a larger amount of
mass than BH-BH accretors, because the tertiary star is typically more massive
in the former case. This is because for the tertiary to fill its Roche lobe
while the inner stars are still on the MS, the initial tertiary star needs to
be evolve faster and hence be more massive than the MS stars. The relative
transferred mass expressed as a fraction of the total mass of the inner binary
(i.e. $M_{\rm transferred}/M_{\rm tot,inner}$) has the same maximum value
($\sim$ 0.5) for both BH-BH and MS-MS accretors (see grey histogram in left
panel of Fig. 4).

####  3.5.5 Formation of circumbinary disc

We discuss how common it is for TMT systems to develop a circumbinary disc at
the onset of the mass transfer episode. As explained in section 2.5, whether a
TMT episode is accompanied by a formation of a circumbinary disc can have
important consequences for the evolution of the inner orbit.

We find that about 63 per cent of all TMT systems develop circumbinary discs
in our moderate metallicity model, while in the rest TMT proceeds in a
ballistic fashion. Systems in which a circumbinary disc is formed during the
TMT phase typically have larger outer pericenters at the onset of the mass
transfer ($a_{\rm p,out}\approx$ 300-6000 $R_{\odot}$) than in those where TMT
proceeds in a ballistic manner ($a_{\rm p,out}\approx$ 100-600 $R_{\odot}$).

TMT with circumbinary disc is more prevalent at lower metallcities. About 74
per cent of all TMT systems develop circumbinary discs at Z = 0.0005. This
occurs because the ratio of the inner and the outer orbital separation
decreases less by the onset of the mass transfer phase due to weaker stellar
winds (see equation 24).

TMT episodes with inner BH-BH binaries are somewhat more likely to occur in a
ballistic fashion than with MS-MS inner binaries. About 45 (23) per cent of
TMT systems with BH-BH inner binaries do not develop circumbinary discs at Z =
0.005 (Z = 0.0005), while 32 (27) per cent of TMT episodes with MS-MS inner
binaries occur in a ballistic fashion. This is mainly because the inner
apocenter to outer pericenter ratios at the onset of TMT are typically higher
for inner BH-BHs than for inner MS-MS binaries (see equation 24). This
difference is due to Wolf-Rayet winds, supernova kicks and possible ZLK
oscillations that BH-BH inner binaries experienced prior to the TMT episode.

####  3.5.6 Three-body dynamics prior to TMT

Three-body dynamics can increase the eccentricities of the inner binary. This
can, for example significantly decrease the coalescence time due to GWs (e.g.
Miller & Hamilton, 2002; Blaes et al., 2002; Wen, 2003; Thompson, 2011).

Three-body dynamics are almost always suppressed during the MS phase of the
inner binaries due to the strong tides (see also section 3.2). Consequently,
the inner orbits of TMT systems with MS-MS inner binaries are always circular
at the onset of the mass transfer episode. On the other hand, this is no
longer the case when the inner stars are in their post-MS. In Fig. 5, we show
the cumulative distribution of the inner binary eccentricities at the onset of
the mass transfer phase of TMT systems with BH-BH accretors at Z = 0.005. We
see that systems without circumbinary discs tend to have eccentric inner
orbits at the onset of mass transfer. The high eccentricities are caused by
ZLK cycles during the post-MS evolution of the inner binary. About 40 per cent
of such triples have $e_{\rm in}\gtrsim 0.4$ at this stage. This is in
contrast with the systems with circumbinary discs; about 90 per cent of the
systems have eccentricities $e_{\rm in}\lesssim 0.1$. The difference is due to
the smaller inner period to outer period ratios that systems without
circumbinary discs have (see equation 15). In our low metallcity model, high
eccentricites at the onset of TMT are much less common (see Fig. 14 in section
of A of Appendix). For these systems the inner period to outer period ratio
does not increase significantly because of the weak stellar winds.

Figure 5: The cumulative distribution of the inner binary eccentricities of
systems experiencing TMT with BH-BH accretors at the onset of the mass
transfer phase at Z = 0.005. The green curve shows the systems with accretions
disc formed during the mass transfer phase. The blue curve shows the systems
where ballistic accretion occurs.

###  3.6 Unbound systems

In this channel, one of the stars in the triple becomes unbound as a result of
core-collapse. We distinguish systems based on whether this occurs via PISN or
via classical core collapse (e.g. Fryer et al., 2012). As shown in Table 2,
PISN does not occur in our moderate metallicity model, whereas at Z = 0.0005,
it becomes quite prevalent; about 84 per cent of the unbound systems occur due
to PISN.

If the triples becomes unbound as a result of a classical core-collapse, we
further distinguish whether it is due to the core-collapse occuring in the
inner binary (97 per cent of all classical core-collapse systems at Z = 0.005
and 99 at Z = 0.0005) or of the tertiary star (3 per cent at Z = 0.005 or 1
per cent at Z = 0.0005). As the inner binary consists of CHE stars, they have
large initial masses (i.e. $M_{\rm ZAMS}\gtrsim 30\,M_{\odot}$) and
furthermore they develop more massive CO cores than their classically evolving
counterparts. Therefore, they get weak (if any) natal kicks when they form BHs
according to our implemented natalk kick prescription Yet weak natal kicks, or
even completely symmetrical instantaneous mass losses due to neutrino losses
(which we assume to be 10 per cent of the pre-collapse mass according to Fryer
et al., 2012) can unbind the tertiary star, if the outer star has high
eccentricities. We find that in systems in which one of the stars becomes
unbound due to the core-collapse in the inner binary, the outer eccentricities
are large, about 70 per cent of them $e_{\rm out}\geq 0.8$. In the vast
majority of the cases (about 99 per cent of such unbound systems), only the
tertiary is ejected, while the inner binary remains bound.

If the triple becomes unbound due to the core-collapse of the tertiary star
and with low outer eccentricity, it almost always occurs as a result of a
strong natal kick. Consequently, most of such unbound systems have initial
tertiary masses of $M_{\rm out,ZAMS}\approx 8$-$25\,M_{\odot}$ (see also
discussion in section 4), as these systems are expecetd to receive the largest
kicks according the supernova prescription of Fryer et al. (2012) .

###  3.7 Systems which become dynamically unstable

These triples typically have very short initial outer pericenters ($a_{\rm
p,out,ZAMS}\approx 70$-$400\,R_{\odot}$) and therefore are very close to the
stability limit at ZAMS. Such systems can transition to non-secular or non-
hierarchical evolution, if $a_{\rm in}/a_{\rm out}$, $e_{\rm out}$ or $q_{\rm
out}$, significantly increases during evolution (see Mardling & Aarseth,
2001). Among CHE triples, there are primarily two processes that can trigger
this change: stellar winds and core collapse.

If the relative wind mass loss rate (e.g. $\dot{M}/M$) in the inner binary is
higher than that of the tertiary star, $a_{\rm in}/a_{\rm out}$ and $q_{\rm
out}$ will increase, which can prompt the triple to experience a dynamical
instability (see also Kiseleva et al., 1994; Iben & Tutukov, 1999; Perets &
Kratter, 2012; Toonen et al., 2022). 30 per cent of the systems of this
channel destabilise due to stellar winds and the destabilisation occurs when
the stars of the inner binary are in their post-MS phase. At this stage, the
inner stars experience strong Wolf-Rayet winds, while the tertiary star is
still on the MS with significantly lower mass loss rates.

In the remaining 70 per cent, the instability sets in due to the collapse of
the core of one of the stars. We find that this only occurs when BH formation
takes place in the inner binary. As noted in section 2.4, CHE stars typically
form BHs via direct collapse, such that $q_{\rm out}$ only increases slightly.
Furthermore, the direct collapse is expected to be accompanied by a weak
Blauw-kick due to neutrino losses such that $a_{\rm in}/a_{\rm out}$ and
$e_{\rm out}$ only increase significantly, if the inner or the outer pre-core-
collapse orbits are eccentric, respectively. The pre-core-collapse inner orbit
is eccentric in 72 per cent of the systems of this channel. And the
eccentricity is caused by ZLK oscillations. In the remaining 28 per cent,
three body-dynamics is not efficient in driving up the eccentricity because
the mutual inclination is outside of the critical Kozai range (see e.g. Naoz,
2016). Therefore, the core collapse occurs when the inner orbit is circular.
These systems still become unstable during the BH formation, because 1) either
$a_{\rm in}/a_{\rm out}$ already increased strongly due to stellar wind mass
losses before the BH formation or 2) the outer orbit is eccentric and the core
collapse occurs while the tertiary star is near the outer pericenter (leading
to a significant increase in $e_{\rm out}$).

The occurrence rate of this channel is strongly dependent on metallicity (3.5
per cent of all CHE triples at Z = 0.005 and 0.7 per cent at Z = 0.0005, see
Table 2). This dependence is due to the reduced strength of stellar winds and
ZLK oscillations (which are responsible for any eccentricity in CHE inner
binaries) at lower metallicities.

##  4 The origin of each evolutionary channel

In this section, we discuss the initial parameters of the triples from each
evolutionary channel introduced in section 3.1. We find that initial
parameters can be used as a proxy to determine the final evolutionary outcome
of CHE triples. In particular, the evolutionary outcome can be parameterised
by the initial mass and orbital separation of the tertiary star. The
parameters of the inner binary play a less important role in this regard, as
the parameter space for CHE inner binaries is already quite reduced. We
illustrate this in the left panel of Fig. 6 by showing an ensemble of CHE
triples at $Z=0.005$, in which the parameters of the inner binary are the
same, but the mass and the orbital separation of the tertiary star are varied
(therefore this grid represents only a small subset of the entire CHE
population discussed in section 3.1). The inner binary consists of two 70
$M_{\odot}$ stars and with a circular initial orbit with $a_{\rm
in,ZAMS}=22.4\,R_{\odot}$ (similarly to the example systems discussed in
section 3.2). The initial tertiary mass ranges from 5 to 100$\,M_{\odot}$,
while $a_{\rm out,ZAMS}$ ranges from 200 to $10^{4}\,R_{\odot}$.

###  4.1 Initial parameters of systems of different evolutionary channels

The majority of the triples shown in the left panel of Fig. 6 experience TMT
episodes. Their initial outer orbital separations are relatively short and
range roughly from 100 to 3300 $R_{\odot}$. The evolutionary phase of the
inner stars at the onset of the TMT episode depends on the initial mass of the
tertiary star. For the systems shown in the left panel of Fig. 6, the inner
binary at the onset of TMT comprise of BHs, if $M_{\rm out,ZAMS}\lesssim
59\,M_{\odot}$, helium stars, if $59\,M_{\odot}\lesssim M_{\rm out,ZAMS}\leq
70\,M_{\odot}$, and MS stars, if $M_{\rm out,ZAMS}\geq 70\,M_{\odot}$. The
majority (53 per cent) of the TMT systems in the left panel of Fig. 6 have a
BH-BH inner binaries. For the entire population of CHE triples presented in
section 3.1, the same percentage is smaller (i.e 31 per cent) at the same
metallicity (see Table 1). As shown in Fig. 7, this quantity (i.e. the ratio
of the number of TMT systems with BH-BH inner binaries and the number of all
TMT system) scales proportionally to the initial mass of the secondary star in
the inner binary. This means that TMT episodes occur more frequently with BH-
BH accretors among CHE triples with more massive inner stars. This is due to
our assumptions about the initial distribution of the triples (section 2.6).
If the TMT occurs towards a BH-BH inner binary, the tertiary has to be
initially the least massive in the triple. With increasing $M_{\rm 2,ZAMS}$,
the fraction of triples for which $M_{\rm 2,ZAMS}>M_{\rm out,ZAMS}$ increases
because of our assumptions of a maximum initial stellar mass of $M_{\rm
ZAMS,max}=100\,M_{\odot}$ and a flat outer mass ratio distribution.

In 15 per cent of the triples shown in the left panel of Fig. 6, the inner
binary merges before BH formation or before a TMT episode occurs. All such
mergers in the grid occur between two helium stars, and are due to ZLK
oscillations that arise when the stars of the inner binary evolve off the MS.
The initial outer orbital separations in this channel are very short, i.e. 200
to 241$\,R_{\odot}$, while the tertiary masses range between $32\leq M_{\rm
out,ZAMS}/M_{\odot}\leq 68$. For lower tertiary masses ($M_{\rm out,ZAMS}<32\
M_{\odot}$), the ZLK oscillations are not strong enough to boost the inner
eccentricity and cause a mass transfer episode in the inner binary. For larger
tertiary masses ($M_{\rm out,ZAMS}>70\ M_{\odot}$), the tertiary typically
fills its Roche-lobe before the stars of the inner binary evolve off the MS.
However, during the main sequence phase of the inner stars, the effects of ZLK
cycles are quenched and consequently no mergers are prompted by three-body
dynamics before the tertiary initiates a TMT episode.

Triples of the no post-MS MT channel in the left panel of Fig. 6 have initial
outer orbits $a_{\rm{out}}\gtrsim 2000$-$3000\ R_{\odot}$. Their initial
tertiary mass is also typically outside of the range of $\sim$8-25
$M_{\odot}$, such that the system does not dissociate due to SN kicks. As we
shown in the next subsection, three-body dynamics are not important for the
evolution of these systems.

In left panel of Fig. 8, we show the inital pericenter ($a_{\rm outer,ZAMS}$)
distribution of the entire CHE triple population for each evolutionary channel
at Z = 0.005. As it can be seen, the range of initial pericenters are in
agreement with those shown in Fig. 6 for all channels except for the unbound
systems (since for the unbound systems, the outer eccentricity plays a crucial
role, as explained in section 3.6, and in the grid we assume circular outer
orbits). This again confirms that the parameters of the tertiary star play the
most important role in determining the evolutionary path of a CHE triple. As
shown in left panel of Fig. 8, the range of $a_{\rm outer,ZAMS}$ of systems
with TMT episodes increases with decreasing metallicity. At lower metallicity,
the stellar winds are weaker and consequently, the outer orbit widens less.
Therefore, the maximum $a_{\rm outer,ZAMS}$ at which the tertiary stars can
still fill their Roche-lobes also increases with decreasing metallicity.

###  4.2 Initial parameters of triples with three-body dynamics

In the right panel of Fig. 6, we show the maximum eccentricities that the
inner binaries reach during their evolution ($e_{\rm in,max}$). About 29 per
cent of the triples shown in the right panel of Fig. 6 reach $e_{\rm
in,max}\geq 0.4$ due to ZLK cycles. In all of these triples, the tertiary star
eventually fills its Roche-lobe (although in some cases, the inner binary
merges first).

For the systems shown in Fig. 6, ZLK cycles are efficient when $a_{\rm
out,ZAMS}\lesssim 1200\,R_{\odot}$ and $M_{\rm out,ZAMS}\lesssim
70\,M_{\odot}$. When the outer orbit is $a_{\rm out,ZAMS}\gtrsim
1200\,R_{\odot}$, the ZLK cycles are quenched by various short range forces
(e.g. precession caused by tides or general relativistic effects). If $a_{\rm
out,ZAMS}\lesssim 1200\,R_{\odot}$ but $M_{\rm out,ZAMS}\gtrsim
70\,M_{\odot}$, the tertiary star fills its Roche-lobe while the stars in the
inner binary are still on the MS. The inner binaries of these triples do not
develop high eccentricities, as ZLK cycles are quenched during MS due to
strong tides (see section 3.2), and TMT episode with MS-MS accretors are
expected to result in the merger of the inner binary (see section 5).

The right panel of Fig. 6 also shows that $e_{\rm in,max}$ does not decrease
smoothly with decreasing outer orbital separations, instead it drops rather
abruptly across $a_{\rm out,ZAMS}\approx 1200\,R_{\odot}$. Triples with
$a_{\rm out,ZAMS}\approx 1200\,R_{\odot}$ reach very large inner eccentricties
($e_{\rm in,max}\approx 0.9$), while at slightly larger orbital separations
(i.e. $a_{\rm out,ZAMS}\approx 1500\,R_{\odot}$) the ZLK cycles are completely
quenched.

These above mentioned effects are qualitatively also true for the entire CHE
triple population presented in section 3.1 (see right panel of Fig. 8). At Z =
0.005, the ZLK oscillations are only efficient, if $a_{\rm p,out,ZAMS}\lesssim
1200\,R_{\odot}$. This implies that three-body dynamics are only relevant for
those triples, in which the tertiary star would eventually fill its Roche-lobe
(compare right and left panel of Fig. 8). Consequently, if the tertiary in a
CHE triple remains detached throughout its evolution, the evolution of the
inner binary will almost always be kinetically decoupled from the tertiary
star. If $a_{\rm p,out,ZAMS}\lesssim 1200\,R_{\odot}$, a wide range of inner
eccentricites are possible ($e_{\rm in,max}=$ 0-0.9) for all $a_{\rm
p,out,ZAMS}$. In this case, the value of $e_{\rm in,max}$ is primarily
determined by the mutual inclination of the triple (see also e.g. Anderson et
al., 2017).

In our low metallicity model ($Z=0.0005$) the maximum initial outer pericenter
at which three-body dynamics are still relevant is lower compared to our
moderate metallicity model (right panel in Fig. 16 in section A of Appendix).
At such low metallicities, stellar winds do not widen the orbit of the inner
binary significantly and thus the timescales of the ZLK cycles do not decrease
as much as at Z = 0.005.

Figure 6: Left panel: we show the evolutionary outcome of collection of
triples with fixed inner binary parameters and different parameters for the
tertiary at a metallicity $Z=0.005$. The parameters of the inner binary are
the same for each system shown in the grid, i.e. $M_{\rm 1,ZAMS}=M_{\rm
2,ZAMS}=70,M_{\odot}$, $a_{\rm in,ZAMS}=22.4\,R_{\odot}$, $e_{\rm in}=0$. The
initial tertiary mass $M_{\rm out}$ ranges from 5 to 100$\,M_{\odot}$ on a
linear scale, while $a_{\rm out,ZAMS}$ ranges from 200 to 1000$\,R_{\odot}$ on
a logarithmic scale. We assume zero outer eccentricity and an initial
inclination of $90^{\circ}$. Right panel: we show the maximum eccentricity of
the inner binary reached during the evolution. The parameters of the inner
binary are the same as in left panel Figure 7: The distribution of the initial
secondary masses of the entire CHE triple population that undergo tertiary
mass transfer at Z = 0.005. We only show those systems which have either BH-BH
or MS-MS accretors. The grey unfilled histogram, with the corresponding
secondary x-axis on the right hand side, shows the number of systems with BH-
BH accretors as a fraction of all systems undergoing tertiary mass transfer
phases.

Figure 8: Left panel: The distribution of the initial pericenter of a few
selected evolutionary types of the entire CHE triple population at Z = 0.005.
The histograms have been normalised to the full population of the CHE evolving
triples. The histograms shown are stacked. Each colour represents a different
evolutionary type. For clarity we do not show all the types introduced in
section 3.1, see text for discussion. Right panel: The distribution of the
initial outer pericenter of the entire CHE triple population at Z = 0.005. We
distinguish systems based on the maximum eccentricity of the inner binary
reached during their evolution. The histograms are normalised to one and
stacked.

##  5 Gravitational waves sources

We now discuss the possible formation channels of GW sources that originate
from CHE triples and their properties. In section 5.1 we predict the merger
rate densities and compare them to that of GW sources from isolated CHE
binaries. For this, we assume two test populations with different stellar
multiplicity fractions. One population is composed of only single and binary
stellar systems (i.e. with stellar multiplicity fractions at ZAMS of $f_{\rm
single}=0.3$, $f_{\rm binary}=0.7$, $f_{\rm triple}=0$), in the other triples
dominate ($f_{\rm single}=0.06$, $f_{\rm binary}=0.21$, $f_{\rm
triple}=0.73$).

In sections 5.2 \- 5.5 we discuss the properties of each GW formation channel
from CHE triples and binaries. These predictions are based on the synthetic
populations discussed previously, and in cases where the simulations are
stopped before the formation of a BH-BH binary, we predict the further
evolution of CHE triples beyond the stopping conditions (Section 2.6) by
applying simple assumptions (as detailed below).

Figure 9: The possible formation channels of merging binary BHs from our CHE
triples population.

The four main identified formation channels of GW sources within our CHE
triple population are (see also Fig. 9):

  * •

Effectively isolated inner binary: For such triples, three-body dynamics is
suppressed by various short-range forces and the tertiary star remains
detached throughout the entire evolution. The inner binary therefore evolves
effectively as an isolated binary and the properties of these GW sources are
indistinguishable from those of the CHE binary channel. There are two ways
these systems can form: i) with the tertiary star bound to the triple (systems
from the no post-MS MT channel, see section 3.3) and ii) systems in which the
tertiary star becomes unbound from the triple (from the unbound channel
discussed in section 3.6). For the latter, we assume that the orbit of the
inner binary is not affected by the tertiary unbinding from the triple system.

  * •

TMT with a BH-BH accretor: This channel comprises systems in which the
tertiary star fills its Roche-lobe when the inner binary is a BH-BH binary.
The inner binary components do not coalesce during the TMT phase, but will
merge afterwards due to GW emission. In these systems, the tertiary star can
affect the evolution of the inner binary in two major ways, via TMT episode
and via three-body dynamics (see section 3.5). In section 2.5, we introduced
our assumptions regarding the evolution of the inner orbiy during a TMT
episode.

  * •

TMT with a MS-MS accretor: In this scenario, there are two sequential mergers
taking place in the system (see also e.g. Stegmann et al., 2022a). First, the
inner binary merges when the stars are still on the MS as a result of mass
transfer from the tertiary to the inner binary. This reduces the triple to a
binary. We assume that the merger product of the inner binary evolves further
in a classical way (as opposed to CHE). Consequently, the merger product
expands and eventually fills its Roche-lobe and transfers mass to the initial
tertiary star. The orbit shrinks due to this second phase of mass transfer and
as a result, a merging double compact object is formed. The second phase of
mass transfer is essential. Systems in which no mass transfer takes place
after the inner binary merger might form detached BH-BH binaries but are too
wide to merge due to GWs within the Hubble time. We note that double MS
mergers among CHE triples typically occur due to TMT episodes as three-body
dynamics are suppressed during the MS phase.

  * •

Dynamical mergers: In the triples of this channel, ZLK oscillations are very
efficient and drive up the inner eccentricities to $e_{\rm in}\approx
0.6$-$0.9$ after the stars of the inner binary have become BHs. Such systems
merge due to GW emission within a few Myrs. The tertiary remains detached
until the inner binary merges and therefore these triples belong to the no
post-MS MT channel. As discussed in section 3.3, these systems are rare.

Table 3: Summary of the statistics of GW sources from the population with
triples (with $f_{\rm single}=0.06$, $f_{\rm binary}=0.21$, $f_{\rm
triple}=0.73$), and from the population without triples (with $f_{\rm
single}=0.3$, $f_{\rm binary}=0.7$, $f_{\rm triple}=0$). The ’of all CHE
systems’ is the number of systems, expressed as a fraction of all systems
containing a binary that contains two, tidally-locked CHE stars. Formation
efficiency gives the number of systems expressed as a faction of all stellar
systems (see equation 31). Merger rate density is the merger rate density in
the local universe (see equation 45).

|  | of all CHE systems
---
at Z = 0.005 [%]
| of all CHE systems
---
at Z = 0.0005 [%]
| Formation efficiency
---
at Z = 0.005
| Formation efficiency
---
at Z = 0.0005
| Merger rate density
---
[$\rm{Gpc^{-3}yr^{-1}}$]
Population with triples
CHE triple channels: | 29.1/31.9 | 23.7/23.7 | $\boldsymbol{6.9\cdot 10^{-7}$/$7.2\cdot 10^{-7}}$ | 9.2 | 12.7/11.8
\- Effectively isolated inner binary | 19.3 | 12 | $4.6\cdot 10^{-7}$ | $3.0\cdot 10^{-7}$ | 8.8
\- TMT with BBH & BA (sc. 1/sc. 2) | 3.8/6.6 | 2.3/2.3 |  $9\cdot 10^{-8}$/$1.2\cdot 10^{-7}$ |  $5.8\cdot 10^{-8}$/ $5.8\cdot 10^{-8}$ | 1.4/0.5
\- TMT with BBH & CBD | 5.9 | 7.8 | $1.4\cdot 10^{-7}$ | $1.9\cdot 10^{-7}$ | 2.4
\- TMT with MS-MS channel | 0.3 | 5 | $6\cdot 10^{-9}$ | $1.2\cdot 10^{-7}$ | 0.2
\- Dynamical mergers | 0.2 | 0.1 | $4\cdot 10^{-9}$ | $2.0\cdot 10^{-9}$ | 0.05
CHE binaries | 8.7 | 12.9 | $\boldsymbol{7.2\cdot 10^{-7}}$ | $\boldsymbol{5.2\cdot 10^{-7}}$ | 11
Population without triples
CHE binaries | 65 | 88.6 | $\boldsymbol{1.2\cdot 10^{-6}}$ | $\boldsymbol{1.2\cdot 10^{-6}}$ | 44.2

Figure 10: Orbital parameters of GW sources that evolve through the TMT
channel with BH-BH accretors and ballistic accretion at Z=0.005. In the upper
panel we show the inner pericenter of the inner BH-BH binary ($a_{\rm p,in}$).
The black line represents $a_{\rm p,in}$ at the onset of TMT. The coloured
lines show the orbital characteristics at the end of the mass transfer for
scenario 2, where the change in the orbit has been estimated using an energy
formalism with three different values of $\alpha_{\rm TMT}\lambda_{\rm TMT}$
(equation 26 and 27). Lower panel: we show the frequency of the GW radiated by
the inner BH-BH binary before and after the tertiary mass transfer phase.
Figure 11: The pericenter of the outer orbit before and at the onset of TMT
for MS-MS inner (upper panel) and BH-BH inner binary accretors (lower panel)
as calculated with equation 30 at Z = 0.005. The grey line shows the ratio of
the pericentre after the TMT and at the before the TMT episode (i.e. $a_{\rm
p,out,afterTMT}/a_{\rm p,out,beforeTMT}$ with the corresponding values shown
in the upper x-axis) We note that for we have normalised these distributions
to $10^{4}$ (such that the numbers on the y-axis are in the order of unity).

We ignore the possibility of GW source forming in a CHE triple through a
stellar merger that do not occur between two MS stars. Such mergers can occur
due to TMT or three-body dynamics with (i) helium star-MS binary or (ii)
double helium star binaries. We justify the omission of the first type, as
they are relatively rare. This type of merger occurs in 0.2-2 per cent of all
CHE triples depending on metallicity. For the second type, the merger product
is a helium star, it is not expected to significantly expand and it is
unlikely to ever fill its Roche-lobe. Without a phase of mass transfer that
leads to orbital shrinkage, the binary remains too wide to merge within a
Hubble time. However, if the merger remnant can accrete matter during the TMT
phase it could regain a hydrogen-rich envelope, and expand later in its
evolution. For simplicity, we neglect this scenario.

###  5.1 Rates of GW mergers

In the population without triples, the predicted merger rate density is
$R_{\rm merger}=44.2\,\rm{Gpc^{-3}yr^{-1}}$ (see Table 3). This is about a
factor of two higher than predicted by Riley et al. (2021), giving a rough
agreement given the simplicity of our rate calculation (see discussion in
Appendix B). The total merger rate density of the population containing
triples is $R_{\rm merger}$ = 23 $\rm{Gpc^{-3}yr^{-1}}$. This is about a
factor of two lower than that of the population without triples. There are two
reasons for this difference. Firstly, stellar mergers frequently occur in CHE
triples, preventing the formation of compact BH-BH binaries. While all CHE
binaries form BH-BH binaries, only about 60 (45) per cent of CHE triples form
(inner) BH-BH binaries at Z = 0.005 (Z = 0.0005).

Secondly, the number of systems formed in the population with triples is
always lower per unit stellar mass formed than in the population without
triples, as triple systems, on average, have larger total masses than binaries
and/or single stars.

In the population with triples, about half of the GW mergers originate from
formation channels involving CHE originate from triples. The role of the
tertiary is negligible for 69 per cent of GW progenitors from CHE triples. In
the remaining 31 per cent, the evolution of the inner binary is affected by
the tertiary star via TMT and/or three-body dynamics.

###  5.2 Isolated binaries

At $Z=0.005$, about 68 per cent of the CHE binary population forms a BH-BH
binary that merges within the Hubble time, while at $Z=0.0005$, all CHE
binaries merge due to GWs within the age of the universe. In our moderate
metallicity model, the delay times of these BH-BH binaries from this
population ranges from 3 to 50 Gyr (and therefore the delay time of GW sources
ranges from 3 to 13.5 Gyr). In our low metallicity model, the delay times are
considerably shorter, ranging roughly from 100 to 600 Myrs. At Z = 0.005, only
those binaries merge which were in contact during their MS phase. At Z =
0.0005, about 97 per cent of all GW progenitors were in contact during their
MS phase. Since we assume such binaries equalise in mass, we predict that the
vast majority of GW sources consist of equal mass black home binaries from
this population (in broad agreement with Marchant et al. 2016). The masses of
the merging binary black holes from this channel range from 20 to 42
$M_{\odot}$ at Z = 0.005 and 33 to 54 $M_{\odot}$ at Z = 0.0005.

###  5.3 Effectively isolated inner binaries

This is the dominant channel among CHE triples with a predicted merger rate
density of 8.8 $\rm{Gpc^{-3}yr^{-1}}$. At Z = 0.005 (Z = 0.0005), about 19
(12) per cent of all CHE systems (e.g CHE binaries and CHE triples, see
section 2.6) are expected to form GW sources via this channel. In 53 per cent
of the GW progenitors of this channel, the tertiary star becomes unbound by
the time both stars in the inner binaries form BHs. This percentage drops to
38 per cent at Z = 0.0005.

The demographics of this channel are nearly indistinguishable from the
isolated binary population. The merger efficiency of this channel, which we
define as the GW sources as a fraction of BH-BH inner binaries formed via a
certain channel, is 68 per cent. Unsurprisingly, this is the same as the
merger efficiency of the isolated CHE binary channel. Similarly to the CHE
binary case, the majority of the inner binaries of these triples were also in
contact during their MS phase and therefore this channel also produces
overwhelmingly equal mass mergers.

###  5.4 TMT with a BH-BH accretor

This is the dominant formation channel in which the evolution of the inner
binary is affected by the tertiary star. The predicted merger rate density is
$R_{\rm merger}$ = 3.8 $\rm{Gpc^{-3}yr^{-1}}$, which accounts for about 16 per
cent of all GW mergers from CHE systems. About 10 per cent of all CHE systems
form merging binary BHs via this channel.

With our simplistic models of TMT (see subsection 2.5), we predict that the
outer orbit widens as a result of the TMT episode for all triples considered
in this study. In the lower panel of Fig. 11, we show how the outer pericenter
changes after the mass transfer phase for triples experiencing TMT with a BH-
BH inner binary accretor for our moderate metallicity model (and in the lower
panel of Fig. 17 for our low metallicity model). The orbital separations widen
typically by a factor 1.5-2.

Even, if the inner orbit remains unchanged due to TMT, the outer orbit widens
so much, such that three body-dynamics become typically negligible after the
TMT episode for the majoirty of these triples. For example, at Z = 0.005, in
those TMT systems, in which ZKL oscillation are effective prior to the mass
transfer event, 70 per cent of the inner binary becomes decoupled from the
tertiary star after the TMT episode. If the evolution of the inner BH-BH inner
binary is decoupled from the tertiary, its orbital evolution is solely
determined by the emission of GWs (and therefore the coalescence time can be
determined according Peters 1964, otherwise, we use equation 18).

As noted in section 2.5, we make different assumptions about the evolution of
the inner orbit based on whether a circumbinary disc is formed during TMT. We
therefore discuss the properties of GW sources from these two subtypes
separately.

####  5.4.1 Accretion through a circumbinary disc

The predicted merger rate of this channel is 2.4 $\rm{Gpc^{-3}yr^{-1}}$. The
merger rate efficiency is just 6 per cent higher than the merger rate
efficiency from isolated binaries. The slight increase is due to the small
number of eccentric inner binaries at the onset of the mass transfer ($\sim$10
per cent of systems undergoing TMT with BH-BH accretors and circumbinary discs
have $e_{\rm in}$>0.4, see Fig. 5). The small difference is not surprising as
we have assumed here that the orbit of the inner binary does not change due to
circumbinary disc accretion. However, if circumbinary disc accretion leads to
a significant increase (decrease) in the inner period, the compact object
merger fraction decreases (increases) significantly as well. Clearly, better
models are required to understand circumbinary accretion of a BH binary from a
mass transferring tertiary star.

####  5.4.2 Ballistic accretion

The properties of these GW sources depend on how the inner binary evolves due
to TMT. If we simplistically assume that that the inner orbit does not change
(i.e. Scenario 1, see section 2.5), then the merger rate density of this
channel in the local universe is $R_{\rm merger}$ = 1.4 $\rm{Gpc^{3}yr^{-1}}$.
In this case about 3.8 (2.3) per cent of all stellar systems containing a CHE
binary form GW sources via this channel at Z = 0.005 (Z = 0.0005). The merger
efficiency of this channel is 75 per cent at Z = 0.005, which is slightly
higher than that of the CHE binary population (68 per cent). As discussed in
section 3.5, a considerable fraction of these sources have high
eccentricities, namely, 48 per cent with $e_{\rm in}\gtrsim 0.4$ at Z = 0.005
and 10 per cent at Z = 0.0005. This results in shorter delay times and more
mergers with respect to the isolated CHE binary channel (top left panel of
Fig. 12).

If the orbital evolution can be described by equation 28 (i.e. Scenario 2, see
section 2.5), then the inner pericenters of BH-BH binaries decrease by 1-3
orders of magnitude due to the TMT episode, depending on the efficiency
parameter $\alpha_{\rm TMT}$. In this case, all inner bineries become
dynamically decoupled from the tertiary star after the TMT episode. As shown
in the left panel of Fig. 10, the peak of the orbital separation distribution
shifts from $32\,R_{\odot}$ to 25, 5 and 1 $R_{\odot}$ with $\alpha_{\rm
TMT}=0.05,5,5$. With such short periods, nearly all (i.e. typically $\gtrsim$
99 per cent) of the inner binaries eventually emerge. However, none of the
inner binaries merge during the mass transfer, in fact they merge due to GW
emission afterwards. In Fig. 12, we show that the typical delay times in
Scenario 2 are also orders of magnitude shorter with respect to that of
isolated CHE binaries. With $\alpha_{\rm TMT}=0.05$, the delay times of these
GW sources is dominated by the stellar evolution. Such timescales could make
TMT episodes relevant in young clusters in which star-formation is still
active. Even when assuming a weaker friction exerted by the transferred mass
(i.e. $\alpha_{\rm{TMT}}\lambda_{\rm TMT}=5$) resulting in the smallest
orbital shrinkage in our models, most of the BHs merge within a few hundred
Myr at Z = 0.005.

Despite the higher merger efficiency, the predicted merger rate density for
Scenario 2 is considerably lower (i.e. $R_{\rm
merger}=0.5\,\rm{Gpc^{-3}yr^{-1}}$) than in Scenario 1. This is due to the
extremely short delay times, implying the progenitor stars must have formed
recently, when the cosmic star formation rate is low (e.g. Madau & Fragos,
2017). As the cosmic star formation rate is expected to increase strongly from
$z=0$ to $z=2$ , we expect the merger rate density of this channel to be
significantly higher at $z\approx 2$ than at $z=0$. This would make these
sources more relevant for third-generation GW detectors.

We mention two interesting aspects of this channel. Firstly, depending on the
efficiency parameter of the TMT episode, these systems could be in the LISA
frequency band (Amaro-Seoane et al., 2022) during the mass transfer phase. In
the right panel panel of Fig. 10, we show the frequency at which the BH-BH
binaries emit GWs after the mass transfer episode. With
$\alpha_{\rm{TMT}}=0.5$, about half, and with $\alpha_{\rm{TMT}}=0.05$, all of
our systems enter the mHZ regime during the mass transfer phase. The evolution
through the LISA frequency range would be primarily driven by gas dynamics
instead of GW emission (see also Renzo et al., 2021). Such sources would be
detectable by LISA, if the corresponding luminosity distances are not larger
than $\sim 10\,\rm{kpc}$ and $\sim 10^{4}\,\rm{kpc}$ in case of
$\alpha_{\rm{TMT}}=0.5$ and $\alpha_{\rm{TMT}}=0.05$, respectively (see e.g.
Fig. 1 in Amaro-Seoane et al., 2022).

Secondly, a TMT episode could be accompanied by a detectable electromagnetic
signal, as the transferred mass is expected to heat up when it reaches the
inner BH binary. If the delay time between this signal and the GW merger is
within the lifetimes of typical observing missions, then the GW merger could
be associated with this electromagnetic counterpart (see also e.g. de Mink &
King, 2017). We find that the time between the end of the TMT episode and the
GW merger in case of $\alpha_{\rm TMT}\lambda_{\rm TMT}=0.05$ is shorter than
a year for 6 per cent of these sources at Z = 0.0005. This implies that in
this case a electromagnetic counterpart could be detected, shortly before the
GW merger. This is in contrast with the possible electromagnetic signatures
associated with BH mergers in AGN discs, where the electromagnetic counterpart
would occur after the GW merger (see e.g. McKernan et al., 2019)

###  5.5 TMT with a MS-MS accretor

This channel has a low merger rate density of $R_{\rm merger}$ = 0.2
$\rm{Gpc^{-3}yr^{-1}}$. Even though 25 per cent of all systems containing a
CHE binary experience a double MS merger in the inner binary at Z = 0.005,
only 1.1 per cent of them form merging binary BHs. This low merging efficiency
is due to two reasons. Firstly, if the mass transfer episode between the
merger product and the tertiary star proceeds in a dynamically unstable way,
the process mostly ends in stellar merger and no double compact binary is
formed. Secondly, if the same mass transfer proceeds instead in a stable way,
the binary BH typically has too wide orbit to merge within the Hubble time. We
note, however, that these predictions are sensitively dependent on uncertain
stellar physics (such as the efficiency of CEE phase, mass-loss radius
exponent and binding energy of stars with $M_{\rm ZAMS}\gtrsim
100\,M_{\odot}$). We also note that the merger efficiency is significantly
higher in our low metallicity model, 12.3 per cent of triples with double MS
merger forms merging binary BHs. As the merger efficiency seems to increase
with decreasing metallicity, and we only calculate the merger rate density
based on two metallicities, it is likely that we underestimate the merger rate
density for this channel (see more detailed explanation in Appendix section
B).

In case of a TMT episode with a MS-MS accretor, we always assume that the
inner binary merges due the mass transfer phase. We justify this assumption by
the fact that that CHE MS-MS binaries tend to be on very close orbits
($\sim$20-30 $R_{\odot}$) compared to their stellar radii
($\sim$5-10$\,R_{\odot}$). A significant fraction of them are already in
contact. Furthermore, these stars may swell up as a result of accretion, this
type mass transfer event is likely to end in merger (e.g. Braudo et al., 2022;
Leigh et al., 2020).

The merger product is a rejuvenated MS star with a mass of
$M_{1+2}=M_{1}+M_{2}$. This means that we neglect any accretion during TMT and
we assume a fully conservative merger without mass outflows. At Z = 0.005, the
mass of the inner binary merger remnant $M_{\rm 1+2}$ ranges from 65 to 188
$M_{\odot}$. The distribution has a peak around $\sim 100\,M_{\odot}$. At Z =
0.0005, the mass of the merger product ranges from 70 $M_{\odot}$ to 190
$M_{\odot}$.

The orbital separations after the TMT episode are shown in the upper panel of
Fig. 11 (and Fig. 17 for our low metallicity model). We can see that the outer
orbit widens typically by a factor of 1.7-2.5 and the orbital separations
range from 150 to 6800 $R_{\odot}$. While the ranges are similar at both
metallicities, at Z = 0.0005, the typical orbital separations are
significantly shorter.

Most of the systems experience a second phase of mass transfer after the TMT
episode (62 per cent at Z = 0.005 and 96 per cent at Z = 0.0005) and typically
the donor star is on the Hertzsprung gap during this second phase of mass
transfer ( about 99 per cent at Z = 0.005, and about 86 per cent at Z =
0.0005). More evolved donor stars are not expected to occur frequently, as the
onset of CHeB occurs at a cooler effective temperature with increasing mass
with and followed by a less significant subsequent radial expansion (Hurley et
al., 2000). In particular for $M_{\rm ZAMS}\gtrsim 100\,M_{\odot}$, stars are
predicted to expand negligibly after the CHeB, even at low metallicities.

Regarding the stability of the mass transfer between the merger remnant and
the initial tertiary, we find that it occurs in an dynamically unstable manner
in 66 (30) per cent of cases at Z = 0.005 (Z = 0.0005). We assume that CE
phases with a donor star on the Hertzsprung gap result in a merger, following
Dominik et al. (2012) (but see also Klencki et al., 2020; Marchant et al.,
2021). At both metallicities, binary BHs are only produced when the second
phase of mass transfer proceeds in a stable manner. Furthermore, in order to
form a GW source, the orbit needs to be compact enough ($a_{\rm out}\lesssim
1000\,R_{\odot}$) at the onset of the second mass transfer event. This only
occurs in about 5 per cent (30 per cent) of systems with stable mass transfer
at Z = 0.005 (Z = 0.0005).

Figure 12: The time delay distribution of the GW sources of CHE triple
population shown by a stacked histogram. As a comparison, we also show the
time delay distribution of the GW sources from CHE isolated binaries at the
same metallicity (gold). The different colours of the stacked histograms refer
to different evolutionary paths of the triples (shown by the legend in the top
of the panel)

This is the only GW formation channel of CHE triples that yields a
significantly different mass and mass ratio distributions than the CHE binary
channel. The masses of the merging binary BHs range from 16 to 27 $M_{\odot}$
at Z = 0.005 and 17 to 54 $M_{\odot}$ at Z = 0.0005. The mass ratios range
from 0.7 to 0.8 at Z = 0.005 and 0.5 to 1.0 at Z = 0.0005. All other channels
produce merging binary BHs with masses that range from 20 to 42 at Z = 0.005
and 33 to 54 at Z = 0.0005. The vast majority ($\gtrsim$ 90 per cent) of these
systems have equal masses, as the inner binaries had been in a contact during
their MS phase.

###  5.6 Dynamical mergers

The merger rate density of these channel is very low, $R_{\rm merger}$ = 0.05
$\rm{Gpc^{-3}yr^{-1}}$. The delay times of these systems are very short and
range from 4 to 20 Myrs. Similarly to the GW progenitors that have experienced
TMT episodes with ballistic accretion, the short delay times imply that the
merger rate density could be about an order of magnitude larger at z $\approx$
2\. About 25 per cent of these systems have eccentricities $e_{\rm in}\gtrsim
10^{-4}$ when the characteristic GW frequency reaches 10 Hz, making
eccentricities detectable by third-generation detectors (Lower et al., 2018).

For all systems, the tertiary star is still on the MS when the inner binary
merges due to GWs with outer pericenters of $a_{\rm{p,out}}\approx$
120-790$\,R_{\odot}$. It is therefore expected that the initial tertiary star
will eventaully fill its Roche-lobe, once it evolves of the MS.

##  6 Conclusion

We studied the evolution of hierarchical triples with CHE stars in the inner
binary with a rapid populations synthesis approach. We performed simulations
with the triple population synthesis code TRES at two representative
metallicities: Z = 0.005 and Z = 0.0005. We showed that the evolution of CHE
stars can be altered by the presence of a tertiary star in several ways. This
can potentially lead to a formation of a number of diverse and unique
astrophysical phenomena, e.g. TMT phases with BH-BH accretors, highly
eccentric mergers of helium stars, and mergers of binary BHs with very short
(few Myrs) delay times.

To summarise our main findings:

  1. 1.

Tertiary mass transfer (TMT) episodes are common among CHE triples: Unlike in
classically evolving hierarchical triples, we predict that TMT phase is very
common among CHE triples. The tertiary star fills its Roche-lobe in about 50
per cent of all triples with CHE inner binaries. The same fraction for
classically evolving systems is predicted to be a few percent at best (see
e.g. de Vries et al., 2014; Toonen et al., 2020; Hamers et al., 2022a; Kummer
et al., 2023). We find that the mass transfer episodes initiated by the
tertiary star typically occurs in a dynamically stable way.

  2. 2.

BH-BH inner binaries that accrete from tertiary star are also common: About 31
(24) per cent of the tertiary-driven mass-transfer episodes occur with BH-BH
accretors at Z = 0.005 (Z = 0.0005). Previous population synthesis studies
suggest that such scenario is probably not possible triples with classically
evolving stars (see e.g. Toonen et al., 2020; Hamers et al., 2022a).
Therefore, mass transfer towards a BH-BH inner binary represents a unique
scenario for triples (or higher-order multiples) with CHE stars in the inner
binaries. An exciting prospects would be a possible EM counterpart from such
an event (e.g. de Mink & King, 2017).

  3. 3.

Importance of three-body dynamics: ZLK oscillations can be effective for CHE
triples, if the stars in the inner binary have evolved off MS (otherwise
precession due to strong tides quench ZLK cycles) and if the initial outer
pericenter is $a_{\rm p,outer,ZAMS}\lesssim 2000\,R_{\odot}$ (otherwise ZLK
cycles are quenched by various short range forces throughout the entire
evolution of the inner binary). ZLK oscillations are only present in those CHE
triples, in which the outer pericenter is so short, such that the tertiary
star would eventually fill its Roche-lobe. The inner eccentricities of these
systems can reach values up to $e_{\rm in,max}\sim 0.9$ (left panel of Fig.
8). The effects of three-body dynamics are negligible for those CHE triples in
which the triple remains detached. In this case, the inner binary evolves
effectively as an isolated binary.

  4. 4.

Three-body dynamics can drive the inner binary to a stellar merger: In about 3
per cent of CHE triples, the inner binary merges before BH-BH formation. The
most common type is a merger of a double helium star binary, that comes into
contact in a highly eccentric orbit (Table 2).

  5. 5.

CHE triples form GW sources efficiently: About 30 (24) per cent of the CHE
triple population forms BH binaries that merge due to GWs within Hubble time
at Z = 0.005 (Z = 0.0005). We predict a merger rate density of GW sources from
CHE triples of $R_{\rm merger}\approx 12\,\rm{Gpc^{-3}yr^{-1}}$ (Table 3). We
also predict that about half of the GW sources from CHE systems originate from
triples. In 69 per cent of all GW sources from CHE triples, the inner binary
evolves effectively as an isolated binary and therefore its properties are
indistinguishable from those of CHE binaries. In the remaining 31 per cent,
the evolution of the GW progenitor is affected by three-body dynamics and/or
TMT episodes.

  6. 6.

Tertiary mass transfer and three-body dynamics could lead to the formation of
BH-BH binaries that merge within Myrs

The vast majority of those GW progenitors of CHE triples, in which the
evolution of the inner binary is not decoupled from the tertiary object,
experience a TMT episode with a BH-BH inner binary. In this case, we model the
evolution of the inner binary during the TMT phase with energy arguments
(following de Vries et al., 2014, see also subsection 2.5) and with different
assumptions on how efficiently the transferred mass shrinks the orbit of the
inner binary. We find typical values for the delay time of these GW sources of
few hundred Myrs and few Myrs in our model variation with the least and the
most orbital shrinkage, respectively.

## Acknowledgements

SdM acknowledges Fabio Antonini, Adrian Hamers and Lieke van Son for
insightful discussions. AD acknowledges travel grant from the HPC3 Europa
programme for providing computational resources at the Snelius supercomputer
in the Netherlands and acknowledges support fro API for allowing an extended
visit. Computational work was performed by the Snelius supercomputer in the
Netherlands and by the University of Birmingham’s BlueBEAR HPC service. ST
acknowledges support from the Netherlands Research Council NWO (VENI
639.041.645 and VIDI 203.061 grants). SdM acknowledges funding by the
Netherlands Organization for Scientific Research (NWO) as part of the Vidi
research program BinWaves with project number 639.042.728.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2019a) Abbott B. P., et al., 2019a, Physical Review X, 9, 031040
  * Abbott et al. (2019b) Abbott B. P., et al., 2019b, The Astrophysical Journal, 882, L24
  * Abbott et al. (2021) Abbott R., et al., 2021, Physical Review X, 11, 021053
  * Amaro-Seoane et al. (2022) Amaro-Seoane P., et al., 2022, Astrophysics with the Laser Interferometer Space Antenna, doi:10.48550/ARXIV.2203.06016, https://arxiv.org/abs/2203.06016
  * Anderson et al. (2017) Anderson K. R., Lai D., Storch N. I., 2017, MNRAS, 467, 3066
  * Antognini et al. (2014) Antognini J. M., Shappee B. J., Thompson T. A., Amaro-Seoane P., 2014, MNRAS, 439, 1079
  * Antoni et al. (2019) Antoni A., MacLeod M., Ramirez-Ruiz E., 2019, ApJ, 884, 22
  * Antonini et al. (2017) Antonini F., Toonen S., Hamers A. S., 2017, ApJ, 841, 77
  * Barkat et al. (1967) Barkat Z., Rakavy G., Sack N., 1967, Phys. Rev. Lett., 18, 379
  * Bartos et al. (2017) Bartos I., Kocsis B., Haiman Z., Márka S., 2017, ApJ, 835, 165