© The Authors 2025

1. Introduction

Massive stars (M ≥ 8 M_⊙) are critical to our understanding of the Universe. During their short but dramatic evolution, they synthesize heavy elements, drive stellar winds and emit intense ultraviolet radiation. Their explosive deaths as supernovae release vast amounts of energy and enrich the interstellar medium of their nucleosynthetic products. These contributions make them play a pivotal role in the evolution of galaxies and the Universe (e.g., Woosley et al. 2002; Crowther 2012; Hayward & Hopkins 2017).

Recent observations indicate that a large fraction of stars possess one or several companions, and that this fraction increases with stellar mass (e.g., Sana et al. 2012, 2014; Moe & Di Stefano 2017; Offner et al. 2023). A significant number of these systems are likely to interact during their evolution (Sana et al. 2012). The fraction of stars in multiple systems increases with the mass, as does the multiplicity order (number of companions). In particular, the observed fraction of stars in triple systems among early B- and O-type stars is larger than 30% (Moe & Di Stefano 2017).

Although the physics of the evolution of binary systems has received extensive attention over the last two decades (e.g., Hurley et al. 2002; Langer 2012; de Mink et al. 2013; Paxton et al. 2015; Belczynski et al. 2016; Marchant et al. 2016; Eldridge et al. 2017; Fragos et al. 2023; Bavera et al. 2023; Andrews et al. 2024), less is currently known about triple systems.

Recent studies have focused on the evolution of hierarchical triple star systems (e.g., Naoz et al. 2016; Toonen et al. 2016, 2020, 2022; Hamers et al. 2021; Stegmann et al. 2022; Preece et al. 2022, 2024; Kummer et al. 2023; Dorozsmai et al. 2024; Generozov & Perets 2024; Bruenech et al. 2025; Shariat et al. 2025; Vigna-Gómez et al. 2025), where the tertiary star is situated at far enough distances from the inner binary to maintain system stability over timescales on the order of or longer than the evolutionary timescale of the stars. One of the main results of these population synthesis studies is that the presence of a third object tends to increase the rate of interactions of the inner binary by triggering von Zeipel-Lidov-Kozai oscillations (ZLK, von Zeipel 1910; Lidov 1962; Kozai 1962). Notably, Kummer et al. (2023) proposed statistical predictions of the main evolutionary pathways of massive hierarchical triple star systems following a population synthesis approach.

The large majority of triple system simulation studies have been performed utilizing rapid stellar codes based on the fitting formulae by Hurley et al. (2000, 2002, hereafter H00, H02). This family of fast stellar codes notably includes SSE (Hurley et al. 2000), SEBA (Portegies Zwart & Verbunt 1996; Toonen et al. 2012), COSMIC (Breivik et al. 2020), COMPAS (Riley et al. 2022), and STARTRACK (Belczynski et al. 2002). One exception is the study by Portegies Zwart & van den Heuvel (2016), who coupled a detailed stellar code to an N-body code to perform simulations of hierarchical triple systems. Despite their outstanding computing performance, fast stellar codes are affected by simplifying assumptions, and some processes are not treated self-consistently (e.g., wind mass-loss and mass transfer). The limitations of the fast stellar codes have been known since their conception, but less well studied is the impact of these simplifications on the evolution of multiple systems. Bavera et al. (2023) showed that by not self-consistently modeling the reaction of the star to mass-loss by stellar winds, fast stellar codes can in some cases overpredict the maximal radial expansion of massive stars by several orders of magnitude, which can have dramatic consequences for the predicted evolution of binary systems. The impact of the discrepancies between stellar tracks predicted by detailed and rapid evolutionary codes on the evolution of triple systems has received little attention.

Although there have recently been new approaches and alternatives to the H00 formulae, which have proven successful in overcoming some of their limitations (e.g., Spera et al. 2015; Eldridge et al. 2017; Kruckow et al. 2018; Agrawal et al. 2023; Fragos et al. 2023), these methods have not yet been applied in triple system population synthesis studies, except in the very recent work by Shariat et al. (2025). In this study, they couple a triple secular code to the POSYDON tracks (Fragos et al. 2023) to investigate the origin of the V404 low-mass X-ray binary. They report differences between POSYDON and the SSE stellar evolution tracks and discuss the impact of these differences on their triple system simulations.

For this study, we coupled the TRiple Evolution Simulation code (TRES, Toonen et al. 2016) to the Modules for Experiment in Stellar Astrophysics (MESA, Paxton et al. 2011, 2013, 2015, 2018, 2019) to perform the first simulations of triple systems combining a triple secular evolutionary code with a detailed on-the-fly stellar code. Our approach differs from that of Shariat et al. (2025) in that in our simulations the detailed stellar evolution is computed simultaneously to the three-body dynamics. This allows for alterations of the stellar physics on a per system basis and provides access to the full structure of the stellar components, enabling the self-consistent modeling of key interactions specific to triple systems (in particular the precession caused by the distortion of the stellar components and its competition with the precession induced by the ZLK mechanism).

We primarily focused on stars with masses M ≥ 50 M_⊙ as they are frequently found in triple configurations and because above this mass threshold the H00 tracks extrapolate the fitting formulae obtained from the Pols et al. (1998, hereafter P98) grid. We followed an approach similar to that used in Bavera et al. (2023), aiming to illustrate the importance of a more self-consistent stellar evolution modeling for predicting the evolutionary pathways of massive triple systems. To this end, we compared the outcome of our simulations to analogous simulations performed with a fast stellar code (SEBA, Portegies Zwart & Verbunt 1996; Toonen et al. 2012).

The structure of this paper is as follows. In Sect. 2 we provide a description of the methods we used and the assumptions of single and triple star physics we made. In Sect. 3 we compare the single star evolution provided by the two different stellar codes we used. In Sect. 4 we present simulations of a set of triple systems for which the two codes predict divergent evolutionary pathways. Finally, we summarize our findings in Sect. 5.

2. Methods

To investigate the evolutionary pathways of massive hierarchical triple systems, we used the TRiple Evolution Simulation code (TRES, Toonen et al. 2016). TRES couples a stellar evolution code to a secular evolution code (SECULARTRIPLE, Toonen et al. 2016) that solves the three-body dynamics of the system following the secular approximation. In this approach, the perturbations occasioned by the third body are assumed to act on timescales much longer than the dynamical timescale of the inner binary, so that only the long-term (secular) evolution is considered, following a perturbative method. By default, the stellar evolution in TRES is handled by SEBA, a rapid population synthesis code based on the fitting formulae of H00.

TRES uses the Astrophysical MUltipurpose Software Environment (AMUSE, Portegies Zwart et al. 2009, 2013; Pelupessy et al. 2013), a library of codes that allows the user to easily exchange between codes according to the necessity of the modeling.

MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) is a 1d detailed stellar evolution code, which is extensively used in the field of single star and binary systems modeling. MESA is interfaced in AMUSE in its version r15140. Utilizing the MESA interface in AMUSE, we coupled MESA to TRES, which allowed us to perform state-of-the-art, detailed simulations of hierarchical triple systems. The code developments required to make TRES fully compatible with MESA have been made publicly available¹.

The approach chosen in this study is to use the versatility of AMUSE to compare the evolutionary outcome of triple systems obtained with a rapid code (SEBA) to that obtained with a detailed code (MESA). For all the presented systems, the evolution was computed both with SEBA and MESA, keeping the exact same physics for the binary and triple interactions.

2.1. Connecting TRES and MESA

The connection between TRES and MESA was done through the interface of MESA in AMUSE. In AMUSE, MESA is treated as a library. Parameters are set with a set_control method, instead of the usual inlist file. Additional getters and setters functions are available, which allow communications with the MESA instances.

Some binary and triple interactions implemented in TRES depend on quantities directly related to the structure of the stars, such as the apsidal motion constant k_AMC and the gyration radius r_g (see Sect. 2.3 for more details). When using TRES with SEBA, these quantities are given rough values which depend mainly on the evolutionary phase of the star. For instance, the apsidal motion constant is assumed to be k_AMC = 0.0144 for main sequence (MS) stars, helium-MS stars, and white dwarfs, k_AMC = 0.260 for neutron stars, k_AMC = 0 for black holes (BHs) and k_AMC = 0.05 for all other types of stars, following the prescriptions of Brooker & Olle (1955), Claret & Gimenez (1992)². In reality, k_AMC and r_g vary throughout the evolution of the stars (e.g., Claret & Gimenez 1989; Claret 1999; Claret & Giménez 2010; Rosu et al. 2020), see also Appendix C. When using TRES with MESA, these quantities are directly retrieved from the structure, which offers a more self-consistent treatment of the corresponding interactions. For that purpose, we added the methods get_apsidal_motion_constant and get_gyration_radius to the MESA interface in AMUSE.

In TRES, the binary and triple interactions are handled by SECULARTRIPLE. For this reason, we used the single star module of MESA instead of the binary module. All the interactions are thus handled by TRES. As in Kummer et al. (2023), we did not follow in detail the mass transfer (MT) events but stopped the simulations whenever a star experiences Roche lobe overflow (RLOF). Therefore in the presented simulations, the stars’ evolution affect the orbital evolution of the system, but the feedback from the interactions to the stellar evolution are not accounted for. As a result, the stellar evolution is identical to that of isolated single stars. The list of stopping conditions are given in Sect. 2.3, together with a description of the binary and triple interactions.

2.2. Single star physics

For the comparisons of the evolutionary pathways obtained with MESA and SEBA to be pertinent, it is crucial to use the two codes with equivalent stellar physics ingredients. We sought to follow this approach as precisely as possible. Some stellar physics ingredients however could not be kept exactly consistent between the two codes, due to intrinsic differences between rapid and detailed codes. In this section we discuss our choice of input physics and disparities between the two stellar codes. Rapid codes such as SEBA are based on the formulae of H00 fitted to the tracks of P98. The P98 grid covers the mass range M ∈ [0.5,  50] M_⊙. Although nothing prevents the use of the fitting formulae for stars outside this range of mass, they get increasingly divergent from tracks computed with a detailed code for increasing masses, as illustrated in Sect. 3.1. The grid of P98 was computed without accounting for mass-loss by stellar winds. In their paper, H00 provide receipes for incorporating mass-loss a posteriori in rapid codes. In their approach, after a star has lost a certain amount of mass due to stellar winds, it shifts to the track of a star of its current mass. The reaction to mass-loss is therefore not self-consistent. In our simulations, we sought to use equivalent winds prescriptions with the two codes, and show that an inconsistent treatment of the reaction to mass-loss can imply divergent evolutionary pathways for certain triple systems. The reactions to mass-loss obtained with the two codes are compared in Sect. 3.2.

2.2.1. Metallicity

We simulated systems at solar metallicity. The grid of P98 was computed with the old value Z_⊙ = 0.02, which we adopted for all simulations. This is also the default value of the solar metallicity in MESA (in particular for the winds scaling, see Sect. 2.2.3 in Andrews et al. 2024).

2.2.2. Rotation

The P98 grid does not account for the effect of rotation on stellar evolution. This is the main reason why we decided to only compute non-rotating MESA models. The second reason is related to the coupling of the stellar rotation to the overall evolution of the system and is discussed in Sect. 2.3.6.

2.2.3. Convection

