© The Authors 2024

1. Introduction

Star systems containing two or more stars are common in the Universe. Observations show that about half of solar-type stars reside in binaries or higher-order multiples, with both the multiplicity fraction and number of companions increasing with the mass of the star (Tokovinin 2014; Moe & Di Stefano 2017; Offner et al. 2023). As a consequence, the majority of massive stars (M ≥ 8 M_⊙) are expected to have at least two companions. A significant fraction of these systems are likely to interact during the lifetimes of the stars, which can drastically affect the evolution of the stellar components (Sana et al. 2012).

Stellar objects in binary systems follow fixed orbits as long as the system does not experience any external perturbations or internal changes. In comparison, the addition of a third body introduces effects that can alter the configuration of the system over secular timescales due to purely gravitational interactions (von Zeipel 1910; Kozai 1962; Lidov 1962; Szebehely & Zare 1977; Naoz et al. 2013). Stellar evolution combined with these three-body dynamics allows greater variation of evolutionary pathways when compared to binary systems. This includes mass transfer from the tertiary onto the inner binary (Portegies Zwart & Leigh 2019; Gao et al. 2022; Dorozsmai et al. 2024; Kummer et al. 2024), triple common envelope (Glanz & Perets 2021), and dynamical destabilisation (Perets & Kratter 2012; Hamers et al. 2022; Toonen et al. 2022), which can lead to either the ejection of a component or a collision between two bodies.

A triple system destabilises when the ratio between the outer and inner semi-major axis decreases, which occurs due to mass loss in the inner binary from stellar winds (Perets & Kratter 2012), mass transfer in the inner binary when the donor becomes less massive than the accretor (Hamers et al. 2022), supernova kicks in either the inner or outer binary, or from the gravitational perturbation of a passing object (Michaely & Perets 2020; Hamers et al. 2022; Stegmann et al. 2024). The destabilisation of triples due to mass loss is known as triple evolution dynamical instability (TEDI; Perets & Kratter 2012), and may be responsible for a large population of Galactic stellar collisions, in addition to being an efficient pathway for producing single stars through the ejection of one of the bodies. Though the potential outcomes of unstable triples are known from theory, a complete statistical overview of physical systems is lacking, with a few recent studied having provided more insight into these numbers. Toonen et al. (2022, herafter referred to as TBZ22) performed N-body simulations of low-mass triple systems on the edge of stability, and find that up to ∼45% of these systems disintegrate into an ejected single star and a remaining binary, with the single star often reaching velocities of tens of kilometres per second. Additionally, between 13% and 24% of the unstable triples produce a collision between two bodies, wherein the vast majority happen between two main sequence (MS) stars. Hamers et al. (2022, hereafter HPTN22) performed a similar study of systems that undergo TEDI, and likewise found that unstable triples can produce stellar collisions and unbound stars at a Galactic rate of 10⁻⁴ yr⁻¹. In both of these studies, the mass of the primary star was sampled from a Kroupa distribution (Kroupa 2001). TBZ22 excluded stars with masses greater than 7.5 M_⊙, while HPTN22 sampled the full mass range. Consequently, a similar study has not yet been produced for triple systems with exclusively massive primaries. As already mentioned, the vast majority of massive stars are found in triples or higher-order multiples, which means that a full consideration of triple dynamics is crucial in the context of understanding the evolutionary pathways of massive stars. Their shorter lifetimes and distinct orbits combined with strong stellar winds means that the orbits of massive multiples are likely to change on shorter timescales than their lower-mass counterparts, potentially resulting in a higher fraction of triples that become unstable within a Hubble time.

Unstable low-mass triples may lead to merger products and ejected single stars and binaries. If the same holds true for massive stars, triples may partly explain the large numbers of observed massive runaway stars (Gies 1987; Stone 1991; Oey et al. 2018; Renzo et al. 2019; Stoop et al. 2023), which cannot be fully explained by supernovae in binaries or dynamic interactions in clusters (de Wit et al. 2005; Gvaramadze et al. 2012). Additionally, as massive stars end their lives as compact remnants – neutron stars (NSs) or black holes (BHs)–, massive triple systems that survive one or more supernova explosions may yet get pushed to the stability limit, after which they disintegrate and produce outcomes, such as runaway BHs and NSs, and mergers or collisions between compact remnants or between a compact remnant and a star, which may lead to phenomena such as tidal disruption events and gravitational wave sources.

In this study, we follow up on TBZ22 by focusing exclusively on triples with a massive primary star. With a combination of numerical techniques for both secular orbit-averaged evolution and N-body dynamics, we follow the evolution of triples starting from a stable hierarchical configuration at the zero-age main sequence (ZAMS), through the regime of dynamical instability, and until the potential disintegration of the system. In the following section, we present a brief overview of three-body dynamics and the topic of stability in hierarchical triples (Sect. 2), followed by an overview of the methodology in which we discuss the choices related to our numerical setup (Sect. 3). We then present the results of our simulations and the statistical properties of the outcomes, including the components involved in collisions and the velocities of ejected stars (Sect. 4). Finally, we discuss our results in the context of observable phenomena (Sect. 5), before presenting our conclusions as a summary of our results along with suggestions for future work (Sect. 6).

2. Three-body dynamics

2.1. Hierarchical triples

As opposed to the case of two gravitationally bound objects, the addition of a third component drastically changes many aspects of the system, one of which is the stability. A two-body system that does not experience any external or internal changes will remain in its orbital configuration indefinitely, and the state vectors of the components can be calculated using analytical formulae. In a three-body system, no analytical solution exists, and the dynamics must therefore be solved using numerical methods. Most three-body systems are chaotic, meaning their outcome is sensitive to the initial conditions, wherein the degree of sensitivity depends on the configuration of the system (Boekholt et al. 2020).

Triples in the universe have so far been found in a configuration consisting of an inner binary whose centre of mass is in a wide orbit with a tertiary companion. This configuration, known as a hierarchical triple, occurs when the separation between the components in the inner binary, denoted as a_in, is significantly smaller than that of the outer binary, denoted by a_out. The larger the ratio a_out/a_in, the stronger the hierarchy and the longer the system can remain stable. The hierarchical configuration allows the system to be described as two nested Keplerian orbits with semi-major axes a_in and a_out, eccentricities e_in and e_out, arguments of periapsis g_in and g_out, and longitudes of ascending nodes Ω_in and Ω_out. The two orbits can be inclined with respect to each other. This mutual inclination, denoted as i, as is defined as the angle between the orbital angular momentum vectors of the inner and outer orbits. The components in a hierarchical triple are denoted as the primary, secondary, and tertiary, with m₁ ≥ m₂ labelling the masses of the primary and secondary, and m₃ the tertiary.

A triple with a stronger hierarchy (higher value of a_out/a_in) can remain in its configuration for longer. However, stability is determined by a combination of orbital and stellar parameters, and not just the separations alone. There have been several proposed prescriptions for estimating the stability of hierarchical triple systems, one of which is the Mardling & Aarseth (1999) criterion, which provides an analytical estimate for the critical semi-major axis ratio, defined such that a triple is classified as unstable if the actual semi-major axis ratio $\frac{a_{\text{out}}}{a_{\text{in}}}$ is smaller than the critical ratio $\frac{a_{\text{out}}}{a_{\text{in}}}{|}_{\mathrm{crit}}$, where

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

Here, q_out is the outer mass ratio defined as q_out ≡ m₃/(m₁ + m₂).

