The Impact of Angular Momentum Loss on the Outcomes of Binary Mass Transfer

We review the analytics for binary MT and define the relevant free parameters. We borrow heavily from the didactic approach of Pols (2018) as well as Soberman et al. (1997). For the simplified case of a circular, corotating binary, the stellar material around each star is gravitationally bound to its host if the stellar radius does not exceed its Roche lobe, the equipotential surface passing through the L1 point. If the star expands beyond its Roche lobe, gas will flow through a nozzle near the L1 point and fall toward the companion star. The volume-equivalent Roche-lobe radius R_L of the donor is

$$\begin{eqnarray}&&{R}_{L}=\displaystyle \frac{0.49a}{0.6+{q}^{2/3}\mathrm{ln}(1+{q}^{-1/3})},\end{eqnarray} \tag{ 1 }$$

where a is the orbital separation and q = M_a /M_d is the binary mass ratio in terms of the mass of the accretor M_a and the donor M_d (Eggleton 1983). Note that q here is inverted from the definition used in Eggleton (1983) and Team COMPAS Team COMPAS et al. (2022). If the donor radius R_d > R_L , MT ensues via Roche-lobe overflow (RLOF).

For a binary in a circular orbit, the orbital AM is

$$\begin{eqnarray}&&J={M}_{a}{M}_{d}\sqrt{\displaystyle \frac{{Ga}}{{M}_{a}+{M}_{d}}},\end{eqnarray} \tag{ 2 }$$

where G is the gravitational constant. If the orbit is eccentric, MT may begin in bursts near the point of closest approach; we assume that MT drives the binary to circularize at periapsis (though see Sepinsky et al. 2007, 2009, 2010; Dosopoulou & Kalogera 2018).

As MT proceeds, some fraction β of the mass lost from the donor will be accreted by the companion

$$\begin{eqnarray}&&\dot{{M}_{a}}=-\beta \dot{{M}_{d}}.\end{eqnarray} \tag{ 3 }$$

The MT efficiency β in general depends on the properties of the binary prior to the MT and changes with time, though it is often assumed to be a constant for a given MT episode. β = 1 corresponds to conservative MT, in which all material lost from the donor is accreted, while β < 1 refers to nonconservative MT, and β = 0 is fully nonconservative.

Any mass that is not accreted is ejected from the binary system, extracting orbital AM. Traditionally, this is represented by the specific AM γ of the ejected material in units of the orbital specific AM. For simplicity, we introduce a new parameterization for the AM loss, the effective decoupling radius a_γ . This can be thought of as the distance from the center of mass at which matter is ejected. The parameters a_γ and γ are related by

$$\begin{eqnarray}\begin{array}{rcl}\gamma & = & \left(\displaystyle \frac{\dot{J}}{\dot{{M}_{a}}+\dot{{M}_{d}}}\right)/\left(\displaystyle \frac{J}{{M}_{a}+{M}_{d}}\right)\\ & = & \left({a}_{\gamma }^{2}{\omega }_{\mathrm{orb}}\right)/\left(\displaystyle \frac{J}{{M}_{a}+{M}_{d}}\right)\\ & = & {\left(\displaystyle \frac{{a}_{\gamma }}{a}\right)}^{2}\displaystyle \frac{{\left(1+q\right)}^{2}}{q},\end{array}\end{eqnarray} \tag{ 4 }$$

where ω_orb is the orbital angular frequency.

As with the efficiency parameter β, γ is generally dependent on time and the attributes of the binary, such as the mass ratio and degree of rotational synchronization (MacLeod & Loeb 2020a, 2020b), but is similarly often simplified to a constant. For matter which is removed from the vicinity of the donor star (the "Jeans" mode for a fast, isotropic wind), a_γ = aq/(1 + q), so γ = q. Material which instead decouples from the binary near the accretor has a_γ = a/(1 + q) and γ = 1/q (the "isotropic reemission" mode). The maximal specific AM loss for corotating matter occurs at the L2 point, a_γ ≈ a_L2 ≲ 1.25a (see Figure 1). ⁵ Larger specific AM loss may be achieved if the ejected matter forms a circumbinary ring which applies a torque on the orbit, though we do not consider this case here.

Zoom In Zoom Out Reset image size

Taking the time derivative of Equation (2) and rearranging with the definitions of β and γ, the relative change in separation is

$$\begin{eqnarray}&&\displaystyle \frac{\dot{a}}{a}=-2\displaystyle \frac{\dot{{M}_{d}}}{{M}_{d}}\left[1-\displaystyle \frac{\beta }{q}+(\beta -1)\left(\gamma +\displaystyle \frac{1}{2}\right)\displaystyle \frac{1}{1+q}\right].\end{eqnarray} \tag{ 5 }$$

2.1.1. Stability Criteria

To predict the outcome of the binary postinteraction, we distinguish MT which leads to CEE from that which does not. Using Equations (1) and (5), we define the response of the donor's Roche lobe to mass loss by the logarithmic derivative

$$\begin{eqnarray}\begin{array}{rcl}{\zeta }_{L} & = & \displaystyle \frac{{\rm{\partial }}{\rm{ln}}{R}_{L}}{{\rm{\partial }}{\rm{ln}}{M}_{d}}\\ & = & \displaystyle \frac{{\rm{\partial }}{\rm{ln}}a}{{\rm{\partial }}{\rm{ln}}{M}_{d}}+\displaystyle \frac{{\rm{\partial }}{\rm{ln}}({R}_{L}/a)}{{\rm{\partial }}{\rm{ln}}q}\displaystyle \frac{{\rm{\partial }}{\rm{ln}}q}{{\rm{\partial }}{\rm{ln}}{M}_{d}}\\ & = & -2\left[1-\displaystyle \frac{\beta }{q}+(\beta -1)\left(\gamma +\displaystyle \frac{1}{2}\right)\displaystyle \frac{1}{1+q}\right]\\ & & +\displaystyle \frac{1}{3}\left(2-\displaystyle \frac{1.2+{\left({q}^{-1/3}+{q}^{-2/3}\right)}^{-1}}{0.6+{q}^{2/3}{\rm{ln}}(1+{q}^{-1/3})}\right)\left(\displaystyle \frac{\beta }{q}+1\right).\end{array}\end{eqnarray} \tag{ 6 }$$

This is purely a function of q, β, and γ, and is explicitly independent of the MT timescale. For most of the parameter space, MT from a high-mass donor to a low-mass accretor will tighten the binary and shrink the Roche lobe, and vice versa. Equation (6) is shown graphically in Figure 1 as a function of these parameters, with γ replaced by the decoupling separation a_γ /a according to Equation (4).

By contrast, the radial response of the donor to mass loss

$$\begin{eqnarray}&&{\zeta }_{* }=\displaystyle \frac{\partial \mathrm{ln}{R}_{d}}{\partial \mathrm{ln}{M}_{d}},\end{eqnarray} \tag{ 7 }$$

strongly depends on the MT timescale, as well as the evolutionary phase of the donor. Following a small amount of mass loss, the stellar structure will react on the short dynamical timescale to restore hydrostatic equilibrium. This initial, adiabatic radial response is parameterized as ζ_* = ζ_ad.

Following Soberman et al. (1997), we label MT as dynamically unstable if the adiabatic response of the donor to mass loss at the beginning of MT is to expand relative to its Roche lobe, ζ_* < ζ_L , leading to an increase in the mass-loss rate and a runaway process on the dynamical timescale. According to this definition, unstable MT leading to CEE sets in only if the MT occurs on the dynamical timescale at the onset of the mass exchange. Otherwise, the MT is stable.

Equation (7) can be solved analytically for some simplified polytropic models (Hjellming & Webbink 1987; Soberman et al. 1997), or numerically using detailed 1D models of stars undergoing mass loss (de Mink et al. 2007; Ge et al. 2020; Klencki et al. 2021). Broadly, the adiabatic response of the star is sensitive to the specific entropy gradient at the surface. For stars with radiative envelopes (intermediate- and high-mass stars on the main sequence (MS) and early on the Hertzsprung gap (HG)), the specific entropy rises steeply near the surface. Mass loss exposes lower-entropy layers which contract rapidly (ζ_ad > 0). As the mass loss continues, in the adiabatic regime (i.e., fixing the entropy profile) the entropy gradient tapers off and the donor contraction is moderated.

