Stellar triples with chemically homogeneously evolving inner binaries

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: 𝑍 = 0 . 005 and 𝑍 = 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 ∼ 1200 𝑅 ⊙ . 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.


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.
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;Abbott et al. 2019aAbbott et al. , 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 rapidlyrotating stars (Maeder 1987;Yoon & Langer 2005;Yoon et al. 2006;Brott et al. 2011;Köhler et al. 2015;Szécsi et al. 2015), 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 buildup of a chemical gradient inside the star and the classical coreenvelope 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;Hastings 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.Finally, in section 5, we use analytical and numerical methods to explore our synthetic population of stellar triples in the context of GW sources.

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), as implemented in the 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 (see e.g.Veras et al. 2011;Debes & Sigurdsson 2002) with mass loss rates 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: =  in,GR +  in,TF +  in,wind  out =  out,GR +  out,TF +  out,wind  in =  in,3b +  in,GR +  in,TF  out =  out,3b +  out,GR +  out,TF  in =  in,3b +  in,GR +  in,tides +  in,rotate  out =  out,3b +  out,GR +  out,tides +  out,rotate ℎ in = ℎ in,3b where , , , ℎ and   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  ≡ cos(), where  is the mutual inclination between the inner and outer orbit, and Ω 1 , Ω 2 , Ω 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 (  in,tides ) and GR effects (  in,GR ) occurs on a much shorter timescale than the precession due to threebody dynamics ( 𝑔 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): The timescale related to the apsidal precession due to tides is (e.g.Smeyers & Willems 2001;Liu et al. 2015): where  am the apsidal motion constant, which we assume to be 0.0144 for MS and helium stars,  in =  ( 1 +  2 ), i.e. the standard gravitational parameter for the inner binary and  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): If  ZLK ≫ min(t GR , t 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 ( 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),  tides ≫  rot and therefore precession due to stellar rotation does not play a role in suppressing three-body dynamics (Liu et al. 2015).
Secondly, the octupole term in the perturbing function of the Hamiltonian (e.g.Naoz et al. 2013) is typically negligible for CHE triples, as the vast majority of the inner binaries are in contact, which leads to equal mass components in our models (see section 2.2).Furthermore, for the relatively rare detached CHE inner binaries, the mass ratio is always  in ≥ 0.7 (see section 2.6).
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): where  GW is the time required for the merger,  GW,Peters is the time to merger based on the relation of Peters (1964),  in,max is the maximum eccentricity reached during ZLK oscillations and  in is the initial inner semimajor axis.The approximation in equation 5 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 5 underestimates the actual merger timescale typically by a factor of 2-3.

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.  star >  CHE,crit .Riley et al. (2021) provides a fit to this critical value based on MESA (Paxton et al. 2011) models of Marchant et al. (2016) at different masses and metallicities .In order to determine whether a star evolves chemically homogeneously, we check whether our simulated star is spinning above  CHE,crit at every timestep.If a star meets this criteria, we determine its radius and luminosity according to the fits of Hurley et al. (2000) for ZAMS stars.We note that the mass of CHE stars are in general affected by stellar winds, therefore, its radius and luminosity do not remain constant during the core-hydrogen burning phase, even with this simplifying assumption (see also Riley et al. 2021).
We assume that the star by the end of core hydrogen burning forms a helium star with a mass  He,ZAMS =  TAMS , where  He,ZAMS is the initial mass of the helium star and  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 stars, 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.

Contact binaries
We follow the implementation of Riley et al. (2021) of modelling contact binaries for rapid population synthesis codes (which is based on the detailed stellar models of Marchant et al. 2016).We assume that contact binaries, i.e. binaries in which both stars fill their Rochelobes, 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.If this mass equalisation indeed occurs in nature, the mass ratio distribution of massive contact binaries would exhibit a prominent peak near one, which is not in an agreement with the observations of massive contact binaries residing in the Magellanic Clouds and in the Milky Way (see Menon et al. 2021;Abdul-Masih et al. 2022).This discrepancy could be due to missing physics in the contact binary models of Marchant et al. (2016), e.g.related to energy transfer between the two stars (see e.g.Abdul-Masih et al. 2021;Fabry et al. 2023).
We follow Marchant et al. (2016) and approximate the L2 point as where  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): where  init ,  fin are the initial and the final orbital separation and  1,init ,  2,init are the initial masses of the primary and the secondary, respectively.The final masses are 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).

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 1).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).
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.With these assumptions, the inner and the outer orbit changes due to stellar winds as: and 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(Huang , 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).We note, however, that the assumptions of fast and spherically symmetric wind might not always be valid for short period binaries (e.g.Brookshaw & Tavani 1993), furthermore, rapidly rotating stars might not have fully symmetric outflows (Georgy et al. 2011).Particularly, stellar winds in certain binaryconfigurations might even lead to orbital shrinking (Schrøder et al. 2021).