2.2. The Kozai-Lidov mechanism

The von Zeipel-Lidov-Kozai (ZKL; Lidov 1962; Kozai 1962) mechanism is one of the most important effects in a hierarchical triple, as it can drastically affect the evolution of the system and its stellar components. ZKL is a secular process that introduces periodic changes in the mutual inclination and inner eccentricity of a hierarchical triple system. The physical mechanism behind the process is described in detail in e.g. Naoz (2016), but can be summarised as a torque between the two orbits, allowing them to exchange angular momentum while conserving energy. This causes the orbits to alter their shapes and orientations. The ZKL process is also known as the ZKL cycle, as it can cause the system to oscillate between high inclination and low inner eccentricity, to a low inclination and high eccentricity, in a periodic cycle.

3. Methodology

We combine two different methods for studying unstable triples: orbit-averaged dynamics with stellar evolution (SE) for long-timescale simulation of orbits and stellar objects, and a direct N-body approach. In the orbit-averaged technique, differential equations describing the time evolution of the orbital elements are solved, as opposed to the positions and velocities of the bodies. As these change on timescales significantly longer than the dynamical timescales, systems can be simulated for long periods at a fraction of the computational cost when compared to direct simulations. In this study, a population of massive triple system are evolved from ZAMS using the triple population synthesis code TRES (Toonen et al. 2016; Kummer et al. 2023). The subset of systems that cross the stability limit (Eq. 1) are then further evolved using a dynamical approach in an N-body code.

At the onset of instability, the orbit-averaged approach of secular evolution can no longer be utilised, as this technique is only valid for strongly hierarchical systems. On the stability limit, perturbations on the inner orbit from the outer star may occur on timescales on the order of the dynamical timescales, meaning that these changes are not captured by the orbit-averaged approach. To further study the evolution of these now-unstable triple systems, they therefore have to be simulated dynamically using an N-body approach.

3.1. Evolving triples with TRES

TRES is a code developed with the goal of rapid evolution of large populations of hierarchical triples. For a full overview of the internals of TRES, see Toonen et al. (2016). In TRES, SE is included as single star evolution implemented using the fast SE code SeBa (Portegies Zwart & Verbunt 1996; Toonen et al. 2012), which is based on fitted SE tracks from Hurley et al. (2000). SeBa provides fast evaluations for stellar parameters as a function of the initial stellar mass and time. Evolution of the orbital parameters is implemented in TRES as a system of ordinary differential equations (ODEs) derived using time-independent perturbations of the Hamiltonian in the semi-major axis ratios (Ford et al. 2000). The set of ODEs includes the inner and outer orbital parameters, in addition to total angular momentum, mutual inclination, and spin frequencies of each of the three stars. TRES includes both the lowest order expansion of the Hamiltonian, the quadrupole-level, and the octupole-level, in addition to gravitational wave emission, precession, and tidal friction. The code evolves triple systems until a specific stopping condition is reached. In the simulations studied in this work, the specific stopping conditions are:

• Primary or secondary RLOF.

• Tertiary RLOF.

• System crosses the stability limit defined by the stability criterion in Eq. (1).

• Inner orbit becomes unbound due to core-collapse supernova in either the secondary or primary.

• Outer orbit becomes unbound due to core-collapse supernova in any of the three components.

• The system reaches Hubble time without any of the aforementioned stopping criteria being triggered.

Here, we present a summarised version of the choice of initial conditions for the triple population that produced the unstable triples, based on Kummer et al. (2023). For a full description, see the aforementioned paper. The initial population of triples are generated by sampling distributions of the component masses and orbital elements. The generated systems consist of co-evolving stars at the ZAMS bound together in an initially dynamically stable configuration. Given that these are systems containing massive stars, there are large uncertainties in the distribution of stellar and orbital parameters due to the low number of observed counterparts. To compensate for this, Kummer et al. (2023) produces sets of varied initial conditions by applying a series of assumptions. Specifically, they define their fiducial model which follows a Kroupa distribution (Kroupa 2001) for the primary mass, while the secondary and tertiary masses follow from the inner and outer mass ratios q_in ≡ m₁/m₂ and q_out ≡ m₃/(m₁ + m₂), which are both sampled from a flat distribution in [0.1, 1]¹. For the orbital parameters, the semi-major axes a_in and a_out are both sampled from various distributions such that systems that are initially overflowing their Roche lobe at the ZAMS are discarded, in addition to systems that are far enough apart for them to be weakly gravitationally bound. This gives limits of a_{in, out} ∈ [5, 5 × 10⁶] R_⊙. The eccentricities e_in and e_out are both sampled from thermal distributions in the range [0, 0.95] (Ambartsumian 1937; Heggie 1975; Moe & Di Stefano 2017; Hwang 2023). The arguments of pericentre g_in and g_out, along with the mutual inclination, i_mut were sampled using uniform distributions in [ − π, π] and [0, π] respectively. Finally, stellar rotational velocities were sampled using the distribution described in Hurley et al. (2000).

3.2. From orbit-averaged to dynamical evolution

As the orbit-averaged approach utilised by TRES has no information about the positions and velocities of each body at any given time, and consequently, we do not have full initial conditions for the N-body integrator a priori. Therefore, when initialising a triple evolved with TRES in the dynamics code, we need to make a choice for the anomaly, i.e. the initial angular positions in the orbits at the start of the simulation. When a triple system is on the edge of a stability, the hierarchical configuration might be lost within a few outer orbits (Toonen et al. 2022), which means that the initial positions of the bodies in the system is likely to play a larger role in the time it takes for the system to disintegrate and consequently the final outcome. Therefore, randomly positioning each system in the orbit is likely to give us only a small sample of the potential outcomes of the population. To mitigate this, we make ten instances of each system, and distribute them according to a uniform distribution of the true anomaly θ ∈ [0, 2π], where an anomaly of 0 corresponds to periapsis. In other words, each system has ten different variations, where the only difference is the initial true anomaly of the outer orbit. TBZ22 performed similar initialisations of the systems in their study. Our final statistics include the outcomes from all realisations of all systems. The reason for varying the outer anomaly as opposed to the inner, is that it is the outer orbit that has the most effect on the stability of the system, as seen by the dependency on e_out in the stability criterion (Eq. 1).

3.3. Simulating dynamics with Syzygy.jl