In stars with convective envelopes (low-mass MS stars and cool giants), energy is efficiently transported in convective layers leading to a flat specific entropy profile within the convective region. Unless the core mass constitutes ≳20% of the total stellar mass, the adiabatic response of the donor to mass loss is to expand (ζ_ad < 0; Soberman et al. 1997). In the extreme case of a fully convective star, the response of the donor follows R ∝ M^−1/3, i.e., ζ_ad = −1/3 (Hjellming & Webbink 1987; Kippenhahn et al. 2012).

However, Ge et al. (2010) argued that comparing ζ_* and ζ_L only at the onset of MT may be a poor predictor of MT stability, since both ζ_* and ζ_L change throughout the MT. They found that while fully convective stars do expand rapidly during adiabatic mass loss, the outermost layers which extend beyond the Roche lobe are very diffuse, so that the actual mass-loss rate is low. Although this is initially a runaway process, the mass-loss rate may never exceed a (model-dependent) critical value, in which case the MT will proceed stably.

Conversely, for some radiative donors, ζ_ad may be initially very high (rapid contraction), but decrease over time until ζ_ad ∼ ζ_L , at which point the binary will experience what is known as the delayed dynamical instability (Hjellming & Webbink 1987; Ge et al. 2010). As with the convective case, this delayed effect cannot be effectively captured from the initial stellar response at the onset of MT, as defined in Soberman et al. (1997), but requires a detailed treatment of the stellar response over the MT episode.

2.1.2. Critical Mass Ratios

We use the COMPAS ⁶ rapid BPS suite (Stevenson et al. 2017; Vigna-Gómez et al. 2018), including the fiducial parameter choices listed in Team COMPAS Team COMPAS et al. (2022), with a few exceptions and recently added prescriptions detailed below. For each distinct model, we evolve 10⁵ binaries from zero-age main sequence (ZAMS). By default, we follow the traditional initial sampling distributions used in many BPS codes. For the more massive primary, we draw its mass M₁ from the Kroupa initial mass function (IMF) $p({M}_{1})\propto {M}_{1}^{-2.3}$ between 5 and 100 M_⊙ (Kroupa 2001). The mass ratio distribution is uniform q ∈ [0.1, 1] (with no additional constraint imposed on the minimum of M₂), and the separation is drawn log-uniformly, $\mathrm{log}(a/\mathrm{AU})\in [-2,2]$, which has a lower upper limit than the COMPAS default 3 (see Team COMPAS Team COMPAS et al. 2022). All single stars are included implicitly as very wide binaries for normalization purposes. Since wide binaries do not interact, extending this upper limit is analogous to increasing the effectively single star fraction. As our focus is on interacting binaries, we omit this noninteracting subset. The initial eccentricity is fixed at 0. We discuss variations to the initial parameter distributions in Appendix.

Each binary is then evolved under the specified evolution model until one of the following termination conditions is reached:

• 1.  
  the binary experiences unstable MT,
• 2.  
  either star experiences a supernova (SN), or
• 3.  
  both components evolve into white dwarfs.

We investigate the binaries at two epochs: immediately following the first MT episode and following the final MT episode prior to the termination condition, henceforth referred to as the "End" state. We classify binaries based on whether or not they have experienced only stable MT up to the specified epoch. We additionally consider the orbital separations following the first MT episode, if the MT was stable. The outcome of CEE is notoriously uncertain and the subject of ongoing research, so for simplicity we do not distinguish common-envelope survival from stellar mergers here (Ivanova 2017; Hirai & Mandel 2022; Lau et al. 2022).

