© The Authors 2023

1. Introduction

A significant fraction of stars exists in multiple systems where the components will interact (Abt & Levy 1976; Abt 1983; Duquennoy & Mayor 1991; Sana et al. 2012; Moe & Di Stefano 2017). As such, mass transfer is a crucial process in the evolution of many stellar systems (e.g. Morton 1960; Paczyński 1971), including for the formation of white-dwarf mergers (Iben & Tutukov 1984), core-collapse supernova diversity (Podsiadlowski et al. 1992), and binary gravitational-wave sources (e.g. Abbott et al. 2016, 2017; Amaro-Seoane et al. 2022). Whether the onset of Roche-lobe overflow (RLOF) in a particular binary leads to stable or unstable mass transfer produces a qualitative difference in the future evolution, as well as in the properties of the final remnants. Where unstable mass transfer and the resulting common-envelope (CE) phase typically result in close binaries or a single, merged object (e.g. Paczynski 1976), stable mass transfer tends to produce wider binaries (e.g. Soberman et al. 1997). Different assumptions regarding mass-transfer stability therefore strongly affect predictions for, for example, double-white-dwarf (DWD) binaries and their associated transients, such as thermonuclear supernovae and gravitational-wave transients. These assumptions also affect transients linked to non-compact stars, such as luminous red novae, which are thought to be an observational manifestation of unstable mass transfer (Tylenda & Soker 2006; Tylenda et al. 2011; Ivanova et al. 2013). Despite its clear importance, only a few systematic studies of the stability of mass transfer, which we discuss in more detail below, have been performed over the years.

As is well known, the response of a star to mass loss depends on the structure of its envelope. The canonical view (e.g. Webbink 1985) is as follows. Donor stars with radiative envelopes tend to shrink on their dynamical timescales in response to a sudden loss of mass. Therefore, mass transfer involving this type of donor star is expected to be relatively stable. If a donor star with a convective envelope loses mass, its radius tends to increase, unless the core-mass fraction exceeds 0.5. On the other hand, its Roche radius typically decreases if the donor star is more massive than the accretor. This means that such a donor star will overfill its Roche lobe by an ever-increasing amount, eventually leading to mass transfer on a dynamical timescale and the formation of a CE. Indeed, classical studies (such as Paczyński 1965; Paczyński et al. 1969; Webbink 1985; Hjellming & Webbink 1987) found that mass transfer from a red giant donor star to a less massive companion should be unstable for nearly the entire giant branch.

However, this classical picture relies on simplifying assumptions, such as the approximation of the stellar structure by simple (composite or condensed) polytropes, and the assumption that the donor star responds to the loss of mass fully adiabatically. It has long been known that this is a serious limitation (see e.g. Paczyński 1965), and these results should be interpreted and applied with care. In fact, more recent theoretical studies (see e.g. Han et al. 2002; Podsiadlowski et al. 2002; Woods & Ivanova 2011; Passy et al. 2012; Pavlovskii & Ivanova 2015) found stable mass transfer for a wider range of mass ratios and orbital periods than the classical results mentioned above. These authors argued that the ability of the outer layers of giant stars to thermally readjust on dynamical timescales should significantly increase the stability of mass transfer compared to the classical results. Unfortunately, these studies typically focus on only a small part of the parameter space. An exception to this is the collection of recent work of Ge et al. (2010, 2015, 2020), which does map the stability boundary systematically for a large set of masses, periods, and mass ratios. While the authors use a realistic equation of state to model the stellar structure, this set of studies still employs the adiabatic approximation mentioned above.

Additionally, over the years, different physically motivated criteria have been proposed to identify unstable mass transfer in both detailed 1D simulations and rapid binary population synthesis (BPS) methods. Indeed, there are multiple physical reasons why mass transfer might become unstable (see e.g. Ivanova et al. 2020). These reasons include dynamical timescale evolution (which has become a common assumption, as in the references above) and the overflowing of an outer lobe (OL) in the Roche potential. As we explain in more detail in Sects. 3 and 4, these criteria typically do not identify unstable mass transfer consistently; mass transfer can still be stable and self-regulating even when these criteria are met.

Not only the definition of reliable criteria for stable mass transfer poses a potential problem for the results mentioned above. As outlined in, for example, the discussion of Podsiadlowski et al. (1992), there are several observed systems that are in tension with the predicted limits for stable mass transfer from polytropic adiabatic calculations. This includes systems that appear to have experienced mass transfer from donors on the red giant branch (RGB; case B) or asymptotic giant branch (AGB; case C) but also have relatively long orbital periods (see e.g. Eggleton & Tout 1989). This suggests that this phase of mass transfer was likely stable, whereas classical results would have predicted otherwise. Additionally, Podsiadlowski et al. (1992) note that several massive systems are observed to be in such a phase of late case B or early case C mass transfer (e.g. VV Cep; Wright 1977, AZ Cas; Cowley et al. 1977, and KQ Pup; Cowley 1965) but have apparently avoided the dynamical instability predicted by the classical results above.

Further populations of observed systems have long been known to indicate problems with classical, simplified results regarding the stability of mass transfer. For example, Tauris & Savonije (1999) and Podsiadlowski et al. (2002) found that the increased mass-transfer stability from their detailed stellar calculations, as compared to classical stability criteria, helps explain the properties of low-mass X-ray binaries and their descendants. Hot subdwarfs formed by stable mass transfer from a giant star onto a lower-mass companion also naturally suggest mass transfer that is more stable than predicted by the classical results (as noted and investigated by e.g. Han et al. 2002, 2003; Vos et al. 2019, 2020). Stars that have become blue stragglers after stably accreting from a giant companion also motivate consideration of increased mass-transfer stability (e.g. Chen & Han 2008). Leiner & Geller (2021) constructed multiple synthetic blue straggler populations using a range of different mass-transfer stability prescriptions. They find that none could reproduce the observed Gaia sample adequately and stressed that a better understanding of the stability of mass transfer is needed to explain observations that point towards more stable mass transfer. Woods et al. (2012) were able to reproduce a subset of the observed DWD population by assuming the first phase of mass transfer is stable. Despite the fact that it has long been clear that assumptions made in rapid BPS methods tend to underestimate the stability of mass transfer, in particular the stability boundaries of Hjellming & Webbink (1987) and Hurley et al. (2002), they still remain widely in use today.

In this work we present a systematic study of mass-transfer stability for RLOF from donor stars with zero-age main-sequence (ZAMS) masses in the range 1 M_⊙ ≤ M_d, ZAMS ≤ 8 M_⊙. We simulate the evolution of 1404 binaries with the Modules for Experiments in Stellar Astrophysics (MESA) code (see Paxton et al. 2011, 2013, 2015, 2018, 2019) in order to obtain the boundary between stable and unstable mass transfer for stars that fill their Roche lobes as post-main-sequence (post-MS) stars, without assuming the donors respond to the loss of mass adiabatically. In our simulations we only find such an adiabatic response when the mass-transfer rate exceeds a critical value, which is several orders of magnitude higher than the rate of mass loss occurring on the Kelvin-Helmholtz (KH) timescale of the donor (see Sect. 3.3) We introduce a criterion for stable mass transfer based on this transition to complete adiabatic evolution (our ‘quasi-adiabatic’ criterion) and use it to identify the boundary between stable and unstable mass transfer. Furthermore, we also find alternative boundaries assuming two commonly used stability criteria, which are based on dynamical-timescale evolution and the overflowing of its OL by the donor star. We treat the accreting star as a point mass and therefore do not follow its evolution.

The contents of this work are organised as follows. In Sect. 2 we describe our numerical method and assumptions. In Sect. 3 we describe the three different physically motivated instability criteria we employ, and argue why we base our main results on the quasi-adiabatic criterion (Sect. 3.3). We present our results in Sect. 4.1, where we focus on evolutionary aspects of the mass transfer, and Sect. 4.2, where we focus on the critical mass ratios for stable mass transfer. We critically analyse the limitations of our simulations and results in Sect. 5, where we also compare our theoretical results to other results from the literature (Sect. 5.1) and to observations of post-mass transfer binaries (Sect. 5.4). We end by summarising our conclusions in Sect. 6.

2. Numerical methods and physical assumptions

We used MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), version r12115, to simulate the evolution of single stars with masses of 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, and 8.0 M_⊙, as well as to compute a grid of binary models. In this section we describe our most important assumptions for the single-star models, followed by a description of our binary evolution calculations. For more details regarding our MESA settings, we refer the reader to the MESA inlists (input files) used in this work¹. The computed single and binary star models are available on request, as are the other data used in this work².

We adopt a slightly super-solar metallicity of Z = 0.02 (see e.g. Asplund et al. 2009; Magg et al. 2022). For the rest of our settings and assumptions, we generally follow the solar-calibrated assumptions described in Choi et al. (2016). We use the Ledoux criterion (Ledoux 1947) for determining whether or not a region inside the stars is stable against convection, and use a value for the mixing length parameter of α_MLT = 1.82 for convective regions. We also include the effects of semi-convection (implemented as in Langer et al. 1983), with α_semi = 0.1. We treat convective overshooting of stellar cores and burning shells as in the solar calibrated models of Choi et al. (2016), which were calibrated by matching models to the main-sequence turnoff in the open cluster M 67. We employ MESA’s default basic nuclear network to model nuclear energy generation inside the donor stars, which is sufficient for our purposes and for the parameter space we consider. We also include mass loss through stellar winds: we assume a modest ‘Reimers’ wind (Reimers 1975) with strength η_R = 0.1 on the RGB, and a ‘Blöcker’ wind (Bloecker 1995) with strength η_B = 0.2 on the AGB. We furthermore increase the resolution of the numerical mesh with respect to the default settings around the boundaries of convective zones, as well as in regions with strong H or He abundance gradients, which typically results in around 1500 zones. We end our single-star simulations when the first thermal pulse (TP) occurs.

We use the resulting single-star tracks to construct a grid of input models for our binary stellar simulations. For each track, we save models at pre-determined radii, which we then use as starting points in our MESA binary simulations. The points along the evolutionary tracks where models are saved are chosen by requiring a minimum radial resolution of Δlog R_d = 0.2 across the entire track. We furthermore ensure that the smallest radius chosen for each evolutionary phase is larger than the largest stellar radius reached during the previous phase. For stars with masses < 2 M_⊙, we do not consider the AGB phase, as we find that these stars do not expand significantly beyond their maximum radius on the RGB before the first TP occurs. The above selection criteria result in about 11 points in total for each donor mass. Figure 1 shows our single-star evolutionary models in the Hertzsprung-Russell diagram (HRD), as well as the aforementioned points where we initialise mass transfer in binaries.

Fig. 1. Single-star evolutionary tracks in the HRD for different stellar masses. Each colour corresponds to a given initial mass, as indicated to the lower left of each line. The coloured open circles indicate where we save single-star models for use in our binary star simulations. We have indicated lines of constant radius using grey dashed lines. Note that we have truncated these tracks to start at the ZAMS, and end at either the first TP or the maximum radius of the star, depending on which occurs first.

For each of those saved models, we initialise 12 binary systems, with mass ratios (q = M_a/M_d; throughout this work, the subscripts ‘a’ and ‘d’ will refer to the accretor and donor, respectively) ranging from 0.1 to 1.2 in increments of 0.1, for a total of 1404 binary star models. For each binary, we thus set M_a = qM_d, where M_d is the actual mass of the donor star (which is not necessarily equal to its initial mass) at each of the pre-determined points, including mass lost due to winds. We purposely include mass ratios > 1 in our grids to ensure we always enclose the stability boundary within our grid points. Previous studies (see the references in Sect. 1) indicate that for a sizeable part of our parameter space, the critical mass ratio (above which mass transfer is stable) could be above unity. We chose the initial orbit of the binary such that the donor star fills its Roche lobe (with equivalent-volume radius R_L1) at the pre-determined stellar radii (see above). To that end, we used the fitting formula from Eggleton (1983, their Eq. (2)) to calculate the size of the initial orbits, given the radii of the donor stars from our saved single-star models:

$$\begin{array}{c}\hfill R_{{\mathrm{L}}_1}/a\approx \frac{0.49q^{-2/3}}{0.6q^{-2/3}+ln(1+q^{-1/3})},\end{array}$$(1)

where a denotes the binary orbital separation. We note that we make the simplification that the initial eccentricity is equal to 0.

We generally use the same settings and configurations for the donor star in the binary system as for the single-star simulations described above. Additionally, we use the ‘Kolb’ scheme (Kolb & Ritter 1990) to calculate the mass-transfer rate. As previously mentioned, we do not evolve the accreting star, and take it to be a point mass. This allows us to study the stability of the donor stars separately from the evolution of the accretor. Additionally, simulating accreting stars in binaries requires special attention to the proper treatment of many quantities, such as the specific entropy of the accreted material and the spin of the accretor, which is outside the scope of this work (but see Sect. 5).

We furthermore assume that all the material lost by the donor is accreted by its companion, that is, the mass transfer phases are fully conservative in terms of both mass and angular momentum. This assumption cannot be expected to hold in all or even most circumstances. The effect of non-conservative RLOF on the stability of mass transfer strongly depends on assumptions regarding the loss of angular momentum. For example, in the case where any material that is not accreted is assumed to leave the system with the specific angular momentum of the accretor (isotropic re-emission), it can be shown analytically that an increase in the fraction of mass leaving the system (β) results in a decrease of dlnR_L1/dlnM_d and dlna/dlnM_d (see e.g. Soberman et al. 1997; Tauris & Savonije 1999; Tauris 1996), making mass transfer more stable. Hence, our results should provide an upper limit (in terms of the critical M_a/M_d) on the location of the boundary between stable and unstable mass transfer in the parameter space, that is, mass transfer is expected to be more stable than found here if mass were lost from the system under the assumption of isotropic re-emission (but see Sect. 5.2 for further discussion).

In this work, we do not model stellar rotation and thus do not include effects of spin-orbit coupling on the evolution of the binaries and the stability of mass transfer (but see Misra et al. 2020, for an extensive study). Hence, we also do not keep track of the possible occurrence of Darwin instabilities (Darwin 1879), which may become relevant for extreme mass ratios M_a/M_d ≲ 0.1 (but see Sect. 5.2).

We focus on the first phase of mass transfer and end our calculations when the donor star has shrunk well inside its Roche lobe, either after its envelope has mostly been exhausted or after it has started He burning in its core. As such, in this work we do not consider subsequent episodes of mass transfer that could occur when partially stripped giants re-expand to fill their Roche lobes after a core He-burning phase.

Having simulated the grid points mentioned above, we systematically decrease the uncertainty on the location of the stability boundary resulting from our main criterion (see Sect. 3) by using a simple bisection method. We execute sufficient bisection iterations such that the relative uncertainty on the critical mass ratio is 5% or less. For the other two criteria, we interpolate linearly in our grid points to estimate the critical mass ratios. Finally, we systematically verify that our results are converged by re-simulating a representative selection of stable and unstable systems with increased resolution in both time and mesh spacing. We discuss these tests in more detail in Appendix A.

3. Criteria for unstable mass transfer

A classical approach (e.g. Soberman et al. 1997) for determining whether or not mass transfer is stable, is to compare the adiabatic response to mass loss of the donor star radius ζ_ad ≡ d log R_d/d log M_d to the response of its Roche lobe (ζ_L ≡ d log R_L/d log M_d). If the radius of the donor star persistently increases relative to its Roche lobe during the mass transfer, a runaway situation develops and the mass transfer is assumed to become unstable. Ultimately, a CE that engulfs both stars may form (Paczynski 1976). However, if the donor star can remain roughly equal in size or even shrink relative to its Roche lobe, the mass transfer is stable and self-regulating. This comes with the caveat that, even if at any point during the mass transfer ζ_ad < ζ_L briefly, mass transfer does not necessarily become unstable, provided that the donor star recovers and shrinks relative to its Roche lobe again. As such, one would need to consider the entire evolution of the mass transfer, and it is not always possible or practical to accurately predict the outcome of the mass transfer at the onset using the ζ method. Well-known examples of this are so-called delayed dynamical instabilities (see e.g. Hjellming & Webbink 1987; Han & Podsiadlowski 2006; Pavlovskii & Ivanova 2015; Ge et al. 2020), where mass transfer is initially stable (according to the ζ method), but only turns unstable after some mass has been lost. Alternatively, it is common practice to define an effective critical mass ratio q ≡ M_a/M_d, which is based on the full mass-transfer history. If the mass ratio at the onset of RLOF is larger than these critical values, mass transfer is stable. Due to it being an effective value, delayed dynamical instabilities, where the masses become more similar before the instability, are naturally taken into account. However, any set of critical mass ratios is only valid under a single set of assumptions regarding mass-transfer efficiency and angular momentum loss.

There are also difficulties of a numerical nature. For example, it is insufficient to assume that numerical breakdown of a 1D stellar-evolution code during mass transfer indicates physical instability. Multiple physically motivated criteria have been proposed to identify unstable mass transfer. Indeed, there are multiple physical reasons why mass transfer might become unstable (see e.g. Ivanova et al. 2020). In this section we describe different criteria that aim to identify unstable mass transfer. We later apply these different criteria to our simulations and, in Sect. 4.2, compare the resulting critical mass ratios.

3.1. Dynamical-timescale mass transfer

A common way to determine whether RLOF is stable is by assuming that when the mass-transfer evolution accelerates to a dynamical timescale, such runaway mass transfer leads to the formation of a CE that engulfs both components of the binary. In this work, we compare the relevant timescales of both the donor star and the system as a whole to the binary period P, which we take as the characteristic timescale of the binary. We then consider a phase in binary evolution to be ‘dynamical’ if the fractional amount by which either the mass of the donor star or the binary orbit changes over one orbital period exceeds a certain threshold. Specifically, we check whether

$$\begin{array}{c}\hfill max(|{\dot{M}}_{\mathrm{d}}/M_{\mathrm{d}}|,|\dot{a}/a|)·P>A_{\mathrm{dyn}},\end{array}$$(2)

where A_dyn is a constant. We adopt A_dyn = 0.05, which is in agreement with the start of the plunge-in phase as described in Ivanova et al. (2020). We then define the critical mass ratio for dynamical mass transfer, q_dyn, as the lowest value of M_a/M_d at the onset of RLOF for which the condition in Eq. (2) is not met. We note that dynamical evolution does not always lead to runaway mass-transfer rates, and hence unstable mass transfer. In many cases, we still find self-regulating mass transfer phases for dynamically evolving binaries (see Sect. 4.1).

We find that our results for q_dyn are sensitive to the precise value of A_dyn, but not in equal measure over the full parameter space. The sensitivity to the value of A_dyn typically increases with decreasing initial donor mass or in more extended giants (see Sect. 4.2 for more details).

However, we consider that it would be misplaced to fixate too strongly on the precise choice of A_dyn. Once the binary is evolving on a dynamical timescale, the assumptions used to derive the Roche geometry themselves break down. The Roche potential is derived in the co-rotating reference frame, while neglecting the Euler force ${\mathit{f}}_{\mathrm{Euler}}\propto \dot{\mathit{\omega }}\times \mathit{r}$, where ω is the angular velocity of rotation of the reference frame. Generally, this is an acceptable simplification. However, in the regime of dynamical evolution, the Euler force will compete with the other fictitious forces included in the potential, and so the Roche approximation is no longer valid. This is most easily recognised by considering the ratio of the Euler force to the centrifugal force:

$$\begin{array}{c}\hfill |{\mathit{f}}_{\mathrm{Euler}}/{\mathit{f}}_{\mathrm{centrifugal}}|\sim \left|\frac{\dot{\omega }}{{\omega }^2}\right|=\frac{3}{4\pi }|\dot{a}/a|·P.\end{array}$$(3)

With A_dyn = 0.05 in Eq. (2), the Euler force contributes by at most about 1% in systems that avoid dynamical evolution according to Eq. (2). Hence, once mass-transferring binaries are evolving dynamically, any results that assume the Roche geometry should be interpreted with care; this includes our own results.

3.2. Overflow of an outer Lagrangian point

Alternatively, one could consider the overflowing of the equipotential surface passing through one of the outer Lagrangian points by the donor star as the onset of a dynamical instability. A similar criterion was adopted in Paczynski (1976), the work commonly cited for the origin of the concept of CE evolution. Significant angular momentum loss through the outer Lagrangian points during mass loss (roughly the specific angular momentum of a circumbinary ring with radius ∼1.2a for mass loss through L₂; see e.g. Pribulla 1998) is thought to dramatically shrink a binary’s orbit and render the system dynamically unstable (see e.g. Nandez et al. 2014; Pavlovskii & Ivanova 2015). However, we stress that the loss of mass through either of the outer Lagrangian points is not guaranteed to always result in unstable mass transfer (see e.g. Bowler 2010, for an observational example of this). For example, under certain approximations, Linial & Sari (2017) find that any mass loss through the L₂ point occurs at a lower rate than conventional mass loss through the L₁ point and that its effects on the overall stability of a binary could be only marginal. Nevertheless, donor stars that extend past their OL, which passes through the L₃ point if M_a ≤ M_d and otherwise through the L₂ point, significantly violate the assumption of spherical symmetry made in 1D evolution codes, which makes reliably simulating this phase challenging.

We do not model mass loss through outer lobe overflow (OLOF), and hence ignore the effects thereof on the evolution of our binaries (but see Marchant et al. 2021; Bobrick et al., in prep.). However, we keep track of whether or not the donor star experiences OLOF by comparing its radius to the equivalent-volume radius of the OL, R_OL. To that end, we numerically calculated R_OL using Fortran routines for −4 ≤ log q ≤ 4, including a damped Newton-Raphson method to find the location of the outer Lagrangian point, and double Simpson method to calculate the volume of the potential surface passing through this point. The results we obtained are converged to better than one part in 10⁸. A Monte Carlo method written in Python was used to extend our calculations of R_OL beyond our Fortran routines to −7 ≤ log q ≤ 7. We fit the ratio of the equivalent-volume radius of the OL to that of the L₁ lobe with the following function:

$$\begin{array}{c}\hfill \frac{R_{\mathrm{OL}}}{R_{{\mathrm{L}}_1}}=1+\frac{0.441q^{-0.325}}{1+0.412q^{-0.8}}·\end{array}$$(4)

This fit is accurate to better than 0.2% over the full range of q values. Analogous to the previous section, we define a critical OLOF mass ratio q_OLOF as the lowest value of M_a/M_d|_RLOF for which the donor star does not experience OLOF during the entire mass transfer. Typically, we find that the overflow of outer Lagrangian points occurs in systems where the evolution becomes dynamical, meaning that q_dyn and q_OLOF have similar values (see Sect. 4.2).

3.3. Thermal readjustment of surface layers