For this study, we utilise the package Syzygy.jl², written by the main author in the Julia programming language³. The package utilises the DifferentialEquations.jl (DiffEq) (Rackauckas & Nie 2017) ecosystem to solve the governing differential equations, providing a combination of high performance, stability, and flexibility, as it allows the users to choose from a large variety of optimised ODE solvers. In this study, we utilise the 8th order explicit adaptive Runge-Kutta-Nyström method with absolute and relative error tolerances of 10⁻¹⁰. We found that this combination of algorithm order and error tolerances provided the optimal balance of speed and long-term stability. We include a set of criteria to stop the simulations in the case of specific events. These stopping conditions are implemented as callbacks, which, in the DiffEq ecosystem, is defined as code that is injected by the user into the solver algorithms. Callbacks allow the solver to handle events, including discontinuities, in a safe and stable way, though in the case of stopping conditions they are simply functions that, at each time step, takes in the state of the integrator, checks one or more criteria and terminates the integration in the case of the criterion being fulfilled. In this study, we implement four stopping conditions: (1) Collision between two bodies: if any two bodies have overlapping radii, we stop the simulation and flag the outcome as a collision with the two components (see Sect. 5.4 for a discussion). (2) Unbound: if a body gets ejected due to a close passage, we stop the simulation and flag the outcome as escaper. To check whether a body, at any point, has become unbound, we implement the constraints from Standish (1971), which consists of three criteria that must be fulfilled: (1) the body is a given distance away from the centre of mass of the other two bodies. We set this distance to be 100 times the oscullating semi-major axis of the remaining two bodies. (2) The body is moving away from the centre of mass of the remaining components, and (3), the body has positive energy. It may occur that the third criteria is not satisfied even when a body has been ejected far from the other components, meaning that a tiny background gravitational potential (as would be present in any astrophysical gravitational system) would unbind the ejected star. As a consequence, we flag the body as a drifter if only the first two criteria are satisfied, and the body is a distance of 1 parsec from the remaining binary centre of mass. (3) To set a limit on the amount of computational resources we use, and to avoid running the simulations for unpractical timescales, we stop the simulation if a system is simulated for longer than a preset maximum CPU time. For the simulations in this study, we set the max CPU time to be 10 hours. (4) The system reaches the maximum allotted simulation time without any of the aforementioned conditions being triggered.

Syzygy.jl also includes flags for checking whether a system remains in a hierarchical configuration, or if it has entered a so-called democratic phase, wherein the hierarchy has broken down. These are not stopping conditions, but rather bookkeeping for later analysis. This is an important inclusion when analysing the systems that eject one of the components. There are two pathways that can cause a body to be ejected from the triple: (1) ZKL effects in triples close to the instability limit can cause extreme eccentricity excitations in the outer orbit, resulting in a small outer periastron distance where the tertiary can reach and exceed the escape velocity of the system due a nudge from the inner binary. (2) The system enters a democratic phase, in which close passages can occur between any pair of stars, potentially accelerating one of the components to escape velocities. For the systems that eject a component, we label the outcome as either a Hierarchical or Democratic-X ejection, depending on whether the democratic flag was raised at any point during the simulation, with X being either 1, 2, or 3, denoting which component was ejected. We implement three different methods for checking whether the triple remains hierarchical:

1. Check whether the pair of bodies with the smallest separation is the initial inner binary. This test is quick, but it might be prone to false positives in systems with highly eccentric outer orbits.

2. Check if the original inner binary is still the binary with the smallest semi-major axis. This check is more computationally expensive, as it requires the calculations of the total energy in each pair to check if they are bound, followed by the calculation of the semi-major axis. It is, however, more robust than the first check.

3. Check if any of the components in the inner binary have hyperbolic orbits (e ≥ 1). This check should in theory only be true if the system is not hierarchical; however, the osculating orbital elements might be prone to instantaneous deviations from the expected values due to numerical inaccuracies.

We label a system as being democratic if any of these three criteria are satisfied. Additionally, we implement the ability to raise flags if a body overflows its Roche lobe. For the tertiary component in a hierarchical triple, its Roche lobe is calculated by approximating the system as a binary, where the companion mass is equal to the total mass of the inner binary. However, unstable systems can lose their hierarchy (Toonen et al. 2022), which begs the question of how to check for Roche lobe overflow (RLOF) if the system is in a democratic phase. To implement this, we include two callbacks to check for RLOF, one for while the system is hierarchical, and one for if the hierarchy has broken down. We use the previously mentioned criteria to check if the democratic flag has been raised, and then use the appropriate check to test for RLOF. For a system in which the democratic flag is raised, we go through each body in the system, find the nearest other body, and use the generalised volume-equivalent Roche lobe radius (Eggleton 1983; Sepinsky et al. 2007):

$$\begin{array}{c}\hfill R_L^{\mathrm{Egg}}(t)=D(t)\frac{0.49q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})},\end{array}$$(2)

where D(t) is the distance between the two bodies at a time t, and q is their mass ratio. If the radius of a star is ever greater than R_L^Egg(t), we store the time and the identities of the mass-transferring components. We emphasise that the RLOF check is simply a flag, and not a stopping condition. Realistically, mass transfer during RLOF would alter the dynamics and likely change the outcome of these simulations. However, eccentric mass transfer is still poorly understood (see Sepinsky et al. 2007; Church et al. 2009; Hamers & Dosopoulou 2019; Dosopoulou et al. 2017), and prescriptions suitable for N-body simulations are nonexistent. Therefore, we simply flag these systems to get an overview of the occurrence of RLOF in unstable massive triple systems, and make a note that these systems might have a different outcome if mass transfer was included. In addition to omitting mass transfer, we also do not include tidal interactions between the stellar objects in the N-body simulations.

3.3.1. Stellar evolution

In this study, we simulate the triple population using two variations of dynamical evolution: a pure N-body approach in which only Newtonian gravity is included, and an N-body + SE model. The latter incorporates the effects of SE through changes in mass and radius. As showed in TBZ22, triple stars that reside on the stability limit may remain hierarchical for thousands and up to tens of millions of crossing times⁴, which means that some systems will need to be simulated for timescales approaching or exceeding the nuclear timescales before destabilisation occurs. As a consequence, the effect of SE should be taken into account to get a more realistic view of the dynamics of the system. Since this study is looking at massive triples, SE is likely to play an even bigger role in the dynamics due to the substantially smaller nuclear timescales and higher mass loss when compared to their low-mass counterparts.

In order to incorporate SE into the direct simulation without adding substantial computational cost, we use SeBa to pre-calculate the single SE for each component in each system. The system is evolved until every component becomes a remnant (white dwarf, NS, or BH). We implement SE into the gravity simulation by constructing an interpolator for each stellar parameter, followed by setting up a callback which takes in the current state of the integrator and interpolates the parameters at each time step. The masses, radii, and stellar types of each object are then updated with the interpolated values. Using a callback allows for the ODE system to remain stable in the case of sudden mass loss due to e.g. a supernova (SN). We also include instantaneous velocity kicks on the bodies that experience supernovae using the SN kick model of Verbunt et al. (2017), combined with a fallback model (Fryer et al. 2012) to reduce the kick velocities of higher mass BHs.

4. Results

4.1. Secular evolution

In a typical population of massive triples, with minimum primary masses of 10 M_⊙, there are several potential evolutionary channels that the systems may evolve through. As shown by Kummer et al. (2023), 67.3% of the systems undergo mass transfer in the inner binary, 29% experience unbinding in either the inner or outer orbit due to a supernova, 1.8% undergo mass transfer from the tertiary onto the inner binary, while 1.7% become dynamically unstable due to winds in the inner binary or supernovae. At the onset of instability, the triple population contains 71% MS stars and 26% BHs, with the remaining percentages consisting of post-MS stars and neutron stars. In the case of our simulations with only gravity, this distribution of stellar types remains fixed, while the inclusion of SE means that the stars are evolving as they are simulated with the N-body code, resulting in a final distribution of stellar types that is distinct from the original (Fig. 1).

Fig. 1. Distributions of stellar types for the systems that become dynamically unstable from secular evolution. The purple bars indicate the stellar types at the end of the simulations in which SE is included.

4.2. Duration of dynamically unstable phase

