Adiabatic Mass Loss in Binary Stars - I. Computational Method
The asymptotic response of donor stars in interacting binary systems to very rapid mass loss is characterized by adiabatic expansion throughout their interiors. In this limit, energy generation and heat flow through the stellar interior can be neglected. We model this response by constructing model sequences, beginning with a donor star filling its Roche lobe at an arbitrary point in its evolution, holding its specific entropy and composition profiles fixed as mass is removed from the surface. The stellar interior remains in hydrostatic equilibrium. Luminosity profiles in these adiabatic models of mass-losing stars can be reconstructed from the specific entropy profiles and their gradients. These approximations are validated by comparison with time-dependent binary mass transfer calculations. We describe how adiabatic mass loss sequences can be used to quantify threshold conditions for dynamical time scale mass transfer, and to establish the range of post-common envelope binaries that are allowed energetically.
In dynamical time scale mass transfer, the adiabatic response of the donor star drives it to expand beyond its Roche lobe, leading to runaway mass transfer and the formation of a common envelope with its companion star. For donor stars with surface convection zones of any significant depth, this runaway condition is encountered early in mass transfer, if at all; but for main sequence stars with radiative envelopes, it may be encountered after a prolonged phase of thermal time scale mass transfer, a so-called delayed dynamical instability. We identify the critical binary mass ratio for the onset of dynamical time scale mass transfer as that ratio for which the adiabatic response of the donor star radius to mass loss matches that of its Roche lobe at some point during mass transfer; if the ratio of donor to accretor masses exceeds this critical value, dynamical time scale mass transfer ensues.
In common envelope evolution, the dissipation of orbital energy of the binary provides the energy to eject the common envelope; the energy budget for this process consists essentially of the initial orbital energy of the binary and the initial self-energies of the binary components. We emphasize that, because stellar core and envelope contribute mutually to each other’s gravitational potential energy, proper evaluation of the total energy of a star requires integration over the entire stellar interior, not the ejected envelope alone as commonly assumed. We show that the change in total energy of the donor star, as a function of its remaining mass along an adiabatic mass-loss sequence, can be calculated either by integration over initial and final models, or by a path integral along the mass-loss sequence. That change in total energy of the donor star, combined with the requirement that both remnant donor and its companion star fit within their respective Roche lobes, then circumscribes energetically possible survivors of common envelope evolution.
Mass transfer is the signature feature of close binary evolution. In relatively wide binaries, this mass transfer can occur via capture of some part of the stellar wind of the donor by the accretor, but in close systems, the dominant mode of mass transfer is by Roche lobe overflow (RLOF). Hydrodynamic estimates of Roche lobe overflow rates (Paczyński & Sienkiewicz, 1972; Eggleton, 2006) show that mass transfer rates are typically of order
where is the mass of the donor star, the orbital period of the binary, the effective Roche lobe radius, which is the radius of a sphere of the same volume as the interior of the Roche lobe, the donor’s radius minus its Roche lobe radius , and the donor’s effective polytropic index. This mass transfer rate is exquisitely sensitive to ; so long as mass transfer occurs on a time scale much longer than the dynamical time scale of the binary, .
We can examine the stability of mass transfer by comparing the responses of donor star radius and Roche lobe radius to perturbations in the mass transfer rate. So long as , hydrostatic equilibrium prevails through the interior of the donor star, and dynamical terms become important only in the vicinity of the inner Lagrangian point. Variations in the radius of the donor can be separated into mass- and time-dependent parts:
The first term on the right-hand side represents not only the evolutionary expansion (or contraction) of the donor star in the absence of mass transfer, but also the expansion or contraction due to thermal relaxation if, for example, mass transfer has driven the donor out of thermal equilibrium. The second term on the right-hand side represents the hydrostatic response of the radius to changing mass, with thermal or chemical diffusion and nuclear evolution suppressed; that is, it represents the adiabatic response of the donor star to mass loss. We will henceforth label this logarithmic derivative by the subscript ‘ad’ instead of , and define
The Roche lobe radius of the donor depends on three quantities (assuming a circular orbit) – its mass , the total mass of the binary , and the orbital angular momentum of the binary :
where is the Roche lobe radius in units of the binary separation, and is a function only of the binary mass ratio, or as written here, of the ratio of donor to total mass. We can thus write
In conservative mass transfer, as that term is commonly used, and are conserved quantities, and therefore the latter two terms on the right-hand side of Eq. (5) vanish. However, more generally, a binary may suffer background mass loss through a stellar wind from the companion star, or angular momentum losses through, for example, stellar winds or gravitational wave radiation. Moreover, systemic mass and angular momentum losses may contain contributions arising from a loss of part of the mass transfer stream between binary components; we can call these consequential losses. We therefore decompose the terms in and into terms linearly independent of and given the subscript , and those consequential losses linearly proportional to and given subscript , by analogy with Eq. (2). We then write
By analogy with Eq. (3) above, we define
Now suppose that the donor star exactly fills its Roche lobe. In order for mass transfer to occur, must increase faster (or decrease slower) than would if it were not for mass transfer; that is, the coefficient of in Eq. (10) must be positive. If the mass transfer rate is self-limiting, the stellar radius will parallel the Roche lobe radius as mass transfer proceeds, i.e., . This condition then requires that the second term on the right-hand side of Eq. (10) be negative. Since , it follows that self-limiting (stable) mass transfer can only occur if . In that case, we can solve directly for the mass transfer rate:
Typically, and are of order unity, and mass transfer is driven on a time scale . If the donor remains in thermal equilibrium, mass transfer proceeds on a nuclear time scale. If, on the other hand, the donor has been driven out of thermal equilibrium by mass loss, it typically expands much more rapidly, on its thermal time scale, and accordingly mass transfer proceeds on that thermal time scale. In some circumstances, such as among the shorter-period cataclysmic variables or double-degenerate binaries, the driving time scale may be that characterizing contraction of the Roche lobe by angular momentum losses.
Dynamical time scale mass transfer occurs if . Both right-hand terms in Eq. (10) are then positive, and the donor star radius is then driven by simple hydrostatic relaxation far beyond its Roche lobe radius. In this case, the mass transfer rate can in principle accelerate up to a truly dynamical rate (cf. Eq. (1)). Three-dimensional simulations of this process (Rasio & Livio, 1996; Ricker & Taam, 2006, 2008) show that the donor envelope is tidally disrupted, engulfing the companion star while tidal torques from the lagging, distorted envelope plunge that companion and the core of the donor into orbits deeply embedded in that envelope. The binary enters common envelope evolution (Paczyński, 1976). The common envelope itself lags in rotation behind the embedded cores, which spiral inward as they lose orbital energy and angular momentum to the envelope. Ultimately, either the dissipation of orbital energy suffices to eject the common envelope, revealing the core of the original donor and its companion star111Because the common envelope has a much higher specific entropy than the surface of the accreting companion, it is much hotter and more tenuous at that interface. The companion can accrete very little matter during common envelope evolution – see Hjellming & Taam (1991)., or the cores merge within the common envelope. Examples of the products of common envelope evolution include planetary nebulae with close binary nuclei, cataclysmic variables, and other short-period binaries containing at least one degenerate component. It is important to recognize that common envelope evolution must be very rapid, lest radiative losses quench common envelope ejection, and that energy sources other than orbital energy dissipation are unlikely to contribute significantly to the energy budget of the system (Webbink, 2007). Adiabatic mass loss sequences can therefore define the energy difference between initial and final states. This energy difference must be supplied by orbital energy dissipation222common envelope systems must, of course, also observe angular momentum conservation (cf. Nelemans et al., 2000, 2001; Nelemans & Tout, 2005), but this is almost always a less stringent constraint on the final system..
Useful insights into the behavior of donor stars under adiabatic mass loss can be gleaned from very simplified models. Hjellming & Webbink (1987) investigated the properties of polytropic models with power-law equations of state, including complete polytropes, composite polytropes ( cores with envelopes and equations of state), and condensed polytropes ( envelopes with point mass cores). These models are relevant for some main-sequence stars and giant-branch stars, and illuminate the qualitative (and, within limits, quantitative) behavior of those stars. They were further employed by Soberman et al. (1997) to explore stability criteria for binary mass transfer. However, these models leave many important evolutionary stages of a donor unexplored, and they leave considerable room for improvement even where they may be applicable (Eggleton & Tout, 1989; Han et al., 2002; Podsiadlowski et al., 2002; Chen, & Han, 2008). A more sophisticated approach, using realistic stellar models, is both desirable and possible, as shown by Hjellming (1989a, b). This paper represents a step toward filling this need.
We describe here a method of construction of stellar adiabatic mass loss sequences. In a given sequence, the initial model corresponds to a full stellar model at an arbitrary state in its evolution. As in the pioneering work of Hjellming (1989a, b), mass is removed from the surface while stellar interior at each point relaxes adiabatically to a new state of hydrostatic equilibrium. In Section 2, we describe the modifications to the normal set of stellar structure equations that are needed to construct adiabatic mass loss sequences; how one can recover the internal luminosity profile, and hence thermal relaxation rate, from the adiabatic sequences; and how we have implemented these concepts in a numerical code. Adiabatic mass loss sequences constructed in this way yield , fixing the critical mass ratio for dynamical time scale mass transfer, but they also reveal thresholds for a delayed dynamical instability arising when the donor is stripped down to a nearly isentropic core. These applications are described in Section 3. In Section 4, we discuss the calculation of the total energy of a donor star as it loses mass adiabatically, and the application of those results to determining in what orbit, if any, a binary can survive common envelope evolution. We conclude with a brief appraisal of the relevance of our results to population synthesis studies in Section 5.
2 Adiabatic Mass Loss
Adiabatic mass loss models rest on two basic approximations: that in the limit of rapid mass loss, the response of the donor star is everywhere locally adiabatic; and that hydrostatic equilibrium nevertheless prevails throughout the bulk of the interior of the donor star. The first of these approximations implies not only that heat flow can be ignored, but that (under the Schwarzschild criterion for convection) convective boundaries are preserved with respect to mass. Both the specific entropy profile and the composition profile of the donor interior are therefore fixed along an adiabatic mass loss sequence. The second approximation implies that fluid flow velocities are negligible (highly subsonic) through the interior of the donor. Solutions for the outflow of gas from a lobe-filling star show that this assumption breaks down only in the vicinity of the inner Lagrangian point (see, e.g., Eggleton (2006)).
To illustrate the strength and limitations of the adiabatic approximation, we show in Fig. 1 the specific entropy profiles of the donor star at five epochs during mass transfer for a binary with initial mass = , initial mass ratio (donor/accretor) = , and initial orbital period = . As mass loss progresses from epoch to epoch , the mass transfer rate rises from about (a thermal timescale) to about (a quasi-dynamical timescale). We see that the entropy profile of the donor remains nearly unchanged during mass transfer on quasi-dynamical time scales (from curve to curve ), but that even for the thermal timescale mass transfer, thermal relaxation is largely concentrated in the outermost layers, which are immediately stripped away. Little heat transport occurs in the deep interior. This pattern of thermal relaxation typifies stars with radiative envelopes, with rising entropy driven by energy absorption in the outer envelope (as opacity rises rapidly with decreasing temperature during decompression), accompanied by falling entropy in the deep interior (as thermal energy generation supplants nuclear sources as they are frozen out by decompression). Among stars with deep surface convection zones, the entire convective envelope slowly declines in specific entropy, as radiative energy losses from the surface exceed the input of energy at the base of the convection zone (these being choked off by the local rise in opacity as the underlying radiative zone undergoes decompression). Thus, the adiabatic approximation not only defines the asymptotic response of a donor star in the limit of faster-than-thermal time scale mass transfer, but serves as a good approximation to the deep interior response even during thermal time scale mass transfer.
2.1 The Equations of Stellar Structure
In principle, we can draw the initial model for an adiabatic mass loss sequence from any point in the evolution of the donor star. Following common practice in close binary evolutionary models, we neglect rotational and tidal effects on the structure of that star, and treat it as spherically symmetric. The initial model then satisfies the usual set of four first-order ordinary differential equations of stellar structure, where four structure variables, i.e., pressure , radius , temperature and luminosity , vary with an independent variable, mass ,
Here density , opacity , adiabatic temperature gradient , nuclear energy generation rate , neutrino loss rate , and specific heat at constant pressure are functions of , , and the abundances of various nuclear species . Subscript means entropy. , and are of course the universal gravitational constant, the speed of light and the radiation constant, respectively. is temperature gradient that would prevail in the absence of convection, while the temperature gradient in a convective zone can be written in a common mixing-length treatment (Eggleton, 1983) as
The characteristic convective velocity is the real root of the cubic equation
where is the pressure scale height, the dimensionless mixing-length parameter, and
Four boundary conditions are required to close this set of equations. At the center (),
and at the surface (),
where the latter surface boundary condition approximates a classical gray atmosphere. The four equations are solved for a given distribution of , which is updated according to nuclear reactions and convective mixing. Stellar mass loss can be included by applying the surface boundary conditions at a moving surface such that . Roche lobe overflow in a binary system can likewise be accommodated either by expressing explicitly in terms of the degree of Roche lobe overflow (see, for example, the Appendix to this paper), or (in the approximation that ) by rewriting the structure equations above in terms of rather than as the independent variable, and applying the surface boundary conditions at .
In adiabatic mass loss models, the equations of energy transport and conservation (Eqs. (2.1) and (15), respectively) are supplanted by algebraic conditions fixing the specific entropy and composition profiles:
where the zero subscripts denote the values of these quantities at the same mass coordinate in the initial model of the sequence. Eqs. (12), (13), (23)and (24), together with appropriate boundary conditions (identified below), are sufficient to describe completely the state of the gas at any point in the interior of the donor. The local temperature, , follows implicitly.
To understand the internal relaxation processes set in motion by adiabatic mass loss, we would also like to recover the luminosity profile through the interior of the donor star. This we can accomplish by writing the local energy flux through the stellar interior in terms of the local entropy gradient, which is a known quantity. Given the local state of the gas, the entropy gradient can be written in terms of the pressure and composition gradients:
where we have invoked the definitions of specific heat, , and of the adiabatic temperature gradient, . As overlying layers are stripped away, successively deeper mass layers are brought toward the stellar surface, pressure in those layers decreases, and increases accordingly. The ambient temperature gradient, , is therefore driven toward the adiabatic temperature gradient, , since and are fixed in Equation (25).
Let us define an entropy gradient, , as
Then the ambient dimensionless temperature gradient, , can be written as
The value of can readily be obtained at each layer of an adiabatic mass loss model, and so is a known quantity. If , then ; such a layer is stable against convection. In this case, , and
If , on the other hand, then . Such a layer is convectively unstable and energy is transported by both radiation and convection. The temperature gradient that would prevail in the absence of convection is (cf. Eqs. 16, 17, and 18):
The total luminosity is then
where is calculated from Eqs. (29), (30) and (31). Thus the additional information needed to reconstruct the luminosity profile in an adiabatic model is already carried in the entropy and composition gradients.
It remains to specify the boundary conditions for the adiabatic mass loss models. At the center of the star (), clearly
as before. At the surface of the star, however, there is a fundamental inconsistency between the assumption of adiabatic relaxation (implying no radiative losses) and the finite effective temperature that finite photospheric optical depth inevitably demands. In the outer atmosphere of a star, at small optical depth, radiative relaxation always proceeds more rapidly than dynamical relaxation, so that purely adiabatic relaxation is impossible. However, at large optical depth, the diffusion approximation remains valid, and it fixes the stellar luminosity in terms of the ambient entropy gradient, as we have seen above. We therefore resolve this inconsistency by insisting on continuity of luminosity at the stellar surface. Substituting the blackbody law into the radiative diffusion equation (Eq. (32)), we obtain for the surface boundary condition
This condition is tantamount to one relating the entropy gradient to the entropy itself at the photosphere.
2.2 Numerical Implementation
The numerical code to solve the equations of the adiabatic mass loss model in section 2.1 is written in , based on the stellar evolution code developed by Eggleton (1971, 1972, 1973) and Paxton (2004). This code has been widely used and is updated with the latest input physics, as described by Han et al. (1994, 2003) and Pols et al. (1995, 1998). The essential features of this code have been retained in our adiabatic mass loss model, namely (a) the use of an adaptive, moving, non-Lagrangian mesh; (b) the treatment of both convective and semiconvective mixing as diffusion processes; and (c) the simultaneous, implicit solution of both the stellar structure equations and (in the evolutionary models) the chemical composition equations, including convective mixing.
Our code uses an equation of state, described by Eggleton et al. (1973) and Pols et al. (1995), that includes pressure ionization, Coulomb interactions and dissociation of molecular hydrogen. The opacity tables are taken from Alexander & Ferguson (1994a, b) in the low-temperature limit (where molecules become important), from Itoh et al. (1983) in the high-density, high-temperature regime (where electron conduction becomes important), and from Rogers & Iglesias (1992) elsewhere. The nuclear reaction rates are interpolated in temperature, mostly from Caughlan & Fowler (1988) and Caughlan et al. (1985). Electron screening theory, including strong and intermediate screening, comes from Graboske et al. (1973). The neutrino energy loss rates resulting from pair, photo-, plasma, bremsstrahlung, and recombination neutrino processes originate from Itoh et al. (1996).
The adaptive, non-Lagrangian mesh automatically distributes mesh points at uniform intervals in a function, , of the local structure variables (Eggleton, 1971). The choice of function is in principle arbitrary, subject only to the constraint that it vary monotonically from center to surface of the stellar model. Its introduction involves rewriting the stellar structure equations in terms of , the mesh point number, as the independent (radial) variable, and introducing an additional ‘structure’ equation for the mesh function:
where is the number of mesh points in the model. In our case, we employ 400 mesh points, distributed according to the mesh function
where . We find that his mesh function provides good resolution of structural features of interest throughout our evolutionary and mass loss sequences. Taking the derivative of Eq. (36), and substituting Eqs. (12, 13, and 2.1), with determined from Eq. (27) in lieu of Eq. (14), we can find . Using Eq. (35), we can then cast the structure equations in terms of rather than as the independent (radial) variable.
In practice, we use (a degeneracy parameter related to the electron chemical potential — see Eggleton et al. (1973); Eggleton (2006)), and as our state variables, and , , and as our global variables. The differential equations for , , , and are written in difference form, as usual, but we include also the algebraic equation in the solution matrix as a convenient tool for satisfying the entropy constraint. The two additional differential equations require boundary conditions. At the surface (),
At the center (), we should have , as before, and , but to avoid singular behavior in that limit, the central mesh point is offset:
As intimated above, the construction of adiabatic mass loss sequences involves the interpolation in this adaptive mesh of the specific entropy at the corresponding mass coordinate in the initial model of the sequence. Near the surface of an initial model, typically varies very rapidly with mass, while the mass increment from one mesh point to the next becomes a very small fraction of the stellar mass. Because of the limitations of numerical accuracy (we typically require that the average correction to structure variables in a model deemed ‘converged’ not exceed one part in ), the situation arises frequently that the nominal error in mass coordinate is comparable with (or even exceeds) the mass difference between one mesh point and the next. In this case, the nominal mass coordinates of mesh points near the surface may not be well-ordered, making the accurate interpolation in mass of the specific entropy and its gradient impossible. We avoid this problem altogether by using not the mass coordinate , but the mass complement, , where is the total mass of the star, over the outer half of our mesh. Our convergence criterion, now applied to , ensures a mesh well-ordered in mass. For similar reasons of numerical accuracy, we tabulate , derived from the local stellar luminosity as described above, for the initial model of each adiabatic sequence, and interpolate in this gradient, rather than numerically differentiating and themselves, in reconstructing the internal luminosity profiles of models in our adiabatic mass loss sequences.
As a demonstration of the ability of our adiabatic mass loss code to reproduce faithfully the response of a donor star undergoing rapid mass loss, we have calculated an adiabatic mass loss sequence, beginning at model () of the time-dependent mass transfer calculation shown in Fig. 1. At this point, the mass loss rate from the donor has grown to roughly a thermal time scale rate. In Fig. 2 we show comparisons of the interior density, luminosity, temperature, and radius profiles in the adiabatic models with those in the time-dependent calculation when both models have reached point () of Fig. 1. Only very minor differences appear in layers very close to the surface, but in general these two calculations are in excellent agreement, giving strong support to the use of adiabatic mass loss sequences to model rapid mass transfer.
3 Thresholds for Dynamical Instability
3.1 Physical Significance of
An important application of adiabatic mass loss sequences is the evaluation of threshold conditions for the onset of dynamical time scale mass transfer. As discussed in Section 1, this condition corresponds formally to ; in that case, a sufficiently rapid perturbation to the donor star (i.e., one faster than a thermal time scale) can lead to unbridled growth of the donor radius beyond its Roche lobe (cf. Eq. 10). On the other hand, if the donor will transfer mass on a slower time scale, on a nuclear or angular momentum loss time scale if the donor thermal equilibrium radius-mass relation () is sufficiently steep (), or on a thermal time scale if thermal relaxation tends to drive the donor beyond its Roche lobe ().
Adiabatic responses of stars to mass loss depend strongly on the variation of specific entropy with mass in their outermost layers. This behavior is illustrated in Fig. 3 for lower main sequence stars. As one progresses down the main sequence, the surface convection zone grows deeper, and, as illustrated here, the initial contraction of stellar radius with decreasing mass becomes less and less severe. In the limit of fully convective stars ( in Fig. 3), their behavior approximates the dependence that follows from the homologous expansion of an isentropic gas sphere supported by classical ideal gas pressure. In contrast, stars on the upper main sequence, with nearly completely radiative envelopes, undergo very rapid contraction with adiabatic mass loss, since their high-entropy outer layers occupy a disproportionately large fraction of their stellar volume. Broadly speaking, then, stars with convective envelopes tend toward dynamically unstable mass transfer, while those with radiative envelopes tend toward stability.
The situation in reality is somewhat more complex. As is evident in Fig. 3, the radius-mass exponents , , and are not constant over the course of mass transfer, but vary continuously, in some cases very rapidly. This is true particularly of in the case of a star with a radiative envelope as it first fills its Roche lobe. Typically, the specific entropy rises extremely rapidly through the surface layers of such a star (see Fig. 1), which therefore decrease rapidly in density. Loss of even a small amount of mass leads to rapid reduction in the stellar radius, so is initially very large and positive. Donor stars with radiative envelopes are therefore typically strongly stable against dynamical time scale mass transfer, at least initially. However, as mass loss proceeds, and the surface layers are stripped away, entropy gradients flatten, and the initially rapid stellar contraction moderates considerably. If the deep interior lacks strong entropy gradients (a circumstance typically prevailing during core hydrogen burning), then as the donor star is stripped down to this region it may begin expanding with further mass loss, much as an (isentropic) convective star expands with decreasing mass. Under these circumstances, the donor star may belatedly become unstable to delayed dynamical mass transfer.
3.2 The Delayed Dynamical Instability
The adiabatic response of a terminal main sequence star, shown in Fig. 4, illustrates the nature of this delayed dynamical instability. This star has a deep radiative envelope, and we see that in the adiabatic limit, its initial response is indeed extremely rapid contraction. Recalling that even in thermal time scale mass transfer, little thermal relaxation occurs in the deep interior of a donor star (see section 2, Fig. 1), we can regard the adiabatic mass loss sequence as approximating the minimum radius of the donor star as a function of its mass. We see that, in this example, a binary with an initial mass ratio , assumed to transfer mass conservatively, never reaches this minimum radius in the course of mass transfer; in this case, the mass transfer rate need never exceed a thermal rate in order for the donor star to track its Roche lobe through the mass transfer episode. In contrast, a binary with an initial mass ratio undergoes such rapid Roche lobe contraction (again in the conservative approximation) that, despite the rapid initial adiabatic contraction of the donor, the Roche lobe radius reaches the minimum radius set by the adiabatic mass loss sequence (intersection of the dashed line with the solid line in Fig. 4) when the donor has been reduced to a fraction 0.86 of its initial mass. In this case, the system may transfer mass at a thermal rate up to that point, but the donor is then no longer able to remain within its Roche lobe, and dynamical time scale mass transfer must ensue. This is the delayed dynamical instability. At an initial mass ratio (dash-dotted curve), the Roche lobe radius is tangent to the adiabatic mass-radius relation (at ), a condition which then marks the threshold for delayed dynamical instability, in the limit that the donor star response is truly adiabatic, and that mass transfer is conservative.
Since the delayed dynamical instability follows a period of slower, thermal time scale mass transfer, in which some measure of thermal relaxation inevitably occurs, the question arises, how good is the adiabatic approximation for the threshold of delayed dynamical instability? We can glean at least a partial answer to this question by comparing the adiabatic results in the example above with a family of time-dependent mass transfer models spanning the same range of initial mass ratios, and again assuming conservative mass transfer, but now with the full set of stellar structure equations. We again assume conservative mass transfer, but now require the donor star to exactly fill its Roche lobe throughout mass transfer. The results of this exercise are shown in Fig. 5, showing mass transfer rates as a function of donor mass on the same scale as Fig. 4 to facilitate comparison. In this figure, the onset of delayed dynamical instability is marked by a rapidly accelerating mass transfer rate as a function of decreasing donor star mass. We see that the mass transfer rates in systems with initial mass ratios and both indeed accelerate well beyond the thermal time scale rate, while the transfer rate for the system with initial mass ratio never accelerates beyond a thermal rate. The mass loss rates for the unstable systems and peak at donor star masses and , respectively, somewhat earlier than predicted by the adiabatic calculations seen in Fig. 4, indicating that the modest thermal relaxation that occurs in the time-dependent calculation drives the donor slightly in the direction of instability. In this case, we may overestimate the critical initial mass ratio for delayed dynamical instability by about , which we may then consider a mark of the accuracy and reliability with which we can define this threshold from adiabatic mass loss sequences alone.
3.3 Superadiabatic Expansion
A different sort of complication arises among luminous giant branch stars, in which convection becomes very inefficient in their outer envelopes. The ambient density there falls so low that an increasing fraction of the stellar luminosity must be carried by radiation, notwithstanding the fact that the local medium is unstable to convection. This condition drives the local temperature gradient, , well above the adiabatic gradient, (see subsection 2.1). The entropy gradient thus becomes negative over a significant fraction of the outer envelope; in effect, the outer envelope forms a cool, dense blanket, suppressing the higher-entropy convective envelope contained within it. Removal of this blanketing superadiabatic layer leads to an extremely rapid initial expansion of the star. As seen in Fig. 6, this superadiabatic expansion can be quite dramatic in luminous giant stars, but in fact it is a feature of any donor star with a surface convection zone. Among less luminous stars (generally, those which have not yet reached the giant branch), the mass and radial extent of the superadiabatic layer is small enough that its effect on threshold conditions for sustained dynamical time scale mass transfer are negligible. However, superadiabatic expansion is a pervasive feature of luminous giant branch stars.
The behavior seen in Fig. 6 reveals certain limitations to the utility of spherically-symmetric adiabatic mass loss models in analyzing the stability of mass transfer. In reality, mass is not removed uniformly over the whole surface of the donor, but from an initially small region about the inner Lagrangian point, . Moreover, sufficiently close to the stellar photosphere, radiative relaxation inevitably becomes rapid compared with convective turnover time scales, so that conditions for adiabatic outflow cannot be satisfied. Indeed, for a giant branch star near Eddington luminosities, the thermal time scale of the entire envelope may approach its dynamical time scale. Without full three-dimensional models, including both advective and radiative energy transport, we cannot be certain how much of this superadiabatic expansion is real, and how much an artifact of our simplifying assumptions. Nevertheless, there exists a real physical basis for this behavior, and the adiabatic models show that it is at least energetically possible for superadiabatic expansion to trigger mass transfer when the donor photosphere still lies well within its Roche lobe. In the convective envelope, hydrostatic equilibrium describes an average state of the gas, but in fact it consists of rising and descending flows, which may differ considerably from each other in specific energy or entropy, and approach sonic velocity where the mean entropy gradient is strongly superadiabatic. That kinetic energy is in principal capable of launching gas through the saddle point in gravitational potential at even as the mean stellar photosphere lies well within the donor star Roche lobe. Rising convective flows may be captured by the companion star, with the result that mass transfer rates can fluctuate wildly on convective turnover time scales333Episodic mass transfer driven by superadiabatic expansion was first postulated by Bath (Bath (1972, 1975); see also Papaloizou & Bath (1975))..
3.4 Critical Mass Ratios
In future papers (Ge et al., in preparation), we will present results of a survey of adiabatic responses of stars throughout the Hertzsprung-Russell diagram. Inasmuch as one of the principal objectives of this survey is to establish the conditions under which a lobe-filling donor star of arbitrary evolutionary state is stable or unstable to dynamical time scale mass transfer, we describe here the issues which arise in characterizing our results for that purpose.
We idealize mass transfer between binary components as fully conservative, by which we mean specifically that the total mass of the binary
and orbital angular momentum of the binary
where the mass ratio is defined, as before, as the ratio of donor star mass to accretor mass, we can write the Roche lobe mass-radius relation directly in terms of the binary mass ratio:
Given from the adiabatic mass loss sequence for a donor star of interest, Eq. (44) then implicitly defines a corresponding critical mass ratio, , satisfying the equation
above which a binary containing that donor star is unstable to dynamical time scale mass transfer.
In reality, mass transfer is never truly conservative, and Eq. (44) is at best an approximation to the dependence of Roche lobe radius on mass ratio. In contrast, to the extent that our spherically-symmetric models are appropriate, is intrinsic to the donor star alone, independent of the mass or nature of the companion star. In this sense, it is the more fundamental parameter by which we should characterize our results. However, as noted above, varies continuously over the course of mass loss, and the decrease in with decreasing donor mass is crucial to the existence of a delayed dynamical instability. Conditions for the onset of that instability therefore depend on the extent and nature of mass and angular momentum losses accompanying mass transfer. For the sake of clarity and simplicity, we here and henceforth implicitly assume conservative mass transfer in quantifying thresholds for dynamical time scale mass transfer, and denote by , specifically the threshold initial mass ratio for dynamical instability, and define , the corresponding critical donor mass-radius exponent marking this same threshold.
The evaluation of is complicated by the fact that the distinction between a prompt dynamical instability and a delayed dynamical instability is to some extent artificial. Even if the donor star satisfies the inequality for dynamical instability () upon filling its Roche lobe, that instability will not be manifested until the mass transfer rate exceeds a thermal rate, and becomes too rapid for thermal relaxation to damp its growth. We have therefore adopted an approach to evaluating threshold conditions for dynamical instability that accounts explicitly, if crudely, for the phase of growth in mass transfer rate up to the point at which donor star response can be reasonably approximated as adiabatic.
As a first approximation of the response of the donor star, we suppose that its structure can still be approximated by its adiabatic mass loss sequence, notwithstanding thermal relaxation. This is the same assumption employed above in illustrating the threshold for a delayed dynamical instability, and is subject to the limitations and shortcomings described there. In a prompt dynamical instability, on the other hand, this should be a much better approximation, as the growth time scale for the mass transfer rate is generally much shorter than a thermal time scale.
The mass loss rate from the donor is a function of the structure of the outer envelope of the donor star, and of the degree to which it overfills its Roche lobe,
Details of the calculation of as adopted in this study can be found in the Appendix. At the start of a mass loss sequence, the donor is assumed to fill its Roche lobe exactly:
As mass transfer proceeds, we assume that the Roche lobe follows the radius-mass relation for conservative mass transfer appropriate to its initial mass ratio, . Defining mass fraction
we can write
where clearly . The donor radius follows the adiabatic mass loss sequence, by assumption:
where follows from conservation of total mass. The desired threshold condition for dynamical time scale mass transfer then consists in finding that initial mass fraction (or, equivalently, initial mass ratio ) such that the Roche lobe radius, at its deepest penetration into the envelope of the donor star, suffices to drive mass transfer at just a thermal (Kelvin-Helmholtz) rate. For that rate, we adopt
where the subscripts “1i” refer to the values of the corresponding variables at the beginning of the mass loss sequence.
Let us define as the radius within the donor star at which
The our desired solution simultaneously satisfies the equations
is fixed by the choice of initial model, so these constitute two equations in two unknowns ( and ). The solution value for then corresponds to the minimum initial mass fraction above which the binary suffers dynamical instability; that for specifies the mass fraction at which the instability actually occurs. The corresponding critical initial mass ratio is , and we define
(see Eq. 44) as the initial mass-radius exponent characterizing this critical case. If the binary is subject to a prompt instability, and differ by a negligible amount, but in a delayed dynamical instability the difference provides an approximate measure of the extent of mass transfer prior to the onset of dynamical instability.
In practice, we have found that this approach works quite satisfactorily, and is insensitive to the vagaries of round-off error that can afflict numerical derivatives like and . However, in the case of donors undergoing pronounced superadiabatic expansion, the ambiguities described in Section 3.3 lead us to characterize our results somewhat differently.
We conjecture that, when superadiabatic expansion becomes important, the true threshold for dynamical instability probably lies between two extremes. If we suppose, on the one hand, that a donor only initiates mass transfer when it fills its Roche lobe, as assumed above, then the rapid superadiabatic expansion that follows presents a most favorable circumstance for dynamical instability. The algorithm detailed above then yields estimates for and , the minimum plausible mass ratio, and corresponding initial mass-radius exponent, unstable to dynamical time scale mass transfer. If superadiabatic expansion is severe (), no simultaneous solution exists for Eqs. (55) and (56) in the case of conservative mass transfer. In this case, , i.e., all mass ratios are dynamically unstable, and we set equal to the minimum value of along the mass loss sequence.
At the opposite extreme, as noted in Section 3.3, it is energetically possible for a donor star to initiate mass transfer while its photosphere still lies well within its Roche lobe, given a sufficiently strong perturbation. This circumstance we consider the least prone to dynamical instability, because sustaining mass transfer at a high rate requires rapid contraction of the donor Roche lobe, and therefore a relatively extreme mass ratio. We characterize the initial rapid superadiabatic expansion by , the logarithmic difference between the Roche lobe radius and the stellar radius at the onset of mass transfer. To quantify , we build a pseudo-envelope model for the donor star in question by assigning to its entire convective envelope a uniform specific entropy equal to that at the base of the convective envelope of the true evolutionary model. By construction, this pseudo-model has no superadiabatic surface layer. As seen in Fig. 6 below, the difference between true and pseudo-models is evident only very near the surface. Then, identifying the radius of this pseudo-model with the Roche lobe radius of the donor, we define
where is the radius of the pseudo-model as a function of mass. is an estimate of the fraction by which a donor giant branch star may underfill its Roche lobe and yet initiate tidal mass transfer to its companion. In turn, the same algorithm described above for finding and , applied now to the pseudo-envelope model, yields estimates for and , the maximum plausible mass ratio, and corresponding initial mass-radius exponent, stable against dynamical time scale mass transfer.
Fig. 6 illustrates the application of the algorithms described above to the determination of threshold conditions for dynamical time scale mass transfer. In this example, the donor is a star near the tip of the giant branch. In response to adiabatic mass loss, the radius of such a donor (heavy solid line) undergoes a rapid initial expansion that will trigger dynamical time scale mass transfer, provided that the mass transfer rate becomes large enough to damp thermal relaxation. In this example, the donor has a very extended envelope and short thermal time scale, and so it must overfill its Roche lobe by a significant factor in order to drive the mass transfer rate up to ; the depth in the donor at which this occurs is shown by the thin solid line. The threshold condition for dynamical time scale mass transfer then occurs when the Roche lobe radius just reaches . In the example at hand, however, expansion is extremely rapid that all mass ratios are unstable: and . In this critical case, reaches at , ; that is, dynamical instability is virtually instantaneous.
Also shown in Fig. 6 is a pseudo-envelope model sequence for the same giant branch star (heavy dashed line). The radius in this model at which is shown by the thin dashed line. In the limiting case, it is energetically possible for this star to sustain Roche lobe overflow by adiabatic expansion when it initially underfills its Roche lobe by . At the threshold for this case the Roche lobe radius (dotted line) again touches the thin dashed line but, because the donor star underfills its Roche lobe, more dramatic contraction of the Roche lobe is needed to drive up to . The corresponding tangent point occurs at , , where and , corresponding to and at .
Our conjecture that the threshold initial mass ratio for dynamical time scale mass transfer, , probably lies within the interval , as these limits are estimated above (and likewise for the critical initial mass-radius exponent, ) must remain a matter of speculation at present. It should be noted that the example illustrated in Fig. 6 is a relatively extreme one, and that given the very short thermal time scale of the donor, even dynamically-stable mass transfer on a thermal time scale could trigger common envelope evolution.
4 Common Envelope Evolution
The family of binary stars contains many examples of non-interacting systems containing a compact component (white dwarfs, neutron stars, or black holes) in an orbit so compact that the immediate progenitor of the compact component could not possibly have fit within its Roche lobe. These systems must have originated through common envelope evolution (Paczyński, 1976). Examples include planetary nebulae with close binary nuclei, pre-cataclysmic white dwarf-main sequence star pairs, close double white dwarfs, and many binary pulsars. These systems are presumably the progenitors to interacting binaries containing compact components, either through evolution of a non-degenerate stellar component, or through angular momentum loss from the binary via a magnetic stellar wind or gravitational radiation, for example.
Modeling common envelope evolution poses one of the most demanding challenges in computational stellar astrophysics. The characteristic dimensions of the common envelope may easily exceed those of a compact core embedded within it by factors of or more, corresponding to a ratio of dynamical time scales (envelope versus core) exceeding . Furthermore, the dissipation of energy within the common envelope is governed by ill-understood energy- and angular momentum-transport processes, characterized by a time scale orders of magnitude in excess of the local dynamical time scale, and possibly approaching a stellar thermal time scale. At present, only the initial immersion of the binary components within a common envelope can be modeled with a semblance of realism (e.g., Ricker & Taam (2008)).
All is not lost, however, in projecting the outcome of common envelope evolution, as the entire common envelope process must still satisfy conservation laws for energy and angular momentum. Considering first the energy equation, we can write the initial energy of the binary, , as
where is the initial mass of the donor (with envelope intact), is of course the initial binary separation, and is the initial total energy (gravitational potential plus internal) of the donor star, a negative quantity. Subscripts ‘2i’ refer to the corresponding initial quantities for the companion star. A similar equation holds for the final, post-common envelope state of the binary (subscripts ‘f’). As a consequence of common envelope evolution, the envelope of the donor star, mass , is shed from the binary. For reasons outlined in the Introduction, we can assume that the companion star survives common envelope virtually unchanged in mass or total energy, and so write and . We require that the final energy, , of the binary be no greater than its initial energy, , which yields the inequality
is the difference between final and initial total energies of the donor star, the initial binding energy of the ejected envelope. Clearly , but it is the fact that the initial binding energy of the envelope usually greatly exceeds the initial orbital energy of the binary in magnitude that is most responsible for dramatic decreases in orbital separation resulting from common envelope evolution.
Considering in turn the conservation of angular momentum, we write the initial angular momentum of the binary, , as
where we explicitly allow for initial orbital eccentricity , but neglect the rotational angular momenta of the binary components. The rotational angular momenta are demonstrably small, typically contributing less than 2% to the total angular momentum, provided the mass ratio of the binary is not too extreme. A parallel expression holds for the final angular momentum, . Then requiring, as for energy, that the final angular momentum of the binary be no greater than its initial angular momentum, we find
In general, dissipative processes prior to, or within, common envelope evolution lead us to expect that . Clearly, save for the possibility that , i.e., that spontaneous ejection by the donor star of its envelope is energetically possible, the energy constraint, Eq. (60), always poses a stronger constraint on post-common envelope orbital separation than does Eq. (63).
In the years since common envelope evolution became the generally-accepted model for the origin of compact binaries with an evolved component, various approximations have been put forward for (e.g., Iben & Tutukov, 1984; Webbink, 1984; de Kool, 1990; Iben & Livio, 1993; Webbink, 2007). A more realistic approach (e.g., Han et al., 1994, 1995; Dewi & Tauris, 2000) has been explicitly to integrate over the putative binding energy of the envelope:
In the limit that this integral extends over the entire mass distribution (), this expression gives the total energy of a star of total mass ; but it is not rigorously correct when parsed over a portion of the star (). The reason is that the total gravitational potential energy of a star is a global property of that star. The gravitational potential at any point in the star (in the usual convention, where that potential vanishes at infinity) depends not only on the mass interior to that point, but also on the mass exterior to it, which determines the shape (and integrated depth) of the potential outside that point.
It may help to clarify the issue to consider an analytic toy model. Consider an isentropic polytrope. Such a polytrope has the equation of state
where , the polytropic constant, is a function only of entropy, by assumption. The specific internal energy of a gas obeying this equation of state is
An polytrope has the analytic solution (cf. Chandrasekhar, 1939)
where is a dimensionless radial coordinate (), and is the th (i.e. first) root of the ratio of local mass density to central mass density (). The surface of the polytrope corresponds to the first zero of , i.e., . The scale for the radial coordinate is, for an polytrope,
and the central density is
Note (cf. Eq. 68) that, in the case of an polytrope, the radius is a function only of the polytropic constant , independent of the mass . The local mass coordinate corresponding to is
The putative gravitational potential energy per unit mass at radial coordinate is then
and the internal energy per unit mass
where and are the total mass and radius of the polytrope, respectively. Integrating these energy terms with respect to mass, one finds
Thus, at the surface of the polytrope (), we find
giving for the total energy of the polytrope,
Now suppose that all matter outside is removed adiabatically from the polytrope. This value of corresponds to mass coordinate and radius coordinate in the initial polytrope. Since the specific entropy of the polytrope remains unchanged (by assumption), the polytropic constant , the radial scale factor , and hence the surface radius all remain unchanged by mass loss. We then have for the final total energy of the polytrope
corresponding to a change in total energy
But if we integrate over only the envelope that was removed, we find
that is, the total energy of the initial polytrope! Note that the work done on the polytrope in removing this envelope increment by increment from its surface is
i.e., precisely , the difference between final and initial total energies of the polytrope. Indeed, an adiabatic change of mass is only possible through work done on the polytrope.
The difference between Eqs. (80) and (79) is not particular to this example, but a direct result of hydrostatic readjustment in the mass-losing star. As mass is stripped away from the surface of the star, a complex readjustment occurs in the interior, as the internal energy of the gas does work against the star’s self-gravity to re-establish hydrostatic and virial equilibrium. Since the surface layers are the coolest part of a star, their loss from the star removes little internal energy directly — none at all for a polytropic model, where at the surface. But the loss of matter impacts directly the self-gravitational potential energy of the star. Re-establishment of virial equilibrium then converts internal energy from the interior into gravitational potential energy, to power expansion to a new state of hydrostatic equilibrium.
The adiabatic mass loss sequences constructed according to the precepts of this paper permit us, for the first time, to calculate rigorously the change in total energy between initial and final states of the donor star, integrated over their entire interiors. In common envelope evolution, this is the change in total energy, , that must be supplied by the dissipation of orbital energy and that fixes the largest remnant orbital separation allowed energetically (Eq. 60).
The practical problem of calculating reliably is somewhat more involved. Common envelope evolution is of greatest interest, and most prevalent, among binary stars with giant branch donors. These stars have very compact cores containing a significant fraction of the total stellar mass, and so their cores utterly dominate the evaluation of total energy of these stars. In evaluating the change in total energy on loss of the extended giant envelope, one is then typically seeking the relatively small difference between large numbers. We estimate that our typical relative uncertainty in calculating is of order . In comparison, the ratio of envelope energy to core energy is (very crudely) of order to ; with , and , this ratio is . Numerical limitations to the brute force calculation of the total stellar energy therefore severely limit the utility of that approach in just the circumstances when common envelope evolution is of greatest interest.
Fortunately, the observation above that one can also calculate the change in energy of the donor by tracking the energy content of matter removed from its surface provides a suitable alternative approach to this problem. Therefore, in addition to calculating by direct integration of the total energy of initial and final states (Eq. 61), we can calculate the change in energy of the donor star by integration along the mass loss sequence:
where the integrand is evaluated at the surface (, ) of the model. This expression differs from Eq. (81) only in that we must account for the finite internal energy of matter at the surface of the model; that contribution vanishes at the surface of a classical polytrope. Provided we remove small enough mass increments along an adiabatic mass loss sequence to track the change in radius smoothly, we achieve a typical relative uncertainty in of order . Thus, Eq. (82) affords the more accurate measure of the change of energy of the donor star early in a mass loss sequence, when , while Eq. (61) gives a more accurate measure later in the sequence. In practice, we find the two estimates are in close agreement at the transition, and use a weighted mean of the two to determine best estimates of and .
Given an initial donor of mass and radius , knowledge of allows us to constrain the range of companion masses, , and post-common envelope orbital separations, , that are energetically allowed. We assume that the initial donor fills its Roche lobe: