© The Authors 2025

1. Introduction

Triple stars are ubiquitous among the observed population of stars. For a detailed overview of stellar multiplicity see Moe & Di Stefano (2017), Offner et al. (2023). Notably, approximately half of solar-type stars are found to possess at least one stellar companion (Tokovinin 2014a), with 10% of all systems existing as triple systems (Tokovinin 2014b). As we move towards higher masses, the companion fraction increases. For young O-type stars, nearly every star is part of a multiple system, with almost 70% comprising a triple or higher order multiple, although this estimate carries some uncertainty (Evans et al. 2005, 2006; Sana et al. 2012, 2014; Kobulnicky et al. 2014; Moe & Di Stefano 2017; Offner et al. 2023). Given these statistics, the influence of stellar multiplicity on the evolution and ultimate fates of stars cannot be understated.

Observed triple stars are found in hierarchical configurations, where a binary is orbited by a more distant tertiary companion, forming a dynamically stable system. In these hierarchical triples, unique dynamical, tidal, and mass transfer interactions emerge through the interplay between the stars, phenomena not encountered in single- or binary-star evolution. For example, angular momentum exchange between the inner binary’s orbit and the tertiary star can induce variations in the inner orbit’s eccentricity and the mutual inclination of the triple over long, secular timescales, a process known as the von Zeipel-Lidov-Kozai (ZLK) mechanism (von Zeipel 1910; Lidov 1962; Kozai 1962; Naoz 2016). Because the effect acts over extended periods, the equations of motion are usually averaged over both the inner and the outer orbit, a method referred to as the secular approximation. The typical secular timescale is on the order of P_out²/P_in ≫ P_out ≫ P_in (Antognini 2015). However, smaller outer-to-inner period ratios could lead to higher eccentricities and merger/collision rates compared to the secular approximation (Grishin et al. 2018; Fragione et al. 2019; Mangipudi et al. 2022). Additionally, tidal deformation of the tertiary star can extract energy from the inner binary (Fuller et al. 2013; Gao et al. 2020). Furthermore, mass transfer from a tertiary proceeds in a fundamentally different manner compared to mass transfer in binary systems (de Vries et al. 2014; Di Stefano 2020; Comerford & Izzard 2020; Glanz & Perets 2021). While conventionally only a small fraction of stars in hierarchical triples undergo tertiary mass transfer (TMT) (Hamers & Thompson 2019; Toonen et al. 2020; Stegmann et al. 2022; Kummer et al. 2023), recent studies by Dorozsmai et al. (2024) have shown that TMT in massive triples occurs in approximately 50% of systems if the stars of the inner binary evolve in a chemically homogeneous manner. The main reason for this high fraction in triples with a chemically homogeneous evolving (CHE) inner binary is the initially compact inner orbits, which permit tight outer orbits, allowing the tertiary to eventually fill its Roche lobe frequently. A comprehensive understanding of the system’s evolution during TMT is therefore necessary to understand the final fate of these complex systems.

Chemically homogeneous evolution, initially presented within the framework of rapidly rotating stars (Maeder 1987), was later proposed to be applicable as well to closely orbiting binaries (de Mink et al. 2009; Song et al. 2016). This mechanism suggests that stellar deformation caused by rapid stellar rotation in tidally locked systems with short orbital periods induces enhanced internal mixing, leading to a homogeneous chemical distribution across the stars. By continually replenishing the core with hydrogen from the envelope during the main sequence (MS), these stars avoid developing a core-envelope boundary, thus preventing significant expansion during their evolution. Consequently, CHE binaries have been proposed as potential progenitors of long gamma-ray bursts (Yoon et al. 2006; Aguilera-Dena et al. 2018), (pulsational) pair-instability supernovae (du Buisson et al. 2020), and gravitational wave (GW) mergers, as they can maintain compact orbits while avoiding mass transfer (de Mink & Mandel 2016; Mandel & de Mink 2016; du Buisson et al. 2020; Hastings et al. 2020; Riley et al. 2021; Vigna-Gómez et al. 2021; Dorozsmai et al. 2024), and they potentially contributed to the reionisation of the Universe (Sibony et al. 2022; Ghodla et al. 2023). Analysis by Tokovinin et al. (2006) indicated that in solar-type systems, short-period binaries are more often accompanied by a tertiary star than wider binaries, highlighting the growing importance of triples in the context of CHE. Despite being a theoretical implication of rapid rotation, observational evidence lends support to the plausibility of CHE, as systems with observed properties consistent with this evolutionary channel have been identified (e.g. Almeida et al. 2015; Shenar et al. 2017; Sharpe et al. 2024). However, Abdul-Masih et al. (2019) did not find clear evidence of efficient mixing due to rotation in the massive overcontact binary VFTS 352 (a promising candidate for CHE evolution), as the observed surface abundances show insufficient helium enrichment.

In this paper, we aim to investigate the impact of stable TMT on the properties of the inner binary, and we explore the population of GW mergers resulting from CHE evolution in triples. To achieve this, we created a rapid analytical model that simulates the evolution of a triple system during stable mass transfer from a tertiary that is applicable to population synthesis simulations of hierarchical triple stars. We applied the model to two population synthesis simulations performed by Dorozsmai et al. (2024) at metallicities of Z = 0.005 and Z = 0.0005.

In Sect. 2, we describe the setup of the TMT model, and this is followed by a discussion on the formation channel of TMT around an inner binary black hole (BBH) and an overview of the initial conditions for the population synthesis simulations in Sect. 3. We then provide a qualitative description of the evolution of an example system during the TMT phase for varying physical assumptions in Sect. 4. In Sect. 5, we present a statistical analysis of the evolution of the simulated populations during the TMT phase and the merger rates and properties of merging black holes (BHs). In Sect. 6, we discuss the caveats and assumptions of our study, and finally, in Sect. 7, we give a summary of our findings.

2. Modelling tertiary mass transfer

In this section we present our analytical framework designed to simulate stable mass transfer originating from a tertiary donor. The analytical nature of our model makes it well suited for integration within population synthesis simulations of triple systems. To provide necessary context, we outline the hierarchical triple star configuration under consideration. Such a system comprises an inner binary, characterised by primary and secondary stars possessing initial zero age MS (ZAMS) masses (m₁, m₂) and angular velocities (Ω₁, Ω₂) assuming rigid-body rotation, and an outer binary consisting of the tertiary star with mass m₃ and angular velocity Ω₃, and the centre of mass of the inner binary. We define the orbital elements of the inner binary as {a_in, e_in, g_in}, where a_in represents the semi-major axis, e_in the eccentricity, and g_in the argument of pericentre. We establish a parallel set for the outer orbit, denoted as {a_out, e_out, g_out}. Lastly, we introduce the mutual inclination between the planes of the inner and outer orbit as i_mut. A sketch of a hierarchical triple system is shown in Fig. 1.

Fig. 1. Illustration of the inner and outer orbits in a hierarchical triple system with masses m1, m2, and m3 along with their orbital elements: semi-major axes, eccentricities, and mutual inclination. The mass streams of a Roche-lobe-filling tertiary that lead to BA (intersecting with the inner orbit) or the formation of a CBD (intersecting with itself) are depicted. We note that the figure is not to scale.

We divide this section in three parts: first, we describe the evolution of the inner binary throughout the TMT phase. Second, we describe the evolution of the outer orbit. Third, we discuss variations to the analytical model. The default model described in this section should be regarded as the most simplistic model, omitting numerous physical processes that are likely to take place. With each model variation, we introduce an extra physical process to the model to assess its significance.

2.1. Evolution of the inner binary

2.1.1. Conditions for disc formation

The trajectory of the mass stream originating from a tertiary star and approaching the inner binary can undergo two distinct evolutionary paths depending upon the stellar and orbital characteristics of the hierarchical triple system. Should the trajectory of the mass stream fail to intersect with the orbit of the inner binary, it is assumed that it will travel around the inner binary, intersect itself, and form a circumbinary disc (CBD) encompassing the binary. Conversely, if the mass stream crosses with the binary orbit, the formation of a CBD is prevented. In this study, we refer to this as ballistic accretion (BA). Figure 1 illustrates the mass flow within a triple system for both scenarios. To determine the minimum distance, R_min, at which the mass stream can approach the accreting binary, we adopted the fitting formulas from Lubow & Shu (1975) and Ulrich & Burger (1976):

$$\begin{array}{c}\hfill R_{\min }=0.425a_{\mathrm{out}}(1-e_{\mathrm{out}})[\frac{1}{q_{\mathrm{out}}}(1+\frac{1}{q_{\mathrm{out}}})],\end{array}$$(1)

where q_out is the mass ratio of the outer binary, defined as q_out = m₃/(m₁ + m₂). If the apocentre of the inner binary, a_apo = a_in(1 + e_in), exceeds R_min, we assume BA occurs, whereas a CBD forms when the apocentre is smaller.

2.1.2. Ballistic accretion

During BA, we assumed that the gas from the stream is redistributed as it intersects with the binary orbit, forming a gas cloud encapsulating the inner binary. This structure bears resemblance to a common envelope phase, as found by de Vries et al. (2014). Recent analytical studies of stable mass transfer in hierarchical triples have treated these configurations as a classic common envelope, parameterised utilising the αλ formalism in order to predict the evolution of the inner binary (Hamers et al. 2022; Gao et al. 2023; Dorozsmai et al. 2024). In contrast, our methodology involved determining the evolution of the inner binary through the calculation of drag forces imparted on the binary components by the surrounding gaseous medium. Gravitational attraction of the wake trailing the objects remove linear momentum and slow down the objects. The amplitude of the force on a perturber mass (m_p) resulting from the gravitational drag force (GDF) was calculated by Ostriker (1999)

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{GDF}}=-\frac{4\pi G^2{m}_p^2{\rho }_g}{v_{\mathrm{rel}}^3}{\mathit{v}}_{\mathrm{rel}}\mathcal{I}(\mathcal{M}),\end{array}$$(2)

with ρ_g the density of the gas surrounding the perturber, v_rel the relative velocity between the gas and the perturber, and ℐ(ℳ) a dimensionless function dependent on the Mach number ℳ = v_rel/c_s, where c_s is the speed of sound in the gas. For lack of a better estimate, we made the assumption that the gas is stationary with respect to the centre of mass (COM) frame. Consequently, the relative velocity is equivalent to the Keplerian velocity of the perturber. The value of ℐ(ℳ) has been determined by Ostriker (1999) for an object that moves on a straight trajectory through a uniform gaseous medium either at subsonic ℳ < 1 or supersonic ℳ ≥ 1 velocities:

$$\begin{array}{c}\hfill I(\mathcal{M})=\{\begin{array}{cc}\frac{1}{2}\ln (\frac{1+\mathcal{M}}{1-\mathcal{M}})-\mathcal{M}\hfill & \text{if}\mathcal{M}<1\hfill \\ \ln (1-\frac{1}{{\mathcal{M}}^2})+\ln (\frac{r_{\max }}{r_{\min }})\hfill & \text{if}\mathcal{M}\ge 1\hfill \end{array},\end{array}$$(3)

where r_max and r_min represent the characteristic sizes encompassed by the gravitational wake. In the case of a binary system, r_max typically approaches the orbital diameter (Kim & Kim 2007; Antoni et al. 2019). Thus we assigned r_max a value of 2a_in. While the exact value of r_min remains uncertain, we adopted the approach by Kim & Kim (2007) and set r_min = a_in/10, which makes ln(r_max/r_min) = ln20 ≈ 3.

The description of the wake formulated in Ostriker (1999) is under the assumption of a single object moving along a linear path through a gaseous medium. However, our situation involves two objects orbiting each other. To address this, Kim & Kim (2007) derived an expression for the gravitational wake of a perturber moving along a circular trajectory. Unlike the symmetrical wake in the former case, the wake has both radial and azimuthal components, with the azimuthal component dominating the orbital decay. To capture this effect, Kim & Kim provided a fitting function for the azimuthal component of the wake:

$$\begin{array}{c}\hfill {\mathcal{I}}_{\varphi }^{\mathrm{K}07}(M)=\{\begin{array}{c}0.7706\ln (\frac{1+\mathcal{M}}{1.0004-0.9185\mathcal{M}})-1.4703\mathcal{M}\hfill \\ \text{if}\mathcal{M}<1.0\hfill \\ [0.5em]\begin{array}{c}\ln (330(a_{\mathrm{in}}/r_{\min })\times {(\mathcal{M}-0.71)}^{5.72}{\mathcal{M}}^{-9.58}).\hfill \\ \text{if}1.0\le \mathcal{M}<4.4\hfill \end{array}\hfill \\ [0.5em]\ln (\frac{a_{\mathrm{in}}/r_{\min }}{0.11\mathcal{M}+1.65})\text{if}4.4\le \mathcal{M}\hfill \end{array}\end{array}$$(4)

In the case of a binary system, the wake generated by one component influences the gravitational drag experienced by the other component. To address this mutual interaction, Kim et al. (2008) introduced an additional term to the gravitational wake expression:

$$\begin{array}{c}\hfill \frac{{\mathcal{I}}_{\varphi }^{\mathrm{K}08}(\mathcal{M})}{{\mathcal{M}}^2}=\{\begin{array}{cc}-0.022(10-\mathcal{M})tanh(3\mathcal{M}/2)\hfill & \text{if}\mathcal{M}<2.97\hfill \\ -0.13+0.07{\tan }^{-1}(5\mathcal{M}-15)\hfill & \text{if}\mathcal{M}\ge 2.97\hfill \end{array}.\end{array}$$(5)

To obtain the total effect of the wake, we summed the contributions from both terms, ${\mathcal{I}}_{\varphi }(\mathcal{M})={\mathcal{I}}_{\varphi }^{\mathrm{K}07}(\mathcal{M})+{\mathcal{I}}_{\varphi }^{\mathrm{K}08}(\mathcal{M})$.

In Appendix B, we derived the change in semi-major axis and eccentricity of the inner orbit over time. For a circular orbit, the eccentricity remains constant, and the semi-major axis decreases due to the dissipation of orbital energy (Grishin & Perets 2016; Rozner & Perets 2022):

$$\begin{array}{c}\hfill \frac{d{a}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{BA}}=-8\pi \sqrt{\frac{G{a}^5}{M}}{\rho }_g\mathcal{I}[\frac{1}{q}{(1+q^{-1})}^2+q{(1+q)}^2],\end{array}$$(6)

where M = m₁ + m₂ is the total mass of the inner binary and q is the mass ratio of the inner binary, defined as q = m₂/m₁.

To simulate the evolution of the inner binary, we had to make assumptions regarding the gas density and local sound speed around the perturbers. However, detailed simulations providing realistic values for these parameters in such systems are currently lacking. In particular, the gas density is anticipated to have a considerable influence on the rate of evolution of the inner binary given that the drag force is directly proportional to ρ_g. To address this uncertainty, we evaluated F_GDF at two distinct gas density values: ρ_g = 10⁻⁸ g cm⁻³ and ρ_g = 10⁻¹⁰ g cm⁻³. We ensured that these densities are less dense than those of typical common envelope scenarios in binaries or triples, which are about 10⁻⁵ − 10⁻⁷ g cm⁻³ (e.g. Ivanova et al. 2013; Glanz & Perets 2021). Moreover, the chosen values are similar to the gas densities found in de Vries et al. (2014). For simplicity, we assumed a constant gas density over time, implying that the mass transfer rate onto the inner binary equals the amount of mass expelled from the binary. Unless otherwise specified, we assumed that none of the mass is accreted by the binary components. Additionally, we adopted a sound speed value of 30 km s⁻¹, corresponding to a gas temperature of approximately 10⁵ K.

2.1.3. Circumbinary discs

The dynamics governing the interactions between a binary system and the surrounding CBD are inherently intricate. For a detailed review on CBDs and their applications, we direct the readers to the comprehensive review by Lai & Muñoz (2023). Various processes can significantly influence the orbital evolution of the binary. On one hand, gravitational torques stemming from resonant interactions between the binary and the disc can extract angular momentum from the binary. Conversely, accretion from the disc onto the binary may transfer angular momentum inward. The balance of these torques dictate how the inner orbit evolves over time. The configuration of the binary, characterised by parameters such as mass ratio (q_in), eccentricity (e_in), and mutual inclination (i_mut), appears to play a pivotal role in determining the binary’s fate (Nixon et al. 2013; Miranda & Lai 2015; Miranda et al. 2017; Muñoz et al. 2019, 2020; Duffell et al. 2020; Zrake et al. 2021; D’Orazio & Duffell 2021; Siwek et al. 2023a,b). Furthermore, alterations in binary properties can lead to consequential shifts in its evolutionary trajectory (Valli et al. 2024).

We utilised findings from Siwek et al. (2023b), who performed hydrodynamical simulations of CBDs around supermassive BH (SMBH) binaries. The results are scale free, making them applicable to stellar-mass BHs as well. Regarding the disc properties, a constant value of 0.1 was assigned to the disc viscosity parameter α. Additionally, the aspect ratio h, representing the disc’s thickness relative to its size, is fixed at 0.1. The study encompasses an extensive parameter exploration using a grid-based approach, wherein the eccentricity, e_in, ranges from 0–0.8, and the mass ratio, q_in, ranges from 0.1–1. For each grid point, time derivatives for the semi-major axis and eccentricity are provided. The rate of change of a_in and e_in is expressed in terms of the mass accretion rate and binary mass as follows:

$$\begin{array}{cc}\hfill \frac{d{a}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{CBD}}& =k_1\frac{{\dot{m}}_{\mathrm{bin}}}{m_{\mathrm{bin}}}{a}_{\mathrm{in}},\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill \frac{d{e}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{CBD}}& =k_2\frac{{\dot{m}}_{\mathrm{bin}}}{m_{\mathrm{bin}}},\hfill \end{array}$$(8)

where k₁ and k₂ are factors dependent on q_in and e_in, and ṁ_bin is always positive as the binary is accreting material. We used linear interpolation to determine the values of k₁ and k₂ for any given q_in and e_in, as illustrated in Fig. 2. A noteworthy observation is the existence of an equilibrium eccentricity around e_in = 0.45 − 0.5 for near-equal mass ratios, likely attributed to a resonant locking mechanism (Artymowicz et al. 1991). For e_in > 0.8, extrapolation is avoided, and instead, the values for e_in = 0.8 were adopted.

Fig. 2. Constants for the time derivatives of the semi-major axis (top) and the eccentricity (bottom) of an inner binary surrounded by a CBD based on the findings of Siwek et al. (2023b). The figure illustrates how certain parts of the parameter space lead to either a decrease or increase of the semi-major axis and the eccentricity. The dashed lines indicate the points where the derivative equals zero.

In scenarios where the orbital configuration permits the formation of a CBD, we assumed that the disc promptly reaches a steady state. Here, the amount of mass accreted onto the inner binary is continuously replenished by the mass transferred from the tertiary star. We adopted the fitting formula proposed by Siwek et al. (2023a) for the preferential accretion onto each binary component, which aligns well with the outcomes of other hydrodynamic simulations of CBDs surrounding SMBHs (Muñoz et al. 2020; Duffell et al. 2020). This relationship is expressed as

$$\begin{array}{c}\hfill \frac{{\dot{m}}_2}{{\dot{m}}_1}=(\frac{m_2}{m_1}{)}^{-0.9},\end{array}$$(9)

resulting in a tendency towards mass equalisation within the binary system.

2.1.4. Gravitational waves

At small separations, angular momentum loss through emission of GWs becomes non-negligible for binary compact objects and can assist the inspiral of the binary towards a merger. The evolution of the semi-major axis and the eccentricity are described by Peters (1964)

$$\begin{array}{cc}\hfill \frac{d{a}_{\mathrm{in}}}{\mathit{dt}}& =\frac{64}{5}\frac{G^3{m}_1{m}_2(m_1+m_2)}{c^5{a}_{\mathrm{in}}^3{(1-e_{\mathrm{in}}^2)}^{7/2}}[1+e_{\mathrm{in}}^2\frac{73}{24}+e_{\mathrm{in}}^4\frac{37}{96}]\hfill \\ \hfill [6pt]\frac{d{e}_{\mathrm{in}}}{\mathit{dt}}& =\frac{304}{15}\frac{G^3{m}_1{m}_2(m_1+m_2)}{c^5{a}_{\mathrm{in}}^4{(1-e_{\mathrm{in}}^2)}^{5/2}}[1+e_{\mathrm{in}}^2\frac{121}{304}].\hfill \end{array}$$(10)

The complete evolution of a_in and e_in will be governed by the summed contribution of the gas drag during BA, CBD dynamics, and GW emission:

$$\begin{array}{cc}\hfill \frac{d{a}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{total}}& =\frac{d{a}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{BA}}+\frac{d{a}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{CBD}}+\frac{d{a}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{GW}}\hfill \\ \hfill [6pt]\frac{d{e}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{total}}& =\frac{d{e}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{CBD}}+\frac{d{e}_{\mathrm{in}}}{\mathit{dt}}{|}_{\mathrm{GW}}.\hfill \end{array}$$(11)

When the TMT phase proceeds via BA, the contribution of the CBD is zero, and vice versa.

2.1.5. Eccentric gas drag

Initially, we considered gas drag effects on objects moving along circular orbits, leading to a constant time derivative of the semi-major axis along the orbit and a constant eccentricity of zero. However, three-body dynamical effects, such as ZLK cycles, can excite the eccentricity of the inner binary prior to the onset of mass transfer. The drag forces become phase dependent, as the velocity of the binary objects relative to the gas is different between pericentre and apocentre. In Appendix B we derived the orbit-averaged time derivatives of the semi-major axis and eccentricity. For binaries embedded in a homogeneous gas cloud, the GDF leads to an increase in eccentricity for supersonic trajectories (e.g. Szölgyén et al. 2022; O’Neill et al. 2024), and a decrease for subsonic ones (e.g. O’Neill et al. 2024). Rozner & Perets (2022) demonstrate that a steep growth of eccentricity towards unity can be achieved, potentially facilitating a merger of the binary.

2.1.6. Isotropic re-emission

So far, we have assumed that the components of the inner binary are able to accrete all material that is transferred from the CBD. However, in reality, the accretion rate is likely to be limited. For stars, accretion efficiency is often tied to the star’s thermal timescale, while for compact objects, it is constrained by their Eddington limit. If the mass transfer rate exceeds the maximum accretion rate, any excess material is expelled from the vicinity of the accreting object. Therefore, we explored a scenario where none of the material is accreted and is instead ejected from the system in the form of a fast isotropic wind or accretion disc wind (e.g. Begelman & McKee 1983; Gallegos-Garcia et al. 2024). The loss of angular momentum due to this expulsion affects the orbital evolution of both the inner and outer orbits, in addition to the evolution resulting from CBD interactions. We assumed the matter is expelled symmetrically with respect to the object and the evolution of the inner binary is similar to that of an isotropic outflow (Hadjidemetriou 1963; Dosopoulou & Kalogera 2016):

$$\begin{array}{cc}\hfill \frac{d{a}_{\mathrm{in}}}{\mathit{dt}}& =-a_{\mathrm{in}}\frac{{\dot{m}}_1^{\text{iso}}+{\dot{m}}_2^{\text{iso}}}{m_1+m_2}\hfill \\ \hfill \frac{d{e}_{\mathrm{in}}}{\mathit{dt}}& =0,\hfill \end{array}$$(12)

where ṁ^iso is the mass lost via an isotropic outflow. It is important to note that a_in increases over time due to the inner binary losing mass, with m₁,m₂ < 0. The consequences of re-emission on the outer orbit are discussed later.

2.1.7. Retrograde torques

Thus far, our model for the CBD interactions assumed that the CBD rotates prograde with respect to the inner binary. However, triples are observed to have mutually inclined inner and outer orbits (Borkovits et al. 2016). If the mutual inclination is larger than 90°, the tertiary star will orbit in retrograde motion with respect to the inner binary, which may result in the formation of a retrograde disc. Nixon et al. (2011) demonstrated that for retrograde discs, orbital resonances are absent. Additionally, the binary components accrete negative angular momentum as the angular momentum vector of the gas particles in a retrograde disc points the opposite direction, leading to a rapid inspiral of the inner binary. Tiede & D’Orazio (2024) performed an extensive parameter study for retrograde discs around SMBH binaries and provided an analytical fit based on their numerical results for the evolution of the semi-major axis and eccentricity:

$$\begin{array}{cc}\hfill \frac{d{a}_{\text{in}}}{\mathit{dt}}& =-10a_{\text{in}}\frac{{\dot{m}}_{\text{bin}}}{m_{\text{bin}}}\hfill \\ \hfill \frac{d{e}_{\text{in}}}{\mathit{dt}}& =30max\{0.1,e_{\text{in}}\}\frac{{\dot{m}}_{\text{bin}}}{m_{\text{bin}}}·\hfill \end{array}$$(13)

2.2. Evolution of the outer binary

The orbital evolution of the outer binary is governed by the exchange and loss of angular momentum within the triple system. Assuming that the tertiary star is sufficiently distant from the inner binary, we can simplify the system by treating the inner binary as a point mass. Under this approximation, the evolution of the semi-major axis of the outer orbit is described by the equation (Soberman et al. 1997; Hamers et al. 2021)

$$\begin{array}{c}\hfill \frac{d{a}_{\mathrm{out}}}{\mathit{dt}}=-2a_{\mathrm{out}}\frac{{\dot{m}}_3}{m_3}[1-\beta \frac{m_3}{m_1+m_2}-(1-\beta )(\gamma +\frac{1}{2})\frac{m_3}{m_1+m_2+m_3}],\end{array}$$(14)

which is analogous to mass transfer in isolated binaries. Here, β is the mass accretion efficiency onto the accretor (i.e. the inner binary), where β = 0 and β = 1 respectively represent completely non-conservative and conservative mass transfer. The parameter γ represents the specific angular momentum of the material that is expelled from the system. For simplicity, we assumed that the eccentricity remains unchanged during the mass transfer.

In our fiducial model, we assumed that all transferred material during BA is eventually expelled isotropically. This material carries a specific angular momentum equal to that of the inner binary. Thus, we set β = 0 and γ = m₃/(m₁ + m₂). If a CBD forms, we assumed that all transferred mass is accreted onto the inner binary. In this scenario, β = 1 (and γ remains unspecified).

In scenarios involving isotropic re-emission within the CBD, the mass transfer is considered entirely non-conservative such that β = 0 and γ = m₃/(m₁ + m₂). Concerning the evolution of the outer orbit, we assumed that all material is ejected at the centre of mass of the inner binary and therefore carries the specific angular momentum of the inner binary. This assumption is justified by the substantial difference in size between the outer and inner orbits.

2.3. Model variations

Here we describe the various models that we investigate later in the paper. A summary of all models is presented in Table 1. First, we start with the most simplistic scenario and only consider BA and GW emission. We call this the Model Simple. The GDF exerted on the binary components is based on the linear perturbation theory (Eqs. (3) + (6)). Estimating these forces requires making assumptions for the gas density and local sound speed around the binary components. To evaluate the implications of the uncertainty in the gas properties, we ran 4 models with varying gas densities and mass transfer rates (defined in Sect. 3.3), comparing the resulting final inner orbits across the models. Additionally, we included the mass-transfer rate as a parameter, as it is expected to influence the gas density around the inner binary. Thus, we explored all four permutations of two gas densities, ρ_gas = 10⁻⁸ g cm⁻³ and ρ_gas = 10⁻¹⁰ g cm⁻³, and two mass-transfer rates, namely the mass transfer that occurs on the thermal and the nuclear timescale of the donor star.

For the next model, we introduced the potential formation of a CBD around the inner binary. The evolution of the semi-major axis and eccentricity in the presence of a CBD are described by Eqs. (7) and (8). We assumed the mass transfer onto the components of the inner binary is conservative. We refer to this as Model Basic. Similar to the previous model, we explored the same four combinations of gas densities and mass transfer rates. However, in Sect. 5.1.1 we demonstrate that in the presence of a CBD, the final distribution of semi-major axes in a population tends to converge with increasing gas density. Therefore, for subsequent models, we focused solely on the combination of ρ_gas = 10⁻⁸ g cm⁻³ with nuclear timescale mass transfer, and ρ_gas = 10⁻¹⁰ g cm⁻³ with thermal timescale mass transfer, representing the two extremes.

The following four models extend Model Basic by incorporating additional physical phenomena, with each model focusing on one aspect at a time to assess its significance. In Model Basic + BinGDF, we integrated the non-symmetric GDF experienced by binary components during BA, as described by Eqs. (4) and (5), originating from the wake generated by two binary components. In Model Basic + EccGDF, we introduced a modification to the GDF under the assumption of eccentric orbits during BA. The evolution of semi-major axis and eccentricity follows Eqs. (B.6) and (B.9). In Model Basic + Iso, we assumed that during the CBD phase, none of the transferred mass accretes onto the inner binary components, instead being isotropically (or symmetrically) ejected near the accreting object. An additional component of the orbital evolution due to the angular momentum loss from the inner binary is described by Eq. (12). Model Basic + Retro incorporates an adjusted formalism for the evolution of semi-major axis and eccentricity during the CBD phase for systems with a retrograde disc, described by Eq. (13). Finally, we merged the Basic + models, incorporating all four additional physical processes described in Sect. 2.1. We refer to this as Model Advanced. We emphasise that Model Advanced is not necessarily the most realistic model, as the physical assumptions are subject to uncertainties.

We introduced additional variations to the Basic and Advanced models, referred to as Model Basic + NoDrag and Model Advanced + NoDrag, where the effects of gas drag are completely ignored. If the total mass of the gaseous cloud enveloping the inner binary is significantly lower than that of the binary components, the formation of an overdense wake trailing the objects becomes inefficient, leading to a minimal inspiral of the binary. In this scenario, we specifically examined the case where the inspiral process is entirely ineffective.

3. Triple population synthesis

3.1. Preceding secular evolution

We first give an overview of the typical evolution of a triple system with a CHE inner binary that leads to TMT and highlight the important intermediate phases. For the system to potentially produce a GW merger, we implement the condition that the components of the inner binary have already evolved to their remnant stage prior to the onset of mass transfer. A cartoon of the evolution is provided in Fig. 3.

Fig. 3. Evolutionary stages of triples with a CHE inner binary that result in TMT with an inner BBH. Depending on the properties of the triple at the onset of TMT, mass transfer occurs via BA or the formation of a CBD. The masses of the three companions are indicated at the ZAMS (MZAMS), MS (MMS), and helium MS (MHeMS) phases. Furthermore, the BH mass (MBH), pre-tertiary mass-transfer mass (Mpre − TMT), envelope mass (Menv), and core mass (Mcore) are indicated. The grey and red crosses denote the centre of mass of the inner binary and outer binary, respectively. The figure is not to scale and has been adapted from van Son et al. (2022).

(I) At the ZAMS, the stars of the inner binary are on a sufficiently short orbit such that tidal dissipation is strong enough to have synchronised the rotational period of the stars with the orbital period, commonly referred to as tidal locking. Consequently, rotation-induced mixing can prevent the star from establishing a chemical gradient (Maeder 1987; Maeder & Meynet 2000; Heger et al. 2000). As opposed to classical evolving stars, these CHE stars do not increase their radius during the MS. On the contrary, by preventing the formation of a chemical gradient within the star and enriching the envelope with helium, the star is expected to contract (Yoon et al. 2006; Mandel & de Mink 2016). Due to the short orbital period, the binary is likely to be in over-contact. During this phase, one of the stars must prevent its second Lagrangian (L2) point from overflowing, as this would lead to a loss of angular momentum through mass outflow, followed by a merger of the binary (Soberman et al. 1997; Marchant et al. 2016). During the MS phase of the inner binary components, any three-body dynamical effects are effectively quenched by strong tidal forces, preventing (e.g. an eccentric merger from occurring). Initially, the tertiary is less massive than either component of the inner binary, ensuring the star will not significantly expand and fill its Roche lobe before the other stars have formed a BH.

(II) Both stars of the inner binary have completed hydrogen fusion. As the interior of the stars has been efficiently mixed during the MS phase, the entire hydrogen content has been fused into helium. The radius of such a helium star is smaller compared to the radius at the ZAMS, and hence the stars detach and stop evolving in a chemically homogeneous manner, with the continued evolution being similar to that of a helium star. As tidal forces are not as strong as before, dynamical interactions between the tertiary and the inner binary can become more pronounced.

(III) The stars of the inner binary undergo core collapse and a BBH is formed. The supernovae might be accompanied by a kick, mainly leading to an increase in eccentricity of both the inner and the outer orbit, but not strong enough to unbind the system. Eventually, the tertiary will evolve off the MS and expand until it fills its Roche lobe, typically when evolving through the Hertzsprung Gap, or after helium has been ignited in the core. We stress that the radial evolution of the most massive stars remains poorly understood, and constitutes a major uncertainty in determining whether the tertiary star can fill its Roche lobe.

(IV) The tertiary star commences mass transfer towards the inner binary. The TMT can proceed in two distinct ways: (i) if the inner orbit is wide and/or the outer orbit is compact, the mass stream intersects with the orbit, resulting in BA. (ii) if the inner orbit is compact and/or the outer orbit is wide, the mass stream misses the orbit and intersects with itself, forming a CBD.

3.2. Initial model assumptions

We applied the TMT model described in Sect. 2 on a large number of hierarchical triple systems with a Roche-lobe filling tertiary star. The prior evolution of each system was simulated with the triple population synthesis code TRES¹ (Toonen et al. 2016). In total, we evolve 6 × 10⁵ systems from the ZAMS, spread over metallicities of Z = 0.005 (0.25 Z_⊙) and Z = 0.0005 (0.025 Z_⊙), identical to the set-up described in Dorozsmai et al. (2024). We only selected systems that follow the evolutionary steps described in the previous section. In the remainder of this section, we provide a description of the initial set-up of the simulations and the conditions that ensure CHE in the inner binary.

The ZAMS properties of the stars and orbits are based on surveys of massive binaries (e.g. Sana et al. 2012; Kobulnicky et al. 2014) for the inner binaries and observations of solar-type stars in hierarchical triples (Tokovinin et al. 2006; Tokovinin 2014b) for the outer binaries. The initial masses of the primary star, m₁, follow a power-law distribution with index α = −2.3 (Kroupa 2001). The initial mass ratios between of the inner orbit, q_in = m₂/m₁, and the outer orbit, q_out = m₃/(m₁ + m₂), both follow a flat distribution. The initial inner and outer semi-major axes follow a flat distribution in log-space. The inner orbits are initially circular and the stellar spins are synchronised with the orbital period, consistent with being tidally locked. The outer eccentricities were drawn from a thermal distribution. The arguments of pericentre of the inner and outer binary were uniformly sampled between −π and π. Lastly, the relative inclination is sampled uniformly from a cosine distribution between 0 and π.

To ensure that the stars of the inner binary are tidally locked and can rotate rapidly enough to sustain efficient mixing within the star, we sampled the initial system properties from a restricted range of values. The initial mass of the primary stars are within the range [20,100] M_⊙ and the semi-major axis of the inner orbits are within the range [14,40] R_⊙. Additionally, following Dorozsmai et al. (2024), to prevent an early stellar merger, the mass ratios of the inner binary were sampled between 0.7 and 1. The ranges for the outer binary are more relaxed: the tertiary masses are in the range [5,100] M_⊙, mass ratios in the range [0.1,1] and the semi-major axes in the range [100,10⁵] R_⊙.

We determined for each sampled system whether CHE takes place in the stars of the inner binary by comparing their rotational angular velocity to a critical angular velocity that is necessary for efficient internal mixing. We used fits produced by Riley et al. (2021) that are based on numerical simulations performed by Marchant et al. (2016) that provide the critical angular velocities for the given set of parameters of the inner binary. Any binary whose stellar components do not have an angular velocity exceeding this critical limit at the ZAMS, or at any later stage of the MS evolution, were discarded. During the CHE, the radius remains constant, apart from the effects of mass loss due to stellar winds. Once a CHE star leaves the MS, it is assumed to directly transition to a helium star. Following the Hurley tracks for MS and helium stars (Hurley et al. 2000), this coincides with a sudden drop in radius of the star. However, detailed simulations of chemically homogeneously evolving systems predict a more gradual decrease of the radius over the MS (see Maeder 1987).

As the initial separation of the inner binary is small, the binary is likely to be in contact. During the contact phase, we assumed that the masses equalise instantaneously (see Riley et al. 2021); although observations suggest that contact binaries are preferably in non-equal mass configurations (Menon et al. 2021; Abdul-Masih et al. 2022), which could be an indication that theoretical models of contact binaries are missing some important physics (Abdul-Masih et al. 2021; Fabry et al. 2023). Due to the redistribution of angular momentum within the binary, the semi-major axis evolves as

$$\begin{array}{c}\hfill \frac{a_f}{a_i}=(\frac{m_{1,i}{m}_{2,i}}{m_{1,f}{m}_{2,f}}{)}^2,\end{array}$$(15)

where the subscripts i and f refer to the initial and final state of the mass equalisation. If either star in the binary overflowed its L2 point, the system was discarded. The location of the L2 point was modelled following Marchant et al. (2016):

$$\begin{array}{c}\hfill \frac{R_{L2,2}-R_{RL,2}}{R_{RL,2}}=0.299{tan}^{-1}(1.84q^{0.397}),\end{array}$$(16)

where R_RL, 2 is the Roche lobe of the secondary star.

3.3. Stability of mass transfer

The stability of the TMT is an important ingredient for the outcome of the mass-transfer phase. CE evolution in triples has been associated with dynamical ejections through destabilisation of the triple (e.g. Comerford & Izzard 2020); the plunge-in timescale of the inner binary towards the core of the tertiary star is short compared to the orbital evolution of the inner binary. Although, other outcomes such as a merger of the inner binary or a merger between the evolved tertiary star and the inner binary are also possible (Glanz & Perets 2021). To that end, we only considered systems where the TMT phase is dynamically stable. Dorozsmai et al. (2024) predicts the vast majority of TMT phases to be dynamically stable, as the mass ratios between the tertiary star and the inner binary (q_out) are typically smaller than 1.

The stability of the mass transfer phase is determined by the evolution of the Roche lobe and the radius of the donor star in response to mass loss. As the star tries to regain hydrostatic and thermal equilibrium, the radius might shrink or expand, depending on the structure of the envelope (Hjellming & Webbink 1987; Soberman et al. 1997). We determined the logarithmic derivatives of the radius with respect to the mass on the dynamical timescale:

$$\begin{array}{c}\hfill {\zeta }_{\text{ad}}={\left(\frac{\text{d}lnR}{\text{d}lnm}\right)}_{\text{ad}}.\end{array}$$(17)

The value of ζ_ad depends on the structure of the envelope and thus on the evolutionary type of the donor star. The values used for each evolutionary type are similar to those described in the binary population synthesis code SeBa (Toonen et al. 2012). The evolution of the Roche lobe during the mass transfer phase is mainly governed by the redistribution of mass within the binary and eventual angular momentum loss. We calculated ζ_RL:

$$\begin{array}{c}\hfill {\zeta }_{\text{RL}}=\left(\frac{\text{d}ln{R}_{\text{L}}}{\text{d}lnm}\right).\end{array}$$(18)

The mass transfer is dynamically stable when ζ_RL ≤ ζ_ad; that is, the star can restore hydrostatic equilibrium.

For dynamically stable mass transfer, the rate of mass loss from the donor star depends on the thermal response of its envelope. If the star is able to regain thermal equilibrium, the mass transfer occurs on the nuclear timescale of the donor star. The corresponding mass transfer rate can be estimated as follows:

$$\begin{array}{c}\hfill {\dot{m}}_{\mathrm{nuc}}=\frac{m}{{\tau }_{\mathrm{nuc}}}\approx \frac{{10}^{-10}}{\varphi }\frac{L}{{\mathrm{L}}_{\odot }}{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(19)

where ϕ is a factor depending on the type of nuclear burning. For core-hydrogen burning, ϕ = 1. For core-helium burning, we assumed ϕ = 0.1. If the star does not regain thermal equilibrium, the mass transfer occurs on its thermal timescale. The corresponding mass transfer rate can be estimated as follows:

$$\begin{array}{c}\hfill {\dot{m}}_{\mathrm{th}}=\frac{m}{{\tau }_{\mathrm{th}}}\approx \frac{2RL}{\text{G}m}\approx 6.67\times {10}^{-8}\left(\frac{{\mathrm{M}}_{\odot }}{m}\right)\left(\frac{R}{{\mathrm{R}}_{\odot }}\right)\left(\frac{L}{{\mathrm{L}}_{\odot }}\right){\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(20)

with τ_th as the thermal timescale and L as the luminosity of the donor star.

In our simulations, if the mass transfer is dynamically stable, we either assumed mass transfer on the nuclear timescale or the thermal timescale for all systems. This way, we can more effectively capture the full range of uncertainties of the mass transfer on the formation of GW mergers.

4. Example TMT model

We applied our model to a representative system drawn from the population outlined in Sect. 3, mainly focusing on the evolution of the inner binary. To highlight the effects from the different TMT models, we qualitatively compare the differences between all the models described in Sect. 2.3. At the onset of TMT, the inner orbit comprised two closely orbiting BHs with orbital properties a_in = 25 R_⊙ and e_in = 0.2. The outer orbit had a_out = 200 R_⊙ and e_out = 0.3. The masses of the components are m₁ = m₂ = 40 M_⊙ and m₃ = 25 M_⊙, with a tertiary envelope mass of 16.7 M_⊙. We chose a gas density of ρ = 10⁻⁸ g cm⁻³ and a mass transfer rate of ṁ₃ = ṁ_nuc. To determine the nuclear timescale of the tertiary, we assigned the initial ZAMS mass and luminosity of the tertiary to m_3i = 27 M_⊙ and 8.5 × 10⁴ L_⊙, giving a timescale of 2.9 × 10⁵ yr (Eq. (19)). The triple has a mutual inclination i_mut = 115°, which implies a retrograde orbit of the tertiary around the inner binary.

In Fig. 4, we show the semi-major axis and eccentricity evolution of the inner binary during the TMT phase for all setups:

Fig. 4. Evolution of the semi-major axis and eccentricity of the inner binary of an example system during the TMT phase and the subsequent GW inspiral. Different models are applied to the same system, as described in Sect. 2.3. The red shaded region indicates that the system undergoes BA. The green shaded region indicates that the system has formed a CBD. The yellow region implies that the mass transfer is non-conservative during the CBD phase, and the material is assumed to be lost via an isotropic outflow. The blue shaded region indicates that the inspiral of the inner binary is dominated by GW emission. Overlapping colours show that more than one of the aforementioned processes apply.

• Model Simple: We have assumed that the complete evolution is driven by BA, precluding CBD formation. The gas surrounding the binary exerts drag forces on the BHs, prompting a merger within 640 years. After 500 years, the GW inspiral starts to dominate over the BA, resulting in a more rapid decay and circularisation of the orbit and an eventual merger.

• Model Basic: The evolution of the system shows significant differences with Model Simple. For roughly the first 10 years, the inner orbit is too wide for a CBD to form with the current orbital configuration, and the mass transfer proceeds in a ballistic fashion. During this period, the inner binary shrinks rapidly from 25 R_⊙ to 10 R_⊙. Afterwards, a CBD is formed and the evolution is directed by angular momentum exchange between the binary and the disc. As a consequence, during the next 1.7 × 10⁵ years, the eccentricity increases quickly to the equilibrium value of about 0.48, while the orbit shrinks modestly, down to 7.5 R_⊙. For the following 1.8 × 10⁵ years, angular momentum loss through GW emission becomes dominant, complementing the CBD torques. As a result, the eccentricity decreases, while the binary orbit shrinks more rapidly. The binary is not able to merge before the tertiary has transferred its entire envelope. At the termination of the TMT phase, the inner binary has a separation of 6.2 R_⊙. Subsequent GW emission brings the system to merger within another 4 × 10⁵ years. Evidently, inspiral via BA is significantly more efficient than via CBD accretion for sufficiently high gas densities. It is noteworthy that a system that starts out in the regime of BA is always expected to transition to a CBD if initially q_out < 1. The reason is that during BA, the outer orbit widens for the assumption of isotropic re-emission, while the inner orbit shrinks, facilitating the formation of a CBD.

• Model Basic + BinGDF: The qualitative evolution of the inner orbit is similar to that of Model Basic. The only notable difference is that the transition from BA to CBD formation occurs later (after 40 years instead of 10 years). This is because the azimuthal drag force exerted by the wake of the companion is negative (it tends to speed up the object instead of slowing it down) and hence slowing the inspiral (Kim et al. 2008). The final semi-major axis and eccentricity are close to identical, as the transition to the CBD formation takes place at a similar inner separation.

• Model Basic + EccGDF: Whereas Model Basic + BinGDF shows results similar to those of Model Basic, this model significantly deviates. During BA, the eccentricity of the inner binary increases drastically, reaching 0.99 after only 12 years. Because the apocentre of the inner binary increases at a faster rate than the semi-major axis decreases, the transition from BA to the formation of a CBD never occurs. Eventually, the eccentricity is sufficiently large, and the semi-major axis sufficiently small, for GW emission to dominate the evolution of the binary, leading to a rapid merger after 13.5 years.

• Model Basic + Iso: The evolution during the CBD phase is altered with respect to Model Basic. Since the conservativeness of mass transfer is only changed in the presence of a CBD, the evolution remains unchanged during BA. As the mass carried away has a specific angular momentum smaller than that of the binary orbit, the orbit effectively widens. The competing effects of the CBD torques, GW emission, and mass outflow leads to a slower decrease of the inner semi-major axis compared to Model Basic. As a result, the separation at the termination of TMT is 8.6 R_⊙ (versus 7.5 R_⊙ for Model Basic).

• Model Basic + Retro: Similar to Model Basic + EccGDF, this model significantly deviates from Model Basic. During the CBD phase, the retrograde torques are more efficient at shrinking the orbit than for the prograde scenario, and the binary merges within 1.7 × 10⁵ years. During the merger, the tertiary is still transferring mass towards the inner binary. During the final phase of evolution, the orbit is almost completely circularised by the GW emission.

• Model Advanced: All features of the Basic + models are combined. However, the evolution of the inner binary is very similar to that of Model Basic + EccGDF. This is because during BA the eccentricity rapidly increases, transitioning to GW dominated inspiral before a CBD can be formed. Therefore, the evolution of the inner binary only depends on the assumptions of the gas drag.

• Model Basic + NoDrag: The inner binary remains unchanged during BA, as no drag forces are acting on its components. After approximately 10⁵ years, the system transitions to a CBD phase. This shift occurs due to the outer orbit’s expansion, caused by the redistribution and/or loss of angular momentum during the mass transfer process. The orbit shrinks to 20.2 R_⊙ once the mass transfer ceases, which is considerably wider than in models where drag forces are included. Nevertheless, GW emission leads to a merger within 60 Myr.

• Model Advanced + NoDrag: The evolution follows a similar pattern to that of Model Basic + NoDrag, but during the CBD phase, retrograde torques cause the orbit to contract more efficiently to 10.2 R_⊙. As a result, the merger occurs significantly faster than in Model Basic + NoDrag, within 8.5 Myr. However, this is still at least a factor 20 longer than for models including GDF.

5. Application to a triple population

Across both metallicities, we find that for about 7.5% of triples in the simulated populations, the inner binary evolves chemically homogeneous for the entire duration of the MS phase. For this subset of systems, the initial masses of the primary stars are at least 30 M_⊙. Among this population, 9.4% of systems at Z = 0.005 and 9.5% at Z = 0.0005 evolve towards TMT with an inner BBH. For these systems, we investigate the effect of TMT on their evolution and their prospects as GW sources.

5.1. Effect of TMT on inner orbits

In this section we examine the impact of TMT on the inner orbits consisting of a BBH, with a primary focus on the semi-major axis and eccentricity. We first discuss the population of high metallicity stars (Z = 0.005), followed by a separate analysis of the role of metallicity.

5.1.1. Inner semi-major axis

Only BA (Model Simple). Under the assumption that a CBD does not form, the semi-major axis at the end of TMT is mainly determined by the mass-transfer rate of the tertiary and the gas density of the cloud encompassing the inner binary during BA. A combination of high gas density and a long mass-transfer timescale results in very efficient inspiral of the binary, whereas low gas density and a short mass-transfer timescale lead to less efficient inspiral. This is evident from the final separations of the inner binary at the end of TMT: the average semi-major axis reduces from 49.0 R_⊙ to 34.4 R_⊙ under conditions of ρ_g = 10⁻¹⁰ g cm⁻³ with ṁ₃ = ṁ_th. Conversely, all the inner orbits merge when ρ_g = 10⁻⁸ g cm⁻³ with ṁ₃ = ṁ_nuc.

BA and CBD (Models Basic(+)/Advanced). If a CBD is allowed to form during TMT, the differences in the final semi-major axis between different assumptions for the gas density and mass transfer rate become less pronounced. In Fig. 5, we show the inner semi-major axes after TMT for the four combinations of the mass transfer rate and gas density described in Sect. 2.3, applied to Model Basic. The average semi-major axes range from 26.7 R_⊙ to 34 R_⊙, indicating more consistent semi-major axes compared to scenarios considering only BA. Moreover, the final separations for all the combinations of gas density and mass-transfer rate, except for ρ_g = 10⁻¹⁰ g cm⁻³ with ṁ₃ = ṁ_th, are nearly identical. This implies that the final semi-major axis of the inner binary is not significantly dependent on the gas properties during BA, especially towards high gas densities. Moving forward, we focus on two specific combinations of gas density and mass transfer rate for the other models that allow for the formation of a CBD: [ρ_g = 10⁻¹⁰ g cm⁻³ with ṁ₃ = ṁ_th] and [ρ_g = 10⁻⁸ g cm⁻³ with ṁ₃ = ṁ_nuc].

Fig. 5. Cumulative distribution of the semi-major axis of the inner binaries at the end of TMT for Model Basic at Z = 0.005. The few systems that have merged during TMT are omitted. The black dashed line denotes the semi-major axis distribution at the onset of mass transfer. Among different assumption for the gas density during BA and the mass transfer rate, the final distributions are fairly similar.

Our findings indicate that the inclusion of eccentric drag forces during BA, the orientation of the CBD, and the conservativeness of mass transfer all have a significant influence on the evolution of the inner orbit. Systems that initially do not form a CBD can attain eccentricities close to unity, leading to an efficient merger through GW emission. However, if the orbit shrinks more rapidly than the eccentricity increases, a CBD might be formed before a high eccentricity is reached. Among the entire population, up to 20% of the inner binaries are expected to merge during TMT due to the inclusion of eccentric drag forces. For systems that do not have a merger during TMT, the final inner semi-major axes are similar to Model Basic. Similarly, if we assume that the triples with a tertiary in a retrograde orbit relative to the inner binary can also form a retrograde CBD, the inner orbits at the end of TMT are more compact. On average, the inner semi-major axis is 12−13 R_⊙ smaller compared to scenarios considering only prograde discs. Isotropic mass outflow from the vicinity of the BHs on the other hand, leads to widening of the orbit through loss of angular momentum. If the TMT is completely non-conservative, the typical semi-major axis at the end of TMT is approximately 10 R_⊙ wider compared to completely conservative mass transfer.

No drag during BA (Models NoDrag). Neglecting drag forces entirely during BA leads to significantly wider final separations of the inner binaries. For systems initiating with BA, the average semi-major axes can be up to four times larger by the end of TMT compared to the models that include drag forces. However, in most systems, a CBD is expected to form immediately at the onset of the mass transfer, and its subsequent evolution remains consistent with models that account for drag forces.

5.1.2. Eccentricity

The evolution of the eccentricity of the inner binary is mainly governed by the interplay between the CBD and GW emission. For prograde discs, the eccentricity evolves towards an equilibrium point (until GW effects take over), as illustrated by the example system in Fig. 4. Consequently, for prograde discs, the typical inner eccentricities increase from 0.22 to 0.47–0.48 at the end of TMT, aligning with the equilibrium eccentricity for near-equal masses. We note that the semi-major axis consistently decreases once this equilibrium eccentricity is reached. This trend holds across the range of mass ratios considered in this study. However, it may not apply for significantly unequal mass ratios (see Fig. 2). In contrast, systems influenced by retrograde disc torques do not reach an equilibrium eccentricity; instead, the eccentricity continually increases towards unity. This high eccentricity leads to rapid loss of orbital energy via GW emission, resulting in fast inspiral of the binary. As GWs also carry away angular momentum, the system does not maintain high eccentricities for long and eventually becomes nearly circularised before merging.

However, if we include eccentric drag forces, the evolution of the eccentricity is most prominent during BA and the GW dominated regime. As demonstrated in the previous section, the eccentricity can increase towards unity during BA for many systems, leading to an efficient merger

5.1.3. Metallicity

In this section we discuss the role of metallicity on the final semi-major axis and eccentricity of the inner binaries. While the evolution of the system during TMT is not directly affected by metallicity in our model, the properties at the onset of TMT can be vastly different.

Firstly, stellar winds are metallicity dependent (Vink et al. 2001; Vink & de Koter 2005) and effectively widen the inner orbit before the stars evolve into BHs, mainly during the post-MS phase of the BH progenitors. Therefore, high-metallicity systems are expected to develop wider orbits, whereas low-metallicity systems tend to remain more compact. The quantified effect of winds on the inner semi-major axes at the onset of TMT is shown in the left panels of Fig. 6. For systems with a metallicity of Z = 0.005, the orbits range from approximately 20–200 R_⊙, with an average of 49.0 R_⊙. Notably, the smallest inner semi-major axes are comparable to the MS radii of the BH progenitors, while the widest orbits significantly exceed the initial upper sampling limit of 40 R_⊙. In contrast, systems with a lower metallicity of Z = 0.0005 experience considerably weaker stellar winds, resulting in a narrower range of semi-major axes at the onset of TMT. The widest orbits are approximately 30 R_⊙, with an average of 20.3 R_⊙. Some BBH systems have semi-major axes smaller than 15 R_⊙, which is less than the combined radii of their MS progenitors. These systems have been subjected to strong three-body dynamical interactions, coupled with either tidal forces or GW emission, depending on the stellar type of the inner binary stars, resulting in energy dissipation from the binary.

Fig. 6. Distribution of the semi-major axis and eccentricity of the inner binaries with a metallicity of Z = 0.005 (top) and Z = 0.0005 (bottom) at the onset of TMT. The purple and red shaded regions correspond to systems that start the TMT phase with the formation of a CBD and BA, respectively. The red bars are plotted on top of the purple bars. We note that the range of the semi-major axis is different between the upper and lower panel.

Secondly, the wider inner orbits of high-metallicity stars also lead to more pronounced three-body dynamical interactions during the post-MS phase of the inner binary components compared to the lower metallicity systems. This is because the timescale of the ZLK cycles scales with the ratio between the outer and inner orbital period, and the separations of the outer orbits that undergo TMT are comparable across both metallicities considered in this study. Therefore, systems with metallicities of Z = 0.005 have higher eccentricities at the onset of TMT than systems with metallicities of Z = 0.0005, as shown in the right panels of Fig. 6. In both cases, the majority of systems have a non-zero eccentricity, with a prominent peak just above 0.1, attributed to the Blaauw kick experienced by BHs in our models due to neutrino losses. However, we predict a tail extending from eccentricities above 0.15 up to 1 only for the high-metallicity systems.

Thirdly, the properties of the inner and outer orbits at the onset of TMT play a crucial role in determining the formation of a CBD. With wider and more eccentric inner orbits at higher metallicities, the likelihood of disc formation decreases. Initially, 57.1% of systems with a metallicity of Z = 0.005 form a CBD, whereas for systems with a metallicity of Z = 0.0005, the percentage increases to 74.4%.

Lastly, stellar winds determine the available envelope mass of the tertiary to be transferred during TMT. A higher metallicity, and thus higher mass-loss rate, would typically be expected to result in a smaller envelope mass of the donor star at the onset of TMT. However, our findings reveal the opposite, as shown in Fig. 7. For systems with metallicities of Z = 0.005, the mass ratios of the outer binary are typically larger than those in lower metallicity systems, as the BHs in the inner binary have lower masses due to higher wind mass loss from their progenitors. As a result, the Roche lobe of the tertiary is smaller, causing the star to fill the Roche lobe at an earlier evolutionary stage. Consequently, these stars frequently engage in mass transfer while still evolving through the Hertzsprung Gap (HG). Conversely, at Z = 0.0005, the majority of tertiary stars are already supergiants, which are much more luminous and have experienced significantly stronger wind mass loss during the pre-TMT evolution. However, the mass-loss rates of massive stars near the Humphreys-Davidson limit that are classified as luminous blue variables (LBVs) are highly uncertain. Current models typically assume the mass-loss mechanism is metallicity independent, which may overestimate the mass lost at low metallicities.

Fig. 7. Distribution of envelope masses of the tertiary star at the onset of the TMT phase at metallicities of Z = 0.005 and Z = 0.0005.

The combined effects of the factors discussed above result in different final properties for the inner binary at the end of TMT. Specifically, in low metallicity systems, the typical inner semi-major axis ranges from 8.7 R_⊙ to 18.1 R_⊙, depending on the physical assumptions of the model. For the high metallicity systems, the separation is, on average, twice as large. Regarding the eccentricities, the distinction between the two populations is much less significant due to the CBD driving the systems towards the same equilibrium eccentricity.

5.1.4. Effect of TMT on outer orbits

Mass transfer from the tertiary to the inner binary inevitably leads to changes in the semi-major axis of the outer orbit due to redistribution or loss of orbital angular momentum. For all systems, TMT occurs in a stable manner (even for donor stars with a convective envelope), largely because of the prevalence of unequal mass ratios in the outer orbit. At metallicities of Z = 0.005 and Z = 0.0005, the mass ratio q_out does not exceed 1 and 0.6, respectively. If the mass ratio were higher, the tertiary would have initially been more massive than the stars of the inner binary, causing it to evolve and fill its Roche lobe before the inner binary became a BBH (Dorozsmai et al. 2024). Since the mass ratio is always less than one, Eq. (14) indicates that the orbit widens throughout the entire mass transfer phase under the assumption of either conservative mass transfer or isotropic re-emission.

5.2. GW sources

Having demonstrated that the inner binary can significantly shrink during TMT, we anticipate efficient formation of GW mergers. Assuming the absence of TMT, already a sizeable fraction of the population would be able to merge solely due to orbital inspiral via GW emission. Specifically, at Z = 0.005 (Z = 0.0005), 44.4% (100%) of the inner binaries are predicted to merge within a Hubble time. When TMT is included, we predict that 85.1–99.9% of systems at Z = 0.005 could form a GW merger within a Hubble time, assuming a CBD is formed. Without a disc, the merger fractions range between 91.5–100%.

Additionally, some mergers occur while TMT is still ongoing, potentially producing a GW signal accompanied by an electromagnetic (EM) counterpart. We discuss this further in Sect. 6.1.1. We predict that < 0.03 − 46.8% and < 0.07 − 29.3% of systems could have an EM signal during the merger with metallicities of Z = 0.005 and Z = 0.0005, respectively. However, if we assume no CBD is formed, up to 100% of systems can produce such mergers. The large uncertainty in these fractions arises from differences in metallicity, mass transfer rate, and physics of the BA and the CBD phase. An overview of all the merger fractions can be found in Table A.1.

5.2.1. Eccentricities with aLIGO

Binary eccentricity can serve as a key indicator to distinguish between different formation channels of GW mergers, as eccentric mergers are often linked to a dynamical history. Various dynamical pathways have been proposed to produce mergers with eccentricities exceeding 0.01–0.1 when entering the aLIGO frequency band at 10 Hz, such as isolated triple stars (Antonini et al. 2017; Liu et al. 2019), galactic nuclei (Antonini & Perets 2012; Takátsy et al. 2019; Gondán & Kocsis 2021; Samsing et al. 2022), and cluster environments (e.g. Antonini et al. 2016; Samsing 2018; Rodriguez et al. 2018a,b; Antonini & Gieles 2020; Dall’Amico et al. 2024). At design sensitivity, current ground-based GW detectors will be sensitive to eccentricities of at least 0.05 (Lower et al. 2018). Future third generation detectors are expected to be sensitive to eccentricities as low as 10⁻³ − 10⁻⁴ (Lower et al. 2018; Saini 2024). We explore whether eccentric BBH mergers can form via the evolutionary channel considered in this study. The GW frequency is determined according to the highest order amplitude (Wen 2003):

$$\begin{array}{c}\hfill f_{\mathrm{GW}}=\frac{\sqrt{G(m_1+m_2)}}{\pi }\frac{{(1+e)}^{1.1954}}{{[a(1-e^2)]}^{1.5}}·\end{array}$$(21)

We have shown that the interactions between the inner binary and the surrounding gas can increase the eccentricity and hence may counteract the circularisation of the orbit due to GW emission, thereby maintaining an eccentric orbit.

As shown in Fig. 8, our predictions indicate a wide range of eccentricities for systems that enter the aLIGO frequency band, with three distinct peaks identified in the distribution when including all additional physical phenomena considered in this study. The largest population, with low eccentricities ranging from 10⁻⁷ − 10⁻⁶, consists of systems whose final phase of the spiral-in is driven purely by GW emission. This peak shows little dependence on metallicity, and is at the same location also characteristic of GW mergers in field binaries and dynamical captures and ejections in globular clusters (e.g. Nishizawa et al. 2017; Zevin et al. 2021a). The second peak, with eccentricities ranging between 10⁻⁵ − 10⁻⁴, arises from systems with retrograde discs that have a merger during the CBD phase. The position of this peak is highly sensitive to the mass-transfer rate, as the time derivative of eccentricity from CBD torques depends on it. Therefore, if the mass transfer takes place on the thermal timescale, the peaks shifts to eccentricities ranging between 10⁻³ − 10⁻². A high eccentricity peak is predicted only for systems with metallicities of Z = 0.005 and high gas densities, reaching eccentricities up to order unity. These systems originate from eccentric gas drag during BA. We find that 19% of these systems reach eccentricities exceeding 0.01. However, we stress that the assumptions for the gas properties during BA are highly simplified and uncertain. When excluding drag forces altogether, the peak at the highest eccentricities disappears.

Fig. 8. Cumulative eccentricity distribution for the GW mergers of Model Advanced with ρg = 10−8 g cm−3 and ṁ3 = ṁnuc when entering the aLIGO frequency band (solid lines) and in the frequency range at which LISA will be most sensitive (dashed and dotted lines).

5.2.2. Eccentricities with LISA

Similar to aLIGO, future space-based GW detectors like LISA could provide valuable insights into the formation history of merging stellar-mass BHs (see Amaro-Seoane et al. 2023). Eccentricity is a key property in this context (Breivik et al. 2016; Nishizawa et al. 2017; Samsing & D’Orazio 2018; Zevin et al. 2019), with predictions indicating that eccentricities of 0.001 should be measurable for 90% of systems within a 5-year observation period, and for 25% of systems within a 2-year period (Nishizawa et al. 2016).

When the inner binaries enter the LISA frequency band at 0.001 Hz, they typically have semi-major axes of 1 − 1.5 R_⊙. By this stage, the systems have generally completed TMT, except for those mergers with a potential EM counterpart. Specifically, systems with an EM counterpart maintain high eccentricities throughout the LISA band because the circularisation due to GW emission is counteracted by CBD torques or eccentric gas drag, which drive the system to elliptical orbits. In Fig. 8, we show the eccentricities at GW frequencies of 0.001 and 0.01 when including all variations of the physical phenomena. We predict similar peaks in the eccentricity distribution as for the aLIGO frequency band, spread across a wide range of eccentricities: 0.0001 − 1. Consequently, the typical eccentricities can be characteristic of all types of evolutionary models, from mergers in isolated binaries to mergers through three-body interactions (see Wang et al. 2024), depending on the physical processes involved in the merger. We note that even without gas drag, high eccentricities (e > 0.1) can be achieved, although it is less common by a factor two compared to cases where eccentric gas drag is included.

5.2.3. Delay times

We define the delay time as the time between the formation of the triple and the merger of the compact objects:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{delay}}={\tau }_{\mathrm{tse}}+{\tau }_{\mathrm{TMT}}+{\tau }_{\mathrm{insp}},\end{array}$$(22)

where τ_tse is the time from the ZAMS up to the onset of TMT, τ_TMT is the duration of the TMT, and τ_insp is the inspiral time of the inner binary when its evolution is only driven by GW emission after the TMT phase. The delay time is an indicator for when a merger occurs and can be translated to a cosmic distance since GWs travel with the speed of light. Given the limited sensitivity of the current ground-based GW detectors, the distance of the merger affects its observability.

Our findings show that the typical delay time shortens due to TMT, regardless of the physical assumptions, as demonstrated in Fig. 9. However, there is a considerable variation in the delay times across different models. Non-conservative mass transfer during the CBD phase leads to delay times around 3–6 times longer compared to conservative mass transfer. Allowing for the formation of retrograde discs shortens delay times by 1.5–2 times compared to only prograde discs. Moreover, systems with low metallicity have significantly shorter delay times compared to high metallicity systems. As discussed in Sect. 5.1.3, the lower metallicity systems typically end up in more compact orbits at the end of TMT. Given the strong dependence of the GW inspiral timescale on the semi-major axis, the delay times at Z = 0.0005 are 1 to 2 orders of magnitude shorter than at Z = 0.005. We assess the impact of these differences in the delay times on the cosmic production rate of GW events in the next subsection.

Fig. 9. Delay time distribution of BBH mergers for different physical assumptions at Z = 0.005 (top) and Z = 0.0005 (bottom). The red histograms show the delay times if TMT had not taken place. The black dashed line denotes the Hubble time of approximately 13.5 Gyr.

5.2.4. Merger rates

Even though the formation efficiency of GW sources is high for the studied channel, triples with a CHE inner binary can only be formed in a limited range of inner semi-major axes, stellar masses, inner mass ratios, and metallicities. We calculate the local merger rates of the inner binaries at redshift z = 0. First, we determine the efficiency at which mergers are formed via this channel within the complete stellar population. This is expressed as

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{eff}}=f_{\mathrm{CHE}}\frac{N_{\mathrm{merge}}}{N_{\mathrm{simulated}}},\end{array}$$(23)

where f_CHE is the fraction of the simulated parameter space relative to the entire parameter space in which stars are born. For details on the underlying sampling distributions, we refer the reader to Sect. 3.2. Assuming a triple fraction of 0.73, a binary fraction of 0.21, and a single star fraction of 0.06 within the simulated mass range, the value of f_CHE is 6.5 × 10⁻⁵. The calculation of f_CHE is identical to that of Dorozsmai et al. (2024). N_merge is the number of GW mergers produced in our simulations, while N_simulated is the total number of triple systems simulated. Next, we translated this merger efficiency to a birth-rate density, which gives the number density of merging systems born per unit star-forming mass:

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{birth}}=\frac{1}{\stackrel{\sim }{M}}{\displaystyle \sum _{Z_i}}\text{SFRd}(Z_i,z_{\text{form}}){ϵ}_{\mathrm{eff}},\end{array}$$(24)

where SFRd(Z_i, z_form) is the star-formation rate density at the redshift at which the system formed following Madau & Dickinson (2014), Madau & Fragos (2017), and $\stackrel{\sim }{M}$ is the average mass per system. The term Z_i is the metallicity of the system, while z_form is the formation redshift. In this study, we considered only two metallicities despite the broad range of metallicities at which mergers can form. This limitation introduces some uncertainty into the calculated rates. Finally, we computed the merger rate by integrating the birth rate over the delay times:

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{merge}}=\frac{1}{\stackrel{\sim }{M}}{\displaystyle \sum _{Z_i}}{\displaystyle {\int }_0^{t_{\mathrm{merge}}}}\text{SFRd}(Z_i,z_{\text{form}}(t_{\mathrm{delay}})){ϵ}_{\text{eff}}({\mathrm{t}}_{\mathrm{delay}})d{t}_{\mathrm{delay}},\end{array}$$(25)

where t_merge is the time at which the merger occurs. For a more detailed description of how the merger rates were calculated, we refer to Dominik et al. (2015), Dorozsmai et al. (2024).

The predicted local (redshift z = 0) merger rates range between 0.69 − 1.74 Gpc⁻³ yr⁻¹ depending on the assumed physics in the model. The variation in these rates is mainly a result of the different delay times for the mergers between the models. To clarify, the range in the predicted rates are solely a result of the different physics assumption during the TMT and does not account for uncertainties in the initial distributions of stellar and orbital properties or the preceding stellar evolution. For instance, Kummer et al. (2023) show that rare evolutionary pathways in hierarchical triples can introduce uncertainties of up to a factor of a few. Interestingly, the high metallicity population contributes significantly more – by a factor 10–30 – towards the local merger rate than the low metallicity population. This discrepancy is due to the short delay times in the low metallicity population, which necessitates the formation of the original triples at a low redshift for them to be observable locally. However, at low redshifts, the star-formation rate (SFR) is relatively low, and stars tend to be born with higher metallicities, resulting in a low merger rate for this population.

Regarding mergers with a potential EM counterpart, the local merger rate varies substantially across the different models. The predicted rates range from as high as 0.62 Gpc⁻³ yr⁻¹ to almost negligible values. Models that predict higher rates typically include retrograde discs, eccentric gas drag, and high gas densities, as these conditions more efficiently facilitate GW mergers during TMT. However, observing these mergers locally is challenging because their delay times are generally short, often not exceeding a few tens of millions of years. An overview of all merger rates can be found in Table A.2.

Figure 10 illustrates the merger rates for both metallicity populations across redshift in our Model Advanced setup. The merger rate for both high and low metallicity populations increases steeply towards z = 2, aligning closely with the peak of the cosmic star-formation rate. Beyond this point, the merger rate for the high metallicity population (Z = 0.005) rapidly declines, while the merger rate for the low metallicity population (Z = 0.0005) decreases more gradually. At redshifts greater than z = 5, most of the mergers originate from the low metallicity stars. This trend results from the combination of shorter delay times and the peak of the metallicity distribution shifting towards lower values at higher redshifts.

Fig. 10. Binary BH merger rate as a function of redshift for Model Advanced at Z = 0.005 (blue) and Z = 0.0005 (orange). The dotted lines refer to the merger rate of triples with a CHE inner binary where the tertiary has not influenced the evolution of the inner binary.

Variations of the SFR and the cosmic metallicity distribution can significantly impact predictions for compact object merger rates (e.g. Neijssel et al. 2019; Chruslinska et al. 2019; Broekgaarden et al. 2021, 2022). This is especially true at high redshifts (z > 2), where the SFR is not well established. Such uncertainty primarily affects local GW mergers from the high metallicity population in our study, as the short delay times of low metallicity mergers indicate that such mergers are only locally observable if the system was formed at a low redshift. Moreover, recent cosmological simulations of galaxy formation and evolution have identified a low metallicity tail at low redshifts (e.g. van Son et al. 2023), which is not captured byg conventional log-normal distributions based solely on empirical evidence. While we do not quantify these uncertainties, they could significantly increase the merger rate for the low metallicity population.

Previous studies have already predicted redshift-dependent BBH merger rates for CHE binaries, but without a tertiary component (Mandel & de Mink 2016; Riley et al. 2021). To compare with their findings, we computed the merger rate for systems where the tertiary star does not significantly influence the evolution of the inner binary, either through mass transfer or dynamical interactions; that is, the inner binary evolves effectively isolated from the tertiary. As shown in Fig. 10, the most pronounced difference with the studied triple channel is that the peak of the merger rate shifts towards lower redshift, particularly for the high metallicity population. This shift occurs because TMT generally shortens the merger delay time. Mandel & de Mink (2016) considered a population at a fixed metallicity of Z = 0.004, similar to our high metallicity population. Their merger rate peaks around a redshift of 0.3–0.4, somewhat higher than our findings, possibly due to minor differences in metallicity. Conversely, Riley et al. (2021) integrated the merger rate over a broader range of metallicities, using 30 evenly spaced metallicity bins, and find that the rate typically peaks at higher redshifts (z = 3 − 4). This discrepancy, however, might be attributed to differences in the sampled metallicity range, the number of different metallicities considered, the SFR model, and redshift-dependent metallicity distribution of stars.

6. Discussion

In this section we address additional key observational signatures associated with GW mergers and examine the limitations of the mass transfer model. Specifically, we discuss the assumptions made related to the CBD and BA.

6.1. Gravitational wave sources

6.1.1. Electromagnetic signal

A merger of two BHs preceded (and accompanied) by strong EM emission would be a smoking gun to pin down the evolutionary origin of the merger. Stone et al. (2017) found that BBHs embedded in an AGN disc could be driven towards mergers assisted due to three-body interactions and interactions with the surrounding disc, potentially producing a luminous EM counterpart upon merger. Presumably, the X-ray emission is constrained by the Eddington mass accretion rate of the accreting object:

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}=1.3\times {10}^{38}(\frac{M}{M_{\odot }})\mathrm{erg}{\mathrm{s}}^{-1}.\end{array}$$(26)

