© The Authors 2023

1. Introduction

Massive stars (M ≳8 M_⊙) play an important role in the dynamical and chemical evolution of galaxies. Stellar feedback in the form of stellar winds, radiation pressure, photo-ionisation, and supernova (SN) explosions controls star formation (e.g. Hayward & Hopkins 2017) and enriches 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 in particular 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 these can more easily be explained with the inclusion of a third companion; for example, 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. 2016, 2019; Shultz et al. 2019). Therefore, complimentary population synthesis studies of triples are essential in order to create 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, three-body 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 orbital plane of the third star with respect to the orbital plane of the inner orbit. For an elaborate review of 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 or unbound, or at the onset of mass transfer, that is, 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 Sect. 2, we give a description of the triple population synthesis code and motivate our choice of initial model conditions. In Sect. 3, we present the predicted incidence rates of each evolutionary channel and discuss a few channels in more detail. In Sect. 4, we compare our results with similar studies and discuss some caveats. Finally, we give a summary of our findings in Sect. 5.

2. Method

To investigate the most common evolutionary channels in massive hierarchical triple stars, we use the triple population synthesis code TRES¹ (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 –, with 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 centre of mass of the inner two stars, forming the outer orbit. In these configurations, the energy of the outer and inner orbits is separately conserved (Heggie 1975). To model stellar evolution, SeBa makes use of analytic fitting formulae from Hurley et al. (2000, 2002), which are based on the single stellar evolution models produced by Pols et al. (1995). These fitting formulae allow for rapid simulation of stellar evolution, which is essential for population synthesis studies.

With TRES, we simulate ten sets of 25 000 initially dynamically stable, coeval, massive hierarchical triple stars that are initialised 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 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.

2.1. Updated wind prescription

By default, SeBa follows the line-driven wind mass-loss model for massive stars of Vink et al. (2000, 2001) in the metallicity range 0.03 Z_⊙  <  Z  <  3 Z_⊙ and those of Nieuwenhuijzen & De Jager (1990) in the range where the Vink models are not applicable. We modified these prescriptions by reducing the mass-loss rates by a factor of 3, following and generalising the findings of Björklund et al. (2021) for the implemented mass range. Additionally, we assume luminous blue variable mass-loss rates following Belczynski et al. (2010), introducing a constant enhanced mass loss of 1.5 × 10⁻⁴ M_⊙ yr⁻¹ for stars exceeding a luminosity of 6 × 10⁵ L_⊙.

2.2. Initial parameters

To simplify the description of a hierarchical triple star, the system is split up into two subsystems: Firstly, we define an inner binary, which comprises the primary and secondary star with masses m₁ and m₂. These stars orbit each other with an inner semi-major axis a_in, eccentricity e_in, and argument of pericentre g_in. Secondly, we have an outer binary, which comprises the centre of mass of the inner binary and the tertiary star, with mass m₃. The outer binary has an outer semi-major axis a_out, eccentricity e_out, and argument of pericentre 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₁ 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₂ and m₃ are specified through the inner mass ratio q_in ≡ m₂/m₁ and the outer mass ratio q_out ≡ m₃/(m₁ + m₂), respectively. Both the inner and outer mass ratios are sampled from a flat distribution in the range [0.1,1]², which is based on mass-ratio observations of massive binary stars (Sana et al. 2012; Kobulnicky et al. 2014). Although the tertiary mass is generally lower than the mass of the inner binary, there are exceptions (e.g. Eisner et al. 2022). The orbital semi-major axis ranges for a_in and a_out both cover [5,5 × 10⁶] 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⁶ 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, 2010). When the sampling from our orbital semi-major axis distribution function yielded a_in > a_out, we decided to swap the values, as, by definition, a_in < a_out. This is a different approach from one of the models in Toonen et al. (2020) for example, where the authors fix the primordial distribution of a_in and accept or reject a_out in accordance with the stability properties of the system. The resulting ZAMS distributions of a_in and a_out might differ between the two methods. The eccentricities e_in and e_out are sampled from a thermal distribution in the range [0,0.95] (Ambartsumian 1937; Heggie 1975; Moe & Di Stefano 2017; Hwang 2023). The arguments of pericentre g_in and g_out are sampled from a uniform distribution in the range [−π,π]. The relative inclination i_rel is sampled from a uniform distribution in cosine in the range [0,π]. Lastly, the initial stellar rotational velocities are drawn from a distribution as described in Hurley et al. (2000). We refer to Toonen et al. (2016) for formulations of the effects of stellar rotation on the tidal forces and orbital precession in TRES.

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), which are 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₁₀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⁴ 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.

