Mass transfer from giant donors
Abstract
The stability of mass transfer in binaries with convective giant donors remains an open question in modern astrophysics. There is a significant discrepancy between what the existing methods predict for a response to mass loss of the giant itself, as well as for the mass transfer rate during the Roche lobe overflow. Here we show that the recombination energy in the superadiabatic layer plays an important and hitherto unaccountedfor role in the donor’s response to mass loss, in particular on its luminosity and effective temperature. Our improved optically thick nozzle method to calculate the mass transfer rate via allows us to evolve binary systems for a substantial Roche lobe overflow. We propose a new, strengthened criterion for the mass transfer instability, basing it on whether the donor experiences overflow through its outer Lagrangian point. We find that with the new criterion, if the donor has a welldeveloped outer convective envelope, the critical initial mass ratio for which a binary would evolve stably through the conservative mass transfer varies from to , which is about twice as large as previously believed. In underdeveloped giants with shallow convective envelopes this critical ratio may be even larger. When the convective envelope is still growing, and in particular for most cases of massive donors, the critical mass ratio gradually decreases to this value, from that of radiative donors.
keywords:
instabilities — methods: numerical — stars: massloss — binaries: close — stars: evolution1 Introduction
Many interacting binaries start their mass exchange when the donor, which has evolved to the giant branch, overfills its Roche lobe. The stability at the start of the mass transfer (MT) differentiates between the binaries that will live long as a mass transferring system and will appear, e.g., as an Xray binary, and ones that will be transformed dramatically by a common envelope event. A clear separation of the possible evolutionary channels is important for our understanding of the formation of compact binaries. E.g., double whitedwarf binaries were thought in the past to be formed via two common envelope events (Nelemans et al., 2000; Nelemans & Tout, 2005), while more recent detailed simulations of MT with lowmass red giant donors show that the first episode of MT was stable (Woods et al., 2012).
There are cases when MT from a main sequence star, while deemed to be stable initially, increases its rate significantly with time, in some cases leading to socalled delayed dynamical instability (for a thorough discussion, see Ge et al., 2010). A delayed rapid growth of MT rate can also occur in systems with giant donors(Woods et al., 2012). Unfortunately, to date there exists no simple criterion for the start of a dynamically unstable MT that would be based just on its rate. It can be defined as MT happening on the dynamical timescale of the donor, but even this definition does not determine the fate of a system (for a review of various problems with the initiation of a common envelope event, see Ivanova et al., 2013).
Standard definition of a dynamicaltimescale MT is that the mass is transferred on the donor’s dynamical timescale (see, e.g., Paczyński & Sienkiewicz 1972). When a star is subjected to a dynamicaltimescale mass loss (ML), stellar codes operate outside of the regime in which they were designed to operate, leading to poorly understood behaviour. Indeed, a conventional stellar code, developed to treat evolutionary changes, typically does not include hydrodynamical terms in structure equations, and hence is unable to treat dynamicaltimescale features of the ML response.
There are then two ways to deal with a dynamical ML in practice, both of them relying on plausiblesounding but arbitrary assumptions.
One way is to preset a threshold for the MT rate close to (but still within) the limits of validity of a conventional stellar code, and assume MT to become unstable once that threshold is exceeded. For example, a dynamical MT has been defined as 10 times the thermal timescale MT (Nelson & Eggleton, 2001).
Another way is to assume that during a very short period of very rapid mass loss (ML), any mass shell in the stellar model has no time to exchange heat with its neighbouring shells. Hence, its entropy remains almost the same as at the start of the ML, in the other words the specific entropy profile of the donor is “frozen”.
The MT rate and stability depend crucially on the response of the donor’s radius to the ML compared to the response of its Roche lobe radius to the same ML (for foundations, see Webbink, 1985a). The governing feature that shaped studies of MT in binary systems with giant donors was the wide acceptance of the theory that a convective donor would intrinsically expand as a reaction to ML. In the framework of this theory, the immediate response of a star to ML takes place on a dynamical timescale (to restore its hydrostatic equilibrium) – this assumes that no thermal effects can take place, and it was hence called an adiabatic response. It is characterised by the quantity known as the massradius exponent
(1) 
With polytropic models (including also composite and condensed polytropes), it was found that a convective donor would have (Paczynski, 1965; Hjellming & Webbink, 1987). The presence of a nonconvective core increases the value of , which becomes positive once the relative mass of the core is more than 20% of the star’s total mass (Hjellming & Webbink, 1987).
It is commonly assumed that MT is dynamically stable only if a donor remains within its Roche lobe, or when , where is the massradius exponent for the Roche lobe reaction to the MT(Webbink, 1985b). Using the Roche lobe approximation as in Eggleton (1983), for conservative MT, when the mass ratio is , and is positive for a larger mass ratio. For completely nonconservative isotropic reemission this critical mass ratio is (Soberman et al., 1997).
Hjellming & Webbink (1987) showed that taking into account the core improves stability, however for all stars with relative mass of the core of the total stellar mass, the first episode of conservative MT will be dynamicall unstable. Later, the results based on polytropic models were reevaluated in the studies of the detailed adiabatic stellar models by Ge et al. (2010). At the same time, detailed codes that traced mass transfer in binaries up to thermal timescale mass transfer, found that critical mass ratio can be up to 1.1 (e.g., Han et al., 2002).
It is crucial to realize that the difference in rates between a nearlythermal timescale MT that can be legitimately calculated with a conventional stellar code and a dynamicaltimescale MT – the regime of adiabatic codes – is several orders of magnitude, and currently there is no adequate treatment to describe the donor response in between these two regimes.
A detailed stellar code can be supplied with hydrodynamical terms in its structure equations and be forced to work in a hydrodynamical regime (please note that adiabatic codes do not include hydrodynamical terms in them). This way, in the recent study of responses of stellar models of giants to fast ML, it was found that this response is a function of the MT rate (Woods & Ivanova, 2011). As a result, it was found that for a large range of ML rates taking place in binaries, giants can also contract, and the nature of the donor’s response was attributed to the behaviour of the giant’s superadiabatic surface layer and its short thermal timescale. In that study, two different stellar codes were used for comparison and they showed the same results. The two codes were: a) STARS/ev code developed by Eggleton (Eggleton (1971, 1972, 1973); Eggleton et al. (1973)) with the most recent update described as in Glebbeek et al. (2008); b) the Heyneytype code developed by Kippenhahn et al. (1967) with the most recent update described in Ivanova & Taam (2004). These results were later confirmed by Passy et al. (2012), with yet another stellar code, MESA.
As can be seen, three different detailed stellar codes appeared to produce the same result – the donor response is very different from the expectations given by a simplified adiabatic model. However, no further analysis of why there is a dramatic difference between an adiabatic approach and a detailed stellar code was provided. The important question remains: which result we should trust more? What could be the possible problems with the adiabatic approach? Could we trust the results of the detailed stellar codes when they are forced to work at dynamicaltimescale MT rates, and when they possibly break down?
It is clear that for further progress in understanding of the response of the giant’s radius to fast ML, one needs to understand better the nonadiabatic processes taking place in a red giant subjected to dynamicaltimescale ML. This is intrinsically related to the nature – the underlying physics and structure – of the superadiabatic layer, which is not yet well understood. A possible influence of the surface layer in giants on MT was realised in the past. For example, Osaki (1970) had suggested that, in U Gem giants, the mechanism for dynamical mass exchange, proposed by Paczynski (1965), can not work, “because of changes occurring on the thermal time scales of the photosphere and transition zone, which are shorter than the dynamical time scale of the star”). However, the physics of the superadiabatic layer had never been studied in detail. On the other hand, the difference in the response in different codes also could be related to numerical effects, which could be either inconsistency in the adopted set of equations, or numerical artefacts.
In this paper, in §2, we systematically analyse the physics of the superadiabatic layer: what difference an adiabatic approach makes, which artefacts might appear if this layer is not numerically modelled properly, what affects the structure of this layer, how massive this layer is in different donors, and finally, whether we can numerically obtain the response of the stellar radius properly with a detailed stellar code, and with what limitations.
All detailed stellar models and numerical experiments in this paper are calculated with MESA (Modules for Experiments in Stellar Astrophysics, http://mesa.sourceforge.net). This set of stellar libraries is described in Paxton et al. (2011, 2013). We chose this modern stellar package, as, due to ease of use, it is becoming widely popular, and may soon become the most commonly used stellar tool. For the purpose of comparison, we clarify that, by default, we use solar metallicity, and employ a standard set of stellar equations: the simple grey atmosphere boundary condition with radiation pressure at zero optical depth taken from Cox & Giuli (1968), mixinglength theory (MLT) as in Cox & Giuli (1968), and radiative opacity tables based on OPAL and Ferguson et al. (2005). If there is any deviation from this default set, it will be clarified for each numerical experiment. While some numerical issues that we will discuss in this paper are relevant only to our choice of the stellar code, most of the issues are relevant to modeling of fast ML with any detailed stellar code. We are aware that MESA was developed to evolve normal stars, in a conventional way, and in this paper we will analyse how MESA, and possibly similar codes, work at an extreme beyond its design.
While the response of the donor’s radius itself to ML is the first important issue, the second foremost is how fast stellar material, overfilling the donor’s Roche lobe, can be transferred to its companion, and hence be effectively lost from the donor. In §3 and §4 we discuss treatment of the ML in detailed stellar models and describe our adopted model of an optically thick stream appropriate to a case of fast ML. In §5 we provide the results of our ML simulations. This includes a comparison with the previously published results, and analysis of the stability of the MT in a binary system.
2 Surface superadiabatic layer
A major feature distinguishing adiabatic calculations of the MT from the results of detailed stellar codes for convective donors is the retention of a thin (in mass sense) layer of lowered entropy near the surface. Understanding the physics of this layer is likely to be key to understanding a realistic stellar response to ML.
2.1 Definition
The term “superadiabatic” makes sense only within the adopted convection theory, MLT.
In MLT, the relative efficiency of heat transfer via two mechanisms (radiative and convective) can be characterized by the ratio of radiative and convective conductivites (see equation 7.12 in Kippenhahn & Weigert, 1994):
(2) 
where is the mixing length and is pressure scale height. is a quantity that characterises the thermal expansion properties of a medium.
In the framework of MLT, convective efficiency directly depends on two free parameters – on how far a blob rises or sinks before it mixes with the environment (which can be seen in Equation 2 directly as ), and on the linear size of convective blobs (used in Equation 2). These free parameters are best calibrated from observations and are not obtained from physical considerations.
In convective zones overadiabaticity, which is the difference between the actual and adiabatic temperature gradients, , is positive. We define the superadiabatic layer as the region where overadiabacity is comparable to the gradients themselves.
In this paper, for the efficiency of convection, we adopt the relation between conductivities of convective and radiative heat transfer denoted by . Larger values of correspond to less efficient convection and accordingly to a larger overadiabaticity. In red giants, due to a combination of factors increases in the vicinity of the stellar surface to the extent that the real gradient of temperature becomes substantially detached from the adiabatic gradient () (see Figure 1).
2.2 Superadiabatic layer and the effect of partial ionisation
Let us consider in detail what affects the behaviour of the efficiency of the convection and of the overadiabaticity . In the region where both quantities have simultaneous slowdown in their almost monotonic decrease towards the center (see Figure 2), hydrogen and helium change their degree of ionisation (see Figure 3). The degrees of ionisation affect and by strongly changing the heat capacity of the medium and the Rosseland mean opacity (see Figure 4).
For the range of densities and temperatures in the envelope, the heat capacity of a partially ionised gas is higher than that of both neutral and fully ionised gases (ionisation consumes a part of the heat transferred to the gas without increasing the kinetic energy of its molecules). The three major jumps in , from highest to lowest temperature, are associated with the following recombinations:

double ionised helium single ionised helium;

single ionised helium neutral helium;

ionised hydrogen neutral hydrogen.
Let us discuss the reason for two overadiabaticity plateaus in the area of interest marked with the solid opacity curve in Figure 4. In this area, partial recombination of elements affects the convective efficiency mainly through , because opacity does not vary much there. Due to the change in monotonicity of at halfionisation, recombination of either hydrogen or helium impedes the growth of overadiabaticity towards the surface while the ionisation fraction is between and , but when the ionisation fraction is lower than the recombination, on the other hand, accelerates the growth of overadiabaticity (see Equation 2 and Figure 4). Because of this, the areas where the growth of overadiabaticity is impeded (plateaus) are shifted with respect to the maxima of . The initial stage of partial recombination of hydrogen is able to completely suppress the growth of overadiabaticity towards stellar surface in the envelope; the initial recombination of helium also impedes it to a large extent. The shape of the superadiabatic layer, hence, is substantially affected by the recombination of hydrogen and helium.
At the final stages of hydrogen recombination, the opacity falls sharply by a few orders of magnitude. As a result, the radiative gradient, which is proportional to opacity, falls below the adiabatic gradient (Figure 1). We call this region the ”surface radiative zone”. It takes place outside the temperature range shown in Figures 24, at . The spatial size of the surface radiative zone is of the order of the mixing length. Hence it is plausible that convective blobs overshoot the boundary between the superadiabatic layer and the surface radiative zone (see also Kuhfuss, 1986).
2.3 Mass of superadiabatic layers in convective donors
For analysis of superadiabatic layers and radiative zones, we introduce the quantity , the mass of these layers. We measure from the stellar surface to the point near the start of the entropy drop, . Since the convection is never fully adiabatic, the entropy profile in the envelope is not flat, and there is no sharp transition between the entropy plateau and the entropy drop (see Figure 5). For integrity between stars of different masses, we therefore choose to define as the point where entropy is
(3) 
Here is the local minimum of entropy in the outer layer of the star and is the (maximum) value of entropy in the convective envelope.
Clearly, it can be seen in Figure 6 that the mass of this shell varies between giants of different masses and radii. From the behaviour of the envelope’s simplified analytical solution (see e.g.Figure 10.2 in Kippenhahn & Weigert, 1994), we expect that overadiabaticity of convective envelopes becomes substantial at about the same pressure, to an order of magnitude. Indeed, in realistic stellar models for stars of different masses, the pressure at the bottom of the superadiabatic layer does not vary much from dyn/cm. In this case, we can estimate by simply considering a pressure drop. The equation of hydrostatic equilibrium over the superadiabatic layer and radiative zone is
(4) 
As the pressure at the top of the superadiabatic layer (photospheric pressure) is much less than at its bottom, we find
(5) 
Indeed, as illustrated in Figure 6, stars with convective envelopes at any evolutionary stage have which can be found using this relation. This approximation has only limited, though still very good, applicability to more massive stars, where is comparable to the total mass of the star’s envelope (we note also that it is well known that in massive stars almost the entire convective envelope could be superadiabatic).
It is important to note that this is the mass of the superadiabatic layer and surface radiative zone in a star unperturbed by ML. Mass loss affects both our assumptions about the pressure drop (as the star is not in hydrostatic equilibrium anymore), and the entropy profile.
We see from the preceding analysis that while the mass of the superadiabatic layer is relatively small in lowmass giants, it grows for more massive stars, increasing substantially as the star is evolving and expanding. Stellar models that have a constant entropy profile are therefore not justified, as the error in the case of massive stars could be enormous. We will also show in §2.4 that the mass of the superadiabatic layer helps estimate its rate of energy generation for each ML rate, and how to find the timesteps that provide valid results for the radial response on the ML.
2.4 Energetics of the superadiabatic layer
We saw a hint in § 2.2 that the formation of a superadiabatic layer is to a certain extent governed by recombination, hence this layer might be energetically important for obtaining the correct ML response. Energetically, recombination enters the structure equations through the gravitational energy term, which is the last term in the energy generation equation and is defined as:
(6) 
The term can be found directly from the Equation 6, using only one temporal derivative of the entropy which is derived from the EOS. Alternatively, can be calculated by subtracting two terms that depend on temporal derivatives of thermodynamic quantities, for example:
(7) 
When mass is removed from the surface of a red giant, the partial ionisation zones of hydrogen and helium discussed in § 2.2 move to areas that were previously fully ionized, deeper in the star. Recombination of these previously ionised areas leads to the release of recombination energy. One can estimate that the recombination energy of pure hydrogen is erg g. For the ML rate of yr the recombination energy release should be of the order of erg s.
Let us consider now this mass loss on an example of a and red giant. Using our Equation 5 for the mass of the superadiabatic layer, we find that its total mass is about , though the part where hydrogen recombination energy is released is usually only a fraction of it, about (see also Figure 7). The contribution of this recombination component alone in in such a superadiabatic layer is expected to be of the order of erg g s.
We have checked what values of we obtain in practice, using the detailed stellar code of our choice, MESA. First, we tested the method where the entropy derivative is found using Equation 6, as it takes into account composition changes, unlike the second method that uses Equation 7. However, this method suffers from an inexact calculation of entropy from the EOS. After performing numerical tests, we found it to not be suitable for the small timesteps necessary for fast mass loss calculations: as the timesteps become smaller, the errors in the entropy values become comparable to their Lagrangian differences between subsequent timesteps. The twocomponent formula 7 always uses structural variables. Entropy, used in the formula 6, is never a structural variable: it is always derived from the structural variables through the tabulated EOS. In addition, composition artefacts can play a role in entropy fluctuations.
Second, we tested the method where the entropy derivative is found using Equation 7, see Figure 7. We found that the expected physically reasonable level of is only reached when the timestep is below yr. We relate it to the inaccuracy of the firstorder numerical differentiation formulae used in practice to calculate the righthand side of (7) at larger timesteps. Note that the righthand side of (7) is effectively proportional to the Lagrangian derivative of entropy^{1}^{1}1If one neglects composition changes., which is calculated indirectly through the corresponding derivatives of and . Therefore, the error in calculation of with formula (7) is, as expected, the most significant in the superadiabatic layer and surface radiative zone, because the second Lagrangian derivative of entropy is the highest there^{2}^{2}2 We’d also like to mention that quite recently a new scheme for the calculation of was implemented in MESA, that follows a mixed LagrangianEulerian approach. We performed a preliminary testing of this scheme on masslosing red giants. Unfortunately, this new scheme is not suitable for our simulations yet because of severe numerical artefacts that it introduces under certain conditions. These artefacts have to be carefully studied and eliminated before we can consider incorporating this new scheme into our simulations. Because of this, we only discuss the standard purely Lagrangian scheme for the calculation of in this paper..
We do not suggest a complete formal procedure to calculate the errors in calculation of . Instead, we resort to a plain comparison of results obtained with various timestep selection approaches. For example, the comparison for one of our models (Figure 7) shows that if it loses not more than of its initial in one timestep, then the calculation of in the superadiabatic layer is affected only marginally. Removing mass in such small timesteps is, however, quite resource intensive. We can foresee that the outlined numerical problem with a recombination zone that moves too fast in mass can be eventually resolved using a technique similar to the one used for calculating thin shell burning as in Eggleton (1967, 1973). However, for the purpose of the studies presented in this paper, we can afford to use a numerically intensive way.
In Section 5.1 we will show how inaccurate calculations of affect the radial response to the mass loss. While the numerical problems that we described in this section are native only to the code we use, when fast mass loss rate calculations are performed with another code, the numerically obtained has to be tested against its expected value that can be found using the method described in this subsection.
3 Fast ML in an onedimensional star
There currently does not exist a comprehensive and selfconsistent hydrodynamical simulation to treat the threedimensional problem of MT in a binary. Even when threedimensional simulations eventually become selfconsistent, they would likely resolve only the question of the initial stability of the Roche lobe overflow (RLOF), but not the longterm MT, which will remain a prerogative of onedimensional stellar codes in the foreseeable future.
In a standard onedimensional ML model, the stellar codes use the regular set of structure equations, but adopting the boundary condition that the total mass decreases with time. This boundary condition is an unavoidable reduction of the threedimensional picture of the ML to one dimension. In other words, represents our best understanding of the stream that is formed in the vicinity of and that carries the donor’s material away from the donor. In this Section we examine the reaction of giants to the ML in onedimensional stellar codes, explaining in particular the nature of the feature observed in the previous MT calculations.
3.1 Understanding the initial contraction of a red giant upon the instantaneous ML in 1D stellar codes
Two completely opposite responses to ML were found in detailed 1D simulations by Passy et al. (2012). According to them, hydrostatic stellar models expand in response to ML, while models with hydrodynamical terms, on the opposite, shrink.
In neither of these two approaches the ML experiment is close to what would happen in Nature, where MT never starts abruptly at some fixed high MT rate. However, it is important to understand what causes the dramatic difference between these two approaches.
Passy et al. (2012) provided the following explanation: ”…some energy that is stored in gravitational form in the hydrostatic models is actually in a kinetic form [in hydrodynamic models], leading to the star contracting instead of expanding”, although how exactly the transformation of gravitational energy into kinetic leads to contraction was not explained.
Instead of the energy argument, we argue that the main reason for the radial response is due to the material being consumed from the surface of a red giant with linear velocity , where is the radius of a red giant, and is the surface density. Note that this term intrinsically decreases the radius.
To validate that this term is dominant, we need to consider the involved stellar equations. A standard way to introduce hydrodynamical treatment into a stellar code is to use a truncated NavierStokes equation in spherical coordinates to calculate the hydrodynamical term of the pressure derivative, for example:
(8) 
where the acceleration term is the ratio of the local Lagrangian acceleration to the local gravitational acceleration:
(9) 
where is the local Lagrangian acceleration. For a star in hydrostatic equilibrium, . In Figure 8 we illustrate how significant the acceleration term can be in the superadiabatic layer and surface radiative zone of a red giant at high ML rates.
Let’s now consider the initial response after some (small compared to the dynamical timescale) time to ML of a red giant of initial mass without artificial viscosity or any other effects that would alter Equation (8).
If is a smooth and bounded function and and are continuous and bounded, then
(10) 
By integrating Equation (10), we get for the velocity at the surface^{3}^{3}3Here we use the ”littleo” notation for the vicinity of zero. A brief definition of this notation is (for details see, e.g., Kevorkian & Cole, 1985).:
(11)  
where . To obtain the complete radial response, this velocity must be combined with the material consumption velocity, which is equal to
(12) 
where is the initial surface density of the star. Their sum must in turn be integrated in time. Thus, the complete radial response of a red giant to ML is given by:
If the original star was in hydrostatic equilibrium, then , and, according to Equation (10), the term, proportional to in Equation (3.1), would be equal to , and hence is vanishing. The total radial response reduces to
(14) 
This result may be directly compared to the decrease in the radius of a stellar model that occurs during the first timestep after the start of the ML, divided by this timestep. This comparison helps to understand the role of the adopted outer boundary condition (see Sections 3.2 and 5.1). We conclude that the contraction observed in the simulations of Passy et al. (2012) is neither a numerical error, nor a conversion of gravitational energy to kinetic, but a natural consequence of a finite pace of expansion of the onedimensional envelope that was forced to experience a fast ML.
3.2 The acceleration term
As shown in Figure 8, the acceleration term is becoming dangerously large in the surface layers during fast ML. In this Section we describe the modification of stellar equations and the boundary condition to take this term into account.
To obtain the pressure at the outer boundary in a simple grey atmosphere approximation, it is common to use the following wellknown formula^{4}^{4}4Here for simplicity we omit radiation pressure at zero optical depth. This term may be important if radiative pressure is a substantial fraction of the total pressure.:
(15) 
Unfortunately, the above formula is based on the equation of hydrostatic equilibrium, hence in our case it is not applicable, not to speak that it is inconsistent with the Equation 8 that is used in many stellar codes. For this reason we took into account the acceleration term not only in the equation of motion (8), but also to find the boundary condition for a simple grey atmosphere, which we use in our experimental version of MESA:
(16) 
Numerical experiments with MESA showed that the choice of the surface boundary condition affects the resulting radial response. The modified boundary condition must be used because in this case we obtain a radial response substantially closer to that predicted by Equation 14 (see Section 5.1).
Similarly, we take into account the acceleration term for the temperature equation
(17) 
for convective conductivity in the Mixing Length Theory, and for the radiative gradient:
(18) 
4 Rate of the mass transfer
We describe below an optically thick model for MT that we have adopted for our studies of hydrodynamic response to rapid MT in binaries.
We follow a conventional way to calculate the MT rate by integrating the mass flow over a “nozzle” crosssection that is taken on a plane perpendicular to the line connecting the centres of the two stars and passing through the point:
(19) 
Here is the velocity that the stream has within the nozzle, and is the density that the stream has at the same position. We note that the use of the nozzle crosssection only in the neighbourhood implies that only the MT rate can be found, and the rate of ML via overflow cannot be calculated.
Furthermore, any scheme that finds an ML boundary condition for stellar onedimensional codes adopts some set of simplifications to find the distribution of the stream’s density and velocity throughout the nozzle, as well as the nozzle geometry. The flow is considered to be steady, which permits the use of the Bernoulli theorem along the streamlines. The standard assumptions are that the nozzle at coincides with the sonic surface of the flow (Lubow & Shu, 1976), and that initial velocities are negligible at the origin of the streamlines. We adopt the same assumptions.
The MT rate then becomes dependent only on these assumptions:

the adopted geometry of the donor and the nozzle;

the streamlines – this assumption, coupled with the adopted evolution of the specific entropy along the streamlines, allows relating the donor’s thermodynamical properties to those of the flow crossing the nozzle.
4.1 RLOF formalism and geometry of the problem
The fundamental simplifying assumption that governs the whole Roche lobe formalism is volume correspondence. More precisely, this is the assumption that thermodynamical parameters (pressure, density, composition, etc.) at a certain radius in a onedimensional stellar model are the same in threedimensional space at a Roche equipotential whose enclosed volume is . Whenever onedimensional stellar evolution is considered in terms of RLOF formalism, the volume correspondence assumption is automatically used. Note that when a donor experiences a substantial RLOF this volume correspondence for thermodynamical parameters can be applied only far from the point, where the donor’s material is almost at rest. The Roche lobe volume radius is usually found using one of two wellknown approximations (Paczyński, 1971; Eggleton, 1983). The area of the nozzle is then usually approximated by the area of an ellipsoid by taking the secondorder term in the Roche potential expansion within the sonic surface near the point.
An additional simplifying assumption is applied when integrating the ML rate over the potential region between the Roche lobe potential and the photospheric Roche potential (note, that ). At this point, to avoid volume integrations, it is a common approach to implicitly break the volume correspondence assumption and replace it with the pressure correspondence assumption.
The pressure correspondence assumption asserts that the thermodynamical parameters (temperature, density and composition) at a point A in threedimensional space are the same as at a point B in the onedimensional model, provided that the pressure at the point A in threedimensional space is the same as the pressure at point B in the onedimensional model. If one applies this pressure correspondence far from the point and assumes, in addition, that
(20) 
where a is the local Lagrangian acceleration, and is the local Roche potential, then it becomes possible to replace the integral over radius with the integral over pressure. The benefit of the pressure correspondence assumption is that one can avoid the painful calculation of the differential , which requires volume integrations. We note that once a star is not in hydrostatic equilibrium, does not hold (see Section 3.1), and condition 20 cannot be satisfied. Hence, pressure correspondence should not be used for rapid mass loss.
RLOF during the evolution of a masslosing red giant can however be very substantial. When an equipotential surface is far from the Roche lobe surface, the expansions obtained in the vicinity of the point break down since the expansions are truncated after the quadratic term. Higherorder terms are no longer negligible. Instead, we conduct the realistic Roche lobe integrations which employ the RungeKuttaFehlberg integrator and the damped NewtonRaphson solver to obtain all geometrical parameters. These integrations have been conducted for 275001 mass ratios from 0.06 to 19.145 with an increment of . For each mass ratio we calculate the nozzle areas and volume radii for 200 equipotentials lying between the potential and the potential. In addition to the use of a more precise relation between the star’s volume and the Roche lobe volume, volume correspondence that is necessary for rapid mass loss and a nonsimplified nozzle shape, this also allows us to track whether the donor overfills the equipotential during the mass transfer.
4.2 Streamlines
For a rapid mass transfer rate, an ”optically thick” approximation is usually used (see Paczyński & Sienkiewicz, 1972; Savonije, 1978; Kolb & Ritter, 1990; Ge et al., 2010, and many others). The flow of matter towards in this approximation is adiabatic and streamlines go along the equipotential surfaces of the Roche potential. The photosphere corresponds to a Roche equipotential which lies outside the Roche lobe. Close to the point, the photosphere turns into the outer boundary of the optically thick flow which flows across the sonic surface along the Roche equipotential surfaces with the local adiabatic speed of sound (Kolb & Ritter, 1990).
To find the stream’s density and sonic velocity at , one is required to adopt the evolution of the specific entropy in the flow along the streamlines. Often the flow is taken to be polytropic – that is, the flow preserves a constant value of along a streamline, where is the adiabatic exponent. In certain cases, the polytropic stratification of the donor itself is adopted. Note that the stream’s velocity, while locally sonic, is not constant across the nozzle. The flow also does not need to be isentropic if the donor is not isentropic – the streamlines can originate at different initial equipotentials.
We also assume that a flow is adiabatic. An adiabatic flow is not, of course, entirely polytropic due to, for example, recombination, that occurs as gas flows towards the sonic surface. It means that varies along an adiabat. For this reason, we employ a realistic equation of state taken from MESA, that among other effects takes into account the ionisation of elements in the mix and the radiative component of pressure. For test purposes, we also can consider the flow to be polytropic. We have verified that for those models that are provided as examples in this paper, the effect of the adopted equations of state on mass flux does not exceed .
We also should mention that in our scheme we do not use any approximation for an optically thin MT stream that originates from above the photosphere. The primary reason for this is that we have not yet developed a model that would allow one to combine both optically thick and optically thin ML schemes when a donor overfills its Roche lobe. Kolb & Ritter (1990) have suggested that during RLOF, the total MT rate can be found by a summation of the MT of an optically thick stream, calculated as a function of the current RLOF, plus the maximum MT rate obtainable for an optically thin photospheric flow in the case of RL underflow, which is at the instant when a star exactly fills its Roche lobe. We are not confident whether this method would estimate correctly the contribution of an isothermal photospheric outflow in the case of a nonnegligible RLOF, considering that both the geometrical crosssection and the potential of the nozzle for the isothermal flow (located now around the nozzle for the adiabatic flow) are changed substantially from those in the neighbourhood. Hence, we do not use this approach.
The degree to which optically thin mass transfer affects the behaviour of stars before their photospheres overfill the Roche lobe depends on the rate of their intrinsic evolutionary expansion when their photospheres approach the Roche lobe. Relatively massive giants (5 and up), cross the interval of radii for which the optically thin mass transfer is dominant in a short time thanks to their fast expansion. Hence, the fraction of the envelope they lose via optically thin transfer is small. On the other hand, less massive giants epxand slowly and can lose a lot of mass through atmospheric ML before the optically thick mechanism kicks in. This decreases the mass ratio at the onset of RLOF and improves stability. We observed an extreme case, where with the Ritter (1988) prescription, a giant could never actually overfill a Roche lobe, and steadily lost the whole convective envelope via optically thin mechanism reaching an optically thin MT rate of ! We therefore warn the reader that the mode of mass transfer from lowmass giants might be sensitive to the very uncertain model of optically thin MT.
We summarise the simplifications eliminated in the existing optically thick schemes and in our case in Table 1.
Reference  GS  PC  PD  PS 

Paczyński & Sienkiewicz (1972)  
Savonije (1978)  
Kolb & Ritter (1990)  
Ge et al. (2010)  
This work 
GS – geometrical simplification for the nozzle, PC – pressure correspondence, PD – polytropic stratification of the donor, PS – along the streamlines is constant. Empty circle – simplified potential is adopted the effect of which is equivalent to the PC assumption.
5 Results of ML simulations in red giants
5.1 Response to constant ML
First, we can compare the initial reaction of a red giant to constant ML to our prediction made in Section 3.1. For example, we take a and giant and subject it to the constant ML rate yr. The initial derivative of radius predicted by formula (14) is about . With modified boundary conditions and with a timestep of yr, the code produces . If instead the default hydrostatic boundary condition is used, the code produces .
Second, to compare our results to the ones published earlier, we calculate the behaviour of the and red giant examined by Passy et al. (2012). As in Passy et al. (2012), we subject this giant to the constant ML rate yr. We evolve one giant using the unmodified MESA code with viscosity and timestep setup described in Passy et al. (2012); we call this the “passya” model. For the other giant, “passyb” model, we take into account the circumstances outlined in Section 2.4. In particular, to obtain correct , we evolve this model using a constant time step of yr. We also remove some pulsation artefacts by disabling composition smoothing at the bottom of the convective envelope.
The initial behaviour of the radius in both models has been explained in Section 3.1. Note, however, that in our model, the stellar radius is always smaller than that in Passy et al. (2012) (see our Figure 9 and their Figure 4). We find that the difference in the value of the radii is mainly due to how accurately the value of is found.
In addition to the initial radius contraction, one can notice ensuing radial oscillations, visible both in the models passya and passyb. We define those nonnumerical radial pulsations, excited either by the start of constant ML or, in the case of evolution in a binary, by the rapid growth of the ML rate, as mass loss induced pulsations (MLIPs). These pulsations might be caused by the sonic rarefaction wave, reflected from the bottom of the envelope, which is theoretically predicted to occur in a fluid as it abruptly expands into vacuum, hence we think that they are not of numerical nature. In the model passya, these pulsations are largely smoothed out because the timestep in this model grows after the start of simulations and exceeds the dynamical timescale.
A crucial parameter that defines the level of importance of MLIP is the pmode damping rate , which characterises the timescale on which a giant roughly attains hydrostatic equilibrium and dynamical oscillations are damped. With contemporary highprecision photometric instruments such as COROT and Kepler, it is possible to obtain highprecision measurements of the profile widths (i.e, damping rates divided by ) of pmodes (see e.g. Baudin et al., 2011; Belkacem et al., 2012). Damping rates of intrinsic pulsations of red giants found from those observations are to .
The damping rates, given that the timestep used is much smaller than the pulsation period, depend hugely on the artificial viscosity coefficients: increasing the artificial viscosity increases the damping rates. The damping rates of MLIPs shown in Figure 9 are comparable, within an order of magnitude, with those obtained directly from observations. Note that the default artificial viscosity used for the models passya and passyb is very moderate (, ) and does not substantially affect the damping rates. Due to the lack of widerange observational calibrations for damping rates across giants of different radii and masses, we do not use an artificial viscosity in our models, except in models passya and passyb, where it is taken into account only for the purpose of comparison with the original paper of Passy et al. (2012).
We have performed several simulations where we subjected giants of several initial masses and radii to constant ML, to determine the realistic stellar response defined as:
(21) 
Note that generally is an implicit function of not only the current MT rate, but also of the previous MT, and here we look at for constant MT rates only. We found that:

At high ML rates, saturates. It becomes mainly a function of the core mass fraction and not of the ML rate (see Figure 10 where we show representative examples);

Nonsaturated values of are usually lower than the saturated .

at the very onset of RLOF is often hard to determine due to MLIPs.
While the radial response seems to be most natural to look at during the MT, we find that other properties of the donor’s surface also depend dramatically on how well is obtained. In Section 2.5 we showed that substantial luminosity can be generated if recombination energy is calculated properly. We can compare the surface luminosities of and giant, subjected to a constant ML yr and evolved with a default timestep adjustment and with a timestep chosen to resolve . An estimate provided in Section 2.5 shows that the energy coming from recombination in a giant subjected to ML is approximately
(22) 
Where does this energy mainly go – is it radiated away as excess luminosity or spent on the star itself, e.g., to increase its thermal or gravitational energy? The unperturbed giant has luminosity only about . The calculations show that a masslosing giant in which in the superadiabatic layer was calculated properly, has a luminosity about 6 times higher than the masslosing giant with a default relatively large timestep! The first giant has also become about 50% hotter. This increase in luminosity indicates that a large portion of energy from recombination is radiated away, rather than being spent on the star itself. In the case of a and giant, a simple prediction for a luminosity increase and an exact calculation would agree fairly well – in a more expanded giant, less recombination energy was spent on the star itself, and most went directly to surface luminosity. Similarly, the luminosities in the passya and passyb models are different by about 6 times, also as a simple estimate would predict. Interestingly, in the case of giants, the radii obtained by the two methods were not much different, except that in giants evolved with a small timestep we can observe MLIPs.
Of course, very fast MT rates are not commonly observed. However, the path of a star through the fast MT is different, and can define both the stability and/or how the donor appears when the fast MT phase is completed. We also can evaluate that this effect might be important (i.e. provides a difference in the luminosity by more than a few per cent) for as long as
(23) 
and hence might be important for lowluminosity giant donors and, in general, for all donors transferring mass at the thermal timescale of their envelope.
Indeed, consider the thermal timescale of the envelope taken as usual
(24) 
For most giants , and their envelope thermal timescale MT, , satisfies condition 23.
5.2 Stability of the Roche lobe mass transfer
The conventional way to determine stability of a binary system to the MT is to compare the initial responses of the donor radius and the Roche lobe radius to the mass loss. As was mentioned in Section 1, such analysis, if it does not involve actual MT calculations in a binary, is usually done by comparing their massradius exponents at the onset of RLOF. In order to do this, the adiabatic radial response of the donor is found from one of the existing approximations (composite polytropes, condensed polytropes, or even simply adopting that for convective donors ). This adiabatic radial response is used instead of the realistic stellar response.
However, we’d like to stress, that even if the realistic stellar response is obtained, there is no reason to assume that if at the start of the MT, and a donor’s relative Roche lobe overflow initially increases, the instability will necessarily occur. Instead of considering just the moment of RLOF, we can trace binary evolution through the MT in detail and determine if the instability eventually takes place, or not. The determination of MT instability in a binary system, while the donor is overflowing its RL, requires an alternative (to a simple comparison of the initial s) criterion to delineate when the MT proceeds stably, and when it results in a common envelope.
To understand the characteristic behaviour of the systems during MT, we performed a large set of simulations with donors of several initial masses (). We considered several values of radii for each donor, taken within the range where the donors have nonnegligible, outer convective envelopes, both at the red giant and asymptotic giant branches. We conducted simulations for several initial mass ratios (between 0.9 and 3.5). We assume that a companion is compact, where the compactness only means that we neglect any possible accretor’s RLOF.
Further, for these detailed MT sequences we adopted fully conservative MT; any nonconservation in the MT would lead to a relative increase of MT stability. Most of our binary systems were evolved using our modifications: obtaining proper , and our method to find MT rates and boundary conditions.
Even for systems with larger initial than that dictated by the conventional adiabatic stability criterion, the MT rate after RLOF smoothly increases, approaches a peak value and then decreases; the peak MT rate is usually much less than the dynamical MT rate. With the increase of , the peak MT rate increases. We define as the smallest initial mass ratio for which the binary experiences overflow during the MT (see Figure 12). We refer to this overflow as overflow, as during the MT, the mass ratio can reverse, and overflow becomes overflow. This is a function of both the initial mass ratio and the donor’s radius. For many initial donor radii, in binary systems with , the mass of the donor encompassed between the donor’s Roche lobe and its lobe is very small and the ML rate when overflow is approached is far from dynamical. Numerically, we are capable of evolving such systems with through the overflow, and many of them will not even approach the dynamical MT rate (see Figure 12). We, however, do not trust the MT rates obtained in this regime because our stream model only considers the nozzle.
We adopt therefore that for as long as overflow of the donor does not happen, the system is stable. In part, this is confirmed by the threedimensional hydrodynamical simulations, in which a common envelope is always associated with a severe, albeit very short in duration, overflow of the donor (Nandez et al., 2014). Note however, that we consider this to be a minimum stability criterion, as the MT in a binary system may remain stable even after the donor’s overflow takes place, as, e.g., SS 433 shows (e.g., Bowler, 2010; Perez M. & Blundell, 2010).
We also compared the surfaceaveraged effective gravitational acceleration over the lobe with the spherically symmetric acceleration that would be expected at its volume radius and found that the difference is about 13%.
In addition to overflow condition, we trace whether the binary orbit is changing rapidly. For the latter condition, we adopt that if , the orbit is not changing rapidly, where is the binary period. Rapid evolution of orbital parameters can invalidate the entire Roche formalism. For example, if angular velocity changes substantially over one period, the Euler term of the fictitious force, which is ignored in the Roche formalism along with the Coriolis term, becomes comparable to the centrifugal term. Note, however, that in Nature, binary systems with a larger might not experience dynamicaltimescale MT.
The key reason behind the unexpected stability of the MT in systems with is that once the MT starts, rises as the donor’s mass decreases (see Figures 10, 11 and 13). At the same time, decreases (Figure 13). A decrease in means that, as the MT proceeds, the stability of the system increases. We define as the critical masspoint, , the mass of the donor at which . When the donor decreases its mass to the critical masspoint, the ratio of the donor radius to its Roche lobe starts to decrease. Therefore, assuming that the relative RL overflow at the equipotential weakly depends on the mass ratio (see below), if there was no before the donor mass decreased to , the binary system is stable with respect to MT. We have searched for the point when overflow starts using our extended set of detailed simulations (see Figure 14). It can be seen that the critical mass ratios for which the MT becomes “unstable” are about twice as large as would be predicted by conventional comparison of massradius exponents at the onset of the MT.
Let us analyse the results and show how the approximate critical mass ratios can be predicted without doing detailed binary calculations.
To predict whether a system experiences overflow during the MT, we need to identify whether at any moment after the start of MT the donor’s radius exceeds – the volume radius at which the donor starts overflow. In order to do this, we can find the ratio for a range of the donor masses. Note that this critical ratio also depends on the mass ratio. However, we find that this dependence is quite weak, for . Nevertheless, we took this dependence into account in our calculations by examining the ratio not only down to donor mass , but also further, almost all the way down to the donor core mass.
The radius of the donor during the MT, when it has shed to mass , is
(25) 
It follows from Section 5.1 that in order to predict the radius of the donor at any moment during fast (saturated) MT, instead of the real one can use , which is obtained from the composite polytrope approximation. “Compact” giants, should be expected to expand less than this estimate would predict, and hence be more stable because their .
Using the definitions of massradius exponents (see also Figure 13), we can find the approximate radius of the donor at any moment after the start of MT as
(26) 
Now we can find whether the system experiences overflow at any moment during the MT and thus produce a semianalytical estimate of – we will refer to this estimate as “overflow condensed polytrope simplification”. It is plotted in Figure 14 for various giants.
We can see that the condensed polytrope simplification works best for those giants that are neither too compact nor too expanded. In both compact and averagesized giants that approach overflow at the critical point, most mass before the critical point is lost at ML rates which are generally much higher than , so it’s fine to use the saturated value of . At the same time, compact giants have higher saturated values of , which makes them more stable than averagesized giants whose saturated approximately follows the condensed polytrope.
In large giants with very rarefied envelopes, overflow occurs at much lower mass loss rates, and the saturated values of become inapplicable. Instead, a lower, nonadiabatic, nonsaturated values should be used. This makes these giants less stable. However, such overflow will be nondynamical and does not have to result in a common envelope.
In Figure 15 we show how the critical mass ratio obtained from the overflow condensed polytrope simplification changes with the initial donor mass and radius. We have also checked and found that nonconservative MT leads, as expected, to the increase in , increasing it by about 0.3.
The detailed simulations (Figure 14) show that is larger when a donor is more compact and the mass fraction of its convective envelope is smaller. In all detailed simulations we did with donors where the outer convective envelope is nonnegligible, is below 3.5, where 3.5 is known as the critical mass ratio leading to a delayed dynamical instability with radiative donors (see, e.g, Ge et al., 2010). It shows therefore that during the development of the outer convective envelope, the convective donors are in the transitional regime. When the convective envelope is well developed, the critical mass ratios are in the range for most donors.
We conclude that the criterion based on overflow is definitely predicting that more binary systems will evolve through their MT in a stable way.
Whether the overflow should necessarily lead to a common envelope is not entirely clear. In principle, MT through the nozzle can be treated in the same simplified way as that through the nozzle, and the material can be considered to leave the system carrying away the specific angular momentum of the donor. It is a subject of our future work.
6 Conclusions
In this paper, we considered a number of theoretical and practical challenges on the path to understanding the stability of the MT in binaries consisting of a convective giant donor and a compact accretor.
We find that in order to obtain the correct response to the mass loss, the in the superadiabatic layer must be calculated properly; depending on which stellar code is used this may require a number of numerical tricks, and at the very least control of the timestep to match the predicted energy generation rate. We provide simple estimates for how to find the mass of this superadiabatic layer and the rate of the energy release, in order to quantify how well an arbitrary stellar code performs under mass loss.
We have shown that the massradius exponents in the giants that are compatible with the compositepolytrope description converge to the prediction in Hjellming & Webbink (1987) for fast mass loss. On the other hand, the giants that cannot be described by a composite polytrope have , and the binary systems with such donors are more stable than the composite polytrope would predict. In our detailed models, we could not find the strong superadiabatic expansion predicted by the adiabatic models in Ge et al. (2010).
In addition to radial response, we examined the changes in surface luminosity and effective temperature of masslosing stars. We find that at fast ML rates, the luminosity can differ by a factor of several times between the models where in the superadiabatic layer is calculated correctly and the ones where it is calculated in the Lagrangian way with large timesteps. We note that in some observed binary systems, the donor’s properties such as effective temperature and mass are determined with a few % precision (e.g., Leahy & Abdallah, 2014). This difference in effective temperature obtained by the two models, and comparable to the observed precision, is expected to occur in systems with a lowluminosity giant donor, and also in thermaltimescale MT systems.
We have enhanced the classical scheme to find the MT rate via an optically thick stream approximation. In particular, our use of the detailed system geometry allows us to find the MT rate in the case of a substantial RLOF. This leads us to a new criterion of MT instability, based on whether the donor starts overflow. We have also found that a binary system does not become immediately unstable as the donor’s envelope becomes convective, but rather the mass ratio at which the instability occurs gradually decreases from the regime predicted for radiative donors. Note that in principle a binary can survive the MT even after overflow, without starting a common envelope. However, we find that even the new overflow criterion warrants that binary systems will proceed with stable conservative MT if their mass ratio is up to twice that given by the conventional criterion, for stars with deep convective zones (and the mass ratio can be even larger for shallow outer convective zones).
Acknowledgements
K.P. acknowledges support from Golden Bell Jar fellowship. N.I. acknowledges support by NSERC Discovery Grants and the Canada Research Chairs Program; this research was supported in part by the National Science Foundation under Grant No. NSF PHY0551164. Authors thank B. Paxton for help with MESA/star and MESA/eos modules, T. Fragos for help with MESA/binary module, C. Heinke for help with the manuscript and the anonymous referee for useful comments. Authors are also grateful to R. Webbink, Ph. Podsiadlowski and S. Justham for useful discussions. This research has been enabled by the use of computing resources provided by WestGrid and Compute/Calcul Canada.
References
 Baudin et al. (2011) Baudin F., Barban C., Belkacem K., Hekker S., Morel T., Samadi R., Benomar O., Goupil M.J., Carrier F., Ballot J., Deheuvels S., De Ridder J., Hatzes A. P., Kallinger T., Weiss W. W., 2011, A&A, 529, A84
 Belkacem et al. (2012) Belkacem K., Dupret M. A., Baudin F., Appourchaux T., Marques J. P., Samadi R., 2012, A&A, 540, L7
 Bowler (2010) Bowler M. G., 2010, A&A, 521, A81
 Cox & Giuli (1968) Cox J., Giuli T., 1968, Principles of stellar structure. GORDON AND BREACH, Science Publishers, Inc.
 Eggleton (1967) Eggleton P. P., 1967, MNRAS, 135, 243
 Eggleton (1971) Eggleton P. P., 1971, MNRAS, 151, 351
 Eggleton (1972) Eggleton P. P., 1972, MNRAS, 156, 361
 Eggleton (1973) Eggleton P. P., 1973, MNRAS, 163, 279
 Eggleton (1983) Eggleton P. P., 1983, ApJ, 268, 368
 Eggleton et al. (1973) Eggleton P. P., Faulkner J., Flannery B. P., 1973, A&A, 23, 325
 Ferguson et al. (2005) Ferguson J. W., Alexander D. R., Allard F., Barman T., Bodnarik J. G., Hauschildt P. H., HeffnerWong A., Tamanai A., 2005, ApJ, 623, 585
 Ge et al. (2010) Ge H., Hjellming M. S., Webbink R. F., Chen X., Han Z., 2010, ApJ, 717, 724
 Glebbeek et al. (2008) Glebbeek E., Pols O. R., Hurley J. R., 2008, A&A, 488, 1007
 Han et al. (2002) Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R., Ivanova N., 2002, MNRAS, 336, 449
 Hjellming & Webbink (1987) Hjellming M. S., Webbink R. F., 1987, ApJ, 318, 794
 Ivanova et al. (2013) Ivanova N., Justham S., Chen X., De Marco O., Fryer C. L., Gaburov E., Ge H., Glebbeek E., Han Z., Li X.D., Lu G., Marsh T., Podsiadlowski P., Potter A., Soker N., Taam R., Tauris T. M., van den Heuvel E. P. J., Webbink R. F., 2013, A&Ar, 21, 59
 Ivanova & Taam (2004) Ivanova N., Taam R. E., 2004, ApJ, 601, 1058
 Kevorkian & Cole (1985) Kevorkian J., Cole J. D., 1985. Springer
 Kippenhahn & Weigert (1994) Kippenhahn R., Weigert A., 1994, Stellar Structure and Evolution. SpringerVerlag Berlin Heidelberg New York
 Kippenhahn et al. (1967) Kippenhahn R., Weigert A., Hofmeister E., 1967, Methods Comput. Phys., 7, 129
 Kolb & Ritter (1990) Kolb U., Ritter H., 1990, A&A, 236, 385
 Kuhfuss (1986) Kuhfuss R., 1986, A&A, 160, 116
 Leahy & Abdallah (2014) Leahy D. A., Abdallah M. H., 2014, ApJ, 793, 79
 Lubow & Shu (1976) Lubow S. H., Shu F. H., 1976, ApJ, 207, L53
 Nandez et al. (2014) Nandez J. L. A., Ivanova N., Lombardi Jr. J. C., 2014, ApJ, 786, 39
 Nelemans & Tout (2005) Nelemans G., Tout C. A., 2005, MNRAS, 356, 753
 Nelemans et al. (2000) Nelemans G., Verbunt F., Yungelson L. R., Portegies Zwart S. F., 2000, A&A, 360, 1011
 Nelson & Eggleton (2001) Nelson C. A., Eggleton P. P., 2001, ApJ, 552, 664
 Osaki (1970) Osaki Y., 1970, ApJ, 162, 621
 Paczynski (1965) Paczynski B., 1965, AcA, 15, 89
 Paczyński (1971) Paczyński B., 1971, ARA&A, 9, 183
 Paczyński & Sienkiewicz (1972) Paczyński B., Sienkiewicz R., 1972, Acta Astronomica, 22, 73
 Passy et al. (2012) Passy J.C., Herwig F., Paxton B., 2012, ApJ, 760, 90
 Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJs, 192, 3
 Paxton et al. (2013) Paxton B., Cantiello M., Arras P., Bildsten L., Brown E. F., Dotter A., Mankovich C., Montgomery M. H., Stello D., Timmes F. X., Townsend R., 2013, ApJs, 208, 4
 Perez M. & Blundell (2010) Perez M. S., Blundell K. M., 2010, MNRAS, 408, 2
 Ritter (1988) Ritter H., 1988, A&A, 202, 93
 Savonije (1978) Savonije G. J., 1978, A&A, 62, 317
 Soberman et al. (1997) Soberman G. E., Phinney E. S., van den Heuvel E. P. J., 1997, A&A, 327, 620
 Webbink (1985a) Webbink R. F., 1985a, Stellar evolution and binaries
 Webbink (1985b) Webbink R. F., 1985b, Stellar evolution and binaries. p. 39
 Woods & Ivanova (2011) Woods T. E., Ivanova N., 2011, ApJ, 739, L48
 Woods et al. (2012) Woods T. E., Ivanova N., van der Sluys M. V., Chaichenets S., 2012, ApJ, 744, 12