By restricting our study to pre-SN and pre-double white dwarf interactions, and by ignoring post-CEE outcomes, we limit the scope of our parameter space and ensure that our conclusions are more robust. Many systems will not interact: this is largely dependent on the initial conditions (and, to a lesser extent, on the stellar winds at early phases). These noninteracting systems are generally ignored here, except where otherwise specified.

2.3.1. Roche-lobe Response

The rate of accretion onto a nondegenerate companion is often assumed to be restricted by the accretor's thermal timescale, $\dot{{M}_{a}}\sim {M}_{a}/{\tau }_{\mathrm{KH},a}$. We calculate the accretion efficiency as

$$\begin{eqnarray}&&\beta =-\displaystyle \frac{\dot{{M}_{a}}}{\dot{{M}_{d}}}=\min \left(\displaystyle \frac{{{CM}}_{a}/{\tau }_{\mathrm{KH},a}}{{M}_{d}/{\tau }_{\mathrm{KH},d}},1\right),\end{eqnarray} \tag{ 8 }$$

where τ_KH is the Kelvin–Helmholtz timescale, and the prefactor C (10 in this study) is included to account for an increase in the maximum accretion rate due to expansion of the accretor (Paczyński & Sienkiewicz 1972; Neo et al. 1977; Hurley et al. 2002; Schneider et al. 2015). During Case A MT—defined when the donor is on the MS—β ≈ 1 if the component masses are roughly comparable. During Case B/C MT—defined for more evolved donors—β is typically much closer to 0. Throughout, we follow the stellar type convention defined by Hurley et al. (2000).

We first consider the impact of uncertainties in the Roche-lobe response to mass loss, ζ_L , by varying the MT efficiency β and the specific AM of nonaccreted matter γ. The default efficiency in COMPAS, β_Comp, uses the ratio of the thermal timescales defined in Equation (8). We also include variations with fixed values of β = {0.0, 0.5, 1.0}. To investigate the effect of different AM loss modes, we introduce a new parameter, f_γ ∈ [0, 1], which increases linearly with the decoupling radius between the accretor a_acc and the L2 point a_L2

$$\begin{eqnarray}&&{a}_{\gamma }={a}_{\mathrm{acc}}+{f}_{\gamma }({a}_{L2}-{a}_{\mathrm{acc}}),\end{eqnarray} \tag{ 9 }$$

such that f_γ = 0 corresponds to the isotropic reemission model and f_γ = 1 corresponds to AM loss from the L2 point. This range in f_γ reflects results from hydrodynamic simulations, which indicate that nonaccreted mass could decouple from the binary between the accretor and the L2 point, depending on the structural properties of the donor (MacLeod et al. 2018b, 2018a; MacLeod & Loeb 2020b). Many rapid BPS codes assume the isotropic reemission model by default.

Since the impact of γ is moderated by β, we consider pairwise combinations of these parameters, noting that γ is irrelevant for fully conservative MT.

2.3.2. Donor Radial Response

We additionally consider variations to the radial response of the donor to mass loss ζ_ad. In our default model both MS and HG stars (in the mass range of interest) have radiative envelopes for the entirety of the phase, with ζ_ad,MS = 2 and ζ_ad,HG = 6.5, respectively (these values are motivated by a constant approximation to the Ge et al. 2020 results; see Team COMPAS Team COMPAS et al. 2022). Giants, by contrast, are assumed to have convective envelopes. By default, we follow the prescription outlined in Soberman et al. (1997). According to this prescription, mass loss leads to radial expansion unless the core mass fraction becomes nonnegligible. We refer to this set of prescriptions encompassing all stellar types as ζ_Comp.

We also include several prescriptions for stability based on the critical mass ratios from Claeys et al. (2014) and Ge et al. (2020). From Claeys et al. (2014) we obtain q_c,Claeys14, which takes on constant values that depend on the donor stellar type and whether or not the accretor is a degenerate star. For giant stars, they use a function of the core mass ratio (Claeys et al. 2014; Hurley et al. 2002). The q_c,Claeys14 prescription is defined based on the initial donor response at the onset of mass loss, and does not capture the delayed stability effects discussed in Section 2.1.1.