Remnant formation
The mass of the compact object remnant is computed based on the delayed supernova model from Fryer et al. (2012).This prescription determines the mass of the stellar remnant as a function of CO core mass at the onset of the core-collapse.The latter is determined in SeBa based on the fits of Hurley et al. (2000).The natal kick velocity for BHs is calculated as where  b is the fallback fraction (Fryer et al. 2012),  NS is the canonical neutron star mass ( NS = 1.4 ⊙ ) and  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 Pĳloo et al. (2012).Models of Fryer et al. (2012) predict that the most massive stars collapse directly (typically  ZAMS ≳ 40  ⊙ ), 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-corecollapse mass of the star.We note that the actual neutrino mass loss is considered to be uncertain, and other population synthesis codes assume considerably smaller losses, e.g. 1 per cent of pre-collapse mass (Belczynski et al. 2020), or a fixed value of 0.1  ⊙ loss (Riley et al. 2022).Recent observations suggest that mass decrements from neutrino emission could indeed be only about a few percent of the pre-collapse mass (Vigna-Gómez et al. 2023).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  core ∼ 3 × 10 9 .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. 2019;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. 2019;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-corecollapse is  HE,pre−SN ≥ 35  ⊙ , 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 ≤  HE,pre−SN ≤ 130  ⊙ , we assume the star undergoes PISN, and leaves no remnant behind.In principle, if  HE,pre−SN ≥ 130  ⊙ , photo-disintegration prevents the pair instability supernova 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,since in our simulation the maximum stellar mass is limited to 100  ⊙ .We note that the above-quoted mass ranges are sensitively dependent on the the poorly constrained reaction rate of 12 C(, ) 16 O (see e.g.Takahashi et al. 2018;Farmer et al. 2019;Woosley & Heger 2021;Costa et al. 2021).In particular, more recent simulations, with updated 12 C(, ) 16 O rates (see e.g.deBoer et al. 2017), find that the lower edge of the BH mass gap is located at a considerably higher value ( ≈ 60  ⊙ ) than previously determined (see e.g.Mehta et al. 2022;Farag et al. 2022).

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. out ≫  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  crit =  donor / 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.  >  crit .We obtain  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  crit can be found in Toonen et al. 2016).These are  crit = 3 and  crit = (1.37+2[donor,core / 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 ∼ 0.45 -0.65, which is fairly typical for massive CHeB stars (Dorozsmai & Toonen 2022), then  crit ≈ 0.7-0.75.This reflects the assumption made by Hurley et al. (2000), that CHeB stars tend to have deep convective envelopes (cf.Klencki et al. 2020), and are therefore more likely to experience unstable mass transfer episodes (see e.g.Hjellming & Webbink 1987, but see Woods & Ivanova 2011).
Stable TMT could be accompanied with the formation of a circumbinary disc or it could occur in a ballistic accretion fashion.These two types of mass transfer phases could lead to significantly different evolution of the inner orbit (de Vries et al. 2014).We assume that TMT occurs via ballistic accretion, if  in (1 +  in ) ≥  cd at the onset of the TMT phase, where  cd is: where we adopted the fitting formulas for mass transferring binaries of Lubow &Shu 1975 andUlrich &Burger 1976 to triples.

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: (i) an inner binary with compact objects and with ballistic accretion, (ii) an inner binary with compact objects and with a circumbinary disc, (iii) 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 -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.Δ trnsf is the mass that is transferred from the tertiary in a timestep Δ.When Δ trnsf ends up encompassing the inner binary, it has binding energy of  bind .As the inner orbit is shrinking due to the friction during the TMT episode, the orbital energy of the inner binary changes by Δ orb .We assume that a fraction ( TMT ) of Δ orb is used to unbind Δ trnsf .We can write an equation expressing the energy balance as: with and where  TMT is a parameter related to the structure of Δ trnsf , parameterising its binding energy,  in,init is the initial orbital separation before Δ trnsf is transferred to the inner binary and  in,fin is the final orbital separation after Δ 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  out,env (but see Laplace et al. 2020).Then assuming a constant  TMT and  TMT , the orbit changes due to the entire TMT episode as: As both  TMT and  TMT are unknown, we combine them and try three different values:  TMT  TMT = 0.05, 0.5, 5.Here  TMT  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  in,new =  in (1 −  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  ≳ 0.95 might retain eccentricities as high as  ∼ 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).

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: where  is the fraction of mass accreted by the inner binary,  is the specific angular momentum lost from the system as a fraction of the specific angular momentum of the triple and  3 is the mass transfer rate from the tertiary star.Equation 16 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  1 +  2 .We assume that eventually all the transferred mass is isotropically expelled from the triple ( = 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 ( =  3 /( 1 +  2 ), see also Hamers et al. 2021, for a similar approach).In this case equation 16 can be expressed as 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  = 0 might no longer be a good approximation.

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 Gyr) 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 ( 1,ZAMS ) follows the power-law mass distribution of Kroupa (2001) for  ZAMS < 0.5  ⊙ .We sample  1,ZAMS from a mass range of 20-100  ⊙ .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. in,ZAMS =  2,ZAMS / 1,ZAMS ) distribution, which is in broad agreement with Sana et al. (2012).We restrict the range of  in,ZAMS to 0.7-1 given that inner binaries in which both of the stars are chemically homogeneously evolving and have  in ≤ 0.7 are in contact and merge early during the MS due to outer Langrange overflow (where we found the lower limit of 0.7 from our simulations).We sample the inner semimajor axis from a loguniform distribution (Öpik 1924; and in broad agreement with Sana et al. 2012) in the range of 16 to 40  ⊙ .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  ⊙ .We assume that the distribution of the outer mass ratio (i.e. out,ZAMS =  out,ZAMS /( 1,ZAMS +  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  ⊙ .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 cos(i ZAMS ), where i ZAMS is the initial inclination.The initial argument of the pericenter is assumed to be uniformly distributed between − and .
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  single = 0.06,  binary = 0.21,  triple = 0.73, see e.g.Moe & Di Stefano 2017), whereas in the second test population, we assume there are no triples and about 70 per cent of massive stars are in binaries (  single = 0.3,  binary = 0.7,  triple = 0, see e.g.Sana et al. 2012).Regarding our first test population, we note that Moe & Di Stefano (2017) finds that 73 per cent of O stars are either in triples or quadrupoles.Therefore  triple = 0.73 should be considered as a rough upper limit.However, 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.We also note that tertiary companions have been detected for many massive contact binary systems (see e.g.Kennedy et al. 2010;Mayer et al. 2013;Lorenzo et al. 2017;Janssens et al. 2021).Regarding, the second test population, we note that 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.

RESULTS OF POPULATION SYNTHESIS SIMULATIONS
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 have an inner binary where both stars evolve chemically homogeneously from ZAMS (CHE at ZAMS triples, see Table 1).We follow the further evolution only for these systems.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  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 can counteract both of these processes, inner binary merger occurs more frequently at low metallicities (i.e. 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).Computational issues arise more often for systems that are close to the dynamical instability (i.e. have low  out / in ).Therefore, our estimates for the occurrence rates of systems, in which stellar merger due to ZLK oscillations or dynamical instability occurs should be considered as a lower limit (see next subsection).

Main evolutionary outcomes
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 such systems.
• 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.

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:  1,ZAMS =  2,ZAMS = 70  ⊙ ,  in,ZAMS = 22.4  ⊙ .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 Gyr).The inner binary is tidally locked and therefore  in,ZAMS = 0.The stars of the inner binary are in contact from ZAMS.The initial mutual inclination is i ZAMS = 90 • 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 ( out ) and outer semi major axes ( 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  tides ∼  −5 , see equation 3) 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  out ≲ 500  ⊙ .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 for CHE triples, if the tertiary star is the initially least massive star in the system), the timescale related to ZLK oscillations will typically decrease.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  out ≲ 2000  ⊙ will develop ZLK oscillations.

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  out,ZAMS = 32.1  ⊙ and a circular outer orbit with  out,ZAMS = 200  ⊙ .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 Myr, 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  in ≈ 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 stop the simulation at this point and assume it leads to a merger.In principle it is possible that the system would enter a stable contact phase, however given the eccentricities, the majority of these systems would probably soon experience outer Langrange overflow, before tides would quench ZLKs and before the inner stars would form BHs.More detailed models are necessary to understand the evolution of a double helium star system in contact, in particular when undergoing three-body dynamical effects.

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:  out,ZAMS = 421  ⊙ .
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 eventually effective, however, because 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,  GR increases by a factor of 5, while  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  in,max = 0.75.Since core-collapse occurs in an eccentric orbit, large range of possbile post-supernova orbits are possible ( in = 42-186  ⊙ ) 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  in = 171  ⊙ 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 ( in = 0.94) BH-BH binary at a highly inclined orbit (i = 71.5 • ).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.

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  out,ZAMS = 1069  ⊙ .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  out,ZAMS = 1720  ⊙ and the tertiary mass is  out = 31.9 ⊙ .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  in = 46.6  ⊙ and  out = 1860  ⊙ and therefore the ZLK oscillations remain quenched.At 6 Myr, the tertiary reaches a radius of 547  ⊙ 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.

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  = 0.0005.However at  = 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. p,out,ZAMS ≳ 2000-3000  ⊙ 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,  p,out,ZAMS ≈ 100-700  ⊙ ), and secondly, the tertiary has a relatively low mass (typi-  cally  out,ZAMS = 10-30  ⊙ ).In these systems, the ZLK oscillations drive the eccentricity of the inner BH-BH binary up to large values (e.g. in ≳ 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).

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 2).
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  ⊙ 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  ⊙ and it does not exceeds 700  ⊙ .The eccentricities of the inner binary at the moment of the merger typically have values of  in ≈ 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.

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  = 0.005 ( = 0.0005) of all CHE triples follow this type of evolutionary path (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  binary = 0.21,  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.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.

Donors of TMT episodes
Here, we discuss the different stellar evolutionary stages of the donor stars at the onset of the mass transfer phase, as it is highly relevant for determining the stability of the mass transfer episode (which has dramatic effect on the outcome of the TMT, compare e.g. de Vries et al. 2014;Glanz & Perets 2021b).
In particular, core-helium-burning or asymptotic giant branch stars tend to develop deep convective envelopes.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).
At  = 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 massive stars 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  out,ZAMS ≲ 8  ⊙ ) or stars on the asymptotic giant branch.
At lower metallicities, CHeB donors are more prevalent.At  = 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.

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  out <  crit , see right panel of Fig. 4 for our moderate metallicity model, and Fig. A1 for our low metallicity model).Typical mass ratios for systems with HG donors are  out = 0.4-0.8,while for CHeB donors, they are  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.

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 of the inner binary, 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 episode is is less likely to lead to merger by itself, since the orbit has to shrink by a much larger factor with respect to double MS systems.In fact, in our models, this process never leads directly to coalescence (see later discussion in section 5.4, but also see models, in which compact binary mergers are predicted in gaseous environemnts Antoni et al. 2019;Rozner & Perets 2022;Lai & Muñoz 2022;Siwek et al. 2023;D'Orazio & Duffell 2021).Nevertheless, a TMT epsiode with a BH-BH inner binary could be a source of (an observable) X-ray emission (e.g.Lewin et al. 1997).
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 orbital separation 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.

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 (i.e. transferred / tot,inner ) 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  ⊙ for BH-BH accretors and between 10-50  ⊙ for MS-MS accretors (see left panel of Fig. 4 for Z = 0.005 and A1 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 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 ( transferred / tot,inner ) has the same maximum value (∼ 0.5) for both BH-BH and MS-MS accretors (see grey histogram in left panel of Fig. 4).

Formation of circumbinary disc
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.In this subsection, we discuss how common it is for TMT systems to develop a circumbinary disc at the onset of the mass transfer episode.
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 ( p,out ≈ 300-6000  ⊙ ) than in those where TMT proceeds in a ballistic manner ( p,out ≈ 100-600  ⊙ ).
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 11).
TMT episodes with inner BH-BH binaries are somewhat more likely to occur in a ballistic fashion than with MS-MS inner binaries.We show the properties of the donor star at the onset of tertiary mass transfer episode.The right panels show the outer mass ratio,  out =  out /(  1 + 2 ), while the left panels show the amount of mass transferred from the tertiary to the inner binary for systems undergoing TMT at Z = 0.005, 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 (see e.g.Table 2) .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.
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 11).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.

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  in ≳ 0.4 at this stage.This is in contrast with the systems with circumbinary discs; about 90 per cent of the systems have eccentricities  in ≲ 0.1.The difference is due to the smaller inner period to outer period ratios that systems without circumbinary discs have (see equation 2).In our low metallcity model, high eccentricites at the onset of TMT are much less common (see Fig. A2 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.

Significance of Humphreys-Davidson limit on the predicted rate of TMT among CHE triples
There are several poorly understood aspects of stellar physics that make our predictions for the occurrence rate of TMT episodes of CHE triples uncertain, such as mixing processes that trigger CHE, stellar winds, and the radial evolution of classically evolving ter- tiary stars.The latter is especially uncertain for stars that eventually reach the so-called Humphreys-Davidson (HD) limit (roughly stars with  ZAMS ≳ 40  ⊙ ).This empirically determined limit represents a boundary beyond which no cool supergiants are observed in the Milky Way and the Magellanic Clouds (see e.g.Humphreys & Davidson 1979;Davies et al. 2018;Davies & Beasor 2020).While there is still no consensus regarding the origin of the HD limit, its existence could imply that the most massive stars remain relatively compact throughout their lives and never become red supergiants or spend only a very short fraction of their CHeB lifetimes as red supergiants.If the former case is true, the maximum radius of the tertiary stars are significantly overestimated in our models, and therefore so is the occurance rate of TMT episodes among CHE triples.
Several possible explanations for the HD limit have been proposed in the recent years.One of the most common hypothesis is that strong (steady-state or eruptive) stellar winds strip the hydrogen envelopes of stars near the HD limit (e.g.Lamers & Fitzpatrick 1988;Smith 2014;Vink & Sabhahit 2023).Other studies suggest that different mixing mechanisms could be responsible for preventing the redward evolution of massive stars (see e.g.Langer & Maeder 1995).More specifically, this could be caused by enhanced energy transport in convective regions near the Eddington limit in evolved stars (Sabhahit et al. 2021), or by efficient semiconvection (Higgins & Vink 2020), or by strong convective overshooting (Gilkis et al. 2021;Schootemeĳer et al. 2019).
Below, we determine the fraction of TMT episodes that occur with donors beyond the HD limit.For this, we assume the following metallicity independent HD limit (see e.g.Hurley et al. 2000): We note that the metallicity dependence of the HD limit is still debated (see e.g.Davies et al. 2018;Davies & Beasor 2020).64 (67) per cent of the TMT episodes occur with donor stars that have already crossed the HD limit at Z = 0.05 (Z = 0.0005).This implies that, if classically evolving stars indeed never evolve beyond the HD limit, only about 20 (17) per cent of CHE triples would initiate TMT in our moderate (low) metallicity model instead of 55 (52) per cent.Furthermore, under such conditions, 32 (31) percent of the TMT episodes would occur with double MS-MS inner binaries, and 62 (60) percent with BH-BH inner binaries at Z = 0.05 (Z = 0.0005).
We do not run additional simulations to predict how CHE triples would evolve, if the radial expansion of the tertiary stars would be restricted by the HD limit.We note, however, that at the onset of the mass transfer episode, ZLK oscillations are not quenched in 22 per cent of CHE triples experiencing TMT episodes with BH-BH inner binaries and with donor stars beyond the HD limit at Z = 0.005.The typical ZLK timescales of these systems range from few tens to a few hundreds of Myr.This suggests, that in these systems, the tertiary star could still affect the evolution of the inner binary, even if the radial expansion of the tertiary star was restricted by the HD limit and therefore the TMT episode did not take place.Specifically, ZLK oscillations could lead to the formation of GW sources with considerably shorter delay times than those formed from isolated CHE binaries.
Whether this would really occur also depends on the exact mechanism responsible for the HD limit.In particular, if strong, steadystate stellar winds are responsible, the mass loss rates have to be  LBV ∼ 10 −3  ⊙ yr −1 , which is about an order of magnitude larger than assumed in our models 1 .As a result of such intense mass loss, the outer orbit would significantly widen soon after the tertiary star evolves off the MS.This, along with the decreased mass of the tertiary star, would lead to a significant increase of the ZLK oscillation timescale, which could potentially quench the three-body dynamics already after few Myr after ZAMS.On the other hand, if the radial expansion is prevented by various mixing processes instead of stellar winds, then the mass loss rates are comparable with that of our models.This could lead a non-negligible systems in which the tertiary remains detached but ZLK oscillations still dominate the evolution of the inner binary.

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 corecollapse, 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 1 This mass loss rate follows from the following considerations.Stellar models of Hurley et al. (2000) predict that at metallicities that are typical for the Milky Way and the Magellanic Clouds (i.e. ≳ 0.005), stars typically reach the HD limit, when they are in their very short Hertzsprung gap phase and initiate core-helium burning well beyond the HD limit, as red supergiants (i.e. eff ≲ 4800 ).Therefore, if steady-state stellar winds are responsible for preventing the formation of red supergiants for stars that cross the HD-limit, the envelope stripping has to occur in a timescale that is similar or shorter than the lifetime of the Hertzsprung gap phase (i.e. ≲ 10 4 yrs).Since the typical hydrogen envelope mass for the most massive stars is a few tens of  ⊙ , the mass loss rate has to be  LBV ∼ 10 −3  ⊙ yr −1 in this case.We have also confirmed this by varying the LBV mass loss rate of single stars and checking the corresponding stellar tracks in the Hertzsprung-Russell diagram.Similar conclusions can also be found in Mennekens & Vanbeveren (2014).
binary consists of CHE stars, they have large initial masses (i.e. ZAMS ≳ 30  ⊙ ) 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  out ≥ 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  out,ZAMS ≈ 8-25  ⊙ (see also discussion in section 4), as these systems are expected to receive the largest kicks according the supernova prescription of Fryer et al. (2012) .

Systems which become dynamically unstable
These triples typically have very short initial outer pericenters ( p,out,ZAMS ≈ 70-400  ⊙ ) and therefore are very close to the stability limit at ZAMS.Such systems can transition to non-secular or non-hierarchical evolution, if  in / out ,  out or  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./) in the inner binary is higher than that of the tertiary star,  in / out and  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 corecollapse in one of the inner stars.As noted in section 2.4, CHE stars typically form BHs via direct collapse, such that  out only increases slightly.Furthermore, the direct collapse is expected to be accompanied by a weak Blauw-kick due to neutrino losses such that  in / out and  out only increase significantly, if the inner or the outer pre-core-collapse orbits are eccentric, respectively.The pre-corecollapse inner orbit is eccentric in 72 per cent of the systems of this channel.This high inner 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 in circular inner orbits.These systems still become unstable during the BH formation, because 1) either  in / 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  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.

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  = 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  ⊙ stars and with a circular initial orbit with  in,ZAMS = 22.4  ⊙ (similarly to the example systems discussed in section 3.2).The initial tertiary mass ranges from 5 to 100  ⊙ , while  out,ZAMS ranges from 200 to 10 4  ⊙ .

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  ⊙ .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  out,ZAMS ≲ 59  ⊙ , helium stars, if 59  ⊙ ≲  out,ZAMS ≤ 70  ⊙ , and MS stars, if  out,ZAMS ≥ 70  ⊙ .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  2,ZAMS , the fraction of triples for which  2,ZAMS >  out,ZAMS increases because of our assumptions of a maximum initial stellar mass of  ZAMS,max = 100  ⊙ 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  ⊙ , while the tertiary masses range between 32 ≤  out,ZAMS / ⊙ ≤ 68.For lower tertiary masses ( out,ZAMS < 32  ⊙ ), 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 ( out,ZAMS > 70  ⊙ ), 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.
Unbound systems shown in Fig. 6 have a specific initial tertiary mass range of 8 ≲  out,ZAMS / ⊙ ≲ 25.Tertiary stars in this mass range receive strong natal kicks during remnant formation, which can lead to unbinding the triple system.Unbound systems from the entire population of CHE triples, however, are not confined to this specific range of initial tertiary masses.In fact, as it was mentioned in section 3.6, the majority of CHE triples becomes unbound due to the core-collapse occurring in the inner binary.The main reason why this does not occur in any of the triples shown in Fig. 6 is because these systems have circular outer orbits and therefore a weak Blauuw kick does not change the outer orbit significantly.
Triples of the no post-MS MT channel in the left panel of Fig. 6 have initial outer orbits  out ≳ 2000-3000  ⊙ .Their initial tertiary mass is also typically outside of the range of ∼8-25  ⊙ , 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 ( 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 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  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  outer,ZAMS at which the tertiary stars can still fill their Roche-lobes also increases with decreasing metallicity.

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 ( in,max ).About 29 per cent of the triples shown in the right panel of Fig. 6 reach  in,max ≥ 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  out,ZAMS ≲ 1200  ⊙ and  out,ZAMS ≲ 70  ⊙ .When the outer orbit is  out,ZAMS ≳ 1200  ⊙ , the ZLK cycles are quenched by various short range forces (e.g.precession caused by tides or general relativistic effects).If  out,ZAMS ≲ 1200  ⊙ but  out,ZAMS ≳ 70  ⊙ , 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  in,max does not decrease smoothly with decreasing outer orbital separations, instead it drops rather abruptly across  out,ZAMS ≈ 1200  ⊙ .Triples with  out,ZAMS ≈ 1200  ⊙ reach very large inner eccentricties ( in,max ≈ 0.9), while at slightly larger orbital separations (i.e. out,ZAMS ≈ 1500  ⊙ ) 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  p,out,ZAMS ≲ 1200  ⊙ .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  p,out,ZAMS ≲ 1200  ⊙ , a wide range of inner eccentricites are possible ( in,max = 0-0.9)for all  p,out,ZAMS .In this case, the value of  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 ( = 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. A4 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.

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  single = 0.3,  binary = 0.7,  triple = 0), in the other triples dominate (  single = 0.06,  binary = 0.21,  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).
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 threebody 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  in ≈ 0.6-0.9 after the stars of the inner binary have become BHs.Such systems merge due to GW emission within a few Myr.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.
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.Table 3. Summary of the statistics of GW sources from the population with triples (with  single = 0.06,  binary = 0.21,  triple = 0.73), and from the population without triples (with  single = 0.3,  binary = 0.7,  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 B1).Merger rate density is the merger rate density in the local universe (see equation B15

Rates of GW mergers
In the population without triples, the predicted merger rate density is  merger = 44.2Gpc −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  merger = 23 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 progen-itors 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.

Isolated binaries
At  = 0.005, about 68 per cent of the CHE binary population forms a BH-BH binary that merges within the Hubble time, while at  = 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 Myr.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  p,out,afterTMT / 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).
holes from this channel range from 20 to 42  ⊙ at Z = 0.005 and 33 to 54  ⊙ at Z = 0.0005.

Effectively isolated inner binaries
This is the dominant channel among CHE triples with a predicted merger rate density of 8.8 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.

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  merger = 3.8 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. A5 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 5).
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.

