Contact tracing of binary stars: Pathways to stellar mergers

Stellar mergers lead to diverse phenomena: rejuvenated blue stragglers, magnetised and peculiar stars, transients and nebulae. Using a grid of about 6000 detailed 1D binary evolution models (initial component masses of 0.5-20$\,\text{M}_{\odot}$ at solar metallicity), we investigate which initial binary-star configurations lead to contact and classical common-envelope (CE) phases and assess the likelihood of a subsequent merger. Considering rotation and tides, we identify five mechanisms leading to contact and mergers: runaway mass transfer, $\text{L}_{2}$-overflow, accretor expansion, tidally-driven orbital decay, and non-conservative mass transfer. At least 40% of mass-transferring binaries with initial primary masses of 5-20$\,\text{M}_{\odot}$ enter contact, with>12% and>19% likely merging and evolving into a classical CE phase, respectively. Classical CE evolution occurs in late Case-B and Case-C binaries for initial mass ratios $q_{\text{i}}$<0.15-0.35, stable mass transfer for larger $q_{\text{i}}$. Early Case-B binaries enter contact for $q_{\text{i}}$<0.15-0.35 and in initially wider Case-A binaries, this occurs for $q_{\text{i}}$<0.35. All initially closest Case-A systems form contact binaries. We predict that binaries entering contact with $q$<0.5 merge or detach on a thermal timescale, while those formed with $q$>0.5 lead to long-lived contact phases. The fact that contact binaries are almost exclusively observed with $q$>0.5 confirms our expectations. Our contact, merger and classical CE incidences are lower limits because the mass transfer in our models is non-conservative. In most binaries, the non-accreted mass cannot be ejected and may settle in disks or lead to contact phases and mergers. Overall, contact binaries are a frequent and fascinating result of binary mass transfer of which the exact outcomes still remain to be understood and explored further.


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. 2020Toonen et al. , 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).
3D simulations of stellar mergers with magnetohydrodynamics (MHD, Schneider et al. 2019) and smoothed particle hydro-dynamics (SPH, Sills et al. 1997Sills et al. , 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 lobe1 , 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 find out 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 done, e.g. 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 ratios2 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 (Ste ¸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.

Methods
We compute a grid of 5957 1D binary evolution models using the binary module of MESA (release 12778; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019)).First, we describe the adopted singlestar 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.

Adopted stellar physics
Each binary component is initialised from a precomputed zeroage main-sequence (ZAMS) model at solar metallicity, that is, Z = 0.0142 and Y = 0.2703 (Asplund et al. 2009).We use 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 allow 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 Ω c = GM/R 3 eq ≃ 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 are hydrostatic, meaning that MESA's implicit hydrodynamic solver is disabled.We use the approx21 nuclear network.

Mixing
Convective mixing is 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).The Ledoux criterion is used to identify regions in the star that are unstable to convection.The efficiency of semi-convective mixing is set to α sc = 10.0 (Schootemeijer et al. 2019).Thermohaline mixing is also included with an efficiency of α th = 1.0 (Marchant et al. 2021).
Convective boundary mixing (CBM) is included via the stepovershoot scheme, in which we allow the convective hydrogenburning 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 use 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. 2023).Because of this, we follow Marchant et al. (2021) and use 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 are included (see for example Heger et al. 2000 for a detailed description).Additional diffusion of angular momentum is taken care of by the Spruit-Tayler dynamo.We scale 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 µ , is set to f µ = 0.1 (see also Pinsonneault et al. 1989).