The prescriptions from Ge et al. (2020) use detailed 1D adiabatic mass-loss models to account for the evolution of stability over the course of the MT episode (see Section 2.1.1), and build off of Ge et al. (2010, 2015) to provide q_c as a function of both donor mass and radius (or analogously, evolutionary state). The parameter q_c,Ge20 refers to their critical mass ratios calculated using adiabatic mass loss from their standard stellar profiles (including superadiabatic regions where relevant), while the ${\tilde{q}}_{c,\mathrm{Ge}20}$ variant refers to the critical mass ratios calculated when the donor envelopes are artificially isentropic, as a work-around for the rapid superadiabatic expansion observed in their standard models (see Ge et al. 2020 for a more thorough explanation).

The ζ_ad prescriptions are shown in Figure 2 as a function of donor radius for a grid of masses. Note that ζ_Ge20 is calculated directly from the critical mass ratio q_c,Ge20, as ζ_Ge20 = ζ_L (q_c,Ge20; β = 1), and similarly for ${\tilde{\zeta }}_{\mathrm{Ge}20}$ and ζ_Claeys14. We include them in Figure 2 to highlight differences with ζ_Comp. In practice, the Ge et al. (2020) values are interpolated in both mass and radius. The vertical lines represent the stellar radii in our models when the stellar type transitions to the named type.

Zoom In Zoom Out Reset image size

For a given binary undergoing RLOF, the MT efficiency β and the specific AM of the nonaccreted material γ completely determine the response of the Roche lobe, as well as the final separation of the binary following stable MT, via Equation (5).

In Figure 3, we show how the post-MT orbital separations of binaries that experienced stable MT depend on the specific AM loss through our proxy parameter f_γ . Here, we follow default assumptions in the other parameters, namely, ζ_* = ζ_Comp and β = β_Comp.

Zoom In Zoom Out Reset image size

Figure 3 shows the separation distribution before (pre-MT, top half) and after (post-MT, bottom half) the first episode of MT. Figure 3 shows the orbital tightening ratio $\mathrm{log}({a}_{\mathrm{post}}/{a}_{\mathrm{pre}})$ for the same systems. Different colors represent different values for f_γ , as indicated in the legend. Only binaries that undergo stable MT are shown, so variations in the normalization for different values of f_γ correspond to differences in the number of binaries which remain stable.

As expected, increasing f_γ leads to a reduction in the number of systems that undergo stable MT, as many of these are driven toward instability. This is particularly prominent for the pre-MT distributions (upper half of Figure 3), in the range $\mathrm{log}(a/{R}_{\odot })\in \sim [2.7,3.5]$, which show a population of systems that undergo stable MT only for very small values of f_γ . The post-MT separation distributions cluster into two peaks near $\mathrm{log}(a/{R}_{\odot })\sim$ 1.5 and 2.5. For f_γ = 0.0, there is also a prominent bump toward large post-MT separations which is not seen in the pre-MT distribution, indicating that these are systems which ultimately widened as a result of the MT.

This can be seen more directly in Figure 3, which shows the fractional amount of tightening or widening that binaries experience. For most values of f_γ shown, the separation ratio distribution is bimodal, with peaks at ${a}_{\mathrm{post}}/{a}_{\mathrm{pre}}\sim 1$ (little net change in the separation) and ${a}_{\mathrm{post}}/{a}_{\mathrm{pre}}\sim 5$ (substantial widening). Interestingly, the peak at ${a}_{\mathrm{post}}/{a}_{\mathrm{pre}}\sim 5$ becomes relatively less prominent with decreasing f_γ , except for f_γ = 0 when that peak dominates over the one at ${a}_{\mathrm{post}}/{a}_{\mathrm{pre}}\sim 1$.

To understand this behavior, we highlight that, for fully conservative MT, the orbital separation decreases until the mass ratio reverses. While this is only approximately true for nonconservative MT, it is nonetheless illustrative. For nearly equal mass ratio binaries undergoing Case A MT, β_Comp ≈ 1, so the mass ratio will reverse after only a modest amount of mass is lost, and the MT will cease when the Roche lobe encompasses a donor radius only moderately smaller than it was at the onset of MT. These binaries thus contribute primarily to the peak at ${a}_{\mathrm{post}}/{a}_{\mathrm{pre}}\sim 1$, independently of the AM loss prescription.