In our MESA simulations, we adopted convection physics as close as possible to that of the P98 grid. We used the same value of α_MLT = 2.0. We also did not account for semi-convection. In order to reduce the superadiabaticity in radiation-dominated convective regions of the most massive stars, we used MLT++ as in Choi et al. (2016). In the P98 grid, they used a modified criterion for convective stability that is not available in the default MESA version. In their modeling, a layer is convective whenever ∇_rad ≥ ∇_ad + δ, where δ is the product of a free parameter δ_ov, their overshooting constant, and a factor which depends on the ratio of radiation to gas pressure ζ. The effect of this prescription corresponds to a step overshoot with a decreasing value of α_ov along the MS³. Instead of implementing the P98 prescription in MESA, we chose to follow a more standard approach by using a step overshoot prescription for the MESA simulations. We calibrated the value of α_ov as to fit the maximum radius reached during the MS by the 50 M_⊙ SEBA star, which corresponds to the upper mass limit of the P98 grid. As will be discussed in Sect. 3.2, the reaction to mass-loss of the H00 tracks is not self-consistent and alters the maximum radius reached during the MS. In order to free ourselves from this issue, we considered necessary to perform the α_ov calibration without accounting for any mass-loss, as the tracks are reliable in this case (at least for the MS evolution and for a mass inside the range of the P98 grid). As to illustrate our point, we show in Sect. 3.2 a comparison of the maximum MS radius reached with MESA and SEBA for different initial masses and mass-loss scaling factors.

Following this approach, we obtained the value α_ov = 0.31, which we adopted for all masses. We note that this value is similar to that used in the Brott et al. (2011) models (they used α_ov = 0.335). As P98 show in their Fig. 1, their prescription corresponds to an increasing Zero Age Main Sequence (ZAMS) value of α_ov with the mass. This calibrated value may therefore be overestimated for the lower mass stars and underestimated for the higher mass stars. However, in the mass range we considered (M ∈ [8, 120] M_⊙), the corresponding ZAMS value of α_ov stemming from P98 convection prescription only moderately increases with the mass, ranging from 0.3 to 0.42. The convection criteria in P98 implies a decreasing value of α_ov along the MS. In order to mimic this treatment, we computed a MESA toy model starting with a value of α_ov higher than the calibrated one and manually reducing it along the MS. We did not find significant differences between the MS evolution of this model and the one with a constant value of α_ov = 0.31 in terms of radial expansion (see Appendix A for more details). Based on these considerations, we judge that the difference in the treatment of overshoot between the two codes does not affect significantly our results.

2.2.4. Stellar winds

In our MESA simulations, we used the Dutch wind combination as in Glebbeek et al. (2009), Choi et al. (2016). It consists of the Vink et al. (2001) prescription for the temperature range 10 000 K ≤ T_eff ≤ 50 000 K, the de Jager et al. (1988) for T_eff ≤ 10 000 K and the Nugis & Lamers (2000) Wolf-Rayet (WR) prescription for stars with a hydrogen surface mass fraction X < 0.4 and T_eff ≥ 10 000 K. Observations by Bouret et al. (2005), Fullerton et al. (2006) and simulations by Sundqvist et al. (2010) indicate that the usual prescriptions in general overestimate the mass-loss rates by a factor of 2 or more. For this reason, we lowered the prescriptions by a factor 3 using Dutch_scaling_factor = 0.333 as done in Kummer et al. (2023). We recall that wind mass-loss rates in the range of masses we consider are still subject to large uncertainties (e.g., Vink et al. 2001; Krtička et al. 2021; Brands et al. 2022; Vink 2022; Björklund et al. 2023). For specific systems whose evolution is strongly affected by mass-loss, we investigated the choice of the wind prescription by comparing the obtained evolutionary pathways according to the wind scaling factor and the chosen stellar code.

In the SEBA simulations, we used the same set of prescriptions except for a small difference in the WR winds prescription. The Nugis & Lamers (2000) prescription in principle requires a value for the surface hydrogen mass fraction, and is metallicity dependent. More precisely, their Eq. (22) depends on the luminosity L, helium mass fraction Y and metallicity Z as:

$$\begin{array}{c}\hfill \dot{M}={10}^{-11}{\left(\frac{L}{{\text{L}}_{\odot }}\right)}^{1.29}{Y}^{1.7}{Z}^{0.5},\end{array}$$(1)

and is applied when X < 0.4. Rapid codes such as SEBA do not follow in detail the evolution of the abundance profile of the star, and therefore do not provide the values of X,  Y,  Z, which are necessary to use this prescription. In order to get around this issue, we decided to use Eq. (25) in Nugis & Lamers (2000) for the SEBA simulations. It reads:

$$\begin{array}{c}\hfill \dot{M}={10}^{-5.73}{\left(\frac{M}{{\text{M}}_{\odot }}\right)}^{0.88},\end{array}$$(2)

and depends only on the mass, which makes it usable for the SEBA simulations. Because the WR mass-loss rates are in principle metallicity dependent, Nugis & Lamers (2000) recommend to use Eq. (1), but given the nature of SEBA, we consider that using Eq. (2) provides an acceptable alternative.

The WR winds treatment between MESA and SEBA thus differ in two aspects. The first is that two different formulas are used. The second is in what activates the wind prescription. In the MESA simulations, the WR winds are activated when the surface hydrogen mass fraction becomes smaller than 0.4. In contrast, in the SEBA simulations, they occur when the star becomes a helium star, i.e., when it has lost its whole envelope⁴. As will be discussed in Sect. 3.2, these two different approaches can in some cases lead to significantly different evolution in the Hertzsprung-Russel diagram (HRD).

In Sect. 3.1, we compare the MESA and SEBA single star MS evolution with our default wind prescription (Dutch combination lowered by a factor 3). We then investigate the impacts of mass-loss by comparing the evolution with stronger winds (standard Dutch combination) in Sect. 3.2.

2.3. Triple systems physics

The physics of the triple system adopted in this study is identical to that in Toonen et al. (2016), Kummer et al. (2023) (apart from the destabilization criterion, see below). We recall here the relevant components.

2.3.1. Stability of the triple system

To ensure stability over a long timescale, the system needs to be hierarchical, in the sense that the outer separation a_out needs to be somewhat larger than the inner separation a_in. In such configuration, the timescale of the perturbation implied by the tertiary component is large compared to the dynamical timescale of the inner binary. When the orbit of the tertiary becomes too close to that of the inner binary, this is not the case anymore, and the secular approximation breaks down. A proper three-body dynamics integration is required to follow the chaotic behavior of the orbits. As shown in (e.g., Hamers et al. 2022; Toonen et al. 2022; Bruenech et al. 2025), the system may in this case dissolve into either a single star and a binary or three single stars, be subject to a collision, or preserve its hierarchy. By default, in TRES the stability criterion is that of Mardling & Aarseth (1999, 2001). In this study, we used the more recent criterion of Vynatheya et al. (2022), based on N-body simulations. In their paper, they provide a publicly available python routine that follows a machine learning approach, which performs better than the criterion by Mardling & Aarseth (2001) for determining the stability of triple systems. They also provide a fitting formula to their result, which reads:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}|}_{\mathrm{crit}}=& 2.4\frac{1+{\stackrel{\sim }{e}}_{\mathrm{in}}}{1-e_{\mathrm{out}}}{\left(\frac{1+q_{\mathrm{out}}}{(1+{\stackrel{\sim }{e}}_{\mathrm{in}})\sqrt{1-e_{\mathrm{out}}}}\right)}^{2/5}\hfill \\ \hfill & \times [\frac{1-0.2{\stackrel{\sim }{e}}_{\mathrm{in}}+e_{\mathrm{out}}}{8}(cos{i}_{\mathrm{mut}}-1)+1],\hfill \end{array}\end{array}$$(3)

where e_out is the outer eccentricity; i_mut the mutual inclination; and $q_{\mathrm{out}}=\frac{M_3}{M_1+M_2}$ where M₁,  M₂, and M₃ are the masses of the components. The term ${\stackrel{\sim }{e}}_{\mathrm{in}}$ is defined in Eq. (3) of Vynatheya et al. (2022). In our simulations we used Eq. (3) as the definition of the critical ratio of $\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}$. Systems are stable when

$$\begin{array}{c}\hfill \frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}>{\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}|}_{\mathrm{crit}}.\end{array}$$(4)

Equation (4) defines a parameter space of stable orbits. In all the presented simulations, the systems are born stable. Some of them may destabilize at some point of their evolution. This is typically the case when the inner binary is subject to significant mass-loss, which makes a_in increase and q_out decrease. The modeling of triple systems in the dynamical unstable regime with a detailed stellar code was beyond the scope of this project. We therefore stop our simulations whenever criterion (4) is violated.

2.3.2. Three-body dynamics and ZLK oscillations

The lowest order (i.e., quadrupole order) perturbation of the third component to the orbital evolution of the inner binary is the ZLK mechanism. It consists of an exchange of angular momentum between the inner and the outer orbit. In some cases, it can induce significant oscillations of the eccentricity of the inner orbit and the mutual inclination.

In the approximation of the test-particle, and for an initially circular orbit, the maximum inner eccentricity reached through ZLK oscillation only depends on the initial mutual inclination:

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

In this approximation, ZLK cycles only occur when the initial inclination is in the range i_mut ∈ [39.2^A°, 140.8^A°] (Ford et al. 2000; Naoz 2016).

For systems with larger values of a_in/a_out, i.e., with a less pronounced hierarchy, additional terms in the secular approximations are needed. The next (octupole) term in the secular approximation can lead to extreme eccentricity oscillations for systems with a non-zero outer eccentricity and an inner binary with unequal stellar masses. The importance of the octupole term with respect to the quadrupole term can be measured with the octupole parameter, defined by:

$$\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}$$(6)

Systems with an octupole parameter |ϵ_oct|≳0.001 are said to be in the eccentric ZLK regime, where the octupole term can play a significant role (e.g., Lithwick & Naoz 2011; Shappee & Thompson 2013).

2.3.3. Precession

Binary and triple systems can be subject to precession of the orbit. In the case of binary systems, precession can affect the orientation of the orbit, but not the overall evolution of the system.

In triple systems, the picture is more complicated as ZLK cycles also induce a precession of the orbit. Other causes of precession (distortion of the star due to its intrinsic rotation or tides, general relativistic effects) can enter in competition with the precession driven by the ZLK oscillations and quench it, which can in turn alter the evolution of other orbital parameters (e.g., eccentricity, inclination) (see Naoz 2016 for a comprehensive review of the ZLK mechanism). The different sources of precession included in TRES are listed in Toonen et al. (2016). Among them, we recall the equation driving the precession caused by the tidal distortion of the star (Smeyers & Willems 2001):

$$\begin{array}{c}\hfill {\dot{g}}_{\mathrm{tides}}=15\frac{k_{\mathrm{AMC}}}{{(1-e^2)}^5}{\mathrm{\Omega }}_{\mathrm{orb}}(1+\frac{3}{2}{e}^2+\frac{1}{8}{e}^4)\frac{M_2}{M_1}{\left(\frac{R}{a}\right)}^5,\end{array}$$(7)

and by intrinsic stellar rotation (Fabrycky & Tremaine 2007):

$$\begin{array}{c}\hfill {\dot{g}}_{\mathrm{rotate}}=\frac{k_{\mathrm{AMC}}}{{(1-e^2)}^2}\frac{{\mathrm{\Omega }}_{\mathrm{spin}}^2}{{\mathrm{\Omega }}_{\mathrm{orb}}}\frac{M_1+M_2}{M_1}{\left(\frac{R}{a}\right)}^5.\end{array}$$(8)