The maximum simulation time for triples is set to 10⁷ × P_{out, init}, where P_{out, init} is the period of the outer orbit at the time of destabilisation. However, the majority of the simulations were terminated before reaching this time due to another stopping condition being triggered. Figure 2 shows cumulative histograms of the simulation time divided by the initial outer period (N ≡ t_end/P_out, init). Half of all the systems have not triggered any stopping condition after N = 1000, and 20% of systems remain stable after a simulation time of more than N = 30 000. These results show that triples are likely to remain hierarchical for long timescales despite being classified as unstable, an observation that coincides with that of TBZ22.

Fig. 2. Duration of the unstable phase for the most prominent outcomes. The x-axis shows the total simulation time expressed as a multiple of the outer orbital period at the time of destabilisation, denoted N. The dotted line shows the distribution of N for all outcomes, while the dashed line shows the duration in units of crossing times. Panels (a) and (b) show N for the pure gravity and SE simulations, respectively.

4.3. All outcomes

For the simulations in which only gravity is taken into account, we find that the systems have different potential outcomes; collision, drifter, escape, CPU time, and numerical issues (labelled as “Numerics” in Fig. 3). When SE is taken into account, we find that the systems experience the same outcomes with the addition of ionisation due to supernova kicks (Fig. 3). We find that 40% of systems experience a collision between two bodies when only gravity is considered. This percentage is slightly lower when we include SE, at 35%. For the system that produce an escaper, 31% eject a body in the pure gravity case, while this rate increases to 40% when SE is also taken into account. In the case of the drifters, we find a large discrepancy between the two simulation suites, with 19% of the systems experiencing this outcome when only gravity is involved, while the inclusion of SE reduces this to 3.5% of systems. Ionisation, in which all three stars become unbound from each other, only occurs when SE is considered and when one or more supernovae take place. We find that 20% of triples undergo ionisation with the inclusion of SE. A smaller fraction of simulations is terminated due to the CPU time stopping condition being triggered. This happens in 8.5% and 0.6% of systems for the gravity and SE suites respectively.

Fig. 3. Fractions of the different outcomes of the dynamical simulations for both the pure N-body and the SE models. The maximum simulation time was 107 × Pout, init.

We also find that only 0.6% of triples reach the CPU time limit when we include SE, compared to only 8.5% in the case of pure gravity. This is likely due to the inclusion of stellar winds. If a system has a large semi-major axis ratio, a_out/a_in, while still being classified as unstable (due to a high outer eccentricity, significant mass ratio, or a mutual inclination close to 90°), the simulation will be limited by the small time step in the inner orbit. Consequently, the code needs a substantial amount of CPU time to complete an outer orbit. In the N-body model, a system in this configuration will remain so until it disintegrates, which, due to the small time step, might take longer than the allotted CPU time. With the inclusion of stellar winds in the SE model, the inner binary may widen due to mass loss, which decreases the semi-major axis ratio and increases the time step used by the ODE solver. These systems will get pushed further into the instability region, and will therefore not only be faster to simulate, but are also likely to disintegrate more rapidly.

4.4. Collisions

Between 35% and 40% of the destabilising systems experience a collision between two bodies. We find that the collisions happen between the primary and the secondary components in 96 − 98% of the colliding systems, showing significant difference between the pure gravity systems and the systems with SE included. The stellar types of the collision pairs vary between the two models, with the majority of collisions in both models taking place between two MS stars (Fig. 4). Approximately 12.5 − 28% of collisions occur between a MS and a post-MS star, with the larger percentage occurring when SE is taken into account. When only gravity is included, we find that about 5% and 3.2% of collisions happen between post-MS + BH and MS + BH respectively, while these percentages are reduced to 0.86 and 0.72 for the SE model.

Fig. 4. Percentages of the stellar types involved in collisions.

Following the collision, we calculate the properties of the newly formed binary using assumptions of a completely inelastic collision with mass conservation. We find that both simulation suites produce binaries with semi-major axes similar to the outer semi-major axis at the onset of instability. (Fig. 5).

Fig. 5. Properties of the newly formed binaries following the collision between two bodies in the original triple. The left figure shows the cumulative distribution of the semi-major axis, while the right shows the eccentricity. The dotted and dashed lines represent the semi-major axes and eccentricities at the onset of instability.

We observe a similar trend in the eccentricities of the newly formed binaries, namely that they tend to have eccentricities similar to the initial outer eccentricity of the triple progenitor. Given that the vast majority of collisions occur between the primary and secondary components, this result is to be expected. The new eccentricity distributions are also found to be sub-thermal⁵, which differs from the results of binary-single scattering experiments (Heggie 1975) by producing binaries that are generally less eccentric. A thermal eccentricity distribution is often used as the initial conditions of binary systems in population synthesis calculations, such as in SeBa. These results show that binaries formed from disintegrating triples are likely to deviate from this distribution.

4.5. Ejections

Between 31% and 41% of the unstable systems result in the ejection of one of the components. When only gravity is included, the majority of the ejected bodies are MS stars and BHs, with percentages of 71% and 26% respectively (Fig. 6). When taking SE into account, these percentages are close to being equal, at 43% and 42% respectively, with the remaining escapers consisting mainly of neutron stars.

Fig. 6. Fractions of the different stellar types that are ejected from the system during the unstable phase.

Ejected bodies are categorised as either escapers or drifters, depending on the criteria as explained in Sect. 3.3. One of the major differing characteristics of these bodies are their velocities (Fig. 7). In both simulation suites, we find that drifters have consistently lower velocities than their escaper counterparts, with a mean velocity of 0.3 km/s. This is to be expected from the definition of the drifters: bodies that have been ejected to substantially large distances away from the remaining binary, but are still loosely bound to the system. When it comes to the escapers, the velocities differ more between the two models, with more high-velocity escapers being produced when SE is taken into account. The tail of higher velocity escapers in the SE model arise from the inclusion of supernova kicks in the simulations, and consequently the highest velocity objects are exclusively neutron stars. However, as seen in Fig. 8, ejected neutron stars in the pure N-body model still reach velocities on the order of tens of kilometres per second without the help of a supernova kick. For the MS stars, the median velocity is 3.6 km/s, with a standard deviation of 12.5 km/s, which shows that dynamical interactions can produce single stars with a wide range of velocities.

Fig. 7. Velocities of the components labelled as drifters (orange) or escapers (blue). The same distribution is visualised as a cumulative frequency plot (left) and two histograms for each model (right).

Fig. 8. Distributions of the velocities of the ejected components as a function of stellar type. The crossbar marks the median value of each distribution, with the whiskers indicating the 1 interquartile range. The markers show the outliers laying outside the whiskers.

We find that the original tertiary is the component that gets ejected in the majority of simulations. Without the inclusion of SE, the tertiary is ejected in 61% of cases, while this percentage increases to 67% when SE is taken into account (Fig. 9). In the case of pure gravity, we observe about half as many primary ejections versus secondary ejections, while this difference is only 3% when we include SE. When we look at the triples that produce drifters, we find that it is the tertiary that becomes a drifter in 82 − 87% of these case, with the smaller percentage coming from the SE model. We expect that the tertiary is more likely to become a drifter due to the already lower binding energy in the outer orbit.

Fig. 9. Bar chart showing which fraction of escaped stars were the primary, secondary, or tertiary in their original configuration. The solid and patterned bars indicate the distributions for escapers and drifters respectively, with the colours indicating the two different models.