Case A MT in unequal mass ratio binaries will initially be nonconservative, leading to a dependence on the AM loss. However, as the donor mass approaches that of the accretor, the MT efficiency β will rise toward unity. If the orbit remains sufficiently wide following the initial phase of nonconservative MT, the mass ratio will reverse and the the orbit will widen from its minimum value. These binaries thus contribute primarily to the extended tail toward low ${a}_{\mathrm{post}}/{a}_{\mathrm{pre}}$ values. If instead the orbit is too tight following this phase, due perhaps to a high f_γ value, the two stars will merge.

Case B systems, by contrast, are nearly always nonconservative, β_Comp → 0. If the binary has a mass ratio close to unity, the orbit will not shrink much before the mass ratio reversal, and is unlikely to become unstable. The mass loss will then act to widen the binary substantially, populating the second peak at ${a}_{\mathrm{post}}/{a}_{\mathrm{pre}}\sim 5$.

The impact of f_γ variations is maximal for low values of the MT efficiency β and reduced as β → 1. We investigate the influence of β and f_γ by varying both parameters simultaneously in Figure 4. In Figure 3, we show how the stable MT fraction during the first MT episode depends on the adopted value of f_γ (shown on the abscissa), assuming the default donor response ζ_* = ζ_Comp. Higher β values (represented by line color) decrease the sensitivity of this f_γ dependence. Figure 3 shows the same analysis for systems evaluated prior to the first SN, which we discuss in more detail in Section 3.4.

Zoom In Zoom Out Reset image size

The height of each line represents the stable MT fraction f_SMT. The remaining fraction, 1 − f_SMT, includes systems that experienced unstable MT during the first interaction, including any stellar mergers. Noninteracting systems are excluded. When β = 0, f_SMT varies over a very broad range from ∼0% – 50%. Even for more realistic values of β = β_Comp, the range is still substantial, closer to ∼20%–50%.

In Figure 5, we show f_SMT across all model variations following the first MT episode (solid, red line) and at the End state (dashed, green line), which we discuss further in Section 3.4. Model variations are listed on the abscissa. The ordering groups together similar models; those on the left have fixed ζ = ζ_Comp and β and γ vary, while those on the right have a varied donor response prescription and β is implicitly fixed to 1. For models which depend on β and f_γ , we include only the extremal values of these parameters (as well as our default β_Comp) for brevity.

Zoom In Zoom Out Reset image size

All of the models which define stability based on critical mass ratios q_c assume fully conservative MT (β = 1). The ζ_Comp model furthest to the right also includes fully conservative MT to facilitate comparison with the critical mass ratio models. The (ζ_Comp, β = 1) model and q_c,Claeys14 both define stability at the onset of MT, i.e., ignoring delayed instabilities.

By contrast, the Ge et al. (2020) models consider the evolution of stability over the course of MT. The aggregate effect is to increase the fraction f_SMT of systems undergoing stable MT during the first MT event from ∼40% for the (ζ_Comp, β = 1) model to ∼70% for the Ge et al. (2020) models. This is primarily due to the reduction in convective envelope donors which experience instability in Ge et al. (2020), since the delayed dynamical instability mechanism in radiative donors acts in the other direction, to increase the number of unstable systems that would previously have been considered stable.

Meanwhile, we see as before that varying just β and γ at fixed ζ_* = ζ_Comp drives significant changes in f_SMT. The maximal variation in f_SMT from all variations in β and γ is ∼50%—seen in the difference between (ζ_Comp, β = 0.0, f_γ = 0.0) and (ζ_Comp, β = 0.0, f_γ = 1.0). For the default COMPAS efficiency β = β_Comp, f_SMT varies by ∼25% across the full range of f_γ . Fluctuations in the stable MT fraction f_SMT due to changes in the efficiency and specific AM loss during MT are comparable in magnitude to the variations from the donor stability prescriptions.

