© The Authors 2024

1. Introduction

Mergers of non-compact stars frequently occur in the Universe (Podsiadlowski et al. 1992; Sana et al. 2012; de Mink et al. 2014). They can be caused by the evolution of the components in a binary system, during which the stars come into contact because of their radial expansion or orbital decay. Such orbital decay is not necessarily a result of the binary evolution itself, but can also be induced by von Zeipel–Kozai–Lidov oscillations (von Zeipel 1910; Lidov 1962; Kozai 1962; Naoz 2016) caused by a third component or a circumbinary disk (e.g. Lubow & Artymowicz 2000; Perets & Fabrycky 2009; Toonen et al. 2020, 2022). A third option is dynamically driven mergers, which occur during close encounters between stars in dense stellar environments (e.g. Hills & Day 1976; Portegies Zwart et al. 1997, 1999, 2004).

The products of stellar mergers can explain a multitude of objects. Examples include blue stragglers (e.g. Rasio 1995; Sills et al. 1997, 2001; Mapelli et al. 2006; Glebbeek et al. 2008; Ferraro et al. 2012; Schneider et al. 2015), some of the most massive stars observed in the Universe (e.g. Portegies Zwart et al. 1999; Banerjee et al. 2012; Schneider et al. 2014; Boekholt et al. 2018), B[e] supergiants (e.g. Podsiadlowski 2006; Wu et al. 2020), OBA stars with large-scale surface magnetic fields (Schneider et al. 2019), highly magnetic white dwarfs and magnetars (e.g. Tout et al. 2004; Ferrario & Wickramasinghe 2005; Wickramasinghe & Ferrario 2005; Wickramasinghe et al. 2014; Shenar et al. 2023), and α-rich young stars (e.g. Chiappini et al. 2015; Martig et al. 2015; Izzard et al. 2018; Hekker & Johnson 2019). Transients linked to stellar mergers include supernovae such as the core-collapse supernova SN 1987A (e.g. Podsiadlowski et al. 1990; Podsiadlowski 1992; Morris & Podsiadlowski 2007), the great eruption of η Car (Frew 2004; Gallagher et al. 1989; Iben et al. 1999; Podsiadlowski 2006, 2010; Morris & Podsiadlowski 2006; Fitzpatrick 2012; Portegies Zwart & van den Heuvel 2016; Smith et al. 2018; Owocki et al. 2019; Hirai et al. 2021), and luminous red novae such as V1309 Sco (Tylenda et al. 2011; Stȩpień 2011), V838 Mon (Soker et al. 2007), and V4332 Sgr (Tylenda et al. 2005).

Three-dimensional simulations of stellar mergers with magnetohydrodynamic (MHD, Schneider et al. 2019) and smoothed particle hydrodynamic (SPH, Sills et al. 1997, 2001; Freitag & Benz 2005; Dale & Davies 2006; Suzuki et al. 2007; Gaburov et al. 2008; Antonini et al. 2011; Glebbeek et al. 2013; Ballone et al. 2023) codes provide useful insights into the merger events and merger products (Schneider et al. 2020; Costa et al. 2022). However, to obtain realistic initial conditions for these computationally expensive simulations, it is crucial to understand in which binary configurations stellar mergers are most commonly expected to occur. Moreover, it is important to characterise the interior structures of stars directly before they enter a merger phase, as these largely determine the merger outcome.

Before stars in a binary system merge, they go through a phase of contact (e.g. Langer 2012). Yet, not every contact phase necessarily leads to a merger. These contact phases can be (over-)contact binaries in which both stars (over-)fill their Roche lobe¹, or classical common-envelope (CE) phases, in which one star is engulfed by the envelope of the other star (e.g. Ivanova et al. 2013). Hence, as a first step towards predicting the occurrence of mergers, it is important to determine in which binary configurations contact phases occur and what the outcomes of these phases are. Mapping the occurrence of contact binaries using detailed 1D binary evolution simulations has been carried out, for example, by Pols (1994), Wellstein et al. (2001), de Mink et al. (2007), Claeys et al. (2011), and Mennekens & Vanbeveren (2017) for massive stars. Marchant et al. (2016) and Menon et al. (2021) additionally computed through contact binary phases, which yields information on their lifetime and stability (i.e. likelihood to merge). More recently, three extended grids of detailed binary evolution models, spanning masses of 0.5 − 300 M_⊙, have been computed as part of the binary population synthesis code POSYDON (Fragos et al. 2023). These contain valuable information about the onset of contact phases. Using rapid binary population synthesis codes, de Mink et al. (2014) and Schneider et al. (2014) evolved entire populations of binary systems and mapped the occurrence of contact phases.

In this work, we use a grid of low-mass and massive binary evolution models with component masses between 0.5 and 20.0 M_⊙ at solar metallicity (Z = 0.0142, Asplund et al. 2009) to trace the occurrence of contact phases over the complete range of initial mass ratios² and orbital separations. We evolve both components and include physical processes such as stellar winds, rotation, and tidal interactions. It is known that binary systems can evolve towards tidal instabilities, which lead to rapid orbital decay and subsequent mergers (Darwin 1879). These instabilities have been proposed to be responsible for the lack of observations of W UMa type contact binaries at low mass ratios (Rasio 1995) and the final spiral-in of the progenitor system of V1309 Sco (Stȩpień 2011).

Including these physical processes allows us to arrive at a picture of the physical mechanisms leading to contact and their likelihood to lead to stellar mergers that is as complete as possible. Moreover, it illustrates the relative importance of, for example, tides, wind-mass loss, and mass-transfer efficiency on the evolution of binaries. We allow for non-conservative mass transfer and expect differences in the occurrence rate of contact binaries compared to works that assume conservative mass transfer (e.g. Pols 1994; Wellstein et al. 2001; Menon et al. 2021) or that use fixed mass-transfer efficiencies (e.g. de Mink et al. 2007; Claeys et al. 2011). Lastly, by including low-mass and massive binaries in our grid, we can compare the onset of contact over a wide mass range.

This paper is structured as follows. In Sect. 2, we describe the computational setup of the grid of binary evolution models. Section 3 covers the physical mechanisms that have been identified to lead to contact phases, their likelihood to result in a merger, and the way in which they are traced throughout the evolution of the binary models. Our findings of the occurrence of contact phases and mergers over the whole mass range, as well as some notable cases, are given in Sect. 4. Our results are discussed in Sect. 5, and summary and conclusions can be found in Sect. 6.

2. Methods

We computed a grid of 5957 1D binary evolution models using the binary module of MESA (release 12778; Paxton et al. 2011, 2013, 2015, 2018, 2019). First, we describe the adopted single-star physics used in the stellar models in Sect. 2.1 before describing the binary star physics in Sect. 2.2. We briefly outline the setup of the grid in Sect. 2.3 and list the stopping conditions of our binary evolution models in Sect. 2.4. In Sect. 2.5, we describe how we compute the birth probabilities for the binaries in our grid, which we use for population studies.

2.1. Adopted stellar physics

Each binary component was initialised from a precomputed zero-age main-sequence (ZAMS) model at solar metallicity, that is, Z = 0.0142 and Y = 0.2703 (Asplund et al. 2009). We used a blend of the OPAL (Iglesias & Rogers 1993, 1996) and Ferguson et al. (2005) opacity tables appropriate for the chemical composition of Asplund et al. (2009). We allowed the stars to rotate using the shellular approximation as implemented in MESA, with a limit on the rotation rate at 97% of the Roche critical rotation rate ${\mathrm{\Omega }}_{\mathrm{c}}=\sqrt{GM/R_{\mathrm{eq}}^3}\simeq \sqrt{8GM/27R^3}$ (Maeder 2009). In this expression, M and R are the mass and radius of the star, respectively, G is the gravitational constant and R_eq is the equatorial radius of a rotationally deformed star. All models were hydrostatic, meaning that MESA’s implicit hydrodynamic solver was disabled. We used the approx21 nuclear network.

2.1.1. Mixing

Convective mixing was handled via the mixing length theory (Böhm-Vitense 1958; Cox & Giuli 1968) and a mixing length parameter of α_mlt = 2.0 (Paxton et al. 2013). We used the Ledoux criterion to identify regions in the star that were unstable to convection. The efficiency of semi-convective mixing was set to α_sc = 10.0 (Schootemeijer et al. 2019). Thermohaline mixing was also included with an efficiency of α_th = 1.0 (Marchant et al. 2021).

Convective boundary mixing (CBM) was included via the step-overshoot scheme, in which we allowed the convective hydrogen-burning core to extend by 0.20 H_P beyond the core boundary set by the Ledoux criterion, with H_P the pressure scale height (Martinet et al. 2021). At the bottom of nonburning convective envelopes, we used step overshoot with a 0.05 H_P extension of the convective region towards the centre. This corresponds to one-half of the upper limit typically inferred for the Sun (Angelou et al. 2020). For convectively burning cores beyond the main sequence (MS), overshooting is not yet understood properly and is known to have a large effect on the final fate of stars (Herwig 2000; Temaj et al. 2024). Because of this, we followed Marchant et al. (2021) and used exponential overshoot extending only 0.005 H_P beyond the edge of the convective region set by the Ledoux criterion.

To account for the rotational mixing of chemical elements and diffusion of angular momentum, the Goldreich–Schubert–Fricke instability, Eddington-Sweet circulation, and the secular and dynamic shear instabilities were included (see for example Heger et al. 2000 for a detailed description). Additional diffusion of angular momentum was taken care of by the Spruit–Tayler dynamo. We scaled the strength of the mixing of chemical elements by a factor f_c = 1/30 as in Heger et al. (2000). The sensitivity of the rotational instabilities to stabilising composition gradients, which is incorporated in the factor f_μ, was set to f_μ = 0.1 (see also Pinsonneault et al. 1989).

2.1.2. Wind mass-loss prescription

Our binary models contain low-, intermediate-, and high-mass stars, and the masses can change significantly over the star’s evolution. Therefore, we employed a wind mass-loss prescription that covered this large mass range and all evolutionary stages.

We considered two distinct regimes, the hot-wind regime with surface temperature³T_surf ≥ 11 kK and the cool-wind regime for stars with T_surf ≤ 10 kK. We used linear interpolation to determine the mass-loss rate in the temperature region between those two values.

Within the hot-wind regime, the mass-loss rate for stars with hydrogen envelopes (surface hydrogen mass fraction X_surf > 0.5) was computed via the Vink et al. (2000) prescription. When X_surf dropped below 0.4, either the prescription for Wolf-Rayet (WR) stars from Sander & Vink (2020) or the prescription for low-mass helium stars from Vink (2017) was used, depending on whether the star’s luminosity L was higher or lower than a certain luminosity L₀ respectively. As described in Sander & Vink (2020), L₀ is the asymptotic limit below which no WR-like wind mass-loss is expected to occur. Its value is metallicity-dependent and obtained from stellar atmosphere models. In the regime where L > L₀, we computed the mass-loss rate with both prescriptions and took the maximum value as the adopted wind mass-loss rate (Vink, priv. comm.). All of the aforementioned prescriptions had a scaling factor of 1.0. For 0.4 ≤ X_surf ≤ 0.5, the mass-loss rate was determined via linear interpolation between the mass-loss rates from both regimes.

When a model reached the cool-wind regime, the distinction between stars expected to become giants or supergiants was made. The cut was made at log₁₀(ℒ/ℒ_⊙) = 3.15, where ℒ is defined as in Langer & Kudritzki (2014),

$$\begin{array}{c}\hfill \mathcal{L}=\frac{1}{4\pi \sigma G}\frac{L}{M}.\end{array}$$(1)

Here, σ is the Stefan-Boltzmann constant. This cut corresponds roughly to a mass of 10 M_⊙ at the base of the (super-)giant branch. Models below the cut used the Reimers (1975) wind prescription on the red giant branch (RGB) and the Bloecker (1995) wind prescription on the asymptotic giant branch (AGB). Following Choi et al. (2016), we used a scaling factor of 0.1 for the former and 0.2 for the latter. Models above the cut in ℒ used the Nieuwenhuijzen & de Jager (1990) prescription with a scaling factor of 1.0.

We increased the scaling factor for the Bloecker (1995) wind to 3.0 at the onset of thermal pulses (TP) during the AGB phase following Choi et al. (2016). This increase aimed to ease the computations through this phase by mimicking the enhanced mass loss during the TP-AGB phase while simultaneously avoiding the TPs themselves (by removing the envelope). Additionally, by removing part of the envelope, we aimed to avoid the Hydrogen Recombination Instability (HRI) and the Fe-Peak Instability (FePI). These instabilities, which lead to envelope inflation over multiple orders of magnitude, can occur when the cold, expanded envelopes of AGB stars are modelled with a hydrostatic code (Rees et al., in prep.). The HRI is caused by the increased dynamical instability of the envelope because of hydrogen recombination (Wagenhuber & Weiss 1994), and the FePI occurs when luminosity at the base of the convective envelope exceeds the Eddington luminosity because of a local iron opacity bump (Lau et al. 2012). The physical mechanism behind these instabilities most likely leads to events of extreme mass loss. The timescales on which this envelope inflation occurs are too short to be captured correctly in MESA’s hydrostatic mode and can lead to numerical issues. The increase of the Bloecker (1995) wind scaling factor was not successful in avoiding numerical issues in each model, especially in those models where the aforementioned instabilities occur. Because of this, we opted to disregard models in which binary mass transfer occurs after the TP-AGB phase. A more elaborate approach to compute through the TP-AGB phase and these instabilities is provided in, for example, Rees et al. (in prep.).

2.2. Adopted binary physics

In our models, we evolved both binary components. This allowed us to consider the behaviour of both the donor and the accretor star for tracing potential contact scenarios (see Sect. 3). Only mass transfer from the initially more massive or ‘primary’ star onto the initially less massive or ‘secondary’ star was considered. Hence, in this work, references to primary (secondary) and donor (accretor) are equivalent. We mostly employ the names primary (subscript ‘1’) and secondary (subscript ‘2’) star. We assumed circular orbits for all binary systems.

2.2.1. Mass transfer and accretion

Whenever the primary star was on the MS, mass transfer was computed using the contact scheme (Marchant et al. 2016). In semi-detached binaries, this scheme uses MESA’s roche_lobe scheme, which ensures that the donor stays within its Roche lobe. When both stars (over-)fill their Roche lobe it switches to a different solver for the mass-transfer rate suitable for contact binaries. In this scheme, only the computation of the mass-transfer rate is handled. Energy transport between the components of the contact binary and the tidal distortion are not taken into account (for the effect of including those, see Fabry et al. 2022, 2023). For systems with post-MS primary masses smaller than 1.3 M_⊙, the Kolb scheme (Kolb & Ritter 1990) was used because this scheme is better suited for envelopes with larger pressure scale heights (Fragos et al. 2023). Both schemes were solved implicitly.