The binaries that are formed following the ejection of one of the components tend to have separations similar to that of the inner semi-major axis at the onset of instability (Fig. 10). This trend is strongest in the case of pure gravity, with the inclusion of SE producing more binaries with semi-major axes larger than the initial inner separation. This result is clearly visible in Fig. 11, which also shows how the new, wider binaries are mainly formed from democratic encounters where one of the stars of the inner binary is ejected. The hierarchical interactions, in which the tertiary is ejected without the democratic flag having been raised, is expected to produce new binaries with similar properties as that of the initially inner. We also find that both models produce binaries with eccentricities closely following a thermal distribution for e ≤ 0.6, while deviating towards a sub-thermal distribution for larger eccentricities.

Fig. 10. Properties of the remaining binary following the ejection of one of the components.

Fig. 11. The semi-major axis of the newly formed binary following an ejection differs from the initial inner semi-major axis. This figure shows the change in separation for the N-body model (top) and the SE model (bottom), with different markers indicating whether the system produced an escaper or drifter, and the colours denoting which component was ejected. The solid line indicates a binary semi-major axis equal to the initial inner semi-major axis.

The value of a_in at the onset of instability is correlated with the type of ejection. In the case of pure gravity, encounters that eject the tertiary without first having a democratic interaction tend to have smaller initial inner separations compared to the other types of ejections (Fig. 12). When SE is considered, we find that the same category of ejection occurs in systems with a substantially higher a_in.

Fig. 12. For the systems that eject one of the components, what are the distributions of the initial inner semi-major axis? This figure shows the cumulative frequency of this quantity for both the N-body model (solid) and the SE model (dashed). The different colours indicate whether the democratic flag was raised during the simulation, and which body was ejected. In the case of hierarchical ejection (blue lines), the tertiary was ejected without the democratic flag being raised, while the other colours indicate that the primary (yellow), secondary (green), or tertiary (black) was ejected after the flag was raised.

The initial mutual inclination also has a significant impact on whether the democratic flag is raised before the ejection of a body. The majority of hierarchical ejections in the case of pure gravity occur in systems with prograde orbits (i ≤ 90°), while non-hierarchical ejections are more evenly distributed over the full range of angles (Fig. 13). When we include SE, the number of systems that eject a body without raising the democratic flag increases. We find that smaller initial inclinations tend to eject the primary (Democratic-1) when SE is included. Retrograde orbits (i ≥ 90°) are generally more stable than prograde orbits (Harrington 1972; Mardling & Aarseth 2001; Vynatheya et al. 2022), but it is possible that the mass loss present in the systems with SE push the triples further into the instability regime, where the effect of prograde versus retrograde orbits becomes less important.

Fig. 13. Similar to Fig. 12 but showing the cumulative distribution of the initial mutual inclination.

In the systems that eject one of the components, the ejected body is the least massive in about 60 − 90% of cases (Fig. 14), with the higher percentage representing ejections of the secondary (Democratic-2). This is similar to results found by TBZ22, who also note that one may naively assume the least massive body will always be the one to get ejected⁶. While this assumption indeed holds true for the systems that eject the secondary, it breaks down when we look at the other ejection channels, wherein the least massive body is ejected in around 65% of cases.

Fig. 14. The ejected star tends to be the least massive. This figure shows a cumulative histogram of the mass of the ejected star divided by the lowest mass in the remaining binary. The different colours indicate the different ejection channels, with the solid lines corresponding to the case of pure gravity, and dashed indicating the inclusion of SE. The vertical, black dashed line indicates a value of 1, so that the left side of the plot corresponds to the ejected body having the lowest mass.

5. Discussion

5.1. Instability rates

As noted by HPTN22, the TEDI is triggered not only by mass loss, but also RLOF in the inner binary. Mass transfer is not included in TRES due to a lack of robust prescriptions for eccentric mass transfer. Therefore, our fractions for triples that become unstable is a lower limit. Including mass transfer in the models can potentially increase the rate of unstable triples by up to a factor 2, as seen in HPTN22. Furthermore, including perturbations caused by close fly-bys of other stars can also push triples over the stability limit, potentially increasing the fraction of TEDI systems by another percent. From Kummer et al. (2023), almost 70% of massive triples undergo mass transfer in the inner binary. If 1% of these lead to an expansion of the inner separation, pushing the system over the stability limit, the fraction of unstable triples is increased by over 40%. As a consequence, we note that the rates calculated for the different outcomes should be considered a lower limit of the real population.

5.2. Galactic rates

To compute galactic rates for the different outcomes, we assume a Galactic star formation rate of 3 M_⊙/yr, a triple star fraction of 0.6 (Offner et al. 2023). The rates for different outcomes are shown in Table 1. We see that the rates are all in the order of 1 event per Myr, with the largest rate originating from collisions in the pure N-body simulation, at 1.3 events/Myr. These rates are dominated by the initial mass function; only about 0.45% of all stars have masses ≥10 M_⊙ (Kroupa 2001). The Galactic rate for the destabilisation of a triple due to stellar winds is estimated to be 3.4 events/Myr. As already mentioned, this is likely a lower limit due to the absence of mass transfer prescriptions in TRES and perturbations from fly-bys. If we assume that these mechanisms could bump up the fraction of systems that become unstable to 5%, as seen in HPTN22, the corresponding Galactic rate would then be 1 event per 10⁵ yr.

5.3. Ejections and runaways

A large fraction of observed massive stars are moving with a high velocity relative to their environment. These stars make up about 5–10% of all observed B-stars, and 10–30% of observed O-stars (Gies 1987; Stone 1991; Oey et al. 2018; Renzo et al. 2019; Stoop et al. 2023), and were dubbed ‘runaway’ stars by Blaauw (1961). These massive runaway stars can have velocities up to 200 km/s, with a dispersion of 30 km/s (Stone 1991). Historically, multiple scenarios have been proposed as mechanisms for producing massive runaways, with two channels being generally regarded as the most dominant; the disruption of a binary due to an internal supernova (Blaauw 1961), and the ejection of a star due to binary-single star interactions in dense stellar clusters (Fujii & Zwart 2011). TBZ22 proposed a new channel for the formation of runaway stars; the destabilisation of a triple system and the subsequent ejection of one of the components. They find that for low-mass stars, the destabilisation of a triple system often leads to the ejection of one of the stellar components, with a velocity distribution of the ejected stars peaking at velocities of a few km/s, and including a tail up to orders of 10² km/s. Additionally, the authors propose estimates of the maximum terminal velocity of ejected massive stars with mass M = 10 M_⊙ and a minimum distance of 10 R_⊙ to be 95 − 222 km/s. Our simulations produce similar results, with a median velocity of 4.25 km/s, and a tail that extends up to a maximum velocity of 116 km/s. With our Galactic model, we estimate rates of 1 − 1.3 ejections per Myr due to the destabilisation of a massive triple. Out of all the ejected bodies, about 8% of them have velocities ≥30 km/s, and around half of these have masses ≥10 M_⊙, which leads to a TEDI-induced massive, high-velocity runaway rate of about 4 − 5.4 × 10⁻⁸ Myr⁻¹. Looking at percentages alone, by assuming a triple rate of 60%, destabilisation rate of 2%, and the results in this work (41% of systems eject a body, 48% of the ejected bodies are stars, 44% of these are massive, and 23% of these are ejected with high velocities), we find that about 0.02% of all massive stars become high-velocity runaways due to triple destabilisation. If we assume that 20% of all massive stars are runaways, then about 0.1% of these originate from triple destabilisation. Even with the assumption that this is a lower limit on the number of massive triples that become unstable, it is likely that the real upper limit is not high enough for this to become a dominant contribution to the population of massive runaways.

