The Main Evolutionary Pathways of Massive Hierarchical Triple Stars

So far, stellar population studies have mainly focused on the evolution of single and binary stars. Recent observations show that triple and higher order multiple star systems are common, especially among massive stars. Introducing three-body dynamical effects can influence the evolution of an individual stellar system and hence affect the predicted rates of astrophysical sources that are a product of stellar evolution. We aim to constrain the main evolutionary pathways of massive hierarchical triple star systems and to quantify the effect of the third star on the evolution of the system. We model the massive triple star population by performing simulations of triple star evolution with the TRES code, which combines stellar evolution with secular evolution of triple systems, and explore how robust the predictions of these simulations are under variations of uncertain initial conditions. We find that interactions are common among massive triple stars. The majority of systems (65-77\%) experience a phase of mass transfer in the inner binary, often with an unevolved donor star. This differs significantly from isolated binary evolution, where mass transfer is less frequent (52.3\% instead of 67\% for our fiducial model) and the donors are typically post-main sequence stars. Initial constraints for dynamical stability as well as eccentricity oscillations driven by the third body facilitate the occurrence of interactions, for instance mass transfer. The requirement of dynamical stability at formation places quite stringent constraints on allowed orbital properties, reducing uncertainties in triple evolution that resort from these initial conditions. Ignoring three-body dynamics during evolution of non-interacting triples leads to triple compact-object systems with stronger eccentricity oscillations and thereby likely over-predicts the merger rate of compact objects in such systems.


Introduction
Massive stars (M 8M ) play an important role in the dynamical and chemical evolution of galaxies.Stellar feedback in the form of stellar winds, radiation pressure, photo-ionization and supernova (SN) explosions control star formation (e.g.Hayward & Hopkins 2017) and enrich the interstellar medium with metals (e.g.Burbidge et al. 1957;Larson 1974;Larson & Dinerstein 1975).Most massive stars are observed to have at least one stellar companion (e.g.Evans et al. 2006;Sana et al. 2014;Kobulnicky et al. 2014;Moe & Di Stefano 2017) and are expected to interact at some stage of their evolution (Sana et al. 2012).These interacting systems can be the progenitors for a great variety of observed astrophysical sources, such as gravitational wave (GW) mergers (e.g.Abbott et al. 2016), X-ray binaries (e.g.Verbunt 1993;Langer 2012), massive equivalents of Algol systems (e.g.Budding 1989) and blue stragglers (e.g.Perets & Fabrycky 2009), and ultra-stripped supernovae (e.g.Tauris et al. 2013).Mass transfer, specifically, can significantly impact the subsequent evolution of both the donor and accretor stars (e.g.Eldridge et al. 2008;Langer 2012;Laplace et al. 2021).
In the past, the vast majority of theoretical work in stellar evolution has focused on single-and binary stellar models.However, more recent advances in observational astronomy reveal that higher order multiples are prevalent, especially among the massive star population.Roughly 10% of the solar-type star systems are triples (Tokovinin 2014;Moe & Di Stefano 2017), while for early B-type and O-type stars this fraction increases to over 30% (Evans et al. 2006;Moe & Di Stefano 2017).Unsurprisingly, stellar products have been observed whose formation poses a challenge from the perspective of single-and binary evolution, but can more easily be explained with the inclusion of a third companion, e.g., low mass X-ray binaries containing a black hole (Podsiadlowski et al. 2003;Naoz et al. 2016), barium stars (Gao et al. 2023), cataclysmic variables with high accretion rates (Knigge et al. 2022) and binary stars with one magnetic component (Alecian et al. 2014;Schneider et al. 2016Schneider et al. , 2019;;Shultz et al. 2019).Therefore, complimentary population synthesis studies of triples are essential toward creating a more comprehensive picture of stellar evolution and interaction.
Introducing a third companion to a binary adds complexity to the system.On top of stellar evolution and interactions, threebody dynamics has to be taken into account to accurately describe the evolution of a triple star system over time.The lowest order manifestation of these three-body interactions in a stable triple system are the Von Zeipel-Lidov-Kozai (ZLK) oscillations (von Zeipel 1910;Lidov 1962;Kozai 1962), which are periodic variations of the eccentricity of the orbit of the inner two stars and the inclination of the third star with respect to the orbital plane of the inner orbit.For an elaborate review on the effects of three-body dynamics, see Naoz (2016).Recent population synthesis studies of triple stars have shown that increased eccentricities induced by three-body dynamical effects such as ZLK oscillations can lead to an enhanced degree of stellar interaction between the components of the inner binary compared to such interactions in isolated binary stars (Antonini et al. 2017;Hamers & Thompson 2019;Toonen et al. 2020;Stegmann et al. 2022a,b;Hamers et al. 2022a,b).
In this paper, we aim to provide an inventory of the most common evolutionary channels within a population of massive triple stars and show how these predictions differ from standard binary evolution.To explore this, we perform simulations of massive hierarchical triple stars with a triple population synthesis code.We stop the computation of an individual stellar and orbital evolution when the system becomes dynamically unstable, unbound or at the onset of mass transfer, i.e., for this subset of evolutionary channels we do not continue until the final remnant stages.The main reasons for this are that the physics of mass transfer in a triple configuration is still poorly understood and awaits more detailed studies, and that it allows us to focus on the effects of a third companion on the initial orbital conditions and early evolution.
In Section 2, we give a description of the triple population synthesis code and motivate our choice of initial model conditions.In Section 3, we present the predicted incidence rates of each evolutionary channel and discuss a few channels in more detail.In Section 4, we compare our results with similar studies and discuss some caveats.Finally, we give a summary of our findings in Section 5.

Method
To investigate the most common evolutionary channels in massive hierarchical triple stars we use the triple population synthesis code TRES1 (Toonen et al. 2016).This is a code that couples a single and binary stellar evolution code, in this case SeBa (Portegies Zwart & Verbunt 1996;Toonen et al. 2012), which models stellar evolution, and a code that solves the orbital motions of three-body systems using the secular approximation (Toonen et al. 2016).In the secular approximation, the system is sufficiently hierarchical such that the tertiary star effectively orbits the center of mass of the inner two stars, forming the outer orbit.In these configurations, the energy of the outer orbit and the orbit of the inner two stars is separately conserved (Heggie 1975).To model stellar evolution, SeBa makes use of analytic fitting formulae from Hurley et al. (2000Hurley et al. ( , 2002)), based on the single stellar evolution models produced by Pols et al. (1995).These fitting formulae allow for rapid simulation of stellar evolution, essential for population synthesis studies.
With TRES, we simulate ten sets of 25000 initially dynamically stable, coeval, massive hierarchical triple stars that are initialized at the zero age main sequence (ZAMS).These systems are then evolved for approximately a Hubble time (13.5 Gyr), unless the systems experience a phase of mass transfer, become dynamically unstable or unbind, after which the simulation is terminated.In a very small number of simulations ( 0.05%) the timescale for dynamical interactions is extremely short compared to the nuclear timescale; we decided to discontinue these models for practical reasons.For a similar number of systems, the secular code could not find an analytical solution for the orbital elements.These too were discarded.
We will first discuss the adopted treatment of mass loss in stellar winds and our choices for the initial stellar and orbital properties of the primordial triple star population.Subsequently, we review a few physical implications of multiplicity that can affect the evolution of a system.

Updated wind prescription
By default, SeBa follows the line-driven wind mass-loss model for massive stars of Vink et al. (2000Vink et al. ( , 2001) ) in the metallicity range 0.03 Z < Z < 3 Z and Nieuwenhuijzen & De Jager (1990) in the range where the Vink models are not applicable.We have modified these prescriptions by reducing the mass-loss rates by a factor 3, following and generalising the findings of Björklund et al. (2021) for the implemented mass range.Additionally, we have assumed luminous blue variable mass-loss rates following Belczynski et al. (2010), introducing a constant enhanced mass loss of 1.5 × 10 −4 M /yr for stars exceeding a luminosity of 6 × 10 5 L .