The mass transfer efficiency β is defined as the effective, overall change in mass of the accretor over the mass transferred from the donor to the accretor, β ≡ −Ṁ₂/Ṁ_trans. In this definition, Ṁ_trans < 0 and Ṁ₂ > 0. Ṁ₂ also includes the wind mass loss of the accretor. During conservative mass transfer in our models, β ≈ 1, because typically |Ṁ_trans| is orders of magnitude larger than the absolute value of the wind mass-loss rate. When the accretor star spun up to its critical rotation rate Ω_c,  2, mass transfer was non-conservative and β < 1. When the accretor star’s accretion timescale τ_acc ≡ Ṁ₂/Ṁ_trans approached the star’s dynamical timescale τ_dyn,  2, defined for a star with sound speed c_S as τ_dyn,  2 ≡ R₂/c_S (Kippenhahn et al. 2012), we limited the accretion rate to that of 0.1 times the star’s Kelvin-Helmholtz (or thermal) timescale τ_KH. This also resulted in non-conservative mass transfer, that is, β < 1. We did this because the accreting models tend not to converge numerically when τ_acc ∼ τ_dyn,  2. Models for which the accretion rate was limited are marked in the results. The accretion of angular momentum during mass transfer was computed following Lubow & Shu (1975) and Ulrich & Burger (1976), which includes accretion through ballistic impact and a Keplerian disk.

2.2.2. Tides and angular momentum loss

Tidal synchronisation was computed uniformly over the components’ structure using the convective synchronisation timescale from Hurley et al. (2002). Upon initialisation of the binary models, the rotation periods of both components were equal to the orbital period. In our models, orbital angular momentum evolved via mass loss from the system (through winds or isotropic re-emission) and spin-orbit coupling. In the former case, the lost mass was assumed to have the specific angular momentum of the star’s orbit in which vicinity it was leaving the system.

2.3. Binary-star grid

In our grid, we varied the initial mass of the primary star M_1,  i (in units of M_⊙), the initial mass ratio q_i = M_2,  i/M_1,  i with M_2,  i the initial mass of the secondary star, and the initial separation a_i (in units of R_⊙). Twenty primary masses were selected between 0.8 and 20.0 M_⊙ with logarithmic mass spacing Δlog₁₀M_1,  i ≈ 0.074. At the high-mass end, we added three additional primary masses (13.1, 15.6 and 18.4 M_⊙) for increased resolution, bringing the total number of primary masses to 23. Values between 0.1 and 0.9 were considered for the mass ratios, with linear spacing and steps of 0.1. We imposed a lower limit on the secondary mass of 0.5 M_⊙ to avoid fully convective companions, which are difficult to converge numerically when they accrete. We added an additional mass ratio at q_i = 0.97 to model twin systems. Lastly, we chose the initial binary separations a_i such that the first phase of mass transfer occurred during all stages of the primary star’s evolution. These stages at which mass transfer occurs are called Case A, B, and C, for when the donor star is on the MS, before core-helium ignition, and after core-helium ignition, respectively⁴. This distinction is related to the phases in which the (primary) star expands, as illustrated in Fig. 1. Case B phases were further divided into early (Case Be) and late (Case Bl) phases. The distinction is based on the presence of a deep convective zone in the envelope of the star⁵.

Fig. 1. Radial expansion of a 10.2 M⊙ star (solid black line). The radial evolution is divided into different cases based on the main phases of expansion (solid horizontal lines). Each dashed horizontal line indicates a sample point at which we put R = RRL.

To ensure sufficient sampling of all these stages, we used the radius evolution of single-star models in combination with the volume-equivalent Roche lobe radius R_RL for a given mass ratio q = M₂/M₁ and separation a, (Eggleton 1983),

$$\begin{array}{c}\hfill \frac{R_{\mathrm{RL}}}{a}=\frac{0.49q^{-2/3}}{0.6q^{-2/3}+ln(1+q^{-1/3})}.\end{array}$$(2)

For a given primary mass M₁ and mass ratio q, we selected a number of points along the radius evolution of the star, indicated by the dashed horizontal lines in Fig. 1. At each of these points, we assumed the star filled its Roche lobe, R = R_RL, and mass transfer ensued. From Eq. (2), we then got the separation a at which this happened. Assuming that the orbit stays approximately constant until mass transfer starts, we took this value for the separation a as our initial value, a_i = a. Although this method of sampling the initial separation space is adequate for closer binaries, it cannot accurately predict the division lines between later mass-transfer cases for a number of reasons. Firstly, the assumption that the orbital separation a stays approximately constant does not hold in systems with strong wind mass loss and tidal interaction prior to mass transfer. Secondly, post-MS donors use MESA’s Kolb scheme for mass transfer (Kolb & Ritter 1990). This scheme can result in significant mass-transfer rates before the donor star formally overfills its Roche lobe, since it takes both the optically thick and thin Roche lobe overflow into account (Paxton et al. 2015).

2.4. Stopping conditions

In principle, we computed each binary model until one of the components, usually the initially more massive primary star, reached the end of core carbon burning. This is defined as the moment when the central mass fractions of He and C fall below 10⁻⁶ and 10⁻², respectively. Systems with MS primaries having predominantly radiative envelopes (M_1,  i > 1.3 M_⊙) used MESA’s contact mass-transfer scheme (see Sect. 2.2.1) and could therefore be modelled through such phases. During a contact phase, the computation terminated when the condition for mass overflow through the second Lagrange point (L₂) or L₂-overflow from Marchant et al. (2016) was met (see Sect. 3.3). For primaries with mostly convective envelopes (post-MS or M_1,  i ≤ 1.3 M_⊙), we used the Kolb mass-transfer scheme. The computation stopped whenever the secondary star overfilled its Roche lobe by more than R_RL,  2. When the mass-transfer rate Ṁ_trans reached a value of 10 M_⊙ yr⁻¹, the evolution was terminated. The timescale on which such fast mass transfer occurs nears typical dynamical timescales, which is not numerically feasible with our MESA setup. Lastly, reverse mass transfer in a semi-detached system from the secondary to the primary star was not considered in this study. The onset of reverse mass transfer was used as a stopping condition.

2.5. Birth probabilities

To describe the models in our grid as a population (Sects. 4.3–4.5), we computed their birth probability p_birth from the initial distribution function, which is the product of the distribution functions of the initial primary mass M_1,  i, mass ratio q_i and period P_i. For the M_1,  i distribution, we used the Kroupa (2001) initial mass function (IMF),

$$\begin{array}{c}\hfill \psi (M_{1,\mathrm{i}})=C_m{M}^{\alpha },\end{array}$$(3)

with α = −2.3 as appropriate for primary-star masses ≥0.5 M_⊙ and C_m a normalisation constant. We assumed the other two initial distribution functions to be uniform in q_i and log₁₀P_i, respectively,

$$\begin{array}{c}\hfill \varphi (q_{\mathrm{i}})=C_q\text{and}\chi ({log}_{10}{P}_{\mathrm{i}})=C_p,\end{array}$$(4)

with C_q and C_p being normalisation constants. The birth probability p_birth of each model was then computed by integrating the product of the initial distribution functions over the parameter size in the model grid,

$$\begin{array}{cc}\hfill p_{\mathrm{birth}}=& {\displaystyle {\int }_{{log}_{10}{P}_{\mathrm{l}}}^{{log}_{10}{P}_{\mathrm{u}}}{\int }_{q_{\mathrm{l}}}^{q_{\mathrm{u}}}{\int }_{M_{\mathrm{l}}}^{M_{\mathrm{u}}}\psi (M_{1,\mathrm{i}})}\hfill \\ \hfill & \times \varphi (q_{\mathrm{i}})\chi ({log}_{10}{P}_{\mathrm{i}})\mathrm{d}{M}_{1,\mathrm{i}}\mathrm{d}{q}_{\mathrm{i}}\mathrm{d}{log}_{10}{P}_{\mathrm{i}}.\hfill \end{array}$$(5)

Here, log₁₀P_u,  l, q_u,  l and M_u,  l are the upper- and lower-boundaries of the parameter size of a model in the grid for log₁₀P_i, q_i, and M_1,  i, respectively. They were chosen as the midpoints between the initial parameter values of the model and its neighbouring models.

2.6. Software for analysis and plotting

For the analysis in this work, we used the following open-source software: PyMesaReader (Wolf & Schwab 2017⁶), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), MPI for Python (Dalcín et al. 2005; Dalcin & Fang 2021), and Astropy (Astropy Collaboration 2013, 2018, 2022). For constructing the figures, we used Matplotlib (Hunter 2007).

3. Physical mechanisms leading to contact

We use the following nomenclature of contact phases traced in our grid, based on the physical picture in Röpke & De Marco (2023) and illustrated in Fig. 2. Systems with MS and/or Hertzsprung-gap (HG) stars (i.e. Case A & Be) that both (over-)fill their Roche lobes or undergo (unstable) runaway mass transfer are referred to as contact binaries (Figs. 2.1b and 2.2b). Systems with a (super-)giant primary (i.e. Case Bl and C) and an MS secondary undergoing (unstable) runaway mass transfer enter a ‘classical’ common-envelope phase (Fig. 2.3b). What sets the onset of a classical CE phase apart from the formation of a contact binary is that now the primary has a clear core-envelope boundary and the radius of the secondary is at least an order of magnitude smaller than that of the primary.

Fig. 2. Schematic representation of the physical mechanisms leading to contact (not to scale). The filled blue circles and purple squares indicate the position of the L1- and L2-points, respectively. The filled grey-blue circles represent the stellar cores in panels 3a and 3b. Panel 1a shows the expansion of the accretor leading to a contact binary (1b). This corresponds to the ‘Accretor expansion’ mechanism described in Sect. 3.1. Subsequent overfilling of the components’ Roche lobes leads to the formation of an overcontact binary (2b). The primary increasingly overfills its Roche lobe (2a) and can eventually fill the secondary’s Roche lobe (2b). This can eventually lead to L2-overflow (2c), which likely results in a stellar merger. The scenarios (1b–2b–2c) and (2a–2b–2c) correspond to the ‘L2-overflow’ mechanism described in Sect. 3.3. In panel 3a, runaway mass transfer from a (super-)giant (left) to an MS star (right) leads to the onset of a classical common-envelope phase 3b, where the cores of both stars revolve in the (super-)giant’s envelope. This corresponds to the ‘Runaway MT’ mechanism described in Sect. 3.5.

In this section, we discuss the different physical mechanisms that can lead to contact phases. As stated before, being in contact does not necessarily mean that the binary will merge. Hence, for each mechanism leading to contact, we discuss the likelihood of a merger or rather a longer-lived contact phase. Additionally, we explain how we trace each mechanism in the binary star models.

3.1. Expansion of accretor

Arguably the most intuitive way to form contact binaries is when the accreting secondary star expands during a mass-transfer phase and (over-)fills its Roche lobe. A schematic representation of this mechanism is shown in Figs. 2.1a–b. The timescale on which the secondary expands, generally defined as τ_R/Ṙ ≡ R/Ṙ with R the stellar radius and Ṙ its time derivative, has implications for the evolution of the subsequent contact phase.

When the mass-transfer timescale τ_trans ≡ |M₁/Ṁ_trans| is shorter than the secondary’s thermal or Kelvin-Helmholtz timescale⁷ defined as (Kippenhahn et al. 2012)

$$\begin{array}{c}\hfill {\tau }_{\mathrm{KH}}=\frac{G{M}^2}{2RL}\approx 1.5\times {10}^7{\left(\frac{M}{M_{\odot }}\right)}^2\frac{R_{\odot }}{R}\frac{L_{\odot }}{L}\mathrm{yr},\end{array}$$(6)

the star is out of thermal equilibrium. In an effort to regain thermal equilibrium, the secondary expands on its thermal timescale (Ulrich & Burger 1976; Kippenhahn & Meyer-Hofmeister 1977; Webbink 1976; Neo et al. 1977; Pols 1994). A contact binary is formed if the increase in radius is sufficient for the secondary to fill its Roche lobe. Such a contact binary can rapidly merge on a thermal timescale, or evolve back to a semi-detached binary when the accretor regains thermal equilibrium and shrinks.

When the secondary is in thermal equilibrium during mass transfer, it expands on a nuclear timescale, defined as (Kippenhahn et al. 2012)

$$\begin{array}{c}\hfill {\tau }_{\mathrm{nuc}}=\frac{E_{\mathrm{nuc}}}{L}\approx {10}^{10}\frac{M}{M_{\odot }}\frac{L_{\odot }}{L}\mathrm{yr}.\end{array}$$(7)

Here, E_nuc is the available nuclear energy. Contact phases driven by the nuclear expansion of the accretor are longer lived since the accretor will not shrink inside its Roche lobe again until the end of the nuclear burning phase. As a result, such contact phases are expected to persist on a nuclear timescale, or until the stars merge.

The onset of contact phases through the expansion of the accretor is traced by checking whether at any point the accretor overfills its Roche lobe. These models are labelled as ‘Accretor expansion’ in the following sections. This is a stopping condition for models using the Kolb mass-transfer scheme. Models using the contact mass-transfer scheme can evolve through these contact phases. Hence, multiple contact phases might occur throughout the evolution of the system. If more than one contact phase occurs in a model using the contact scheme, only the onset of the first contact phase is considered for a better comparison with models using the Kolb scheme.

An example of a system forming a contact binary because of the thermal expansion of the accreting secondary star (shown up to the point of contact) is shown in Fig. 3. From the Hertzsprung–Russell diagram (HRD, Fig. 3a), it can be seen that, as a result of Case-A mass transfer, the primary and secondary become under- and over-luminous compared to their single-star counterparts, respectively. The single-star tracks are shown up to their terminal-age main sequence (TAMS). As a result of mass accretion, the secondary star expands rapidly and fills its Roche lobe (Fig. 3b). The expansion timescale τ_R/Ṙ,2 becomes about an order of magnitude shorter than the secondary’s global thermal timescale τ_KH,  2 due to the mass-transfer timescale τ_trans becoming comparable to τ_KH,  2 (Fig. 3c). An example of a system forming a contact binary because of the nuclear expansion of the accretor is provided in Appendix A. The accretor expansion timescales of all systems forming contact binaries through accretor expansion are shown in Figs. D.1–D.2.