Ejected stars tend to escape with velocities on the order of the orbital velocity, and, consequently, triples with more compact inner orbits will produce ejected stars with larger velocities⁷. Combining this with the observations that massive stars in binaries tend to have compact orbits (Sana et al. 2012; Moe & Di Stefano 2017), TBZ22 proposes that destabilised triples may play an important role in answering the questions about the origins of OB runaways. In this study we are able to produce high velocity massive runaways, though we cannot conclude nor dismiss the idea that triples account for a significant fraction of observed OB runaways due to our initial conditions. Our sample of destabilised triples that eject a star are heavily biased towards wider inner orbits (Fig. 10), meaning that we are not investigating the outcome of destabilised triples with more compact orbits, which would be the progenitors of the highest velocity stars. There is a combination of factors that cause our population of unstable triples to have wide inner orbits. As seen in Kummer et al. (2023), the distribution of the inner semi-major axis from which the triple systems are sampled peaks at ∼10³ R_⊙ with an exponentially decaying tail towards higher separations. The outer semi-major axis peaks at 5 × 10⁶ R_⊙, with a similar exponential decaying tail to lower values. When generating a triple system from these distributions, there are four general categories that can be initialised: (1) small a_in and small a_out, (2) large a_in and small a_out, (3) large a_in and large a_out, and (3) small a_in and large a_out. Any of the cases with a small a_in are more likely to experience inner RLOF before any other interaction occurs. Thus, the remaining candidates for producing unstable triples from winds alone are the ones with a large inner a_out. A consequence of the large typical outer semi-major axis means that it is statistically more likely that systems with a wide outer orbit are generated, and as a result, the triples that become unstable from winds must have extremely wide inner orbits if the stellar winds are to expand the orbit enough for instability to occur.

5.3.1. Runaways from RLOF-induced instability

Given that the velocity of an ejected star is strongly anti-correlated with the initial inner separation, we investigate what values of a_in are required to produce stars with high velocities (≥30 km/s). Out of the 25 000 triples simulated by Kummer et al. (2023), 16 833 experienced RLOF in the inner binary as the first interaction. As mentioned in Sect. 5.1, this mechanism might destabilise triples through expansion of a_in. To get an idea how much the inner separation can expand due to mass transfer, we look at two cases: conservative mass transfer and completely non-conservative mass transfer with isotropic reemission. For these two cases, the change in the orbital configuration is given by (Soberman et al. 1997; Pols et al. 1998)

$$\begin{array}{c}\hfill \frac{a_f}{a_i}={\left(\frac{M_{\mathrm{d},\mathrm{i}}}{M_{\mathrm{d},\mathrm{f}}}\frac{M_{\mathrm{a},\mathrm{i}}}{M_{\mathrm{a},\mathrm{f}}}\right)}^2\end{array}$$(3)

and

$$\begin{array}{c}\hfill \frac{a_f}{a_i}=\frac{M_{\mathrm{d},\mathrm{i}}+M_{\mathrm{a}}}{M_{\mathrm{d},\mathrm{f}}+M_{\mathrm{a}}}{\left(\frac{M_{\mathrm{d},\mathrm{i}}}{M_{\mathrm{d},\mathrm{f}}}\right)}^2\exp (-2\frac{M_{\mathrm{d},\mathrm{i}}-M_{\mathrm{d},\mathrm{f}}}{M_{\mathrm{a}}}),\end{array}$$(4)

respectively. Here, subscript i and f denote the initial and final values respectively, while d and a denote the donor star and the companion. We assume the entire envelope is entirely stripped from the donor, which means that in the case of conservative mass transfer, the final mass of the accretor is equal to its initial mass plus the initial envelope mass of the donor. We then check whether a system would undergo stable or unstable mass transfer using the critical mass ratio q_crit from Claeys et al. (2014). Out of the systems that undergo inner mass transfer, 47% of these will be stable. From these systems, 55% of them will increase their separation such that the triple crosses over the instability limit. Finally, 23% of these will destabilise with an inner separation a_in ≤ 10³ R_⊙. We find no significant difference in the results in the case of conservative versus non-conservative mass transfer. If we assume that out of these systems, 40% eject one of the components, then we find that 1.6% of all massive triples can produce high-velocity runaway stars due to destabilisation from inner RLOF, which is 80 times higher than the equivalent rate from stellar winds alone. We emphasise that these results are from back-of-the-envelope calculations, and thus do not encompass the full picture of inner mass transfer in a triple system. For example, we have assumed that the outer orbit does not change during the process. In the case of completely non-conservative mass transfer, the outer orbit is likely to widen due to mass lost from the system. If both the inner and outer orbit widens, the triple might never reach instability, though this depends on the timescales of the expansions, which is beyond the scope of this study. Consequently, these results give an upper limit for the number of system that become unstable due to inner mass transfer.

5.3.2. Ejected BHs

We find that 26–42% of the ejected stellar objects are BHs. These tend to move with velocities similar to that of the ejected stars (Fig. 8), with slightly higher velocities in the SE model due to the addition of a supernova kick. Single BHs can potentially be observed as a microlensing events (Wambsganss 2006). From our Galactic rates, we estimate that the ejection of a BH from an unstable triple to happen at a rate of 0.3 − 0.5 per Myr.

5.4. Collisions

A significant fraction of the destabilising massive triples end up with a collision between two MS stars. These events are thought to be the progenitors of luminous red novae (Kulkarni et al. 2007; Thompson et al. 2009; Tylenda et al. 2011). The rates of massive triples that become unstable and produce a collision correspond to an estimated Galactic rate of 1.1 − 1.3 events/Myr. The majority of these involve the collision of two MS stars, with the second largest population being that of post-MS + MS collisions. Collisions tend to happen in orbits with a small initial inner semi-major axis (Fig. 15). In the case of pure gravity, the 80th percentile of collision systems is around 10⁵ R_⊙, while this is about one order of magnitude lower when we include SE. In other words, including SE tends to favour smaller inner separations when collisions occur, while the pure N-body model is more uniformly distributed over a larger range of a_in. One explanation for this is the increase in radii due to the evolution of the stars; with a larger radii, the systems would not need to be kicked up to extreme eccentricites before collision occur, and would therefore produce collisions at earlier times. This is reflected in Figs. 2a and 2b, with more collisions occurring after fewer outer orbits in the SE models.

Fig. 15. Cumulative histogram of the inner orbital separation at the onset of instability, colour-coded according to the final outcome of the system.

5.4.1. Properties of the newly formed binary

The properties of the newly formed binary following the collision of two components (Fig. 5) are calculated using two main assumptions: (1) perfect inelastic collision without mass loss, and (2) each collision, regardless of the impact parameter, merges the two components. For the first assumption, hydrodynamical simulations have shown that only a small amount of mass is actually lost during a MS-MS merger (Lombardi et al. 2002; Dale & Davies 2006; Glebbeek et al. 2008). As the majority of the collisions in this study are between two MS stars, this assumption is likely to hold for most of the systems. For systems involving a compact object and a star, the star is more likely to be heavily disrupted due to strong tidal forces (Regev & Shara 1987), potentially forming a disc around the compact object before the two components merge. With a potentially higher mass loss during the collision, it is more likely that the newly formed object is not bound to the remaining component, meaning that the collision does not form a new binary, but rather two single objects, one of which is a merger product. The second assumption is less realistic, as contact between two stars may lead to three potential outcomes depending on the impact parameter⁸.