Initial parameters
To simplify the description of a hierarchical triple star, the system is split up into two subsystems: (i) The inner binary, that comprises the primary and secondary star with masses m 1 and m 2 .These stars orbit each other with an inner semi-major axis a in , eccentricity e in and argument of pericenter g in ; (ii) The outer binary, that comprises the centre of mass of the inner binary and the tertiary star, with mass m 3 .The outer binary has an outer semi-major axis a out , eccentricity e out and argument of pericenter g out .We define a relative inclination angle i rel between the orbital planes of the inner and outer orbit.
The incidence rates of transients depend on the ZAMS conditions of a stellar population (e.g.Abate et al. 2015;de Mink & Belczynski 2015;Stevenson et al. 2022).Therefore, we would ideally like to possess a complete understanding of how the initial stellar and orbital parameters of massive triples are distributed.Unfortunately, the observed sample size of massive stars is limited compared to low-mass stars and our understanding is fragmented, as certain parts of the parameter space (e.g., low mass ratios and wide orbits) are difficult to probe.Moreover, observations suggest there is a degeneracy in constraining certain parameter combinations (Moe & Di Stefano 2017).To account for these uncertainties, we produce a set of model variations with a diverse set of assumptions on the primordial parameter distributions.
Our fiducial model follows a Kroupa IMF distribution (Kroupa 2001) for the stellar mass m 1 in the range [10,100] M .This distribution is a power-law function with index α = −2.3.Above a mass of 100 M , extrapolation of the Hurley stellar evolution fitting formulae becomes less accurate.The stellar masses m 2 and m 3 are specified through the inner mass ratio q in ≡ m 2 /m 1 and the outer mass ratio q out ≡ m 3 /(m 1 + m 2 ), respectively.Both the inner and outer mass ratios are sampled from a flat distribution in the range [0.1,1]2 , based on mass-ratio observations of massive binary stars (Sana et al. 2012;Kobulnicky et al. 2014).The orbital semi-major axis range for a in and a out both cover [5,5×10 6 ] R , and the orbital distribution functions (see below) are sampled uniformly in logarithmic space, following (Öpik 1924;Kobulnicky & Fryer 2007).Below separations of 5 R , the stars of the inner binary are likely to touch at the ZAMS.Above 5 × 10 6 R , the system is only very weakly bound through gravity and can easily be disrupted by small perturbations, such as stellar flybys (Kouwenhoven et al. 2007(Kouwenhoven et al. , 2010)).When the sampling from our orbital semi-major axis distribution function yielded a in > a out , we decide to swap the values, as, by definition, a in < a out .This is a different approach from, e.g., one of the models in Toonen et al. (2020), who fixed the primordial distribution of a in and accepts or rejects a out in accordance with the stability properties of the system.The resulting ZAMS distributions of a in and a out might differ between both methods.The eccentricities e in and e out are sampled from a thermal distribution in the range [0,0.9](Ambartsumian 1937 ;Heggie 1975;Moe & Di Stefano 2017;Hwang 2023).The arguments of pericenter g in and g out are sampled from a uniform distribution in the range [-π,π].Lastly, the relative inclination i rel is sampled from a uniform distribution in cosine in the range [0,π].
In each population model we vary only one (two in a single case) of the aforementioned distributions, in an effort to isolate the impact each parameter has on the evolution of the population.Five model variations assume the distributions for the semi-major axes and eccentricities as described by Sana et al. (2012), based on an observed population of O-type multiple stars in nearby galactic open clusters.The orbital periods are sampled from a power law in the range 0.15 < log 10 P < 8.5, with index π = −0.55 and are then converted to semi-major axes using Kepler's third law.The eccentricities are sampled from a power law with index η = −0.45.In two model variations we assume a flat distribution of e in and e out (Kobulnicky et al. 2014).In one variation, we assume a power law with index γ = −1.5 for the outer mass ratio, q out , which seems more appropriate for early B-and O-type primaries with wide companions (Moe & Di Stefano 2017).For the final variation, we assume i rel = 0.This model is motivated by Tokovinin (2017), who showed that for outer orbital separations smaller than about 10 4 R the inner and outer orbits appear to be more aligned.Table 1 provides a complete overview of the parameter distributions, sampling ranges and model variations.

Dynamical stability
To ensure long-term stability of a three-body system, a hierarchical configuration is generally required, i.e., the outer binary has a considerably larger orbit than the inner binary.In that case, the timescale at which the tertiary perturbs the system is long compared to the dynamical timescale of the inner binary.When the tertiary component orbits the inner binary too closely, the two timescales become comparable, breaking down the hierarchical structure and rendering the secular approximation invalid.We define the critical ratio for hierarchical stability following Mardling & Aarseth (1999, 2001) as a criterion for hierarchical stability: . (1) A system is dynamically stable if a out /a in > (a out /a in ) crit .For instance, for a circular orbit, equal mass system (q out =0.5), with the inner and outer orbit in the same plane, (a out /a in ) crit = 3.3.
Applying this stability criterion to our primordial population of triple stars affects the shape of the initial stellar and orbital parameter distributions.Fig. 1 shows how the sampled distributions of a in and a out are altered by rejecting systems that are dynamically unstable or Roche-lobe filling at birth.While the fiducial model follows a flat sampling distribution in the logarithm of both a in and a out , the population of hierarchical triple stars that satisfies the stability criterion (black line) is clearly more biased against large (small) inner (outer) semi-major axes.The most compact inner binaries (a in < 20 R ) are often in contact at initialisation and are discarded.However, the number of systems discarded due to contact is much smaller for models with lower initial inner eccentricities.The important thing to notice is that the semi-major axis distributions of the ZAMS population of models that vary a in and/or a out are quite similar to the fiducial model, even though their sampling distributions are vastly differing.

Three-body dynamics
The perturbations induced by the tertiary star onto the inner binary can be described through a mutual torque acting between the inner-and outer orbit.This torque allows angular momentum to be transported from one orbit to the other and vice versa, leading to an oscillatory behavior of the inner eccentricity and the relative inclination.The lowest order manifestation of this mechanism, the quadruple regime or ZLK effect, acts on a timescale that is mostly dependent on the ratio of orbital periods, P 2 out /P in , and e out .An approximate expression for this timescale is (Kinoshita & Nakai 1999;Antognini 2015): In the test-particle approximation, large oscillations of the inner eccentricity occur only at initial relative inclinations i init between 39.2 • − 140.8 • for initially non-eccentric inner orbits (Ford et al. 2000;Naoz et al. 2013;Naoz 2016).The maximum inner eccentricity reached through ZLK oscillations in the test-particle regime can be described analytically for circular orbits (Innanen et al. 1997;Grishin et al. 2017): Close to dynamical destabilisation and at non-zero outer eccentricities, higher order terms of the secular approximation treatment become important, such as the octupole term (Ford et al. 2000).Besides reaching extreme eccentricity amplitudes, in this regime three-body dynamics can provoke a flip in relative inclination (from prograde to retrograde or vice versa) (Naoz et al. 2011;Li et al. 2014) and inner eccentricity variations can occur at a wider range of relative inclinations than 39.2 • −140.8 • (Ford et al. 2000).The octupole term is relevant when the octupole parameter, is of the order | oct | 0.001-0.01(Lithwick & Naoz 2011;Shappee & Thompson 2013).We include a new function in the code that tracks the amplitude of the inner eccentricity oscillations during the entire evolution of a system.Per time step of the evolutionary code, the largest eccentricity amplitude is stored.This information can be used as a tool to quantify the impact of a third stellar component on the dynamical evolution of the system.In Fig. 2 we show the maximum change in inner eccentricity during each time step of the evolutionary code for an example system.The system has initial stellar masses m 1 = 20 M , m 2 = 15 M and m 3 = 30 M , initial semi-major axes a in = 2 × 10 4 R and a out = 1.5 × 10 5 R , initial eccentricities e in = 0.2 and e out = 0.4, and an initial relative inclination i rel = π/2.Until about 4.5 Myr, three-body dynamics persistently induce variations of the inner eccentricity up to 0.8.As the octupole parameter of this system is 0.009, the octupole term is important, which is manifested as the sinusoidal variations on timescales of about 0.8 Myr.With the analytical expression for the octupole timescale from Antognini (2015), we derive a comparable approximate oscillation time of 1.25 Myr.