Fig. 3. Example of a M1,  i = 10.2 M⊙, qi = 0.4, and ai = 12.4 R⊙ binary system forming a contact binary during the MS of the primary (Case A) because of the thermal expansion of the secondary (accretor) star. Panel a shows an HRD with the evolutionary tracks of the primary and secondary stars (solid lines). The dashed black lines show the evolutionary tracks of single stars with the same initial masses as the binary components and initial rotation rates of Ω/Ωc = 0.25. Panel b shows the evolution of the radius R (solid lines) and Roche lobe radius RRL (dashed lines) of both components. The timescales governing the evolution of the binary and the mass-transfer rate (black dashed line) are shown in panel c as a function of the decreasing primary star mass. The secondary’s nuclear (τnuc,  2) and thermal (τKH,  2) timescales are shown with a dashed green and pink line, respectively. The expansion (τR/Ṙ,2) and mass-transfer (τtrans) timescales are shown with a solid blue and orange line, respectively.

3.2. Non-conservative mass transfer

As demonstrated in Packet (1981) and, for example, more recently in Ghodla et al. (2023), a star only has to accrete between ∼2 and ∼10% of its own mass to be spun up to the Roche critical rotation rate Ω_c. When a star is rotating at (near) this critical rotation rate, it may not be able to accrete any more material. It is then often assumed that the excess mass is lost through an enhanced stellar wind, or even instantaneously through a so-called ‘fast’ wind from the accretor, resulting in non-conservative mass transfer. In our models, instantaneous wind mass loss is invoked to expel the excess matter for accretors rotating near the critical rotation rate and for models in which the accretion timescale τ_acc is restricted to 10% of the Kelvin-Helmholtz timescale τ_KH. It is assumed that there is an energy source that can expel this excess mass. Let us consider the energy per unit mass ε required to drive away the mass to infinity, ε = GM₂/2R₂. This is under the simplifying assumption that the primary’s gravitational potential can be ignored and that the mass is lost from the surface of the accretor. Equating ε to (L₁ + L₂/Ṁ_max), we eventually find a maximum mass-loss rate Ṁ_max the accretor can have under the condition that all this mass is driven away to infinity by the combined luminosity (L₁+L₂) of the binary (Marchant 2017),

$$\begin{array}{c}\hfill \frac{{\dot{M}}_{\max }}{M_{\odot }{\mathrm{yr}}^{-1}}={10}^{-7.19}{f}_{\mathrm{eff}}\frac{R_2}{R_{\odot }}\frac{M_{\odot }}{M_2}\frac{(L_1+L_2)}{L_{\odot }}.\end{array}$$(8)

The factor f_eff is a free parameter added to take into account the uncertainty on the exact radius at which the non-accreted mass is expelled, the fact that only the gravitational potential of the accretor is accounted for, the unknown fraction of the total luminosity that can be used to expel the non-accreted mass, and the unknown energy of the ejected mass at infinity.

Since the non-accreted mass in systems with Ṁ₂ > Ṁ_max cannot be expelled to infinity, it must remain in or around the binary system. Although it is unclear what exactly happens, the non-accreted mass can potentially lead to a contact phase, for example, if it fills the secondary’s Roche lobe or introduces a drag on the binary components. Alternatively, the non-accreted matter can form an accretion disk. In this work, we merely flag such models and discuss the potential consequences for their evolution in Sect. 5.3.

In the models, we trace the failure to eject non-accreted matter in post-processing by using Eq. (8) and assuming f_eff = 1.0. Models for which the mass-loss rate of the accretor, defined as Ṁ_ej = (1 − β)Ṁ_trans, at one point exceeds Ṁ_max are labelled in the following sections as ‘Non-conservative MT + cannot eject’.

An example of a system in which the accretor cannot eject the non-accreted matter during non-conservative mass transfer is shown in Fig. 4. The system undergoes Case-Be and -C mass transfer. The primary is partially stripped in the former mass-transfer phase, and the remaining H+He envelope amounts to ≈20% of the mass of the post-mass-transfer primary star. After core-He exhaustion, the star expands again, initiating a Case-C mass-transfer phase. During Case-Be mass transfer, the secondary star becomes over-luminous because it is out of thermal equilibrium. However, its luminosity decreases again when its rotation rate Ω reaches the critical rotation rate Ω_c (Fig. 4c). At this point (M₁ ≈ 9.5 M_⊙), the mass-transfer efficiency β decreases to almost zero. During the non-conservative Case-Be mass transfer, the combined luminosity is insufficient to expel the non-accreted matter to infinity. This holds both for f_eff = 0.1 and f_eff = 1.0 (Fig. 4b).

Fig. 4. Example of a M1,  i = 10.2 M⊙, qi = 0.9, and ai = 39.2 R⊙ binary star model in which non-conservative mass transfer occurs and for which the non-accreted matter cannot be driven to infinity. Panel a is the same as Fig. 3a (‘ign.’ = ‘ignition’ and ‘exh.’ = ‘exhaustion’). In panel b, the mass-loss rate of the secondary (solid black line) is compared to the maximum mass-loss rate Ṁmax (dashed lines; pink for feff = 1.0, green for feff = 0.1) set by Eq. (8). The surface rotation rates (solid blue and orange lines) and mass transfer efficiency (dashed black line) are shown in panel c as a function of the primary mass.

3.3. Outflow from second Lagrange point

In the Roche potential, there are three equilibrium points located on the line connecting the centres of the binary components. From lowest to highest potential, these are L₁ (through which mass transfer occurs), L₂ and L₃. The L₂-point is always located on the side of the less massive star in the binary, and L₃ on the side of the more massive star. In Case-A and -Be binaries, mass loss from the L₂-point (Fig. 2.2c) takes away a significant amount of specific angular momentum from the system, which can lead to rapid orbital shrinkage and a subsequent merger. The rate of orbital shrinkage and hence the time until the merging event depends on the mass-loss rate and the outflow velocity through L₂ (Marchant et al. 2021). In Case-Bl and -C binaries, L₂ indicates the onset of a classical CE.

We trace L₂-overflow in our models when either one (semi-detached) or both (contact) components overfill(s) their Roche lobe. In the latter case, the radii of the components are compared to L₂ volume-equivalent radii from Marchant et al. (2016). This is done during the computation of the model and the moment of L₂-overflow is a stopping condition (Sect. 2.4). In semi-detached systems, we compare the radius of the primary star to the volume-equivalent radius R_L2 and the distance to the L₂-point D_L2 from Misra et al. (2020; similar fitting formulae for R_L2 are derived in Ge et al. 2020b). This is done in post-processing, meaning that the outflow from L₂ is not modelled and does not affect the evolution of the model. We refer to, for example, Marchant et al. (2021) to see the effect of these outflows. In both cases, if L₂-overflow occurs, the models are labelled ‘L₂-overflow’.

The evolution of a binary system experiencing L₂-overflow is shown in Fig. 5. The primary and secondary become under- and over-luminous, respectively, during Case-Be mass transfer (Fig. 5a). As in the previous example, mass transfer is non-conservative and early on during the mass transfer the non-accreted matter cannot be driven to infinity (Fig. 5b). As the binary orbit shrinks during mass transfer, also the Roche lobe radius R_RL,  1 and the volume-equivalent radius R_L2,  1 decrease (Fig. 5c). After about 1.5 M_⊙ is lost from the primary star, its radius decreases slower than R_RL,  1 and R_L2,  1. As a consequence, the star increasingly overfills its Roche lobe until it reaches the point of L₂-overflow. The subsequent loss of mass and angular momentum from the L₂-point will result in an accelerated orbital shrinkage, which is not accounted for in this model. The onset of L₂-overflow is shown schematically in Fig. 2. In the evolutionary scenario depicted in Figs. 2.1a–b and 2.2b–c, a contact binary is formed because of the expansion of the accretor. The components continuously overfill their Roche lobes (Fig. 2.2b) until the overcontact binary fills the L₂-lobe (Fig. 2.2c). In the second scenario, from Figs. 2.2a–c, the primary overfills its own Roche lobe (Fig. 2.2a) and later fills the Roche lobe of the secondary, forming an overcontact binary (Fig. 2.2b). Eventually, this overcontact binary fills the L₂-lobe and L₂-overflow occurs (Fig. 2.2c).

Fig. 5. Example of a binary model with M1,  i = 10.2 M⊙, qi = 0.2, and ai = 175.7 R⊙ going through a phase of (delayed) runaway mass transfer. Panel a is the same as Fig. 3a (‘ign.’ = ‘ignition’). Panel b shows the mass-transfer rate Ṁtrans (solid black line), thermal mass-transfer rate ṀKH (solid green line) and dynamical mass-transfer rate Ṁdyn (solid gold line) for the primary star on the left axis. The right axis shows the component mass evolution as a function of age (dashed blue and orange lines). In panel c, the radius R (solid blue line), Roche lobe radius RRL (dashed light-blue line), L2-volume-equivalent radius RL2 (dahsed green line, Misra et al. 2020, Eq. (15)) and orbital separation aorb (dashed black line) evolution for the primary are shown as a function of the decreasing primary mass. The blue cross indicates the moment of L2-overflow.

3.4. Tidally driven contact

When a single star ascends the (super-)giant branch, its moment of inertia I increases, such that its surface rotation velocity decreases. In a binary and if the tidal synchronisation timescale τ_sync (Hurley et al. 2002, Eq. (27)) is shorter than the expansion timescale τ_R/Ṙ, the star is tidally locked to the orbit and does not spin down. The transfer of angular momentum from the orbit to the spin of the star shrinks the orbit. Under certain circumstances, tidally driven orbital decay can lead to contact phases, which demonstrates the importance of considering the effect of tides in binaries.

In the most extreme case, the binary system becomes Darwin unstable (Darwin 1879). This instability arises if

$$\begin{array}{c}\hfill \frac{L_{\mathrm{orb}}}{(I_1+I_2)\mathrm{\Omega }}<3,\end{array}$$(9)

where L_orb is the orbital angular momentum, Ω is the orbital angular rotation velocity, and I₁ and I₂ are the moment of inertia of the primary and secondary star, respectively. This condition is derived under the assumption of a circular (e = 0, with e the eccentricity), coplanar (spins of components are aligned), and synchronised (ω₁ = ω₂ = Ω) system. Here, ω₁ and ω₂ are the angular rotation rates of the primary and secondary, respectively. Moreover, solid-body rotation is assumed for the individual components, allowing the spin angular momenta S_1,  2 to be written as S_1,  2 = ω_1,  2I_1,  2. Especially for close or contact binaries, where τ_sync is typically of the order of a few years, the Darwin instability can lead to a dynamical inspiral of the components, resulting in a merger.

Models labelled ‘Tidally driven contact’ experience orbital decay caused by tides before the onset of contact, while orbital widening is expected from mass transfer (q > 1). We trace this condition in post-processing. This label is not used for systems which are at one point during their evolution Darwin unstable. Although these systems will experience orbital decay, this can happen on timescales longer than the evolutionary timescales of the system. However, most of the systems that are driven into contact by tides do eventually become Darwin unstable.

An example of a system in which tides lead to the onset of contact is shown in Fig. 6. In this low-mass binary, the primary overfills its Roche lobe near the end of the MS. During mass transfer, the primary star reaches the TAMS, but does not become a red giant because of the continuous stripping of its envelope. During the mass-transfer phase, the secondary star overtakes in evolution, leaves the MS and turns to the giant branch. At this point, the mass ratio of the system has already reversed, and mass transfer is almost conservative, widening the orbit. Two important changes occur when the secondary becomes a giant. First, there is the development of a deep convective envelope and an increase in radius, which lead to increased tidal coupling of the star: τ_sync,  2 drops by several orders of magnitude (Fig. 6b). Because of this, the orbital and rotation period of the secondary synchronise. Secondly, the increase in radius and density redistribution (the envelope now has a deep convective zone) both increase the moment of inertia I₂ of the star (Fig. 6b). Following the conservation of angular momentum, the secondary spins down, yet it is immediately spun up again by tides. Hence, the spin angular momentum of the secondary S₂ increases and the orbital angular momentum L_orb decreases (Fig. 6a). Eventually, the system becomes Darwin unstable, as it fulfils the condition⁸L_orb/(S₁ + S₂)≲3. Because of the decrease in L_orb, the orbit, which was widening from mass transfer, starts shrinking again. As a consequence of the shrinking Roche lobe, the primary star is stripped of its outer envelope increasingly rapidly, leading to a decrease in radius and an increase in mass-transfer rate (not shown). The shrinking of the Roche lobes eventually leads to the formation of a contact binary once the secondary also fills its Roche lobe (Fig. 6c). At the point when the contact binary is formed, the helium core mass M_He has increased beyond that of the primary star. This means that the secondary is closer to core-helium ignition, that is, it has overtaken the primary in evolution.

Fig. 6. Example of a M1,  i = 1.1 M⊙, qi = 0.8 and ai = 4.1 R⊙ binary system in which tides lead to contact. In Panel a we show the exchange between the orbital angular momentum Lorb (solid black line) and the secondary’s spin angular momentum S2 (solid orange line), which coincides with the decrease in Lorb/(S1 + S2) (dashed green line). The primary’s spin angular momentum S1 is shown with a solid blue line. Panel b shows the evolution of the moment of inertia (solid lines) and tidal synchronisation timescale (dashed lines) for both components. While the moment of inertia of the primary decreases around M1 = 0.3 M⊙, it sharply increases for the secondary. At the same time, the secondary star’s tidal synchronisation timescale decreases by approximately two orders of magnitude. Panel c shows the evolution of the primary’s and secondary’s radii (solid lines) and Roche lobe radii (dashed lines). While the former fills its Roche lobe, tides cause orbital shrinkage, which results in the secondary also filling its Roche lobe.

3.5. Unstable mass transfer and CE phases

During mass transfer, the donor star’s radius and Roche lobe can both either shrink or expand. Mass transfer is generally stable when the donor star shrinks faster or expands slower than its Roche lobe. However, when the donor star overfills its Roche lobe by increasing amounts as a reaction to mass loss, the binary enters an unstable runaway situation (e.g. Soberman et al. 1997). This happens, for example, when the Roche lobe shrinks while the donor expands (e.g. during mass transfer with q < 1 and/or because of orbital angular momentum loss). The responses of the donor’s radius and Roche lobe radius are quantified in the mass–radius exponents ξ_R1 ≡ dlnR₁/dlnM₁ and ξ_RL ≡ dlnR_RL,  1/dlnM₁, respectively (Webbink 1984). If ξ_R1 < ξ_RL, mass transfer is unstable, and the mass-transfer rate keeps increasing, potentially reaching values larger than 1 M_⊙ yr⁻¹. During runaway (unstable) mass transfer, the primary star increasingly overfills its Roche lobe, and the secondary’s expansion timescale becomes orders of magnitude lower than its thermal timescale. When the primary star is an MS or HG star (Case-A or -Be), a contact binary forms and merge on a timescale shorter than the primary’s thermal timescale because of the runaway expansion of both components. In Case-Bl and -C binaries, runaway mass transfer leads to the engulfment of the secondary in the primary’s envelope, that is, a classical CE phase. A classical CE phase can lead to a merger or leave a close binary behind (see, e.g. Röpke & De Marco 2023). The onset of a classical CE through runaway mass transfer is shown schematically in Figs. 2.3a,b.