2.3. Dynamical stability

To ensure the long-term stability of a three-body system, a hierarchical configuration is generally required, that is, 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:

$$\begin{array}{c}\hfill \frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}{|}_{\mathrm{crit}}=\frac{2.8}{1-e_{\mathrm{out}}}(1-\frac{0.3i_{\mathrm{rel}}}{\pi })[\frac{(1+q_{\mathrm{out}})(1+e_{\mathrm{out}})}{\sqrt{1-e_{\mathrm{out}}}}{]}^{2/5}.\end{array}$$(1)

A system is dynamically stable if a_out/a_in > (a_out/a_in)_crit. For instance, for an equal-mass system (q_out = 0.5) with a circular orbit, 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. Figure 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 relatively similar to the fiducial model, even though their sampling distributions are vastly different.

Fig. 1. Initial distributions of the inner- (blue) and outer (red) semi-major 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 ain & aout-Sana. The bottom panel shows the initial sampling distributions of the fiducial model and the Sana distribution.

2.4. Three-body dynamics

2.4.1. Secular effects

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 behaviour 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_{\mathrm{out}}^2/P_{\mathrm{in}}$, and e_out. An approximate expression for this timescale is (Kinoshita & Nakai 1999; Antognini 2015):

$$\begin{array}{c}\hfill {\tau }_{\mathrm{ZLK}}\approx \frac{8}{15\pi }(1+\frac{m_1}{m_3})\left(\frac{P_{\mathrm{out}}^2}{P_{\mathrm{in}}}\right){(1-e_{\mathrm{out}}^2)}^{3/2}.\end{array}$$(2)

In the test-particle approximation, large oscillations of the inner eccentricity occur only at initial relative inclinations i_init between 39.2° and 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):

$$\begin{array}{c}\hfill e_{\max }=\sqrt{1-\frac{5}{3}{\cos }^2(i_{\mathrm{init}})}.\end{array}$$(3)

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, three-body dynamics can provoke a flip in relative inclination in this regime (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,

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{oct}}=\frac{m_1-m_2}{m_1+m_2}\frac{a_{\mathrm{in}}}{a_{\mathrm{out}}}\frac{e_{\mathrm{out}}}{1-e_{\mathrm{out}}^2},\end{array}$$(4)

is of the order |ϵ_oct|≳0.001–0.01 (Lithwick & Naoz 2011; Shappee & Thompson 2013).

2.4.2. Eccentricity oscillation tracking

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. Figure 2 shows 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₁ = 20 M_⊙, m₂ = 15 M_⊙, and m₃ = 30 M_⊙, initial semi-major axes a_in = 2 × 10⁴ R_⊙ and a_out = 1.5 × 10⁵ 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 in the inner eccentricity by up to 0.8. As the octupole parameter of this system is 0.009, the octupole term is important, which manifests as 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.

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: Evolution of the eccentricity of the inner orbit during the first 0.5 Myr.

2.5. Tidal interaction

Stars in a binary system are subjected to tidal forces that act to synchronise the stellar rotation with the orbital period and circularise 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 parameterised 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 pericentre passage (e.g. Fabrycky & Tremaine 2007). After complete circularisation, the semi-major axis has decreased by a factor of $(1-e_{\text{in}}^2)$.