However, if the mass transfer rate exceeds the Eddington accretion rate significantly, the unaccreted mass is speculated to be expelled, obstructing the emitted X-ray radiation in most directions, thereby directing the emission towards the polar regions above and below the accretion disc (e.g. King et al. 2023). This leads to an observed X-ray luminosity in excess of 10³⁹ erg s⁻¹ for neutron stars and stellar-mass BHs, classified as ultraluminous X-ray (ULX) sources. To date, the observed number of ULX candidates ranges in the thousands (Walton et al. 2022), with many originating from low-metallicity environments (Soria et al. 2005; Mapelli et al. 2009; Prestwich et al. 2013; Basu-Zych et al. 2016; Kovlakas et al. 2020).

If the inner binaries of our triples manage to coalesce while the TMT phase is ongoing, this pathway could give rise to BBH mergers accompanied by EM signals. The BHs in our sample typically have Eddington mass accretion rates of 10⁻⁷ − 10⁻⁶ M_⊙ yr⁻¹. The mass-transfer rate of massive tertiary stars, especially on the thermal timescale, far exceeds Eddington limit, leading to high X-ray luminosities. It is important to emphasise that the frequency of mergers with a disc still present heavily depends on the underlying model assumptions (see Sect. 5.2).

6.1.2. BH spin