We look for runaway mass transfer in post-processing and label such models ‘Runaway MT’ given the following three criteria. Firstly, the mass-transfer rate Ṁ_trans needs to exceed the thermal-timescale mass-transfer rate, set by Ṁ_KH = M₁/τ_KH,1. Secondly, the second time derivative of log₁₀Ṁ_trans has to be positive. When this is the case, the rate of change in Ṁ_trans is increasing, which is indicative of a runaway situation. Moreover, this condition also ensures that Ṁ_trans does not decrease again, as is observed for stable Case-C mass transfer (see Sect. 4.2). Thirdly, the condition for unstable mass transfer described above needs to be fulfilled: ξ_R1 < ξ_RL. For semi-detached models using MESA’s contact mass-transfer scheme, ξ_R1 = ξ_RL by definition. Hence, in this case, only the first two conditions are evaluated. We note that L₂-overflow is not a necessary condition for the onset of runaway mass transfer. However, it does often occur in systems with runaway mass transfer.

An example of a binary system experiencing delayed runaway mass transfer is shown in Fig. 5 and was previously discussed to demonstrate L₂-overflow in Sect. 3.3. The evolution of the mass-transfer rate Ṁ_trans as a function of the binary system’s age is shown in Fig. 5b, where it is also compared to the thermal and dynamical mass-transfer rates. Shortly after Ṁ_trans exceeds the former, the runaway nature of the evolution is observed, during which Ṁ_trans starts nearing the dynamical mass-transfer rate. Simultaneously, we see in Fig. 5c that the slope of R₁ as a function of the decreasing M₁ becomes larger than that of R_RL,  1, which means that the primary increasingly overfills its Roche lobe. This continues up to the point where L₂-overflow occurs, after which an accelerated orbital shrinkage is expected (see Sect. 3.3). In this example, the binary does not enter a runaway phase immediately at the onset of mass transfer. Therefore, it is an example of a delayed runaway mass-transfer phase (e.g. Hjellming & Webbink 1987; Han & Podsiadlowski 2006; Pavlovskii & Ivanova 2015; Ge et al. 2015, 2020a). Since in this particular example the primary star is not yet a supergiant at the onset of mass transfer, this system is not expected to enter a classical common-envelope phase, but rather become a contact binary that will most likely eventually merge.

4. Occurrence of contact phases

In this section, we present our contact tracing results for three initial primary masses (Sect. 4.1). Then, we focus on a particular region of the initial binary parameter space in which Case-Bl and -C mass transfer are found to be stable (Sect. 4.2). Next, we look at the incidence of the different physical mechanisms leading to contact in a population of binary systems (Sect. 4.3), and put lower limits on the stellar merger and classical CE fractions of mass-transferring binaries (Sect. 4.4). Lastly, we compare the properties of Case-A contact systems found in our grid with those of observed systems (Sect. 4.5).

4.1. Contact phases for different initial primary masses

4.1.1. Initial primary mass of 10.2 M_⊙

In Fig. 7, we show our binary models in the initial mass ratio–separation (q_i − log₁₀a_i) plane for a fixed initial primary mass M_1,  i of 10.2 M_⊙. Systems with M_1,  i = 20.0 and 1.6 M_⊙ are shown in Fig. 8, and described in Sects. 4.1.2 and 4.1.3, respectively. The results for the other M_1,  i are shown in Appendix B. A table with the evolutionary outcomes of all computed MESA models as well as the model input and output files are available online⁹. An extract of this table can be found in Appendix G.

Fig. 7. Occurrence of contact phases for models with initial primary masses M1,  i = 10.2 M⊙ on the initial mass ratio–separation plane. Models marked with a dot are those for which accretion was limited to 0.1 times the secondary’s global thermal timescale τKH (see Sect. 2.2.1). The other ones are marked with a star symbol. The dark-blue quasi-horizontal lines indicate the initial mass transfer cases, which can be read from the right side. Systems on the left of the solid black line are Darwin unstable at the onset of mass transfer according to Eq. (9) and assuming R1 = RRL,  1.

Fig. 8. Same as Fig. 7, but for M1,  i = 20.0 M⊙ (left) and M1,  i = 1.6 M⊙ (right).

The coloured regions in Fig. 7 indicate the different physical mechanisms leading to contact (‘Accretor expansion’, ‘Non-conservative MT + cannot eject’, ‘L₂-overflow’, ‘Tidally driven contact’, and ‘Runaway MT’, see Sect. 3). They have been constructed using nearest neighbour interpolation. Models labelled ‘No contact’ undergo at least one phase of mass transfer but manage to avoid any form of contact until the end of our computations (see Sect. 2.4). In some models, reverse mass transfer occurs from the initially less massive secondary to the primary, but we do not consider this here for contact tracing. These models are also labelled ‘No contact’. The potential fates of reverse mass transfer systems are briefly described later in this section. Models labelled ‘Numerical issues’ were numerically unable to reach the desired stopping conditions of our computations (see Sect. 2.4). Lastly, models labelled ‘MT after TP’ are those for which Case-C mass transfer occurs after the TP-AGB phase. Their final fate is not further interpreted (see Sect. 2.1.2).

In some regions of the initial binary parameter space, more than one of the above conditions are met, and we indicate this by the corresponding hatching (ancillary outcome). For the background colour (principal outcome), priority is given in the following order (highest to lowest priority): ‘Tidally driven contact’, ‘Runaway MT’, ‘Accretor expansion’, ‘L₂-overflow’ ‘No contact’, and ‘Numerical issues’. Regions in which ‘No contact’ is indicated by hatching contain models that make it to the end of core-He burning and are, therefore, strong candidates for avoiding contact but encountered numerical difficulties before they could reach the end of core-C burning (see Sect. 2.4). We use ‘Non-conservative MT + cannot eject’ exclusively as an ancillary outcome (hatching). Although this is a viable mechanism leading to contact (Sect. 3.2), it is uncertain to what evolutionary outcome it leads. For example, a system with ‘No contact’ as its principal outcome and ‘Non-conservative MT + cannot eject’ as its ancillary outcome can either avoid contact when, for example, an accretion or circumbinary disk is formed, or form a contact binary when the non-accreted matter fills the secondary’s Roche lobe.

For systems with M_1,  i = 10.2 M_⊙, contact binaries form because of the expansion of the accreting secondary star for q_i = 0.15 − 1.00 in the initially closest Case-A binary systems (‘Accretor expansion’; Fig. 7). For q_i = 0.15 − 0.45 with log₁₀(a_i/R_⊙) ≲ 1.23 and q_i = 0.75 − 1.00, these contact systems additionally experience L₂-overflow, which leads to orbital angular momentum loss and subsequent stellar mergers. For q_i = 0.45 − 0.75, the contact binaries do not expand beyond the L₂-lobe. This behaviour is found consistently throughout the range of M_1,  i = 2.6 − 20.0 M_⊙ for q_i between 0.45 and 0.65 − 0.75, where the upper boundary decreases with decreasing M_1,  i. By comparing with the expansion timescales of the accretors (Figs. D.1–D.2) in systems avoiding L₂-overflow, we find that systems with q_i ≳ 0.5 likely produce longer-lived (τ_nuc,  1) contact binaries, which may be observable. The contact binaries with q_i ≲ 0.5 form through the thermal-timescale expansion of the accretor and merge or detach again on the accretor’s thermal timescale.

The initially wider Case-A (log₁₀(a_i/R_⊙) ≳ 1.23) systems with q_i = 0.15 − 0.35 that avoid L₂-overflow reach mass-transfer rates Ṁ_trans orders of magnitude larger than their thermal mass-transfer rate Ṁ_KH. Ṁ_trans eventually reaches the stopping condition value of 10 M_⊙ yr⁻¹ (see Sect. 2.4). In the Case-Be models labelled ‘Accretor expansion’, the primary is an HG star and the secondary is expanding on a timescale orders of magnitude shorter than its thermal timescale. The resulting contact binary will likely experience L₂-overflow and merge.

The Case-A binary systems with q_i ≤ 0.15 and log₁₀(a_i/R_⊙) ≲ 1.37 experience a phase of runaway mass transfer (‘Runaway MT’). In addition, the mass-accretion rate of the secondaries is limited to 10 Ṁ_KH,2, which brings the mass-transfer efficiency β almost down to zero right after the onset of mass transfer. The lack of further accretion quenches initially rapid expansion (τ_R/Ṙ,2 <τ_KH,2) of the secondary. Assuming that in reality this rapid expansion continues (the accretion rate in these models is limited for numerical reasons), the two MS stars are likely to form contact binaries.

In the region of q_i = 0.35 − 1.00 and log₁₀(a_i/R_⊙) ≈ 1.2 − 1.5, we find Case-A systems that avoid contact (‘No contact’). In these systems, the primary stars are stripped in Case-A, Case-AB, and in some models even Case-ABC mass-transfer phases. Eventually, they reach or are expected to reach (horizontal hatching) the end of core carbon burning, or enter a phase of reverse mass transfer. The latter occurs for some of the binary systems with q_i = 0.65 − 1.00 close to the border between forming contact binaries and avoiding contact. Here, the secondary (over-)fills its Roche lobe after the primary has detached, leading to a phase of accretion onto a stripped primary star. Other models that experience reverse mass-transfer phases are found at q_i = 0.97 and log₁₀(a_i/R_⊙) ≈ 1.52 − 2.30. Such phases are not considered further in our contact tracing.

For Case-Be systems, contact is likely avoided for q_i ≳ 0.25 − 0.35. The primary stars in these systems all follow the same evolutionary pathway as the example described in Fig. 4: the HG primaries are stripped and reach core-C exhaustion. The secondary stars are spun up by accretion and mass transfer becomes non-conservative. Tidal synchronisation timescales are longer than the spin-up timescales, hence tides are not able to prevent spin-up to the critical rotation rate as is the case in the closer-orbit Case-A systems. Virtually all systems fail to eject the non-accreted matter at one point during this phase of non-conservative mass transfer (‘Non-conservative MT + cannot eject’).

At q_i ≤ 0.25 − 0.35, Case-Be binaries experience runaway mass transfer. All systems experience non-conservative mass transfer and fail to eject the non-accreted matter. Except for systems at log₁₀(a_i/R_⊙) ≲ 1.78 and q_i = 0.15 − 0.35, all of them also experience L₂-overflow. The secondary stars are still on the MS. The primaries are HG stars without a clear core-envelope boundary (see Figs. C.1–C.3 for the evolutionary state of all models at contact and/or termination). Hence, these systems are expected to evolve into contact binaries and not result in classical CEs.

We find runaway mass transfer in Case-Bl systems with q_i ≤ 0.25 and log₁₀(a_i/R_⊙) ≲ 2.92, and q_i ≤ 0.35 and log₁₀(a_i/R_⊙) ≳ 2.92 as well as a small set of Case-C systems with q_i = 0.15 − 0.35. The secondary stars in these systems are MS stars. The difference compared to the Case-Be systems experiencing runaway mass transfer is that all primary stars have now evolved into supergiants with a clear core-envelope boundary by the time mass transfer starts (this can be observed from, for example, the steeper drop in density and binding energy around the core-envelope boundary). As a consequence, all these systems are expected to enter a classical CE phase.

Following Rasio (1995), we have computed the mass ratios below which binary systems are Darwin unstable at the onset of mass transfer (black solid line in Fig. 7). Systems with supergiant donors with deep convective envelopes (Case Bl and C) are Darwin unstable at the onset of mass transfer for q_i ≲ 0.1. In Case-A and -Be binaries, this is the case for q_i ≲ 0.05.

At q_i = 0.15 − 0.45, we find Case-C systems (and one Case-Bl system) with runaway mass-transfer where the primary stars overfill their L₂-lobes by a factor of more than one. The primaries have (almost) engulfed their companion and we witness the beginning of a classical CE phase.

For q_i ≥ 0.45, Case-Bl and -C binary systems evolve through a stable phase of mass transfer. The Case-Bl systems are actually found to go through two stable mass-transfer phases, Case-Bl and Case-BC. Thanks to the stability of these mass-transfer phases, contact is avoided. The stability of Case-Bl and -C mass transfer is described further in Sect. 4.2.

The region around the transition between early and late Case-B systems contains models for which MESA does not converge numerically because of the rapid spin-up of the accretor. These issues with the spin-up of the accretor occur in systems with initial primary masses down to 8.6 M_⊙. For systems with initial primary masses of ≥13.1 M_⊙, models with Case-Be mass transfer are also affected. This can be seen in Figs. 8 and B.1.

4.1.2. Initial primary mass of 20.0 M_⊙

For M_1,  i = 20.0 M_⊙ (left panel in Fig. 8), the general picture is similar to that of models with M_1,  i = 10.2 M_⊙ presented above. However, there are some notable differences.

Similarly to the models with M_1,  i = 10.2 M_⊙, the initially closest Case-A systems form contact binaries through the expansion of the secondary star during mass transfer. For the twin systems (q_i ≈ 0.97), this even happens for all the Case-A systems. Binary systems with q_i ≤ 0.15 and log₁₀(a_i/R_⊙) ≲ 1.37 enter contact through accretor expansion, whereas binary systems in the equivalent region for M_1,  i = 10.2 M_⊙ do so through a phase of runaway mass transfer. In general, it is found that the former holds for models with M_1,  i = 12.6 − 20.0 M_⊙ and the latter for M_1,  i = 5.2 − 10.2 M_⊙. For M_1,  i < 5.0 M_⊙, models with q_i = 0.1 are not computed, so information about this region is not available.

Just as for the systems with M_1,  i = 10.2 M_⊙, contact binaries with q_i ≤ 0.45 formed through accretor expansion experience L₂-overflow (Fig. 8, left). However, at q_i ≥ 0.75, L₂-overflow is found to be largely avoided during the first phase of contact, contrary to what is found for systems with M_1,  i = 10.2 M_⊙. This, however, does not imply that these contact phases are long-lived (τ_contact ∼ τ_nuc,  1) since the secondary has been found to shrink again after regaining thermal equilibrium (τ_contact ∼ τ_KH,  2). In almost all of these systems, a second contact phase follows later in the evolution because of the nuclear-timescale expansion of the secondary star. The primary stars are at his point either MS or post-MS stars. In the former case, the models have again been computed through the contact phase, which almost exclusively results in L₂-overflow. In the latter case, the evolution is terminated when the accretor fills its Roche lobe (see Sect. 2.4).