5.4.2. Collisions between two MS stars

A collision between two massive MS stars is likely to lead to a merger, with the merger timescale and final properties of the product depending on several factors, including the components masses and the impact parameter. It is possible that an impact parameter close to unity (grazing encounter) might lead to a contact binary-like system. Several properties of contact binaries are still not well understood, including the effects of tides, the transport of energy between the two stars, and the internal structure (Shu et al. 1980). As the stars in a contact binary continue to evolve and expand, they might eventually overflow the L2 equipotential, resulting in an outflow of mass and loss of angular momentum, causing a rapid coalescence and merger (Pejcha et al. 2016). Merger products from coalescence might have peculiar properties, such as being rapidly rotating (Dale & Davies 2006), having high helium abundance in the envelope (Ivanova & Podsiadlowski 2002), or being highly magnetised (Marchant & Bodensteiner 2024). Merger products from two MS stars can also have the property of being too luminous and blue in comparison to other MS stars in its local environment (Chatterjee et al. 2013 and references within).

Dale & Davies (2006) performed hydrodynamical simulations of collisions and close encounters of massive main-sequence stars. They find that close encounters are likely to lead to a rapid merger, with timescales of a few tens stellar free-fall times. Encounters with an impact parameter close to unity (grazing encounters) were found to produce a common envelope system, which may result in the envelope being expelled and potentially leading to a merger. The most massive stars were found to be more likely to expel the envelope without merging of the two cores, resulting in a compact binary consisting of stripped stars. Additionally, energy deposited into the envelope of the merger product might cause them to expand their radii by up to two orders of magnitude. Such a swell-up of the merger product can precipitate further interactions in the newly formed binary following the merger, potentially resulting in mass transfer and/or another merger, this time between one of the original triple components and the merger product. Collision products are likely to be rapidly rotating, due to the orbital angular momentum of the progenitors being deposited into the spin of the remnant or in the mass that is lost during the merger. If little mass is lost, the product will be rapidly rotating, making them possible progenitors for hypernovae or gamma-ray bursts (Podsiadlowski et al. 2004).

5.4.3. Collisions between an MS and an evolved star

In addition to MS-MS collision, we also produce a high rate of collisions between an evolved star and a MS star. The post-MS stars involved in collisions have a mean mass of 25 M_⊙ with a standard deviation of 8.6 M_⊙. The MS counterparts have a similar spread, but with a mean mass of 12 M_⊙. The outcome of a collision between a MS star and an evolved star depends on factors such as the mass and the exact stellar type of the post-MS components, and the impact parameter. As noted by TBZ22, there are three potential scenarios that may occur: (1) the MS star may pass through the envelope of the evolved star, stripping some of the mass; (2) the two stars merge, and; (3) the two stars form a contact binary through the process of a common envelope scenario, stripping the envelope in the process. In the case that a system such as this does merge, the hydrogen rich core of the MS star rejuvenates the evolved component, effectively extending its lifetime. These remaining hard binaries might become potential progenitors for gravitational wave sources, as the small separation can lead to a rapid coalescence once the stars become compact objects.

5.5. Low-mass vs. high-mass unstable triples

As this study focuses exclusively on triple systems with a high-mass primary, we present here an overview of the most important distinctions between our results and those from TBZ22 and HPTN22, who mainly looked at low-mass stars.

Our population of high-mass triples on the edge of stability consists of 79% MS stars, with a smaller percentage – 6.6% – post-MS stars (Fig. 1), while TBZ22 find that 30 − 40% of low-mass triples contain a giant (post-MS) primary star at the onset of instability. Low mass stars have weaker winds while on the MS compared to high-mass stars, resulting in limited widening of the inner orbit before the primary evolves off the MS. High mass stars are therefore more likely to widen the inner orbit to the point of instability before any of the components become post-MS stars, consequently resulting in a higher fraction of MS stars at the onset of instability.

35 − 40% of our systems end up in a collision between two components (Fig. 3). This is significantly higher than the collision rate for low-mass systems from TBZ22, at 13 − 24% (Fig. 1 in TBZ22), and 14% in HPTN22⁹. In both high-mass and low-mass systems, collisions tend to happen between two MS stars (Fig. 4 and Fig. 1 in TBZ22). The greater collision percentage in the high-mass population is likely a combination of two main factors: the larger radii of high-mass stars, and the initial orbital separation. The inner orbital separation of the high-mass systems are comparable to that of the low mass systems. Having larger radii in the same orbit results in a higher likelihood for collisions to occur due to a greater collisional cross-section.

5.6. Model caveats

As with any model, our model has several shortcomings that should be considered when discussing the results. We know what certain systems experience RLOF, either during the hierarchical phase in an eccentric inner binary, or during a close passage during a democratic interaction. In these simulation we only flag these systems for bookkeeping, but it is likely that many of these occurrences would result in either mass transfer, mass loss, or both. RLOF in eccentric orbits is still poorly understood, and thus it is difficult to state whether or not the inclusion of accurate mass transfer/loss prescriptions would significantly affect the outcomes of the systems where one of the stars overflows its Roche lobe. It is possible that some of these interactions are similar to partial tidal disruption events, in which a star is partially disrupted following a close encounter with a compact object (Wang et al. 2021).

Another obvious caveat is how we include SE. By pre-calculating the evolution of each star individually, we are assuming they evolve in isolation, an assumption that breaks down for the systems that experience RLOF at any point during the simulations, and also for the systems that do not, as some portion of the mass lost from winds would likely be accreted by the companion(s). Mass transfer and wind accretion is likely to change the outcome of some of the triples due to the chaotic nature of systems with strong three-body effects. However, whether the cumulative effect of these changes would significantly affect the overall statistics is unknown.

The choice of initial conditions for the population of stable triples, as described in Sect. 3.1, should also be discussed. In this work, the eccentricities of these systems were drawn from a thermal distribution function. However, as observed by e.g. Moe & Di Stefano (2017), Tokovinin (2020), Hwang et al. (2022) the eccentricity distribution of wide binaries is correlated with the binary semi-major axis, where binaries with separation of ∼10⁴ R_⊙ seem to have uniformly distributed eccentricities, while even wider binaries (a ≈ 10⁵ R_⊙) tend to have eccentricities that follow a super-thermal distribution. Given that the massive triple systems that become unstable due to stellar winds and supernovae kicks tend to be very wide at the onset of instability, one might argue that these systems should be even more eccentric. If this applies to the outer binary, more triples could become unstable within a Hubble time due to the stability of the system being dependent on the outer eccentricity.

Finally, our N-body implementation does not include tidal effects. The inclusion of tides would introduce another energy dissipation term, which could potentially change the outcome of several systems. Given that our sample mostly contains wider binaries, tidal effects would only be significant during close passages in a democratic phase or during the pericentre passage of a highly eccentric triple. It is likely that the inclusion of tidal forces would produce more compact binaries instead of colliders or escapers by shrinking and circularising the inner orbit (Orlov & Petrova 1996 and references within), effectively moving the systems back to the region of dynamical stability. Though further evolution of the system may again widen the inner orbit and trigger another instability. However, the overall effect on the statistics is unclear and beyond the scope of this study.

6. Conclusion