Accretion through a circumbinary disc
The predicted merger rate of this channel is 2.4 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 (∼10 per cent of systems undergoing TMT with BH-BH accretors and circumbinary discs have  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.

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  merger = 1.4 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  in ≳ 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 15 (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  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  ⊙ to 25, 5 and 1  ⊙ with  TMT  TMT = 0.05, 0.5, 5.With such short periods, nearly all (i.e.typically ≳ 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  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. TMT  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. merger = 0.5 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  = 0 to  = 2 , we expect the merger rate density of this channel to be significantly higher at  ≈ 2 than at  = 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  TMT = 0.5, about half, and with  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 ∼ 10 kpc and ∼ 10 4 kpc in case of  TMT = 0.5 and  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  TMT  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)

TMT with a MS-MS accretor
This channel has a low merger rate density of  merger = 0.2 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  ZAMS ≳ 100  ⊙ ).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 (∼20-30  ⊙ ) compared to their stellar radii (∼5-10  ⊙ ).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  1+2 =  1 +  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  1+2 ranges from 65 to 188  ⊙ .The distribution has a peak around ∼ 100  ⊙ .At Z = 0.0005, the mass of the merger product ranges from 70  ⊙ to 190  ⊙ .
The orbital separations after the TMT episode are shown in the upper panel of Fig. 11 (and Fig. A5 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  ⊙ .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  ZAMS ≳ 100  ⊙ , 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 ( out ≲ 1000  ⊙ ) 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).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  ⊙ at Z = 0.005 and 17 to 54  ⊙ 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 (≳ 90 per cent) of these systems have equal masses, as the inner binaries had been in a contact during their MS phase.

Dynamical mergers
The merger rate density of these channel is very low,  merger = 0.05 Gpc −3 yr −1 .The delay times of these systems are very short and range from 4 to 20 Myr.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 ≈ 2.About 25 per cent of these systems have eccentricities  in ≳ 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  p,out ≈ 120-790  ⊙ .It is therefore expected that the initial tertiary star will eventaully fill its Roche-lobe, once it evolves of the MS.

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 Myr) delay times.
To summarise our main findings: (i) 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.
(ii) BH-BH inner binaries that accrete from tertiary star are also common: About 31 (24) per cent of the tertiary-driven masstransfer 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 for 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).
(iii) 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  p,outer,ZAMS ≲ 2000  ⊙ (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 Rochelobe.The inner eccentricities of these systems can reach values up to  in,max ∼ 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.
(iv) 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).
(v) 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  merger ≈ 12 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.
(vi) Tertiary mass transfer and three-body dynamics could lead to the formation of BH-BH binaries that merge within Myr 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 Myr (ii) Birth rate density (equation B8) (iii) Merger rate density (equation B15) for each identified evolutionary channels.In this section, we discuss in detail how we determine these quantities.
(i) Formation efficiency: The formation efficiency expresses the number of ZAMS stellar systems formed that will evolve according to a specific evolutionary channel as a fraction of all ZAMS stellar systems formed.We calculate this quantity as: where  channel is the number simulated systems that evolves according to the channel of interest,  simulated the total number of sampled systems, and  pm is the portion of the simulated parameter space with respect to the complete parameter space, that is: where  triple is the assumed triple fraction,  M 1,ZAMS is the fraction of the simulated parameter space of primary masses: where we assumed that the absolute minimum stellar mass is  ZAMS,min = 0.08  ⊙ and the absolute maximum stellar mass is  ZAMS,max = 100  ⊙ , and as explained in section 2.6, we sample primary masses in the range of 20-100  ⊙ .The fraction of the simulated parameter space of inner mass ratios is: since the distribution of (inner and outer) mass ratios is assumed to be uniform.In equation B4, we assume that inner mass ratios of hierarchical triples have an interval of (0,1] and we sample from the interval of [0.7,1].The fraction of the simulated parameter space of outer mass ratios is  q,out = 1.0 − 0.1 1.0 − 0.0 where we assume that outer mass ratios triples have an interval of (0,1] and we sample from the interval of [0.1,1].The fraction of the simulated parameter space of inner semimajor axis is: to be uniform in a logarithmic space.We assume that inner mass semimajor axes of all triples range from 14  ⊙ to 10 5  ⊙ and we sample from the interval of [14,40]  ⊙ .Finally, the fraction of the simulated parameter space of outer semimajor axis is:  a,out = log 10 (10 5 R ⊙ ) − log 10 (10 2 R ⊙ ) log 10 (10 5 R ⊙ ) − log 10 (10 2 R ⊙ ) where we assume that inner mass semimajor axes of all triples range from 10 2  ⊙ to 10 5  ⊙ and we sample from the enitre interval.Equation B2 for channels involving isolated binaries reduces to  pm =  binary •  M 1,ZAMS •  q,in •  a,in .
ii) Birth rate density: The birth rate density gives the number density of ZAMS stellar systems in the local universe (that is at redshift  ≈ 0), which will evolve according to a specific channel.We calculate the birth rate of systems in a certain channel as: where  low and  high are 0.0015 (10 −10 ) and 0.01 (0.0015), respectively, for our model with Z = 0.005 (Z = 0.0005).Here, Z = 0.0015 is the midpoint between Z = 0.005 and Z = 0.0005 in logarithmic space, Z = 0.01 is the highest metallicity at which CHE binaries can still form GW sources at appreciable numbers and Z = 10 −10 is an arbitrarily chosen, extremely low metallicity value.In equation B9, SFRd(z) is the star formation rate density, and we use the model from Madau & Dickinson (2014): SFRd() = 0.01 • (1 + ) 2.6 1 + ((1 + )/3.2) 6.2  ⊙ yr −1 Mpc −3 , (B10) and  met (, ) is the metallicity distribution of the stellar mass formed.This quantity is also redshift dependent and assumed to follow a log-normal distribution (Madau & Fragos 2017): where  IMF is the normalised, piecewise continuous initial mass function of Kroupa (2001),  single and  binary are the single and binary fractions, respectively.We neglect higher order systems, such that  triple = 1 −  single −  binary .
We note that we also assume that the binary and triple fractions are independent on the primary mass of the system (which is clearly not consistent with observations, see e.g.Moe & Di Stefano 2017, but commonly assumed in population synthesis studies as a simplification).Assuming flat mass ratio distributions for both the inner and outer binary, equation B12 becomes: The term SFRd * (  ,  ZAMS = 0)/ M in equation B8 then gives the average number of stars formed at redshift  = 0 in a metallicity range of  i,low ≤  ≤  i,high .Multiplying this term with  formation gives the number of systems formed in a given formation channel as a fraction of all systems formed in the above mentioned metallicity range for a given star formation history model.Summing these values over all of our metallicity bins therefore yields the total birth rate of systems in a specific channel.
iii) Merger rate density: The merger rate density gives the rate density of a given astrophysical event (such as GW transients from coalescing double compact objects) in the local universe.The main difference between the birth and merger rate is due to the considerable delay time between the formation of the stellar system and the occurrence of the GW merger.For example, if the delay time for a GW source at  = 0 is  delay = 10.5 Gyr, then the redshift at ZAMS of its progenitor systems is  ZAMS ≈ 2, at which the star formation rate density is an order of magnitude higher with respect to its value at  = 0 (e.g.see models of Madau & Dickinson 2014).We determine the merger rate density at  = 0 as: where  ZAMS is the redshift at which the progenitor of a given astrophysical event is formed (and therefore it is a function of delay time), ε is the number of astrophysical events occurring at z = 0 with a delay time of  delay as a fraction of all ZAMS stellar systems formed at z =  ZAMS .We determine  ZAMS for a given delay time via the standard relation for lookback time: where  () = √︁ Ω  (1 + ) 3 + Ω  , with Ω  = 0.3, Ω  = 0.7 and  0 = 70kms −1 Mpc −1 .
We note that our merger rate density should be only considered as an order of magnitude estimate at best.This imprecision is due to several uncertainties in stellar physics and, notably, the limited density of our metallicity grid.We performed simulations only at two metallicities to determine the merger rate density.However, the formation efficiency and delay times of GW sources originating from CHE systems is expected to be sensitively dependent on metallicity.
In particular, we overestimate the delay times for GW sources formed at 0.001 <  ≤ 0.005, which in turn leads to an overestimation of the merger rate density at z = 0.This is because, we represent all systems formed in this metallicity range with our models at Z = 0.005, at which the stellar winds are stronger and therefore lead to wider BH-BH binaries.The longer time delays imply that GW sources merging at z = 0 are predicted to have formed at a larger redshift, at which the star formation rate is higher.In particular, Madau & Dickinson 2014 predicts that the cosmic star formation rate montonically increases up to  ∼ 2. This could also explain why our merger rate is a factor of two higher than predicted by (Riley et al. 2021).Similarly, we underestimate the delay times for GW sources formed at 0.0005 <  ≤ 0.001, and therefore we might underestimate the merger rate densities for such systems.In particular, this could mean that the merger rate density of the TMT with a MS-MS accretor channel (discussed in section 5.5) could be significantly higher than predicted (shown in Table 3).This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure1.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  1 ,  2 and inner separation  in as indicated in the top left of each panel) but with varying tertiary masses,  out (x-axis) and outer semimajor axes,  out (y-axis, logscale).The parameters of the inner binary in each panel are the same as that of an isolated CHE binary with  1,ZAMS =  2,ZAMS = 70  ⊙ ,  ZAMS = 22.4  ⊙ 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 i = 90 • .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 2, 4, and 3).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 .
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:  1,ZAMS =  2,ZAMS = 70  ⊙ ,  in,ZAMS = 22.4  ⊙ , i ZAMS = 90 • ,  out,ZAMS = 32.1  ⊙ ,  out,ZAMS = 200  ⊙ ,  out = 0.The outer eccentricity remains  out ≲ 0.01 throughout the evolution.
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.
Figure 4. We show the properties of the donor star at the onset of tertiary mass transfer episode.The right panels show the outer mass ratio,  out =  out /(  1 + 2 ), while the left panels show the amount of mass transferred from the tertiary to the inner binary for systems undergoing TMT at Z = 0.005, 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 (see e.g.Table2) .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.

Figure 5 .
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.

Figure 6 .
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  = 0.005.The parameters of the inner binary are the same for each system shown in the grid, i.e.  1,ZAMS =  2,ZAMS = 70,  ⊙ ,  in,ZAMS = 22.4  ⊙ ,  in = 0.The initial tertiary mass  out ranges from 5 to 100  ⊙ on a linear scale, while  out,ZAMS ranges from 200 to 1000  ⊙ on a logarithmic scale.We assume zero outer eccentricity and an initial inclination of 90 • .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 .
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 .
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.

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

Figure 12 .Figure A1 .
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)

Figure A2 .Figure A3 .
Figure A2.The same figure as Fig. 5 but for our model at Z = 0.0005.

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.

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.See equation B8 and the accompanying discussion in appendix B for the definition of birth rate density. ).