Initially wider Case-A systems, with log₁₀(a_i/R_⊙) = 1.37 − 1.68 and q_i ≤ 0.35, go through a phase of runaway mass transfer and thus likely form contact binaries. Compared to the systems with M_1,  i = 10.2 M_⊙, this region extends to higher values of q_i. At similar initial separations and q_i = 0.35 − 0.55, MESA does not converge numerically. As in the equivalent region for M_1,  i = 10.2 M_⊙ models, the solver fails to find a suitable solution for the accretor such that the rotation rate remains below the critical rotation rate.

Primary stars in Case-A systems with q_i = 0.550 − 0.935 and log₁₀(a_i/R_⊙) = 1.38 − 1.77 are stripped in Case-A, Case-AB, and Case-ABC mass-transfer phases and avoid or are expected to avoid contact. Contrarily to similar systems with M_1,  i = 10.2 M_⊙, none experience a phase of reverse mass transfer. There is, however, a set of models of twin systems with q_i = 0.97 and log₁₀(a_i/R_⊙) ≈ 1.78 − 2.80 where reverse mass transfer does occur, as in the binaries with M_1,  i = 10.2 M_⊙.

In Case-Be binaries, the differences with respect to the equivalent M_1,  i = 10.2 M_⊙ systems are more pronounced. Firstly, a few systems only experience L₂-overflow and thus avoid runaway mass transfer. Closer inspection shows that in these systems, the radius of the primary star exceeds R_L2 by ≲20% for ≲10³ yr. It is uncertain whether these contact binaries (see Fig. 2.2c) are stable, given that the phase of L₂-overflow is short and might not lead to sufficient angular momentum loss to cause significant orbital decay and hence a merger. Moreover, the primary star shrinks again rapidly afterwards, causing the binary to become semi-detached again. Secondly, Case-Be models have a harder time converging numerically. In these models, the primary stars are stripped in one or more mass-transfer phase(s), and continue core-He burning and core-C burning as stripped stars with a thin hydrogen layer (≲1 M_⊙). At these masses, MESA runs into numerical difficulties and the timesteps of the simulations drop well below one year. This is a known issue for this kind of stripped stars (Götberg, priv. comm.) and has prevented us from computing the evolution to the end of core-C exhaustion for all models. A select number of models are computed with small timesteps until core-C exhaustion. In these models, the hydrogen surface layers expand during core-C burning and drive a stable and short-lived (∼10³ yr) Case-C mass-transfer phase. Mass transfer is stable because the donor stars shrink again when nearing core-C exhaustion and the core-mass fraction is > 0.5, which typically signals stability (Temmink et al. 2023). At lower initial primary masses, such as for models with M_1,  i = 18.4 M_⊙ (see Fig. B.1), the aforementioned numerical difficulties are less severe, allowing most stripped primaries to reach core-C exhaustion. Most of them avoid Case-C mass transfer, such that we do not expect contact phases.

From the left panel of Fig. 8, we see that for Case-Be, -Bl, and -C systems with q_i ≤ 0.25 and log₁₀(a_i/R_⊙) ≲ 2.5, and q_i ≤ 0.15 and log₁₀(a_i/R_⊙) = 2.50 − 3.42, contact is reached again through runaway mass transfer. Among these, Case-Bl and Case-C systems are likely to enter classical CE phases. We also see two models at q_i = 0.2 with log₁₀(a_i/R_⊙) ≈ 1.72 and log₁₀(a_i/R_⊙) ≈ 2.57, respectively, that enter contact through accretor expansion. However, the accretor expansion is driven by numerical difficulties with finding accretors that spin below critical. Contact appears likely in these two models, albeit for another reason (unstable mass transfer or L₂-overflow).

Just as for the M_1,  i = 10.2 M_⊙ systems, there is a small region in the initial binary parameter space at q_i ≲ 0.35 and log₁₀(a_i/R_⊙) ≈ 3.4 − 3.5 where runaway mass transfer does not occur, but a classical CE phase is expected from L₂-overflow. For q_i ≳ 0.35, a similar region of stable Case-Bl and Case-C mass transfer as for the M_1,  i = 10.2 M_⊙ systems is found (see Sect. 4.2).

In the same way as the 10.2 M_⊙ initial primary mass models, the spin-up of the accretor leads to numerical convergence problems in both the Case Bl and -Be regions. At these higher masses, the problems persist down to initial separations of a few 100 R_⊙.

4.1.3. Initial primary mass of 1.6 M_⊙

At lower initial primary masses, more specifically M_1,  i = 1.6 M_⊙, the situation is different than at higher masses (right panel in Fig. 8)¹⁰. Although the primary stars have radiative envelopes during the MS, none of our systems undergo Case-Be mass transfer, because of the resolution in a_i¹¹.

Case-A models at q_i = 0.35 − 0.65 encounter numerical issues when the accreting secondary star reaches critical rotation. Also in Case-Bl systems with mass ratios q < 0.65 − 0.75, the models do not converge numerically. As can be seen from Figs. B.2–B.3, this is an issue that occurs for M_1,  i ≤ 2.2 M_⊙. Although this prevents us from including binary systems with M_1,  i ≤ 2.2 M_⊙ in further population analysis, specific regions of the initial binary parameter space can still be described and hold valuable information.

The expansion of the accretor star is leading to the formation of contact binaries in Case-A systems with log₁₀(a_i/R_⊙) ≲ 0.57 and q_i = 0.35 − 0.55, and log₁₀(a_i/R_⊙) = 0.57 − 0.79 and q_i = 0.55 − 1.00 (Fig. 8). Similar to contact systems with M_1,  i = 10.2 M_⊙, systems with q_i ≳ 0.80 experience L₂-overflow.

Case-A systems with q_i = 0.65 − 0.94 and log₁₀(a_i/R_⊙) = 0.73 − 0.94 evolve through similar pathways as those in equivalent regions of the initial parameters space with M_1,  i = 10.2 M_⊙ and M_1,  i = 20.0 M_⊙. The primaries are stripped in Case-A and Case-B mass-transfer phases, after which reverse mass transfer occurs.

In between the Case-A systems at q_i = 0.65 − 1.00 that form contact binaries through accretor expansion and those that avoid contact, contact binaries form by tidal interaction. In our grid, such tidally driven contact in Case-A systems can be found for initial primary masses between 0.8 and 1.8 M_⊙. For initial primary masses between 0.8 and 1.1 M_⊙ tidally driven contact also occurs in the initially closest Case-B systems with q_i = 0.97 (Fig. B.3).

Case-Bl binary systems at q_i = 0.65 − 1.00 and log₁₀(a_i/R_⊙) = 0.91 − 1.45 evolve into contact binaries because of the expansion of the accretor. The donor stars in these binaries have deep convective envelopes, resulting in mass-transfer rates ≳ 10⁻⁵ M_⊙ yr⁻¹. Because of the relatively high mass transfer rates, the accretors are out of thermal equilibrium and expand with τ_dyn,2 < τ_R/Ṙ,2 < τ_KH,2 and fill their Roche lobes. In comparison, the Case-A binaries at log₁₀(a_i/R_⊙) < 0.91, which avoid contact, have mass transfer rates ≲ 10⁻⁷ M_⊙ yr⁻¹, resulting in a nuclear-timescale expansion of the accretor star. The three models at q_i = 0.85 − 1.00 not labelled ‘Non-conservative mass transfer + cannot eject’ have conservative mass transfer because they manage to reach contact before the secondary star rotates at its critical rotation rate.

At larger initial separations (log₁₀(a_i/R_⊙) = 1.07 − 2.64), Case-Bl systems with q_i > 0.65 − 0.75 avoid contact. The primary stars have their envelopes removed during the red giant (RG) phase and end up as (partially) stripped stars. Mass transfer is stable because of the relatively high initial mass ratios (the mass ratio becomes greater than one early on during mass transfer). The mass transfer stops when the entire hydrogen envelope is removed. A pure helium WD remains if the hydrogen-burning shell is stripped before the helium core grows to a mass above roughly 0.45 M_⊙. If this is not the case, central helium ignition occurs and the primary ends up as a carbon-oxygen white dwarf. Systems with log₁₀(a_i/R_⊙) = 2.42,  2.27,  2.28,  2.12 do not ignite helium in the centre, whereas the initially wider systems do. We find reverse mass transfer in all cases once the secondary star becomes an RG and fills its Roche lobe. Given that the primary stars are WDs, these systems might be observable as symbiotic binaries.

Case-Bl systems with q_i = 0.35 − 0.45 and log₁₀(a_i/R_⊙) = 2.38 − 2.56 go through a phase of runaway mass transfer and experience L₂-overflow. Since the primary star is an RG with a deep convective envelope and the secondary star an MS star at the onset of mass transfer, we expect a classical CE phase. Initially slightly wider binary systems (q_i = 0.35 − 0.65 and log₁₀(a_i/R_⊙) ≈ 2.60) avoid runaway Case-Bl mass transfer. After core-He exhaustion, Case-BC mass transfer starts and the model quickly fails to converge numerically. While there is no clear indication for future runaway mass transfer, we cannot rule it out. Binary systems with q_i = 0.45 − 0.65 and log₁₀(a_i/R_⊙) = 2.39 − 2.60 only go through Case-Bl mass transfer. The models stop when the primary stars evolve into a pure helium WD, and the secondaries climb the giant branch.

Finally, we find that in the initially widest Case-C systems, contact phases are avoided. Here, mass transfer is stable because the primary star is in its TP-AGB phase, during which it experiences enhanced mass loss in our models (Sect. 2.1.2). Moreover, the star’s radius periodically decreases again, preventing a runaway situation. Models that failed to converge numerically experience Case-C mass transfer before or early on during the TP-AGB phase of the primary. Their outcomes are uncertain.

4.2. Stable Case-Bl and -C mass transfer

For all initial primary masses considered in this work, Case-C mass transfer has been found to be stable over a wide range of initial mass ratios, even down to q_i ≈ 0.2 in some cases (e.g. for M_1,  i = 15.5 M_⊙, see Fig. B.1). This also applies to the initially widest Case-Bl systems, but the exact extent of this region in the initial binary parameter space is unknown because of the aforementioned numerical difficulties.

The response of the donor star’s radius to mass loss has been studied using models with polytropic equations of state in Hjellming & Webbink (1987) and detailed adiabatic mass loss computations such as in Ge et al. (2010, 2015, 2020a). It is found that donors with radiative or convective envelopes with core-mass fractions¹² greater than 0.5 shrink and hence have stable mass transfer (Temmink et al. 2023)¹³. Donors with convective envelopes and core-mass fractions below 0.5 typically expand in response to mass loss and cause unstable mass transfer.

In our models, we identify two stabilising effects during Case-Bl and -C mass transfer. First, the primary star’s envelope is partially lost due to winds prior to mass transfer. This increases the core-mass fraction. At core-helium ignition, values of the core-mass fraction in our stable Case-C models are 0.22, 0.15, 0.19, 0.26, and 0.31 for M_1,  i = 1.9,  4.4,  8.6,  14.3, and 20.0 M_⊙, respectively. At core-helium exhaustion, the values for the core-mass fraction are 0.26, 0.22, 0.28, 0.36, and 0.44, respectively. Only for the most massive primary stars with M_1,  i = 20.0 M_⊙ at core-helium exhaustion, the core-mass fraction approaches the stabilising value of 0.5. Hence, the increased core-mass fraction alone is insufficient to explain the mass transfer stability. However, stellar wind mass loss increases the mass ratio q towards unity before the onset of mass transfer. With q’s approaching or exceeding one, orbits will (soon) widen, stabilising mass transfer.

Second, orbital widening can stabilise mass transfer. Because the mass transfer in these systems is non-conservative, it results in less orbital shrinkage ȧ/a than in the conservative case. Furthermore, this also causes orbital widening already at mass ratios q < 1, whereas this only occurs for q ≥ 1 in conservative mass transfer (Tauris & van den Heuvel 2006). How orbital widening stabilises mass transfer can be seen in the example of a M_1,  i = 10.2 M_⊙, q_i = 0.6, and a_i = 1398.2 R_⊙ system undergoing stable Case-C mass transfer in Fig. 9. The orbital separation drops sharply before the onset of mass transfer because of an enhanced spin-orbit coupling when the primary becomes a supergiant. After the primary loses about 0.1 M_⊙, the mass-transfer efficiency β becomes zero when the secondary star reaches critical rotation. The orbit starts to widen when the mass ratio reaches a value of 0.67 and the primary star’s Roche lobe radius stays nearly constant, preventing a situation of runaway mass transfer as described in Sect. 3.5. Further on, the primary star’s Roche lobe radius increases faster than the star’s radius, such that R₁ < R_RL,  1 at q = 1.21. At M₁ = 7.2 M_⊙ (q = 0.87, age = 25.8815 Myr) the mass-transfer rate reaches its peak, after which it decreases sharply. This model ends its evolution when the primary reaches core-C exhaustion, which happens during mass transfer as in most of our stable Case-C models.

Fig. 9. Example of a M1,  i = 10.16 M⊙, qi = 0.6, and ai = 1398.2 R⊙ binary model with stable Case-C mass transfer. Panel a is the same as Fig. 3a (‘ign.’ = ‘ignition’). Panel b is the same as Fig. 5b. Panel c is the same as Fig. 5c, with the addition of a solid orange line indicating the evolution of the mass ratio.

4.3. Population properties

To understand what fraction of systems avoid or end up in contact and through which mechanism, we count the systems in each category (‘No contact’, ‘Accretor expansion’ etc.) and weigh them with their birth probability p_birth (see Sect. 2.5). We show the results for initial primary mass ranges of [4.8; 12.6) M_⊙ and [12.6; 20.8] M_⊙ in so-called sunburst charts in Fig. 10. Figure 11 shows the results for the whole mass range of M_1,  i  ∈  [4.8;20.8] M_⊙¹⁴. We only consider binaries with M_1,  i ≥ 4.8 M_⊙ because models with lower masses have not been computed for the entire q_i = 0.1 − 0.97 range (Sect. 2.3). The sunburst charts in the top row show the incidences of the physical mechanisms leading to contact for mass-transferring binaries. The inner level of these charts shows the fraction of all mass-transferring binary systems in the specified mass range that end their evolution with the indicated outcome. In other words, they represent the principal outcomes, ‘Accretor expansion’, ‘Runaway MT’, ‘L₂-overflow’, and ‘No contact’ (see Sect. 4.1.1), which correspond to the outcomes indicated by the background colours in Figs. 7 and 8. Models expected to avoid contact but which encountered numerical issues before the primaries reached core-C exhaustion (‘No contact’ hatching in Figs. 7 and 8) are included in the principal ‘No contact’ outcomes. The outer level of the sunburst charts shows the incidence of ancillary outcomes, that is, the hatching in Figs. 7 and 8 (‘L₂-overflow’ and ‘Non-conservative MT + cannot eject’). We assign systems that experienced numerical issues a likely evolutionary outcome based on their initial mass ratio and first mass-transfer case. For more information, we refer to Appendix E. These outcomes are shown at the intermediate level of the charts in the ‘Numerical issues’ slice. On the bottom row of Fig. 10, we show the incidences of the principal outcomes per initial mass-transfer case. The outer level of these charts shows the lower limits of the incidences for stellar mergers and classical CEs (see Sect. 4.4).