In this study, we have not considered the spin distribution of merging BHs. The effective spin χ_eff of the merging binary, which contains both the magnitude and orientation of the spin vector for the binary system, is a property that can be inferred directly from the GW signal (Abbott et al. 2023). Disentangling the individual spin components poses a significant challenge. Previous studies of CHE in isolated binaries predict BH spins that stand out from other GW channels in terms of their high spin magnitudes retained from rapid rotation during the CHE phase (e.g. Zevin et al. 2021b), although Marchant et al. (2024) predict that the majority of locally detected BBHs from this channel have χ_eff < 0.5. In our sample, interactions between the disc and the binary could potentially alter the spin properties of the merging BHs. For instance, in inclined accretion discs, the BHs are driven to alignment by the Lense-Thirring effect (e.g. Volonteri et al. 2005), although the spins of the BHs are already expected to be aligned with respect to each other. An additional consequence of Lense-Thirring is the alignment between the angular momentum vector of the binary and the disc (and thus the tertiary), which could suppress dynamical interactions within the triple. This would prevent misalignment due to three-body interactions during the successive inspiral of the binary. Furthermore, accretion of material could result in spin-up of the objects. However, Muñoz et al. (2020) find that the time-averaged spin torques are only a few percent of the rate of change in orbital angular momentum. A more detailed investigation of the spin evolution in this channel is necessary to accurately determine the final spin distribution, but we anticipate the spins to be relatively similar to CHE in isolated binaries.