In Eqs. (7) and (8), Ω_orb stands for the orbital angular velocity of the binary system, Ω_spin the angular velocity of the star, R its radius. These two equations depend on the apsidal motion constant of the star, which can be directly determined from the structure of the stars in the MESA simulations, offering a more self-consistent modeling of the precession of the triple system, and hence of its overall orbital evolution (see Sect. 2.3.6).

2.3.4. Tidal interactions

Unless stated otherwise, the chosen tides formalism in the presented simulations is that of Hurley et al. (2002, hereafter H02). We recall that under this formalism, both the equilibrium and dynamical tides are incorporated into the formalism of Hut (1981).

In Sciarini et al. (2024, hereafter S24), we demonstrated that following this procedure, the tidal evolution obtained in the case of the dynamical tides is not consistent with the original formalism by Zahn (1977, hereafter Z77). In this letter, we obtained an updated set of equations consistent with the original formulation by Z77. Yet, for this study we chose in most cases to stick to the H02 formalism. There are mainly two reasons for this choice.

The first is that the approach followed in this study is to primarily investigate the effect of the detailed stellar physics onto the evolution of the triple system, keeping the interactions unchanged. The present study complements that of Kummer et al. (2023), who investigated the evolutionary pathways of massive hierarchical triple systems. We decided to keep the same physics for the interactions as in their study.

The second is that the derivation proposed by Z77, which we complemented in S24, is an expansion keeping only low-order eccentricity terms. It is well suited for systems with moderate eccentricity (say e ≲ 0.5), but likely to be less appropriate for systems with large eccentricities, which are often found among triple and higher order systems due to the ZLK oscillations. The equilibrium tides formalism of Hut (1981) assumes that tides make the star reach an equilibrium shape with a constant time lag between the tidal bulge and the line of centers of the two bodies, an assumption not well justified in the case of the dynamical tides where the star oscillates. Its strength lies in its applicability to systems with any orbital eccentricity. Under the formalism of H02, the functional dependence of the eccentricity and the angular velocity in the case of the dynamical tides is exactly the same as in the case of the equilibrium tides formalism of Hut (1981). It might be more suitable for systems with very high eccentricities, even though its functional dependence (in particular the angular velocity dependence) intrinsically comes from the formalism of Hut (1981), i.e., a formalism of equilibrium tides.

In Sect. 4.7, we compare the evolution obtained with the H02 formalism to that obtained with the S24 prescription (updated Z77) for a couple systems subject to ZLK oscillations and tidal interactions.

2.3.5. Mass transfer

Because triple systems are subject to ZLK oscillations, mass transfer frequently occurs in eccentric orbits (Toonen et al. 2020). Mass transfer in eccentric orbits is still poorly understood (see, e.g., Sepinsky et al. 2007; Davis et al. 2013; Hamers & Dosopoulou 2019). Whenever a star fills its Roche lobe, the simulation is stopped and mass transfer (either in the inner binary or from the tertiary star) is considered as the evolutionary pathway of the system. It is therefore appropriate to use the value of the Roche lobe radius at periastron to determine if a star fills its Roche lobe. In this case, the expression of Eggleton (1983) needs to be corrected with an eccentricity dependent term

$$\begin{array}{c}\hfill R_{\mathrm{L},i}=a_j\frac{0.49q_i^{2/3}}{0.6q_i^{2/3}+ln(1+q_i^{1/3})}(1-e_j),\end{array}$$(9)

where $q_1=\frac{M_1}{M_2}$, $q_2=\frac{M_2}{M_1}$, $q_3=\frac{M_3}{M_1+M_2}$, j refers to the inner (i = 1, 2) or outer orbit (i = 3). As such, for the inner binary Eq. (9) neglects distortions to the Roche geometry caused by the third star.

2.3.6. Orbital evolution and coupling to stellar evolution

The set of equations governing the orbital evolution of the triple system is given in Toonen et al. (2016). In TRES, it is solved within the SECULATRIPLE code. Notably, the evolution of the stars’ angular velocity is included in the set. The rotation of the star is thus treated outside the stellar evolution code. The main reason is that by default TRES uses SEBA as stellar evolution code, which does not account for the effect of rotation on stellar evolution.

During a simulation with TRES, the stellar and secular evolution codes are called successively until the next computed age (next_age). In order to reach the next_age, different timesteps are generally used by the stellar and the secular evolution codes. The main reason is that depending on the orbital configurations, the secular timescale can differ significantly from the timescale of stellar evolution. The angular velocities being included inside the set of equations governing the orbital evolution, their evolution is solved for together with the quantities whose evolution is altered by the dynamics of the system. For this reason, the consistent coupling of the rotation of the star to the overall evolution of the system is not straightforward. In order to account for the effect of rotation onto stellar evolution, stellar and orbital evolution would need to share the same timesteps, which might be very constraining for the stellar evolution when the ZLK timescales are short with respect to the stellar evolution timescales.

For these reasons, we followed a simpler approach by neglecting the effects of stellar rotation on stellar evolution. This is the second reason why there are no feedback from the interactions to the stellar evolution in our simulations.

2.3.7. Evolutionary pathways

The possible evolutionary pathways of the triple systems in this study are the same as in Kummer et al. (2023). They consist of stopping conditions of the presented simulations:

1. Mass transfer in the inner binary when the primary or the secondary fills its Roche lobe.

2. Mass transfer from the tertiary to the inner binary when the tertiary fills its Roche lobe.

3. Dynamical destabilization when criterion (4) is violated.

4. Dynamical unbinding of one of the orbit as a result of core collapse supernova of one of the three components.

5. None of the above interactions occur within a Hubble time. In this case the systems are considered as non-interacting.

3. Single star evolution

One of the main driver of the evolution of massive hierarchical triple systems is the stellar evolution of its components. The final state of the system depends on the interactions the stars experience during their lives, which are strongly impacted by their radial expansions and mass lost by stellar winds. In our simulations we do not account for feedback from the interactions onto stellar evolution. Thus, it is sufficient to get a satisfactory understanding of the differences in the single star evolution provided by MESA and SEBA to understand how they can imply divergences in the predicted evolutionary pathways of triple systems. In this section we aim at comparing and explaining the obtained differences in the single star evolution provided by MESA and SEBA under the choice of stellar physics ingredients presented in Sect. 2.2.

3.1. Grid comparison

The HRD of massive stars with initial masses in the range M ∈ [8,  120] M_⊙ obtained with MESA and SEBA and the stellar physics ingredients described in Sect. 2.2 is shown in Fig. 1a. Comparing MESA and SEBA tracks, we observe that:

Fig. 1. HRD obtained with MESA (solid lines) and SEBA (dashed lines). Isoradius lines are represented by gray dashes. (a) HRD for stars in the mass range M ∈ [8,  120] M⊙. (b) HRD for stars in the mass range M ∈ [50,  120] M⊙.

1. The chosen calibration of α_ov offers a satisfactory match between the MS width of MESA and SEBA models among the less massive stars. The location of the MS hook for stars in the mass range M ∈ [8,  50] M_⊙ (i.e., for stars inside the range of the P98 grid) is relatively consistent.

2. The late evolution of MESA and SEBA models is relatively consistent in the mass range M ∈ [8, 50] M_⊙, although the maximum radius reached during the evolution is slightly higher with SEBA.

3. MESA and SEBA tracks exhibit increasing discrepancies with mass, especially for M ≥ 50 M_⊙, i.e., outside the range of the P98 grid. In particular, the maximum radius differs by about an order of magnitude in this mass range.

In order to get a closer look at the most massive stars (outside the range of the P98 grid), we show in Fig. 1b a denser HRD for stars in the range M ∈ [50,  120] M_⊙, which allows to observe the following differences:

1. The MESA MS width in this range of mass and with this choice of wind scaling factor is larger than that of SEBA. MESA stars expand more during their MS, therefore the MESA MS width is larger among massive stars. The simplest explanation for this difference is that stars of these masses are outside the range of the P98 grid, and the fitting formulae proposed by H00 do not generalize well for stars outside the range of the grid. There is a second explanation for these discrepancies: at high masses, stars are subject to significant stellar winds (even when they are moderated by a scaling factor of 0.333), and the H00 tracks do not react consistently to mass-loss. This point will be discussed in more details in Sect. 3.2.

2. The maximum radius reached during the evolution differ by one order of magnitude or more. No matter the initial mass, SEBA tracks reach the red supergiant (RSG) phase, and extend up to several thousands solar radii. In contrast, the most massive MESA models stop at some point during their crossing of the Hertzsprung gap and start going back to the blue. As discussed in Sect. 3.2, this is a consequence of the more consistent reaction to mass-loss of the MESA models. Because of this, the maximum extent of the most massive MESA stars is significantly smaller than that of the SEBA models. This has important consequences for the evolution of triple systems, as the occurrence of some interactions (in particular mass transfer) strongly depend on the radial expansion of the stars.

3. SEBA models with M ≥ 70 M_⊙ end their evolution as WR stars. The corresponding MESA models do not reach the WR phase and end their evolution as blue supergiants.

4. Even with a moderate wind scaling factor, the inconsistent reaction to mass-loss of the H00 tracks can be noticed among the more massive stars. A distinct inflection point is observed during the MS evolution of the SEBA models. They are marked with a dot on each track. After this inflection point, the luminosity of the stars starts decreasing. The reason of the inflection points and the subsequent decrease in luminosity is a direct consequence of H00 tracks reaction to mass-loss. They correspond to the moment the models reach the Vink bistability jump, which makes the mass-loss increase. The greater mass-loss changes the terminal age main sequence (TAMS) luminosity of the star, which is a simple function of the current mass in the H00 tracks⁵. As a result, the luminosity of the star decreases as its mass decreases, and the track contains a direct imprint of the bistability jump. Although MESA models’ evolution also react to mass-loss, they do not show a decrease in luminosity when reaching the bistability jump, which arises in the SEBA models when the star shifts towards the track of a star of lower mass.

3.2. Impacts of mass-loss

In this section, we focus on stars in the upper mass range (M ∈ [50,  120] M_⊙) and on the models’ respective reactions to mass-loss. The HRD for stars in the same mass range as in Fig. 1b but with Dutch_scaling_factor=1 is shown in Fig. 2.

Fig. 2. Same as Fig. 1b, but with a wind scaling factor three times higher.

As expected in this high-mass range where stellar winds play an important role, increasing the wind scaling factor strongly affects the evolution of the stars. However, the reaction to mass-loss is very different between the MESA and the SEBA tracks.

As with the lower wind scaling factor, the maximum radial extent significantly differs between MESA and SEBA tracks, but this difference is even more pronounced. Even with this high mass-loss rate, SEBA models reach the RSG phase and expand up to several thousands solar radii. In contrast, the maximal radial extent of the MESA models is strongly reduced in this case, as a result of the significant mass-loss.

With the high wind prescription, MESA and SEBA models all end their evolution as WR in this mass range, but the evolution prior to the WR phase differs. Due to the strength of the winds, the more massive MESA stars already reach the WR phase during their MS and never expand more than a few tens of solar radii (in the presented models this is the case for the M = 100 M_⊙ and the M = 120 M_⊙ stars). This can be explained as follows: when the winds are strong enough, the mass-loss timescale becomes smaller than the nuclear timescale. In this case, winds are able to strip the whole hydrogen envelope, which induces the WR phase before the end of the MS. Once the helium enriched layers are revealed, the surface opacity changes which makes the star contract (see Bavera et al. 2023 for a similar discussion). Among the less massive models, the winds are weaker and the stars reach the WR phase later. After the end of their MS, they start crossing the Hertzsprung gap, but come back to the blue and do not become RSGs (for the 85 M_⊙ this happens shortly after the MS and the star’s Hertzsprung gap crossing is minor). Their maximal radial extent is lower than with moderate winds.