In Figure 3, we consider the impacts of the β−γ variations at the End State. We find that ∼0%–20% of interacting binaries will undergo only stable MT throughout their evolution.

As with the first MT episode, the stable MT fraction at the End state is strongly dependent on the adopted AM loss prescription f_γ (for values of β < 1), but here the impact is more concentrated at low f_γ values. For f_γ ≳ 0.7, f_SMT drops to a few percent with only a weak dependence on f_γ , which suggests that if f_γ is universally high, even binaries which remain stable during the first MT will eventually experience unstable MT.

We show f_SMT at the End state across all model variations in Figure 5 (dashed, green line). Here, the differences between models are qualitatively similar to the differences following the first MT episode, though not quite as pronounced. Interestingly, although f_SMT is still sensitive to β−γ variations at the end of the evolution, the impact of different donor responses (for fully conservative MT, β = 1) is greater at this epoch, with f_SMT ranging from ∼15% to 45%. Taken at face value, this suggests that uncertainties in the donor response are more influential for the long-term evolution than uncertainties in the Roche-lobe response, though we emphasize that it is difficult to estimate without models that vary both the donor response and the Roche-lobe response simultaneously.

The fraction of systems that undergo stable MT during the first MT event, as well as the fraction of systems that will experience only stable MT throughout their evolution, strongly depends on the specific AM γ carried away by nonaccreted material (or equivalently, our proxy parameter f_γ ), as shown in Figure 5. BPS codes often assume that nonaccreted matter follows the "isotropic reemission" model, taking with it the the specific AM of the accretor, i.e., f_γ = 0, in all cases. However, both models and observations suggest that this may not be the case.

Using detailed hydrodynamic simulations, MacLeod et al. (2018a, 2018b) find that f_γ should not be constant, but depends sensitively on properties of the system, such as the evolutionary state of the donor. Kalogera & Webbink (1996) argue that f_γ ∼ 0 if the accretor is an NS or BH on the basis that gas ejection is driven by accretion energy deposited very close to the compact object. However, Goodwin & Woods (2020) find that f_γ < 0 (mass lost near the L1 point) is required to explain the large period derivative of SAX J1808.4-3658. Moreover, there may exist mechanisms to keep ejected matter synchronized to the binary orbit beyond the L2 point, such as magnetic fields, collisions in the ejecta streams, or continuous torquing from circumbinary disks. Thus f_γ may not be a priori limited to f_γ ≤ 1.

We find that the fraction of systems which undergo stable MT during the first MT event is sensitive to the value of f_γ across the entire range f_γ ∈ [0, 1], but that by the end of the evolution, most systems experience CEE at some point if f_γ is universally high. Under the default β_Comp efficiency, fewer than 15% of interacting binaries will experience only stable MT before the first SN if f_γ > 0.5, with that fraction dropping quickly with increasing β or f_γ (see Figure 4). If very high f_γ values are more representative of reality, this would restrict the viability of the stable MT channel for the formation of binaries containing a compact object (in agreement with van Son et al. 2022 in the context of binary BH formation). This could potentially have a strong impact on the predicted formation rates of these binaries, and thus, conversely, could be used as a constraint on f_γ . However, such impacts are outside of the scope of this study as they invariably require modeling of the CEE, SNe, and other physics which we have explicitly ignored here.

MT stability prescriptions based on the delayed onset of MT instability in both radiative and convective envelopes, such as those from Ge et al. (2020), show that stable MT is significantly more common than is suggested by prescriptions that are restricted to the onset of MT. However, these more nuanced prescriptions are computed based on the assumption of fully conservative MT. To the extent that this is not an accurate representation of most binary interactions, this assumption underpredicts the likelihood of instability.

In Figure 2, we show the prescriptions for ζ_ad derived from the detailed models from Ge et al. (2020), ζ_Ge20 and ${\tilde{\zeta }}_{\mathrm{Ge}20}$, as well as the default ζ_Comp used in COMPAS and ζ_Claeys14 from Claeys et al. (2014). ζ_Ge20 and ${\tilde{\zeta }}_{\mathrm{Ge}20}$ are generally very similar, but both are higher than ζ_Comp and ζ_Claeys14 for virtually all masses and radii, indicating that stable MT may be more abundant than is currently estimated, variations to β and γ notwithstanding.