3. 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 criterion 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 core-collapse 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 merger of two stars as a stopping condition, as such events are preceded by mass transfer. We note that we do not include tidal effects or three-body dynamical effects such as ZLK oscillations as stopping conditions, while, strictly, 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 Sect. 3.1 we give an overview of the predictions of the evolutionary channels for our simulated populations. We then discuss two channels in more detail; the systems that experience a phase of mass transfer initiated by the primary star in Sect. 3.2 and the systems that do not undergo any interaction during their entire evolution in Sect. 3.3. Next, we dive into a few of the model variations in more detail in Sect. 3.4. Finally, we compare the predictions for massive triple star evolution with simulations of isolated binary stars in Sect. 3.5.

3.1. Overview evolutionary channels

Figure 3 shows 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%) experience 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 uncertainties follow a Poisson distribution by approximation. A relatively large fraction of systems lose 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)³. In cases where 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 fraction of a few percent of the systems experience 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.

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.

3.1.1. Model variations

In addition to statistical uncertainties, the predicted contributions of evolutionary channels are dependent on the uncertainties related to the formation of the stellar system; that is, initial stellar and orbital parameters. Figure 4 shows the predicted ranges of the contributing fraction of each evolutionary channel resulting from the model variations discussed in Sect. 2.2 relative to the fiducial model. Overall, the incidence rate of each evolutionary channel is relatively robust across the model variations; the rates only twice differ from the fiducial model by more than a factor 2. The most significant outliers mainly originate from one model variation, Model q_out-Moe. We discuss some of the outliers in more detail in Sect. 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 or 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.

Fig. 4. 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 decrease and increase, 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 corresponds to Model qout-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 is given as a percentage.

In conclusion, the general evolution of a massive triple star population leading up to the first interaction is not seriously affected by the differing assumptions on the initial conditions. This can partially be explained by the initial stability check described in Sect. 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 is identical in both populations. We terminate the simulation at the onset of the first interaction, but the outcome of certain interactions, such as mass transfer, and therefore 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.

3.1.2. 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.

Fig. 5. Initial inner and outer orbital periods covered by each evolutionary channel. The triangular shape of the population is a direct consequence of the requirement on hierarchical structure and dynamical stability. The black dashed line shows where Pin = Pout. There are no systems near this line because of the dynamical stability criterion. There is some overlap between colours, as the evolution of a system is determined by more properties than just the orbital periods.

– 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⁴ days, extreme inner eccentricities induced by ZLK oscillations can provoke the primary star to fill its Roche lobe during pericentre passage. For a few systems, the mass transfer 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, and therefore 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 of larger than ≳10⁴ 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³ and 10⁵ − 10⁶ days, respectively, for this channel.

– Orbital unbinding: The initial conditions for systems whose inner (outer) orbit becomes unbound as a result of a SN kick are restricted to inner (outer) orbital periods of larger than a few thousand days. Their subsequent evolution is difficult to predict without proper n-body calculations, as a phase of non-secular 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 additional parameters, 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³ days. Smaller periods would most likely result in a phase of mass transfer within the inner binary.

3.1.3. Time-evolution of evolutionary channels

The stellar and orbital properties of a system are affected when interactions between stars or orbits take place. As 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.

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

At very early times (≲4 Myr), we predict only dynamical destabilisation of the triple and mass transfer in the inner binary. The extent 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 rates of interactions are low.

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. At 15 Myr, 78% and 46% of the systems that unbind the outer orbit and the inner orbit, respectively, have interacted. 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 stars with the highest masses 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₁ has a sharp cut-off at 10 M_⊙, while m₃ 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 increases to approximately 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 would like to know what the triple population looks like for a realistic stellar population. We can gain insight here by convolving the results from Fig. 6 with a star formation history. Figure 7 demonstrates the appearance of the intrinsic population of massive triples 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%) have 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 emphasises the importance of mass transfer in multiple stars. A significantly smaller fraction of the systems (15.9%) have 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 one-third of post-mass-transfer binary systems result in a merger. Our findings may therefore indicate that a large fraction of systems born as massive triple stars will be reduced to a population of single and binary stars at the moment of observation.

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

3.2. Mass transfer initiated by the primary