Fig. 10. Sunburst charts displaying the fractions of evolutionary outcomes for mass-transferring binary systems in the grid over initial primary mass ranges of [4.8; 12.6) M⊙ (left) and [12.6; 20.8] M⊙ (right). Wedges with a percentage < 1% are not labelled. Top row: the inner level shows the principal outcome of the evolution (‘Accretor expansion’, ‘Runaway MT’, ‘L2-overflow’, and ‘No contact’), while the outer level shows the ancillary outcome (‘L2-overflow’, ‘Non-conservative MT + cannot eject’). Ancillary ‘No contact’ outcomes are incorporated in the inner ‘No contact’ category. Models in the ‘Numerical issues’ category are assigned a likely evolutionary outcome based on their initial mass ratio qi and mass-transfer case, displayed on the sunburst chart’s intermediate level. The labels ‘A’, ‘Be’, and ‘Bl/C’ refer to Case-A, Case-Be, and Case-Bl or -C mass transfer, respectively. Bottom row: the inner level shows the percentage of Case-A, -Be, -Bl, and -C systems. We show the principal outcome per case on the middle level. The outer level shows the lower limits of the total fraction of stellar mergers and classical CEs per mass-transfer case.

We find that the fractions for M_1,  i  ∈  [4.8;12.6) M_⊙ and M_1,  i  ∈  [12.6;20.8] M_⊙ are relatively similar (Fig. 10). The most noticeable differences are found for the systems going through a phase of runaway mass transfer (‘Runaway MT’) and systems failing to eject non-accreted matter (‘Non-conservative MT + cannot eject’). The lower fraction of systems that simultaneously avoid contact and fail to eject non-accreted matter for higher primary masses can be traced back to the mass–luminosity relation. In general, L ∼ M^α, where the exponent α > 1 (Kippenhahn et al. 2012). So, even though the gravitational potential increases linearly with increasing mass, the luminosity increase is steeper since α > 1. Hence, it is easier for higher mass systems to expel non-accreted matter.

The fraction of systems with runaway mass transfer in the lower mass range is 5% higher than in the higher mass range. At the same time, the fraction of systems with numerical issues is 7% higher in the higher mass range. However, the total fraction of systems with numerical issues that have runaway mass transfer as their most likely outcome is 4 − 6% in both mass ranges. The higher fraction of systems with runaway mass transfer in the lower mass range is thus not caused by increased numerical difficulties at higher masses but is physical.

Excluding the systems with numerical issues, we find that ≥41% of binaries with M_1,  i ∈ [4.8;12.6) M_⊙ and ≥34% of binaries with M_1,  i  ∈  [12.6;20.8] M_⊙ enter a contact phase (Fig. 10). Over the whole initial primary mass range of M_1,  i ∈ [4.8;20.8] M_⊙, the percentage of binaries entering a contact phase is ≥40% (Fig. 11). These are lower limits because, in addition to excluding the systems with numerical issues, we do not take into account the binaries that avoid contact in our models but fail to eject non-accreted matter.

Fig. 11. Sunburst charts displaying the fractions of evolutionary outcomes for mass-transferring binary systems in the grid over initial primary mass ranges of [4.8; 20.8] M⊙. The left and right charts are equivalent to those in the top and bottom row of Fig. 10, respectively.

4.4. Stellar merger and classical CE incidence

We assume that all binaries experiencing runaway mass transfer and/or L₂-overflow merge or enter a classical CE phase. Following the physical picture described in Sect. 3, we make the distinction based on the structure of the primary star. This means that Case-A and -Be binaries with runaway mass transfer and/or L₂-overflow lead to stellar mergers, and Case-Bl and -C binaries to classical CEs. Based on this, we compute lower limits on the incidences of stellar mergers and classical CEs, shown in the bottom row of Fig. 10 for initial primary mass ranges of [4.8; 12.6) M_⊙ and [12.6; 20.8] M_⊙, and for each initial primary mass in Table 1. This table also lists the critical mass ratios q_crit for which binaries with q_i < q_crit merge or enter classical CE phases. The stellar merger and classical CE incidences are ≥13% and ≥21%, respectively, for the primary mass range of [4.8; 12.6) M_⊙, and ≥12% and ≥11%, respectively, for the primary mass range of [12.6; 20.8] M_⊙. Figure 11 shows similar charts for the total initial primary mass range of [4.8; 20.8] M_⊙, and we find that ≥12% of mass-transferring binaries merge and ≥19% evolve towards a classical CE phase. The stellar merger incidence for Case-A binaries of 8% (Fig. 11) is similar to the incidence of massive stars with strong, large-scale magnetic fields of ∼10% (Donati & Landstreet 2009; Fossati et al. 2015; Grunhut et al. 2017), which are likely formed by (pre-)MS stellar mergers (Schneider et al. 2019).

The incidences reported in Figs. 10 and 11 and Table 1 are lower limits since the criteria given here do not take into account binary systems that merge as a result of a classical CE phase, binary systems that experience L₂-overflow in later contact phases (so far if more than one contact phase occurred, we only considered the first one, see Sect. 3.1), and potential mergers among the models that had numerical issues. Furthermore, a certain fraction of models avoiding contact but failing to eject non-accreted matter (‘Non-conservative MT + cannot eject’) might in reality also merge (see discussion in Sect. 5.3). If we now assume that all contact binaries formed through the expansion of the accretor eventually merge (except those which survive contact and reach core-C exhaustion in either component), we find stellar merger incidences of ≥18%, ≥20%, and ≥19% for initial primary masses of [4.8;12.6) M_⊙, [12.6;20.8] M_⊙, and [4.8;20.8] M_⊙, respectively. For comparison, the merger incidence for O-type stars (M₁ ≳ 15 M_⊙) found by interpreting observations in the context of binary evolution and reported in Sana et al. (2012) is 20 − 30%.

We find that the incidences of stellar mergers of ∼12% are similar for all initial primary masses (bottom row of Fig. 10 and Table 1). The classical CE incidence varies more significantly, from 12% at M_1,  i = 5.2 M_⊙, to 33% at M_1,  i = 10.2 M_⊙ and 9% at M_1,  i = 20.0 M_⊙. The decrease in classical CE incidence from 10.2 M_⊙ to 20.0 M_⊙ is linked to the decrease in Case-C systems (bottom row of Fig. 10). Our values of q_crit = 0.15 − 0.35, correspond reasonably well with those for HG star donors from Temmink et al. (2023).

4.5. Comparison with observed contact binaries

In Fig. 12, our simulated population of Case-A contact binaries with initial primary masses of 4.8 − 20.8 M_⊙ is compared to observed near-contact¹⁵ and contact systems in the Milky Way (MW), Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC) from the compilation in Menon et al. (2021). The sample consists of MS O+O, B+O and B+B massive contact systems, which is why we have chosen to compare only to Case-A binaries from our grid with initial primary masses ≥4.8 M_⊙ (this is as low as we can go in terms of the initial primary mass while still having models at all q_i). The mass ratios of the observed systems are compared to the probability distribution function (PDF) of the models as a function of the mass ratio q at contact. For models that enter a contact phase through the expansion of the accretor, contact is defined as the moment when both the primary and the secondary overfill their respective Roche lobes. When runaway mass transfer is responsible for the onset of a contact phase, the moment at which Ṁ_trans > Ṁ_KH,1 is taken as the moment of contact. To have a meaningful comparison with observed contact systems, the mass ratios at contact of the observed systems are defined here as the mass of the currently less massive component over the mass of the currently more massive one (whereas before q has always been defined as the mass of the initially less massive component over the mass of the initially more massive one). The probabilities used to construct the PDF are the birth probabilities p_birth computed via Eq. (5) in Sect. 2.5. We highlight the following contributions to the PDF: ‘accretor expansion’ systems are those that enter contact because of accretor expansion but do not experience L₂-overflow, while the ‘accretor expansion + L₂-overflow’ do. The third contribution comes from runaway mass transfer systems (‘runaway MT’). The ‘accretor expansion’ and ‘accretor expansion + L₂-overflow’ systems are further divided into systems where the accretor expands on a nuclear and a thermal timescale prior to filling its Roche lobe. This division is made by comparing the mean of the expansion timescale τ_R/Ṙ from the onset of mass transfer to the formation of a contact system with the mean of the thermal and nuclear timescales respectively (see Appendix D).

Fig. 12. Probability density function (PDF) of contact systems formed in Case-A binaries with M1,  i = 7.9 − 20.8 M⊙ as a function of the mass ratio q at the onset of contact. Data points (filled squares) show the observed mass ratio and primary mass (right axis) for all observed MW, LMC and SMC contact and near-contact systems compiled in Menon et al. (2021), including the uncertainties on their values. Systems without reported uncertainties on q are indicated with a dash symbol, and those without reported uncertainties on q and M1 are indicated with a star symbol.

There is a striking difference between the PDF from the model systems and the observations. From the models, we expect two to three times more contact systems at q < 0.5 than at q ≥ 0.5. Contrarily, observations show a dearth of contact binaries at q < 0.5. There are only a few models at q ≥ 0.5 experiencing L₂-overflow, while such systems are dominant at q < 0.5. For q < 0.45, there is an additional contribution of contact systems formed through runaway mass transfer. Among the contact binaries that do not experience L₂-overflow or runaway mass transfer (‘accretor expansion’), it can be seen that at q < 0.5 contact is formed mostly on a thermal timescale, while for q ≥ 0.5 contact is formed on the longer, nuclear timescale. Overall, 29% of the Case-A contact binaries in the mass range considered here form through the nuclear-timescale expansion of the accretor. 59% of those (17% of the total) do not yet experience L₂-overflow and are expected to be longer-lived.

We find a number of ways in which this (apparent) discrepancy can be resolved. Firstly, it should be noted that the sample of observed systems is small, especially if one disregards the near-contact systems, which have been considered in this comparison. The absence of observed systems at mass ratios q < 0.5 might be because they have simply not been observed yet. As argued in Abdul-Masih et al. (2022), a larger sample size, which requires a dedicated search for these massive contact binary systems, should shed light on whether this is the case.

A second reason for the discrepancy could be that contact systems rapidly evolve to equal-mass systems, that is, to q = 1, and are therefore not observed as unequal-mass systems. However, following the same argument as in Abdul-Masih et al. (2022), the observed systems are uniformly distributed for q ≥ 0.5. Should contact systems all equalise in mass, the bulk of them would be expected to be observed at q ≈ 1. Furthermore, they find that based on the observed stability of the orbit, the systems in their sample will continue to evolve as unequal-mass contact binaries on a nuclear timescale.

Based on the different contributions to the PDF in Fig. 12, a third reason for the lack of contact systems at q < 0.5 can be found. The largest contribution to the PDF at these mass ratios comes from contact binaries that experience L₂-overflow. As described in Sect. 3.3, this can lead to considerable orbital angular momentum loss and thereby likely stellar mergers. Case-A contact binaries formed through runaway mass transfer, which we find at q < 0.45, also lead to mergers (see Sect. 3.5). In other words, contact systems are almost exclusively observed at q ≥ 0.5 because those at q < 0.5 merge shortly after they get into contact.

Considering that binaries with L₂-overflow or runaway mass transfer merge quickly, we are left with a PDF that at q < 0.5 is dominated by contact systems formed by the thermal-timescale expansion of the accretor. As explained in Sect. 3.1, such contact binaries either detach or merge quickly. In the former case, these systems would be observed as semi-detached binaries. In the latter case, we argue that if the thermal expansion of the accretor continues during contact, it eventually leads to L₂-overflow and a merger. Ge et al. (2023) find that the critical mass ratio q_crit below which binaries experience runaway mass transfer increases with decreasing metallicities. Therefore, at lower metallicities, we expect an even more significant fraction of binaries to experience runaway mass transfer and a subsequent quick merger. In conclusion, we find that contact systems are almost exclusively observed at q ≥ 0.5 because those at q < 0.5 merge or detach shortly after they get into contact.

5. Discussion

In this section, we put into perspective the influence of mass-transfer efficiency on our contact tracing results (Sect. 5.1), the occurrence of stable mass transfer from (super-)giant donors (Sect. 5.2), and the potential fates of binaries with non-ejected matter (Sect. 5.3).

5.1. The role of mass-transfer efficiency

The mean mass-transfer efficiencies $\overline{\beta }$ of binaries in our grid (Fig. F.1) are relatively low for Case-B and -C mass transfer. In these systems, the secondary usually reaches critical rotation after accreting M_accreted/M_2,  i < 3% of its initial mass (with M_accreted the accreted mass), which is lower than the ∼5 − 10% from Packet (1981) and consistent with the ∼2% found by Ghodla et al. (2023). It should be noted that only a few percent of the secondary’s outer envelope in terms of its mass is rotating near the critical rotation rate. With more efficient angular momentum transport, the relative mass accretion fraction would be higher than 3% because of a delay of critical surface rotation.

Because the secondary stars in our models accrete less matter than those in more conservative models, they also expand less. Compared to models with higher mass-transfer efficiencies, such as those from Pols (1994), Wellstein et al. (2001), de Mink et al. (2007), Claeys et al. (2011), and Menon et al. (2021), we find fewer contact systems from accretor expansion. The difference is particularly noticeable for Case-B binaries, where virtually none of such systems are found in our models. Our results generally agree better with the model grids from Fragos et al. (2023). This is expected since similar assumptions with regard to mass-transfer efficiency are made in their models.

From observations of Be X-ray binaries in the SMC, Vinciguerra et al. (2020) infer a mass transfer efficiency of ∼30% for the progenitor systems (i.e. the systems before the primary has formed a compact object). Potential progenitors of Be X-ray binaries in our grid, those with post-CHeB (core-helium burning) and MS components (see Figs. C.1–C.3), have a wide range of mass-transfer efficiencies (see Fig. F.1). For Case-A binaries with M_1,  i = 6.1 − 8.6 M_⊙ and Case-C binaries with M_1,  i = 5.2 − 14.2 M_⊙ the mass-mass transfer efficiency is consistent with the observationally constrained value of ∼30%.