Unlike MESA models, SEBA models all experience a RSG phase before becoming WR stars. This is because the condition for SEBA models to become WR (helium stars) is to lose their whole envelope, i.e., to reach M_env = M_tot − M_core = 0. This never happens during the MS as the core mass is set to 0 during this stage. At the end of the MS, a core mass is defined, which does not account for the mass already lost during the MS (see Sect. 5.1.2 in H00). As a result, the envelope mass is overestimated for stars that have undergone significant mass-loss due to stellar winds, which delays their transition to the WR phase. The inflection point during the MS evolution of the most massive SEBA models, which is caused by the Vink bistability jump and the models’ reaction to mass-loss is even more pronounced in this case, as the mass-loss is larger.

The maximum MS radius and the maximum radius play an important role in determining the evolutionary pathways of triple systems, as will be discussed in Sect. 4. The interactions often depend on whether the stellar components fill their Roche lobe at some point of the evolution, which is mainly determined by the maximum radius they reach. We therefore show in Fig. 3 the maximum radius reached by MESA and SEBA models during the MS and the whole evolution with different values for the wind scaling factor (for this comparison we also computed MS models without mass-loss; we therefore have three values for the scaling factor for the maximum MS radius: 0, 0.333, and 1).

Fig. 3. Upper panel: Maximum radius reached during the MS by the MESA and SEBA models for different values of the wind scaling factors (0, 0.333, and 1). Lower panel: Same as upper panel, but for the maximum radius reached during the whole evolution and for wind scaling factors of 0.333 and 1.

As can be observed in the upper panel of Fig. 3, the maximum MS radius R_max, MS of the MESA and SEBA models show a very different mass and wind scaling factor dependence. The maximum MS radius of the SEBA models show a monotonic dependence with the mass, which can be easily understood considering that following the H00 formulae, the TAMS radius is a simple increasing function of the mass⁶. More interestingly, one notices that for any initial mass, the maximum MS radius of the SEBA models decreases when more mass-loss is accounted for. This can be simply understood noting that in a SEBA run accounting for mass-loss, the value of the TAMS radius is updated with the current value of the mass, as is done for the TAMS luminosity. After a star has lost some amount of mass, it shifts towards the track of a lower mass star, and therefore ends up with a smaller TAMS radius, which explains the mass-loss dependence of R_max, MS for the SEBA models.

Detailed calculations with MESA show that the picture is more complicated. Although models without mass-loss show a mass dependence of R_max, MS consistent with SEBA (i.e., increasing with the mass), the monotonicity of R_max, MS with the initial mass breaks down when mass-loss is accounted for. With a moderate mass-loss (Dutch/3), R_max, MS peaks around 70 M_⊙, while with a higher mass-loss (Dutch), the peak is around 50 M_⊙. At high masses, the mass-loss is so strong that it prevents the stars from expanding.

It is interesting to have a look at the dependence of R_max, MS on the wind mass-loss for a fixed mass. At 50 M_⊙, the SEBA and the MESA models show an opposite behavior. When no winds are accounted for, the MESA and SEBA models have the same R_max, MS (which is due to the of the calibration we performed, see Appendix A). With greater mass-loss, R_max, MS of the MESA model increases, whereas R_max, MS of the SEBA model decreases. At a higher initial mass (e.g. 70 M_⊙), the mass-loss dependence of R_max, MS of the MESA model is non-monotonic: the model with moderate winds has a higher maximum MS radius than the model without mass-loss, whereas the model with high winds expands less than the model without winds. It is worth noting the convergence of the SEBA and the MESA models for masses smaller than 50 M_⊙, i.e., inside the range of the P98 grid. The excellent match of R_max, MS of all MESA and SEBA models with moderate or no mass-loss in the range M ∈ [8, 50] M_⊙ validates the performed calibration of α_ov.

The maximum radius R_max reached by the MESA and SEBA models is shown in the lower panel of Fig. 3. For any model in the considered mass range, the maximum radial expansion predicted by MESA and SEBA substantially differ. The discrepancy between the MESA and SEBA models increases with the mass and the wind scaling factor. This directly results from the inconsistent reaction to mass-loss of the H00 tracks. It is more pronounced at high masses as the mass-loss is more efficient in this case. As explained above, SEBA models experience a RSG phase before becoming WR stars, which explains their large maximum radial extent. MESA models in contrast react self-consistently to mass-loss. Their maximal radial extent is sensitive to the initial mass and to the efficiency of the winds. For the considered range of masses and wind scaling factors, the maximum radial extent strongly decreases with increasing mass and wind scaling factor. In particular, when a wind scaling factor of 1 is applied, we observe a discrepancy of 2 orders of magnitude or more of the maximum radius of MESA and SEBA models with M ≥ 85 M_⊙. In the mass range M ∈ [8, 50] M_⊙, the MESA and SEBA models show small discrepancies in R_max, which decrease with the mass. They can be attributed to the winds, which are not completely negligible in this mass range, and increase with the initial mass.

4. Triple evolution

The initial conditions of the different systems presented in this section are given in Table 1. The last row indicates whether the evolutionary pathways predicted by MESA and SEBA differ, following the nomenclature adopted in Kummer et al. (2023)⁷.

4.1. Validation with η Carinae

η Carinae is an observed binary system with two massive stars (M₁ ∼ 90 M_⊙ and M₂ ∼ 30 M_⊙) in a highly eccentric orbit (e = 0.9) and period P = 5.5 yr (Damineli et al. 1997; Davidson 1997). It is famous for undergoing the “Great Eruption” in the 19th century, during which event the luminosity of the system increased by orders of magnitudes, until it became the second brightest star of the sky (de Vaucouleurs & Eggen 1952). The Homunculus Nebula, now surrounding the η Carinae system, is an evidence of this tragic event (Humphreys et al. 1999). Up to this day, there is no clear consensus on the cause of the Great Eruption. Different scenarios have been proposed, which include mass transfer from η Carinae B at periastron passage (Kashi & Soker 2010), Luminous Blue Variables (LBV) outbursts (e.g., Humphreys & Davidson 1994), or the merger of a close inner binary in a hierarchical triple system, leaving behind a binary in a wide orbit (Portegies Zwart & van den Heuvel 2016; Hirai et al. 2021).

Toonen et al. (2016) reproduced the evolution of the favored model of Portegies Zwart & van den Heuvel (2016) with TRES. In their simulation, SEBA was used for the stellar evolution. The system being composed of very massive components (in their simulation the primary is of initial mass M₁ = 110 M_⊙), the stellar evolution provided by MESA is expected to differ substantially from that of SEBA (see discussion in Sect. 3.1).

We replicated the system using the same initial conditions as in Toonen et al. (2016), i.e., the favored model of Portegies Zwart & van den Heuvel (2016). The system is labeled “η Carinae” in Table 1. In Fig. 4, the radius expansion of the primary and secondary are shown, together with the Roche lobe radius of the primary. The evolution of the inner eccentricity and of the apsidal motion constant are also presented. The evolution provided by SEBA is compared to that of MESA. We find that for this specific system, the evolutionary pathway of the system obtained with MESA is consistent with SEBA, i.e., the primary fills its Roche lobe, which is expected to lead to a common envelope event and eventually a merger of the inner binary, given the high mass ratio of the binary (Toonen et al. 2016).

Fig. 4. Evolution of the η Carinae system simulated with MESA and SEBA. Upper panel: Inner eccentricity. Middle panel: Radius expansion of the primary (dashed line) and secondary (dash-dotted line). The Roche lobe radius of the primary is shown as a solid line. Lower panel: Evolution of the apsidal motion constant of the primary.

The system begins in a phase of eccentric ZLK oscillations (initially the octupole parameter is ϵ_oct = 0.0048), which are quenched after ∼1.3 Myr in the SEBA simulation, ∼1.5 Myr in the MESA simulation. The quenching of the ZLK oscillations is due to the precession caused by to the distortion and rotation of the stars, which overtakes the precession caused by the ZLK mechanism (Toonen et al. 2016). The visible oscillations of the Roche lobe radius (see Eq. (9)) are a consequence of the ZLK cycles: the Roche lobe radius is taken at periastron and therefore directly depends on the eccentricity. Subsequently, tides circularize the system (around t = 1.6 Myr) and the primary fills its Roche lobe in a circular orbit.

The main differences between the outcome of the MESA and SEBA simulations are found in the timescale of the evolution and the radius expansion of the primary. Because the MESA model expands more during its MS (Sect. 3.2), the star fills its Roche lobe after only 2.2 Myr (2.7 Myr in the SEBA run). We further observe that the amplitude of the ZLK oscillations is larger in the MESA simulation. This difference can be attributed to the fact that in the MESA run, the apsidal motion constant is computed consistently with the stellar structure, whereas it is assigned a fixed value in the SEBA run. As shown in the lower panel of Fig. 4, the value at ZAMS of the apsidal motion of the MESA model is relatively consistent with the value of the SEBA model (k_AMC = 0.0144), but at the end of the simulation, it has decreased by more than an order of magnitude. The decrease of k_AMC reduces the precession caused by the distortion of the star (Eqs. (7) and (8)), allowing the precession due to the ZLK mechanism to play a more prominent role, which enables the ZLK oscillations to perdure longer. We note that this aspect per se may imply divergences in the evolutionary pathways of triple systems computed with MESA and SEBA. In triple systems, the ZLK oscillations typically increase the rate of interactions between the components of the inner binary compared to such interactions in isolated binaries (e.g., Antonini et al. 2017; Hamers & Dosopoulou 2019; Toonen et al. 2020; Stegmann et al. 2022; Kummer et al. 2023). By obtaining the apsidal motion constant from the structure, the relative effect of the precession due to ZLK cycles is enhanced, which should increase the number of systems that interact as a result of the ZLK oscillations. In Sect. 4.6 and in Appendix B, we present examples of triple systems representative of this aspect.

The differences in the stellar evolution provided by MESA and SEBA are not sufficient to imply a divergent evolutionary pathway for this system, nevertheless it consists of an interesting science case and a proof of concept of our coupling of MESA to TRES.

4.2. Divergent evolution: Inner versus tertiary mass transfer

In binary systems, because the primary star is more massive and evolves faster, mass transfer in general occurs from the primary to the secondary (although stellar winds can complicate the picture among very massive stars; see e.g., Bavera et al. 2023). In the case of triple systems, stable configurations can exist where the tertiary is more massive than each component of the inner binary. One observed system with this property is TIC 470710327 (Eisner et al. 2022). In such systems, mass transfer can theoretically be initiated from the tertiary to the two stars of the inner binary.

The system presented in this section is labeled “Tertiary MT” in Table 1. It is composed of three massive stars, with a tertiary star more massive (M₃ = 80 M_⊙) than the two stars of the inner binary (M₁ = M₂ = 70 M_⊙). This system was chosen to represent a class of systems for which a divergent evolutionary pathway is obtained whether MESA or SEBA is used as stellar code and is not unique. We furthermore chose a system not subject to ZLK oscillations (i_mut = e_out = 0) to simplify the description.

Given the initial masses of the system, it is not clear whether an inner or tertiary mass transfer will occur. The tertiary being more massive, it evolves faster and has a shorter MS lifetime. It is expected to expand significantly over its MS (and also after the MS in the case of the SEBA model, as shown in Fig. 1b), which may trigger a mass transfer episode from the tertiary to the inner binary. But since the components of the inner binary are closer, they also possess a smaller Roche lobe radius, and might therefore interact if the tertiary does not expand enough to fill its Roche lobe.