The majority of systems (64.6–76.6%) have 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) phases. 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%). As 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 therefore similar initial masses. Nonetheless, these systems can lead to interesting events, such as a double-core common envelope (Dewi et al. 2006).

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.

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 mentioned above, not all stars evolve through the FGB phase. Over 75% of the primary stars are initially more massive than about 13 M_⊙ and move directly to the CHeB phase. Similarly, stars initially more massive than about 38 M_⊙ 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, and therefore reaches the onset of mass transfer early in its evolution. Second, three-body dynamical interactions can expedite the onset of mass transfer by decreasing the pericentre 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 by 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 of between 0.05 and unity, respectively. 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 100, because additional parameters (e_out, q_out and i_rel) affect the timescale and magnitude of the oscillations.

Fig. 9. Circularisation radius (the semi-major axis if the orbit was circularised by tides) of the inner and outer orbit of each system at the onset of mass transfer initiated by the primary star. Colour-coded are the maximum amplitudes of the inner eccentricity induced by three-body dynamical interactions. Their distribution is also presented as a histogram. A value of $\mathrm{\Delta }{e}_{\mathrm{in}}^{\max }=1$ indicates that, during a single oscillation, the inner eccentricity has varied between 0 and 1.

3.3. 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 its entire evolution. Above, we show that these systems have initially relatively wide inner orbits, preventing a phase of mass transfer, and are therefore also required to have wide outer orbits to ensure a stable hierarchical configuration. However, at these wide separations, the weakly bound tertiary star is prone to becoming unbound because of a SN kick, justifying the small contribution of this channel.

Nevertheless, systems devoid of any interaction are an interesting target for 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 2017; Rodriguez & Antonini 2018; Fragione & Loeb 2019; Martinez et al. 2022). In this scenario, the inner binary forms a double compact object (DCO) that is generally too wide to merge within a Hubble time. The tertiary star can drive the inner two objects significantly closer through three-body dynamics over long timescales, eventually merging the inner binary. Alternatively, in systems with a very wide outer orbit (> 1000 AU), flybys can induce a gravitational wave merger by exciting the eccentricity of the outer orbit (Michaely & Perets 2020).

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 of the non-interacting systems at TCO formation for the fiducial model. Most of the properties show a relatively dissimilar distribution between the two populations. The TCO population evidently favours more moderate inner- and outer eccentricities. Highly eccentric systems are more likely to interact and hence prevent the formation of a TCO that has had no interaction during its evolution. Interestingly, the models e_in-Sana and e_out-Sana, which 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₃/m₁) are strongly favoured. 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⁴ for the TCO systems. The minimum inner orbital period to avoid mass transfer in the inner binary is around 6 × 10³ 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 disfavour relative inclinations of 90°. At these inclinations, the eccentricity oscillations in the inner binary are the largest, complicating the formation of a non-interacting 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).

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.

3.4. 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 favour 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 rates of 73.4% and 75.1%, respectively, compared to the fiducial model rate of 67.3  ±  0.9%. These models have more compact inner orbits, resulting in smaller Roche radii of the inner binary components, 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 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 that unbind 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 SN 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 SN 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 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 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 SN 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 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 pericentre 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 MS donors increases to 49.3% and 51.7%, respectively, compared to the fiducial model value of 38.1%. However, there are slightly 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 pericentre of the inner binary does not decrease, and the onset of mass transfer is expedited.

3.5. 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 found 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; that is, the systems are initiated and evolved without the presence of a tertiary star. We refer to this model as 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 pericentre at the ZAMS for the systems from model Bin-Skew. Interestingly, these differences are much bigger than the differences between the triple star models 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.

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.

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 might expect 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 did 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 might 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.

Fig. 12. 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 (grey).

However, the variations in the inner eccentricity due to three-body dynamics can provide some insight 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 increase in stellar radius during the MS is marginal compared to that in the post-MS evolution. Therefore, eccentricity oscillations are expected to play an important role in driving mass transfer for MS donors. Additionally, the MS lifetime is relatively long, which increases the opportunity for the quadruple and 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 and for a large range of initial masses the radial expansion is minor.