Wind mass-loss prescription
Because our binary models contain low-, intermediate-and highmass stars, and the masses can change significantly over the star's evolution, we employ a wind mass-loss prescription that covers this large mass range and all evolutionary stages.
We consider two distinct regimes, the hot-wind regime with surface temperature3 T surf ⩾ 11 kK and the cool-wind regime for stars with T surf ⩽ 10 kK.We use 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) is computed via the Vink et al. (2000) prescription.When X surf drops below 0.4, either the prescription for Wolf-Rayet (WR) stars from Sander & Vink (2020) or the prescription for lowmass helium stars from Vink (2017) is used, depending on whether the star's luminosity L is higher or lower than a certain luminosity L 0 respectively.As described in Sander & Vink (2020), L 0 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 0 , we compute the mass-loss rate with both prescriptions and take the maximum value as the adopted wind mass-loss rate (J.Vink, 2021, priv. comm.).All of the aforementioned prescriptions have a scaling factor of 1.0.For 0.4 ⩽ X surf ⩽ 0.5, the mass-loss rate is determined via linear interpolation between the mass-loss rates from both regimes.
When a model reaches the cool-wind regime, the distinction is made between stars expected to become giants or supergiants.The cut is made at log 10 (L /L ⊙ ) = 3.15, where L is defined as in Langer & Kudritzki (2014), 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 use 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 use a scaling factor of 0.1 for the former and 0.2 for the latter.Models above the cut in L use the Nieuwenhuijzen & de Jager (1990) prescription with a scaling factor of 1.0.
We increase 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 aims 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 aim 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 is not successful in avoiding numerical issues in each model, especially in those models where the aforementioned instabilities occur.Because of this, we opt 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.).

Adopted binary physics
In our models, we evolve both binary components.This allows 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 is considered.Hence, in this work, references to primary (secondary) and donor (accretor) are equivalent.We will mostly employ the names primary (subscript "1") and secondary (subscript "2") star.We assume that all binaries are on circular orbits.

Mass transfer and accretion
Whenever the primary star is on the MS, mass transfer is 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. 2022Fabry et al. , 2023)).For systems with post-MS primary masses smaller than 1.3 M ⊙ , the Kolb scheme (Kolb & Ritter 1990) is used, because this scheme is better suited for envelopes with larger pressure scale heights (Fragos et al. 2023).Both schemes are 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, β ≡ − Ṁ2 / Ṁtrans .In this definition, Ṁtrans < 0 and Ṁ2 > 0. Ṁ2 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 spins up to its critical rotation rate Ω c, 2 , mass transfer is non-conservative and β < 1.When the accretor star's accretion timescale τ acc ≡ M 2 / Ṁtrans approaches the star's dynamical timescale τ dyn, 2 , defined for a star with sound speed c S as τ dyn, 2 ≡ R 2 /c S (Kippenhahn et al. 2012), we limit the accretion rate to that of 0.1 times the star's Kelvin-Helmholtz (or thermal) timescale τ KH .This also results in non-conservative mass transfer, that is, β < 1.We do this because the accreting models tend not to converge numerically when τ acc ∼ τ dyn, 2 .Models for which the accretion rate is limited are marked in the results.The accretion of angular momentum during mass transfer is computed following Lubow & Shu (1975) and Ulrich & Burger (1976), which includes accretion through ballistic impact and a Keplerian disk.

Tides and angular momentum loss
Tidal synchronisation is 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 are equal to the orbital period.Orbital angular momentum in our models evolves via mass loss from the system (through winds or isotropic reemission) and spin-orbit coupling.In the former case, the lost mass is assumed to have the specific angular momentum of the star's orbit in which vicinity it is leaving the system.

Binary-star grid
In our grid, we vary 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 are selected between 0.8 and 20.0 M ⊙ with logarithmic mass spacing ∆ log 10 M 1, i ≈ 0.074.At the high-mass end, we add 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 are considered for the mass ratios, with linear spacing and steps of 0.1.We impose 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 add an additional mass ratio at q i = 0.97 to model twin systems.Lastly, we choose the initial binary separations a i such that the first phase of mass transfer occurs 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, respectively4 .This distinction is related to the phases in which the (primary) star expands, as illustrated in Fig. 1.Case B phases are 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 star5 .
To ensure sufficient sampling of all these stages, we have 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 2 /M 1 and separation a, (Eggleton 1983), For a given primary mass M 1 and mass ratio q, we select 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 assume the star fills its Roche lobe, R = R RL , and mass transfer ensues.From Eq. (2), we then get the separation a at which this happens.Assuming that the orbit stays approximately constant until mass transfer starts, we can take 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 (see, for example, the dividing line between Case Bl and Case C in Fig. 7).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).