We simulated this system with MESA and SEBA and compared the obtained evolutionary pathways. The radius of the primary (tertiary) and its Roche lobe radius evolution are shown in the upper (lower) panel of Fig. 5. In this case, the evolution obtained with the two codes differ significantly. In the SEBA simulation, the MS expansion of the components of the inner binary is not sufficient for them to fill their Roche lobe. As a result, the tertiary completes its MS while the primary and secondary are still far from interacting. It subsequently expands to several thousands solar radii (the maximum radius of the 80 M_⊙ SEBA model is R_max, SeBa ∼ 3500 R_⊙), becoming a RSG and fills its Roche lobe (R_L ∼ 2000 R_⊙), which triggers a case B MT to the inner binary.

Fig. 5. Evolution of the system labeled “Tertiary MT” in Table 1. Upper panel: Radius expansion of the primary (dashed line) and Roche lobe radius (solid line). Lower panel: Same as upper panel, but for the tertiary.

The system simulated with MESA shows a different behavior. After the tertiary completes its MS, it keeps expanding for a while, before starting to contract. It reaches a maximum radius of only R_max, MESA ∼ 650 R_⊙, which is insufficient for it to fill its Roche lobe. Thus, the two stars of the inner binary end up filling their Roche lobe, which leads to a contact phase.

4.3. Divergent evolution: Mass transfer versus destabilization

The system presented in this section is labeled “Dyn. Destab” in Table 1. It is composed of three massive stars, with the two components of the inner binary significantly more massive (M₁ = 100 M_⊙, M₂ = 90 M_⊙) than the tertiary (M₃ = 25 M_⊙). Initially, it is stable as it has $\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}=10$ and ${\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}|}_{\mathrm{crit}}\sim 5.8$. This system represents a class of systems for which divergent evolutionary pathways (inner mass transfer or dynamical destabilization) can be obtained depending whether MESA or SEBA is used as stellar code.

Being close to dynamical destabilization, the evolution of the system is sensitive to the amount of mass lost through stellar winds. Stellar winds in this mass range being furthermore very uncertain, we simulate the evolution of the system with two winds prescription, the standard Dutch prescription and the Dutch prescription lowered by a factor 3. The use of various winds prescription allows a more robust comparison between the outcome of simulations with MESA and SEBA.

The evolution of the relevant quantities $(\mathrm{mass},\mathrm{radius},\mathrm{Roche}\mathrm{lobe}\mathrm{radius}\mathrm{of}\mathrm{the}\mathrm{primary},\mathrm{ratios}\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}$ and ${\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}|}_{\mathrm{crit}})$ of the system obtained with MESA and SEBA with the two winds prescriptions are shown in Fig. 6. We note that the obtained evolution of the triple system in the SEBA simulations are not strongly altered by the wind prescription, as can be expected considering the discussion in Sect. 3.2. In both cases, the primary fills its Roche lobe after the MS while crossing the Hertzsprung gap (case B MT).

Fig. 6. Evolution of the system labeled “Dyn. Destab” in Table 1. Upper panel: Mass evolution of the primary. Middle panel: Radius and Roche lobe radius evolution of the primary. Lower panel: Evolution of the ratios $\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}$ and ${\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}|}_{\mathrm{crit}}$.

The MESA simulations are more sensitive to the winds prescription. This can be simply understood looking at the stellar tracks in Sect. 3.2. With the standard Dutch prescription, the winds are strong enough for the 100 M_⊙ MESA star to become a WR already during the MS. Its radius starts shrinking around t ∼ 2.7 Myr and never extends beyond ∼20 R_⊙. The star eventually experiences WR winds and the system reaches the dynamical destabilization. With the lowered Dutch winds prescription, the star does not lose enough mass to reach the WR phase during the MS. In this case, the radius expansion is sufficient for the star to fill its Roche lobe during the MS, which leads to case A MT.

The system presented in this section is interesting in several aspects. Firstly, it illustrates how the uncertainties in the winds prescriptions can imply uncertainties in the evolution of triple systems. Secondly, it illustrates how an inconsistent reaction to mass-loss can impact the predicted evolutionary pathway of a triple system, leading it in one case to a mass transfer, in the other case to a dynamical destabilization. In Bavera et al. (2023), the same single star physics aspects were discussed, together with the consequences they imply for the evolution of binary systems. They showed that above a certain mass and with the Dutch winds prescriptions, stars computed with MESA reach the WR phase early in the evolution and as a result do not expand more than a few tens of solar radii. Similar stars computed with H00 tracks based codes (such as SEBA) become RSGs and expand by several thousands of solar radii. They showed the effect of these important differences on the evolution of binary systems. The natural consequence of this discrepancy is that systems simulated with H00 tracks based codes interact more often at high masses given that the stars expand more.

4.4. Divergent evolution: Mass transfer versus unbinding of the orbit

The initial conditions of the system presented in this section are quite close to those of the system of the previous section, except that the chosen masses are smaller and equal (M₁ = M₂ = M₃ = 85 M_⊙) and the tertiary is situated at a larger distance. It is labeled “Orbit unbound” in Table 1. The initial masses were chosen to be identical for simplicity, but a wide range of initial masses would lead to equivalent results given the largely discrepant values of R_max obtained between the two stellar codes (see Fig. 3). The evolution is shown in Fig. 7, similarly to Fig. 6 with in addition the Roche lobe radius of the tertiary.

The main difference between this system and that of Sect. 4.3 is that the present one is initiated far from dynamical destabilization. In this case, stellar winds are not sufficient to destabilize the system. The systems simulated with SEBA or with MESA and moderate winds interact as the primary expands enough to fills its Roche lobe. In contrast, the components of the system simulated with MESA and high winds do not expand beyond ∼35 R_⊙. They do not fill their Roche lobe and thus do not interact before the supernova, similarly to what is found by Bavera et al. (2023) for binary systems. In this case, because the tertiary is on a wide and highly eccentric orbit, it gets unbound from the inner binary as it undergoes the supernova explosion⁸.

4.5. Divergent evolution: Mass transfer versus non-interacting system

The tertiary of the system presented in Sect. 4.4 is in a wide outer orbit, which makes it likely to get unbound by the supernova explosion. In the study by Kummer et al. (2023), only a very small fraction of the simulated systems (0.1–0.4%) do not experience any of the interactions accounted for in TRES (see Sect. 2.3.7) and end up as triple compact objects systems (TCOs). The reason this category is only scarcely populated is due to a combination of two competing factors. On the one hand, the inner orbit needs to be relatively wide to avoid RLOF. They find that systems with P_in ≲ 6 × 10³ days cannot avoid mass transfer. This is a direct consequence of the maximal radial expansion of the stellar models (see Fig. 3). To ensure the stability of the system, this implies that the outer period also needs to be wide, which makes it likely for the tertiary to become unbound due to the SN kick.

Non-interacting triple systems are of particular interest in the context of gravitational wave (GW) events as TCO systems are considered to be a main formation channel of mergers (e.g., Silsbee & Tremaine 2017; Antonini et al. 2017; Rodriguez & Antonini 2018; Fragione & Antonini 2019; Martinez et al. 2022; Trani et al. 2022; Dorozsmai et al. 2024; Su et al. 2024). In this scenario, ZLK oscillations caused by the presence of a third object can significantly increase the eccentricity of the inner double compact object (DCO), which can lead it to merge through GW emission.

According to the discussion in Sect. 3.2, the maximum radial expansion of the MESA models at solar metallicity and with the standard Dutch combination is smaller by one to two orders of magnitudes than that predicted by SEBA in the mass range M ∈ [50, 120] M_⊙. Given Kepler 3rd law (P² ∝ a³) and the proportionality between the semi-major axis and the Roche lobe radius in Eq. (9), we can estimate the limiting period for avoiding RLOF in MESA simulations to be up to three orders of magnitude lower than what is predicted by SEBA (i.e., P_min, MESA ∼ 6 days whereas P_min, SeBa ∼ 6 × 10³ days), substantially increasing the size of the parameter space where TCOs can form at solar metallicity. Such compact systems are also more likely to stay bound after the SN explosion, as low inner separations allow for low outer separations, i.e., more gravitationally bound systems.

The system presented in this section is representative of this category and is labeled “TCO” in Table 1. It is almost identical to that of Sect. 4.4, except that it is much more compact, with a_in = 1 au and a_out = 10 au, corresponding to orbital periods P_in = 28 days and P_out = 723 days. In this case, the system was only simulated with Dutch_scaling_factor = 1. The evolution of the radii and Roche lobe radii of the components are shown in Fig. 8. Since the stars are of the same mass, the radii expansions are the same as in Fig. 7. Only the Roche lobe radii differ (they are smaller in this case since the system is more compact). Even in this close inner orbit, the stellar components do not fill their Roche lobe in the MESA simulations, as the maximum radius of the M = 85 M_⊙ model is R_max ∼ 35 R_⊙ and the inner Roche lobe radius stays larger than 80 R_⊙ throughout the simulation. Given the compactness of the system, it does not get unbound by the SN and forms a TCO system.

Fig. 7. Same as Fig. 6 for the system “Orbit unbound” in Table 1.

Fig. 8. Same as middle panel of Fig. 7 for the system “TCO” in Table 1.

The presented system was initiated in a circular orbit and was not subject to ZLK oscillations (i_mut = e_out = 0). It is therefore unlikely to merge by GW emission, but serves as an illustrative example. The formation of GW sources through this channel would be interesting to investigate. Close systems simulated with MESA can more easily avoid mass transfer and are also more likely to survive the SN explosion, which should increase the fraction of formed TCOs (and therefore GW progenitors) compared to what is predicted by Antonini et al. (2017), Rodriguez & Antonini (2018), Trani et al. (2022), Kummer et al. (2023). This aspect will be investigated in a follow-up study.

4.6. Divergent evolution: The importance of self-consistently modeling the precession

The system presented in this section is labeled “ZLK & precession” in Table 1. It is composed of a tertiary more massive than either star of the inner binary (M₁ = 25 M_⊙, M₂ = 10 M_⊙ and M₃ = 30 M_⊙). The tertiary thus evolves faster, but is situated too far from the inner binary to fill its Roche lobe. The system is subject to ZLK oscillations. We simulated the system with MESA and SEBA and included a model variation where the stellar evolution is performed by MESA but the apsidal motion constant is fixed to k_AMC = 0.0144.

The evolution of the relevant quantities (radii of the stars, Roche lobe radii of the primary and tertiary, apsidal motion constant) are shown in Fig. 2.3.3. The system simulated with MESA and a self-consistent modeling of the apsidal motion constant is subject to maximal ZLK oscillations, as the decrease of k_AMC over time makes the precession induced by the ZLK mechanism increasingly dominant. Therefore, the inner eccentricity oscillations stay more pronounced with this model, whereas they are quenched more efficiently with the SEBA model and the MESA model variation.

Fig. 9. Evolution of the system labeled “ZLK & precession” in Table 1. Upper panel: Radius expansion of the primary (dashed line) and Roche lobe radius (solid line). Middle panel: Same as upper panel, but for the tertiary. Lower panel: Evolution of the apsidal motion constant of the primary.

As a result, the primary fills its Roche lobe in the case of the MESA simulation, which leads to case A eccentric mass transfer. In the MESA model variation and the SEBA simulation, the quenching of the oscillations prevents the primary from filling its Roche lobe. The tertiary eventually completes its evolution, undergoes a SN explosion and gets unbound from the inner binary. It is worth noting that in this case, the chosen initial masses lie inside the range of the P98 grid, and the stellar evolution predicted by SEBA is very consistent with that of MESA. The divergence in the evolutionary pathway can thus be solely attributed to the apsidal motion modeling.

4.7. ZLK cycles and tidal evolution

In this section we aim to illustrate the impact of the dynamical tides prescription proposed in S24 (adapted from Z77) on the evolution of triple systems and show how it differs from that of H02. For this purpose we present four triple systems where the inner binary is in an eccentric orbit close enough to experience substantial tides. Two systems reach an inner eccentricity e_in, max ≳ 0.8 during the evolution, either due to ZLK oscillations or because it is initially high, whereas the inner eccentricity of the other two verifies e_in, max ≲ 0.6 throughout the whole evolution. In each category, one system is subject to ZLK oscillations, the other is not. The systems are labeled according to the maximum inner eccentricity they reach and whether they experience ZLK oscillations. The initial conditions of these systems are given in Table 2.

All systems are simulated with SEBA. For each set of initial conditions, two independent simulations are performed, one time with the default H02 prescription, the other with the S24 prescription. A detailed description of the latter model is given in Appendix D for the interested reader. The main conclusion is that in highly eccentric orbits (e ≳ 0.7), this tidal model predicts the eccentricity of the system to increase, which goes against the common picture that tides always act towards the circularization of the system⁹.

The evolution of the two systems whose inner eccentricity remains lower than 0.6 throughout the evolution is shown in Fig. 10. They correspond to the two first systems in Table 2. One of them is subject to ZLK oscillations, the other is not. As most of the tidal evolution occurs near the end of the simulation, the largest part of the plot focuses on the last 0.8 Myr of evolution.

Fig. 10. Comparison of the evolution of systems in moderately eccentric orbits as predicted by the H02 and the S24 prescriptions. One system is subject to ZLK oscillations (blue and magenta lines), the other is not (violet and orange lines). The scale of the horizontal axis changes at t = 2.6 Myr and is chosen so as to focus on the part of the evolution where tidal effects are important. Upper panel: inner eccentricity evolution. Lower panel: Radius and Roche lobe radius evolution of the primary.

For these systems, the evolutions predicted by the two prescriptions are relatively similar¹⁰. Both predict that the tides circularize the system before the primary fills its Roche lobe, which implies that the system experiences mass transfer in a circular orbit. This is because in this low eccentricity regime, the S24 prescription also predicts that dynamical tides tend to reduce the eccentricity (for any pseudo-synchronized star).

The evolution of the two systems whose inner eccentricity reach 0.8 or more during the evolution is shown in Fig. 11. They correspond to the last two systems in Table 2. For these systems, the tidal evolution predicted by the S24 prescription differs substantially from that predicted by the H02 prescription. In this regime, S24 predicts that the dynamical tides make the eccentricity increase¹¹. The systems in this case do not circularize, therefore they experience an eccentric mass transfer episode as the primary fills its Roche lobe. In contrast, the systems simulated with the H02 prescription circularize before the primary fills its Roche lobe and experience mass transfer in circular orbit.

Fig. 11. Same as Fig. 10, but for systems with ein, max ≳ 0.8. Additionally the evolution of the inner separation is shown in the lower panel.

The question of whether the tidal evolution of largely eccentric systems as predicted by S24 (i.e., where the dynamical tides make the eccentricity increase) is justified physically is beyond the scope of this study. We recall that this result directly stems from Z77 formalism, which is an expansion keeping only low-order eccentricity terms. It is possible that pushing the expansion to higher order terms may alter its prediction, in particular the shape of the function f_e, which is crucial for determining the eccentricity evolution (see Appendix D). When combined with the physics of triple systems (in particular ZLK oscillations), the use of the S24 prescription obviously tends to increase of the number of interacting systems, as any system reaching an inner eccentricity e_in ≳ 0.7 through ZLK oscillations and subject to dynamical tides would not circularize and be therefore more likely to interact, according to Eq. (9).

5. Conclusion

In this study we presented the first simulations of triple systems performed with a detailed, on-the-fly stellar code, MESA and a triple system secular evolution code, TRES. The coupling of TRES to MESA was done through AMUSE. This coupling was relatively straightforward due to the fact that TRES is natively developed within AMUSE and thanks to the versatility of this environment. These simulations offer new insights into the physics of triple systems.

The approach we adopted was to primarily investigate the impact of the single star physics on the evolution of triple systems. We kept the same physics for the interactions and investigated how the differences between the stellar evolution predicted by a detailed code such as MESA and that of a fast stellar code based on the H00 fitting formulae (SEBA) can lead to divergences in the predicted evolutionary pathways of triple systems.

We primarily focused on stars with masses M ≥ 50 M_⊙, i.e., outside the range of the P98 grid, on which H00 tracks are based. Above this limit the fitting formulae become increasingly less accurate with increasing masses as they extrapolate the tracks. We compared single star evolutionary tracks obtained with MESA to those obtained with SEBA in the mass range M ∈ [8, 120] M_⊙. We showed that SEBA and MESA tracks deviate significantly at the upper end of our mass range. In particular, the maximum radii reached during the evolution, which are crucial for determining whether components of multiple systems interact, already differ by an order of magnitude or more among the most massive stars with a moderate mass-loss prescription.

Given the importance and uncertainties of stellar winds, especially for the most massive stars, we investigated the impact of the strength of the stellar winds on our results. We showed that the tracks predicted by MESA and SEBA become increasingly divergent with increasing mass-loss, which is a direct consequence of the fact that the H00 tracks do not offer a self-consistent reaction to mass-loss. In particular, with a standard winds combination (i.e., de Jager et al. 1988; Nugis & Lamers 2000; Vink et al. 2001), the predicted maximum radii reached during the evolution differ by up to two orders of magnitude depending on the stellar code in the mass range we considered.

We also provided illustrative examples of triple systems for which the stellar evolution discrepancies between the two codes imply divergences in the evolutionary pathway of the system. We were able to identify examples of divergent pathways for each of the possible evolutionary pathways considered in TRES. These results illustrate the importance of a consistent stellar evolution modeling for predicting the evolutionary pathways of triple systems.

We showed that these discrepancies can be primarily attributed to the different radial expansion predicted by the codes, which is a consequence of the different reaction to mass-loss they offer and the fact that at masses higher than 50 M_⊙, H00 fitting formulae extrapolate the tracks of lower mass stars. Since MESA models expand less at high masses, the size of the parameter space where RLOF can be avoided is increased. At solar metallicity, the limiting period to avoid RLOF with SEBA models is P_min, SeBa ∼ 6 × 10³ days whereas with MESA it is reduced to P_min, MESA ∼ 6 days. This has important consequences for the predicted evolutionary pathways of triple systems. In particular, we demonstrated that the parameter space where triple compact objects can form is increased, which can have significant consequences in the context of the progenitors of GW events. The determination of the BH merger rate requires carrying out a population synthesis study, which was beyond the scope of this project, but could be investigated in a follow-up study.

Using a detailed stellar evolution code also allowed us to self-consistently model some important interactions between the stellar components, in particular the precession of the orbits due to the tidal and rotational distortion of the stars. This precession plays an important role for the evolution of triple systems, as it competes with the precession driven by the ZLK oscillations. In our MESA simulations, we directly obtained the apsidal motion constant from the structure, which allowed us to self-consistently model the precession due to the distortion of the star and its competition with that driven by the ZLK oscillations. In SEBA, in contrast, the apsidal motion constant remains fixed and is updated only when the star transitions to a new evolutionary stage¹². We showed that MS models with self-consistent precession modeling tend to experience more pronounced ZLK oscillations, making them more likely to interact. This is because the apsidal motion constant decreases by orders of magnitude during the MS evolution, reducing the effect of precession caused by stellar distortion, and thereby increasing the influence of ZLK oscillations.

Finally, we illustrated the impact of using another prescription for the dynamical tides, as proposed by S24, on the evolution of a small set of triple systems, and discussed the implication of this prescription when applied to largely eccentric systems.

In summary, using a more self-consistent modeling of the single star evolution alters the predicted evolutionary pathways of triple systems. The great advantage of stellar codes such as SEBA is their computing efficiency, which makes them an essential tool for population synthesis studies. The present study complements that of Kummer et al. (2023), who provided predictions of the main evolutionary pathways of massive triple systems through a population synthesis approach. Quantifying the effect of these discrepancies at the population level requires carrying out a population synthesis study, which would be computationally expensive if MESA were to be used given the large dimensionality of the parameter space of triple systems, but could be the subject of a follow-up study.

Acknowledgments

We sincerely thank the anonymous referee for their constructive feedback and valuable suggestions, which enabled us to improve the manuscript. LS thanks Prof. Georges Meynet (UNIGE), Dr. Max Briel (UNIGE) and Dr. Patrick Eggenberger (UNIGE) for useful discussions, and Dr. Sophie Rosu (UNIGE) for a valuable feedback. LS and SE have received support from the SNF project No 212143. LS has received support from the Netherlands Research School for Astronomy (NOVA). SE and SR have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grand agreement No 833925, project STAREX). EF is support by SNF grant number 200020_212124. ST acknowledges support from the Netherlands Research Council NWO (VIDI 203.061 grants).

References

Appendix A: Calibration of α_ov

In order to perform our α_ov calibration, we computed several MESA models of 50 M_⊙ stars with various α_ov values and compared them to the SEBA model (all models were performed without accounting for mass-loss by stellar winds). As an illustration, we show in Fig. A.1 the end MS evolution of MESA models with α_ov = {0.3, 0.31, 0.335} (the latter value corresponds to that used in (Brott et al. 2011) which was calibrated using data of the FLAMES survey of massive stars) and compare them to the SEBA model.

Fig. A.1. End of MS radial expansion of the 50 M⊙ MESA model with different step-overshoot prescriptions compared with the SEBA model (all models are without stellar winds). The value αov = 0.31 offers the best fit to the SEBA model radial expansion. All simulations were computed without mass-loss. A toy model with decreasing values of αov over time so as to mimic the effect of the δ prescription was also computed (orange curve).

We find that the value α_ov = 0.31 offers the best fit to the SEBA model in terms of the maximal MS radial expansion. We also computed a model variation were we manually reduced the value of α_ov along the MS as to mimic the effect of the δ prescription used in P98. With this model, we started the MS evolution with α_ov = 0.4, which corresponds to the ZAMS value stemming from the δ prescription for this mass (see Fig. 1 in P98) and subsequently reduced it. More precisely, we set the overshooting parameter value to α_ov = {0.4, 0.35, 0.3, 0.25} when the central hydrogen mass fraction reached $X_{\mathrm{c}}=\{X_0,\frac{3}{4}{X}_0,\frac{1}{2}{X}_0,\frac{1}{4}{X}_0\}$ with X₀ = 0.7 the initial hydrogen mass fraction. The MS radial evolution of this model is almost undistinguishable from that of the model with α_ov = 0.31.

The mass of the helium core formed during the MS evolution of the MESA models directly depends on the performed calibration of α_ov. It is interesting to compare the helium core mass resulting from this calibration to that of the SEBA models. Following the H00 fitting formulae, the core mass at helium ignition is given by their Eq. (44), replacing $M_{\mathrm{c}}[L_{\mathrm{BGB}...}(M_{\mathrm{HeF}}){]}^{†}$ with M_c[L_HeI...(M_HeF)]. In the mass range considered in this work, this formula is very well approximated by

$$\begin{array}{c}\hfill M_{\mathrm{c},\mathrm{HeI}}=0.098M^{1.35},\end{array}$$(A.1)

where M is the current mass of the star in solar mass (see comments below their Eq. (44) and (65)). We show in Fig. A.2 a comparison of the He core masses at He burning ignition of MESA and SEBA models as a function of initial mass and with the three different wind scaling factors (0, 0.333, and 1).

Fig. A.2. He core masses at He burning ignition of MESA and SEBA models as a function of mass and with the three different wind scaling factors (0, 0.333, and 1).

We note that when no winds are accounted for, the SEBA and MESA core masses show a consistent trend with the mass, which is expected. An offset between the core mass of the MESA and SEBA models can be noticed, which increases with the mass. The core masses only match for models with M ≤ 20 M_⊙. Notably, the MESA and SEBA core masses of the 50 M_⊙ model, which is the one we used for calibrating α_ov, differ by about 10%. This discrepancy can primarily be attributed to the different treatment of convection between the two approaches. Secondly, the calibration was performed matching the maximum radius, and not the luminosity, which is more directly linked to the mass of the core. One can note looking at the HRD (Fig. 1) that the MESA models in general have a slightly higher luminosity when crossing the Hertzsprung gap, which is consistent with them having slightly larger cores.

Although the core masses are relatively consistent when no winds are accounted for, the discrepancy increases with increasing mass-loss. The reason is that the mass used in Eq. A.1 is the current mass, not the initial mass. At high masses, the stars have already lost a significant portion of their mass at helium ignition, which is why Eq. (A.1) predicts much smaller core masses. Finally, the drop of M_c of the 100 and 120 M_⊙ MESA models is explained by the fact that they loose all of their hydrogen envelope before the end of the MS. At this point, they start losing mass in regions enriched in helium, which reduces the size of their final helium core.

Appendix B: η Carinae*

In this appendix we present a triple system with initial conditions very similar to those of the system simulated in Sect. 4.1, apart from the inner semi-major axis, which is slightly increased (a_in = 1.1 au) compared to the original η Carinae system (a_in = 1 au). It corresponds to the system labeled “η Carinae*” in Table 1. Given the larger inner separation, the octupole term for this system is initially larger than for the original system (ϵ_oct = 0.052), therefore the ZLK oscillations are more pronounced. We simulated this system with SEBA and MESA, and also included a model variation where the stellar evolution is performed with MESA, but the apsidal motion constant is attributed the same fixed value k_AMC = 0.0144 as in the SEBA simulation. The same quantities as in Fig. 4 for this system are shown in Fig. B.1, which illustrates the impact of the precession.

Fig. B.1. Same as Fig. 4, but with a larger inner separation (system labeled “η Carinae*” in Table 1). A model variation (MESA with constant kAMC) is shown in orange.

The MESA model with a self-consistent value of k_AMC is subject to maximal ZLK oscillations, as the precession due to distortion are weaker compared to the precession due to the ZLK mechanism (for t ≥ 0.2 Myr, the apsidal motion constant of the primary retrieved from the MESA structure gets smaller than the constant value used in the other simulations). The MESA model variation shows an evolution more consistent with the system computed with SEBA, even though differences due to the stellar evolution are still present.

Because the ZLK oscillations are larger in the default MESA system, the primary star fills its Roche lobe after only ∼​0.7 Myr, when the inner binary is still in an eccentric orbit. The system in this case experiences an eccentric mass transfer episode. In contrast, ZLK oscillations in the MESA model variation and the system simulated with SEBA are quenched more efficiently due to a higher k_AMC, and the systems in this case avoid RLOF during the eccentric phase. The systems experience mass transfer later (∼2.1 Myr for the MESA model variation and ∼2.5 Myr for the SEBA model), when they have already circularized.

In Table 1, the system of this section is not labeled as strictly divergent, as the nomenclature we followed is that of the paper by (Kummer et al. 2023) (in all cases an inner MT occurs), but given the above discussion, the evolutionary pathways of this system, as predicted by SEBA and MESA, still differ.

Appendix C: Apsidal motion constant evolution of the MESA models

The MS evolution of the apsidal motion constant k_AMC is represented in Fig. C.1 for the whole range of masses of the grid (M ∈ [8, 120] M_⊙), and with Dutch_scaling_factor=0.333. The constant MS value k_AMC = 0.0144 used in the simulations with SEBA is shown as a reference. It can be noted that although the ZAMS value of the models is relatively consistent with the value used by default in TRES, the apsidal motion constant of the MESA models decreases by more than one order of magnitude along the MS evolution. In terms of mass dependence, we note that lower mass models at ZAMS have lower k_AMC than higher mass ones, but a transition occurs around the middle of the MS and at TAMS the higher mass models have a lower k_AMC, consistently with the models of (e.g., Claret 2004). This can be explained by the fact that higher mass models have larger convective core, which results in a larger portion of the total mass being converted to helium at the TAMS and hence a higher density contrast (lower k_AMC).

Fig. C.1. MS evolution of the apsidal motion constant of the MESA models for masses in the range M ∈ [8, 120] M⊙.

Appendix D: S24 model of dynamical tides in eccentric orbits

In this appendix we provide a detailed description of the dynamical tide model used as an alternative to the H02 formalism in the simulations of Sect. 4.7. The prescription we used is the Equation (9) in S24, which consists of a complement to the Z77 formalism. The prescription reads:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\frac{\text{d}}{\text{d}t}\left(I{\mathrm{\Omega }}_{\mathrm{spin}}\right)|}_{\mathrm{Dyn}}=\frac{3}{2}\frac{G{M}^2}{R}{q}^2{E}_2{\left(\frac{R}{a}\right)}^6·\{s_{22}^{8/3}\text{sgn}\left(s_{22}\right)\hfill \\ \hfill & +e^2(\frac{1}{4}{s}_{12}^{8/3}\text{sgn}\left(s_{12}\right)-5s_{22}^{8/3}\text{sgn}\left(s_{22}\right)+\frac{49}{4}{s}_{32}^{8/3}\text{sgn}\left(s_{32}\right))\}\hfill \\ \hfill & {\frac{\text{d}e}{\text{d}t}|}_{\mathrm{Dyn}}=-\frac{3}{4}e{\left(\frac{\mathit{GM}}{R^3}\right)}^{1/2}q{(1+q)}^{1/2}{E}_2{\left(\frac{R}{a}\right)}^{13/2}·\hfill \\ \hfill & (\frac{3}{2}{s}_{10}^{8/3}\text{sgn}\left(s_{10}\right)-\frac{1}{4}{s}_{12}^{8/3}\text{sgn}\left(s_{12}\right)-s_{22}^{8/3}\text{sgn}\left(s_{22}\right)+\frac{49}{4}{s}_{32}^{8/3}\text{sgn}\left(s_{32}\right)).\hfill \end{array}\end{array}$$(D.1)

Equation (D.1) predicts the time evolution of the eccentricity and rotational velocity of a star subject to dynamical tides with radiative damping, i.e., gravity modes excited by the presence of a companion of mass M₂ = qM, where M is the mass of the considered star and q the mass ratio. We recall that the coefficients s_lm are defined as

$$\begin{array}{c}\hfill s_{\mathit{lm}}\equiv (l{\mathrm{\Omega }}_{\mathrm{orb}}-m{\mathrm{\Omega }}_{\mathrm{spin}}){(R^3/GM)}^{1/2}.\end{array}$$(D.2)

When the two components of the binary system are massive enough to be subject to this type of tides (i.e., when they both possess a convective core and a radiative envelope), the equation can be applied to each component, and the evolution of the semi-major axis and eccentricity of the system can be obtained as $\dot{a}={\dot{a}}_1+{\dot{a}}_2$ and $\dot{e}={\dot{e}}_1+{\dot{e}}_2$, adding up the contributions of the effect of the dynamical tides on each component.

Imposing conservation of the total angular momentum${\dot{\mathcal{L}}}_{\mathrm{tot}}={\dot{\mathcal{L}}}_{\mathrm{spin}}+{\dot{\mathcal{L}}}_{\mathrm{orb}}$ and using Eq. (D.1), the equation governing the evolution of the semi-major axis is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\frac{\text{d}a}{\text{d}t}|}_{\mathrm{Dyn}}=-& 3{\left(\frac{\mathit{GM}}{R}\right)}^{1/2}q{(1+q)}^{1/2}{E}_2{\left(\frac{R}{a}\right)}^{11/2}·\hfill \\ \hfill \{\frac{3}{4}& s_{10}^{8/3}\text{sgn}\left(s_{10}\right)·\frac{e^2}{1-e^2}\hfill \\ \hfill +& s_{12}^{8/3}\text{sgn}\left(s_{12}\right)·[\frac{1}{4}\frac{e^2}{\sqrt{1-e^2}}-\frac{1}{8}\frac{e^2}{1-e^2}]\hfill \\ \hfill +& s_{22}^{8/3}\text{sgn}\left(s_{22}\right)·[\frac{1}{\sqrt{1-e^2}}(1-5e^2)-\frac{1}{2}\frac{e^2}{1-e^2}]\hfill \\ \hfill +& s_{32}^{8/3}\text{sgn}\left(s_{32}\right)·[\frac{49}{4}\frac{e^2}{\sqrt{1-e^2}}+\frac{49}{8}\frac{e^2}{1-e^2}]\}.\hfill \end{array}\end{array}$$(D.3)

We implemented Equations (D.1) and (D.3) in TRES and compared the evolutions of the orbital parameters of a few systems they provide to those obtained with the H02 prescription.

In order to interpret the results Eq. (D.1) and (D.3) provide, it is useful to rewrite Eq. (D.1) as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\frac{\text{d}}{\text{d}t}\left(I{\mathrm{\Omega }}_{\mathrm{spin}}\right)|}_{\mathrm{Dyn}}& =\frac{3}{2}\frac{G{M}^2}{R}{q}^2{E}_2{\left(\frac{R}{a}\right)}^6·f_{\mathrm{\Omega }}({\mathrm{\Omega }}_{\mathrm{spin}},{\mathrm{\Omega }}_{\mathrm{orb}},e)\hfill \\ \hfill {\frac{\text{d}e}{\text{d}t}|}_{\mathrm{Dyn}}& =-\frac{3}{4}e{\left(\frac{\mathit{GM}}{R^3}\right)}^{1/2}q{(1+q)}^{1/2}{E}_2{\left(\frac{R}{a}\right)}^{13/2}·f_e({\mathrm{\Omega }}_{\mathrm{spin}},{\mathrm{\Omega }}_{\mathrm{orb}}),\hfill \end{array}\end{array}$$(D.4)

where the functions f_Ω and f_e are defined as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill f_{\mathrm{\Omega }}({\mathrm{\Omega }}_{\mathrm{spin}},{\mathrm{\Omega }}_{\mathrm{orb}},e)& \equiv s_{22}^{8/3}\text{sgn}\left(s_{22}\right)\hfill \\ \hfill & +e^2(\frac{1}{4}{s}_{12}^{8/3}\text{sgn}\left(s_{12}\right)-5s_{22}^{8/3}\text{sgn}\left(s_{22}\right)+\frac{49}{4}{s}_{32}^{8/3}\text{sgn}\left(s_{32}\right))\hfill \\ \hfill f_e({\mathrm{\Omega }}_{\mathrm{spin}},{\mathrm{\Omega }}_{\mathrm{orb}})& \equiv \frac{3}{2}{s}_{10}^{8/3}\text{sgn}\left(s_{10}\right)-\frac{1}{4}{s}_{12}^{8/3}\text{sgn}\left(s_{12}\right)-s_{22}^{8/3}\text{sgn}\left(s_{22}\right)\hfill \\ \hfill & +\frac{49}{4}{s}_{32}^{8/3}\text{sgn}\left(s_{32}\right).\hfill \end{array}\end{array}$$(D.5)

