Mass Transfer from Giant Donors
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 15 September 2018 (MN L A TEX style file v2.2)
Mass transfer from giant donors
K. Pavlovskii, N. Ivanova, University of Alberta, Dept. of Physics, 11322-89 Ave, Edmonton, AB, T6G 2E7, Canada
15 September 2018
ABSTRACT
The stability of mass transfer in binaries with convective giant donors remains anopen question in modern astrophysics. There is a significant discrepancy betweenwhat 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 showthat the recombination energy in the superadiabatic layer plays an important andhitherto unaccounted-for role in the donor’s response to mass loss, in particular onits luminosity and effective temperature. Our improved optically thick nozzle methodto calculate the mass transfer rate via L allows us to evolve binary systems for asubstantial Roche lobe overflow. We propose a new, strengthened criterion for themass transfer instability, basing it on whether the donor experiences overflow throughits outer Lagrangian point. We find that with the new criterion, if the donor has awell-developed outer convective envelope, the critical initial mass ratio for which abinary would evolve stably through the conservative mass transfer varies from 1 . .
2, which is about twice as large as previously believed. In underdeveloped giantswith shallow convective envelopes this critical ratio may be even larger. When theconvective 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.
Key words: instabilities — methods: numerical — stars: mass-loss — binaries: close— stars: evolution
Many interacting binaries start their mass exchange whenthe donor, which has evolved to the giant branch, overfillsits Roche lobe. The stability at the start of the mass trans-fer (MT) differentiates between the binaries that will livelong as a mass transferring system and will appear, e.g. ,as an X-ray binary, and ones that will be transformed dra-matically by a common envelope event. A clear separationof the possible evolutionary channels is important for ourunderstanding of the formation of compact binaries. E.g.,double white-dwarf binaries were thought in the past to beformed via two common envelope events (Nelemans et al.2000; Nelemans & Tout 2005), while more recent detailedsimulations of MT with low-mass red giant donors show thatthe 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 sig-nificantly with time, in some cases leading to so-called de-layed dynamical instability (for a thorough discussion, seeGe et al. 2010). A delayed rapid growth of MT rate can alsooccur in systems with giant donors(Woods et al. 2012). Un-fortunately, to date there exists no simple criterion for thestart of a dynamically unstable MT that would be based juston its rate. It can be defined as MT happening on the dy- namical timescale of the donor, but even this definition doesnot determine the fate of a system (for a review of variousproblems with the initiation of a common envelope event,see Ivanova et al. 2013).Standard definition of a dynamical-timescale MT is thatthe mass is transferred on the donor’s dynamical timescale(see, e.g. , Paczy´nski & Sienkiewicz 1972). When a star issubjected to a dynamical-timescale mass loss (ML), stellarcodes operate outside of the regime in which they were de-signed to operate, leading to poorly understood behaviour.Indeed, a conventional stellar code, developed to treat evolu-tionary changes, typically does not include hydrodynamicalterms in structure equations, and hence is unable to treatdynamical-timescale features of the ML response.There are then two ways to deal with a dynamical MLin practice, both of them relying on plausible-sounding butarbitrary assumptions.One way is to pre-set a threshold for the MT rate closeto (but still within) the limits of validity of a conventionalstellar code, and assume MT to become unstable once thatthreshold is exceeded. For example, a dynamical MT hasbeen defined as 10 times the thermal timescale MT (Nelson& Eggleton 2001).Another way is to assume that during a very short pe-riod of very rapid mass loss (ML), any mass shell in the c (cid:13) a r X i v : . [ a s t r o - ph . S R ] M a r Pavlovskii & Ivanova stellar model has no time to exchange heat with its neigh-bouring shells. Hence, its entropy remains almost the sameas at the start of the ML, in the other words the specificentropy profile of the donor is “frozen”.The MT rate and stability depend crucially on the re-sponse of the donor’s radius to the ML compared to theresponse of its Roche lobe radius to the same ML (for foun-dations, see Webbink 1985a). The governing feature thatshaped studies of MT in binary systems with giant donorswas the wide acceptance of the theory that a convectivedonor would intrinsically expand as a reaction to ML. Inthe framework of this theory, the immediate response of astar to ML takes place on a dynamical timescale (to restoreits hydrostatic equilibrium) – this assumes that no thermaleffects can take place, and it was hence called an adiabaticresponse. It is characterised by the quantity known as themass-radius exponent ζ ad ζ ad ≡ (cid:18) ∂ log R∂ log M (cid:19) ad . (1)With polytropic models (including also composite andcondensed polytropes), it was found that a convective donorwould have ζ ad < ζ ad , which becomes positive once the relative mass ofthe core is more than 20% of the star’s total mass (Hjellming& Webbink 1987).It is commonly assumed that MT is dynamically sta-ble only if a donor remains within its Roche lobe, or when ζ L ≤ ζ ad , where ζ L is the mass-radius exponent for the Rochelobe reaction to the MT(Webbink 1985b). Using the Rochelobe approximation as in Eggleton (1983), for conservativeMT, ζ L = 0 when the mass ratio is q crit ≈ . ≈ . < ∼ .
45 of the total stellarmass, the first episode of conservative MT will be dynamicallunstable. Later, the results based on polytropic models werere-evaluated in the studies of the detailed adiabatic stellarmodels by Ge et al. (2010). At the same time, detailed codesthat traced mass transfer in binaries up to thermal timescalemass transfer, found that critical mass ratio can be up to1.1 ( e.g. , Han et al. 2002).It is crucial to realize that the difference in rates be-tween a nearly-thermal timescale MT that can be legiti-mately calculated with a conventional stellar code and adynamical-timescale MT – the regime of adiabatic codes –is several orders of magnitude, and currently there is no ad-equate treatment to describe the donor response in betweenthese two regimes.A detailed stellar code can be supplied with hydrody-namical terms in its structure equations and be forced towork in a hydrodynamical regime (please note that adiabaticcodes do not include hydrodynamical terms in them). Thisway, in the recent study of responses of stellar models of gi-ants to fast ML, it was found that this response is a functionof the MT rate (Woods & Ivanova 2011). As a result, it wasfound that for a large range of ML rates taking place in bina- ries, giants can also contract, and the nature of the donor’sresponse was attributed to the behaviour of the giant’s su-peradiabatic surface layer and its short thermal timescale.In that study, two different stellar codes were used for com-parison and they showed the same results. The two codeswere: a) STARS/ev code developed by Eggleton (Eggleton(1971, 1972, 1973); Eggleton et al. (1973)) with the mostrecent update described as in Glebbeek et al. (2008); b) theHeyney-type 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 codesappeared to produce the same result – the donor responseis very different from the expectations given by a simplifiedadiabatic model. However, no further analysis of why thereis a dramatic difference between an adiabatic approach anda detailed stellar code was provided. The important questionremains: which result we should trust more? What could bethe possible problems with the adiabatic approach? Couldwe trust the results of the detailed stellar codes when theyare forced to work at dynamical-timescale MT rates, andwhen they possibly break down?It is clear that for further progress in understanding ofthe response of the giant’s radius to fast ML, one needs tounderstand better the non-adiabatic processes taking placein a red giant subjected to dynamical-timescale ML. This isintrinsically related to the nature – the underlying physicsand structure – of the superadiabatic layer, which is not yetwell understood. A possible influence of the surface layer ingiants on MT was realised in the past. For example, Osaki(1970) had suggested that, in U Gem giants, the mechanismfor dynamical mass exchange, proposed by Paczynski (1965),can not work, “because of changes occurring on the thermaltime scales of the photosphere and transition zone, which areshorter than the dynamical time scale of the star”). How-ever, the physics of the superadiabatic layer had never beenstudied in detail. On the other hand, the difference in theresponse in different codes also could be related to numericaleffects, which could be either inconsistency in the adoptedset of equations, or numerical artefacts.In this paper, in §
2, we systematically analyse thephysics of the superadiabatic layer: what difference an adia-batic approach makes, which artefacts might appear if thislayer is not numerically modelled properly, what affects thestructure of this layer, how massive this layer is in differentdonors, and finally, whether we can numerically obtain theresponse of the stellar radius properly with a detailed stellarcode, and with what limitations.All detailed stellar models and numerical experimentsin this paper are calculated with
MESA (Modules for Experi-ments 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 easeof use, it is becoming widely popular, and may soon becomethe most commonly used stellar tool. For the purpose ofcomparison, we clarify that, by default, we use solar metallic-ity, and employ a standard set of stellar equations: the sim-ple grey atmosphere boundary condition with radiation pres-sure at zero optical depth taken from Cox & Giuli (1968),mixing-length theory (MLT) as in Cox & Giuli (1968), andradiative opacity tables based on OPAL and Ferguson et al. c (cid:13) , 000–000 ass transfer from giant donors (2005). If there is any deviation from this default set, it willbe clarified for each numerical experiment. While some nu-merical issues that we will discuss in this paper are relevantonly to our choice of the stellar code, most of the issues arerelevant to modeling of fast ML with any detailed stellarcode. We are aware that MESA was developed to evolve nor-mal stars, in a conventional way, and in this paper we willanalyse how
MESA , and possibly similar codes, work at anextreme beyond its design.While the response of the donor’s radius itself to MLis the first important issue, the second foremost is how faststellar material, overfilling the donor’s Roche lobe, can betransferred to its companion, and hence be effectively lostfrom the donor. In § § § A major feature distinguishing adiabatic calculations of theMT from the results of detailed stellar codes for convectivedonors is the retention of a thin (in mass sense) layer oflowered entropy near the surface. Understanding the physicsof this layer is likely to be key to understanding a realisticstellar response to ML.
The term “superadiabatic” makes sense only within theadopted convection theory, MLT.In MLT, the relative efficiency of heat transfer via twomechanisms (radiative and convective) can be characterizedby the ratio of radiative and convective conductivites (seeequation 7.12 in Kippenhahn & Weigert 1994): σ rad σ conv = 3 acT c P ρ κl (cid:114) H P gδ , (2)where l m is the mixing length and H P is pressure scaleheight. δ = − ( ∂ ln ρ/∂ ln T ) P is a quantity that charac-terises the thermal expansion properties of a medium.In the framework of MLT, convective efficiency directlydepends on two free parameters – on how far a blob rises orsinks before it mixes with the environment (which can beseen in Equation 2 directly as l m ), and on the linear size ofconvective blobs (used in Equation 2). These free parametersare best calibrated from observations and are not obtainedfrom physical considerations.In convective zones overadiabaticity, which is the dif-ference between the actual and adiabatic temperature gra-dients, ∇ − ∇ ad , is positive. We define the superadiabaticlayer as the region where overadiabacity is comparable tothe gradients themselves.In this paper, for the efficiency of convection, we adoptthe relation between conductivities of convective and ra-diative heat transfer denoted by σ conv /σ rad = 1 /U . Largervalues of U correspond to less efficient convection and ac-cordingly to a larger over-adiabaticity. In red giants, due to -5 -5 -5 -5 tot - m) [M ☉ ] Figure 1.
Interplay of temperature gradients ∇ , ∇ rad and ∇ ad inthe vicinity of the hydrogen ionisation boundary of a 70 R (cid:12) , 5 M (cid:12) red giant. The radiative gradient is shown with a solid grey line,adiabatic gradient – solid black line, actual gradient – dashed line; m is the mass coordinate and M tot is the total giant’s mass. Nearthe surface, where ∇ rad < ∇ ad , this giant has a surface radiativezone. A superadiabatic zone is where ∇ rad > ∇ ad but convectionis not efficient enough, and the actual gradient ∇ exceeds ∇ ad bya value comparable to the gradient itself. a combination of factors U increases in the vicinity of thestellar surface to the extent that the real gradient of tem-perature becomes substantially detached from the adiabaticgradient ( ∇ > ∇ ad ) (see Figure 1). Let us consider in detail what affects the behaviour of theefficiency of the convection 1 /U and of the overadiabaticity ∇−∇ ad . In the region where both quantities have simultane-ous slow-down in their almost monotonic decrease towardsthe center (see Figure 2), hydrogen and helium change theirdegree of ionisation (see Figure 3). The degrees of ionisationaffect U and ∇ − ∇ ad by strongly changing the heat capac-ity of the medium c P and the Rosseland mean opacity (seeFigure 4).For the range of densities and temperatures in the en-velope, the heat capacity of a partially ionised gas is higherthan that of both neutral and fully ionised gases (ionisationconsumes a part of the heat transferred to the gas withoutincreasing the kinetic energy of its molecules). The threemajor jumps in c P , from highest to lowest temperature, areassociated with the following recombinations:(i) double ionised helium → single ionised helium;(ii) single ionised helium → neutral helium;(iii) ionised hydrogen → neutral hydrogen.Let us discuss the reason for two overadiabaticityplateaus in the area of interest marked with the solid opac-ity curve in Figure 4. In this area, partial recombination ofelements affects the convective efficiency mainly through c P ,because opacity does not vary much there. Due to the changein monotonicity of c P at half-ionisation, recombination ofeither hydrogen or helium impedes the growth of overadi-abaticity towards the surface while the ionisation fractionis between 100% and 50%, but when the ionisation fractionis lower than 50% the recombination, on the other hand, c (cid:13)000
Interplay of temperature gradients ∇ , ∇ rad and ∇ ad inthe vicinity of the hydrogen ionisation boundary of a 70 R (cid:12) , 5 M (cid:12) red giant. The radiative gradient is shown with a solid grey line,adiabatic gradient – solid black line, actual gradient – dashed line; m is the mass coordinate and M tot is the total giant’s mass. Nearthe surface, where ∇ rad < ∇ ad , this giant has a surface radiativezone. A superadiabatic zone is where ∇ rad > ∇ ad but convectionis not efficient enough, and the actual gradient ∇ exceeds ∇ ad bya value comparable to the gradient itself. a combination of factors U increases in the vicinity of thestellar surface to the extent that the real gradient of tem-perature becomes substantially detached from the adiabaticgradient ( ∇ > ∇ ad ) (see Figure 1). Let us consider in detail what affects the behaviour of theefficiency of the convection 1 /U and of the overadiabaticity ∇−∇ ad . In the region where both quantities have simultane-ous slow-down in their almost monotonic decrease towardsthe center (see Figure 2), hydrogen and helium change theirdegree of ionisation (see Figure 3). The degrees of ionisationaffect U and ∇ − ∇ ad by strongly changing the heat capac-ity of the medium c P and the Rosseland mean opacity (seeFigure 4).For the range of densities and temperatures in the en-velope, the heat capacity of a partially ionised gas is higherthan that of both neutral and fully ionised gases (ionisationconsumes a part of the heat transferred to the gas withoutincreasing the kinetic energy of its molecules). The threemajor jumps in c P , from highest to lowest temperature, areassociated with the following recombinations:(i) double ionised helium → single ionised helium;(ii) single ionised helium → neutral helium;(iii) ionised hydrogen → neutral hydrogen.Let us discuss the reason for two overadiabaticityplateaus in the area of interest marked with the solid opac-ity curve in Figure 4. In this area, partial recombination ofelements affects the convective efficiency mainly through c P ,because opacity does not vary much there. Due to the changein monotonicity of c P at half-ionisation, recombination ofeither hydrogen or helium impedes the growth of overadi-abaticity towards the surface while the ionisation fractionis between 100% and 50%, but when the ionisation fractionis lower than 50% the recombination, on the other hand, c (cid:13)000 , 000–000 Pavlovskii & Ivanova U x ( ∇ - ∇ a d ) x log (T [K]) Figure 2.
Overadiabaticity ∇ − ∇ ad (solid) and U (dashed) inthe outer layers of a 5 M (cid:12) and 50 R (cid:12) red giant. ( ∇ - ∇ a d ) x A v g . c h a r g e [ e ] log (T [K]) Figure 3.
Overadiabaticity ∇ − ∇ ad (solid) and average chargesof helium (dotted) and hydrogen (dashed) in the outer layers ofa 5 M (cid:12) and 50 R (cid:12) red giant. Note that behaviour of overadia-baticity and partial ionisation zones is strongly coupled. accelerates the growth of overadiabaticity (see Equation 2and Figure 4). Because of this, the areas where the growthof overadiabaticity is impeded (plateaus) are shifted withrespect to the maxima of c P . The initial stage of partial re-combination of hydrogen is able to completely suppress thegrowth of overadiabaticity towards stellar surface in the en-velope; the initial recombination of helium also impedes it toa large extent. The shape of the superadiabatic layer, hence,is substantially affected by the recombination of hydrogenand helium.At the final stages of hydrogen recombination, the opac-ity falls sharply by a few orders of magnitude. As a result,the radiative gradient, which is proportional to opacity, fallsbelow the adiabatic gradient (Figure 1). We call this regionthe ”surface radiative zone”. It takes place outside the tem-perature range shown in Figures 2-4, at lg T < .
73. Thespatial size of the surface radiative zone is of the order of themixing length. Hence it is plausible that convective blobsovershoot the boundary between the superadiabatic layerand the surface radiative zone (see also Kuhfuss 1986). ( ∇ - ∇ a d ) x ; l o g ( κ [ c m g - ]) l o g ( C p [ e r g K - g - ]) log (T [K]) Figure 4.
Overadiabaticity ∇ − ∇ ad (solid black), Rosselandmean opacity κ (grey solid in the area of interest and grey dashedelsewhere) and specific heat capacity at constant pressure c P (dashed black) in the outer layers of a 5 M (cid:12) and 50 R (cid:12) redgiant.
20 25 30 35 40 3.5 4 4.5 5 5.5 6 S p e c i fi c e n t r o p y [ k / N A ] log (T [K]) Figure 5.
Specific entropy in the outer layers of a 14 . M (cid:12) and500 R (cid:12) red giant, the cross denotes the bottom of the superadi-abatic layer calculated as per Equation 3. For analysis of superadiabatic layers and radiative zones, weintroduce the quantity m sad , the mass of these layers. Wemeasure m sad from the stellar surface to the point near thestart of the entropy drop, m . Since the convection is neverfully adiabatic, the entropy profile in the envelope is not flat,and there is no sharp transition between the entropy plateauand the entropy drop (see Figure 5). For integrity betweenstars of different masses, we therefore choose to define m as the point where entropy is S ( m ) = S min + 0 . S conv − S min ) . (3)Here S min is the local minimum of entropy in the outer layerof the star and S conv is the (maximum) value of entropy inthe convective envelope.Clearly, it can be seen in Figure 6 that the mass ofthis shell varies between giants of different masses and radii.From the behaviour of the envelope’s simplified analyti-cal solution (see e.g. Figure 10.2 in Kippenhahn & Weigert1994), we expect that overadiabaticity of convective en-velopes becomes substantial at about the same pressure, c (cid:13) , 000–000 ass transfer from giant donors -10-5 0 5 0 4 8 12 l o g ( M s a d [ M ☉ ]) log (R /M tot [R ☉ M ☉ -1 ])1 M ☉ ☉ ☉
10 M ☉
15 M ☉
20 M ☉ Figure 6.
Total mass of superadiabatic layers and surface radia-tive zones in stars with convective envelopes of different massesand ages. to an order of magnitude. Indeed, in realistic stellar mod-els for stars of different masses, the pressure at the bot-tom of the superadiabatic layer does not vary much from P ( m ) ≈ dyn/cm . In this case, we can estimate m sad by simply considering a pressure drop. The equation of hy-drostatic equilibrium over the superadiabatic layer and ra-diative zone is dPdm sad = G ( m + m sad ( r ))4 πr . (4)As the pressure at the top of the superadiabatic layer (pho-tospheric pressure) is much less than at its bottom, we find m sad ( R ) ≈ πR P GM tot ≈ . × − M (cid:12) ( R/R (cid:12) ) ( M tot /M (cid:12) ) . (5)Indeed, as illustrated in Figure 6, stars with convective en-velopes at any evolutionary stage have m sad which can befound using this relation. This approximation has only lim-ited, though still very good, applicability to more massivestars, where m sad is comparable to the total mass of thestar’s envelope (we note also that it is well known that inmassive stars almost the entire convective envelope could besuperadiabatic).It is important to note that this is the mass of the su-peradiabatic layer and surface radiative zone in a star un-perturbed by ML. Mass loss affects both our assumptionsabout the pressure drop (as the star is not in hydrostaticequilibrium anymore), and the entropy profile.We see from the preceding analysis that while the massof the superadiabatic layer is relatively small in low-massgiants, it grows for more massive stars, increasing substan-tially as the star is evolving and expanding. Stellar modelsthat have a constant entropy profile are therefore not jus-tified, as the error in the case of massive stars could beenormous. We will also show in § We saw a hint in § -1 0 1 2 3 -6-4.5 ε g [ e r g · g - · s - ] log (M-m) [M ☉ ] Figure 7.
Values of (cid:15) g in the same 5 . M (cid:12) , 50 R (cid:12) red giant sub-jected to the same constant ML rate of 10 − M (cid:12) yr − shown atthe same moment of time: 0 . yr since the start of ML. Dash-dotted, dotted, dashed and solid lines correspond to constanttimesteps 10 − , 10 − , 10 − yr and 10 − yr respectively. Theobserved discrepancy in the values of (cid:15) g is due to the inaccuracyof the first order formula used to calculate Lagrangian derivativesusing the Equation 6 (see Section 2.4 for details). hence this layer might be energetically important for obtain-ing the correct ML response. Energetically, recombinationenters the structure equations through the gravitational en-ergy term, which is the last term in the energy generationequation and is defined as: (cid:15) g = − T (cid:18) ∂s∂t (cid:19) m (6)The (cid:15) g term can be found directly from the Equation 6,using only one temporal derivative of the entropy which isderived from the EOS. Alternatively, (cid:15) g can be calculated bysubtracting two terms that depend on temporal derivativesof thermodynamic quantities, for example: (cid:15) g = − T c P (cid:20) (1 − ∇ ad χ T ) 1 T (cid:18) ∂T∂t (cid:19) m − ∇ ad χ ρ ρ (cid:18) ∂ρ∂t (cid:19) m (cid:21) (7)When mass is removed from the surface of a red gi-ant, the partial ionisation zones of hydrogen and heliumdiscussed in § W ∼ erg g − . For the ML rate of˙ M = 10 − M (cid:12) yr − the recombination energy release shouldbe of the order of W ˙ M ∼ × erg s − .Let us consider now this mass loss on an example ofa 5 . (cid:12) and 50 R (cid:12) red giant. Using our Equation 5 forthe mass of the superadiabatic layer, we find that its totalmass is about 10 − M (cid:12) , though the part where hydrogenrecombination energy is released is usually only a fractionof it, about 10 − M (cid:12) (see also Figure 7). The contribu-tion of this recombination component alone in (cid:15) g in sucha superadiabatic layer is expected to be of the order of3 × erg g − s − .We have checked what values of (cid:15) g we obtain in practice ,using the detailed stellar code of our choice, MESA . First, wetested the method where the entropy derivative is found us- c (cid:13)000
Values of (cid:15) g in the same 5 . M (cid:12) , 50 R (cid:12) red giant sub-jected to the same constant ML rate of 10 − M (cid:12) yr − shown atthe same moment of time: 0 . yr since the start of ML. Dash-dotted, dotted, dashed and solid lines correspond to constanttimesteps 10 − , 10 − , 10 − yr and 10 − yr respectively. Theobserved discrepancy in the values of (cid:15) g is due to the inaccuracyof the first order formula used to calculate Lagrangian derivativesusing the Equation 6 (see Section 2.4 for details). hence this layer might be energetically important for obtain-ing the correct ML response. Energetically, recombinationenters the structure equations through the gravitational en-ergy term, which is the last term in the energy generationequation and is defined as: (cid:15) g = − T (cid:18) ∂s∂t (cid:19) m (6)The (cid:15) g term can be found directly from the Equation 6,using only one temporal derivative of the entropy which isderived from the EOS. Alternatively, (cid:15) g can be calculated bysubtracting two terms that depend on temporal derivativesof thermodynamic quantities, for example: (cid:15) g = − T c P (cid:20) (1 − ∇ ad χ T ) 1 T (cid:18) ∂T∂t (cid:19) m − ∇ ad χ ρ ρ (cid:18) ∂ρ∂t (cid:19) m (cid:21) (7)When mass is removed from the surface of a red gi-ant, the partial ionisation zones of hydrogen and heliumdiscussed in § W ∼ erg g − . For the ML rate of˙ M = 10 − M (cid:12) yr − the recombination energy release shouldbe of the order of W ˙ M ∼ × erg s − .Let us consider now this mass loss on an example ofa 5 . (cid:12) and 50 R (cid:12) red giant. Using our Equation 5 forthe mass of the superadiabatic layer, we find that its totalmass is about 10 − M (cid:12) , though the part where hydrogenrecombination energy is released is usually only a fractionof it, about 10 − M (cid:12) (see also Figure 7). The contribu-tion of this recombination component alone in (cid:15) g in sucha superadiabatic layer is expected to be of the order of3 × erg g − s − .We have checked what values of (cid:15) g we obtain in practice ,using the detailed stellar code of our choice, MESA . First, wetested the method where the entropy derivative is found us- c (cid:13)000 , 000–000 Pavlovskii & Ivanova ing 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 entropyfrom the EOS. After performing numerical tests, we foundit to not be suitable for the small timesteps necessary forfast mass loss calculations: as the timesteps become smaller,the errors in the entropy values become comparable to theirLagrangian differences between subsequent timesteps. Thetwo-component 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 thetabulated EOS. In addition, composition artefacts can playa role in entropy fluctuations.Second, we tested the method where the entropy deriva-tive is found using Equation 7, see Figure 7. We found thatthe expected physically reasonable level of (cid:15) g is only reachedwhen the timestep is below ≈ − yr. We relate it to theinaccuracy of the first-order numerical differentiation for-mulae used in practice to calculate the right-hand side of(7) at larger timesteps. Note that the right-hand side of (7)is effectively proportional to the Lagrangian derivative ofentropy , which is calculated indirectly through the corre-sponding derivatives of T and ρ . Therefore, the error in cal-culation of (cid:15) g with formula (7) is, as expected, the mostsignificant in the superadiabatic layer and surface radiativezone, because the second Lagrangian derivative of entropyis the highest there .We do not suggest a complete formal procedure to cal-culate the errors in calculation of (cid:15) g . Instead, we resort to aplain comparison of results obtained with various timestepselection approaches. For example, the comparison for oneof our models (Figure 7) shows that if it loses not morethan ≈ / m sad in one timestep, thenthe calculation of (cid:15) g in the superadiabatic layer is affectedonly marginally. Removing mass in such small timesteps is,however, quite resource intensive. We can foresee that theoutlined numerical problem with a recombination zone thatmoves too fast in mass can be eventually resolved using atechnique similar to the one used for calculating thin shellburning as in Eggleton (1967, 1973). However, for the pur-pose of the studies presented in this paper, we can afford touse a numerically intensive way.In Section 5.1 we will show how inaccurate calculationsof (cid:15) g affect the radial response to the mass loss. While thenumerical problems that we described in this section arenative only to the code we use, when fast mass loss rate cal-culations are performed with another code, the numericallyobtained (cid:15) g has to be tested against its expected value thatcan be found using the method described in this subsection. If one neglects composition changes. We’d also like to mention that quite recently a new schemefor the calculation of (cid:15) g was implemented in MESA , that followsa mixed Lagrangian-Eulerian approach. We performed a prelim-inary testing of this scheme on mass-losing red giants. Unfor-tunately, this new scheme is not suitable for our simulations yetbecause of severe numerical artefacts that it introduces under cer-tain conditions. These artefacts have to be carefully studied andeliminated before we can consider incorporating this new schemeinto our simulations. Because of this, we only discuss the stan-dard purely Lagrangian scheme for the calculation of (cid:15) g in thispaper. There currently does not exist a comprehensive and self-consistent hydrodynamical simulation to treat the three-dimensional problem of MT in a binary. Even when three-dimensional simulations eventually become self-consistent,they would likely resolve only the question of the initialstability of the Roche lobe overflow (RLOF), but not thelong-term MT, which will remain a prerogative of one-dimensional stellar codes in the foreseeable future.In a standard one-dimensional ML model, the stellarcodes use the regular set of structure equations, but adopt-ing the boundary condition that the total mass decreaseswith time. This boundary condition is an unavoidable re-duction of the three-dimensional picture of the ML to onedimension. In other words, ˙ M represents our best under-standing of the stream that is formed in the vicinity of L and that carries the donor’s material away from the donor.In this Section we examine the reaction of giants to the MLin one-dimensional stellar codes, explaining in particular thenature of the feature observed in the previous MT calcula-tions. Two completely opposite responses to ML were found indetailed 1D simulations by Passy et al. (2012). According tothem, 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 experimentis close to what would happen in Nature, where MT neverstarts abruptly at some fixed high MT rate. However, it isimportant to understand what causes the dramatic differ-ence between these two approaches.Passy et al. (2012) provided the following explanation: ”...some energy that is stored in gravitational form in thehydrostatic models is actually in a kinetic form [in hydrody-namic models], leading to the star contracting instead of ex-panding” , although how exactly the transformation of grav-itational energy into kinetic leads to contraction was notexplained.Instead of the energy argument, we argue that the mainreason for the radial response is due to the material beingconsumed from the surface of a red giant with linear velocity˙ m/ (4 πr ρ ), where r is the radius of a red giant, and ρ is thesurface density. Note that this term intrinsically decreasesthe radius.To validate that this term is dominant, we need to con-sider the involved stellar equations. A standard way to in-troduce hydrodynamical treatment into a stellar code is touse a truncated Navier-Stokes equation in spherical coordi-nates to calculate the hydrodynamical term of the pressurederivative, for example: (cid:20) ∂P∂m (cid:21) t = − Gm πr (1 + Q ) , (8)where the acceleration term Q is the ratio of the local La- c (cid:13) , 000–000 ass transfer from giant donors -5 -4
20 24 28 32 Q s / ( k N A ) (M init - m) [M ☉ ] Figure 8.
A red giant of initial mass 5.0 M (cid:12) and radius 50 R (cid:12) subjected to ML of 0.1 M (cid:12) yr − . The entropy profile is shownwith a thin line, the ratio Q of acceleration and gravitationalterms is shown with a thick line. It is clear that Q is far from beingnegligible compared to unity within the superadiabatic layer andsurface radiative zone. grangian acceleration to the local gravitational acceleration: Q = (cid:20) ∂ r∂t (cid:21) m / Gmr = a ( m, t ) GM/r , (9)where a ( m, t ) is the local Lagrangian acceleration. For a starin hydrostatic equilibrium, Q = 0. In Figure 8 we illustratehow significant the acceleration term can be in the supera-diabatic layer and surface radiative zone of a red giant athigh ML rates.Let’s now consider the initial response after some (smallcompared to the dynamical timescale) time τ to ML ˙ m < M without artificial viscosityor any other effects that would alter Equation (8).If P ( m, t ) is a smooth and bounded function and m ( t )and r ( t ) are continuous and bounded, thenlim t → a ( m, t ) = 4 πR (cid:20) ∂P∂m (cid:21) t =0 + GM R (10)By integrating Equation (10), we get for the velocity v atthe surface : v ( t ) = v + (cid:18) πR (cid:20) ∂P∂m (cid:21) t =0 + GM R (cid:19) t (11)+ o ( t ) , where v = v (0). To obtain the complete radial response,this velocity must be combined with the material consump-tion velocity, which is equal to v c ( t ) = ˙ m πρ R + o ( t ) , (12)where ρ is the initial surface density of the star. Their summust in turn be integrated in time. Thus, the complete radial Here we use the ”little-o” notation for the vicinity of zero.A brief definition of this notation is f ( t ) = o ( g ( t )) ⇔ lim t → f ( t ) /g ( t ) = 0 (for details see, e.g., Kevorkian & Cole1985). response of a red giant to ML is given by: R ( t ) = R + (cid:18) ˙ m πR ρ + v (cid:19) t (13)+ (cid:18) πR (cid:20) ∂P∂m (cid:21) t =0 + GM R (cid:19) t + o ( t )If the original star was in hydrostatic equilibrium, then v = 0, and, according to Equation (10), the term, propor-tional to t in Equation (13), would be equal to o (1) t = o ( t ), and hence is vanishing. The total radial response re-duces to R ( t ) = R + ˙ m πR ρ t + o ( t ) (14)This result may be directly compared to the decreasein the radius of a stellar model that occurs during the firsttimestep after the start of the ML, divided by this timestep.This comparison helps to understand the role of the adoptedouter boundary condition (see Sections 3.2 and 5.1). Weconclude that the contraction observed in the simulations ofPassy et al. (2012) is neither a numerical error, nor a conver-sion of gravitational energy to kinetic, but a natural conse-quence of a finite pace of expansion of the one-dimensionalenvelope that was forced to experience a fast ML. As shown in Figure 8, the acceleration term is becomingdangerously large in the surface layers during fast ML. Inthis Section we describe the modification of stellar equationsand the boundary condition to take this term into account.To obtain the pressure at the outer boundary P bc in asimple grey atmosphere approximation, it is common to usethe following well-known formula : P bc = 23 GMR κ (15)Unfortunately, the above formula is based on the equa-tion of hydrostatic equilibrium, hence in our case it is notapplicable, not to speak that it is inconsistent with the Equa-tion 8 that is used in many stellar codes. For this reason wetook into account the acceleration term not only in the equa-tion of motion (8), but also to find the boundary conditionfor a simple grey atmosphere, which we use in our experi-mental version of MESA : P bc = 23 GMR κ (1 + Q ) (16)Numerical experiments with MESA showed that the choice ofthe surface boundary condition affects the resulting radialresponse. The modified boundary condition must be usedbecause in this case we obtain a radial response substantiallycloser to that predicted by Equation 14 (see Section 5.1).Similarly, we take into account the acceleration term for Here for simplicity we omit radiation pressure at zero opticaldepth. This term may be important if radiative pressure is a sub-stantial fraction of the total pressure.c (cid:13)000
A red giant of initial mass 5.0 M (cid:12) and radius 50 R (cid:12) subjected to ML of 0.1 M (cid:12) yr − . The entropy profile is shownwith a thin line, the ratio Q of acceleration and gravitationalterms is shown with a thick line. It is clear that Q is far from beingnegligible compared to unity within the superadiabatic layer andsurface radiative zone. grangian acceleration to the local gravitational acceleration: Q = (cid:20) ∂ r∂t (cid:21) m / Gmr = a ( m, t ) GM/r , (9)where a ( m, t ) is the local Lagrangian acceleration. For a starin hydrostatic equilibrium, Q = 0. In Figure 8 we illustratehow significant the acceleration term can be in the supera-diabatic layer and surface radiative zone of a red giant athigh ML rates.Let’s now consider the initial response after some (smallcompared to the dynamical timescale) time τ to ML ˙ m < M without artificial viscosityor any other effects that would alter Equation (8).If P ( m, t ) is a smooth and bounded function and m ( t )and r ( t ) are continuous and bounded, thenlim t → a ( m, t ) = 4 πR (cid:20) ∂P∂m (cid:21) t =0 + GM R (10)By integrating Equation (10), we get for the velocity v atthe surface : v ( t ) = v + (cid:18) πR (cid:20) ∂P∂m (cid:21) t =0 + GM R (cid:19) t (11)+ o ( t ) , where v = v (0). To obtain the complete radial response,this velocity must be combined with the material consump-tion velocity, which is equal to v c ( t ) = ˙ m πρ R + o ( t ) , (12)where ρ is the initial surface density of the star. Their summust in turn be integrated in time. Thus, the complete radial Here we use the ”little-o” notation for the vicinity of zero.A brief definition of this notation is f ( t ) = o ( g ( t )) ⇔ lim t → f ( t ) /g ( t ) = 0 (for details see, e.g., Kevorkian & Cole1985). response of a red giant to ML is given by: R ( t ) = R + (cid:18) ˙ m πR ρ + v (cid:19) t (13)+ (cid:18) πR (cid:20) ∂P∂m (cid:21) t =0 + GM R (cid:19) t + o ( t )If the original star was in hydrostatic equilibrium, then v = 0, and, according to Equation (10), the term, propor-tional to t in Equation (13), would be equal to o (1) t = o ( t ), and hence is vanishing. The total radial response re-duces to R ( t ) = R + ˙ m πR ρ t + o ( t ) (14)This result may be directly compared to the decreasein the radius of a stellar model that occurs during the firsttimestep after the start of the ML, divided by this timestep.This comparison helps to understand the role of the adoptedouter boundary condition (see Sections 3.2 and 5.1). Weconclude that the contraction observed in the simulations ofPassy et al. (2012) is neither a numerical error, nor a conver-sion of gravitational energy to kinetic, but a natural conse-quence of a finite pace of expansion of the one-dimensionalenvelope that was forced to experience a fast ML. As shown in Figure 8, the acceleration term is becomingdangerously large in the surface layers during fast ML. Inthis Section we describe the modification of stellar equationsand the boundary condition to take this term into account.To obtain the pressure at the outer boundary P bc in asimple grey atmosphere approximation, it is common to usethe following well-known formula : P bc = 23 GMR κ (15)Unfortunately, the above formula is based on the equa-tion of hydrostatic equilibrium, hence in our case it is notapplicable, not to speak that it is inconsistent with the Equa-tion 8 that is used in many stellar codes. For this reason wetook into account the acceleration term not only in the equa-tion of motion (8), but also to find the boundary conditionfor a simple grey atmosphere, which we use in our experi-mental version of MESA : P bc = 23 GMR κ (1 + Q ) (16)Numerical experiments with MESA showed that the choice ofthe surface boundary condition affects the resulting radialresponse. The modified boundary condition must be usedbecause in this case we obtain a radial response substantiallycloser to that predicted by Equation 14 (see Section 5.1).Similarly, we take into account the acceleration term for Here for simplicity we omit radiation pressure at zero opticaldepth. This term may be important if radiative pressure is a sub-stantial fraction of the total pressure.c (cid:13)000 , 000–000
Pavlovskii & Ivanova the temperature equation (cid:20) ∂T∂m (cid:21) t = − Gm πr (1 + Q ) TP ∇ , (17)for convective conductivity in the Mixing Length Theory,and for the radiative gradient: ∇ rad = 3 κP L πacT GM (1 + Q ) . (18) We describe below an optically thick model for MT that wehave adopted for our studies of hydrodynamic response torapid MT in binaries.We follow a conventional way to calculate the MT rateby integrating the mass flow over a “nozzle” cross-sectionthat is taken on a plane perpendicular to the line connectingthe centres of the two stars and passing through the L point: ˙ m = (cid:90) nozzle v flow ρ flow dS (19)Here v flow is the velocity that the stream has within thenozzle, and ρ flow is the density that the stream has at thesame position. We note that the use of the nozzle cross-section only in the L neighbourhood implies that only the L MT rate can be found, and the rate of ML via L /L overflow cannot be calculated.Furthermore, any scheme that finds an ML boundarycondition for stellar one-dimensional codes adopts some setof simplifications to find the distribution of the stream’s den-sity and velocity throughout the nozzle, as well as the nozzlegeometry. The flow is considered to be steady, which per-mits the use of the Bernoulli theorem along the streamlines.The standard assumptions are that the nozzle at L coin-cides with the sonic surface of the flow (Lubow & Shu 1976),and that initial velocities are negligible at the origin of thestreamlines. We adopt the same assumptions.The MT rate then becomes dependent only on theseassumptions:(i) the adopted geometry of the donor and the nozzle;(ii) the streamlines – this assumption, coupled with theadopted evolution of the specific entropy along thestreamlines, allows relating the donor’s thermodynam-ical properties to those of the flow crossing the nozzle. The fundamental simplifying assumption that governs thewhole Roche lobe formalism is volume correspondence . Moreprecisely, this is the assumption that thermodynamical pa-rameters (pressure, density, composition, etc.) at a certainradius r in a one-dimensional stellar model are the same inthree-dimensional space at a Roche equipotential whose en-closed volume is (4 / πr . Whenever one-dimensional stellarevolution is considered in terms of RLOF formalism, the vol-ume correspondence assumption is automatically used. Notethat when a donor experiences a substantial RLOF this vol-ume correspondence for thermodynamical parameters can be applied only far from the L point, where the donor’smaterial is almost at rest. The Roche lobe volume radius R L is usually found using one of two well-known approxi-mations (Paczy´nski 1971; Eggleton 1983). The area of thenozzle is then usually approximated by the area of an ellip-soid by taking the second-order term in the Roche potentialexpansion within the sonic surface near the L point.An additional simplifying assumption is applied whenintegrating the ML rate over the potential region betweenthe Roche lobe potential Φ L and the photospheric Rochepotential Φ phot (note, that Φ phot > Φ L ). At this point, toavoid volume integrations, it is a common approach to im-plicitly break the volume correspondence assumption andreplace it with the pressure correspondence assumption.The pressure correspondence assumption asserts thatthe thermodynamical parameters (temperature, density andcomposition) at a point A in three-dimensional space are thesame as at a point B in the one-dimensional model, providedthat the pressure at the point A in three-dimensional spaceis the same as the pressure at point B in the one-dimensionalmodel. If one applies this pressure correspondence far fromthe L point and assumes, in addition, that a · ∇ Φ (cid:28) | ∇ Φ | , (20)where a is the local Lagrangian acceleration, and Φ is the lo-cal Roche potential, then it becomes possible to replace theintegral over radius with the integral over pressure. The ben-efit of the pressure correspondence assumption is that onecan avoid the painful calculation of the differential d (Φ( r )),which requires volume integrations. We note that once a staris not in hydrostatic equilibrium, Q (cid:28) L point break downsince the expansions are truncated after the quadratic term.Higher-order terms are no longer negligible. Instead, we con-duct the realistic Roche lobe integrations which employ theRunge-Kutta-Fehlberg integrator and the damped Newton-Raphson solver to obtain all geometrical parameters. Theseintegrations have been conducted for 275001 mass ratiosfrom 0.06 to 19.145 with an increment of 7 · − . For eachmass ratio we calculate the nozzle areas and volume radiifor 200 equipotentials lying between the L potential andthe L /L potential. In addition to the use of a more pre-cise relation between the star’s volume and the Roche lobevolume, volume correspondence that is necessary for rapidmass loss and a non-simplified nozzle shape, this also allowsus to track whether the donor overfills the L /L equipo-tential during the mass transfer. For a rapid mass transfer rate, an ”optically thick” approx-imation is usually used (see Paczy´nski & Sienkiewicz 1972;Savonije 1978; Kolb & Ritter 1990; Ge et al. 2010, and manyothers). The flow of matter towards L in this approximation c (cid:13) , 000–000 ass transfer from giant donors is adiabatic and streamlines go along the equipotential sur-faces of the Roche potential. The photosphere correspondsto a Roche equipotential which lies outside the Roche lobe.Close to the L point, the photosphere turns into the outerboundary of the optically thick flow which flows across thesonic surface along the Roche equipotential surfaces withthe local adiabatic speed of sound (Kolb & Ritter 1990).To find the stream’s density and sonic velocity at L ,one is required to adopt the evolution of the specific entropyin the flow along the streamlines. Often the flow is takento be polytropic – that is, the flow preserves a constantvalue of P/ρ γ along a streamline, where γ is the adiabaticexponent. In certain cases, the polytropic stratification ofthe donor itself is adopted. Note that the stream’s velocity,while locally sonic, is not constant across the nozzle. Theflow also does not need to be isentropic if the donor is notisentropic – the streamlines can originate at different initialequipotentials.We also assume that a flow is adiabatic. An adiabaticflow is not, of course, entirely polytropic due to, for example,recombination, that occurs as gas flows towards the sonicsurface. It means that γ varies along an adiabat. For thisreason, we employ a realistic equation of state taken from MESA , that among other effects takes into account the ioni-sation of elements in the mix and the radiative componentof pressure. For test purposes, we also can consider the flowto be polytropic. We have verified that for those modelsthat are provided as examples in this paper, the effect ofthe adopted equations of state on mass flux does not exceed ± L neighbourhood. Hence, wedo not use this approach.The degree to which optically thin mass transfer affectsthe behaviour of stars before their photospheres overfill theRoche lobe depends on the rate of their intrinsic evolution-ary expansion when their photospheres approach the Rochelobe. Relatively massive giants (5 M (cid:12) and up), cross theinterval of radii for which the optically thin mass transferis dominant in a short time thanks to their fast expansion.Hence, the fraction of the envelope they lose via opticallythin transfer is small. On the other hand, less massive gi-ants epxand slowly and can lose a lot of mass through at-mospheric ML before the optically thick mechanism kicksin. This decreases the mass ratio at the onset of RLOF Table 1.
Simplifications eliminated in the optically thick masstransfer schemesReference GS PC PD PSPaczy´nski & Sienkiewicz (1972)Savonije (1978)Kolb & Ritter (1990) • Ge et al. (2010) ◦ •
This work • • • •
GS – geometrical simplification for the nozzle, PC – pressure cor-respondence, PD – polytropic stratification of the donor, PS –along the streamlines
P/ρ γ is constant. Empty circle – simplifiedpotential is adopted the effect of which is equivalent to the PCassumption. Figure 9.
Models passy-a and passy-b. A 102 R (cid:12) , 0 . M (cid:12) red giant similar to the one considered by Passy et al. (2012) issubjected to constant ML of 0 . M (cid:12) yr − . The solid line repre-sents the result we reproduced following the method described byPassy et al. (2012), which is very similar to the one shown by ared line in their Figure 4. The dashed line represents the resultwe obtained by compensating for the limitations discussed in Sec-tion 2.4, and after removing some pulsation artefacts by disablingcomposition smoothing at the bottom of the convective envelope. and improves stability. We observed an extreme case, wherewith the Ritter (1988) prescription, a 0 . M (cid:12) giant couldnever actually overfill a 100 R (cid:12) Roche lobe, and steadily lostthe whole convective envelope via optically thin mechanismreaching an optically thin MT rate of 10 − . M (cid:12) yr − ! Wetherefore warn the reader that the mode of mass transferfrom low-mass giants might be sensitive to the very uncer-tain model of optically thin MT.We summarise the simplifications eliminated in the ex-isting optically thick schemes and in our case in Table 1. First, we can compare the initial reaction of a red giantto constant ML to our prediction made in Section 3.1. Forexample, we take a 5 M (cid:12) and 50 R (cid:12) giant and subject itto the constant ML rate ˙ M = 10 − . M (cid:12) yr − . The initialderivative of radius predicted by formula (14) is about − . · c (cid:13)000
Models passy-a and passy-b. A 102 R (cid:12) , 0 . M (cid:12) red giant similar to the one considered by Passy et al. (2012) issubjected to constant ML of 0 . M (cid:12) yr − . The solid line repre-sents the result we reproduced following the method described byPassy et al. (2012), which is very similar to the one shown by ared line in their Figure 4. The dashed line represents the resultwe obtained by compensating for the limitations discussed in Sec-tion 2.4, and after removing some pulsation artefacts by disablingcomposition smoothing at the bottom of the convective envelope. and improves stability. We observed an extreme case, wherewith the Ritter (1988) prescription, a 0 . M (cid:12) giant couldnever actually overfill a 100 R (cid:12) Roche lobe, and steadily lostthe whole convective envelope via optically thin mechanismreaching an optically thin MT rate of 10 − . M (cid:12) yr − ! Wetherefore warn the reader that the mode of mass transferfrom low-mass giants might be sensitive to the very uncer-tain model of optically thin MT.We summarise the simplifications eliminated in the ex-isting optically thick schemes and in our case in Table 1. First, we can compare the initial reaction of a red giantto constant ML to our prediction made in Section 3.1. Forexample, we take a 5 M (cid:12) and 50 R (cid:12) giant and subject itto the constant ML rate ˙ M = 10 − . M (cid:12) yr − . The initialderivative of radius predicted by formula (14) is about − . · c (cid:13)000 , 000–000 Pavlovskii & Ivanova cm s − . With modified boundary conditions and with atimestep of 10 − yr, the code produces − . · cm s − . Ifinstead the default hydrostatic boundary condition is used,the code produces − . · cm s − .Second, to compare our results to the ones publishedearlier, we calculate the behaviour of the 0 . M (cid:12) and102 R (cid:12) red giant examined by Passy et al. (2012). As inPassy et al. (2012), we subject this giant to the constantML rate ˙ M = 0 . M (cid:12) yr − . We evolve one giant using theunmodified MESA code with viscosity and timestep setup de-scribed in Passy et al. (2012); we call this the “passy-a”model. For the other giant, “passy-b” model, we take intoaccount the circumstances outlined in Section 2.4. In par-ticular, to obtain correct (cid:15) g , we evolve this model using aconstant time step of 10 − yr. We also remove some pul-sation artefacts by disabling composition smoothing at thebottom of the convective envelope.The initial behaviour of the radius in both models hasbeen explained in Section 3.1. Note, however, that in ourmodel, the stellar radius is always smaller than that in Passyet al. (2012) (see our Figure 9 and their Figure 4). We findthat the difference in the value of the radii is mainly due tohow accurately the value of (cid:15) g is found.In addition to the initial radius contraction, one cannotice ensuing radial oscillations, visible both in the modelspassy-a and passy-b. We define those non-numerical radialpulsations, excited either by the start of constant ML or,in the case of evolution in a binary, by the rapid growthof the ML rate, as mass loss induced pulsations (MLIPs).These pulsations might be caused by the sonic rarefactionwave, reflected from the bottom of the envelope, which istheoretically predicted to occur in a fluid as it abruptly ex-pands into vacuum, hence we think that they are not of nu-merical nature. In the model passy-a, these pulsations arelargely smoothed out because the timestep in this modelgrows after the start of simulations and exceeds the dynam-ical timescale.A crucial parameter that defines the level of importanceof MLIP is the p-mode damping rate η p , which characterisesthe timescale on which a giant roughly attains hydrostaticequilibrium and dynamical oscillations are damped. Withcontemporary high-precision photometric instruments suchas COROT and Kepler, it is possible to obtain high-precisionmeasurements of the profile widths (i.e, damping rates di-vided by π ) of p-modes (see e.g. Baudin et al. 2011; Belka-cem et al. 2012). Damping rates of intrinsic pulsations ofred giants found from those observations are η p ∼ − π to10 π .The damping rates, given that the timestep used ismuch smaller than the pulsation period, depend hugely onthe artificial viscosity coefficients: increasing the artificialviscosity increases the damping rates. The damping rates ofMLIPs shown in Figure 9 are comparable, within an order ofmagnitude, with those obtained directly from observations.Note that the default artificial viscosity used for the modelspassy-a and passy-b is very moderate ( l = 0 . l = 0) anddoes not substantially affect the damping rates. Due to thelack of wide-range observational calibrations for dampingrates across giants of different radii and masses, we do notuse an artificial viscosity in our models, except in modelspassy-a and passy-b, where it is taken into account only for ζ ln(M tot ) [M ☉ ]-1 0 1 2 3 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 Figure 10. ζ for a giant of initial mass 5.0 M (cid:12) and radius50 R (cid:12) subjected to various rates of mass loss (10 − , 10 − and10 − M (cid:12) yr − (solid, dashed and dot-dashed lines respectively). ζ Core mass fraction
Figure 11. ζ for giants at MT rate 10 − M (cid:12) yr − . 5.0 M (cid:12) R (cid:12) (solid line), 1.0 M (cid:12) R (cid:12) (dashed line), 5.0 M (cid:12) R (cid:12) (dot-ted line). ζ comp for a composite polytrope ( n = 3 / the purpose of comparison with the original paper of Passyet al. (2012).We have performed several simulations where we sub-jected giants of several initial masses and radii to constantML, to determine the realistic stellar response defined as: ζ ≡ (cid:18) ∂ log R∂ log M (cid:19) . (21)Note that generally ζ is an implicit function of not only thecurrent MT rate, but also of the previous MT, and here welook at ζ for constant MT rates only. We found that:(i) At high ML rates, ζ saturates. It becomes mainly afunction of the core mass fraction and not of the MLrate (see Figure 10 where we show representative ex-amples);(ii) In “well-expanded” giants, the saturated ζ ap-proaches the behaviour of a composite polytrope ζ comp ,as in Hjellming & Webbink (1987), with n = 3 / ζ is higher thanwhat is expected from the composite polytrope.(iii) Non-saturated values of ζ are usually lower than thesaturated ζ . c (cid:13) , 000–000 ass transfer from giant donors (iv) ζ at the very onset of RLOF is often hard to determinedue to MLIPs.While the radial response seems to be most natural tolook at during the MT, we find that other properties ofthe donor’s surface also depend dramatically on how well (cid:15) g is obtained. In Section 2.5 we showed that substantialluminosity can be generated if recombination energy is cal-culated properly. We can compare the surface luminosi-ties of 5 M (cid:12) and 50 R (cid:12) giant, subjected to a constant ML˙ M = 0 . M (cid:12) yr − and evolved with a default time-step ad-justment and with a time-step chosen to resolve (cid:15) g . An esti-mate provided in Section 2.5 shows that the energy comingfrom recombination in a giant subjected to ML is approxi-mately W rec ≈ L (cid:12) ˙ MM (cid:12) yr − . (22)Where does this energy mainly go – is it radiated awayas excess luminosity or spent on the star itself, e.g. , to in-crease its thermal or gravitational energy? The unperturbedgiant has luminosity only about L donor ∼ L (cid:12) . The cal-culations show that a mass-losing giant in which (cid:15) g in thesuperadiabatic layer was calculated properly, has a luminos-ity about 6 times higher than the mass-losing giant with adefault relatively large time-step! The first giant has alsobecome about 50% hotter. This increase in luminosity indi-cates that a large portion of energy from recombination isradiated away, rather than being spent on the star itself. Inthe case of a 5 M (cid:12) and 200 R (cid:12) giant, a simple prediction fora luminosity increase and an exact calculation would agreefairly well – in a more expanded giant, less recombinationenergy was spent on the star itself, and most went directly tosurface luminosity. Similarly, the luminosities in the passy-aand passy-b models are different by about 6 times, also as asimple estimate would predict. Interestingly, in the case of5 M (cid:12) giants, the radii obtained by the two methods were notmuch different, except that in giants evolved with a smalltime-step we can observe MLIPs.Of course, very fast MT rates are not commonly ob-served. However, the path of a star through the fast MT isdifferent, and can define both the stability and/or how thedonor appears when the fast MT phase is completed. Wealso can evaluate that this effect might be important (i.e.provides a difference in the luminosity by more than a fewper cent) for as long as˙ M > ∼ − M (cid:12) yr × L donor L (cid:12) (23)and hence might be important for low-luminosity giantdonors and, in general, for all donors transferring mass atthe thermal timescale of their envelope.Indeed, consider the thermal timescale τ KH of the enve-lope taken as usual τ KH ≈ × yr (cid:18) MM (cid:12) (cid:19) R (cid:12) R L (cid:12)
L . (24)For most giants
R/R (cid:12) (cid:29) M/M (cid:12) , and their envelope ther-mal timescale MT,
M/τ KH , satisfies condition 23. The conventional way to determine stability of a binary sys-tem to the MT is to compare the initial responses of thedonor radius and the Roche lobe radius to the mass loss.As was mentioned in Section 1, such analysis, if it does notinvolve actual MT calculations in a binary, is usually doneby comparing their mass-radius exponents at the onset ofRLOF. In order to do this, the adiabatic radial response ofthe donor ζ ad is found from one of the existing approxima-tions (composite polytropes, condensed polytropes, or evensimply adopting that for convective donors ζ ad = 0). Thisadiabatic radial response is used instead of the realistic stel-lar response.However, we’d like to stress, that even if the realisticstellar response is obtained, there is no reason to assumethat if ζ < ζ L at the start of the MT, and a donor’s relativeRoche lobe overflow initially increases, the instability willnecessarily occur. Instead of considering just the momentof RLOF, we can trace binary evolution through the MTin detail and determine if the instability eventually takesplace, or not. The determination of MT instability in a bi-nary system, while the donor is overflowing its RL, requiresan alternative (to a simple comparison of the initial ζ s) cri-terion to delineate when the MT proceeds stably, and whenit results in a common envelope.To understand the characteristic behaviour ofthe systems during MT, we performed a large setof simulations with donors of several initial masses(1 , , , , , , M (cid:12) ). We considered several values ofradii for each donor, taken within the range where thedonors have non-negligible, > ∼ . M donor outer convectiveenvelopes, both at the red giant and asymptotic giantbranches. We conducted simulations for several initial massratios (between 0.9 and 3.5). We assume that a companionis compact, where the compactness only means that weneglect any possible accretor’s RLOF.Further, for these detailed MT sequences we adoptedfully conservative MT; any non-conservation in the MTwould lead to a relative increase of MT stability. Most ofour binary systems were evolved using our modifications:obtaining proper (cid:15) g , and our method to find MT rates andboundary conditions.Even for systems with larger initial q than that dic-tated by the conventional adiabatic stability criterion, theMT rate after RLOF smoothly increases, approaches a peakvalue and then decreases; the peak MT rate is usually muchless than the dynamical MT rate. With the increase of q , thepeak MT rate increases. We define as q L the smallest initialmass ratio for which the binary experiences L /L -overflowduring the MT (see Figure 12). We refer to this overflow as L /L -overflow, as during the MT, the mass ratio can re-verse, and L -overflow becomes L -overflow. This q L is afunction of both the initial mass ratio and the donor’s ra-dius. For many initial donor radii, in binary systems with q = q L , the mass of the donor encompassed between thedonor’s Roche lobe and its L /L lobe is very small andthe ML rate when L /L -overflow is approached is far fromdynamical. Numerically, we are capable of evolving such sys-tems with q > q L through the L /L -overflow, and manyof them will not even approach the dynamical MT rate (seeFigure 12). We, however, do not trust the MT rates obtained c (cid:13)000
M/τ KH , satisfies condition 23. The conventional way to determine stability of a binary sys-tem to the MT is to compare the initial responses of thedonor radius and the Roche lobe radius to the mass loss.As was mentioned in Section 1, such analysis, if it does notinvolve actual MT calculations in a binary, is usually doneby comparing their mass-radius exponents at the onset ofRLOF. In order to do this, the adiabatic radial response ofthe donor ζ ad is found from one of the existing approxima-tions (composite polytropes, condensed polytropes, or evensimply adopting that for convective donors ζ ad = 0). Thisadiabatic radial response is used instead of the realistic stel-lar response.However, we’d like to stress, that even if the realisticstellar response is obtained, there is no reason to assumethat if ζ < ζ L at the start of the MT, and a donor’s relativeRoche lobe overflow initially increases, the instability willnecessarily occur. Instead of considering just the momentof RLOF, we can trace binary evolution through the MTin detail and determine if the instability eventually takesplace, or not. The determination of MT instability in a bi-nary system, while the donor is overflowing its RL, requiresan alternative (to a simple comparison of the initial ζ s) cri-terion to delineate when the MT proceeds stably, and whenit results in a common envelope.To understand the characteristic behaviour ofthe systems during MT, we performed a large setof simulations with donors of several initial masses(1 , , , , , , M (cid:12) ). We considered several values ofradii for each donor, taken within the range where thedonors have non-negligible, > ∼ . M donor outer convectiveenvelopes, both at the red giant and asymptotic giantbranches. We conducted simulations for several initial massratios (between 0.9 and 3.5). We assume that a companionis compact, where the compactness only means that weneglect any possible accretor’s RLOF.Further, for these detailed MT sequences we adoptedfully conservative MT; any non-conservation in the MTwould lead to a relative increase of MT stability. Most ofour binary systems were evolved using our modifications:obtaining proper (cid:15) g , and our method to find MT rates andboundary conditions.Even for systems with larger initial q than that dic-tated by the conventional adiabatic stability criterion, theMT rate after RLOF smoothly increases, approaches a peakvalue and then decreases; the peak MT rate is usually muchless than the dynamical MT rate. With the increase of q , thepeak MT rate increases. We define as q L the smallest initialmass ratio for which the binary experiences L /L -overflowduring the MT (see Figure 12). We refer to this overflow as L /L -overflow, as during the MT, the mass ratio can re-verse, and L -overflow becomes L -overflow. This q L is afunction of both the initial mass ratio and the donor’s ra-dius. For many initial donor radii, in binary systems with q = q L , the mass of the donor encompassed between thedonor’s Roche lobe and its L /L lobe is very small andthe ML rate when L /L -overflow is approached is far fromdynamical. Numerically, we are capable of evolving such sys-tems with q > q L through the L /L -overflow, and manyof them will not even approach the dynamical MT rate (seeFigure 12). We, however, do not trust the MT rates obtained c (cid:13)000 , 000–000 Pavlovskii & Ivanova
70 80 90 100 110 120 0 0.05 0.1 0.15 0.2-10-8-6-4-2 0 R a d i u s [ R ☉ ] l o g ( M . [ M ☉ / y r ]) Time since M. > 10 -4 M ☉ /yr [yr]1.6 1.5 1.4 Figure 12.
The MT rate (green solid lines), and the evolutionof the donor’s radius (dashed lines) and the donor’s Roche loberadius (solid black lines), shown for the mass ratio with whichthe binary system does not start L /L -overflow during the MT( q = 1 . L / L q = 1 . , . L /L overflow is indicated with solid black colour on the MTcurves. The donor is a 1 M (cid:12) and 100 R (cid:12) giant at the onset of theMT. The rate of a dynamical MT for this system is ∼ M (cid:12) yr − . in this regime because our stream model only considers the L -nozzle.We adopt therefore that for as long as L /L -overflowof the donor does not happen, the system is stable. In part,this is confirmed by the three-dimensional hydrodynamicalsimulations, in which a common envelope is always asso-ciated with a severe, albeit very short in duration, L /L overflow of the donor (Nandez et al. 2014). Note however,that we consider this to be a minimum stability criterion, asthe MT in a binary system may remain stable even after thedonor’s L /L -overflow takes place, as, e.g., SS 433 shows(e.g., Bowler 2010; Perez M. & Blundell 2010).We also compared the surface-averaged effective grav-itational acceleration over the L /L lobe with the spheri-cally symmetric acceleration that would be expected at itsvolume radius and found that the difference is about 13%.In addition to L /L -overflow condition, we tracewhether the binary orbit is changing rapidly. For the lattercondition, we adopt that if | ˙ a/a | T < /
50, the orbit is notchanging rapidly, where T is the binary period. Rapid evo-lution of orbital parameters can invalidate the entire Rocheformalism. For example, if angular velocity changes substan-tially over one period, the Euler term of the fictitious force,which is ignored in the Roche formalism along with the Cori-olis term, becomes comparable to the centrifugal term. Note,however, that in Nature, binary systems with a larger | ˙ a/a | T might not experience dynamical-timescale MT.The key reason behind the unexpected stability of theMT in systems with q ζ < q < q L is that once the MTstarts, ζ rises as the donor’s mass decreases (see Figures 10,11 and 13). At the same time, ζ L decreases (Figure 13). Adecrease in ( ζ L − ζ ) means that, as the MT proceeds, the sta-bility of the system increases. We define as the critical mass-point , m cp , the mass of the donor at which ζ = ζ L . When thedonor decreases its mass to the critical mass-point, the ra-tio of the donor radius to its Roche lobe starts to decrease.Therefore, assuming that the relative RL overflow at the L /L equipotential weakly depends on the mass ratio (see ζ ; ζ L ln(M tot ) [M ☉ ] 3.0 1.8-1 0 1 2 3 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.60.9 Figure 13. ζ for a giant of initial mass 5.0 M (cid:12) and radius 50 R (cid:12) subjected to a constant mass loss of 10 − M (cid:12) yr − (dot-dashedline). Solid lines show ζ L for the conservative MT at various initialmass ratios. Intersections ζ = ζ L define the critical mass-points m cp , at which mass transfer starts to decrease. The shaded areais equal to the right part of the Equation 26 for q = 1 . m = m cp . below), if there was no L /L before the donor mass de-creased to m cp , the binary system is stable with respect toMT. We have searched for the point when L /L -overflowstarts using our extended set of detailed simulations (seeFigure 14). It can be seen that the critical mass ratios q L for which the MT becomes “unstable” are about twice aslarge as would be predicted by conventional comparison ofmass-radius exponents at the onset of the MT.Let us analyse the results and show how the approx-imate critical mass ratios can be predicted without doingdetailed binary calculations.To predict whether a system experiences L /L -overflow during the MT, we need to identify whether at anymoment after the start of MT the donor’s radius exceeds R L23 – the volume radius at which the donor starts L /L -overflow. In order to do this, we can find the ratio R L23 /R L for a range of the donor masses. Note that this critical ra-tio also depends on the mass ratio. However, we find thatthis dependence is quite weak, ln( R L23 ) /R L ≈ . ± . . ≤ q ≤
4. Nevertheless, we took this dependence intoaccount in our calculations by examining the ratio R L23 /R L not only down to donor mass m cp , but also further, almostall the way down to the donor core mass.The radius of the donor during the MT, when it hasshed to mass m , isln R ( m ) = ln R + (cid:90) mM ζ ( ˙ M ) dm . (25)It follows from Section 5.1 that in order to predict theradius of the donor at any moment during fast (saturated)MT, instead of the real ζ one can use ζ comp , which is ob-tained from the composite polytrope approximation. “Com-pact” giants, should be expected to expand less than thisestimate would predict, and hence be more stable becausetheir ζ > ζ comp .Using the definitions of mass-radius exponents (see alsoFigure 13), we can find the approximate radius of the donorat any moment after the start of MT as c (cid:13) , 000–000 ass transfer from giant donors I n i t i a l m a ss r a t i o Radius at RLOF [R ☉ ] C o n v e n t i o n a l c o n d e n s e d p o l y t r o p e t h r e s h o l d L /L -overflow condensed polytrope threshold 0.5 1 1.5 2 2.5 30 40 50 60 70 80 90 100 I n i t i a l m a ss r a t i o Radius at RLOF [R ☉ ]Conventional condensed polytrope thresholdL /L -overflow condensed polytrope threshold 0.5 1 1.5 2 2.5 150 200 250 300 350 400 450 I n i t i a l m a ss r a t i o Radius at RLOF [R ☉ ]Conventional condensed polytrope thresholdL /L -overflowcondensed polytropethreshold Figure 14.
The outcomes of the MT sequences for donors of dif-ferent initial masses. The top panel is for a 1 M (cid:12) donor, the middlepanel is for a 2 M (cid:12) donor and the bottom panel is for a 10 M (cid:12) donor. A square symbol indicates that | ˙ a | T/a > /
50, otherwise,the symbol is a circle. If a symbol is filled, then no L /L overflowis experienced, if a symbol is empty, then the system experiences L /L -overflow. Our L /L -overflow condensed polytrope sim-plification is shown with a dashed line, and q ζ – the conventionalcondensed polytrope threshold – is shown with a solid line. ln (cid:18) R ( m ) R L (cid:19) = (cid:90) mM ( ζ comp − ζ L ) d ln M . (26)Now we can find whether the system experiences L /L -overflow at any moment during the MT and thus producea semi-analytical estimate of q L – we will refer to thisestimate as “ L /L -overflow condensed polytrope simplifi-cation”. It is plotted in Figure 14 for various giants.We can see that the condensed polytrope simplifica-tion works best for those giants that are neither too com-pact nor too expanded. In both compact and average-sized I n i t i a l m a ss r a t i o log (R[R ☉ ]) at RLOF1 2 3 5 10 15 30 Figure 15.
Critical mass ratios obtained from the L /L -overflow simplified composite polytropes for the donors of dif-ferent initial masses and radii. giants that approach L /L -overflow at the critical point,most mass before the critical point is lost at ML rates whichare generally much higher than ∼ − M (cid:12) yr − , so it’sfine to use the saturated value of ζ . At the same time, com-pact giants have higher saturated values of ζ , which makesthem more stable than average-sized giants whose saturated ζ approximately follows the condensed polytrope.In large giants with very rarefied envelopes, L /L -overflow occurs at much lower mass loss rates, and the satu-rated values of ζ become inapplicable. Instead, a lower, non-adiabatic, non-saturated values should be used. This makesthese giants less stable. However, such L /L -overflow willbe non-dynamical and does not have to result in a commonenvelope.In Figure 15 we show how the critical mass ratio ob-tained from the L /L -overflow condensed polytrope sim-plification changes with the initial donor mass and radius.We have also checked and found that non-conservative MTleads, as expected, to the increase in q L , increasing it byabout 0.3.The detailed simulations (Figure 14) show that q L islarger when a donor is more compact and the mass fractionof its convective envelope is smaller. In all detailed simula-tions we did with donors where the outer convective envelopeis non-negligible, q L is below 3.5, where 3.5 is known as thecritical mass ratio leading to a delayed dynamical instabil-ity with radiative donors (see, e.g, Ge et al. 2010). It showstherefore that during the development of the outer convec-tive envelope, the convective donors are in the transitionalregime. When the convective envelope is well developed, thecritical mass ratios are in the range 1 . − . L /L overflowis definitely predicting that more binary systems will evolvethrough their MT in a stable way.Whether the L /L overflow should necessarily lead toa common envelope is not entirely clear. In principle, MTthrough the L /L nozzle can be treated in the same simpli-fied way as that through the L nozzle, and the material canbe considered to leave the system carrying away the specificangular momentum of the donor. It is a subject of our futurework. c (cid:13)000
Critical mass ratios obtained from the L /L -overflow simplified composite polytropes for the donors of dif-ferent initial masses and radii. giants that approach L /L -overflow at the critical point,most mass before the critical point is lost at ML rates whichare generally much higher than ∼ − M (cid:12) yr − , so it’sfine to use the saturated value of ζ . At the same time, com-pact giants have higher saturated values of ζ , which makesthem more stable than average-sized giants whose saturated ζ approximately follows the condensed polytrope.In large giants with very rarefied envelopes, L /L -overflow occurs at much lower mass loss rates, and the satu-rated values of ζ become inapplicable. Instead, a lower, non-adiabatic, non-saturated values should be used. This makesthese giants less stable. However, such L /L -overflow willbe non-dynamical and does not have to result in a commonenvelope.In Figure 15 we show how the critical mass ratio ob-tained from the L /L -overflow condensed polytrope sim-plification changes with the initial donor mass and radius.We have also checked and found that non-conservative MTleads, as expected, to the increase in q L , increasing it byabout 0.3.The detailed simulations (Figure 14) show that q L islarger when a donor is more compact and the mass fractionof its convective envelope is smaller. In all detailed simula-tions we did with donors where the outer convective envelopeis non-negligible, q L is below 3.5, where 3.5 is known as thecritical mass ratio leading to a delayed dynamical instabil-ity with radiative donors (see, e.g, Ge et al. 2010). It showstherefore that during the development of the outer convec-tive envelope, the convective donors are in the transitionalregime. When the convective envelope is well developed, thecritical mass ratios are in the range 1 . − . L /L overflowis definitely predicting that more binary systems will evolvethrough their MT in a stable way.Whether the L /L overflow should necessarily lead toa common envelope is not entirely clear. In principle, MTthrough the L /L nozzle can be treated in the same simpli-fied way as that through the L nozzle, and the material canbe considered to leave the system carrying away the specificangular momentum of the donor. It is a subject of our futurework. c (cid:13)000 , 000–000 Pavlovskii & Ivanova
In this paper, we considered a number of theoretical andpractical challenges on the path to understanding the sta-bility of the MT in binaries consisting of a convective giantdonor and a compact accretor.We find that in order to obtain the correct responseto the mass loss, the (cid:15) grav in the superadiabatic layer mustbe calculated properly; depending on which stellar code isused this may require a number of numerical tricks, and atthe very least control of the time-step to match the predictedenergy generation rate. We provide simple estimates for howto find the mass of this superadiabatic layer and the rate ofthe energy release, in order to quantify how well an arbitrarystellar code performs under mass loss.We have shown that the mass-radius exponents in thegiants that are compatible with the composite-polytrope de-scription converge to the prediction in Hjellming & Webbink(1987) for fast mass loss. On the other hand, the giantsthat cannot be described by a composite polytrope have ζ > ζ comp , and the binary systems with such donors aremore stable than the composite polytrope would predict. Inour detailed models, we could not find the strong superadi-abatic expansion predicted by the adiabatic models in Geet al. (2010).In addition to radial response, we examined the changesin surface luminosity and effective temperature of mass-losing stars. We find that at fast ML rates, the luminositycan differ by a factor of several times between the modelswhere (cid:15) grav in the superadiabatic layer is calculated correctlyand the ones where it is calculated in the Lagrangian waywith large timesteps. We note that in some observed binarysystems, the donor’s properties such as effective tempera-ture and mass are determined with a few % precision (e.g.,Leahy & Abdallah 2014). This difference in effective tem-perature obtained by the two models, and comparable tothe observed precision, is expected to occur in systems witha low-luminosity giant donor, and also in thermal-timescaleMT systems.We have enhanced the classical scheme to find the MTrate via an optically thick stream approximation. In particu-lar, our use of the detailed system geometry allows us to findthe MT rate in the case of a substantial RLOF. This leadsus to a new criterion of MT instability, based on whether thedonor starts L /L overflow. We have also found that a bi-nary system does not become immediately unstable as thedonor’s envelope becomes convective, but rather the massratio at which the instability occurs gradually decreases fromthe regime predicted for radiative donors. Note that in prin-ciple a binary can survive the MT even after L /L overflow,without starting a common envelope. However, we find thateven the new L /L -overflow criterion warrants that binarysystems will proceed with stable conservative MT if theirmass ratio is up to twice that given by the conventional cri-terion, for stars with deep convective zones (and the massratio 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 sup-ported in part by the National Science Foundation underGrant No. NSF PHY05-51164. Authors thank B. Paxton forhelp with MESA/star and MESA/eos modules, T. Fragos forhelp with MESA/binary module, C. Heinke for help with themanuscript and the anonymous referee for useful comments.Authors are also grateful to R. Webbink, Ph. Podsiadlowskiand S. Justham for useful discussions. This research has beenenabled by the use of computing resources provided by West-Grid and Compute/Calcul Canada.
REFERENCES
Baudin F., Barban C., Belkacem K., Hekker S., Morel T.,Samadi R., Benomar O., Goupil M.-J., Carrier F., BallotJ., Deheuvels S., De Ridder J., Hatzes A. P., Kallinger T.,Weiss W. W., 2011, A&A, 529, A84Belkacem K., Dupret M. A., Baudin F., Appourchaux T.,Marques J. P., Samadi R., 2012, A&A, 540, L7Bowler M. G., 2010, A&A, 521, A81Cox J., Giuli T., 1968, Principles of stellar structure. GOR-DON AND BREACH, Science Publishers, Inc.Eggleton P. P., 1967, MNRAS, 135, 243Eggleton P. P., 1971, MNRAS, 151, 351Eggleton P. P., 1972, MNRAS, 156, 361Eggleton P. P., 1973, MNRAS, 163, 279Eggleton P. P., 1983, ApJ, 268, 368Eggleton P. P., Faulkner J., Flannery B. P., 1973, A&A,23, 325Ferguson J. W., Alexander D. R., Allard F., BarmanT., Bodnarik J. G., Hauschildt P. H., Heffner-Wong A.,Tamanai A., 2005, ApJ, 623, 585Ge H., Hjellming M. S., Webbink R. F., Chen X., Han Z.,2010, ApJ, 717, 724Glebbeek E., Pols O. R., Hurley J. R., 2008, A&A, 488,1007Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R.,Ivanova N., 2002, MNRAS, 336, 449Hjellming M. S., Webbink R. F., 1987, ApJ, 318, 794Ivanova 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., TaamR., Tauris T. M., van den Heuvel E. P. J., Webbink R. F.,2013, A&Ar, 21, 59Ivanova N., Taam R. E., 2004, ApJ, 601, 1058Kevorkian J., Cole J. D., 1985. SpringerKippenhahn R., Weigert A., 1994, Stellar Structure andEvolution. Springer-Verlag Berlin Heidelberg New YorkKippenhahn R., Weigert A., Hofmeister E., 1967, MethodsComput. Phys., 7, 129Kolb U., Ritter H., 1990, A&A, 236, 385Kuhfuss R., 1986, A&A, 160, 116Leahy D. A., Abdallah M. H., 2014, ApJ, 793, 79Lubow S. H., Shu F. H., 1976, ApJ, 207, L53Nandez J. L. A., Ivanova N., Lombardi Jr. J. C., 2014, ApJ,786, 39Nelemans G., Tout C. A., 2005, MNRAS, 356, 753Nelemans G., Verbunt F., Yungelson L. R., Portegies ZwartS. F., 2000, A&A, 360, 1011Nelson C. A., Eggleton P. P., 2001, ApJ, 552, 664Osaki Y., 1970, ApJ, 162, 621 c (cid:13) , 000–000 ass transfer from giant donors Paczynski B., 1965, AcA, 15, 89Paczy´nski B., 1971, ARA&A, 9, 183Paczy´nski B., Sienkiewicz R., 1972, Acta Astronomica, 22,73Passy J.-C., Herwig F., Paxton B., 2012, ApJ, 760, 90Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P.,Timmes F., 2011, ApJs, 192, 3Paxton 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, 4Perez M. S., Blundell K. M., 2010, MNRAS, 408, 2Ritter H., 1988, A&A, 202, 93Savonije G. J., 1978, A&A, 62, 317Soberman G. E., Phinney E. S., van den Heuvel E. P. J.,1997, A&A, 327, 620Webbink R. F., 1985a, Stellar evolution and binariesWebbink R. F., 1985b, Stellar evolution and binaries. p. 39Woods T. E., Ivanova N., 2011, ApJ, 739, L48Woods T. E., Ivanova N., van der Sluys M. V., ChaichenetsS., 2012, ApJ, 744, 12 c (cid:13)000