Using observations of the classical Algol system δ Librae, Dervişoǧlu et al. (2018) found a mass-transfer efficiency of ∼100%. Based on the component masses derived for δ Librae, the initial primary mass would have been between ≈2.6 M_⊙ (q_i = 1.0) and ≈4.1 M_⊙ (q_i = 0.25, stable mass transfer). Our Case-A mass-transfer efficiency for M_1,  i ≈ 2.6 M_⊙ is ∼95%, which is consistent with the observationally constrained value. For M_1,  i ≈ 4.1 M_⊙, we find a Case-A mass-transfer efficiency of ∼ 25%, which is inconsistent with the value inferred for δ Librae.

Sen et al. (2022) computed a grid of Case-A systems with M_1,  i = 10 − 40 M_⊙ at LMC metallicity (Z = 0.0047) and similar assumptions as ours with regards to rotation, tides and accretion of angular momentum. For thermal-timescale mass-transfer phases in systems avoiding contact on the MS, they find relatively low (time-averaged) mass-transfer efficiencies (peak around ∼0.06). During the slower, nuclear-timescale mass transfer phases, the mass-transfer efficiency peaks around ∼0.94. This is consistent with our models, in which the rapid spin-up of the accretor during thermal-timescale mass transfer results in non-conservative mass transfer. For M_1,  i > 10 M_⊙, the time-averaged mass-transfer efficiencies in our models range from ∼0.4 to ∼0.65 (Fig. F.1). The thermal timescale is several orders of magnitude shorter than the nuclear timescale (Eqs. (6) and (7)). From estimating the time-averaged mass-transfer efficiencies over the whole Case-A phase (thermal- and nuclear-timescale mass-transfer phases) in the models of Sen et al. (2022), we see that their values of the mass-transfer efficiency are closer to unity. The slightly different definitions of the mass-transfer efficiency and the different metallicities (and, therefore, wind mass-loss rates) between our grids likely explain this difference.

Our models’ current assumption that stars cannot accrete anymore once they are at breakup velocities might be up for debate. Interactions with an accretion disk have been proposed to spin down the surface of the accretor and allow for additional accretion of material (Paczynski 1991; Popham & Narayan 1991). This contributes to the fact that the incidences of contact systems reported in Sect. 4.4 are lower limits because more accretion allows the accretor to expand more and potentially fill its Roche lobe. This could in turn increase the stellar merger incidences.

5.2. The stability of Case-B and -C mass transfer

In Sect. 4.2 we have described how non-conservative mass transfer and stellar wind mass loss from the donor prior to mass transfer can lead to stable mass transfer from (super-)giant donors (Case Bl and C) with initial mass ratios down to 0.1. Following the discussion in Sect. 5.1, these outcomes might also change with different assumptions regarding the mass-transfer efficiency. Just as for the formation of contact binaries, an increase in the mass-transfer efficiency leads to a higher incidence of classical CEs.

The presence of stable Case-Bl and -C mass transfer in our models agrees with what was found by Chen & Han (2008) for stars with M_1,  i ≲ 8 M_⊙, who also note that the effect of wind mass loss prior to the onset of mass transfer only has a minor effect on the stability.

We find that our critical mass ratios q_i,  crit for Case-Bl and -C mass transfer (Table 1) are higher than those from the adiabatic mass loss computations of Ge et al. (2010) for M_1,  i > 10 M_⊙, while being in relatively good agreement for stars with lower initial primary masses. These differences are likely caused by the fact that contrary to Ge et al. (2010), we also take the response of the accretor and orbit into account. Another difference is that Ge et al. (2010, 2015, 2020a) derive the critical mass ratios under the assumption of fully conservative mass transfer. Our critical mass ratio ranges agree relatively well with similar simulations with fully non-conservative mass transfer (Ge, priv. comm.).

Picco et al. (2024) use the adiabatic mass–radius exponents ξ_R,  ad (see Sect. 3.5) from Ge et al. (2020a) to determine the stability of mass transfer in their detailed binary evolution models. Our critical mass ratios agree relatively well with the ones from their simulations of fully non-conservative mass transfer in binaries with M_1,  i = 8.0 M_⊙.

The adiabatic mass–radius exponents ξ_R,  ad are also used by Li et al. (2023) in their binary population synthesis computations. For binaries with M_1,  i = 8.0 M_⊙ with non-conservative (β = 0 − 0.5) Case-Bl and -C mass transfer, they find values for q_i,  crit ≈ 0.37 − 0.66, which is in broad agreement with our values.

Similar stable Case-Bl and -C mass transfer has been found by Ercolino et al. (2023) for M_1,  i = 12.6 M_⊙. The reported critical mass ratio q_crit = 0.525 − 0.625 is slightly higher than what is found in our models for stars with similar masses (Table 1). This difference can be attributed to the different stability criteria for mass transfer used in their work.

5.3. The fate of binaries with non-ejected matter

In the majority of our models with β < 1, the non-accreted mass cannot be ejected to infinity. This raises the question of where this excess material resides. Should the excess matter remain in the accretor’s Roche lobe or the binary’s L₂-lobe, it can fill it up and lead to a situation similar to a contact binary. In this case, our grid would contain significantly more contact systems. For the orbital evolution, this scenario would essentially correspond to that of systems with higher mass-transfer efficiencies. Should the excess matter overfill the L₂-lobe and hydrodynamic drag becomes significant, a classical CE could form. Alternatively, the matter could settle in a circumbinary disk or shell, which can significantly influence the further evolution of the binary (Wei et al 2023).

Alternatively, a circumstellar disk may form, through which the secondary can potentially continue to accrete mass (Paczynski 1991; Popham & Narayan 1991). In this scenario, systems might show signs of circumstellar disks such as Hα excess and emission features in Be/Oe stars. This could explain the fast-rotating Hα emitters in the extended MS turnoff observed in young clusters (Milone et al. 2018). Moreover, Bodensteiner et al. (2021) have found a significant close-binary fraction of ${34}_{-7}^{+8}\%$ for one of the clusters analysed in Milone et al. (2018), NGC 330.

5.4. The onset of contact through runaway mass transfer and L₂-overflow

In Sect. 4.1, we find that runaway mass transfer and/or L₂-overflow are responsible for the onset of contact for a significant fraction of the Case-A, -Be, -Bl, and -C binaries. In Case-A binaries, L₂-overflow is of lesser importance for the onset of contact but leads to stellar mergers in contact binaries formed through the expansion of the accretor. Runaway mass transfer in Case-A systems leads to the formation of unstable contact binaries and subsequent stellar mergers. The same happens in Case-Be binaries, but now in unison with L₂-overflow. The added effect of the orbital angular momentum loss associated with L₂-overflow contributes to the instability of the contact binary formed through runaway mass transfer. In Case-Bl and -C binaries, the onset of runaway mass transfer and L₂-overflow often occur quasi simultaneously. In these systems, the L₂-overflow serves as an additional indication that the secondary is being engulfed by the rapidly expanding envelope of the (super-)giant primary star during runaway mass transfer. Lastly, some Case-Bl and C systems do not experience runaway mass transfer but do have primary stars that extend far beyond the L₂-lobe. In both cases, we expect the onset of a classical CE phase.

6. Summary and conclusions

Using a grid of ∼6000 detailed binary evolution models including rotation, tidal interactions, the evolution of both components, and with component masses between 0.5 and 20.0 M_⊙, we examine in which regions of the initial binary parameter space we expect contact phases, such as contact binaries and classical common-envelope (CE) phases, to occur. We identify five mechanisms that lead to contact: the expansion of the accretor, runaway mass transfer, L₂-overflow, orbital decay because of tides, and non-conservative mass transfer.

We find that accretor expansion, L₂-overflow, and runaway mass transfer lead to the formation of contact binaries in Case-A and -Be systems, and L₂-overflow and runaway mass transfer to the onset of classical CEs in Case-Bl and -C systems. This distinction stems from the fact that primary stars in Case-Bl and Case-C systems have extended envelopes with a clear core-envelope boundary, which engulfs the more compact MS companion upon contact. Case-A binaries with initial primary masses below 2 M_⊙ also form contact binaries because of the orbital decay caused by tides. Overall, the incidences of mass-transferring binaries forming contact binaries or entering a classical CE phase are ≥41% and ≥34% for M_1,  i ∈ [4.8;12.6) M_⊙ and M_1,  i ∈ [12.6;20.8] M_⊙, respectively. Over the entire mass range of M_1,  i ∈ [4.8;20.8] M_⊙, the incidence is ≥40%. These numbers are lower limits because they do not take into account the binaries that might enter a contact phase because of the interaction with non-accreted matter. Such systems are fairly common in our grid, and could alternatively result in systems with circumstellar disks and Hα emission features. Potential observational counterparts could be found in the extended MS turnoff in young stellar clusters.

We find that mass transfer is non-conservative in a large part of the initial binary parameter space, which is caused by the spin-up of the accretors to critical rotation after accreting < 3% of their own mass. The mass transfer efficiencies are 15 − 65%, 5 − 25%, and 25 − 50% in Case-A,-B, and -C mass transfer, respectively, for primary-star masses above 3 M_⊙. The incidence of systems entering a contact phase is sensitive to the mass-transfer efficiency. Given the relatively low mass-transfer efficiencies in our models, we might be underestimating this incidence. Another consequence of non-conservative mass transfer is that Case-Bl and -C mass transfer is stable for mass ratios ≥0.15 − 0.35.

By assuming that systems Case-A and -Be systems with M_1,  i  ∈  [4.8;20.8] M_⊙ that experience L₂-overflow and/or runaway mass transfer merge, and Case-Bl and -C systems experiencing this enter a classical CE, we find stellar merger and classical CE incidences among mass-transferring binaries of ≥12% and ≥19%, respectively. If we relax this assumption and assume that all contact binaries, except those for which either component reaches core-C exhaustion, eventually merge, the stellar merger incidence increases to ≥19%. Just as for the incidence of mass-transferring binaries reaching a contact phase, the stellar merger and classical CE incidences are lower limits, which are sensitive to the mass-transfer efficiency.

Lastly, we compare our population of massive (M_1,  i ≳ 5 M_⊙) Case-A contact systems with observed (near-)contact systems. We find from our models that approximately two to three times as many contact binaries form with mass ratios q < 0.5 upon contact than with q > 0.5. Moreover, the majority of the systems with q < 0.5 form contact binaries because of runaway mass transfer or the thermal-timescale expansion of the accretor, with subsequent L₂-overflow in more than half of the cases. Therefore, most of these systems quickly merge or detach. This is in agreement with the observations since almost no (near-)contact systems are observed with mass ratios below 0.5. Out of all systems forming contact binaries, 17% do so on the accretor’s nuclear timescale and avoid L₂-overflow, and are expected to be stable and long-lived. The majority of these systems have q > 0.5 upon contact, which is again in excellent agreement with the observations.

Mergers and classical CEs are some of the most complex and fascinating outcomes of binary evolution. While we have identified the physical mechanisms that can lead to contact, the outcome of mergers and classical CEs still need to be explored in more detail and leave many open questions.

Acknowledgments

We wish to thank the referee, Hongwei Ge, for his valuable comments and suggestions, which have helped us further improve this work. We thank S. Hekker, Ph. Podsiadlowski, D. Wei, V. Bronner, P. Marchant, and C. Wang for the meaningful discussions, suggestions, and comments. The authors acknowledge support from the Klaus Tschira Foundation. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 945806). This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1-390900948 (the Heidelberg STRUCTURES Excellence Cluster).

References

Appendix A: Nuclear timescale expansion of the accretor

Here, we present an example of a model forming a contact binary from the nuclear timescale expansion of the secondary (accretor) star. The evolution of the Case-A system with M_1,  i = 10.2 M_⊙, q_i = 0.9, and a_i = 12.5 R_⊙ is shown in Fig. A.1. After regaining thermal equilibrium, the primary star detaches and later fills its Roche lobe again because of nuclear expansion. The detachment phase is short compared to the evolutionary timescale and only visible by a slight wiggle in the HRD (Fig. A.1a). During the first, thermal-timescale mass-transfer phase, the secondary star expands on its thermal timescale (Fig. A.1b). However, during the second, nuclear-timescale mass-transfer phase, the secondary remains in thermal equilibrium and its radius evolves on a nuclear timescale (Fig. A.1c) and increases with mass accretion. Because R₂ grows faster in time than R_RL,  2, the secondary eventually fills its Roche lobe simultaneously with the primary, and a contact binary is formed.

Fig. A.1. Same as Fig. 3 but for a M1,  i = 10.2 M⊙, qi = 0.9, and ai = 12.5 R⊙ binary.

Appendix B: Contact tracing results for other M_1,  i

Figs. B.1–B.3 contain the contact tracing results for the initial primary masses M_1,  i = 0.8, 0.9, 1.1, 1.3, 1.9, 2.2, 2.6, 3.1, 3.7, 4.3, 5.2, 6.1, 7.2, 8.6, 12.0, 13.1, 14.2, 15.6, 16.9, and 18.4 M_⊙.

Fig. B.1. Contact tracing results for M1,  i = 7.2 − 18.4 M⊙. Models with M1,  i = 7.2 M⊙ and initial separations larger than the Case-C systems experience numerical issues after the TP-AGB phase (see Sect. 2.1.2). Only the model with qi = 0.1 avoids these issues and does not initiate mass transfer.

Fig. B.2. Contact tracing results for M1,  i = 1.9 − 6.1 M⊙. Models with M1,  i = 3.7 − 6.1 M⊙ experience numerical issues after the TP-AGB phase (see Sect. 2.1.2). This explains the unexpected onset of mass transfer at initial separations larger than those of systems avoiding mass transfer. For M1,  i = 1.9 − 2.2 M⊙ we find certain models where mass transfer starts when the primary is on the WD cooling track and experiences sudden radial expansion.

Fig. B.3. Contact tracing results for M1,  i = 0.8 − 1.3 M⊙. For M1,  i = 1.1 − 1.3 M⊙ we find certain models where mass transfer starts when the primary is on the WD cooling track and experiences sudden radial expansion.

Appendix C: Evolutionary states at contact or termination of models

In Figs. C.1–C.3, we show the evolutionary state of the binary components for all systems in our grid at contact or at termination for systems that avoid contact.

Fig. C.1. Evolutionary state of the primary (left triangle) and secondary (right triangle) at contact or termination for M1,  i = 8.6 − 20.0 M⊙. The pink, blue and grey background colours indicate systems that get into contact, avoid contact and fail to converge numerically, respectively. Post-MS stars are those that have exhausted hydrogen but have not yet ignited helium in their core. ‘CHeB’ stands for core helium burning.