Notably, in response to Ge et al. (2010), Woods & Ivanova (2011) argue that the thermal timescale in the superadiabatic surface layers is comparable to the dynamical timescale, thus invalidating the adiabatic assumption of Ge et al. (2010, 2015, 2020). They predict that the thermal contraction and restructuring of the superadiabatic layer occurs faster than the mass can be removed, except for unrealistically high mass-loss rates. This prediction is supported by Pavlovskii & Ivanova (2015), who find no evidence for rapid expansion when recombination energy is included in the superadiabatic layers. They argue instead that matter flowing through the nozzle around the L1 point is constricted, and that MT instability only sets in once stellar material begins to overflow through either the L2 or L3 points (Pavlovskii & Ivanova 2015; Pavlovskii et al. 2017). Lu et al. (2023) explored the threshold conditions for mass loss from the L2 point, and find that this occurs only for fairly high donor mass-loss rates ≳ 10⁻⁴ M_⊙ yr⁻¹, though this is dependent on the geometry of the system and the assumed cooling rate of the accretion stream. This is marginally consistent with Ge et al. (2010), who find that MT instability sets in around ∼10⁻⁵ M_⊙ yr⁻¹ but quickly rises above ∼10⁻⁴ M_⊙ yr⁻¹ (see their Figure 5), though a full population synthesis study to compare the models properly is warranted.

In a recent paper, Temmink et al. (2023) compared several different stability definitions, including a critical mass-loss rate at which the donor response becomes adiabatic, overflow through the L2 point, and rapid contraction of the orbit. They found that MT is generally more stable than is predicted in the prescriptions commonly in use in BPS codes, although they also computed their stability criteria assuming fully conservative MT, and their study focused on a lower mass regime than is considered here. Crucially, they conclude that none of the definitions for stability are self-consistent—or particularly reliable—for very evolved stars, which suggests that the difference between stable and unstable MT itself may not be well defined, or perhaps that the transition between the two regimes is not as discrete as has been assumed historically.

In practice, the convective–radiative nature of an envelope is not directly aligned with its stellar type as defined in Hurley et al. (2000). Convective instabilities develop in cool envelopes where the opacity is sufficiently high to make radiation transport inefficient. We implemented an improved prescription in which a giant's envelope is radiative if the surface effective temperature is greater than a given threshold temperature, and convective otherwise. We considered several values for this threshold temperature, including $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})=3.73$ from Belczynski et al. (2008) and $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})=3.65$ from Klencki et al. (2021). However, we found that these variations had only a marginal impact on MT stability, and thus did not pursue them further. This may suggest that the convective–radiative boundary is not so uncertain as to impact the stable MT fraction substantially, however population synthesis with more detailed stellar structures and responses to mass loss will certainly provide some clarity here.

Observations of interacting binaries are crucial in constraining the parameters β and γ. Binaries engaged in ongoing MT can provide vital constraints on the instantaneous values of β and γ in a given mass-exchange episode, via measurable changes to the orbital period. In tandem, the properties of observed populations of specific classes of postinteraction binaries, or their associated transients, can also be informative in determining the values and universality of β and γ in different contexts, though these are subject to modeling constraints.

For an in-depth review of the possible classes of interacting binary observations, we refer the reader to De Marco & Izzard (2017), particularly their Table 1 and the references therein. However, we emphasize that the most useful observations will be those which involve binaries which have experienced stable MT recently in their evolution (e.g., Algols), to reduce modeling uncertainties. By contrast, binaries which are believed to have survived an SN or common-envelope event, must necessarily be convolved with models for these more uncertain phases, reducing the viability of any claimed constraints. The choice taken in this paper to study the stable MT rate only up to the first SN reflects this preference to reduce the modeling complexity in order to derive more robust conclusions.