6.1.3. Pulsational pair-instability supernovae

Stellar rotation decreases the initial mass necessary to form helium and carbon-oxygen (CO) cores of a given mass (Chatzopoulos & Wheeler 2012), an effect that is particularly relevant for stars that spin so fast that they evolve chemically homogeneous (Woosley 2017). For helium and CO core masses of about 30 M_⊙ or more, this leads to pulsational instabilities in the very final phases of evolution. These may lead to very substantial mass shedding, reducing the mass and angular momentum of the BH that forms in the (pulsational pair-instability) supernova that ensues (see also, Marchant et al. 2019; Stevenson et al. 2019).

For the explored parameter space (i.e. initial masses of the primary component of the inner binary between 20 and 100 M_⊙) 7.5% of systems evolves chemically homogeneous, out of which 9.5% lead to mass transfer from the tertiary to an inner BBH. From these, for Z = 0.005, 50% of the inner binary components ultimately have helium and CO core masses in the range 35 − 60 M_⊙. For Z = 0.0005 this is 99%. So, the majority of our systems may be affected by pulsational pair instability effects not accounted for in our simulations. The primary effect the pair instability driven mass loss is anticipated to have is that it drives the inner binary towards larger separations. The implications of this shift on the merger process remain ambiguous. On the one hand, the inspiral time should increase as the binary widens. On the other hand, three-body dynamical interactions become stronger, leading to more eccentric orbits. Additionally, the final BH masses and angular momenta would be reduced. The former influences the evolution of both the inner and outer orbits throughout the entire TMT phase and subsequent GW inspiral. Furthermore, the mass transfer is more prone to be unstable.