4. Discussion

4.1. 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) address a similar scientific question to that in the present work, using the same code, but instead focus on the evolution of low-mass triple stars (1 M_⊙ < m₁ < 7.5 M_⊙). We 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. These authors 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¹⁰ days as 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, these latter authors 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 SN kick model.

4.1.1. 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 non-zero SN kick prescription; that is 20% fewer systems than the lower boundary on our predictions. The most likely explanation for this discrepancy is that the initial orbits used by these latter authors extend to longer periods. Intuitively, with less strongly bound orbits, we would expect a SN kick to unbind the system more easily. This is supported by the findings of Hamers & Thompson (2019), who show an increase of over 20% in the number of systems where the inner and outer orbit become unbound. Furthermore, with lower typical masses, eventual SN kicks are stronger on average and are therefore 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, as these authors 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 is 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%. As 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 to 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. As 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 HG phase.

4.1.2. Non-interacting systems

Previous studies have shown that the formation rate of TCO systems that have not experienced an interaction is 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 fractions of 0.12% and 0.4% non-interacting systems, respectively, which are within the range of our model predictions. However, when neglecting the SN kick, this fraction becomes substantially larger. Naively, one would expect the compact-object merger rates to 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.

4.1.3. 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 of 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 three-body 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.

4.2. 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 than one may expect (see Sect. 2.3). Besides the initial properties, we address the other main sources of uncertainty in the following section.

4.2.1. Massive star evolution

As mentioned in Sect. 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 50 M_⊙, as this is the upper mass limit of the Hurley stellar tracks (Hurley et al. 2000) and one needs to extrapolate to higher initial masses.

Furthermore, many aspects of stellar physics are still poorly understood. We highlight the two physical assumptions most relevant to this study: mass loss throughout evolution and natal kick properties of compact objects. Uncertainties and intricacies in mass loss from hot stars were recently summarised by Vink (2022). Scatter in empirical rates as a 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⁵ L_⊙; e.g. Brands et al. 2022); for cooler stars, uncertainties are even larger. For luminous stars that evolve towards the Humphrey-Davidson limit and experience a luminous blue variable phase, accurate theoretical and observation-based mass-loss 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 SN event. We assume 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. 2012, 2017; Mandel 2016). The uncertainty in SN kicks on the occurrence rates of certain evolutionary channels can be significant (see Sect. 4.1.2).

4.2.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 have experienced an appreciable increase in the inner eccentricity, indicating three-body dynamical interactions have been at play. Only a handful of systems maintain a nearly constant inner eccentricity, as the systems were initialised 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, because three-body dynamics is not effective at low inclinations.

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 cannot be excluded that the components of the inner binary would interact on the pre-MS, 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 that 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.

4.2.3. 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 by 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 non-negligible fraction of systems located at small a_out/a_in even experience 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-MS phase.

Fig. 13. Circularisation radius of the inner and outer orbit of each system at the onset of mass transfer initiated by the primary star. The colour code refers to 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-coloured systems ($e_{\mathrm{in}}^{\mathrm{fin}}/e_{\mathrm{in}}^{\mathrm{init}}\ll 1$) at small circularisation radii of the inner orbit have been completely circularised due to tidal dissipation. The dark-coloured systems ($e_{\mathrm{in}}^{\mathrm{fin}}/e_{\mathrm{in}}^{\mathrm{init}}>1$) at comparable inner and outer circularisation 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.

4.2.4. 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:

$$\begin{array}{cc}\hfill \frac{1}{a_{\mathrm{in}}}\frac{\mathrm{d}{a}_{\mathrm{in}}}{\mathrm{d}t}& =2.22\times {10}^{-8}{\mathrm{yr}}^{-1}\frac{4q_{\mathrm{in}}}{{(1+q_{\mathrm{in}})}^2}(\frac{R_3}{100{\mathrm{R}}_{\odot }}{)}^{5.2}\hfill \\ \hfill & \times (\frac{a_{\mathrm{in}}}{0.2\mathrm{AU}}{)}^{4.8}(\frac{a_{\mathrm{out}}}{2\mathrm{AU}}{)}^{-10.2}(\frac{\tau }{0.534\mathrm{yr}}{)}^{-1.0},\hfill \end{array}$$(5)