Tidal interaction
Stars in a binary system are subjected to tidal forces that act to synchronize the stellar rotation with the orbit and circularize the orbit by the dissipation of energy.Orbital properties, such as the semi-major axis, the eccentricity and the angular velocity are thereby affected and alter the further evolution and possibly final fate of the system.The tidal forces are strongly dependent on the ratio between the stellar radii and the semi-major axis, and hence only significantly impact compact binary systems.We adopt the equilibrium tidal model described by Hut (1980) for stars with convective envelopes.For stars with a radiative envelope, we adopt the dynamical tidal model described by Zahn (1977).Both the equilibrium and dynamical tidal models are parameterized by Hurley et al. (2002) and implemented in TRES.
In orbital configurations where both the tidal dissipation and secular timescales are short compared to the evolutionary timescale of the stars, interplay between tides and three-body dynamics can efficiently reduce the semi-major axis of the inner binary (Mazeh & Shaham 1979;Kiseleva et al. 1998); ZLK oscillations increase the inner eccentricity, allowing tides to efficiently dissipate energy during pericenter passage (e.g.Fabrycky & Tremaine 2007).After complete circularization, the semimajor axis has decreased by a factor (1 − e 2 in ).

Results
We briefly describe the conditions at which we stop individual stellar evolution simulations.We differentiate between the following channels: -Mass transfer in the inner binary when either the primary or secondary star fills its Roche lobe.-Mass transfer from the tertiary star onto the inner binary, where the material may be accreted by one or both components.-Dynamical destabilisation due to stellar evolution (e.g., orbital widening due to mass loss) violating stability criterium Eq. 1. -Dynamical unbinding of the inner orbit when either the primary or secondary star experiences its core-collapse supernova.-Dynamical unbinding of the outer orbit as a result of corecollapse of one of the three system components.-None of the above interactions occur within a Hubble time.
At this point all three components are compact remnants.We refer to this channel as non-interacting systems.We do not explicitly identify the merging of two stars as a stopping condition as such events are preceded by mass transfer.Note that we do not include tidal effects and three-body dynamical effects such as ZLK oscillations as stopping conditions, while technically these could be described as orbital interactions.An overview of the evolutionary channels leading up to a model discontinuation is presented in Table 2.The remainder of this section is structured as follows: in Section 3.1 we give an overview of the predictions of the evolutionary channels for our simulated populations.Then we discuss two  Non-Interacting channels in more detail; the systems that experience a phase of mass transfer initiated by the primary star in Section 3.2 and the systems that do not undergo any interaction during their entire evolution in Section 3.3.Next, we dive into few of the model variations in more detail in Section 3.4.Finally, we compare the predictions for massive triple star evolution with simulations of isolated binary stars in Section 3.5.

Overview evolutionary channels
In Fig. 3 we show the predicted contribution of each evolutionary channel to the total population of massive triple stars for our fiducial model.Similar to massive binary star systems (e.g.Sana et al. 2012), the evolution of massive triple star systems is dominated by both stellar and orbital interactions.The vast majority of systems (67.3±0.9%)experiences an episode of mass transfer in the inner binary.This includes mass transfer initiated by both the primary and the secondary star, although the contribution of the latter is only very minor ( 0.1%).The error bars are the 3σ statistical sampling uncertainties obtained by bootstrapping the data (i.e., re-sampling the data with replacement).These uncer- tainties follow a Poisson distribution by approximation.A fairly large fraction of systems loses a companion as a result of a SN kick, unbinding either the inner-or outer orbit (10.0 ± 0.6% and 18.9 ± 0.7% respectively)3 .In case the velocity kick is not strong enough to unbind the orbit, it can still affect the system by altering the semi-major axis and the eccentricity (Pijloo et al. 2012;Lu & Naoz 2019).Only a few percent of the systems experiences mass transfer initiated by the tertiary star (1.8 ± 0.2%) or reach dynamical destabilisation (1.7 ± 0.2%).By far the smallest fraction of systems (0.2 ± 0.1%) has not engaged in one of the aforementioned interactions within a Hubble time.A detailed overview of the rates of all population models can be found in Appendix A.

Model variations
Apart from statistical uncertainties, the predicted contributions of evolutionary channels are dependent on the uncertainties related to the formation of the stellar system, i.e. initial stellar and orbital parameters.In Fig. 4, we show the predicted ranges on the contributing fraction of each evolutionary channel resulting from the model variations discussed in Section 2.2 relative to the fiducial model.Overall, the incidence rate of each evolutionary channel is fairly robust across the model variations; The rates only twice differ over a factor 2 compared to the fiducial model.The most significant outliers mainly originate from one model variation, Model q out -Moe.We discuss some of the outliers in more detail in Section 3.4.Moreover, the predicted rates for the fiducial model are close to the median for most of the evolutionary channels, confirming that the fiducial model is a good representation of an average model.For the more uncommon channels, such as tertiary mass transfer, non-interacting and dynamically destabilised systems, the fiducial model predictions agree less well with the median.Although, we can not rule out that this deviation is explained by the enhanced statistical uncertainties for the latter two channels.
In conclusion, the general evolution of a massive triple star population leading up to the first interaction is not seriously af- fected by the differing assumptions of the initial conditions.This can partially be explained by the initial stability check described in Section 2.3, reducing the differences in the ZAMS population between models with vastly differing initial parameter distributions.This is fortunate, as the initial parameters are empirically not well constrained.We stress that robust predictions for the incidence rates of evolutionary channels between two model variations do not necessarily imply that the complete evolution of stellar systems in both populations are identical.We terminate the simulation at the onset of the first interaction, but the outcome of certain interactions, such as mass transfer, and thus the subsequent evolution of the system, can be highly dependent on the stellar and orbital properties of the system at the ZAMS and similarly at the onset of the interaction.This is specifically important for the formation of more uncommon sources, such as gravitational wave mergers.