6.2. Assumptions regarding circumbinary discs

6.2.1. Viscosity parameter

The evolution of the inner binary in the presence of a CBD in our model are based on simulations tailored for SMBHs (Siwek et al. 2023b), though our investigation primarily concerns stellar mass objects, which are significantly lower in mass. Notably, the viscosity parameter α in discs around stellar mass systems are typically around an order of magnitude lower than those around SMBHs (Hartmann et al. 1998). Dittmann & Ryan (2022) found that for low aspect ratios, H/r < 0.1, the transfer rate of specific angular momentum between the disc and the binary is highly dependent on the viscosity of the disc. However, at H/r = 0.1, this dependency more or less disappears, resulting in consistent semi-major axis evolution regardless of the disc’s viscosity. This result is in agreement with the findings of Muñoz et al. (2020). Given that stellar-mass discs typically have aspect ratios closer to 0.1 (see Lai & Muñoz 2023), it suggests that the binary’s evolution may not deviate significantly.

6.2.2. Steady-state disc

In our simulations, we assumed that a steady state disc is reached instantaneously, wherein the mass accretion rate onto the inner binary matches the supply rate at the outer boundary of the disc. However, this assumption overlooks the so-called ‘transient phase’, characterised by the cavity’s filling within the CBD and the formation of circumstellar discs encircling the individual binary components. This phase may involve loss of angular momentum from the inner binary as the net torque is provided solely by gravitational interactions (Muñoz et al. 2020), leading to an additional decrease in the semi-major axis. Although the duration of the transient phase is dependent on the initial disc properties, we performed a back-of-the-envelope calculation to assess the significance of this phase. Following Muñoz et al. (2020), if a steady state is reached after a few hundred binary orbits, its duration could extend to approximately 10 years. However, for stronger viscous discs, resembling discs around stellar-mass objects instead of SMBHs, the transient phase is anticipated to last even longer (e.g. Miranda et al. 2017). Considering the thermal timescale of massive evolved stars, estimated at about 10–100 years, it becomes apparent that a significant portion of the mass transfer phase could be governed by the transient phase, and therefore the orbital shrinkage could be underestimated in the current model.