where R₃ is the radius and τ the viscoelastic relaxation time of the tertiary star. We set τ to an optimistically small value of 10⁻⁴ yr. We did not find a noticeable impact of the tertiary tides on the evolution of our fiducial population for the chosen value of τ. We expect tertiary tides to play a more significant role in the post-mass-transfer population.

5. Conclusions

In this work, we studied the main evolutionary pathways of massive hierarchical triple stars up to the first moment of interaction. Additionally, we investigated the impact a tertiary companion has on the evolution of the system. To this end, we performed 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 summarised 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 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 cases, the donor is still a MS star, and has evolved to the Hertzsprung gap in 38%–50% of the systems. We 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, which favours 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 in 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 also show that ignoring three-body dynamics before compact object formation results in TCOs with stronger eccentricity oscillations and thereby likely leads to over-prediction of 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, that is, 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 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 interaction 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.

In conclusion, the evolution of massive hierarchical triple stars includes a wealth of interactions, and differs from the evolution of isolated massive binaries on a broader population level. Future studies should focus both on formation pathways and rates of specific evolutionary (end) products of triple star systems.

Acknowledgments

The authors acknowledge support from the Netherlands Research Council NWO (VENI 639.041.645 and VIDI 203.061 grants). The data and scripts necessary to reproduce the figures in this paper are publicly available on Zenodo: https://zenodo.org/record/8183545.

References

Appendix A: Overview incidence rates

All Tables

All Figures

Fig. 1. Initial distributions of the inner- (blue) and outer (red) semi-major 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 ain & aout-Sana. The bottom panel shows the initial sampling distributions of the fiducial model and the Sana distribution.

In the text

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: Evolution of the eccentricity of the inner orbit during the first 0.5 Myr.

In the text

	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.
In the text

Fig. 4. 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 decrease and increase, 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 corresponds to Model qout-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 is given as a percentage.

In the text

Fig. 5. Initial inner and outer orbital periods covered by each evolutionary channel. The triangular shape of the population is a direct consequence of the requirement on hierarchical structure and dynamical stability. The black dashed line shows where Pin = Pout. There are no systems near this line because of the dynamical stability criterion. There is some overlap between colours, as the evolution of a system is determined by more properties than just the orbital periods.

In the text

	Fig. 6. Rate of interactions as a function of time for the fiducial model. Left: Cumulative normalised distribution per evolutionary channel. The black dashed line combines the contribution of all channels. Right: Total number of interactions occurring at each moment in time, summed across all channels.
In the text

	Fig. 7. Incidence of pre- and post interacting massive triples assuming that stars formed at a constant star-formation rate. The contribution of binary stars is not included.
In the text

	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.
In the text

Fig. 9. Circularisation radius (the semi-major axis if the orbit was circularised by tides) of the inner and outer orbit of each system at the onset of mass transfer initiated by the primary star. Colour-coded are the maximum amplitudes of the inner eccentricity induced by three-body dynamical interactions. Their distribution is also presented as a histogram. A value of $\mathrm{\Delta }{e}_{\mathrm{in}}^{\max }=1$ indicates that, during a single oscillation, the inner eccentricity has varied between 0 and 1.

In the text

	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.
In the text

	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.
In the text

	Fig. 12. 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 (grey).
In the text

Fig. 13. Circularisation radius of the inner and outer orbit of each system at the onset of mass transfer initiated by the primary star. The colour code refers to 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-coloured systems ($e_{\mathrm{in}}^{\mathrm{fin}}/e_{\mathrm{in}}^{\mathrm{init}}\ll 1$) at small circularisation radii of the inner orbit have been completely circularised due to tidal dissipation. The dark-coloured systems ($e_{\mathrm{in}}^{\mathrm{fin}}/e_{\mathrm{in}}^{\mathrm{init}}>1$) at comparable inner and outer circularisation 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.

In the text