The signs of the functions f_Ω and f_e are crucial as they determine whether the corresponding quantities (Ω and e) increase or decrease with time.

With this set of equations, pseudo-synchronization is reached for a given eccentricity when f_Ω(Ω_spin,Ω_orb,e) = 0. Because the equation governing the evolution of the angular velocity is directly obtained from Z77 formalism without being combined to the (Hut 1981) formalism (as is is done following the H02 derivation), it provides different pseudo-synchronization values than what is obtained with the H02 prescription (where the pseudo-synchronization values are those predicted by Hut 1981).

In Fig. D.1, f_Ω and f_e are represented as functions of the ratio Ω_spin/Ω_orb (for a star of star of mass M = 15 M_⊙, R = 4.7 R_⊙). As f_Ω also depends on the eccentricity, we show its shape for different values of the eccentricity (e = {0, 0.25, 0.5, 0.7, 0.9}). Although the functions f_Ω and f_e depend on both Ω_spin and Ω_orb, their zeros only depend on the ratio Ω_spin/Ω_orb. For this reason we only show the functions for a fixed period of P = 1.2 days, corresponding to Ω_orb = 6.06 × 10⁻⁵ rad/s.

Fig. D.1. Upper panel:fΩ as a function of the ratio Ωspin/Ωorb for different values of the eccentricity. The dashed gray line corresponds to Ωspin = Ωorb. Lower panel:fe as a function of the ratio Ωspin/Ωorb. The dashed gray line corresponds to the value of Ωspin/Ωorb where fe changes sign.

Among the important features in the shape of f_Ω and f_e, we can note the following:

1. The zeros of f_Ω (i.e., the values of the ratio of Ω_spin/Ω_orb for which the star is pseudo-synchronized) increase with the eccentricity, which is the expected behavior for pseudo-synchronization. For circular orbits, the zero of f_Ω is found when Ω_spin/Ω_orb = 1 as expected.

2. The function f_e is not positive everywhere. It changes sign at Ω_spin/Ω_orb ≈ 2.007. This result is independent of the eccentricity as f_e only depends on Ω_spin and Ω_orb and was already obtained by Z77 (see Fig. 2 in the paper). It implies that, for a star rotating above this threshold, the dynamical tides makes the eccentricity increase. This goes against the standard picture that tides always act towards the circularization of the system, but directly stems from Z77 formalism.

3. One can note looking at the curves of f_Ω that for systems with high eccentricities (e ≳ 0.7), the pseudo-synchronization value is larger than the threshold value of Ω_spin/Ω_orb = 2.007 at which f_e changes sign. This implies that for any pseudo-synchronized star in a system with e ≳ 0.7, the dynamical tides makes the eccentricity of the system increase.

It is useful to obtain the value of the ratio Ω_spin/Ω_orb for which the star is pseudo-synchronized as a function of the eccentricity. We can thus define the corresponding function of the eccentricity f_sync(e) in the following way:

$$\begin{array}{c}\hfill f_{\mathrm{sync}}(e)\equiv max{\frac{{\mathrm{\Omega }}_{\mathrm{spin}}}{{\mathrm{\Omega }}_{\mathrm{orb}}}|}_{f_{\mathrm{\Omega }}({\mathrm{\Omega }}_{\mathrm{spin}},{\mathrm{\Omega }}_{\mathrm{orb}},e)=0}.\end{array}$$(D.6)

The maximum notation is used as there is a small interval of eccentricities for which f_Ω has several zeros. As can be seen in Fig. D.1, the graph of f_Ω is similar to that of a polynomial function of degree 3, which may have one or three zeros. It turns out that for eccentricities in the interval e ∈ [0.524, 0.566], the function f_Ω has three zeros. In order to correctly define the function f_sync for any eccentricity, we take the maximum of the three values in this interval.

Because the dependence in R/a is larger in the eccentricity equation than in the angular velocity equation, pseudo-synchronization usually occurs faster than circularization. It is therefore useful to determine whether the action of the tides on a pseudo-synchronized star will tend to make the eccentricity increase or decrease. To this purpose, we define the function f_circ(e, Ω_orb) as:

$$\begin{array}{c}\hfill f_{\mathrm{circ}}(e,{\mathrm{\Omega }}_{\mathrm{orb}})\equiv {f_e({\mathrm{\Omega }}_{\mathrm{spin}},{\mathrm{\Omega }}_{\mathrm{orb}})|}_{\frac{{\mathrm{\Omega }}_{\mathrm{spin}}}{{\mathrm{\Omega }}_{\mathrm{orb}}}=f_{\mathrm{sync}}(e)}.\end{array}$$(D.7)

f_circ corresponds to the value of the function f_e for any pseudo-synchronized star. It is eccentricity dependent as the value of the ratio Ω_sync/Ω_orb for which pseudo-synchronization is reached depends on the eccentricity. f_circ still depends on Ω_orb (or Ω_spin) as the pseudo-synchronization condition only constrains the ratio Ω_sync/Ω_orb, but does not constrain both Ω_sync and Ω_orb. Since the zeros of both f_e and f_Ω only depend on the ratio Ω_sync/Ω_orb, it can be shown that the zero of f_circ is the same for any value of Ω_orb, and is reached for a single value of the eccentricity of e ≈ 0.659.

We show in Fig. D.2f_sync and f_circ as functions of the eccentricity. We represent f_circ for different values of the orbital angular velocity.

Fig. D.2. Upper panel:fsync as a function of the eccentricity. Lower panel:fcirc as a function of the eccentricity for different values of Ωorb. The dashed gray line corresponds to the value of the eccentricity where fcirc changes sign.

The following comments can be made regarding Fig. D.2:

1. The visible discontinuity of f_sync at e ≈ 0.566 can be imparted to the fact that in the interval e ∈ [0.524, 0.566], f_Ω possesses three pseudo-synchronization values, and we take the maximum of the three. f_sync shows an expected behavior: it increases with the eccentricity and reaches unity at e = 0, which is interpreted as standard synchronization.

2. f_circ also presents a discontinuity at e ≈ 0.566, which is not surprising given that it is defined with respect to f_sync.

3. The most striking feature of f_circ is that it changes sign at e ≈ 0.659. The interpretation of this result is that according to this model, dynamical tides tend to increase the eccentricity of any pseudo-synchronized star in an orbit with e ≳ 0.659. This can in practice never make the eccentricity reach unity as the increase of e also implies an increase of a by conservation of angular momentum, which will at some point quench the effect of the tides, given the strong dependence in a in the eccentricity evolution.

4. Given this interpretation of f_circ, we can understand that the tidal evolution of systems in highly eccentric orbits as predicted by S24 may differ significantly from what is predicted following the H02 formalism. In the case of triple systems where the inner binary is subject to large eccentric ZLK oscillations, where the inner eccentricity can reach extreme values, the dynamical tides in this case will not tend to make the eccentricity decrease but increase even more, which can have important consequences for the evolution of the system (e.g., RLOF in an eccentric orbit).

All Tables

All Figures

	Fig. 1. HRD obtained with MESA (solid lines) and SEBA (dashed lines). Isoradius lines are represented by gray dashes. (a) HRD for stars in the mass range M ∈ [8,  120] M⊙. (b) HRD for stars in the mass range M ∈ [50,  120] M⊙.
In the text

	Fig. 2. Same as Fig. 1b, but with a wind scaling factor three times higher.
In the text

	Fig. 3. Upper panel: Maximum radius reached during the MS by the MESA and SEBA models for different values of the wind scaling factors (0, 0.333, and 1). Lower panel: Same as upper panel, but for the maximum radius reached during the whole evolution and for wind scaling factors of 0.333 and 1.
In the text

	Fig. 4. Evolution of the η Carinae system simulated with MESA and SEBA. Upper panel: Inner eccentricity. Middle panel: Radius expansion of the primary (dashed line) and secondary (dash-dotted line). The Roche lobe radius of the primary is shown as a solid line. Lower panel: Evolution of the apsidal motion constant of the primary.
In the text

	Fig. 5. Evolution of the system labeled “Tertiary MT” in Table 1. Upper panel: Radius expansion of the primary (dashed line) and Roche lobe radius (solid line). Lower panel: Same as upper panel, but for the tertiary.
In the text

	Fig. 6. Evolution of the system labeled “Dyn. Destab” in Table 1. Upper panel: Mass evolution of the primary. Middle panel: Radius and Roche lobe radius evolution of the primary. Lower panel: Evolution of the ratios $\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}$ and ${\frac{a_{\mathrm{out}}}{a_{\mathrm{in}}}|}_{\mathrm{crit}}$.
In the text

	Fig. 7. Same as Fig. 6 for the system “Orbit unbound” in Table 1.
In the text

	Fig. 8. Same as middle panel of Fig. 7 for the system “TCO” in Table 1.
In the text

	Fig. 9. Evolution of the system labeled “ZLK & precession” in Table 1. Upper panel: Radius expansion of the primary (dashed line) and Roche lobe radius (solid line). Middle panel: Same as upper panel, but for the tertiary. Lower panel: Evolution of the apsidal motion constant of the primary.
In the text

Fig. 10. Comparison of the evolution of systems in moderately eccentric orbits as predicted by the H02 and the S24 prescriptions. One system is subject to ZLK oscillations (blue and magenta lines), the other is not (violet and orange lines). The scale of the horizontal axis changes at t = 2.6 Myr and is chosen so as to focus on the part of the evolution where tidal effects are important. Upper panel: inner eccentricity evolution. Lower panel: Radius and Roche lobe radius evolution of the primary.

In the text

	Fig. 11. Same as Fig. 10, but for systems with ein, max ≳ 0.8. Additionally the evolution of the inner separation is shown in the lower panel.
In the text

Fig. A.1. End of MS radial expansion of the 50 M⊙ MESA model with different step-overshoot prescriptions compared with the SEBA model (all models are without stellar winds). The value αov = 0.31 offers the best fit to the SEBA model radial expansion. All simulations were computed without mass-loss. A toy model with decreasing values of αov over time so as to mimic the effect of the δ prescription was also computed (orange curve).

In the text

	Fig. A.2. He core masses at He burning ignition of MESA and SEBA models as a function of mass and with the three different wind scaling factors (0, 0.333, and 1).
In the text

	Fig. B.1. Same as Fig. 4, but with a larger inner separation (system labeled “η Carinae*” in Table 1). A model variation (MESA with constant kAMC) is shown in orange.
In the text

	Fig. C.1. MS evolution of the apsidal motion constant of the MESA models for masses in the range M ∈ [8, 120] M⊙.
In the text

	Fig. D.1. Upper panel:fΩ as a function of the ratio Ωspin/Ωorb for different values of the eccentricity. The dashed gray line corresponds to Ωspin = Ωorb. Lower panel:fe as a function of the ratio Ωspin/Ωorb. The dashed gray line corresponds to the value of Ωspin/Ωorb where fe changes sign.
In the text

	Fig. D.2. Upper panel:fsync as a function of the eccentricity. Lower panel:fcirc as a function of the eccentricity for different values of Ωorb. The dashed gray line corresponds to the value of the eccentricity where fcirc changes sign.
In the text