Directly following a sudden loss of mass, the donor star loses both its hydrostatic and thermal equilibria. In the simplified picture, the initial response of the donor star occurs on the dynamical timescale, when it attempts to restore its hydrostatic equilibrium. After that, mass transfer is driven by the slower thermal readjustment of the donor, as the donor star attempts to recover the thermal equilibrium radius appropriate for its new mass on the global (KH) thermal timescale:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{KH}}=\frac{1}{2}G\frac{M_{\mathrm{d}}^2}{\mathit{RL}}\sim 1.5\times {10}^7\frac{M_{\mathrm{d}}^2}{M_{\odot }^2}\frac{R_{\odot }}{R}\frac{L_{\odot }}{L}\mathrm{yr}.\end{array}$$(5)

A very simple estimate for the rate of mass transfer occurring on a donor star’s global thermal timescale is then given by

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{KH}}\equiv -\frac{M_{\mathrm{d}}}{{\tau }_{\mathrm{KH}}}\sim -6.7\times {10}^{-7}\frac{M_{\odot }}{M_{\mathrm{d}}}\frac{R}{R_{\odot }}\frac{L}{L_{\odot }}{M}_{\odot }{\mathrm{yr}}^{-1}.\end{array}$$(6)

More specifically, in the simplified picture it is assumed that the donor star initially responds to the loss of mass adiabatically. With this adiabatic assumption, the new sub-surface layers in a radiative star (shown as the red dashed curve in Fig. 2a) will have lower entropy compared to the sub-surface layers before the loss of mass. Hence, after adiabatic decompression, these layers will occupy a smaller volume than the layers that have been removed, which stabilises the mass transfer and leads to relatively low critical mass ratios. In the envelope of a convective giant, approximated as isentropic, the new sub-surface layers (shown as the red dashed curve in Fig. 2b) will have similar entropy as before the loss of mass. This results in the sub-surface layers expanding along with the star as a whole in their response to the loss of mass (we note that even in the context of polytropic models, the conclusion that an isentropic envelope expands upon mass loss depends on the fractional core mass; Hjellming & Webbink 1987). This behaviour of these convective giants is reflected in their generally comparatively large critical mass ratios. Typically, the location of the boundary between stable and unstable mass transfer used in prescriptions used in rapid BPS studies (see e.g. Toonen et al. 2012; Claeys et al. 2014; Riley et al. 2022) is based on model calculations (such as Hjellming & Webbink 1987; Ge et al. 2010) assuming (part of) these simplifications.

Fig. 2. Schematic illustration of the entropy (S) profile versus mass coordinate m in various types of stars. Panel a depicts a representation of a star with a radiative envelope, whilst panel b shows a model for a giant with a convective, isentropic envelope. A more realistic representation of the envelope of a convective giant is shown in panel c, which includes the superadiabatic sub-surface layers (highlighted in green). In each panel, the black line represents the entropy profile of the undisturbed star. Dashed red lines represent the entropy profile after sudden adiabatic (see text) loss of ΔM of mass. The dashed blue line in panel c shows the entropy profile if the mass is lost at a rate below the critical mass-loss rate (Eq. (9)).

However, even on timescales much shorter than the global thermal timescale, stars do not necessarily react to the loss of mass fully adiabatically, since surface layers respond more rapidly than the star as a whole (see e.g. Osaki 1970; Podsiadlowski et al. 2002; Kippenhahn et al. 2012). As emphasised by Woods & Ivanova (2011) and Passy et al. (2012), the thermal readjustment of the surface layers of the donor star increases mass-transfer stability as compared to adiabatic models (see more below). The different thermal timescales on which the different sub-surface layers of a donor star will react depend on the local properties of these layers, for example the efficiency of convection, the degree of ionisation, and the specific energy content (see also Woods & Ivanova 2011; Passy et al. 2012; Pavlovskii & Ivanova 2015). Approximating the local thermal timescale at mass coordinate m as the time it would take a star to radiate away all the thermal energy content of the layers above m at the rate of its surface luminosity, the local thermal timescale of a given layer can be calculated as (e.g. Kippenhahn et al. 2012)

$$\begin{array}{c}\hfill {\tau }_{\mathrm{th}}(m)=\frac{1}{L}{\displaystyle {\int }_m^M}{c}_P(\stackrel{\sim }{m})T(\stackrel{\sim }{m})\mathrm{d}\stackrel{\sim }{m}.\end{array}$$(7)

We note that the above equation implicitly assumes that one can divide the stellar interior into an outer and an inner zone at any mass coordinate, and that the outer zone can thermally readjust independently from the inner zone. This assumption might be inappropriate, especially in a convective region, where the timescale of convective overturn could be shorter than the local thermal timescale, which could stabilize thermal perturbations on a shorter timescale.

We present runs of the local thermal timescale with the fractional exterior mass for different moments in the evolution of a 1.5 M_⊙ donor star in the upper panel of Fig. 3. Similar to the global (KH) thermal timescale, at any mass coordinate relative to the surface the local thermal timescale tends to decrease as the star expands and becomes more luminous. Additionally, the local thermal timescale decreases monotonically with the mass coordinate m. However, there is a part of the sub-surface layers where the stellar interior is unstable to convection, but energy transport via convection is inefficient. This results in temperature gradients that are significantly larger than the adiabatic temperature gradient, and these regions are hence known as superadiabatic regions (here defined as regions where ∇ − ∇_ad ≥ 0.1∇; with $\mathrm{\nabla }=\frac{\partial logT}{\partial logP}$). Within these superadiabatic sub-surface layers (highlighted by black shading in Fig. 3), the steep temperature gradients result in a relatively sharp decrease in the local thermal timescales with m in the upper parts of these regions, compared to above and below the superadiabatic layers. Donor stars do not respond effectively adiabatically to mass loss until a certain critical mass-transfer rate is reached, determined largely by these relatively steep slopes of τ_th(m) in the upper part of the superadiabatic regions.

Fig. 3. Thermal properties of a 1.5 M⊙ donor star at various stages in its post-MS evolution. The colours of the lines correspond to the moments the star reaches certain radii (as shown in Fig. 1), whose values are shown in the legend. Regions with grey shading are significantly superadiabatic (‘sad’; here defined as regions where ∇ − ∇ad ≥ 0.1∇; with $\mathrm{\nabla }=\frac{\partial logT}{\partial logP}$). The top panel shows the run of the local thermal timescale (see Eq. (7)) with the logarithm of the fractional exterior mass, such that the stellar centre is located on the left end of the x axis, and the surface is located to the right. The bottom panel shows the run of the associated thermal mass-loss rate (see Eq. (8)) on the same horizontal axis.

We can make this more explicit by examining the thermal mass-loss rate Ṁ_th(m) as a function of mass coordinate. This is the minimum mass-loss rate required to strip away the outer layers of a star down to this mass coordinate m faster than these layers can thermally readjust. This mass-loss rate can be estimated in the following way:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{th}}(m)=-\frac{M-m}{{\tau }_{\mathrm{th}}(m)},\end{array}$$(8)

which is plotted in the bottom panel of Fig. 3. It can be seen that Ṁ_th(m) generally reaches its maximum value

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{th},\mathrm{crit}}=max{\dot{M}}_{\mathrm{th}}(m)\end{array}$$(9)

at the top of the superadiabatic layer. As long as the mass-loss rate does not exceed Ṁ_th,crit, there is always a part of the outer layers that can thermally readjust on a timescale shorter than the mass loss timescale. Conversely, if the mass-loss rate should exceed this critical value, no layer would be able to thermally adjust on short enough timescales, and the mass transfer can be expected to take place fully adiabatically.

In the outer sub-surface layers, the local thermal timescales are shorter than the donor star’s global and local dynamical timescales. Hence, for mass-transfer rates up to Ṁ_th,crit, the initial response of a convective giant is dictated by the readjustment of the local thermal structure. In this scenario, the new sub-surface layers (shown as the blue dashed curve in Fig. 2c) have a lower entropy than before the mass loss. Thus, similar to what happens in a radiative star, the sub-surface layers occupy a smaller volume after decompression. This stabilises mass transfer over a significantly wider range of initial mass ratios than when the possibility of thermal relaxation is ignored. However, if the mass loss becomes too rapid (Ṁ_d > Ṁ_th,crit), the donor star is not able to readjust to the loss of mass on its local thermal timescale, and its response becomes adiabatic. In that scenario, the new sub-surface layers (shown as the red dashed curve in Fig. 2c) will have a higher entropy, compared to the ones before the mass loss. The sub-surface layers thus expand relative to the star as a whole, raising the mass-transfer rate even more and driving the mass transfer to instability.

As such, the transition to fully adiabatic response of the donor star is an indication for the onset of unstable mass transfer. We thus keep track of Ṁ_th,crit during the evolution of our binary systems, and whether or not this value is exceeded by Ṁ_d. We define a last, ‘quasi-adiabatic’ critical mass ratio q_qad as the lowest value of M_a/M_d|_RLOF for which we can simulate the complete mass transfer phase with MESA and the mass-transfer rate does not exceed Ṁ_th,crit.

In Fig. 4 we show how Ṁ_th,crit changes during the evolution of donor stars with different ZAMS masses. Like the global KH mass-transfer rate, Ṁ_th,crit typically increases with donor radius. This is a direct consequence of the local thermal timescale at any depth decreasing with the stellar radius. At a given donor radius, Ṁ_th,crit typically increases with donor mass, whilst Ṁ_KH typically decreases with donor mass. For the least evolved stars, Ṁ_KH is up to about three orders of magnitude lower than Ṁ_th,crit. However, since Ṁ_KH is a steeper function of R_d than is Ṁ_th,crit, the two rates grow more similar with increasing donor radius.

Fig. 4. Critical mass-loss rate as a function of stellar radius from thermal considerations (see Eq. (9)). Solid lines show the critical mass-loss rate based on local thermal properties of the sub-surface layers of the donor star (see Eq. (8) and Fig. 3), whilst dashed lines show the critical mass-loss rate based on global thermal properties (see Eq. (6)). The linear portions of the curves correspond to the giant phases. The colours of the lines correspond to different initial donor masses, as indicated to the lower left of each dashed line.

As briefly mentioned earlier, we do not find that dynamical or OLOF evolution always leads to runaway, unstable mass transfer. In fact, mass transfer is self-regulating in most of the systems with dynamical or OLOF evolution within the context of our models, provided they did not first encounter the quasi-adiabatic criterion. By contrast, we find that our quasi-adiabatic criterion, based on local thermal properties of the donor star, does consistently pick out the systems with runaway mass transfer. Hence, in this work, we base our main results on the quasi-adiabatic criterion described in Sect. 3.3. However, we also show the effects of choosing different criteria, as described in Sects. 3.1 and 3.2.

4. Results

In this section we present the results from our binary simulations. In Sect. 4.1, we focus on the evolution of the mass-transfer rates in our binaries, whilst we present the critical mass ratios for stable mass transfer in Sect. 4.2. Our main results are tabulated in Appendix B.

4.1. Mass-transfer rates and evolution

For most of the systems that we simulated, the overall evolution of the mass-transfer rate can be divided into two distinct parts (see Fig. 5 for a few examples). The first is a relatively brief ‘thermal’ phase, when the donor star attempts to recover thermal equilibrium and where the mass-transfer rate peaks. This is followed by a significantly longer phase, driven by the natural expansion of the donor star on its current evolution timescale, where the mass transfer occurs at lower rates. The mass-transfer rates depend on the parameters of the systems at the onset of RLOF, in that both more evolved and more massive donor stars will typically lose mass at higher rates. This trend is a result of the global thermal timescale decreasing with increasing luminosity and radius (Eq. (5)).

Fig. 5. Time evolution of mass-transfer rates in systems with Md = 2 M⊙ (left column) and Md = 5 M⊙ (right column). Each row corresponds to RLOF starting in a different phase of the donor stars’ lives (at the radius indicated in the top right of each panel). More specifically, panels a and b show HG donors, panels c and d show donors on their RGB, and panels e and f show donors on their AGB. The colours of the lines correspond to different mass ratios at the onset of RLOF, as indicated in the legend. We have indicated where, if appropriate, our quasi-adiabatic-, dynamical-, and OLOF-based criteria are met with open circles, squares, and triangles, respectively. The zero point of the time axis is defined as the first moment when the mass-transfer rate exceeds the rate at which wind mass loss occurs.

In the following, we discuss our results for Hertzsprung gap (HG) donors separately from those for giant donor stars. In this work, we use the term HG to refer to all stars evolving from the terminal-age main sequence to the base of the RGB. These stars typically have radiative envelopes and evolve on their thermal timescales. Unless explicitly stated otherwise, we also include in that term low-mass subgiants that expand on a much slower, near-nuclear timescale.

4.1.1. Hertzsprung gap donors

Due to the structure of a radiative envelope and a relatively long global thermal timescale, the initial response of a HG donor to mass loss is to shrink. Because such a donor star needs to expand to comparatively large radii to regain its thermal equilibrium, the mass transfer is driven on the thermal timescale of the donor.

Our interests lie predominantly in the first, rapid part of the mass transfer phase, since any potential instabilities will occur there. We use the mass-averaged mass-loss rate

$$\begin{array}{c}\hfill {⟨{\dot{M}}_{\mathrm{d}}⟩}_{\mathrm{M}}=\frac{1}{M_{\mathrm{d},\mathrm{i}}-M_{\mathrm{d},\mathrm{f}}}{\displaystyle {\int }_{M_{\mathrm{d},\mathrm{i}}}^{M_{\mathrm{d},\mathrm{f}}}}{\dot{M}}_{\mathrm{d}}(M_{\mathrm{d}})\mathrm{d}{M}_{\mathrm{d}},\end{array}$$(10)

as a representative mass-transfer rate. The subscripts ‘f’ and ‘i’ correspond to quantities evaluated at the final and initial moments in the mass-transfer evolution (here defined as the moments where Ṁ_d ≥ 10⁻¹⁵ M_⊙ yr⁻¹). The rapid thermal part is short compared to the subsequent mass transfer on the evolutionary timescale of the donor star, and hence the time-averaged mass-transfer rate would be dominated by the latter. The mass-averaged mass-loss rate, as defined above, does not suffer from this issue, as most mass is lost during this rapid phase.

Figure 6 compares ⟨Ṁ_d⟩_M to the KH rate (as given by Eq. (6)) for 2 M_⊙ and 5 M_⊙ donor stars. It is clear that donors that fill their Roche lobes as HG stars typically lose mass at rates that are significantly below Ṁ_KH, by up to a factor of 10. This weakly depends on the mass ratio, and is almost independent of the orbital period for the more massive donors that expand significantly during the HG. Similarly, the peak mass-transfer rates, whilst typically closer, are often below Ṁ_KH, though higher than ⟨Ṁ_d⟩_M by a factor of a few.

Fig. 6. Mass-transfer rates from our calculations compared to the global KH rate (as given by Eq. (6)). The top panel shows results for a 2 M⊙ donor star, whilst the bottom panel shows results for a 5 M⊙ donor star. In each panel, the coloured lines show the difference between the mass-transfer rate predicted by Eq. (6), (which assumes the rate is governed by the global KH timescale at the onset of mass transfer) and the mass-averaged mass-loss rate ⟨Ṁd⟩M (see Eq. (10)) from our calculations. Vertical dashed black lines indicate the base and the tip of the RGB.

For HG donors, the evolution of the mass transfer can be fairly complicated. Often, the rapid phase of mass transfer consists of two peaks instead of one. At the start of mass transfer, the initially radiative donor star shrinks slightly as the radius follows the Roche radius. However, shortly after the first peak in the mass-transfer rate, the convective zone in the outer layers of the donor star extends increasingly far down into the stellar interior. The mass-transfer rate climbs to a second peak because these convective outer layers can supply the required material faster than a radiative envelope. Typically, the second peak coincides with the moment the convective envelope has grown to its largest extent (see also e.g. Harmanec 1970; Van der Linden 1987).

For many M_d = 2 − 3 M_⊙ donors with q ≳ 0.7 that start transferring mass before the base of the giant branch, the aforementioned shrinkage is significant enough to temporarily suspend mass transfer altogether (see Fig. 5a). Typically, such interruptions last of the order of a few megayears, equivalent to a few per cent of the total duration of the mass transfer, or a few times the global thermal timescale of the donor stars at that point.

In several M_d ≤ 2 M_⊙ donors, the convective part of the envelope extends down deeply enough to leave a discontinuity in the abundance profile at the deepest point of the convective zone. When a star’s H-burning shell reaches this discontinuity, the star’s luminosity and radius decrease slightly. This process also occurs in single stars (e.g. Thomas 1967; Iben 1968; Christensen-Dalsgaard 2015), where it is well-known to result in the so-called ‘luminosity bump’ in star clusters (e.g. King et al. 1985; Salaris et al. 2002). For the aforementioned low-mass donors, the associated shrinkage can be significant enough to suspend the mass transfer until the donor re-expands and once more fills its Roche lobe (as already shown by Kippenhahn et al. 1967), of the order of tens of megayears later.

Thus, several of the 2 M_⊙ HG donors experience two interruptions in their mass transfer, as shown in Fig. 5a (see also e.g. Han et al. 2000). By contrast, for the more massive stars (M ≥ 4 M_⊙), mass transfer is not suspended at any moment during the mass transfer in our models. These stars hardly become convective during mass loss (the mass of the convective part of the outer layers is typically only ∼0.1 M_⊙), and therefore they do not undergo significant shrinkage induced by either structural re-arrangement or by dredge-up/mixing effects.

In Fig. 5 we have indicated with different symbols when during the evolution of our binaries our three instability criteria are met, if applicable. For HG donor stars, none of the systems we simulated experienced dynamical or OLOF evolution before Ṁ_th,crit was reached. By contrast, we do find that our quasi-adiabatic criterion (see Sect. 3.3) is consistently met for systems where the mass-transfer rate runs away. Furthermore, we find that mass transfer does not become unstable immediately after the onset of RLOF, but that the binaries instead experience a delayed instability (see e.g. Hjellming & Webbink 1987; Han & Podsiadlowski 2006; Pavlovskii & Ivanova 2015; Ge et al. 2020). The amount of time it takes for mass transfer to become unstable is generally a function of the mass ratio, with the typical delay time being shorter in systems with more extreme mass ratios.

We sometimes find systems with anomalous mass transfer behaviour along the boundary between stable and unstable mass transfer. During the rapid thermal part of the mass transfer, Ṁ_d oscillates by up to three orders of magnitude, with an oscillation period similar to the typical critical thermal timescale during that phase. During several of the peaks in the mass-transfer rate, the rate at which the donor star loses mass exceeds the critical thermal rate (see Sect. 3.3). In these moments, the outermost layers are stripped away quasi-adiabatically (as the surface layers readjust on a longer timescale than the timescale of mass loss), and the mainly radiative regions below the convective outer layers briefly approach the surface. The donors recover as the mass-transfer rate decreases and the convective part of the outer layers regrows, restarting the cycle anew. After a few of these cycles, the mass-transfer rate stabilises and follows the typical evolution described in this section.

4.1.2. Giant donors

Donor stars that fill their Roche lobes as giants typically have convective outer layers. Convective envelopes are effectively isentropic (excepting the non-negligible superadiabatic region) and therefore tend to expand relative to their Roche lobe upon mass loss, resulting in higher, but shorter-lasting peaks in the mass-loss rates. More evolved giant donors have shorter evolutionary timescales and therefore higher mass-transfer rates, as can be seen in Fig. 5. For the majority of giant donors in our grid, we find that general properties of the mass-transfer evolution (such as duration and typical or peak mass-transfer rates) are poorly captured by Eq. (6). This is because this type of formula does not depend on properties of the binary systems, such as the mass ratio or the amount by which the donor overflows its Roche lobe, but only on global properties of the donor star (see also e.g. Langer et al. 2000, who find a similar mismatch).

Unlike for HG donor stars, ⟨Ṁ_d⟩_M for giant donors is strongly dependent on the initial mass ratio and period of a binary (see Fig. 6). This is especially true for the low-mass giants (M_d ≤ 2.5 M_⊙), where ⟨Ṁ_d⟩_M can range from over two orders of magnitude below to almost two orders of magnitude above the global thermal rate of the donors. For the more massive giants, the dependence on the initial mass ratio is less extreme. The peak mass-transfer rate follows a very similar trend with mass ratio and period, being about an order of magnitude higher than ⟨Ṁ_d⟩_M for unevolved giants, with the difference decreasing for more evolved giant donors. Over the full mass range, the mass-averaged mass-transfer rate tends to agree better with the global thermal rate, and depends less strongly on mass ratio, with increasing donor radius. As we describe in more detail below, even in the systems for which ⟨Ṁ_d⟩_M ≫ Ṁ_KH, mass transfer is still self-regulating as long as the critical mass-transfer rate (see Sect. 3.3) is not exceeded.

As for the less evolved donors discussed in Sect. 4.1.1, we have marked in Fig. 5 when the binaries with giant donors meet our instability criteria, if applicable. As Fig. 5f clearly illustrates, neither dynamical-timescale nor OLOF evolution necessarily lead to unstable, runaway mass transfer. As we noted for the HG donors, the mass transfer does typically turn unstable when the critical thermal rate (as given in Eq. (9)) is exceeded. As such, the dynamical and OLOF criteria should be interpreted more as guidelines for when our 1D simulations start becoming unrealistic, rather than actual physical stability criteria. This is discussed further in Sect. 4.2.

Typically, stable mass transfer for giant donors continues until nearly the entire envelope has been stripped. However, for relatively unevolved RGB stars where helium is ignited under non-degenerate conditions (in our grid this corresponds to M_d ≥ 2.5 M_⊙; log T_eff, d ≳ 3.6) with mass ratios that are not too extreme (M_a/M_d ≳ 0.8), the donor star can start burning helium in its core during mass transfer when there is still a significant amount of mass left in the envelope. As a result, after shrinking within their Roche lobes, which suspends the mass transfer, they evolve through the HRD in a bluewards loop within the HG region, without becoming hotter than the main-sequence band. The donor stars can re-expand enough to refill their Roche lobes as partially stripped AGB stars. The typical amount of mass left in the envelope is around 30−40% of the total donor mass at that point.

Lastly, we find that the most evolved AGB donors in our grid can experience TPs whilst transferring mass. As a result of recurring He shell flashes, these donor stars experience cycles of expansion and contraction. Consequently, also the mass-transfer rate increases and decreases significantly during these cycles, as visible in Figs. 5e and f. For many donor stars, the shrinkage is significant enough to completely suspend mass transfer temporarily. Additionally, for several systems, we find that the expansion is large enough to drive mass transfer to higher rates than during the rapid peak after the onset of RLOF.

4.2. Critical mass ratios

In this section we present our results for the critical mass ratios introduced in Sect. 3. For mass ratios greater than these critical ratios (M_a/M_d > q_crit), mass transfer was found to be stable according to one of the three criteria we applied. In Fig. 7, we present the critical mass ratios resulting from these criteria for two representative ZAMS masses 1 M_⊙ and 4 M_⊙.

Fig. 7. Critical mass ratios (above which mass transfer is stable) for 1 M⊙ and 4 M⊙ donor stars as a function of donor radius, under different criteria for stable mass transfer (see Sect. 3). The solid black line shows the critical mass ratios resulting from our quasi-adiabatic criterion, with the grey shading indicating the parameter space where mass transfer would be unstable. The red lines show the critical mass ratios resulting from the dynamical evolution criterion with Adyn = 0.01, 0.05, 0.25, as indicated to the right end of the lines. The solid orange line shows the critical mass ratios resulting from the OLOF-based criterion, and the solid blue line and hatching show where the Darwin instability could be encountered (see Sect. 5.2). The solid horizontal grey line shows where Ma/Md = 1, and the dashed vertical grey lines show the end of the HG phase (panel a and left line in panel b) and the tip of the RGB (panel b).

For HG donor stars, we find that any mass transfer instabilities set in well before their envelopes have become even marginally convective. Hence, the response of these stars to mass loss is essentially that of a donor star with a completely radiative envelope. Therefore, mass transfer from HG donors is relatively stable (see Sect. 3.3), with the typical quasi-adiabatic critical mass ratio being around q_qad = 0.25, which agrees well with values typically assumed for such donors (see e.g. Claeys et al. 2014). Our 1 M_⊙ model with the smallest Roche radius, which is classified as HG, behaves more like a convective donor star (discussed in more detail below), due to the presence of a non-negligible surface convection zone. As a result, this subgiant donor has comparatively larger q_qad than the HG donor stars discussed above. We were unable to determine q_dyn and q_OLOF for any HG donor. In our simulations, the mass-transfer rate and corresponding orbital evolution of the system never reach the thresholds for dynamical mass transfer, as described in Sect. 3.1. Similarly, our model HG stars never expand enough relative to their Roche lobes to overflow their OL in these systems.

Also visible in Fig. 7b is the sharp transition between the HG stars and giants, where the critical mass ratio abruptly increases by a factor of roughly 3−4. This coincides with the region where the structure of the outer layers of the donor stars changes significantly as a function of stellar radius. A significant part of the envelope of giant stars is already convective prior to the onset of mass transfer. As discussed in Sect. 4.1, stars with convective outer layers tend to expand in their initial response to mass loss (note that even for polytropic models this behaviour depends on the fractional core mass; see Hjellming & Webbink 1987) and hence transfer mass occurs at higher rates compared to stars with radiative envelopes. This results in significantly larger values for the critical mass ratio (i.e. less stable mass transfer) for the least evolved giants compared to the most evolved HG stars.

For giant stars, there is no visible discontinuity in q_qad values between donor stars filling their Roche lobes on the RGB, and those that do so on the AGB. Whilst there are structural differences between stars in both phases, these are mainly in the star’s core and burning shell(s). The outer layers are very similar in structure and deeply convective on both giant branches, and hence the trends in response to mass loss are very similar. As donor stars become more evolved giants, their core mass fraction increases (except for dredge-up events) and the local thermal timescales in the outer layers decrease (raising Ṁ_th,crit). Both have a stabilising effect on the mass transfer (see Sect. 3.3) and q_qad decreases with increasing donor radius. For those donor stars that fill their Roche lobes near the thermally pulsing part of the AGB, the local thermal timescales in the outer layers relevant for the response to mass transfer are short enough that q_qad can be even smaller than for HG donor stars.

With increasing donor radius on the RGB and AGB, the dynamical timescale decreases and the mass-transfer rates increase, while the donor stars extend further beyond their Roche lobes. Therefore, both q_dyn and q_OLOF typically follow a trend opposite to q_qad; these critical mass ratios generally increase with both log R and M_d. In the systems with giant donors for which we could not determine the critical mass ratios q_dyn and q_OLOF, the mass transfer turned unstable and our MESA simulations ended before the corresponding criteria could be reached, under conditions in accordance with our quasi-adiabatic criterion (see Sect. 3.3).

For our choice of A_dyn = 0.05, the dynamical and OLOF-based criteria result in similar critical mass ratios at the largest radii (log R/R_⊙ ≳ 1.5). At these initial donor radii, both of these criteria are typically met during the mass transfer in binaries where M_d ≥ 2.5 M_⊙. We stress again, as discussed in Sect. 4.1 that in this regime (where typically q_qad < q < {q_dyn, q_OLOF}) we do not find the mass transfer to run away and become unstable. Instead, the mass transfer is self-regulating and proceeds in a stable manner down to q = q_qad. However, in this part of the parameter space, multiple assumptions in our calculations (such as the Roche geometry and spherical symmetry) break down, and our results start to become unreliable. As such, the OLOF and dynamical criteria should be considered as indicators of where 1D-based calculations break down, rather than criteria for when mass transfer turns unstable.

Also shown in Fig. 7 is the effect of choosing different values for A_dyn (see Eq. (2)). The choice of A_dyn is somewhat arbitrary, as, to the best of our knowledge, there does not exist an obviously most suitable value for it. Adopting a larger value of A_dyn means that the evolution of binaries must become more rapid in order to be considered dynamical. This happens when the mass ratios are more extreme, such that larger values of A_dyn lead to lower critical mass ratios.

We summarise our results for the critical mass ratios for the full parameter space for all three stability criteria in Fig. 8. In general, the critical mass ratios depend on the mass and radius of the donor star when it fills its Roche lobe. However, the critical mass ratios have comparable values for donor stars in similar evolutionary stages. This holds especially well for the low-mass (M ≤ 2 M_⊙) stars in our simulations. As described in more detail above, the critical mass ratios vary with stellar radius and mass in systematic manners: q_qad is roughly constant for HG stars, increases by a factor of 3−4 at the transition to the giant phase, during which q_qad typically decreases with log R_d. The maximum value of q_qad generally decreases with M_d. The other critical mass ratios, q_dyn and q_OLOF follow the opposite trend and tend to increase with log R_d and M_d. For most of the parameter space studied in this work, we find that q_qad is smaller than 1, except for a small region where the maximum critical mass ratio we find is 1.05.

Fig. 8. Critical mass ratios (qcrit, where q ≡ Ma/Md) for all our binary models as interpolated contours. The three panels show critical mass ratios derived using different criteria (outlined in Sect. 3), as indicated at the top of each panel. The colour shows the value of the critical mass ratio, following the colour bar at the top of the figure. The results shown in panel b are for Adyn = 0.05. The green dashed lines show the donor radii corresponding to the end of the main sequence (TAMS), the base of the RGB (BRGB), and the tip of the RGB (TRGB), from bottom to top.

5. Discussion

In this section we compare our results to prior literature, and to commonly used critical mass ratios in BPS works. Furthermore, we critically examine sources of uncertainty in our results, and discuss several simplifying assumptions we made. We end by outlining the implications of our results for CE evolution, and a comparison of our results to observations.

5.1. Comparison to previous works

The classic work of Hjellming & Webbink (1987, hereafter HW87) provides systematic values for the critical mass ratio and adiabatic mass-radius exponent ζ_ad for a wide parameter space. The authors approximated the structure of giant stars as condensed polytropes that respond to mass loss adiabatically. By comparing the response of the stellar radii to the response of the Roche lobe, the authors derived critical mass ratios for dynamically unstable mass transfer as a function of the core mass fraction. We fitted their results (provided in Table 3 of HW87) by the following function, which is accurate to within 0.5% over the full range:

$$\begin{array}{c}\hfill q_{\mathrm{cr},\mathrm{HW}87}=1.577·\frac{{(1-m_{\mathrm{c}}^{3.533})}^{0.924}}{1+m_{\mathrm{c}}^{0.929}},\end{array}$$(11)

where m_c denotes the core mass fraction and q = M_a/M_d, as before. We assume the core mass fraction to be the same function of initial donor mass and radius as in our single-star models, and compare q_cr, HW87 to our results in Figs. 9a and b. We find that our critical mass ratios are significantly lower over the full parameter range studied in this work (i.e. we find mass transfer to be more stable than predicted by HW87), in qualitative agreement with earlier work. The difference increases with the radius at which the donors fill their Roche lobes. Our calculations do not assume that the donor star responds adiabatically and we find that our donor stars can find a state of self-regulating mass transfer down to significantly lower mass ratios. Whereas virtually all critical mass ratios of HW87 for the parameter space studied here are above unity, our results indicate that q_qad < 1 (see Sect. 4.2).

Fig. 9. Comparison of our results to the default sets of critical mass ratios (Ma/Md) implemented in the SeBa and binary_c rapid BPS codes and those found in detailed theoretical modelling efforts as a function of donor star radius. Each line corresponds to the source with the same style in the legend. As in Fig. 7, the solid black line and grey shading indicate the part of the parameter space that would lead to unstable mass transfer according to our quasi-adiabatic criterion.

Rapid BPS codes typically rely on simplified prescriptions for either the critical mass ratios or ζ_ad. The resulting boundaries for stable mass transfer are a collection of several literature values and fits, which are not necessarily consistent with one another in terms of the assumptions made in their derivation. We compare our results to the prescriptions for the critical mass ratios used in the binary_c code (see Izzard et al. 2004, 2018, and references therein) as given in Table 2 of Claeys et al. (2014). These critical mass ratios are based on Hurley et al. (2002), who use approximate fits to single stellar models to obtain ζ_eq (the response of a donor star’s equilibrium radius) for giants, and assume that ζ_ad = ζ_eq. This is the default prescription in binary_c, but we note that other options are available to users. Additionally, we compare our results to the prescriptions for the ζ_ad values used in the SeBa code (Portegies Zwart & Verbunt 1996; Toonen et al. 2012, as listed in Table A.1 of the latter). The SeBa code relies on the results from HW87 for giants with deeply convective envelopes, and employs a constant ζ_ad = 4 for HG stars, as well as giants with shallow convective envelopes. As before, we use the core mass fraction and donor mass as a function of donor radius from our own single-star simulations. We convert the aforementioned expressions for ζ_ad to critical mass ratios using the relations presented in Soberman et al. (1997), under the assumption of conservative mass transfer, appropriate for comparison to our results.

For HG donor stars, the critical mass ratios of binary_c are similar to our results, whilst those from SeBa are systematically larger, as can be seen in Figs. 9a and b. The transition from radiative to convective stars occurs at different radii between the two BPS codes and our models. In part, this is because SeBa employs a different prescription to determine beyond which radius the giants respond as convective stars during mass loss. For the giant stars, the critical mass ratios are lower in SeBa than in binary_c. This difference results from the assumption made in binary_c that ζ_ad = ζ_eq, whereas in reality ζ_ad > ζ_eq. Both sets of values, however, are significantly larger than ours, by a factor of ∼2−3. Whereas we find a strong decrease in the critical mass ratio along the giant branches, such a decrease with radius is much milder or absent in both BPS codes (except for the occasional jump due to dredge-ups). This suggests that predictions based on these popular rapid BPS models not only assume mass transfer to be less stable than we conclude, but also miss an important systematic effect in mass-transfer stability as a function of the donor radius at the onset of RLOF.

We next compare to results from various model calculations that do not implicitly assume an adiabatic response of the donor. By computing binary evolution models with the STARS/ev code (Eggleton 1971; Pols et al. 1995; Glebbeek et al. 2008), Chen & Han (2008, hereafter CH08) find the critical mass ratios for a similar parameter space as ours. Their results for giant stars exhibit a trend mostly opposite to our quasi-adiabatic critical mass ratios, more akin to our OLOF and dynamical criteria, and are typically larger than q_qad (see Figs. 9c and d). This might be due to the different way the mass-transfer rate is computed in their evolution code, which we discuss in more detail in Sect. 5.2.3. For their most evolved donor stars, the critical mass ratios found by CH08 exceed those from HW87 and those used in SeBa.

The same STARS/ev code, along with a ‘standard Henyey-type code’ (more details in Ivanova & Taam 2004) was used by Woods & Ivanova (2011, hereafter WI11) to perform a similar analysis. These authors focused on demonstrating the importance of the superadiabatic surface layer for the response of giant donors to mass loss and, as such, the parameter space they study is limited. The critical mass ratios they find follow a trend similar to those found by CH08, but offset to lower mass ratios. This offset could be explained by the difference in choice for the free parameter C in the expression for the mass-transfer rate (C = 1 in WI11 whilst C = 500 in CH08; see these works for more details), which according to the authors affects mass-transfer stability (but see Sect. 5.2.3). For low-mass stars, their critical mass ratios are significantly larger than the ones presented in this work, but they are more comparable for intermediate-mass stars, at least in the limited parameter space for which they were computed.

Pavlovskii & Ivanova (2015, hereafter PI15) studied the response of giant stars to mass loss in a similar manner as this work. These authors used a modified version of MESA for their simulations, with a set-up and physical parameters comparable to our simulations. The main differences between their simulations and ours are their inclusion of hydrodynamic terms in the stellar structure and evolution equations, and their modifications to the Kolb & Ritter (1990) mass-transfer scheme (see Sect. 5.2.3). We compare our results to the critical mass ratios of Pavlovskii & Ivanova (2015) as approximated by their ‘L₂/L₃-overflow condensed polytrope simplification’ results, which we extracted from their Fig. 15. Because their criterion for unstable mass transfer is based on L₂/L₃-overflow, we would expect it to be quantitatively similar to our OLOF criterion (see Sect. 3.2). Our values for q_OLOF follow a similar trend to the critical mass ratios of Pavlovskii & Ivanova (2015), but with a substantial offset to larger values by Δ(M_a/M_d)∼0.15 or more.

However, as can be observed in Fig. 14 of PI15, their L₂/L₃-overflow condensed polytrope simplification estimate of the critical mass ratios is not always consistent with their detailed calculations. In the case of a 1 M_⊙, 100 R_⊙ giant, our critical mass ratio for OLOF agrees very well with their detailed result, even though the difference between our result and their simplified approximation is significant. For this specific case their calculations for the mass-transfer rate and amount of overflow (as shown in their Fig. 12) also agree very well with ours. On the other hand, for less evolved 1 M_⊙ donor stars, their simplified approximation is closer to the results from their detailed calculations. Likewise, for a 2 M_⊙ donor star shown in their Fig. 14, the difference between their detailed calculations and their simplified approximation depends on the donor radius. Hence, the difference between our results and theirs cannot always be explained by this simplification. Lastly, these authors do not appear to encounter the limitation set by the critical thermal rate, although it is unclear why.

In a series of works focusing on the systematic study of mass-losing stars in binaries, which are similar in spirit to HW87 but do not rely on polytropic approximations, Ge et al. (2010, 2015, 2020) infer critical mass ratios from simulations of single stars. Under the assumption that the donor stars respond to mass loss adiabatically, the critical mass ratios are determined as the smallest mass ratio for which the mass-transfer rate does not exceed the global KH limit. They calculate the mass-transfer rate using a prescription very similar to to our Kolb scheme. The resulting critical mass ratios (q_ad in Ge et al. 2010, 2015, 2020) are very similar to our q_qad for HG stars, but different for giants. For giant donors, both sets of results follow a similar trend, but our critical mass ratios are systematically smaller by an almost constant difference Δ(M_a/M_d)∼0.2, meaning we typically find mass transfer to be more stable. Since we do not impose an adiabatic response on our donor stars, mass transfer can remain stable even when the mass-transfer rate greatly exceeds the global KH rate, as shown in Sect. 4.1. Instead, we find that the limiting mass-transfer rate is best taken to be one set by the shortest local thermal timescale in the envelope’s superadiabatic layers. To investigate whether this could be contributing to the difference between our results and those of Ge et al. (2010, 2015, 2020), we re-define our critical mass ratios in a similar manner to theirs; by setting the limiting mass-transfer rate to the global KH rate instead. We find that this leads to very similar critical mass ratios for low-mass stars (M_d ≤ 2 M_⊙), although our critical mass ratios deviate more from their results for larger initial donor masses.

In addition to their models with a realistic envelope structure, Ge et al. (2010, 2015, 2020) constructed adiabatic mass loss sequences involving stars with artificially isentropic envelopes, thus neglecting the superadiabatic layers. This is intended to mimic the effect of thermal relaxation of the surface layers (as shown in Fig. 2 of WI11 and our Fig. 2). The resulting critical mass ratios based on these isentropic envelope models (${\stackrel{\sim }{q}}_{\mathrm{ad}}$; the red dashed lines in Figs. 9c and d) are much more similar to our set of q_qad values than the Ge et al. (2020) values based on more realistic envelopes (q_ad, red solid lines in Figs. 9c and d). Apparently the effects of the two approximations made by Ge et al. (neglect of thermal relaxation in the outer layers, and neglect of the superadiabatic envelope structure) nearly compensate each other to give critical mass ratios very similar to our more realistic assumptions. We also note that, in Fig. 9d, the results of Ge et al. (2010, 2015, 2020; both q_ad and ${\stackrel{\sim }{q}}_{\mathrm{ad}}$) are shifted to somewhat larger radii compared to our q_qad values. This is probably due to differences in the physical assumptions made in the calculations regarding the amount of overshooting, the adopted mixing length and/or wind mass loss.

5.2. Limitations of the model

Here, we outline the most important simplifications we made in our calculations, and discuss how varying these might affect our results.

5.2.1. Treatment of the accretor and mass-transfer efficiency

We have ignored the evolution of the accretor in this work, treating it as an all-accreting point mass. However, at sufficiently high mass-transfer rates, accreting stars are unable to readjust their structure on timescales shorter than the mass transfer timescale, such that the accretor will be significantly brought out of thermal equilibrium. For those accreting main-sequence stars that have radiative envelopes (M_a ≳ 1.5 M_⊙), this causes substantial expansion of the accretor star beyond its main-sequence radius (Kippenhahn & Meyer-Hofmeister 1977; Neo et al. 1977; Packet & De Greve 1979). As a result, the secondary may fill its own Roche lobe, leading to the formation of a contact binary (Pols 1994; Wellstein et al. 2001). Although such a contact configuration does not necessarily lead to a CE, it could lead to OLOF from the accretor, and thus to non-conservative mass transfer (Marchant et al. 2016). Nonetheless, some stellar systems have been explained as a consequence of orbital instability after the accreting star expands into contact, including some apparently post-merger stars (e.g. Langer & Heger 1998; Justham et al. 2014), and close hot-subdwarf binaries (Justham et al. 2011).

Such a contact configuration may be avoided if stellar rotation is included. During mass transfer, the accreting star not only gains mass, but also angular momentum, which causes it to rotate faster. As the accretor spins up, the accretor’s surface rotation rate can eventually exceed the rotational break-up rate. It has been demonstrated (Packet 1981) that very little mass is required to spin up an accreting star to this break-up speed. Binary evolution models in which this effect is accounted for (Petrovic et al. 2005; de Mink et al. 2013) predict almost fully non-conservative mass transfer in massive binaries (however, also see Deschamps et al. 2013, for a different view). Additionally, rapid accretion could potentially produce strong winds and/or jets from an accreting main-sequence star (as seen in some binary post-AGB stars; see e.g. Van Winckel 2018).

Observational evidence supports a wide range of mass-transfer efficiencies, from almost conservative to almost fully non-conservative (see e.g. the references above and (De Greve & Linnell 1994; Nelson & Eggleton 2001; de Mink et al. 2007; Schootemeijer et al. 2018; Vinciguerra et al. 2020). Additionally, different forms of non-conservative mass transfer would have very different associated AM loss, which would affect the orbital evolution of the binaries. These remain as-of-yet unsolved problems in the evolution of interacting binary systems. Non-conservativeness during mass transfer, whether for the reasons just described or otherwise, affects the orbital evolution. In this work, we have assumed all of the mass transferred from the donor to be accreted by its companion. As mentioned in Sect. 2, under the often-made assumption that mass leaving the system carries with it the specific orbital angular momentum of the accreting star, the fully conservative case is expected to be the least stable one, based on orbital evolution arguments. As such, our results would provide an upper limit on the critical mass ratios separating stable and unstable mass transfer. We investigated the effects of non-conservative mass transfer by re-simulating a selection of systems representative of the full span of our parameter range with different (constant) mass-transfer efficiencies. We find that our completely conservative mass transfer calculations indeed provide a lower limit on the stability boundary.

5.2.2. Tidal interactions

Tidal interactions can act to affect the orbital evolution and therefore, potentially, the stability of mass transfer. For example, in close binaries, tidal interaction with the companion can transfer spin angular momentum back into the orbit, possibly preventing critical rotation. The effect of spin-orbit coupling on mass transfer was studied in detail by Misra et al. (2020), in the context of low- and intermediate-mass X-ray binaries containing neutron star accretors. They found that the inclusion of tidal interactions only had a minor effect on the stability of mass transfer. Since, in most of our systems, the main-sequence accretors are relatively compact compared to their post-MS companions, we conclude that the inclusion of spin-orbit coupling would not significantly alter our results or main conclusions. However, at sufficiently small mass ratios, tidal spin up of the donor star can lead to an ever growing orbital angular frequency, decreasing the orbital separation in a spiral-in manner. This phenomenon, known as the Darwin instability (see e.g. Darwin 1879; Counselman 1973; Hut 1980), occurs when the spin angular momentum of one of the binary components exceeds one third of the binary orbital angular momentum. Our calculations do not take this instability into account, but we estimate in which part of the parameter space it is expected to occur and where tidal effects may affect our results significantly. To that end, the donor radius at the moment it fills its Roche lobe should not exceed

$$\begin{array}{c}\hfill \frac{R_{\mathrm{d}}}{a}={\left[3r_{\mathrm{g}}^2(1+q^{-1})\right]}^{-\frac{1}{2}}\end{array}$$(12)

in order to avoid a Darwin instability (e.g. Darwin 1879; Eggleton & Kiseleva-Eggleton 2001; MacLeod et al. 2017). Here, r_g is the donor stars’ dimensionless radius of gyration ($r_{\text{g}}^2$ is defined as its moment of inertia $I=\frac{8\pi }{3}{\int }_0^{M_{\mathrm{d}}}\rho (r)r^4\mathrm{d}r$, divided by $M_{\text{d}}{R}_{\text{d}}^2$). We equate this expression to Eq. (1), and solve for q using our single-star models that form the initial configuration of our binary simulations. The resulting critical mass ratio for the Darwin instability is shown as a blue line in Fig. 7. As can be seen in that figure, the Darwin instability occurs only for mass ratios below our critical values for stable mass transfer, except perhaps in the most evolved intermediate-mass giants. Therefore, we do not expect the inclusion of tidal affects would alter our main results and conclusions.

5.2.3. The mass-transfer scheme

The comparison of our results to the literature in Sect. 5.1 suggests that the method used to calculate the rate of mass transfer can have strong effects on the stability of the RLOF, at least for evolved giant donors. In this work we used the prescription given in Kolb & Ritter (1990) as implemented in MESA (Paxton et al. 2015; the Kolb scheme). This method treats mass transfer as resulting from a laminar, subsonic flow along equipotential surfaces in the layers of the donor star that exceed the Roche lobe, towards the nozzle around the L₁ point. For the high mass-transfer rates that are relevant here, these layers are optically thick and the flow is assumed to be adiabatic with a constant adiabatic index. The vertical stratification of these layers is assumed to be hydrostatic and taken from the 1D stellar (MESA) model, so that the method explicitly takes into account the density structure in the overflowing layers.

By contrast, in some of the literature we compared with (in particular the studies of CH08 and WI11) the mass-transfer rate is taken to be a function of only the over-extension of the photospheric radius of the donor star beyond its Roche radius (i.e. Ṁ_d ∝ [(R_d − R_L)/R_L]³). Such a relation is expected if the outer layers have a polytropic structure with index n = 1.5 (Paczyński & Sienkiewicz 1972), which would be appropriate for an adiabatically stratified convective envelope. However, the proportionality factor is taken to be a constant free parameter (C) rather than calculated from the physical properties of the outer layers. This approach is sufficient for relatively unevolved donor stars with mass transfer on a thermal or nuclear timescale, where the resulting mass-transfer rate is insensitive to the precise method of computation (Kolb & Ritter 1990). On the other hand, for evolved red giants this simple approach can be problematic, for two reasons. Firstly, the outer envelopes of an extended convective giant are not adiabatically stratified but substantially superadiabatic. Secondly, the layers close to the photosphere have very low density and a pressure scale height that can be a substantial fraction of the stellar radius. This means that these stars must typically extend significantly past their Roche lobes for mass to be transferred at or above the KH rate. This is not taken into account in the simplified prescription used by for example, CH08 and WI11, which could therefore yield mass-transfer rates that are too high. This could explain why such studies result in larger critical mass ratios than our results, and why the difference grows with increasing radius (Figs. 9c and d). Conversely, the similar trend of critical mass ratio with stellar radius found by Ge et al. (2020) could result from their use of a method to compute the mass-transfer rate very similar to our Kolb scheme.

More recent work by Pavlovskii & Ivanova (2015) and Marchant et al. (2021) improved the Kolb & Ritter (1990) method for calculating the mass-transfer rate by taking into account the detailed shape of the Roche potential around L₁ and not approximating overflowing optically thick layers as an ideal gas, computing P/ρ from their detailed stellar models instead. These modifications are not expected to significantly change our results, as ideal-gas pressure always dominates in the outer layers of a red giant, and the effect of the improved L₁ geometry is small for mass ratios that are not too extreme (see Appendix A of Marchant et al. 2021). Indeed, for a 1 M_⊙, 100 R_⊙ giant donor we find almost the same mass-transfer rate and RLOF factor as Pavlovskii & Ivanova (2015, their Fig. 12).

Nevertheless, we should consider that other assumptions made in the Kolb & Ritter (1990) method may break down for evolved, very extended giants transferring mass at high rates. The assumed laminar flow may become turbulent and the structure of the overflowing layers may no longer be in hydrostatic equilibrium, even at locations far away from L₁. A 3D hydrodynamical treatment of the flow may be required to compute the mass-transfer rate in this case, and may lead to a different relation between Ṁ and the degree of Roche-lobe overfilling, potentially affecting the resulting critical mass ratios for all three criteria we consider in this paper.

5.2.4. Convection

Lastly, an appropriate description of convection is essential to determine the structure of a red giant. A reliable treatment of superadiabatic convection would require 3D hydrodynamical and radiation transport calculations (e.g. Freytag & Salaris 1999; Trampedach et al. 2014; Magic et al. 2015; Goldberg et al. 2022). A widely adopted solution within current 1D stellar evolution codes is different variants of the mixing-length theory (MLT; Böhm-Vitense 1958), a simplified analytical formulation of the problem. Whilst effective and widely in use, the MLT might not always be an appropriate treatment, especially in the scenarios considered here. Even for single-star models, studies such as the aforementioned Goldberg et al. (2022) show that the more realistic near-surface structure of giants can differ from those calculated using the MLT framework. Moreover, MLT implicitly assumes that the convective region is in a steady state, or at least that the timescale of change is much longer than the convective turnover time. In our simulations, as the evolution of the donor approaches the dynamical timescale, the timescale of the evolution becomes similar to or shorter than the convective turnover time when subsonic convective velocities are assumed. In this regime the approximations in MLT might be so inappropriate as to diminish the validity of our results.

5.3. Implications for common envelope evolution

In this work we have studied mass transfer in low- and intermediate-mass binaries and calculated three different critical mass ratios for stable mass transfer, corresponding to the three criteria described in Sect. 3. For systems with q < q_qad, unstable, runaway mass transfer appears to be guaranteed and will almost inevitably lead to a CE. This is broadly consistent with observational constraints from post-CE binaries (see Sect. 5.4).

However, our models cannot always provide a definite prediction of the outcome of mass transfer in binaries with q > q_qad. This is partly because of the limitations discussed in Sect. 5.2. The outcome of mass transfer in systems with q_qad < q < {q_dyn, q_OLOF} is particularly uncertain. In our models we find self-regulated mass transfer in this regime, albeit at very high rates, and we were able to simulate the entire mass transfer for most of these systems. One has to keep in mind, however, that we did not simulate OLOF and its effects on the orbital evolution of the binaries in this work, and the MESA code (in part due to its 1D nature) is not designed to simulate the types of rapid binary interactions that occur in systems with q ≲ q_dyn. Additionally, and perhaps more importantly, the outcome of mass transfer in such cases will crucially depend on the reaction of the accretor, which we did not take into account. The very high mass-transfer rates suggest that very little transferred material can be directly accreted by the companion star, and the excess material may rapidly fill the Roche lobe of the companion. Whether this leads to the onset of a CE followed by a spiral-in or to some other, less dramatic form of non-conservative mass transfer cannot be predicted from our models.

Finally, systems that can be classified as stable by all three of our criteria have the best prospects of avoiding a CE altogether. This by itself would significantly enlarge the parameter space for stable mass transfer compared to what is currently assumed in many popular rapid BPS codes.

5.4. Comparison with observed systems

To check that our results are consistent with expectations from observations, we compare our theoretical predictions to the sample of 62 post-common-envelope binaries (PCEBs) constructed by Zorotovic et al. (2011). The pre-CE properties of the binaries were inferred by the authors as follows. For each binary in the sample, the mass of the white dwarf was assumed to be equal to the pre-CE core mass of the donor star, and the mass of the companion star was assumed to be unchanged by the CE event. The single-star evolution code from Hurley et al. (2000) and the envelope binding energy prescription from Zorotovic et al. (2010), where the energy used to unbind the envelope is assumed to be a fraction α of the change in orbital energy during the spiral-in, were employed to calculate the pre-CE mass and radius of the donor star. These were then used to calculate the mass ratios at the onset of the ultimately unstable mass transfer (q_obs), which are shown in Fig. 10 with purple crosses. We note that the assumption made by Zorotovic et al. (2011) is that all the systems in their catalogue indeed experienced a CE. For IK Peg, we are not convinced this is indeed the most plausible formation scenario, given its wide current period of 21.7 days.

Fig. 10. Comparison of our results to observationally derived mass ratios. The purple ‘x’ markers indicate mass ratios of systems at the onset of mass transfer, inferred from observations of post-CE binaries. Different shades of grey show mass ratio ranges leading to unstable mass transfer based on the different criteria studied in this work, as indicated in the legend, interpolated to the properties of each system. The coloured symbols correspond to interpolated critical mass ratios based on results from a variety of sources, as indicated in the legend using the same abbreviations for literature sources as in the text and as used in Fig. 9.

Also shown in Fig. 10 by dark grey shading are our theoretical predictions for the range of mass ratios that lead to a CE, based on our q_qad for the assumed pre-CE donor mass and radius. The typical inferred progenitor of these systems is an evolved giant, in some cases on the thermally pulsing phase on the asymptotic giant branch (TPAGB), which means that we cannot always compare their inferred properties to our predictions, because we do not model mass transfer that begins during the TPAGB phase of our donor stars. Of those observed PCEBs for which comparison to our predictions is possible, we find that 33 out of 41 systems (roughly 80%) have q_obs ≤ q_qad, that is, they are consistent with our predictions for when a CE is unavoidable. However, only for two of the eight systems with q_obs > q_qad (EG UMa and IN CMa) we find this difference to be significant. Because uncertainties on the inferred pre-RLOF mass ratios are not provided in Zorotovic et al. (2011), we estimated them in the following manner. We use the provided uncertainties on the observed white-dwarf mass and companion mass, as well as the range of values for the α parameter that Zorotovic et al. (2010) claim can reproduce all the systems in this sample: 0.2 ≤ α ≤ 0.3 (though we note that for each individual system the range of α values the authors allow is much larger). Then, varying these parameters over their respective ranges of uncertainty individually and independently, we solved Eqs. (3)–(6) in Zorotovic et al. (2010) to derive an estimate of the uncertainty on q_obs. For six of the eight systems with q_obs > q_qad the difference with our results is within this uncertainty. All eight systems furthermore have q_obs ≤ {q_OLOF, q_dyn}, meaning that these systems are all consistent with stability limits derived from those criteria. However, as we stress in Sect. 5.3, this by itself does not necessarily mean a CE is unavoidable based on our results. All things considered, our predictions for the stability of mass transfer agree reasonably well with the inferred properties of PCEBs.

Figure 10 also shows critical mass ratios from literature, similar to Fig. 9. As it can be seen, the inferred pre-CE mass ratios are in all cases smaller than the critical mass ratios derived in these previous studies. This is not surprising given that our q_qad values tend to be smaller than those from the literature.

6. Conclusions

We have performed a systematic theoretical investigation of the critical mass ratios for stable mass transfer in low- and intermediate-mass binaries. We used the MESA stellar evolution code to simulate the mass-transfer evolution of 1404 binaries, with post-MS donor stars ranging from 1 M_⊙ up to and including 8 M_⊙ and with a wide range of mass ratios and orbital periods.

We explored three different criteria for identifying unstable mass transfer. Two are well-known and commonly used criteria, based on dynamical-timescale evolution and OL overflow by the donor star. Additionally, we used a criterion based on the transition to an effectively adiabatic response by the donor star, which takes the thermal response of the outermost layers into account when the timescale of mass transfer is so short that the relevant near-surface layers can no longer radiate away significant entropy before being lost. We find that this third criterion most consistently corresponds to the boundary between the stable, self-regulated mass transfer and unstable, runaway mass transfer shown by our mass-transfer simulations. For many cases of mass transfer from red giant donors where either the evolution accelerates to a dynamical timescale or the donor star swells far enough beyond its Roche lobe to fill its OL, or both, we find that the mass transfer can still be self-regulating and stable in our calculations. Rather than robust criteria for unstable mass transfer, these latter two should instead be employed to estimate where 1D stellar evolution codes might become unreliable in this context.

Therefore, our results should be interpreted with care. It is safe to assume that systems that can be classified as stable by all three criteria we employed will avoid a CE due to the response of the donor star. Similarly, binaries with q < q_qad can surely be expected to encounter unstable mass transfer, leading to a CE. However, for binaries with q_qad < q < {q_dyn, q_OLOF}, the outcome is still uncertain, owing to the limitations of our simulations. These include the fact that we did not evolve the accretors and only considered conservative mass transfer using the 1D MESA code.

In summary, our results confirm those of previous studies: typical assumptions of polytropic donor star structure and an adiabatic response of the donor star, on which prescriptions used in many rapid BPS models are based, severely underestimate the stability of mass transfer from giant donors. Additionally, compared to studies that assume the donor always responds to the loss of mass fully adiabatically, we find mass transfer to be more stable. Because we do not impose an adiabatic response of the donor star, the sub-surface layers of giants with convective envelopes are allowed to respond to mass loss on their short thermal timescales. This postpones a fully adiabatic response and allows stable mass transfer up to a critical mass-transfer rate, which is several orders of magnitude higher than the global KH rate.

We compared our theoretical limits for stable mass transfer to the Zorotovic et al. (2011) sample of PCEBs with inferred pre-RLOF properties, and we find them to be in good agreement. Our results have important consequences for interacting binaries. Typically, the critical mass ratio we find is significantly lower than currently assumed in population calculations, meaning an increase in the number of systems predicted to avoid a CE can be expected. Additionally, we find that the critical mass ratio decreases with donor radius along the giant branch much more strongly than in the typical assumptions made in rapid BPS codes. This would change the properties of post-mass-transfer populations in a rather complicated manner that is important to explore in more detail in the future.

Acknowledgments

The authors express gratitude to the anonymous referee, whose comments helped improve the quality of this work. K.D.T. furthermore sincerely thanks Jakub Klencki for many informative discussions and valuable advice regarding proper MESA usage and binary stellar simulations in general, and Evan Bauer for much-appreciated help understanding MESA logs and fixing errors. Lastly, K.D.T. thanks Rob Farmer, Glenn Oomen, Ylva Götberg and Mathieu Renzo for advising on more specific issues. Parsing of the data from MESA simulations was done using Rob Farmer’s mesadata module. K.D.T. acknowledges support from NOVA. S.J. acknowledges funding from the Netherlands Organisation for Scientific Research (NWO), as part of the Vidi research program BinWaves (project number 639.042.728, PI: de Mink). A.G.I. acknowledges support from the Netherlands Organisation for Scientific Research (NWO). This research made use of the following Python (Python Software Foundation, https://www.python.org/">https://www.python.org/) software packages and tools: matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), IPython (Pérez & Granger 2007), NumPy (Van Der Walt et al. 2011) and Pandas (McKinney 2010). The software references above were gathered from the Astronomy Acknowledgement Generator (http://astrofrog.github.io/acknowledgment-generator/).

References

Appendix A: Convergence study

In order to make sure that our results are numerically converged, we re-ran a selection of simulations of both stable and unstable mass transfer, varying the spatial resolution of our calculations as well as the constraints on the maximum time step size systematically.

More specifically, to assess dependence of our results on spatial resolution, we increase the value of the log_tau_function_weight parameter from 10 to 100 (we use a value of 50 in our grid calculations). This is a way to assign ‘meshing weight’ to regions in the star based on the logarithm of the optical depth. Qualitatively, this leads to smaller allowed jumps in mesh functions and hence higher resolution in the outermost layers of the star. Whilst increasing the number of mesh-points, especially in the outermost layers, generally increases the numerical stability and accuracy of the calculations, increasing them beyond our default leads to impractically long computation times without changing our conclusions about the stability of the mass transfer or results regarding its detailed evolution.

To test the dependence of our results on the size of the time-steps, we vary the fr (our default is 0.01) and delta_lg_star_mass_limit (our default is 10⁻⁴) parameters in the ranges 0.001-0.1 and 5 ⋅ 10⁻⁵ − 5 ⋅ 10⁻³, respectively, and adjust their associated ‘limit’ parameters accordingly. These parameters control the maximum allowed change in the donor star’s relative overflow past its Roche lobe and donor mass per time-step. We find that varying these settings has minimal effects on the simulated binaries’ evolution, besides affecting the calculation time. Another related setting we varied is the fr_dt_limit option. This parameter is the lower limit (in years) to the timestep size based on the donor star’s relative overflow past its Roche lobe, and thus limits the effectiveness of the parameters mentioned before. We varied this parameter between 100 and 0.01 (the MESA default is 10). We found that decreasing fr_dt_limit below its default could lead to increased numerical convergence and stability, reflected mostly in the number of time steps required to simulate challenging parts of the mass transfer as well as the number of retries and solver iterations. However, this behaviour was not monotonic and highly case-dependent. Typically, a value of fr_dt_limit = 1 provided the best balance between numerical stability and practical running times. Hence, this is our default value.

Furthermore, we experimented with two different expressions of the energy equation in the stellar structure and evolution equations (see Paxton et al. 2018, Sect. 8 for all built-in options in MESA). We found that there is no discernible difference in the evolution of the mass-transfer rate Ṁ_d(t) between the two (use_dedt_form_of_energy_eqn equal to .true. or .false.) for binaries that have mass ratios well above or below q_qad. Hence, for these systems, our conclusions are not affected by the choice of energy equations. For systems nearer the boundary between stable and unstable mass transfer, we often found it necessary to ‘hand-pick’ the most appropriate expression of the energy equation to ensure numerical stability and physical accuracy. In our grid of calculations, it did not appear that one of the two expressions we considered was consistently more appropriate than the other, and many times variation or use of other settings was necessary to ensure convergence.

Appendix B: Supplementary tables

Here we provide detailed tabulations of our main results. Table B.1 lists for all our binary models (except for the bisection models) key aspects of the mass-transfer-rate evolution, which can prove useful as a benchmark for rapid BPS models. Table B.2 enumerates for each of our grid points (the open circles in Fig. 1) the critical mass ratios based on the three criteria described in Sect. 3 as well as the critical thermal mass-transfer rate Ṁ_th,crit at the start of mass transfer. Specific values of and uncertainties on q_qad were obtained by performing bisecting MESA calculations until the uncertainty decreased below 5 per cent. Grid points with larger relative uncertainty correspond to models where the donor star expanded past the size of the binary orbit, which were not bisected farther.

All Tables

All Figures

Fig. 1. Single-star evolutionary tracks in the HRD for different stellar masses. Each colour corresponds to a given initial mass, as indicated to the lower left of each line. The coloured open circles indicate where we save single-star models for use in our binary star simulations. We have indicated lines of constant radius using grey dashed lines. Note that we have truncated these tracks to start at the ZAMS, and end at either the first TP or the maximum radius of the star, depending on which occurs first.

In the text

Fig. 2. Schematic illustration of the entropy (S) profile versus mass coordinate m in various types of stars. Panel a depicts a representation of a star with a radiative envelope, whilst panel b shows a model for a giant with a convective, isentropic envelope. A more realistic representation of the envelope of a convective giant is shown in panel c, which includes the superadiabatic sub-surface layers (highlighted in green). In each panel, the black line represents the entropy profile of the undisturbed star. Dashed red lines represent the entropy profile after sudden adiabatic (see text) loss of ΔM of mass. The dashed blue line in panel c shows the entropy profile if the mass is lost at a rate below the critical mass-loss rate (Eq. (9)).

In the text

Fig. 3. Thermal properties of a 1.5 M⊙ donor star at various stages in its post-MS evolution. The colours of the lines correspond to the moments the star reaches certain radii (as shown in Fig. 1), whose values are shown in the legend. Regions with grey shading are significantly superadiabatic (‘sad’; here defined as regions where ∇ − ∇ad ≥ 0.1∇; with $\mathrm{\nabla }=\frac{\partial logT}{\partial logP}$). The top panel shows the run of the local thermal timescale (see Eq. (7)) with the logarithm of the fractional exterior mass, such that the stellar centre is located on the left end of the x axis, and the surface is located to the right. The bottom panel shows the run of the associated thermal mass-loss rate (see Eq. (8)) on the same horizontal axis.

In the text

Fig. 4. Critical mass-loss rate as a function of stellar radius from thermal considerations (see Eq. (9)). Solid lines show the critical mass-loss rate based on local thermal properties of the sub-surface layers of the donor star (see Eq. (8) and Fig. 3), whilst dashed lines show the critical mass-loss rate based on global thermal properties (see Eq. (6)). The linear portions of the curves correspond to the giant phases. The colours of the lines correspond to different initial donor masses, as indicated to the lower left of each dashed line.

In the text

Fig. 5. Time evolution of mass-transfer rates in systems with Md = 2 M⊙ (left column) and Md = 5 M⊙ (right column). Each row corresponds to RLOF starting in a different phase of the donor stars’ lives (at the radius indicated in the top right of each panel). More specifically, panels a and b show HG donors, panels c and d show donors on their RGB, and panels e and f show donors on their AGB. The colours of the lines correspond to different mass ratios at the onset of RLOF, as indicated in the legend. We have indicated where, if appropriate, our quasi-adiabatic-, dynamical-, and OLOF-based criteria are met with open circles, squares, and triangles, respectively. The zero point of the time axis is defined as the first moment when the mass-transfer rate exceeds the rate at which wind mass loss occurs.

In the text

Fig. 6. Mass-transfer rates from our calculations compared to the global KH rate (as given by Eq. (6)). The top panel shows results for a 2 M⊙ donor star, whilst the bottom panel shows results for a 5 M⊙ donor star. In each panel, the coloured lines show the difference between the mass-transfer rate predicted by Eq. (6), (which assumes the rate is governed by the global KH timescale at the onset of mass transfer) and the mass-averaged mass-loss rate ⟨Ṁd⟩M (see Eq. (10)) from our calculations. Vertical dashed black lines indicate the base and the tip of the RGB.

In the text

Fig. 7. Critical mass ratios (above which mass transfer is stable) for 1 M⊙ and 4 M⊙ donor stars as a function of donor radius, under different criteria for stable mass transfer (see Sect. 3). The solid black line shows the critical mass ratios resulting from our quasi-adiabatic criterion, with the grey shading indicating the parameter space where mass transfer would be unstable. The red lines show the critical mass ratios resulting from the dynamical evolution criterion with Adyn = 0.01, 0.05, 0.25, as indicated to the right end of the lines. The solid orange line shows the critical mass ratios resulting from the OLOF-based criterion, and the solid blue line and hatching show where the Darwin instability could be encountered (see Sect. 5.2). The solid horizontal grey line shows where Ma/Md = 1, and the dashed vertical grey lines show the end of the HG phase (panel a and left line in panel b) and the tip of the RGB (panel b).

In the text

Fig. 8. Critical mass ratios (qcrit, where q ≡ Ma/Md) for all our binary models as interpolated contours. The three panels show critical mass ratios derived using different criteria (outlined in Sect. 3), as indicated at the top of each panel. The colour shows the value of the critical mass ratio, following the colour bar at the top of the figure. The results shown in panel b are for Adyn = 0.05. The green dashed lines show the donor radii corresponding to the end of the main sequence (TAMS), the base of the RGB (BRGB), and the tip of the RGB (TRGB), from bottom to top.

In the text

Fig. 9. Comparison of our results to the default sets of critical mass ratios (Ma/Md) implemented in the SeBa and binary_c rapid BPS codes and those found in detailed theoretical modelling efforts as a function of donor star radius. Each line corresponds to the source with the same style in the legend. As in Fig. 7, the solid black line and grey shading indicate the part of the parameter space that would lead to unstable mass transfer according to our quasi-adiabatic criterion.

In the text

Fig. 10. Comparison of our results to observationally derived mass ratios. The purple ‘x’ markers indicate mass ratios of systems at the onset of mass transfer, inferred from observations of post-CE binaries. Different shades of grey show mass ratio ranges leading to unstable mass transfer based on the different criteria studied in this work, as indicated in the legend, interpolated to the properties of each system. The coloured symbols correspond to interpolated critical mass ratios based on results from a variety of sources, as indicated in the legend using the same abbreviations for literature sources as in the text and as used in Fig. 9.

In the text