Similarly, we assumed that the disc promptly disperses as soon as the mass-transfer phase is terminated. However, for a finite disc in quasi-steady state (where the mass-accretion rate scales with the rate at which the disc viscously spreads), the mass accretion onto the binary diminishes with time, following a time dependency of approximately ṁ_acc ∝ t^−4/3 (Muñoz et al. 2020; Siwek et al. 2023a). As a result, it is conceivable that our estimation of the disc’s lifetime is underestimated. Increased disc lifetimes, coupled with strong GW emission, may lead to an increase in the number of mergers in the presence of a CBD, consequently leading to a higher occurrence of mergers accompanied by an EM signal. This could be particularly pronounced at lower metallicity, where delay times can be less than 10 Myr, and GW emission plays a significant role in the orbital evolution. This scenario is especially relevant for massive post-MS tertiary stars, as their evolutionary timescale typically spans 0.1–1 Myr for the tertiary stars considered in this study.

6.2.3. Comparison with other studies

Similar to our approach, Valli et al. (2024) use an interpolation method to analyse changes in semi-major axis and eccentricity, referencing the work of Siwek et al. (2023b), and compare these results with those from Zrake et al. (2021) and D’Orazio & Duffell (2021). Notably, disagreements arise regarding the stability of the eccentricity in initially circular orbits with equal mass ratios. While Siwek et al. (2023b) suggest the configuration is unstable, leading to an increase in eccentricity and ultimately a decreasing orbital separation, the two previous studies indicate an additional equilibrium eccentricity exists at e = 0 characterised by orbital widening. At the onset of TMT, fewer than 1% (2.5%) of the inner binaries have eccentricities below 0.01 for the Z = 0.005 (Z = 0.0005) population. Predominantly, inner binaries at both metallicities have eccentricities around 0.1 upon TMT onset, largely attributed to the Blaauw kick during BH formation. However, we assumed 10% of the mass is lost as neutrinos, which might be an overestimation (Vigna-Gómez et al. 2024) and lead to an overestimation of the number of eccentric orbits. Nevertheless, many triples are influenced by three-body dynamical effects, increasing the eccentricity regardless of the occurrence of a supernova kick. Therefore, the potential equilibrium eccentricity at e = 0 is not of great importance to our systems.

6.3. Assumptions regarding ballistic accretion

6.3.1. Properties of gas around the inner binary

It is important to recognise that uncertainties related to the gas properties could lead to considerable variations in the evolution of the binary during TMT, potentially influencing the findings presented in this study.

For the duration of the BA phase, we relied on simplified assumptions regarding the gas enveloping the inner binary. These simplifications stem from our limited knowledge, largely due to the scarcity of detailed studies into stable mass transfer from tertiary stars. To encompass a wide range, we chose two significantly different gas densities (ρ_gas = 10⁻⁸ g cm⁻³ and ρ_gas = 10⁻¹⁰ g cm⁻³), which are in reasonable agreement with densities observed in previous research (de Vries et al. 2014). Additionally, we imposed a constraint that the density remains several orders of magnitude lower than that of a typical common envelope structure. In Sect. 5.1 we demonstrated that the final orbital separations of the inner binary at the end of TMT are more similar when considering the formation of a CBD, even for different densities of the gas during BA. The choice for the sound speed and gas temperature is uncertain. However, given the high orbital velocities of the inner binary components (hundreds of km/s), their motion remains highly supersonic under any realistic assumption, meaning the GDF is not significantly impacted. If the sound speed were much higher, inner orbits that become highly eccentric due to eccentric gas drag could transition to subsonic motion at apocentre before GW emission dominates the evolution, halting the eccentricity growth (e.g. O’Neill et al. 2024). This could reduce the predicted fraction of GW mergers with high eccentricities in the aLIGO frequency band.

Additionally, in our model we assumed that the density of gas surrounding the inner binary remains constant during the BA phase. This assumption implies a balance between material ejected from the inner binary and mass transferred from the tertiary. If we assume the gas density increases in time, starting from a high initial gas density, the final properties of the inner binary are not expected to change much, as these are mainly determined by the CBD phase. Conversely, starting from a low initial gas density, the final orbital properties might be different if most or the entire mass transfer takes place during BA. Alternatively, we consider a scenario where the gas density around the inner binary decreases over time. This leads to depletion of the gas cloud, potentially completely halting the inspiral. However, the cloud continuously receives mass from the tertiary, maintaining some gas around the binary during the TMT phase. Characterising the properties and evolution of such gas clouds would require more detailed studies of stable TMT.

Another assumption we have made is that the density of the gas is completely homogeneous and has no dependence on the distance from the COM, resulting in efficient inspiral of the binary. Antoni et al. (2019) argue that if the COM is at rest relative to the gas, the density near the COM becomes enhanced, leading to a Bondi density profile where ρ_g ∝ a^−3/2. Consequently, the inspiral becomes more efficient for supersonic motion and less efficient for subsonic motion. Since we expect the compact objects in our study to typically move at supersonic velocities, the transition time from BA to a CBD would become shorter. Furthermore, a centrally concentrated mass distribution would lead to a less efficient increase in eccentricity, or even a decrease for a sufficient high radial dependence of the gas density (Szölgyén et al. 2022). This would affect the subsequent CBD evolution, but its impact on merger frequency and properties is not evident.

The GDF experienced by binary objects is significantly influenced by the relative velocity between the gas and the inspiraling objects. This dependence arises because $F_{\mathrm{drag}}\propto v_{\mathrm{rel}}^{-2}$ for supersonic motion relative to the gas. This relation originates from the accretion radius of the binary objects, which shares the same proportionality with the relative velocity. In our simulations, we have assumed the gas to be at rest with respect to the COM of the inner binary. This is not a realistic assumption, as the gas stream enters the vicinity of the inner binary with an initial velocity relative to the COM and is subsequently disturbed upon intersecting with the binary orbit. A prograde (retrograde) direction of the gas velocity would lead to more (less) efficient binary inspiral during BA.

The findings of Ostriker (1999), Kim & Kim (2007) on the GDF were derived under the assumption of an infinite background. In our scenario, the binary is embedded within a finite amount of gas during BA. A back-of-the-envelope calculation suggests that, given the gas densities under consideration, the total mass of the surrounding gas cloud is likely much lower than the mass of the inner binary components. Although the precise impact of this difference on the results of Ostriker; Kim & Kim is unclear, it is expected that the efficiency of the GDF is reduced, resulting in a less effective inspiral of the binary

6.3.2. Hydrodynamic drag during ballistic accretion

Computing the orbital evolution of the binary embedded in a gaseous medium in our model, we have only considered the effect of the GDF that arises through the interaction between the overdense wake and the orbiting bodies. Additional drag forces can be imparted by direct impact of the gas on the bodies, resulting in a loss of angular momentum. This hydrodynamic drag can be expressed as

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{drag}}=-\frac{1}{2}{C}_D\pi R^2{\rho }_g{v}_{\mathrm{rel}}{\mathbf{v}}_{\mathrm{rel}},\end{array}$$(27)

where C_D is the drag coefficient. The balance between both types of drag force can be estimated in terms of the compactness, M/R, of the objects (Grishin & Perets 2015; Szölgyén et al. 2022). If the compactness exceeds a critical value, the GDF dominates over the hydrodynamic drag force, and vice versa. In this study, we applied the TMT model to inner binary objects consisting of BHs, where we found the hydrodynamic drag force to be negligible as the accretion radius of the BHs, R_acc = 2Gm/v_rel², is significantly larger than their Schwarzschild radius. Therefore, considering only hydrodynamic drag, we find a similar evolution of the systems as when drag is completely neglected. However, if the mass transfer phase is initiated while the stars of the inner orbit are still on the MS, as can happen when the tertiary star is initially the most massive object in the triple (e.g. Toonen et al. 2020), the contribution of the hydrodynamic drag is not evident.

We estimated the parameter space where the hydrodynamic drag is important. Ignoring the constants that are approximately equal to one, the gravitational drag is dominant for each of the components if the following condition is satisfied (Szölgyén et al. 2022):

$$\begin{array}{c}\hfill \frac{M_i}{R_i}\gtrsim \frac{v_{i,\mathrm{rel}}^2}{2G}·\end{array}$$(28)

For MS stars, the compactness increases towards larger masses, as the mass and radius relate as approximately R ∝ M^0.7. Furthermore, the orbital velocity of the more massive component in a binary is lower than the velocity of its companion. Therefore, we only compared the drag forces for the secondary star. If the relative velocity is equal to the Keplerian velocity, we can express the equation in terms of semi-major axis and binary mass:

$$\begin{array}{c}\hfill \frac{M_2}{R_2}\gtrsim \frac{M_{\mathrm{bin}}}{2a(1+\frac{M_2}{M_1}{)}^2}·\end{array}$$(29)

We considered a 1 M_⊙ secondary and calculated the ratio between the competing drag forces for a grid of mass ratios in the range 0.01–1 and semi-major axes in the range 0.01 − 10 R_⊙. We find that the hydrodynamical drag force becomes important towards unequal mass ratios and small separations, as shown in Fig. 11. However, when only considering detached systems, the hydrodynamical drag is non-negligible only at separations of at least about 10 R_⊙ and mass ratios close to 0.1. For any other orbital configurations, the GDF dominates the total drag. However, if gravitational drag is absent, hydrodynamic drag could still lead to changes in the orbital properties. If the gas were to move move prograde (retrograde) with respect to the binary, the hydrodynamic drag becomes less (more) important.

Fig. 11. Dominant type of drag force (gravitational versus hydrodynamic) as a function of the semi-major axis and the mass ratio for a secondary with M2 = 1 M⊙. The red dashed line indicates the radius of the primary (R1) at the ZAMS. The dotted and dash-dotted lines refer to the Roche radii of the primary at a semi-major axis of 3 R⊙ and 8 R⊙, respectively. If the radius exceeds the size of the Roche lobe, the system should not exist.

7. Conclusions

We have presented a rapid analytical model designed to probe the evolution of triple star systems undergoing stable mass transfer from the tertiary component to the inner binary. This model serves as a tool for conducting large-scale population simulations. Our analysis describes two primary modes of mass transfer evolution: ballistic accretion, which is characterised by the inner binary engulfed within a gaseous cloud and governed by the GDF, and CBD evolution, wherein gravitational and accretion torques dictate orbital changes. We applied this model to a simulated population of triples with a Roche-lobe filling tertiary star transferring mass onto an inner binary consisting of two BHs at metallicities of Z = 0.005 and Z = 0.0005. Both components of the inner binary experienced chemically homogeneous evolution throughout their MS phase. We summarise our key findings as follows:

• For triples with a CHE inner binary, tertiary mass transfer is typically stable due to the mass of the tertiary being lower than the mass of the inner binary. Among the simulated systems, 42.9% (25.6%) are not able to form a CBD at the onset of TMT at metallicities of Z = 0.005 (Z = 0.0005), mainly due to small ratios between the outer and inner orbital separation, and start out the evolution with BA.

• The evolution of the inner binary during BA invariably leads to a decrease in the semi-major axis of the inner binary on a timescale much shorter than that of the CBD evolution if gas drag is efficient. Additionally, for initially eccentric inner binaries, drag forces can rapidly increase the eccentricity towards unity, leading to efficient formation of GW mergers. CBD evolution tends to shrink the orbit for the simulated range of near-equal mass ratio binaries. However, if the mass transfer is completely non-conservative, it may lead to widening of the inner orbit during this phase. Eccentricity typically converges towards an equilibrium value of approximately 0.48 under CBD interactions in prograde orbits. However, in cases where retrograde discs are formed, the eccentricity increases towards unity.

• While many systems start out with BA, the evolution of the orbital properties typically leads to the formation of a CBD. Whether a CBD can form primarily depends on the balance between the rate at which the inner semi-major axis shrinks and the outer orbit expands (promoting CBD formation) versus the rate at which the inner eccentricity increases (hindering CBD formation). If a CBD does form, the importance of the efficiency of the binary inspiral during BA – dictated by the properties of the gas surrounding the inner binary – diminishes, as the CBD regulates the inspiral.