Fig. C.2. Same as for Fig. C.1 but for M1,  i = 1.9 − 7.2 M⊙.

Fig. C.3. Same as for Fig. C.1 but for M1,  i = 0.8 − 1.6 M⊙.

Appendix D: Expansion timescales for Case-A contact systems

In Figs. D.1–D.2, the Case-A region of the initial binary parameter space is shown for all initial primary masses of the grid. The contact systems formed through the expansion of the accretor and those which experience L₂-overflow are indicated with the same colour scheme as in Fig. 7. Each model is marked with a red or blue square based on whether the mean of the accretor’s expansion timescale is of the order of its mean thermal or nuclear timescale, respectively.

Fig. D.1. Expansion timescales for Case-A contact systems formed through the expansion of the accretor star with M1,  i = 5.2 − 20.0 M⊙. In models marked with filled red (blue) squares, the accretor expands on its thermal (nuclear) timescale prior to the onset of contact. The colour scheme is the same as in Fig. 7. The dark blue solid line indicates the division between Case-A and -Be systems.

Fig. D.2. Same as Fig. D.1 but for M1,  i = 0.8 − 4.3 M⊙.

Appendix E: Assignment criteria for evolutionary outcomes of models with numerical issues

Table E.1 lists the criteria used to determine the evolutionary outcome of models that experience numerical issues. The outcome assignment criteria for these models are based on the outcomes of neighbouring models and information from equivalent systems at different initial primary masses.

Appendix F: Mass-transfer efficiency

We show the mean mass-transfer efficiency $\overline{\beta }$ for Case-A, -B, and -C mass transfer as a function of the initial primary mass M_1,  i in Fig. F.1. We use the birth probabilities p_birth (Sect. 2.5) of each model as weights for the computation of these mean values. The value of the mass-transfer efficiency for each individual model is averaged over time. We only consider mass transfer in systems that avoid contact.

Fig. F.1. Mean mass-transfer efficiency $\overline{\beta }$ per mass-transfer case for each initial primary mass M1,  i.

First, we see that for all mass-transfer cases, the mean mass-transfer efficiency $\overline{\beta }$ is relatively high for log₁₀(M_1,  i) < 0.25 − 0.45. For log₁₀(M_1,  i) > 0.25 − 0.45, $\overline{\beta }$ has moderate to low values. The high values of $\overline{\beta }$ for binaries with lower initial primary masses should be regarded as an approximation since for all mass-transfer cases the models avoiding contact are sparse due to numerical issues (Fig. B.2–B.3). The models that avoided numerical issues have high mass-transfer efficiencies.

Case-A mass transfer has low to moderate mean mass-transfer efficiencies for log₁₀(M_1,  i) > 0.45. The initially closest Case-A binaries go through conservative mass transfer before they form contact binaries and are, therefore, not taken into account for the computation of $\overline{\beta }$. Case-A binaries that avoid contact are on initially wider orbits. The accretors in these systems have longer tidal synchronisation times, so tides are unable to prevent them from rotating critically. As discussed in Sect. 5.1, Case-A mass transfer consists of a short thermal-timescale phase, which is mostly non-conservative, and a more conservative, longer nuclear-timescale phase. This results in values of $\overline{\beta }$ between ∼ 0.15 and ∼ 0.65.

Case-B mass transfer in binaries with log₁₀(M_1,  i) > 0.25 have relatively low values of $\overline{\beta }$, with values between ∼ 0.05 and ∼ 0.25. In virtually all systems with Case-B mass transfer, the accretors reach critical rotation after accreting a few percent of their own mass, which quenches accretion (Sect. 4.1 and 5.1).

Case-C mass transfer is generally highly non-conservative (Sect. 4.2). However, the contribution of Case-C systems with relatively short mass-transfer phases before core-C exhaustion in which the accretor is not spun up to critical rotation increases the value of $\overline{\beta }$ for log₁₀(M_1,  i) > 0.25.

The contribution of both conservative and non-conservative mass transfer in Case-A and -C mass transfer cases, for the former in terms of mass-transfer phases (thermal and nuclear) and for the latter in terms of systems (highly non-conservative models and conservative models), results in somewhat similar mean mass-transfer efficiencies.

Appendix G: Table with contact tracing results

Table G.1 contains an extract of the table containing the contact tracing results for all our models.

All Tables

All Figures

	Fig. 1. Radial expansion of a 10.2 M⊙ star (solid black line). The radial evolution is divided into different cases based on the main phases of expansion (solid horizontal lines). Each dashed horizontal line indicates a sample point at which we put R = RRL.
In the text

Fig. 2. Schematic representation of the physical mechanisms leading to contact (not to scale). The filled blue circles and purple squares indicate the position of the L1- and L2-points, respectively. The filled grey-blue circles represent the stellar cores in panels 3a and 3b. Panel 1a shows the expansion of the accretor leading to a contact binary (1b). This corresponds to the ‘Accretor expansion’ mechanism described in Sect. 3.1. Subsequent overfilling of the components’ Roche lobes leads to the formation of an overcontact binary (2b). The primary increasingly overfills its Roche lobe (2a) and can eventually fill the secondary’s Roche lobe (2b). This can eventually lead to L2-overflow (2c), which likely results in a stellar merger. The scenarios (1b–2b–2c) and (2a–2b–2c) correspond to the ‘L2-overflow’ mechanism described in Sect. 3.3. In panel 3a, runaway mass transfer from a (super-)giant (left) to an MS star (right) leads to the onset of a classical common-envelope phase 3b, where the cores of both stars revolve in the (super-)giant’s envelope. This corresponds to the ‘Runaway MT’ mechanism described in Sect. 3.5.

In the text

Fig. 3. Example of a M1,  i = 10.2 M⊙, qi = 0.4, and ai = 12.4 R⊙ binary system forming a contact binary during the MS of the primary (Case A) because of the thermal expansion of the secondary (accretor) star. Panel a shows an HRD with the evolutionary tracks of the primary and secondary stars (solid lines). The dashed black lines show the evolutionary tracks of single stars with the same initial masses as the binary components and initial rotation rates of Ω/Ωc = 0.25. Panel b shows the evolution of the radius R (solid lines) and Roche lobe radius RRL (dashed lines) of both components. The timescales governing the evolution of the binary and the mass-transfer rate (black dashed line) are shown in panel c as a function of the decreasing primary star mass. The secondary’s nuclear (τnuc,  2) and thermal (τKH,  2) timescales are shown with a dashed green and pink line, respectively. The expansion (τR/Ṙ,2) and mass-transfer (τtrans) timescales are shown with a solid blue and orange line, respectively.

In the text

Fig. 4. Example of a M1,  i = 10.2 M⊙, qi = 0.9, and ai = 39.2 R⊙ binary star model in which non-conservative mass transfer occurs and for which the non-accreted matter cannot be driven to infinity. Panel a is the same as Fig. 3a (‘ign.’ = ‘ignition’ and ‘exh.’ = ‘exhaustion’). In panel b, the mass-loss rate of the secondary (solid black line) is compared to the maximum mass-loss rate Ṁmax (dashed lines; pink for feff = 1.0, green for feff = 0.1) set by Eq. (8). The surface rotation rates (solid blue and orange lines) and mass transfer efficiency (dashed black line) are shown in panel c as a function of the primary mass.

In the text

Fig. 5. Example of a binary model with M1,  i = 10.2 M⊙, qi = 0.2, and ai = 175.7 R⊙ going through a phase of (delayed) runaway mass transfer. Panel a is the same as Fig. 3a (‘ign.’ = ‘ignition’). Panel b shows the mass-transfer rate Ṁtrans (solid black line), thermal mass-transfer rate ṀKH (solid green line) and dynamical mass-transfer rate Ṁdyn (solid gold line) for the primary star on the left axis. The right axis shows the component mass evolution as a function of age (dashed blue and orange lines). In panel c, the radius R (solid blue line), Roche lobe radius RRL (dashed light-blue line), L2-volume-equivalent radius RL2 (dahsed green line, Misra et al. 2020, Eq. (15)) and orbital separation aorb (dashed black line) evolution for the primary are shown as a function of the decreasing primary mass. The blue cross indicates the moment of L2-overflow.

In the text

Fig. 6. Example of a M1,  i = 1.1 M⊙, qi = 0.8 and ai = 4.1 R⊙ binary system in which tides lead to contact. In Panel a we show the exchange between the orbital angular momentum Lorb (solid black line) and the secondary’s spin angular momentum S2 (solid orange line), which coincides with the decrease in Lorb/(S1 + S2) (dashed green line). The primary’s spin angular momentum S1 is shown with a solid blue line. Panel b shows the evolution of the moment of inertia (solid lines) and tidal synchronisation timescale (dashed lines) for both components. While the moment of inertia of the primary decreases around M1 = 0.3 M⊙, it sharply increases for the secondary. At the same time, the secondary star’s tidal synchronisation timescale decreases by approximately two orders of magnitude. Panel c shows the evolution of the primary’s and secondary’s radii (solid lines) and Roche lobe radii (dashed lines). While the former fills its Roche lobe, tides cause orbital shrinkage, which results in the secondary also filling its Roche lobe.

In the text

Fig. 7. Occurrence of contact phases for models with initial primary masses M1,  i = 10.2 M⊙ on the initial mass ratio–separation plane. Models marked with a dot are those for which accretion was limited to 0.1 times the secondary’s global thermal timescale τKH (see Sect. 2.2.1). The other ones are marked with a star symbol. The dark-blue quasi-horizontal lines indicate the initial mass transfer cases, which can be read from the right side. Systems on the left of the solid black line are Darwin unstable at the onset of mass transfer according to Eq. (9) and assuming R1 = RRL,  1.

In the text

	Fig. 8. Same as Fig. 7, but for M1,  i = 20.0 M⊙ (left) and M1,  i = 1.6 M⊙ (right).
In the text

	Fig. 9. Example of a M1,  i = 10.16 M⊙, qi = 0.6, and ai = 1398.2 R⊙ binary model with stable Case-C mass transfer. Panel a is the same as Fig. 3a (‘ign.’ = ‘ignition’). Panel b is the same as Fig. 5b. Panel c is the same as Fig. 5c, with the addition of a solid orange line indicating the evolution of the mass ratio.
In the text

Fig. 10. Sunburst charts displaying the fractions of evolutionary outcomes for mass-transferring binary systems in the grid over initial primary mass ranges of [4.8; 12.6) M⊙ (left) and [12.6; 20.8] M⊙ (right). Wedges with a percentage < 1% are not labelled. Top row: the inner level shows the principal outcome of the evolution (‘Accretor expansion’, ‘Runaway MT’, ‘L2-overflow’, and ‘No contact’), while the outer level shows the ancillary outcome (‘L2-overflow’, ‘Non-conservative MT + cannot eject’). Ancillary ‘No contact’ outcomes are incorporated in the inner ‘No contact’ category. Models in the ‘Numerical issues’ category are assigned a likely evolutionary outcome based on their initial mass ratio qi and mass-transfer case, displayed on the sunburst chart’s intermediate level. The labels ‘A’, ‘Be’, and ‘Bl/C’ refer to Case-A, Case-Be, and Case-Bl or -C mass transfer, respectively. Bottom row: the inner level shows the percentage of Case-A, -Be, -Bl, and -C systems. We show the principal outcome per case on the middle level. The outer level shows the lower limits of the total fraction of stellar mergers and classical CEs per mass-transfer case.

In the text

	Fig. 11. Sunburst charts displaying the fractions of evolutionary outcomes for mass-transferring binary systems in the grid over initial primary mass ranges of [4.8; 20.8] M⊙. The left and right charts are equivalent to those in the top and bottom row of Fig. 10, respectively.
In the text

Fig. 12. Probability density function (PDF) of contact systems formed in Case-A binaries with M1,  i = 7.9 − 20.8 M⊙ as a function of the mass ratio q at the onset of contact. Data points (filled squares) show the observed mass ratio and primary mass (right axis) for all observed MW, LMC and SMC contact and near-contact systems compiled in Menon et al. (2021), including the uncertainties on their values. Systems without reported uncertainties on q are indicated with a dash symbol, and those without reported uncertainties on q and M1 are indicated with a star symbol.

In the text

	Fig. A.1. Same as Fig. 3 but for a M1,  i = 10.2 M⊙, qi = 0.9, and ai = 12.5 R⊙ binary.
In the text

	Fig. B.1. Contact tracing results for M1,  i = 7.2 − 18.4 M⊙. Models with M1,  i = 7.2 M⊙ and initial separations larger than the Case-C systems experience numerical issues after the TP-AGB phase (see Sect. 2.1.2). Only the model with qi = 0.1 avoids these issues and does not initiate mass transfer.
In the text

Fig. B.2. Contact tracing results for M1,  i = 1.9 − 6.1 M⊙. Models with M1,  i = 3.7 − 6.1 M⊙ experience numerical issues after the TP-AGB phase (see Sect. 2.1.2). This explains the unexpected onset of mass transfer at initial separations larger than those of systems avoiding mass transfer. For M1,  i = 1.9 − 2.2 M⊙ we find certain models where mass transfer starts when the primary is on the WD cooling track and experiences sudden radial expansion.

In the text

	Fig. B.3. Contact tracing results for M1,  i = 0.8 − 1.3 M⊙. For M1,  i = 1.1 − 1.3 M⊙ we find certain models where mass transfer starts when the primary is on the WD cooling track and experiences sudden radial expansion.
In the text

Fig. C.1. Evolutionary state of the primary (left triangle) and secondary (right triangle) at contact or termination for M1,  i = 8.6 − 20.0 M⊙. The pink, blue and grey background colours indicate systems that get into contact, avoid contact and fail to converge numerically, respectively. Post-MS stars are those that have exhausted hydrogen but have not yet ignited helium in their core. ‘CHeB’ stands for core helium burning.

In the text

	Fig. C.2. Same as for Fig. C.1 but for M1,  i = 1.9 − 7.2 M⊙.
In the text

	Fig. C.3. Same as for Fig. C.1 but for M1,  i = 0.8 − 1.6 M⊙.
In the text

	Fig. D.1. Expansion timescales for Case-A contact systems formed through the expansion of the accretor star with M1,  i = 5.2 − 20.0 M⊙. In models marked with filled red (blue) squares, the accretor expands on its thermal (nuclear) timescale prior to the onset of contact. The colour scheme is the same as in Fig. 7. The dark blue solid line indicates the division between Case-A and -Be systems.
In the text

	Fig. D.2. Same as Fig. D.1 but for M1,  i = 0.8 − 4.3 M⊙.
In the text

	Fig. F.1. Mean mass-transfer efficiency $\overline{\beta }$ per mass-transfer case for each initial primary mass M1,  i.
In the text