Effect of initial conditions
We provide an overview of a subset of the initial conditions that generally lead to a specific evolutionary channel in Fig. 5, based on the triple systems in our fiducial model.
-Mass transfer: Mass transfer initiated by the primary star can occur over a wide range of the inner and outer orbital period (1 − 10 6.5 days for P in and 20 − 10 8.5 days for P out ), as the maximum radius of massive stars can be large.Even at P in > 10 4 days, extreme inner eccentricities induced by ZLK oscillations can provoke the primary star to fill its Roche lobe during pericenter passage.For a few systems, the mass trans-Fig.5. Initial inner-and outer orbital periods covered by each evolutionary channel.The triangle shape of the population is a direct consequence of the requirement on hierarchical structure and dynamical stability.The black dashed line shows where P in = P out .There are no systems near this line due to the dynamical stability criterion.There is some overlap between colors, as the evolution of a system is determined by more properties than just the orbital periods.
fer in the inner binary is instead initiated by the secondary star.This can be counter-intuitive as the nuclear timescale of the primary star is shorter than that of the secondary star, hence the primary will evolve faster.A possible way to prevent the primary star from transferring mass first is by stripping the envelope via stellar winds before the star expands sufficiently to fill its Roche lobe.This does require the primary star to be very massive and is susceptible to the choice of wind model.Systems in which the secondary star transfers mass as their first interaction are restricted to initial inner orbital periods larger than 10 4 days.As the triple systems are hierarchical, and thus the outer orbit is significantly wider than the inner orbit, the first phase of mass transfer will only be initiated by the tertiary star if it is initially more massive than either component of the inner binary.Therefore, the distribution of q out is a key parameter in determining the contribution of this channel and P in and P out are restricted within about 10 3 and 10 5 − 10 6 days, respectively, for this channel.-Orbital unbinding: The initial conditions for systems whose inner (outer) orbit become unbound as a result of a SN kick are restricted to inner (outer) orbital periods larger than a few thousand days.Their subsequent evolution is difficult to predict without proper n-body calculations, as a phase of nonsecular and possibly chaotic evolution ensues.However, it is not unlikely to expect at least one of the bodies to be ejected, either directly as a result of the SN kick, or through subsequent dynamical interactions (Pijloo et al. 2012).-Dynamical destabilisation: most systems that become dynamically unstable do not have very large initial ratios of P out /P in (generally within a hundred, apart from a few exceptions).The location of instability is dependent on other parameters as well, such as e out , q out and i rel , resulting in a quite generous spread of initial period ratios where systems experience a dynamical destabilisation.We identify several processes responsible for the transition from dynamically stable to unstable, the first being systems born close to the instability limit.Early on in their evolution, when the primary star is still on the main sequence (MS), three-body dynamical interactions can increase e out or decrease i rel , shifting the critical semi-major axis ratio towards larger values (Toonen et al. 2020).Furthermore, once the primary star has evolved off the MS, either the increased mass-loss rate, resulting in the widening of the inner orbit, or an eventual SN kick, can push the system across the instability limit.-Non-interacting: Systems that do not have any form of interaction after a Hubble time of evolution are restricted to initial inner orbital periods of 6×10 3 days.Smaller periods would most likely result in a phase of mass transfer within the inner binary.

Time-evolution of evolutionary channels
The stellar and orbital properties of a system are affected when interactions between stars or orbits take place.Since interactions occur at different times for each system, the observational properties of a stellar population are expected to change over time.In Fig. 6, we present an overview of the typical stellar ages where each channel dominates.The non-interacting systems have been excluded, simply because their definition comprises the absence of any interaction.At very early times ( 4 Myr), we predict solely dynamical destabilisation of the triple and mass transfer in the inner binary to take place.The extend of this epoch corresponds to roughly the MS lifetime of the most massive stars in our simulations.Both channels show a steep rise initially (< 1 Myr), comprising systems that have short secular timescales and/or are initialised near the instability limit.Later on, the rate of interactions are small.
At later times ( 4 Myr), due to the fast radial expansion of post-MS stars and an eventual SN kick, the number of interactions increases rapidly within each channel and peaks just after 6.9 Myr.Interestingly, the systems in which the outer orbit becomes unbound typically interact at earlier times than systems in which the inner orbit becomes unbound.This may suggest that the former channel is biased towards high tertiary masses, as those stars evolve into a compact object on a shorter timescale.This is generally not the case.The explanation is that the inner binaries are generally more compact and have a higher binding energy compared to the systems in which the inner orbit becomes unbound.Therefore, for the inner binary the low kick velocities of the highest mass stars in our simulation are simply not sufficient to unbind their orbits.The rate of systems that undergo orbital unbinding of the outer binary or mass transfer from the tertiary onto the inner binary drops quickly after about 15 Myr.The reason for this is that m 1 has a sharp cut-off at 10 M , while m 3 extends to masses as low as 1.2 M .These low-mass stars do not experience a SN at the end of their lives, and reach different maximum radii compared to their massive counterparts.The exclusion of low-mass primary stars leads to a sharp drop in the number of systems that undergo a phase of mass transfer in the inner binary after 25 Myr.
Overall, we predict that within 3.7 Myr about 10% of systems have had an interaction.This fraction rapidly increase to ∼50% after 10.2 Myr.After that, the rate of interactions slows down, despite the majority of stars evolving on timescales longer than 10 Myr as a consequence of the initial mass function.This trend follows the declining rate of inner binary mass-transfer systems, as that constitutes the most frequent type of interaction.
Lastly, we ask ourselves how the triple population looks like for a realistic stellar population.This can be done by convolving the results from Fig. 6 with a star formation history.In Fig. 7, we show how the intrinsic population of massive triples looks like after 30 Myr, which is the typical lifetime of star-forming giant molecular clouds, assuming a constant star-formation rate, rather than a single burst of star formation.The majority of systems (62.3%) has already interacted.Most systems (45.1%) have undergone a phase of mass transfer, initiated by any of the three stars.The high incidence of mass transfer is consistent with predictions for massive binary stars (de Mink et al. 2014) and emphasizes the importance of mass transfer in multiple stars.A significantly smaller fraction of the systems (15.9%) has ejected at least one of the stars, while few systems (1.3%) have destabilised as a result of three-body dynamics.de Mink et al. (2014) showed that almost a third of post-mass-transfer binary systems result in a merger.Our findings could thus indicate that a large fraction of systems born as massive triple stars would be reduced to a population of single and binary stars at the moment of observation.

Mass transfer initiated by primary
The majority of systems (64.6-76.6%)has a phase of mass transfer in the inner binary.As we stop our simulations at the first moment of interaction, the primary star will be the donor and the secondary star the accretor in most cases.One of the key ingredients in determining the outcome of mass-transfer systems is the stability of mass transfer.When a massive donor star approaches the giant branch, its envelope will transition from radiative to convective, and the mass transfer becomes more prone to instabilities (Klencki et al. 2021).Therefore, the evolutionary phase of the donor is a useful property to explore.We show the evolutionary phases of the donor and accretor stars at the onset of mass transfer for all model variations in Fig. 8.Given in chronological order, the phases at which the donor transfers mass are the Main Sequence (MS), Hertzsprung Gap (HG), First Giant Branch (FGB), Core Helium Burning (CHeB) and Asymptotic Giant Branch (AGB) phase.The CHeB phase is defined as the moment when helium is ignited in the core, which can occur before the star reaches the FGB phase for sufficiently high masses.These systems move directly from the HG phase to the CHeB phase, avoiding the FGB phase.We only discuss systems with a MS accretor, as the fraction of systems in which the secondary is an evolved star is extremely small ( 0.5%).Since the evolutionary phases after the MS are short lived compared to the total stellar lifetime, this can only occur when the donor and the accretor have similar evolutionary timescales, and thus similar initial masses.Nonetheless, these systems can lead to interesting events, such as a double core common envelope (Dewi et al. 2006).
In the vast majority of systems, the donor star is a MS star (32.8-51.7%)or a HG star (37.8-50.5%).Only 7.9-12.7% of the donor stars are in the CHeB phase, while the FGB and AGB donors contribute merely a few percent.As discussed before, not all stars evolve through the FGB phase.Over 75% of the primary stars are initially more massive than about 13M and move directly to the CHeB phase.Similarly, stars initially more massive than about 38M do not evolve through the AGB phase for our choice of wind mass-loss prescription and metallicity.Instead, they transition to the Wolf-Rayet phase while burning helium in the core.Additionally, for stars that do evolve through the AGB phase, the fractional increase in radius is generally marginal compared to the HG and CHeB phases.
We discuss two effects that explain the large number of systems that undergo mass transfer during an early evolutionary phase of the donor.First, as shown in Fig. 1, the ZAMS distribution peaks at small inner semi-major axes as a result of satisfying the criterion for dynamical stability, disfavouring large inner semi-major axes.At small separations, the donor star does not need to expand much to fill its Roche lobe, hence reaching the onset of mass transfer early in its evolution.Second, three-body dynamical interactions can expedite the onset of mass transfer by decreasing the pericenter of the inner binary through oscillations of the inner eccentricity.
To quantify the impact of three-body dynamics on the systems that go through a phase of mass transfer initiated by the primary star, we show the maximum amplitude of the eccentricity oscillations for each system of the fiducial model in Fig. 9.In 9.1% of systems the inner eccentricity is increased between 0.05 and 0.96.The population model affected the most and the least by three-body dynamics are model e in -Sana, with low initial inner eccentricities, and model i rel -const, with initial relative inclinations of zero, respectively.For these models, 13.7% and 1.1% of the systems that go through a phase of mass transfer initiated by the primary star have experienced eccentricity oscillations with amplitudes between 0.05 and unity.The affected systems generally have small initial ratios of P out /P in , ensuring relatively short timescales on which the three-body dynamical interaction acts.However, there is a generous spread towards orbital period ratios up to over a hundred, since additional parameters (e out , q out and i rel ) affect the timescale and magnitude of the oscillations.