• Even though the formation efficiency of systems that evolve through the evolutionary channel investigated in this study is low (< 1% for systems with initial masses of the primary component of the inner binary between 20 and 100 M_⊙ and separations of the outer orbit up to 10⁵ R_⊙), the channel is highly efficient in producing GW mergers. For systems with a metallicity of Z = 0.005, we predict that 85.1–100% of inner binaries will merge within a Hubble time, depending on the specific physical model assumptions. At a metallicity of Z = 0.0005, the efficiency reaches 100%. This higher efficiency at Z = 0.0005 is primarily due to minimal binary widening caused by stellar winds during pre-supernova evolution, resulting in inner orbits that are, on average, 30 R_⊙ less wide compared to systems with a metallicity of Z = 0.005.

• The total local merger rate spans a range of 0.69–1.74 Gpc⁻³ yr⁻¹, with the primary contribution originating from higher metallicity systems, despite their lower efficiency in forming GW mergers. This behaviour arises due to a combination of factors, including the delay time distribution and the metallicity-specific star-formation rate. Notably, low-metallicity binaries have delay times only up to a few hundred million years, suggesting that they predominantly merge at high redshifts, where the formation of low metallicity stars is favoured. The predicted merger rates only account for uncertainties in the physics assumptions during TMT.

• Our predictions suggest there are two types of GW mergers: those occurring during TMT and those taking place after the TMT phase has finished. The former type of mergers has the potential to produce an EM signal before and during the merger event. However, the occurrence varies significantly depending on the model assumptions. Considering only prograde discs and ignoring eccentric gas drag, we predict an EM signal in up to 5.6% of mergers at low metallicity, whereas including retrograde discs and eccentric gas drag could yield rates as high as 46.8%. Even when ignoring gas drag completely during BA, we predict up to 19.2% of mergers occurring during TMT. However, these rates are subject to large uncertainties, mainly due to the assumptions made for the initial triple orientations and the gas properties during BA.

• Our study reveals that triples with CHE inner binaries share several properties with other formation pathways but also present distinct differences. Compared to CHE in isolated binaries, these systems are expected to have similar mass ratios and effective spins. However, unlike the isolated binary case, our findings indicate that systems with high eccentricities can be formed in the aLIGO and LISA frequency band, under the assumption that retrograde discs and eccentric gas drag play a role in the evolution of the population. Unlike dynamical formation pathways involving classically evolving triples or dynamical interactions in clusters, the effective spin in triples with a CHE inner binary is expected to be typically non-zero. Nevertheless, for most systems, the predicted eccentricities closely align with those found in and around cluster environments.

In conclusion, TMT plays an important role in shaping the evolution of the inner binary within hierarchical triple star systems. In particular, for triples with a CHE inner binary, the mass transfer facilitates the formation of GW sources, and can result in mergers with potentially unique properties. However, the mass transfer model is subject to numerous uncertainties, and further detailed studies are essential for a complete understanding of the evolution of these systems.

Data availability

The data and scripts necessary to run the model and reproduce the figures in this paper are publicly available on Zenodo: 10.5281/zenodo.13610538.

Acknowledgments

FK would like to thank Alejandro Vigna-Gómez, Magdalena Siwek, Mor Rozner, Morgan Mcleod and Nathalie Degenaar for the helpful discussions. The authors also acknowledge support from the Netherlands Research Council NWO (VENI 639.041.645 and VIDI 203.061 grants).

References

Appendix A: Overview merger rates

Appendix B: Orbit averaged evolution due to gravitational drag

In this section we derive the orbit-averaged derivatives of the semi-major axis and eccentricity for a binary subjected to GDF. Built upon the expression for a single orbiting object (Murray & Dermott 1999; Grishin & Perets 2016), the change in the semi-major axis over time can be described by

$$\begin{array}{c}\hfill \frac{\mathit{da}}{\mathit{dt}}=2\frac{a^{3/2}}{\mu \sqrt{GM(1-e^2)}}[\mathrm{\Delta }{F}_{r}esin\nu +\mathrm{\Delta }{F}_{\varphi }(1+ecos\nu )],\end{array}$$(B.1)

where μ = m₁m₂/M is the reduced mass, M = m₁ + m₂ is the total mass of the binary, ΔF_r and ΔF_ϕ are the radial and azimuthal components of the GDF resulting from the differential acceleration between the binary objects, and ν is the true anomaly of the orbit. The force is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathit{F}=\frac{A_0{m}_1^2}{v_1^3}{\mathit{v}}_1-\frac{A_0{m}_2^2}{v_2^3}{\mathit{v}}_2=A_0{M}^2(\frac{1}{q^2}+q^2)\frac{{\mathit{v}}_b}{v_b^3},\end{array}$$(B.2)

where A₀ = 4πG²ρ_gℐ, and the velocities are expressed in terms of the binary velocity: v₁ = (m₂/M)v_b and v₂ = (m₁/M)v_b. The minus sign in the expression becomes positive because the velocity vector of the two objects point in opposite direction. The velocity vector in the COM frame is

$$\begin{array}{c}\hfill {\mathit{v}}_b=\frac{\mathit{na}}{\sqrt{1-e^2}}[esin\nu \widehat{r}+(1+ecos\nu )\widehat{\varphi }],\end{array}$$(B.3)

where n = 2π/P and the velocity magnitude is defined as $v_b=\sqrt{v_{b,r}^2+v_{b,\varphi }^2}$. Assuming the gas is at rest with respect to the COM, substituting these expressions into Eq. B.1 yields

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{\mathit{da}}{\mathit{dt}}=\frac{2a^{3/2}}{\mu \sqrt{GM(1-e^2)}}\frac{A_0{M}^2(1-e^2)}{n^2{a}^2}(\frac{1}{q^2}+q^2)\frac{1}{\sqrt{1+2ecos\nu +e^2}}\\ \hfill =\frac{2A_0{M}^2\sqrt{1-e^2}}{\mu n^3{a}^2}(\frac{1}{q^2}+q^2)\frac{1}{\sqrt{1+2ecos\nu +e^2}}.\end{array}\end{array}$$(B.4)

Here, we used $\sqrt{\mathit{GM}}=a^{3/2}/n$. To calculate the average derivative over one orbital period, we have

$$\begin{array}{c}\hfill \langle \frac{\mathit{da}}{\mathit{dt}}\rangle \equiv \frac{1}{P}{\displaystyle {\int }_0^P}\frac{\mathit{da}}{\mathit{dt}}dt=\frac{n}{2\pi }{\displaystyle {\int }_0^{2\pi }}\frac{\mathit{da}}{\mathit{dt}}\frac{{(1-e)}^{3/2}}{n{(1+ecos\nu )}^2}d\nu ,\end{array}$$(B.5)

where dt = dν(1 − e)^3/2/[n(1 + e cos ν)²]. Combining Eqs. B.4 and B.5, we obtain

$$\begin{array}{c}\hfill \langle \frac{\mathit{da}}{\mathit{dt}}\rangle =\frac{A_0{(1-e^2)}^2}{\pi n^3{a}^2}\frac{M^2}{\mu }(\frac{1}{q^2}+q^2){\displaystyle {\int }_0^{2\pi }}\frac{\mathit{dv}}{{(1+ecos\nu )}^2\sqrt{1+2ecos\nu +e^2}}.\end{array}$$(B.6)

For convenience, we kept ℐ out of the integral, though it depends on the true anomaly. This assumption might influence the accuracy of the calculated values for the semi-major axis derivative. In the case of a circular binary, this expression simplifies to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \langle \frac{\mathit{da}}{\mathit{dt}}\rangle =\frac{2A_0}{n^3{a}^2}\frac{M^2}{\mu }(\frac{1}{q^2}+q^2)=\frac{8\pi G^2}{n^3{a}^2}{\rho }_g\mathcal{I}\frac{M^2}{\mu }(\frac{1}{q^2}+q^2)\\ \hfill =8\pi \sqrt{\frac{G{a}^5}{M}}{\rho }_g\mathcal{I}\frac{M}{\mu }(\frac{1}{q^2}+q^2)\\ \hfill =8\pi \sqrt{\frac{G{a}^5}{M}}{\rho }_g\mathcal{I}[\frac{1}{q}{(1+q^{-1})}^2+q{(1+q)}^2].\end{array}\end{array}$$(B.7)

The factor between the brackets is similar to the one presented in Rozner & Perets (2024).

Using the same approach, we find the orbit-averaged derivative of the eccentricity. The change in the eccentricity over time is

$$\begin{array}{c}\hfill \frac{\mathit{de}}{\mathit{dt}}=\frac{1}{\mu }\sqrt{\frac{a(1-e^2)}{\mathit{GM}}}[\mathrm{\Delta }{F}_{r}sin\nu +\mathrm{\Delta }{F}_{\varphi }(cos\nu +cosE)],\end{array}$$(B.8)

where E is the eccentric anomaly, defined as cos E = (e + cos ν)/(1 + e cos ν). Following the same steps as for the semi-major axis, the final expression is

$$\begin{array}{c}\hfill \langle \frac{\mathit{de}}{\mathit{dt}}\rangle =\frac{A_0{(1-e^2)}^3}{\pi n^3{a}^3}\frac{M^2}{\mu }{k}_q{\displaystyle {\int }_0^{2\pi }}\frac{(e+cos\nu )dv}{{(1+ecos\nu )}^2{(1+2ecos\nu +e^2)}^{3/2}},\end{array}$$(B.9)

where we have defined $k_q=(\frac{1}{q^2}+q^2)$.

All Tables

All Figures

Fig. 1. Illustration of the inner and outer orbits in a hierarchical triple system with masses m1, m2, and m3 along with their orbital elements: semi-major axes, eccentricities, and mutual inclination. The mass streams of a Roche-lobe-filling tertiary that lead to BA (intersecting with the inner orbit) or the formation of a CBD (intersecting with itself) are depicted. We note that the figure is not to scale.

In the text

Fig. 2. Constants for the time derivatives of the semi-major axis (top) and the eccentricity (bottom) of an inner binary surrounded by a CBD based on the findings of Siwek et al. (2023b). The figure illustrates how certain parts of the parameter space lead to either a decrease or increase of the semi-major axis and the eccentricity. The dashed lines indicate the points where the derivative equals zero.

In the text

Fig. 3. Evolutionary stages of triples with a CHE inner binary that result in TMT with an inner BBH. Depending on the properties of the triple at the onset of TMT, mass transfer occurs via BA or the formation of a CBD. The masses of the three companions are indicated at the ZAMS (MZAMS), MS (MMS), and helium MS (MHeMS) phases. Furthermore, the BH mass (MBH), pre-tertiary mass-transfer mass (Mpre − TMT), envelope mass (Menv), and core mass (Mcore) are indicated. The grey and red crosses denote the centre of mass of the inner binary and outer binary, respectively. The figure is not to scale and has been adapted from van Son et al. (2022).

In the text

Fig. 4. Evolution of the semi-major axis and eccentricity of the inner binary of an example system during the TMT phase and the subsequent GW inspiral. Different models are applied to the same system, as described in Sect. 2.3. The red shaded region indicates that the system undergoes BA. The green shaded region indicates that the system has formed a CBD. The yellow region implies that the mass transfer is non-conservative during the CBD phase, and the material is assumed to be lost via an isotropic outflow. The blue shaded region indicates that the inspiral of the inner binary is dominated by GW emission. Overlapping colours show that more than one of the aforementioned processes apply.

In the text

Fig. 5. Cumulative distribution of the semi-major axis of the inner binaries at the end of TMT for Model Basic at Z = 0.005. The few systems that have merged during TMT are omitted. The black dashed line denotes the semi-major axis distribution at the onset of mass transfer. Among different assumption for the gas density during BA and the mass transfer rate, the final distributions are fairly similar.

In the text

Fig. 6. Distribution of the semi-major axis and eccentricity of the inner binaries with a metallicity of Z = 0.005 (top) and Z = 0.0005 (bottom) at the onset of TMT. The purple and red shaded regions correspond to systems that start the TMT phase with the formation of a CBD and BA, respectively. The red bars are plotted on top of the purple bars. We note that the range of the semi-major axis is different between the upper and lower panel.

In the text

	Fig. 7. Distribution of envelope masses of the tertiary star at the onset of the TMT phase at metallicities of Z = 0.005 and Z = 0.0005.
In the text

	Fig. 8. Cumulative eccentricity distribution for the GW mergers of Model Advanced with ρg = 10−8 g cm−3 and ṁ3 = ṁnuc when entering the aLIGO frequency band (solid lines) and in the frequency range at which LISA will be most sensitive (dashed and dotted lines).
In the text

	Fig. 9. Delay time distribution of BBH mergers for different physical assumptions at Z = 0.005 (top) and Z = 0.0005 (bottom). The red histograms show the delay times if TMT had not taken place. The black dashed line denotes the Hubble time of approximately 13.5 Gyr.
In the text

	Fig. 10. Binary BH merger rate as a function of redshift for Model Advanced at Z = 0.005 (blue) and Z = 0.0005 (orange). The dotted lines refer to the merger rate of triples with a CHE inner binary where the tertiary has not influenced the evolution of the inner binary.
In the text

Fig. 11. Dominant type of drag force (gravitational versus hydrodynamic) as a function of the semi-major axis and the mass ratio for a secondary with M2 = 1 M⊙. The red dashed line indicates the radius of the primary (R1) at the ZAMS. The dotted and dash-dotted lines refer to the Roche radii of the primary at a semi-major axis of 3 R⊙ and 8 R⊙, respectively. If the radius exceeds the size of the Roche lobe, the system should not exist.

In the text