In this work, we studied the final outcomes of massive triple systems that become dynamically unstable due to stellar winds. Using a direct N-body approach coupled with a SE code, we simulated a population of triple systems that reside on the edge of stability. We considered two suites of models: one that only takes into account the gravitational interaction between the bodies, and one suite of models that also include SE. Our results can be summarised as follows:

• Even though all the systems simulated in this work are exactly on the limit of dynamical stability, as defined by Eq. (1), a substantial fraction take up to hundreds and even thousands of outer orbital periods to disintegrate. This is in line with the simulated results from TBZ22. This result not only serves as a justification for including SE in the simulation, but also emphasises the importance of simulating unstable triples for longer timescales before categorising them as remaining hierarchical.

• Between 35% and 40% of the simulated triples end in a collision. The collisions happen overwhelmingly between the primary and secondary components, and occur mostly between two MS stars (between 70% and 78% of collisions). These collisions are likely to merge the two components, producing a star with potentially peculiar properties, such as rapid rotation or high magnetism. After MS–MS collisions, the most common scenario involves collisions between an evolved star and a MS star. The evolved stars in these collisions are mostly very massive, and the collisions tend to be grazing encounters. These are likely to result in a common envelope-like event, which might either strip the envelope and produce a compact binary, or merge the two components. In the case of a compact binary, these might be progenitors for future gravitational wave sources, while a merger between these two components can rejuvenate the evolved star, producing a blue straggler-like star.

• 31–40% of systems eject one of the components, leaving behind a newly formed binary and a single stellar object. The distribution of velocities for these ejected bodies peaks at around 4 − 6 km/s, with a tail towards larger velocities. About 4% of the ejected bodies are massive stars with velocities of ≥30 km/s, putting them firmly in the regime of massive runaways. We estimate a lower limit on massive runaways produced via the unstable triple channel of about 0.02%, which would correspond to about 0.1% of all the massive runaways.

• The velocity of an ejected body is strongly anti-correlated with the inner separation at the time of destabilisation. The population of triples studied in this work have mainly wide inner orbits, and it is therefore still possible that triples might contribute more significantly to high-velocity runaway stars. Other mechanisms for destabilising triples might favour more compact orbits, such as mass transfer in the inner binary, which in turn can produce ejected stars with substantially higher velocities. By looking at extreme cases of completely conservative and non-conservative mass transfer, we use analytical formulae to estimate the final properties of the massive triples that undergo mass transfer in the inner binary. We find that this can mechanism can produce up to ten times as many unstable triples than winds alone, with a significant fraction having compact inner orbits, meaning they may be potential progenitors for high-velocity runaways.

• When SE is considered, about 20% of the systems completely ionise, meaning that all three stars become unbound from each other due to multiple internal supernovae. These cases do not produce a remaining binary, but rather three single stellar objects, at least two of which will be stellar remnants.

• With our Galactic model, we estimate Galactic rates for collisions, ejections, drifters, and ionised triples. All the rates are on the order of 1 Myr⁻¹, which is predominantly due to the low number of massive stars.

• Collisions and ejections both leave behind a binary system. As collisions tend to happen between the components in the inner binary, the properties of the newly formed binaries follow the distributions of the initial outer orbit, and thus produce wide binaries. Similarly, the ejection channel tends to expel the tertiary star, leaving behind the inner binary. Consequently, ejections generally produce tighter binaries with orbital parameters closely resembling those of the inner binary.

• Due to the stronger stellar winds present in massive stars, triples with massive primaries tend to destabilise before evolving off the MS.

• Compared to low-mass triples, systems with a massive primary are up to a factor of 2.7 more likely to end up in a collision. The orbital separations of unstable triples differ little between low-mass and high-mass systems, and so the increased radii of massive stars result in a larger collisional cross-section.

Data availability

The data and code necessary to reproduce the figures in this paper are publicly available on Zenodo: 10.5281/zenodo.13884425

Acknowledgments

ST acknowledges support from the Netherlands Research Council NWO (VIDI 203.061 grant). TB is supported by an appointment to the NASA Postdoctoral Program at the NASA Ames Research Center, administered by Oak Ridge Associated Universities under contract with NASA.

References

All Tables

All Figures

	Fig. 1. Distributions of stellar types for the systems that become dynamically unstable from secular evolution. The purple bars indicate the stellar types at the end of the simulations in which SE is included.
In the text

Fig. 2. Duration of the unstable phase for the most prominent outcomes. The x-axis shows the total simulation time expressed as a multiple of the outer orbital period at the time of destabilisation, denoted N. The dotted line shows the distribution of N for all outcomes, while the dashed line shows the duration in units of crossing times. Panels (a) and (b) show N for the pure gravity and SE simulations, respectively.

In the text

	Fig. 3. Fractions of the different outcomes of the dynamical simulations for both the pure N-body and the SE models. The maximum simulation time was 107 × Pout, init.
In the text

	Fig. 4. Percentages of the stellar types involved in collisions.
In the text

	Fig. 5. Properties of the newly formed binaries following the collision between two bodies in the original triple. The left figure shows the cumulative distribution of the semi-major axis, while the right shows the eccentricity. The dotted and dashed lines represent the semi-major axes and eccentricities at the onset of instability.
In the text

	Fig. 6. Fractions of the different stellar types that are ejected from the system during the unstable phase.
In the text

	Fig. 7. Velocities of the components labelled as drifters (orange) or escapers (blue). The same distribution is visualised as a cumulative frequency plot (left) and two histograms for each model (right).
In the text

	Fig. 8. Distributions of the velocities of the ejected components as a function of stellar type. The crossbar marks the median value of each distribution, with the whiskers indicating the 1 interquartile range. The markers show the outliers laying outside the whiskers.
In the text

	Fig. 9. Bar chart showing which fraction of escaped stars were the primary, secondary, or tertiary in their original configuration. The solid and patterned bars indicate the distributions for escapers and drifters respectively, with the colours indicating the two different models.
In the text

	Fig. 10. Properties of the remaining binary following the ejection of one of the components.
In the text

Fig. 11. The semi-major axis of the newly formed binary following an ejection differs from the initial inner semi-major axis. This figure shows the change in separation for the N-body model (top) and the SE model (bottom), with different markers indicating whether the system produced an escaper or drifter, and the colours denoting which component was ejected. The solid line indicates a binary semi-major axis equal to the initial inner semi-major axis.

In the text

Fig. 12. For the systems that eject one of the components, what are the distributions of the initial inner semi-major axis? This figure shows the cumulative frequency of this quantity for both the N-body model (solid) and the SE model (dashed). The different colours indicate whether the democratic flag was raised during the simulation, and which body was ejected. In the case of hierarchical ejection (blue lines), the tertiary was ejected without the democratic flag being raised, while the other colours indicate that the primary (yellow), secondary (green), or tertiary (black) was ejected after the flag was raised.

In the text

	Fig. 13. Similar to Fig. 12 but showing the cumulative distribution of the initial mutual inclination.
In the text

Fig. 14. The ejected star tends to be the least massive. This figure shows a cumulative histogram of the mass of the ejected star divided by the lowest mass in the remaining binary. The different colours indicate the different ejection channels, with the solid lines corresponding to the case of pure gravity, and dashed indicating the inclusion of SE. The vertical, black dashed line indicates a value of 1, so that the left side of the plot corresponds to the ejected body having the lowest mass.

In the text

	Fig. 15. Cumulative histogram of the inner orbital separation at the onset of instability, colour-coded according to the final outcome of the system.
In the text