Non-interacting systems
Across all model variations only 0.1-0.4% of the systems form a bound triple compact object (TCO) system that has had no previous interaction during their entire evolution.We have shown that these systems have initially relatively wide inner orbits to prevent a phase of mass transfer, and are therefore also required to have wide outer orbits to ensure a stable hierarchical configuration.At these wide separations however, the weakly bound tertiary star is prone to become unbound due to a SN kick, justifying the small contribution of this channel.
Nevertheless, systems devoid of any interaction are an interesting target to study, as they are considered to be a main formation channel for producing gravitational wave mergers in hierarchical triple stars (Antonini et al. 2017;Silsbee & Tremaine    In view of these possible compact object merger futures, we discuss the properties of the systems from the non-interacting channel.In Fig. 10, we show the distribution of initial parameters that are of importance for driving the eccentricity oscillations.We compare these initial parameter distributions to the properties at TCO formation of the non-interacting systems for the fiducial model.Most of the properties are distributed quite dissimilar between the two populations.The TCO population evidently favors more moderate inner-and outer eccentricities.Highly eccentric systems are more likely to interact and hence prevent forming a TCO that has had no interaction during its evolution.Interestingly, the models e in -Sana and e out -Sana, that peak at much smaller initial inner-and outer eccentricities respectively, have eccentricities at TCO formation that are very similar to those in the fiducial model.For all model variations, near-equal mass ratios between the primary and tertiary compact object (m 3 /m 1 ) are strongly favored.The system has to survive three subsequent SN kicks in order to form a bound TCO.For our choice of SN kick-velocity model, each star of the triple system needs to be massive to prevent break-up.Consequently, we predict a dearth of initial primary and tertiary masses below 35 M .At high masses, all stars tend to evolve towards similar final black-hole masses due to high rates of wind mass loss.The orbital period ratios P out /P in are almost completely limited to ratios smaller than 10 4 for the TCO systems.The minimum inner orbital period to avoid mass transfer in the inner binary is around 6 × 10 3 days (see Fig. 5).As the outer orbital periods are confined to 10 8.5 days, the orbital period ratio never reaches extreme values.The TCO orbits slightly disfavor relative inclinations of 90 • .At these inclinations, the eccentricity oscillations in the inner binary are the largest, complicating the formation of a noninteracting TCO.The discussed properties set the timescale for the ZLK oscillations, which partially determines how effectively GW mergers can be produced due to the presence of a third compact object.The ZLK timescales (Eq.2) for the TCO systems are significantly longer than for the complete initial population and range from 100 kyr to 1 Gyr.However, ZLK oscillations and possibly higher order three-body dynamical interactions can still affect the system over the span of a Hubble time (∼ 13.5 Gyr).

Differences in the incidence of evolutionary channels and onset of mass transfer between model variations
The rates at which specific interactions occur vary quite significantly from the fiducial model for a few of the model variations (Fig. 4).In particular, models that favor small inner semi-major axes, such as model a in -Sana and a in &a out -Sana, or low tertiary masses, such as model q out -Moe, are responsible for the most significant differences.Model a in -Sana and a in &a out -Sana show an increase in the incidence rate of systems that undergo mass transfer in the inner binary, with a rate of 73.4% and 75.1%, respectively, compared to the fiducial model's rate of 67.3 ± 0.9%.These models have more compact inner orbits, resulting in smaller Roche radii of the inner binary components and thus increasing the likelihood of mass transfer.Model a in &a out -Sana also has more systems where the tertiary star transfers mass onto the inner binary, with a rate of 3.2%, compared to the fiducial model's rate of 1.8±0.2%.In this case, the outer orbits are more compact, making it more likely for the tertiary star to fill its Roche lobe.
Conversely, both models predict a decrease in the incidence rate of systems which unbound their inner orbit, with a rate of 5.9% and 5.1%, respectively.In the fiducial model, the inner orbits were unbound in 10.0 ± 0.6% of all systems.We identify two processes that contribute to this decrease.Firstly, the inner orbits are typically more compact, resulting in a larger binding energy, which means that the range of supernova kick velocities and directions required to unbind the orbit becomes smaller.Secondly, at smaller semi-major axes, the system is more likely to undergo mass transfer before the supernova occurs that ultimately unbinds the orbit.
Similar to models a in -Sana and a in &a out -Sana, our predictions for model q out -Moe indicate an increase in the incidence rate of systems undergoing mass transfer in the inner binary, with a rate of 76.6%.However, we also predict a significant decrease for systems transferring mass from the tertiary star onto the inner binary, with a rate of 0.45%, compared to the fiducial model's rate of 1.8 ± 0.2%.Additionally, we observe a large decrease in the incidence rate of systems experiencing an unbinding of the outer orbit.The predicted rate is 6.9%, compared to the fiducial model's rate of 18.9 ± 0.7%.These differences arise due to the fact that the population peaks at small initial outer mass ratios, resulting in lower tertiary masses.When the tertiary star is less massive than either component of the inner binary, its evolutionary timescale is longer, allowing the inner binary to initiate mass transfer or unbind the inner orbit via a supernova kick before the tertiary star has evolved sufficiently to interact.Moreover, a larger fraction of tertiary stars in model q out -Moe will not experience a SN kick at all, but will instead end their lives as white dwarfs.
For a range of models we predict a significant relative increase in the incidence rate of systems that do not interact, with values ranging from 0.38% to 0.42%, compared to the fiducial model's rate of 0.24%.This includes all model variations with more moderate inner or outer eccentricities (e in -flat, e out -flat, e in -Sana & e out -Sana) and a low relative inclination (i rel -const).At low eccentricities, the pericenter distance as well as the specific binding energy between stars is larger, making it less likely for the system to initiate mass transfer or unbind.At low relative inclination, three-body dynamical effects are less important and will not be able to induce eccentricity oscillations that drive the system toward interaction.
In the models a in -Sana and a in &a out -Sana, donors are usually less evolved, and the percentage of main sequence (MS) donors increases to 49.3% and 51.7%, respectively, compared to the fiducial model's 38.1%.However, there are sligincludely fewer donors in later evolutionary phases, which is expected because a smaller inner semi-major axis results in a smaller Roche lobe for the donor star.For model i rel -const, we expect fewer MS donor stars and more HG donors.This is due to the fact that three-body dynamics are not effective at low relative inclinations, resulting in negligible inner eccentricity oscillations.Consequently, the pericenter of the inner binary does not decrease, and the onset of mass transfer is expedited.

Comparison with isolated binaries
In this section we quantitatively investigate how the evolution of the inner binary is affected by a tertiary star.Previously, we have come across two mechanisms that could alter the evolution of a system with a tertiary component: (1) the initial ZAMS properties are affected by constraints of triple star formation and (2) three-body dynamical effects alter the orbital elements.In order to study the implications of both mechanisms on the evolution of the systems, we present two additional simulations of a population of isolated binary stars based on the same initial sampling distributions as the fiducial model.The first simulation represents a population of true isolated binary stars; i.e. the systems are initiated and evolved without the presence of a tertiary star.We call this model Bin.The systems in the second simulation are initiated as a hierarchical triple (identical to the fiducial model), forcing the initial distributions of the inner orbital parameters to become skewed to ensure a stable hierarchical triple configuration.Subsequently, the tertiary star is removed.Therefore, when we model the evolution of this system from the ZAMS onward, there are only two stars in the system.We call this model Bin-Skew.The inclusion of the latter model allows us to disentangle the contribution of tertiary induced dynamical effects from the selection effects introduced at star formation.
We first investigate how significantly the general evolution of an isolated binary population is affected by skewing the orbital parameter distributions due to the presence of a tertiary star.To this end, we compare the predicted incidence rates of evolutionary channels for the models Bin and Bin-Skew (Fig. 11).There is an enormous discrepancy in the predicted number of systems that undergo mass transfer and systems whose orbit becomes unbound between both models.The model Bin-Skew has about one and a half times as many systems that experience an episode of mass transfer, 80.7% compared to 52.3%, and over a factor of two decrease in the number of orbits that become unbound, 18.0% compared to 45.0%.Furthermore, the number of systems that do not have any interaction during their entire evolution decreases by over a factor of two, 1.2% compared to 2.7%.The differences between these models are primarily a result of a smaller average pericenter at the ZAMS for the systems from model Bin-Skew.Interestingly, these differences are much bigger than the differences between the triple star models, as were shown in Fig. 4. In short, dynamical constraints imposed by a third star during the initialisation of a system can alter the interaction history for a large number of systems.
Next, we investigate the effects of three-body dynamics on the evolution of the inner binary.To this end, we compare the occurrence of mass transfer from the primary onto the secondary star at different evolutionary phases of the donor star between the fiducial model and Bin-Skew (Fig. 12).For the fiducial model, we find an increase of 7.6% in systems that initiate mass transfer on the MS compared to model Bin-Skew.Conversely, during later evolutionary stages, we predict a decrease of 23%-46%.Naively, one would have expected dynamical interactions to boost the total number of systems experiencing mass transfer in the inner binary.In only 67.3% of the systems from the fiducial model the inner binary components undergo mass transfer, while for model Bin-Skew this is 80.7%.This discrepancy follows from the fact that there are a few additional types of interactions unique to triple stars (Fig. 11).The systems in the fiducial model that undergo these unique triple-star interactions migincludegraphics still undergo mass transfer in the inner binary at a later stage of evolution, but their contribution is excluded here as the simulations are terminated at the onset of the first interaction.Therefore, a direct comparison between the predictions of the binary and triple simulations would be undesirable.
However, the variations in the inner eccentricity due to threebody dynamics can give some insigincludegraphics into the difference in incidence rates of mass-transfer systems between the binary and triple simulations.For the triple population, the maximum oscillation amplitude of the inner eccentricity is most pronounced for the MS donors, with a measurable increase of more than 0.05 (0.01) in 19.1% (27.6%) of the systems.The MS lifetime is relatively long, which increases the opportunity for the quadruple terms as well as octupole terms to become relevant.At most other evolutionary phases the contribution of the tertiary star is significantly lower: 2.0% (6.6%) for HG donors, 0.27% (1.9%) for FGB donors, and 1.4% (6.0%) for AGB donors.The exception are CHeB donors, where the eccentricity increase is 15.4% (28.0%).Similar to the MS, the timescale of the CHeB phase is relatively long compared to the timescale of the other evolutionary phases.

Comparison with other work
Hamers & Thompson (2019) and Stegmann et al. (2022a) both study the intermediate and final evolutionary stages of massive hierarchical triple stars in the field using the respective triple stellar evolution codes TSE and SECULARMULTIPLE (Hamers & Portegies Zwart 2016;Hamers 2018).Toonen et al. (2020) addresses a similar scientific question as in this work, using the same code, but instead focus on the evolution of low mass triple stars (1 M < m 1 < 7.5 M ).We will compare the results of these studies with our results and comment on dissimilarities.
Hamers & Thompson (2019) investigate the merger rates of double neutron star and black hole-neutron star mergers in hierarchical triple systems.They sample the initial stellar and orbital properties from similar distributions, but the sampling boundaries differ somewhat from our study.First of all, the initial primary masses are sampled between 8 and 50 M , resulting in lower average masses.Second, their choice of the maximum initial inner orbital period is significantly longer compared to our study: 10 10 days opposed to 10 8.5 days.Stegmann et al. (2022a) investigate the evolution of a massive triple star population from the ZAMS until the formation of a compact object.Similar to Hamers & Thompson (2019), the initial primary masses reach down to 8 M .However, instead of sampling each initial property from an independent distribution, they use correlated distributions as presented in Moe & Di Stefano (2017).Both studies do not vary the initial properties, but do vary physics, such as the supernova kick model.

Mass transfer initiated by primary
Hamers & Thompson (2019) find that mass transfer initiated by the primary star occurs in 44% of systems when including a nonzero SN kick prescription.That is 20% fewer systems than the lower boundary on our predictions.The most likely explanation for this discrepancy is that their initial orbits extend to longer periods.Naturally, with less strongly bound orbits, we would expect a SN kick to unbind the system more easily.This is supported by an increase of over 20% in the number of system where the inner and outer orbit becomes unbound in Hamers & Thompson (2019).Furthermore, with lower typical masses, eventual SN kicks are stronger on average and are hence more likely to unbind the orbit.It is evident that changes in the initial orbital parameter ranges can significantly alter the interaction history of a triple star population.Stegmann et al. (2022a) find that at solar metallicity 74% of the systems undergo mass transfer in the inner binary.However, since they do not terminate the simulation at the onset of an interaction, this population includes systems that experience the unbinding of the outer orbit prior to the mass transfer event.If the tertiary star was originally more massive than either component of the inner binary, mass transfer in the inner binary could be achieved even after the SN kick of the tertiary star has unbound the outer orbit.Unsurprisingly, our results indicate that for most models a lower fraction of inner binaries initiate mass transfer, with a median of 67%.Since we do not have predictions for the number of systems that experience mass transfer in the inner binary after the first interaction, we refrain from drawing conclusions about the comparison between the incidence rates of mass-transfer systems.
The frequency of mass transfer episodes that are initiated by the primary star in low-mass hierarchical triple populations is similar as in our high mass population (Toonen et al. 2020).However, the mass transfer is dominated by more evolved donor stars, namely by red giant branch stars, while in our simulations the mass transfer is mainly initiated by donor stars on the MS or HG.Opposed to massive stars, low mass stars barely expand during the MS phase and avoid filling their Roche lobe.Additionally, low-mass stars leave the MS phase much closer to the Hayashi track, resulting in a shorter Hertzsprung gap phase.

Non-Interacting systems
Previous studies have shown that the formation rate of TCO systems that have not experienced an interaction are highly dependent on the assumed SN kick model (Antonini et al. 2017;Silsbee & Tremaine 2017;Rodriguez & Antonini 2018;Fragione & Loeb 2019).With lower kick velocities, the triple orbits are more likely to remain bound and the probability of producing a TCO becomes higher.With non-zero kick velocities, Stegmann et al. (2022a) and Hamers & Thompson (2019) predict a fraction of 0.12% and 0.4% non-interacting systems, respectively, which is within the range of our model predictions.However, when neglecting the SN kick, this fraction becomes substantially larger.Naively, one would expect that the compact-object merger rates would be drastically lower when higher kick velocities are implemented.Interestingly, Antonini et al. (2017) showed that this is not the case for binary black-hole mergers, as bound systems that have received a stronger SN kick have larger eccentricities on average, making three-body dynamics more effective, which results in shorter merger delay times.

GW sources
Studies exploring the effect of three-body dynamics on the merger rate of compact objects in non-interacting triple stars have often neglected three-body dynamics prior to TCO formation.We explore the consequences of neglecting three-body dynamics during the stellar evolution phase on the population of non-interacting triple stars at TCO formation by excluding the quadrupole and octupole terms.Consistent with other studies, we still include the stability criterion.Our results indicate that neglecting three-body dynamics during the stellar evolution phase has an impact on the properties of the TCO population.For example, the fraction of systems with relative inclinations between 40 • and 140 • (60 • and 120 • ) is 80% (51%) for the model variation compared to 70% (34%) for the fiducial model.The preference for aligned orbits is consistent across the other models.In addition, the fraction of systems in the model variation with outer eccentricities larger than 0.5 (0.6) is 59% (46%) compared to 45% (27%) for the fiducial population.However, not all other models show a clear preference towards lower outer eccentricities.The impact on the ZLK timescales is negligible, but the maximum amplitudes of the inner eccentricity due to threebody dynamical interactions are smaller for the fiducial population (Eq.3).Specifically, 62% of systems in the fiducial population have a maximum ZLK amplitude above 0.5, while this is 78% for the model variation.This suggests that previous studies ignoring three-body dynamics during stellar evolution have overestimated the rate of gravitational wave mergers resulting from non-interacting triple stars.These results are not surprising, as stellar systems with strong dynamics are expected to interact before TCO formation.Furthermore, we predict no discernible difference in the final masses, mass ratios or orbital periods of the inner binary.

Main caveats
Simulations with binary population synthesis codes rely on many assumptions that unavoidably introduce uncertainties to the predicted population.Adding a third star into the equation complicates matters further.In particular, the uncertainties in the initial orbital properties of the triples are prominent.Fortunately, these uncertainties are less important then one may expect (see Section 2.3).Besides the initial properties, we address the other main sources of uncertainty in the following section.

Massive star evolution
As mentioned in Section 2, the approximate nature of the single stellar evolution fitting formulae can lead to imprecise inferred final properties of the stars, especially for stars above 50M , as this is the upper mass limit of the Hurley stellar tracks (Hurley et al. 2000) and one needs to extrapolate at higher initial masses.
Furthermore, many aspects of stellar physics are still poorly understood.We highlight the two physical assumption most relevant for this study: mass loss throughout evolution and natal kick properties of compact objects.Uncertainties and intricacies in mass-loss from hot stars have recently been summarized by Vink (2022).Scatter in empirical rates as function of luminosity are typically at least a factor of three, while discrepancies between theory and observations may exceed that amount in certain parameter ranges (especially at luminosities below ∼ 10 5 L ; e.g., Brands et al. 2022).For cooler stars, uncertainties are even larger.For luminous stars that evolve toward the Humphrey-Davidson limit and experience a Luminous Blue Variable phase, accurate theoretical and observation-based massloss relations are lacking altogether (Smith 2014), forcing crude assumptions to be made.For Red Supergiants, which possibly suffer from multiple mechanisms causing the loss of gas (e.g., Cannon et al. 2021;Montargès et al. 2021), empirical prescriptions of mass loss deviate by up to an order of magnitude or more (e.g., Mauron & Josselin 2011;Beasor et al. 2020).These uncertainties are particularly concerning as stellar mass loss plays a crucial role in governing supernova progenitor properties, such as the core mass, and can affect the outcome of binary evolution.Higher wind mass-loss rates lead to more pronounced orbital widening and can prevent systems from engaging in mass transfer.
Of key relevance as well are the natal kick properties imparted on the compact object formed during a supernova event.
We have assumed the kick velocity model from Verbunt et al. (2017), which is based on the motions of young galactic pulsars, and scaled down for black holes according to the fallback prescription of Fryer et al. (2012) and the momentum of the object.However, the proper shape of the velocity distribution is still under debate.Some studies favour a single peaked Maxwellian distribution (Lyne & Lorimer 1994;Hansen & Phinney 1997;Hobbs et al. 2005), while others one that is double peaked (Fryer et al. 1998;Arzoumanian et al. 2002;Faucher-Giguère & Kaspi 2006;Verbunt et al. 2017).For black holes, the lack of observed natal kicks complicates discrimination between different models even further (Fragos et al. 2009;Repetto et al. 2012Repetto et al. , 2017;;Mandel 2016).The uncertainty in SN kicks on the occurrence rates of certain evolutionary channels can be significant (see Section 4.1.2).

Early interactions
Roughly 4% of the systems interact within 0.1 Myr from the ZAMS.These are mainly systems that engage in mass transfer and a small population of systems that become dynamically unstable.The vast majority of these systems has experienced an appreciable increase in the inner eccentricity, indicating threebody dynamical interactions have been at play.Only a handful of systems maintain a nearly constant inner eccentricity, as the systems were initialized extremely close to Roche-lobe filling or dynamical destabilisation.The number of systems that interact early on is on the order of 3-6% across all model variations, apart from the model i rel -const.At initial relative inclinations of zero (coplanar orbits), we predict only 0.8% of systems to interact within 0.1 Myr, since at low inclinations three-body dynamics is not effective.
We address a few points that could complicate the physical significance of these fast interactions.First, the initial inner semi-major axis of the systems experiencing mass transfer are typically not very large (few tens to few hundreds of solar radii).Before thermal equilibrium sets in on the ZAMS, the stellar radii are larger and the possibility can not be excluded that the components of the inner binary would interact on the pre-main sequence, likely resulting in a stellar merger (Tokovinin & Moe 2020).
Second, observed low velocity dispersion distributions in young stellar clusters suggest that the minimum separation of inner binary orbits at formation might be larger than assumed in this study (Sana et al. 2017;Ramírez-Tannus et al. 2021), and that the Sana distribution is only established ∼ 1 Myr after formation.With a larger minimum separation, we expect a decrease in the number of systems that interact early on, as systems with primary stars initially close to Roche-lobe filling are avoided and three-body dynamical effects need to be strong in order to drive wider systems to transfer mass.
Third, observations of solar mass triple stars suggest that tight inner orbits, a in < 2000 R , are associated with relative inclinations smaller than 40 • (Borkovits et al. 2016;Tokovinin 2017).At low inclinations, three-body dynamical interactions are strongly reduced and likely diminish the contribution of these fast interactions.

Eccentric mass transfer
An important challenge in population synthesis codes is the lack of understanding of how mass transfer progresses in eccentric orbits.Usually, this problem is circumvented by assuming the system circularises by the onset of mass transfer or initialising the orbits as circular.However, recent studies find that the tidal forces are generally inefficient in circularising the orbit of binary systems before mass transfer is initiated (Eldridge 2009;Vigna-Gómez et al. 2020;Vick et al. 2021).Moreover, in triple star systems, if the three-body dynamical interactions are strong enough, complete circularisation due to tidal forces can be prevented (e.g.Toonen et al. 2020).In Fig. 13, we show the degree of circularisation of the inner orbit for all systems in which the primary star fills its Roche lobe for our fiducial model.While many of the most compact inner orbits are completely circularised, the majority of systems retain a nearly constant eccentricity.A nonnegligible fraction of systems, located at small a out /a in even experiences a significant boost in eccentricity as three-body dynamical effects become more important.We stress that Fig. 13 is based on a singular tidal model, where (1) the tidal timescales are uncertain and (2) tidal processes, such as efficient dissipation in highly eccentric orbits (Moe & Kratter 2018;Generozov et al. 2018), are ignored.As a result, we are likely to underestimate the degree and timescale of circularisation in the inner binary.This is especially valid for inner orbits with high eccentricities, but also for small semi-major axes, where tides could be more efficient during the pre-main sequence phase.

Tertiary tides
Besides tidal forces between the primary and secondary star, the tertiary star can also dissipate energy from the inner binary through tidal interactions (Fuller et al. 2013).To investigate the importance of this effect, we apply the tertiary tidal model of Gao et al. (2020) on our fiducial population.As a result of tertiary tides, the semi-major axis of the inner binary shrinks as: with R 3 the radius and τ the viscoelastic relaxation time of the tertiary star.We set τ to 10 −4 yr in accordance with Gao et al. (2020).We did not find a noticeable impact of the tertiary tides on the evolution of our fiducial population for the chosen value of τ.

Conclusions
In this work, we have studied the main evolutionary pathways of massive hierarchical triple stars up to the first moment of interaction.Additionally, we have investigated what impact a tertiary companion has on the evolution of the system.We have done this by performing large-scale simulations of massive triple star evolution, beginning at the ZAMS up to the first point of stellar or orbital interaction.To account for uncertain formation properties in the simulated population, we included predictions for several model variations of the initial parameter distributions.Our key findings can be summarized as follows: -The initial parameters of stable hierarchical triple stars are highly constrained by the combined requirement of dynamical stability and avoidance of contact of the inner binary, which restrict the period distribution.Triple initiation results in less wide inner orbits, leading to a high incidence of interactions such as mass transfer.The inner orbit period distribution strongly differs from a flat distribution and deviates substantially different from a Sana distribution.-The vast majority of systems (65%-77%) have a phase of mass transfer initiated by the primary star as their first interaction.This occurs mainly at early evolutionary stages of the donor.In 32-50% of the cases the donor is still a mainsequence star, and it has evolved to the Hertzsprung gap in 38-50% of the systems.We have identified two main processes responsible for the large number of systems that initiate mass transfer at an early evolutionary phase.First, as a result of the dynamical stability criterion at formation, the inner orbits are typically more compact initially, and therefore the stars need to expand less in order to fill their Roche lobe, thus favoring mass transfer.Second, three-body dynamics can drive up the maximum eccentricity of the inner orbit, shrinking the eccentric Roche lobe.In the fiducial model, 10% of all systems that undergo mass transfer initiated by the primary star experience an increase of the inner eccentricity by at least 0.05.-Across all triple models, we predict that fewer than 0.5% of the systems do not engage in any interaction after being evolved for a Hubble time, such that all three stars effectively evolve as single stars.These systems are known targets as progenitors for gravitational wave mergers.We have shown that systems with strong three-body dynamics tend to interact before a triple compact object (TCO) is formed and hence reduce the likelihood of producing compact objects that merge with the help of von Zeipel-Lidov-Kozai (ZLK) oscillations.However, the typical ZLK timescales are still a few orders of magnitude shorter than a Hubble time and can accelerate the inspiral of the compact objects.We have also shown that ignoring three-body dynamics before compact object formation results in TCOs with stronger eccentricity oscillations and thereby likely over-predicts the merger rate of compact objects in such systems.-The predicted incidence rates of systems evolving through each evolutionary channel is most sensitive to variations in the inner/outer semi-major axis and the outer mass ratio distribution.The incidence rate of inner mass-transfer systems is dominant compared to other evolutionary channels leading to interaction.Though, the absolute differences between population models for these less common channels are not important, the relative differences can be relevant.For only two evolutionary channels, mass transfer from the tertiary star onto the inner binary and the unbinding of the outer orbit, the incidence rate differs by over a factor of two from the fiducial model's rate (see Fig. 4).-The evolution of an isolated massive binary star population up to the first interaction differs significantly from that of massive hierarchical triple stars.The criterion of dynamical stability and eccentricity oscillations due to dynamical inter-action ensure that at least an additional 15% of the population initiates mass transfer in the inner binary.Moreover, the donor stars are generally at an earlier evolutionary phase at the onset of mass transfer.

F
. Kummer et al.: The Main Evolutionary Pathways of Massive Hierarchical Triple Stars

Fig. 1 .
Fig. 1.Initial distributions of the inner-(blue) and outer (red) semimajor axis for our fiducial model and three model variations.The solid black lines represent the ZAMS population, consisting exclusively of initially stable, non-interacting hierarchical triple stars.The lighter colours correspond to the population before excluding dynamically unstable systems and systems that are Roche-lobe filling at birth.To enhance the features in the plot, we have cut off the top part of the distribution for model a in&out -Sana.The bottom panel shows the initial sampling distributions of the fiducial model and the Sana distribution.

Fig. 2 .
Fig. 2. Example of a system with large inner eccentricity oscillations induced by three-body dynamical effects.Top: for each evolutionary time step, the maximum increase in the inner eccentricity is shown.During the first few time steps the ZLK amplitudes are small, because at the start of the simulation the evolutionary time steps are very short.Bottom: the evolution of the eccentricity of the inner orbit during the first 0.5 Myr.

Fig. 3 .
Fig.3.Percentage of systems evolving through each evolutionary channel for our fiducial population of massive hierarchical triple stars.For the outer orbit unbound, tertiary mass transfer and dynamical destabilisation channels, the tertiary star is directly involved.

Fig. 4 .
Fig. 4. The predicted number of systems that evolve through each evolutionary channel for all model variations with respect to the fiducial model (blue diamonds).A factor 0.5 and 2 indicate a factor two increase and decrease, respectively.The predicted rates of all model variations are combined into the median (red stars), and the largest outliers per channel are represented by the error bars.The red part of the error bars correspond to Model q out -Moe, the model for which the predictions deviate most significantly from the fiducial model.The statistical uncertainties (3σ) on the predictions for the fiducial model are also included (orange bars).The relative contribution of each evolutionary channel for the fiducial are given in percentages.

Fig. 6 .
Fig. 6.The rate of interactions as a function of time for the fiducial model.Left: cumulative normalized distribution per evolutionary channel.The black dashed line combines the contribution of all channels.Right: The total number of interactions occurring at each moment in time, summed across all channels.

Fig. 7 .
Fig. 7.The incidence of pre-and post interacting massive triples assuming that stars formed at a constant star-formation rate.The contribution of binary stars are not included.

Fig. 8 .
Fig. 8. Evolutionary phase of the system at the onset of mass transfer in the inner binary.On the x-axis, the stellar types of the donor (primary) and accretor (secondary) are shown respectively.The symbols have the same meaning as those in Fig. 4.

Fig. 9 .
Fig.9.The circularization radius (the semi-major axis if the orbit was circularized by tides) of the inner and outer orbit of each system at the onset of mass transfer initiated by the primary star.Color-coded are the maximum amplitudes of the inner eccentricity induced by threebody dynamical interactions.Their distribution is also presented as a histogram.A value of ∆e max in = 1 indicates that during a single oscillation the inner eccentricity has varied between 0 and 1.

Fig. 10 .
Fig. 10.Cumulative step function of stellar and orbital properties for the fiducial model that are important in driving three-body dynamical interactions.In blue, we show the initial ZAMS conditions for all systems evolved in the simulation.In orange, the conditions at triple compact object formation are shown for the non-interacting systems.

Fig. 11 .
Fig. 11.Percentage of systems evolving through each evolutionary channel for the isolated binary models Bin and Bin-Skew, and the fiducial model for reference.The systems in model Bin are initialised and evolved as isolated binary stars.The systems in model Bin-Skew are initialised as hierarchical triple stars, but evolved as isolated binary stars after removing the tertiary star.

Fig. 12 .
Fig. 12.The number of systems per evolutionary phase that evolve through the (inner) mass transfer channel for the fiducial model (orange), the binary model with skewed initial parameters (blue) and the true isolated binary model (red).The hatched regions correspond to the number of inner binaries with maximum eccentricity amplitudes larger than 0.05 (black) and 0.01 (gray).

Fig. 13 .
Fig. 13.The circularization radius of the inner and outer orbit of each system at the onset of mass transfer initiated by the primary star.Colorcoded is the ratio between the final and initial eccentricity of the inner orbit.Their distribution is also presented as a histogram in the smaller plot.The light-colored systems (e fin in /e init in 1) at small circularization radii of the inner orbit have been completely circularized due to tidal dissipation.The dark-colored systems (e fin in /e init in > 1) at comparable inner-and outer circularization radii have been subjected to three-body dynamical effects.The evolution of the other systems has hardly been affected by either tidal effects or three-body dynamics.

Table 1 .
The initial sampling ranges and distributions of the triple star properties for the fiducial model and all model variations.

Table 2 .
Overview and classification of the different evolutionary channels considered in this study.