Stopping conditions
In principle, we compute each binary model until one of the components, usually the initially more massive primary star, reaches the end of core carbon burning.This is defined as the moment when the central mass fractions of He and C fall below 10 −6 and 10 −2 , respectively.
Systems with MS primaries having predominantly radiative envelopes (M 1, i > 1.3 M ⊙ ) use MESA's contact mass-transfer scheme (see Sect. 2.2.1) and can therefore be modelled through such phases.During a contact phase, the computation terminates when the condition for mass overflow through the second Lagrange point (L 2 ) or L 2 -overflow from Marchant et al. (2016) is met (see Sect. 3.3).
For primaries with mostly convective envelopes (post-MS or M 1, i ≤ 1.3 M ⊙ ), which use the Kolb mass-transfer scheme.The computation stops whenever the secondary star overfills its Roche lobe by more than R RL, 2 .
When the mass-transfer rate Ṁtrans reaches a value of 10 M ⊙ yr −1 , the evolution is 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, is not considered in this study.The onset of reverse mass transfer is a stopping condition.

Birth probabilities
To describe the models in our grid as a population (Sect.4.3-4.5),we compute their birth probability p birth from the initial distribution function, which here 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 use the Kroupa (2001) initial mass function (IMF), with α = −2.3 as appropriate for primary star masses ≥ 0.5 M ⊙ and C m a normalisation constant.We assume the other two initial distribution functions to be uniform in q i and log 10 P i , respectively, ϕ(q i ) = C q and χ(log with C q and C p being normalisation constants.The birth probability p birth of each model is then computed by integrating the product of the initial distribution functions over the parameter size in the model grid, × ϕ(q i ) χ(log 10 P i ) dM 1, i dq i d log 10 P i . (5) Here, log 10 P u, l , q u, l and M u, l are the upper-and lowerboundaries of the parameter size of a model in the grid for log 10 P i , q i and M 1, i , respectively.They are chosen as the midpoints between the initial parameter values of the model and its neighbouring models.

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 (Fig. 2.1b and 2.2b).Systems with a (super-)giant primary (i.e.Case Bl & 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.
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.

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 Fig. 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 1 / Ṁtrans | is shorter than the secondary's thermal or Kelvin-Helmholtz timescale 7 defined as (Kippenhahn et al. 2012) the star will be 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 semidetached binary when the accretor regains thermal equilibrium and shrinks.When the secondary is in thermal equilibrium during mass transfer, it will expand on a nuclear timescale, defined as (Kippenhahn et al. 2012) 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 7 Technically one should compare with the local thermal timescale in the surface layers of the accreting star, defined in Kippenhahn et al. (2012) (see also Temmink et al. 2023) (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 L 2 -overflow (2c), which likely results in a stellar merger.The scenarios (1b-2b-2c) and (2a-2b-2c) correspond to the "L 2 -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.
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

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 to ∼10 percent 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 2 /2R 2 .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 1 + L 2 ) / Ṁ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 1 + L 2 ) of the binary (Marchant 2017) 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 Ṁ2 > Ṁ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 Sec.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 masstransfer 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 1 ≈ 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).

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 1 (through which mass transfer occurs), L 2 and L 3 .The L 2 -point is always located on the side of the less massive star in the binary, and L 3 on the side of the more massive star.In Case-A and -Be binaries, mass loss from the L 2 -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 2 (Marchant et al. 2021).In Case Bl and -C binaries, L 2 indicates the onset of a classical CE.
L 2 -overflow is traced 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 2 volume-equivalent radii from Marchant et al. (2016).This is done during the computation of the model and the moment of L 2 -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 L 2 and the distance to the L 2 -point D L 2 from Misra et al. (2020) (similar fitting formulae for R L 2 are derived in Ge et al. 2020).This is done in post-processing, meaning that the outflow from L 2 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 2 -overflow occurs, the models are labelled "L 2 -overflow".
The evolution of a binary system experiencing L 2 -overflow is shown in Fig. 5.The primary and secondary become underand over-luminous, respectively, during Case-Be mass transfer (Fig. 5a).As in the previous example, mass transfer is non- 0.9 and a i = 39.2R ⊙ 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 f eff = 1.0, green for f eff = 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.conservative and early on during the mass transfer the nonaccreted 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 L 2 , 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 L 2 , 1 .As a consequence, the star increasingly overfills its Roche lobe until it reaches the point of L 2 -overflow.The subsequent loss of mass and angular momentum from the L 2 -point will result in an accelerated orbital shrinkage, which is not accounted for in this model.The onset of L 2 -overflow is shown schematically in Fig. 2. In the evolutionary scenario depicted in Fig. 2.1a-1b-2b-2c, 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 2 -lobe (Fig. 2.2c).In the second scenario, from Fig. 2.2a-2c, 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 2 -lobe and L 2 -overflow occurs (Fig. 2.2c).

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 where L orb is the orbital angular momentum, Ω is the orbital angular rotation velocity, and I 1 and I 2 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 (ω 1 = ω 2 = Ω) system.Here, ω 1 and ω 2 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, 2 I 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 Fig. 5. Example of a binary model with M 1, i = 10.2 M ⊙ , q i = 0.2 and a i = 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 R RL (dashed light-blue line), L 2 -volume-equivalent radius R L 2 (dahsed green line, Misra et al. 2020, Eq. 15) and orbital separation a orb (dashed black line) evolution for the primary are shown as a function of the decreasing primary mass.The blue cross indicates the moment of L 2 -overflow.
increase in radius and density redistribution (the envelope now has a deep convective zone) both increase the moment of inertia I 2 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 2 increases and the orbital angular momentum L orb decreases (Fig. 6a).Eventually, the system becomes Darwin unstable, as it fulfils the condition8 L orb /(S 1 + S 2 ) ≲ 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.

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 ξ (Webbink 1984).
If ξ R 1 < ξ RL , mass transfer is unstable and the mass-transfer rate keeps increasing, potentially reaching values larger than 1 M ⊙ yr −1 .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 will form 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 Fig. 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 1 /τ KH, 1 .Secondly, the second time derivative of log 10 Ṁ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: ξ R 1 < ξ RL .For semi-detached models using MESA's contact mass-transfer scheme, ξ R 1 = ξ RL by definition.Hence, in this case, only the first two conditions are evaluated.Note that L 2 -overflow is not a necessary condition for the onset of runaway mass transfer.However, it does often occur in systems with runaway mass transfer.
Fig. 6.Example of a M 1, i = 1.1 M ⊙ , q i = 0.8 and a i = 4.1 R ⊙ binary system in which tides lead to contact.In Panel (a) we show the exchange between the orbital angular momentum L orb (solid black line) and the secondary's spin angular momentum S 2 (solid orange line), which coincides with the decrease in L orb /(S 1 + S 2 ) (dashed green line).The primary's spin angular momentum S 1 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 M 1 = 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.
An example of a binary system experiencing delayed runaway mass transfer is shown in Fig. 5 and was previously discussed to demonstrate L 2 -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 masstransfer rate.Simultaneously, we see in Fig. 5c that the slope of R 1 as a function of the decreasing M 1 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 2overflow 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. 2015Ge et al. , 2020)).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.

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 has been 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).

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 ratioseparation (q i -log 10 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 Sect.4.1.2and 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 online9 .An extract of this table can be found in Appendix G.
The coloured regions in Fig. 7 indicate the different physical mechanisms leading to contact ("Accretor expansion", "Nonconservative MT + cannot eject", "L 2 -overflow", "Tidallydriven 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 0.0 0.2 0.4 0.6 0.8 1.0 background colour (principal outcome), priority is given in the following order (highest to lowest priority): "Tidally driven contact", "Runaway MT", "Accretor expansion", "L 2 -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.45with log 10 (a i /R ⊙ ) ≲ 1.23 and q i = 0.75-1.00,these contact systems additionally experience L 2 -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 2 -lobe.This behaviour is found consistently throughout the range of M 1, i = 2.6-20.0M ⊙ 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 2 -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 10 (a i /R ⊙ ) ≳ 1.23) systems with q i = 0.15-0.35that avoid L 2 -overflow reach mass-transfer rates Ṁtrans orders of magnitude larger than their thermal masstransfer rate ṀKH .Ṁtrans eventually reaches the stopping condition value of 10 M ⊙ yr −1 (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 2 -overflow and merge.
The Case-A binary systems with q i ≤ 0.15 and log 10 (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 masstransfer 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.00and log 10 (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.00close 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 10 (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 10 (a i /R ⊙ ) ≲ 1.78 and q i = 0.15-0.35,all of them also experience L 2 -overflow.The secondary stars are still on the MS.The primaries are HG stars without a clear core-envelope boundary (see Fig . C.1-C.3 in Appendix C for the evolutionary state of all models at contact/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 10 (a i /R ⊙ ) ≲ 2.92, and q i ≤ 0.35 and log 10 (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 2 -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.

Initial primary mass of 20.0 M ⊙
For M 1, i = 20.0M ⊙ (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 10 (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.0M ⊙ and the latter for M 1, i = 5.2-10.2M ⊙ .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 2 -overflow (Fig. 8, left).However, at q i ≥ 0.75, L 2 -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 2 -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 10 (a i /R ⊙ ) = 1.37-1.68and 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.935and log 10 (a i /R ⊙ ) = 1.38-1.77are stripped in Case-A, Case-AB and Case-ABC mass-transfer phases and avoid/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 10 (a i /R ⊙ ) ≈ 1.78-2.80where 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 2 -overflow and thus avoid runaway mass transfer.Closer inspection shows that in these systems, the radius of the primary star exceeds R L 2 by ≲ 20% for ≲ 10 3 yrs.It is uncertain whether these contact binaries (see Fig. 2.2c) are stable, given that the phase of L 2 -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 (Y.Götberg, 2022, 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 3 yrs) 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. 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 10 (a i /R ⊙ ) ≲ 2.5, and q i ≤ 0.15 and log 10 (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 10 (a i /R ⊙ ) ≈ 1.72 and log 10 (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 2 -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 10 (a i /R ⊙ ) ≈ 3.4-3.5where runaway mass transfer does not occur, but a classical CE phase is expected from L 2 -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 ⊙ .

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) 10 .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 i11 .
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 10 (a i /R ⊙ ) ≲ 0.57 and q i = 0.35-0.55,and log 10 (a i /R ⊙ ) = 0.57-0.79and 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 2 -overflow.
Case-A systems with q i = 0.65-0.94and log 10 (a i /R ⊙ ) = 0.73-0.94evolve 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.0M ⊙ .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.00that 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.00and log 10 (a i /R ⊙ ) = 0.91-1.45evolve into contact binaries because of the expansion of the accretor.The donor stars in these binaries have deep convective envelopes, resulting in masstransfer rates ≳ 10 −5 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 10 (a i /R ⊙ ) < 0.91, which avoid contact, have mass transfer rates ≲ 10 −7 M ⊙ /yr, resulting in a nuclear-timescale expansion of the accretor star.The three models at q i = 0.85-1.00not 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 10 (a i /R ⊙ ) = 1.07-2.64),Case-Bl systems with q i > 0.65-0.75avoid 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 10 (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.45and log 10 (a i /R ⊙ ) = 2.38-2.56go through a phase of runaway mass transfer and experience L 2 -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 10 (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 10 (a i /R ⊙ ) = 2.39-2.60only 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.

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. (2010Ge et al. ( , 2015Ge et al. ( , 2020)).It is found that donors with radiative or convective envelopes with core-mass fractions12 greater than 0.5 shrink and hence have stable mass transfer (Temmink et al. 2023) 13 .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 coremass 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.0M ⊙ 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.2R ⊙ 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 1 < R RL, 1 at q = 1.21.At M 1 = 7.2 M ⊙ (q = 0.87, age = 25.8815Myr) 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.

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.Fig. 11 shows the results for the whole mass range of M 1, i ∈ [4.8; 20.8] M ⊙ 14 .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.97range (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 2 -overflow" and "No contact" (see Sect. 4.1.1),which correspond to the outcomes indicated by the background colours in Figs.7-8.Models expected to avoid contact but which encountered numerical issues before the primaries reached core-C exhaustion ("No contact" hatching in Figs.7-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-8 ("L 2overflow" 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 14 To compute p birth , the initial mass function is integrated from M l to M u (see Sect. 2.5).For the models with M 1, i = 20.0M ⊙ , M u = 20.8M ⊙ .This explains why the initial primary mass range extends to 20.8 M ⊙ .
shows the lower limits of the incidences for stellar mergers and classical CEs (see Sect. 4.4).
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.

Stellar merger and classical CE incidence
We assume that all binaries experiencing runaway mass transfer and/or L 2 -overflow merge or enter a classical CE phase.Following the physical picture described in Sect.3, we make the 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", "L 2 -overflow" and "No contact"), while the outer level shows the ancillary outcome ("L 2 -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 q i 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.
distinction based on the structure of the primary star.This means that Case-A and -Be binaries with runaway mass transfer and/or L 2 -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-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 2 -overflow in later contact phases (so far if more than one contact phase occurred, we only consid-  ered 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 nonaccreted 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 1 ≳ 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.0M ⊙ .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).

Comparison with observed contact binaries
In Fig. 12, our simulated population of Case-A contact binaries with initial primary masses of 4.8-20.8M ⊙ is compared to observed near-contact 15 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 15 Near-contact systems are systems in which for both components R/R RL ≥ 0.9 (Menon et al. 2021).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 2 -overflow, while the "accretor expansion + L 2 -overflow" do.The third contribution comes from runaway mass transfer systems ("runaway MT").The "accretor expansion" and "accretor expansion + L 2 -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).  12. Probability density function (PDF) of contact systems formed in Case-A binaries with M 1, i = 7.9-20.8M ⊙ 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 M 1 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 2 -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 2 -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 2 -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 2 -overflow.As described in Sect.3.3, this can lead to considerable orbital angu-lar 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 2 -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 2 -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.

Discussion
In this section, we put into perspective the influence of masstransfer 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).

The role of mass-transfer efficiency
The mean mass-transfer efficiencies β of binaries in our grid (Fig. F.1 in Appendix F) 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 Table 1.Stellar merger and classical CE incidences, and critical mass ratios q i, crit for mass-transferring binaries with M 1, i ∈ [4.8; 20.8] M ⊙ .For each initial primary mass, the incidence fractions are given per mass-transfer case.a) Case A and -Be systems with q i < q i, crit form contact binaries through runaway mass transfer and/or L 2 -overflow, and merge.Case-Bl and -C systems form classical CEs through runaway mass transfer and/or L 2overflow.
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 in Appendix C), 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.2M ⊙ the mass-mass transfer efficiency is consistent with the observationally constrained value of ∼30%.
Using observations of the classical Algol system δ Librae, Dervis ¸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 (Eq.6 and 7).From estimating the time-averaged mass-transfer efficiencies over the whole Case-A phase (thermal-and nucleartimescale 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.

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. (2010Ge et al. ( , 2015Ge et al. ( , 2020) ) 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 (H.Ge, 2023, priv. comm.).Picco et al. (2023) use the adiabatic mass-radius exponents ξ R, ad (see Sect. 3.5) from Ge et al. (2020) 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.

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 2 -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 2 -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) has found a significant close-binary fraction of 34 +8 −7 % for one of the clusters analysed in Milone et al. (2018), NGC 330.

The onset of contact through runaway mass transfer and L 2 -overflow
In Sect.4.1, we find that runaway mass transfer and/or L 2overflow 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 2 -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 2 -overflow.The added effect of the orbital angular momentum loss associated with L 2 -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 2 -overflow often occur quasi simultaneously.In these systems, the L 2 -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 2 -lobe.In both cases, we expect the onset of a classical CE phase.

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 2 -overflow, orbital decay because of tides, and non-conservative mass transfer.We find that accretor expansion, L 2 -overflow, and runaway mass transfer lead to the formation of contact binaries in Case-A and -Be systems, and L 2 -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 blow 2 M ⊙ also form contact binaries because of the orbital decay caused by tides.Overall, the incidences of masstransferring 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 masstransfer 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 2 -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 masstransferring 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 2 -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 2 -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.Contact tracing results for M 1, i = 1.9-6.1 M ⊙ .Models with M 1, 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 M 1, i = 1.9-2.2M ⊙ we find certain models where mass transfer starts when the primary is on the WD cooling track and experiences sudden radial expansion.Secondary's evolutionary state at contact or termination.

Fig. 1 .
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 = R RL .
as τ local KH (m) = M m c P (m ′ )T (m ′ )dm ′ , with c P (m) the heat capacity at constant pressure P and T (m) the temperature at mass coordinate m.In practice, timescale comparisons should therefore be made at an order-of-magnitude level.

Fig. 2 .
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 L 1 -and L 2 -points, respectively.The filled grey-blue circles represent the stellar cores in Panel (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 L 2 -overflow (2c), which likely results in a stellar merger.The scenarios (1b-2b-2c) and (2a-2b-2c) correspond to the "L 2 -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.

Fig. 3 .
Fig.3.Example of a M 1, i = 10.2 M ⊙ , q i = 0.4 and a i = 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 R RL (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.

Fig. 9 .
Fig.9.Example of a M 1, i = 10.16M ⊙ , q i = 0.6 and a i = 1398.2R ⊙ 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.

Fig. 11 .
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.
Fig. 12. Probability density function (PDF) of contact systems formed in Case-A binaries with M 1, i = 7.9-20.8M ⊙ 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 inMenon 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 M 1 are indicated with a star symbol.
Fig. A.1.Same as Fig.3but for a M i = M ⊙ , q i = 0.9 and a i = 12.5 R ⊙ binary.
Fig. B.2. Contact tracing resultsfor M 1, i = 1.9-6.1 M ⊙ .Models with M 1, 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 M 1, i = 1.9-2.2M ⊙ 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 M 1, i = 0.8-1.3M ⊙ .For M 1, i = 1.1-1.3M ⊙ findcertain models where mass transfer starts when the primary is on the WD cooling track and experiences sudden radial expansion.
thermal pulses( j)    Numerical issues (k) Primary's evolutionary state at contact or termination.MS: before core-H exhaustion.Post-MS: after core-H exhaustion and before core-He ignition.CHeB: after core-He ignition and before core-He exhaustion.Post-CHeB: after core-He exhaustion.(l) Occurrence of contact phases for models with initial primary masses M 1, 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 R 1 = R RL, 1 .
Table G.1.Extract of the table with the contact tracing results of all 5957 binary MESA models.The full table is available online at https://zenodo.org/doi/10.5281/zenodo.10148634.