ANALYTICAL EXPRESSIONS FOR THE ENVELOPE BINDING ENERGY OF GIANTS AS A FUNCTION OF BASIC STELLAR PARAMETERS

We compute our stellar models using a version of the binary stellar-evolution code ev (also known as STARS or TWIN),⁵ developed by Eggleton (Eggleton 1971, 1972; Yakut & Eggleton 2005, and references therein) and updated as described in Pols et al. (1995). In the code, convective mixing is modeled by a diffusion equation with a ratio of mixing length to pressure scale height of l/H_p = 2.0. Convective overshooting in the core is taken into account on the main sequence (MS) for stars with M > 1.2 M_☉ and on the horizontal branch (HB) for stars of all masses. On the MS, we use an overshooting parameter δ_ov  =  0.12 for stars with M  >  2.0 M_☉, which corresponds to an overshooting length of about 0.3 H_p, and 0.0 < δ_ov < 0.12 increasing linearly with mass for stars with 1.2 M_☉ < M < 2.0 M_☉. On the HB, we use δ_ov = 0.12 for all models. The code cannot evolve a model through the helium flash, the violent ignition of helium in a low-mass star with a degenerate core, but it automatically replaces the model at the moment of helium ignition with a tailored model of the same total mass and core mass in which helium has just ignited. For stars above a certain mass (M ≳ 2.1 M_☉ for Z = 0.02), helium ignites non-degenerately and this intervention is not needed.

We define the helium-core mass M_c as the mass coordinate below which the hydrogen abundance does not exceed 10%. We compute the binding energy of the hydrogen-rich envelope, E_bind, by integrating the gravitational and internal energies over the mass coordinate of the model, from the core–envelope boundary to the surface of the star M_s:

$$\begin{equation} E_\mathrm{bind}= \int _{M_\mathrm{c}}^{M_\mathrm{s}} \left(E_\mathrm{int}(m) - \frac{G m}{r(m)}\, \right) {d}m. \end{equation} \tag{ 1 }$$

The term E_int is the internal energy per unit of mass, which contains terms such as the thermal energy of the gas and the radiation energy, but not the recombination energy. More details regarding these assumptions are provided in van der Sluys et al. (2006). In addition, we provide a separate fit for the recombination-energy term, so that it can be included if the user so chooses.

In reality, Equation (1) is a simplification of the actual situation. First, the core–envelope boundary is not uniquely defined, especially for high-mass stars (see the discussion in Section 5.3). Second, formally one would need to compute the total binding energies of the pre-CE donor and post-CE remnant, and subtract the latter from the first. The difference between this exercise and Equation (1) lies in the expansion of the remnant during the CE, which would yield energy that can contribute to the ejection of the envelope (Ge et al. 2010; Deloye & Taam 2010). This may be especially important for donors that fill their Roche lobes of the Hertzsprung gap (HG).

In this study, we compute grids of models with six different metallicities: Z = 10⁻⁴, Z = 0.001, Z = 0.010, Z = 0.015, Z = 0.020, and Z = 0.030. The initial hydrogen and helium abundances (X and Y, respectively) of our model stars are a function of the metallicity Z, and are given by X = 0.76 − 3.0Z and Y = 0.24 + 2.0Z.

For three of these metallicities, Z = 10⁻⁴, Z = 0.020, and Z = 0.030, we calculate additional grids that include mass loss via stellar winds. We use a prescription that was inspired by Reimers (1975):

$$\begin{equation} \dot{M} = - \eta \times \min \left\lbrace \begin{array}{@{}l@{}}3.16 \times 10^{-14}\,M_\odot \,\mathrm{yr}^{-1}\, \left(\!\frac{M}{M_\odot }\!\right) \left(\!\frac{L}{L_\odot }\!\right) \left(\!\frac{E_\mathrm{bind}}{10^{50}\,\mathrm{erg}}\!\right)^{-1} \\ \ms 9.61 \times 10^{-10}\,M_\odot \,\mathrm{yr}^{-1}\, \left(\!\frac{L}{L_\odot }\!\right) \end{array} \right., \end{equation} \tag{ 2 }$$

where we set η = 0.2 for this study. This wind prescription dominates the stellar winds for lower-mass stars on the giant branches, and the upper prescription in Equation (2) is used in most of those models, except where |E_bind| is small, e.g., for stars near the tip of the AGB. For massive, luminous stars, we use the wind prescription by de Jager et al. (1988).

We computed several grids of single-star models for a range of metallicities and for each metallicity, for a range of initial masses. For these grids, we assume that there is no mass loss due to stellar winds. However, for three metallicities (Z = 0.02 (solar) and the two extremes Z = 10⁻⁴ and Z = 0.03), we computed a grid of models where wind mass loss is included, in order to determine how the accuracy of our fits is affected by stellar winds and how we can correct for it (see Sections 3.4 and 4.4).

We selected 73 zero-age main-sequence (ZAMS) masses between 0.8 M_☉ and 100 M_☉ in each grid, distributed uniformly in log M:

$$\begin{equation} \log M_i \simeq -0.038660 + 0.029124 \times i, \end{equation} \tag{ 3 }$$

resulting in M ≈ 0.80, 0.86, 0.91, ..., 93.5, 100.0 M_☉. We ignore stars that have a main-sequence lifetime longer than 13 Gyr for a given metallicity. We use the models with an odd value of i to perform the fits, whereas the models with an even value for i are used as verification, to determine the accuracies of our fits for models that they are independent of.

The full grids of detailed stellar-evolution models are computed for all metallicities considered. For the three stellar-wind grids, only odd-i models are calculated, since this is all that is required to determine the mass-loss correction factor described in Section 3.4. A Hertzsprung–Russell diagram (HRD) for selected stellar models from our grid with Z = 0.02 is shown in Figure 1.

Zoom In Zoom Out Reset image size

After some experimentation we find that the binding energy is best described by the stellar mass and radius, as opposed to other basic parameters like age or core mass. Such a relationship of the binding energy with the mass and radius of a star was also found by Dewi & Tauris (2000). Because the ranges of E_bind, R, and M may span several orders of magnitude, we use the (base-10) logarithm of these variables.

We do our fits separately for each metallicity and do not provide fits that take into account metallicity as a parameter. The reason for this choice is that with each added fitting parameter, subtle features in the models that were well described by the other parameters are washed out when the new dimension is added. However, we compare the changes due to the metallicity in the log E_bind–log R plane for models with a given mass, and find that these changes are relatively smooth, see Figure 3. This suggests that the user could compute the envelope binding energy for the two metallicities that bracket the desired value of Z and interpolate linearly in log E_bind in order to estimate the binding energy for the given star. In practice, however, this grid may be sufficiently dense that the user can take the value of Z that is closest to the desired value.

A CE phase can only be initiated by a star that is at its maximum radius so far in its evolution. Because of this we consider only evolutionary phases in which a star expands (almost) monotonically. This results in the division of the evolutionary tracks into four different groups, and we develop fits for each group separately.

First, we divide our set of stellar-evolution models into a low-mass and a high-mass group. The dividing mass between the two groups lies around 12 M_☉, depends lightly on metallicity, and is given in Table 1.

The low-mass stars in our grid (M ≲ 12 M_☉) evolve through a clear horizontal-branch (HB) phase between the red giant branch (RGB) and AGB. In this phase, the star burns helium in its core and the star's radius is smaller than it was on the RGB, so that this phase can be ignored for CE evolution. In addition, most of these low-mass stars undergo a first dredge-up on the RGB, during which the hydrogen-burning shell meets less-processed material and the star as a whole shrinks slightly. We define the beginning of the RGB as the moment where the hydrogen abundance in the stellar core becomes zero. As a consequence, the envelope binding energy in the HG is also described by the RGB fits for most stars (see Figure 1). We define the end of the RGB as the point of maximum radius before the core helium abundance drops below half of its absolute maximum value over the stellar lifetime. The beginning of the AGB is taken to be the point where the core helium abundance drops below approximately 10% (we vary this value slightly for different stars). We find that this roughly defines the point where the stellar radius begins to increase again after the HB. For all of our models, the "AGB" begins well before the radius surpasses its maximum on the RGB, so that we cover more than the range of radii for which a CE can be initiated. As the tip of the AGB, we take the model where the star reaches its maximum radius in its entire evolution. We find that the accuracies of our fits are greatly improved if we treat the tracks before and after this dredge-up phase separately. Since the dredge-up occurs at different radii for models with different masses, we include a simple function to describe this radius as a function of the stellar mass:

$$\begin{equation} R_\mathrm{LMR} = \sum _{i=0}^4 a_i \left[\log \left(\frac{M}{M_\odot }\right)\right]^i, \end{equation} \tag{ 4 }$$

with the coefficients a_i given in Table 2.

As a result, we split the tracks from the low-mass models into three groups, LMR1 for the low-mass early RGB (before the dredge-up), LMR2 for the low-mass late RGB, and LMA for the low-mass AGB.

For the models with higher masses (M ≳ 12 M_☉) there is no HB and no first dredge-up, and we keep these tracks intact and put them in a fourth group, which we call HM. Hence, these tracks are defined between the moment hydrogen is exhausted in the core near the terminal-age main sequence and the maximum radius ever achieved in the evolution of each star. The lowest-mass model (whose mass is ∼12 M_☉, see Table 1) included in this high-mass group is also included in the low-mass RGB and AGB groups to ensure continuity between the regions (we use the phenomenological parallels in the log E_bind–log R curves, see Figure 2, of lower-mass models to determine where to divide this track between RGB and AGB).

Zoom In Zoom Out Reset image size

We interpolate the evolutionary tracks from each stellar model for each group in the log E_bind–log R plane, in order to give each model the same number (200) of data points, which gives them an equal weight when fitting along the mass axis, and with a constant density along their tracks in this plane, which prevents biasing the fits toward particular regions in these curves. The only exception are the groups LMR1 and LMR2, which are given 200 data points together. Figure 2 shows the model tracks in the LMR, LMA, and HM groups for a selection of the models in our Z = 0.02 grid.

As described above, we choose to find a description for log E_bind as a function of log M and log R. As a general functional form, we select the polynomial of an as yet undefined order in both the instantaneous mass M and radius R. Hence, the general form of our fitting function is

$$\begin{eqnarray} \log \left(\frac{-E_\mathrm{bind}}{\mathrm{erg}}\right) & \approx & E_0 + \Lambda \left(M_0,M\right) \times \ \sum _{m,r} \alpha _{m,r} \left[\log \left(\frac{M}{M_\odot }\right)\right]^m \nonumber \\ && \times \left[\log \left(\frac{R}{R_\odot }\right)\right]^r, \end{eqnarray} \tag{ 5 }$$

with E₀ = 33.29866, α_{m, r} the fitting coefficients, and m and r the integer indices and exponents. The factor Λ depends on the initial (ZAMS) and instantaneous mass and corrects for effects due to wind mass loss in high-mass stars. It is of course equal to unity for stars that do not experience (significant) mass loss and is discussed in more detail in Section 3.4.

Finding a good and accurate prescription for the binding energy then consists of two parts. First, we must determine the orders of each of the polynomials in Equation (5). In fact, we choose to use a liberal definition of the word polynomial, and rather than selecting the orders, we vary the ranges of m and r for which α_{m, r} has non-zero values in order to find the narrowest ranges that still describe E_bind with good accuracy, possibly including negative values for m and r. Second, for each of our choices of the ranges of m and r, we must find the optimal set of fitting coefficients α_{m, r} and quantify the accuracy of the fit, in order to compare to fits that use different polynomials.

For the actual fitting of the data points, we use only the odd-numbered mass models from the grid described in Section 2.2 (the masses with odd values of i in Equation (3)). For each selection of the range of m and r, we fit the data points, i.e., determine the best set of α_{m, r}, using the χ²-minimization-based fitting procedures available in Wolfram Mathematica 8 (Wolfram Research, Inc. 2010). Instead of using Mathematica's "goodness-of-fit parameter," we choose three criteria to compare the accuracies of the different fits in a quantitative way.

First, we use the newly generated fit function to predict the binding energies for all mass models and all data points that are used for the fit (i.e., the odd-numbered stellar-evolution models), and determine for which fraction of these data points the following condition holds:

$$\begin{equation} \left| \frac{E_\mathrm{bind,fit} - E_\mathrm{bind,mdl}}{E_\mathrm{bind,mdl}} \right| < f, \end{equation} \tag{ 6 }$$

where we choose the value f = 0.1 for our criterion (note that here we do not use the logarithm of E_bind). The fraction of data points that fulfills this criterion for f = 0.1 is called $\Delta _{10\%}$. More specifically, we demand that for at least 90% of the data points, the fit has an accuracy of 10% or better; in other words, $\Delta _{10\%} \ge 90\%$.

Second, we use these fits to compute the envelope binding energies for the even-numbered evolution models in our grid, which were not used to create the fits. Again, we use the condition in Equation (6) for f = 0.1 to determine the accuracy of the fit, and demand that the results are similar to those from the first criterion.

Third, we select the simplest polynomial, i.e., the narrowest ranges for m and r in Equation (5), that still pass the first two criteria.

For some groups of data these criteria turn out to be too strict, i.e., there are no reasonable ranges of m and r that could meet the first criterion while still being sufficiently well behaved to also meet the second. For these groups we are satisfied if the fit predicts at least 90% of the data points within 15% of the model value, i.e., we increase f from 0.1 to 0.15 for these groups and demand that $\Delta _{15\%} > 0.9$. For this reason we report both $\Delta _{10\%}$ and $\Delta _{15\%}$ for all groups in the next section.

So far, we have considered conservative stellar-evolution models only. However, we find that when stars lose a significant amount of mass due to stellar winds (i.e., stars with masses ≳ 20–30 M_☉), our fits lose some of their accuracy. Since wind mass-loss rates are uncertain, and because the assumption of evolution without mass loss is unreasonable for the more massive stars in our grids, we introduced the correction factor Λ in Equation (5). The factor is based on a comparison of the fit using conservative models, as discussed in the previous section, to grids of models which experience mass loss as described in Section 2.1. For models in the low-mass grid (LMR/LMA), this correction is not necessary, and we use Λ = 1. For the high-mass models (HM), we find that a correction factor based on the relative amount of mass lost since the ZAMS gives a reasonably good prescription for most models, restoring the accuracy found for conservative models (see Figure 6):

$$\begin{equation} \Lambda _\mathrm{HM}(M_0,M) = 1 + \frac{1}{4} \left(\frac{M_0-M}{M_0} \right)^2, \end{equation} \tag{ 7 }$$

where M₀ is the ZAMS mass and M is the instantaneous mass of the star. Note that Λ_HM is independent of metallicity and that it is valid across the range of metallicities considered here.

The fitting procedure described in this section was carried out for each of the four groups mentioned in Section 3.2 and each of the six metallicities listed in Section 2.2. The results of our fits, their accuracies, and their validity for the case of moderate wind mass loss are discussed in Section 4.

We carried out the fitting of our data points to the fitting function in Equation (5) as described in Section 3. Here, we present the polynomials that best describe the data points, i.e., the ranges of m and r we used in Equation (5), and the values for α_{m, r} that best describe the data for those polynomials. We do this for each of the four groups of our data set, LMR1, LMR2, LMA, and HM, as described in Section 3.2 and defined by Tables 1 and 2 and the equations in Table 7 and the accompanying supplementary material and for each metallicity we consider in this study.

The ranges for m, r are variable, since for some groups of models, the binding energy can be described in fewer terms than for others. The extreme values used are 0 and 20 for m and −5 and 20 for r. We provide the table of coefficients in electronic form, listing m, r, and α_{m, r} for all cases where α_{m, r} is non-zero, as well as the contents of Tables 1 and 2. In addition, we provide Fortran routines which can read these data files and compute the envelope binding energy as a function of metallicity, mass, radius, and evolutionary phase of the star (RGB or AGB) in Loveridge et al. (2011).

Figure 2 shows a comparison between the detailed stellar-evolution models and fits for selected models with a range of masses from our Z = 0.02 grid. Figure 3 shows the dependence of the envelope binding energy on metallicity, by comparing models with different metallicities for three different masses.

Zoom In Zoom Out Reset image size

Table 3 lists the accuracies of our fits for the envelope binding energy for each of the metallicities considered, and for each of the model groups defined in Section 3.2. To express the accuracy of each fit, we list the percentages of the fitted data points that fall within 10% and 15% of the model values ($\Delta _{10\%}$ and $\Delta _{15\%}$).

Table 3. Accuracies of Our Fits, for Each Metallicity Considered and Each Group Defined in Section 3.2

Z	Group	$\Delta _{10\%}$		$\Delta _{15\%}$
		Odd	Even	Odd	Even
10−4	LMR	97.6%	97.7%	98.9%	99.1%
10−4	LMA	86.6%	86.4%	94.8%	94.8%
10−4	HM	94.8%	93.6%	96.4%	94.6%
10−3	LMR	92.0%	90.4%	95.0%	93.6%
10−3	LMA	85.3%	84.3%	95.2%	92.5%
10−3	HM	91.1%	90.0%	95.5%	93.4%
0.01	LMR	95.7%	93.7%	97.1%	96.5%
0.01	LMA	86.3%	87.2%	95.9%	95.0%
0.01	HM	92.4%	90.0%	96.2%	96.1%
0.015	LMR	95.0%	95.9%	96.6%	97.0%
0.015	LMA	91.4%	90.6%	97.8%	96.4%
0.015	HM	90.0%	90.3%	95.5%	95.8%
0.02	LMR	94.3%	92.7%	96.4%	96.4%
0.02	LMA	97.1%	91.9%	99.3%	96.9%
0.02	HM	92.0%	91.7%	96.6%	96.1%
0.03	LMR	95.0%	94.7%	97.1%	96.4%
0.03	LMA	97.0%	90.5%	98.7%	92.5%
0.03	HM	91.1%	90.4%	96.3%	95.8%

Note. For an explanation of the Δs, see Section 3.3 and Equation (6).

Download table as:  ASCIITypeset image

We separately list the accuracies for the odd-numbered stellar-evolution models, which were used to produce the fits, and the even-numbered models, which have masses that lie between those of the odd-numbered models and are used only for verification. Note that the percentages refer to the actual (absolute) value of the binding energy, not its logarithm. The same results are presented in Figure 4, as cumulative histograms of the fraction of data points that lie outside a given accuracy.

Zoom In Zoom Out Reset image size

As these data demonstrate, all of our fits meet the criteria outlined in the previous section, for f = 0.15 in Equation (6). In fact, all but three of the fits also meet the stricter criteria for f = 0.1. For Z ⩽ 0.01, the features in the LMA models vary irregularly from one mass model to the next, so that only for our low-metallicity AGB models, the f = 0.1 criteria are not satisfied. Politano et al. (2010) show that, for the cases where a common envelope (CE) leads to a merger, only 17% of the CEs are initiated on the AGB, against 83% on the RGB, so that the lower accuracy for our LMA results probably has a somewhat smaller impact than suggested by the numbers in the table alone.

In the fits presented so far we consider only the gravitational and internal (e.g., thermal, radiation) energy of the envelope. In what follows we provide separate fit prescriptions for the recombination energy, which may be tapped for unbinding the envelope. This term includes the ionization energy, but not the dissociation energy of molecular hydrogen, which is a factor of ∼5–10 times smaller than the ionization energy. For this, we use the models and methods that were used to fit the actual binding energy above. Because the recombination energy is a simpler function, we do not need to divide the tracks and we can use one prescription per metallicity. As with the binding energy, we provide the coefficients in the supplementary material and include them in the Fortran routines separately. We provide the accuracies for the different metallicity grids in Table 4 and see that they meet the criteria we used for the envelope binding energy. A comparison of our fits with the original models can be found in Figure 5. Since the recombination-energy term is generally an order of magnitude lower than the rest of the envelope binding energy, the accuracy of the prescription will be dominated by the accuracy of the binding energy in most cases. We find that the recombination-energy term is independent of mass loss for the parameter space our prescription applies to, so that we can use Λ = 1 in all cases for the recombination-energy formulae.

Zoom In Zoom Out Reset image size

Table 4. Accuracies of Our Fits for the Recombination Energy Term, for Each Metallicity Considered

Z	Group	$\Delta _{10\%}$		$\Delta _{15\%}$
		Odd	Even	Odd	Even
10−4	All	81.3%	83.5%	91.4%	93.5%
10−3	All	83.7%	84.6%	93.7%	94.4%
0.01	All	91.5%	91.9%	96.2%	96.2%
0.015	All	91.3%	89.5%	96.1%	95.0%
0.02	All	91.5%	91.1%	96.3%	95.8%
0.03	All	91.5%	91.8%	96.1%	96.1%

Note. For an explanation of the Δs, see Section 3.3 and Equation (6).

Download table as:  ASCIITypeset image

The accuracies in the previous section apply to the prediction of the binding energies of stars that have not suffered mass loss due to stellar winds. Here, we compare the outcome of our fits, based on conservative models, to grids of stellar models that undergo wind mass loss as described in Section 2.1. Table 5 lists the accuracies for this comparison for the two extreme metallicities in our grid (Z = 10⁻⁴ and 0.03) as well as "solar" metallicity (Z = 0.02), expressed as $\Delta _{10\%}$ and $\Delta _{15\%}$ (see Section 3.3). We see that when the factor Λ in Equation (5) is ignored, the LMR grid is hardly affected compared to Table 3, whereas the accuracies in the LMA and especially the HM grids suffer appreciably. As described in Section 3.4, the correction factor Λ defined in Equation (7) improves the accuracies for the high-mass models significantly. This improvement is also shown in Figure 6. For the LMA group, we found no simple correction factor that is independent of metallicity. Since the drop in accuracy due to winds is not as large as for the LM models, we decided to leave these fits uncorrected.

Zoom In Zoom Out Reset image size

Table 5. Accuracies of Our Fits as in Table 3, but Applied to Models with Stellar Wind, with Λ = 1 and Λ = Λ_HM (Equation (7))

Group	Z	Λ = 1		Λ = ΛHM
		$\Delta _{10\%}$	$\Delta _{15\%}$	$\Delta _{10\%}$	$\Delta _{15\%}$
LMR	10−4	98.8%	99.3%
	0.02	94.3%	96.4%
	0.03	93.7%	96.4%
LMA	10−4	77.3%	87.4%
	0.02	75.5%	87.5%
	0.03	68.3%	76.6%
HM	10−4	69.7%	86.4%	92.1%	96.1%
	0.02	54.7%	63.1%	87.2%	89.5%
	0.03	60.1%	70.7%	87.0%	89.8%

Download table as:  ASCIITypeset image

The advantage of the Λ factor is that our prescription for the binding energy can be used with arbitrary mass-loss prescriptions and still provide the accuracies listed in Table 5. However, we find that for very strong mass loss, even this correction factor can no longer describe our models accurately. Our results show that when a star loses more than ∼20% of its ZAMS mass, our fits become increasingly inaccurate. For each metallicity, this implies an upper limit in mass for which these fits can be used. We did not include the results of more massive models in Table 5 and discuss these limitations in more detail in Section 4.5.

Here we note the limitations on the applicability of the fit functions. Outside these limits we do not provide fits and hence cannot estimate accuracies.

In the absence of mass loss, Table 3 displays the accuracies of our fits for any star with a mass 0.8 M_☉ ⩽ M ⩽ 100 M_☉, the full range of masses we computed. However, when mass loss is present, we note that the parameter Λ_HM is only effective when stars have lost less than ∼20% of their ZAMS masses. For the wind prescriptions used in this study, this limit can be translated to an upper mass limit of approximately 50 M_☉ for Z = 0.02 and Z = 0.03, and about 75 M_☉ for Z = 10⁻⁴.

In terms of radius, we have considered each model from the end of the MS to the star's maximum radius on the AGB, excepting only phases of radial shrinkage—most notably on the horizontal branch, so that any radius that may correspond to the initiation of a CE in a binary is covered in our study.

The values of the metallicity considered here are discrete, and we do not provide accuracies for Z < 10⁻⁴ and Z > 0.03. Analogous to metallicities that lie between those used in this study, for metallicities outside this range the reader can either choose to use the value for Z that lies closest to the desired value, or to use an extrapolation scheme.

Figures 2 and 3 show examples of the envelope binding energies of a selection of detailed stellar-evolution models (black dotted lines) compared with the results of our fits (colored lines). These figures indicate that the fitting function we defined in Equation (5) describes the binding energy fairly well. A full picture of the accuracy of our fits, when compared to our models, can be read from Table 3 and Figure 4, and from Table 5 and Figure 6, for models with conservative evolution and models that undergo wind mass loss, respectively.

We want to emphasize here the limitation of our results for stars that undergo substantial wind mass loss. When stars lose more than ∼20% of their ZAMS mass, the binding energy is no longer easily related to that of a star of the same mass but without a mass-loss history. This is the result of an interplay between the different effects of wind mass loss on both the mass and structure of the envelope. Strong mass loss ultimately reduces the envelope mass, which results in a less tightly bound envelope compared to a conservative star of the same total mass and radius. However, our models indicate that the envelope structure is also affected by mass loss, in a way that counteracts the first effect. We find that this results in a more tightly bound envelope than that of a conservative star, unless the mass loss becomes very significant. Hence, our compensation factor Λ_HM in Equation (7), which is greater than unity, gives reasonably good results for small amounts of mass loss, but becomes gradually less accurate when more than ∼20% of the ZAMS mass of the star is lost. We note that for our choice of wind mass-loss prescription, this happens for stars with masses larger than ≳ 50 M_☉ for Z = 0.02, 0.03 and for M ≳ 75 M_☉ for Z = 10⁻⁴, hence the most massive stars in our grids. These upper limits will be decreased when stronger winds are used in the stellar-evolution models. However, as we shall discuss in Section 5.3, the dominant source of uncertainty for these high-mass stars is the definition of the core–envelope boundary.

The choice of the mixing length relative to the pressure scale height will directly influence the radius of a model star with a deep convective envelope. Dewi & Tauris (2000) looked at the effect of this parameter and found a weak dependence for the binding energy, with a higher mixing-length ratio translating into a slightly lower binding energy on the AGB. We compared our fits based on our standard assumption l/H_p = 2.0 to a few test models which were computed with l/H_p = 1.5 and found a similar effect: a drop in accuracy almost exclusive to the AGB, resulting in ∼15%–20% more data points outside the 10% accuracy range.

The binding energy of a stellar envelope is computed by integrating the different energy sources (thermal, gravitational, etc.) from the core–envelope boundary to the surface of the star (Equation (1)). Hence, the value of the binding energy will depend on the definition of the core mass M_c, which we chose to define as the mass coordinate at which the hydrogen abundance reaches 10%. Tauris & Dewi (2001) consider several reasonable choices for this definition, one of which is identical to the definition used in our study (the central region of the star with X < 0.1). If we follow these authors and discard the extreme cases (first and last row) in their Table 1 (their λ_b is almost identical to our λ_env defined in Equation (8)), we see that while for low-mass stars different definitions result in similar values for M_c and hence the binding energy, for more massive stars the results can vary quite appreciably. For a 20 M_☉ star, the variation in E_bind can even be more than an order of magnitude. They also note that "our" choice of the definition of the core gives the lowest value for λ_env, hence the highest value for the binding energy (and so the most tightly bound envelope), among the three options considered.

Naively, we would therefore expect that the definition of the core mass would have little or no influence on our results for low-mass stars, with masses ≲ 7–10 M_☉. Since these stars provide the majority of CE donors in binaries (van der Sluys et al. 2010), most instances of CEs in nature would be little affected, as would the results for low-mass compact binaries, such as double white dwarfs, CVs, and low-mass X-ray binaries. On the other hand, for high-mass stars, we expected the binding energy to depend strongly on how the core–envelope boundary is defined.

In order to test this, we compute the binding energy for models with two alternative definitions for this boundary, and determine how accurately the binding energies are described by our original fit functions. We use the odd-numbered stellar-evolution models (see Section 2.2) for this test, and solar metallicity. As previously stated, our original definition of the core–envelope boundary is the mass coordinate where the hydrogen abundance equals 10% (X = 0.10). Here, we include the definitions X = 0.15 and X = 0.20 as well.

Our results in Table 6 show, contrary to expectation, that the application of our results to models with a different definition of the core–envelope boundary yields significantly lower accuracies than before, even though the formulae remain reasonable approximations at 25% accuracy. This effect is most apparent for low-mass stars (M ≲ 11 M_☉), even though for these stars the core (and hence envelope) masses do not change much. The explanation here is that even though the envelope mass changes by only a small fraction, this is the most tightly bound part of the envelope and hence the effect on the total binding energy of the envelope is relatively large. For high-mass stars, the lower accuracy is expected, because a small change in the definition of the core–envelope boundary is known to result in larger differences in core mass.

Table 6. Accuracies of Our Fits as in Table 3, but Applied to Models with Different Definitions of the Core–Envelope Boundary

Group	X = 0.15			X = 0.20
	$\Delta _{10\%}$	$\Delta _{15\%}$	$\Delta _{25\%}$	$\Delta _{10\%}$	$\Delta _{15\%}$	$\Delta _{25\%}$
LMR	65.6%	76.7%	87.8%	54.8%	60.9%	70.9%
LMA	79.5%	93.7%	99.6%	71.3%	80.4%	92.6%
HM	82.6%	85.9%	88.0%	71.8%	74.9%	77.5%

Note. All data are for Z = 0.02.

Download table as:  ASCIITypeset image

In reality, the exact remnant mass after the CE will depend on the response of both the donor star and the orbit to the spiral-in, and hence also on the properties of the binary companion (Deloye & Taam 2010; Ivanova 2011). Therefore, for massive stars no unique core mass can be defined and detailed stellar-structure and binary-evolution models are needed to determine a self-consistent outcome of a CE. The uncertainty thus introduced is likely to be much larger than the uncertainties due to the definition of the core mass and those due to wind mass loss, as described in Section 5.1.

Finally, we would like to remark that our choice of the core–envelope-boundary definition is identical to that used for the fits that resulted in the stellar-evolution prescriptions in the SSE and BSE codes (Hurley et al. 2000, 2002), which form the basis of a number of population-synthesis codes, such as StarTrack (Belczynski et al. 2008) and the latest version of SeBa (Portegies Zwart & Verbunt 1996; Nelemans et al. 2001; S. Toonen et al. 2011, in preparation). Hence, even if there is no unique definition of the core mass, our fits will provide consistent prescriptions when implemented in the most commonly used population-synthesis codes.

The purpose of this study is to provide fits of the binding energy of stellar envelopes, based on basic stellar parameters, so that they can be used to treat CEs in population-synthesis simulations, where detailed stellar-structure models are not available. So far, such codes have often used the so-called envelope-structure parameter, λ_env, defined by

$$\begin{equation} E_\mathrm{bind}= - \frac{G M M_\mathrm{env}}{R \lambda _\mathrm{env}}, \end{equation} \tag{ 8 }$$

to compute the binding energy (Webbink 1984; de Kool et al. 1987), instead of Equation (1). In such a case, one needs to assume a value for λ_env. For example, Nelemans et al. (2000) and Hurley et al. (2002) use λ_env = 0.5, while Belczynski et al. (2008) choose α_CEλ_env = 0.5 and 1.0. Since λ_env has been used extensively, it is useful to give an indication here of our results expressed in terms of this parameter, and to compare them to the choices above. In Figure 7, we show the values of λ_env for three different stellar-evolution models and the results of our fits, converted to λ_env. This figure alone indicates that λ_env is far from constant and will vary for stars of different masses and different evolutionary phases. Dewi & Tauris (2000) and van der Sluys et al. (2006) give more extensive examples of the variation of this parameter. In addition, van der Sluys et al. (2010) show that, for Z = 0.02, the assumption of α_CEλ_env = 1.0 implicitly assumes that α_CE > 1 for roughly 60% of the CEs in a typical stellar population, while for some cases α_CE > 10 is implied. These examples strongly indicate that a better prescription for the envelope binding energy is needed.

Zoom In Zoom Out Reset image size

The reason why we chose to develop an expression for the envelope binding energy, and not λ_env, may now become clear. The parameter was introduced in order to estimate the binding energy when it could not be computed otherwise. From now on, a reasonably accurate prescription for E_bind is available and there is no longer any need to use the parameter λ_env.

Previous studies that provided a simple prescription for the binding energy have been limited to tabulations of λ_env (Dewi & Tauris 2000) for stars of different masses and different evolutionary stages. However, while our paper was in preparation, a paper was published by Xu & Li (2010), largely with the same idea as our study: to create simple and fast prescriptions that can provide the envelope binding energy in order to treat CEs in population-synthesis codes.

The main points where our study improves on the latter paper concern the range of applicability. First, whereas Xu & Li (2010) provide fits for a selection of 14 discreet masses between 1 M_☉ and 20 M_☉, we use the stellar mass as a continuous parameter so that it can take any value between 0.8 M_☉ and 100 M_☉. Hence, the user does not have to interpolate between masses, which may not be a trivial process. For example, while stars of 1.0 M_☉ and 1.5 M_☉ reach a radius of 100 R_☉ on the RGB, a 2.0 M_☉ star does not. It may therefore be impossible to obtain the binding energy for a 1.5 M_☉ star at that radius on the RGB by interpolating the results of the 1.0 M_☉ and 2.0 M_☉ models in a discrete grid. Second, they provide two choices of metallicity (Z = 0.001 and Z = 0.02), whereas we include six different, albeit still discrete, values, ranging from Z = 10⁻⁴ to Z = 0.03. We believe that our grid of metallicities is sufficiently dense that the nearest value can be used without too much loss in accuracy. Third, whereas their study includes a fixed choice of stellar winds, our fits allow for the user's choice of wind prescription (although the accuracy we present is valid only within certain limits, as discussed in Section 5.1).

In addition to having wider applicability, our paper presents a more extensive error analysis. As we have shown in Section 5.1, our accuracies are dominated by the uncertainty in wind mass loss for massive stars, by the uncertainty in the separation between core and envelope at the end of a CE, or even by the intrinsic problems in determining this separation a priori for stars with masses of ∼10–20 M_☉ and up (Section 5.3). These systematic errors are discussed by Xu & Li (2010) as well. However, for the majority of CEs the donor mass is ≲ 10 M_☉, so that these external uncertainties are less important. For these cases the uncertainty in the fits will dominate, which we present explicitly in this paper.

On the other hand, the study by Xu & Li (2010) has two advantages over ours. First, their work covers the HB for all masses they consider. While we include all stages of stellar evolution where a CE could be initiated, knowledge of the envelope binding energy on the HB for low-mass stars, which cannot initiate a CE, may be useful for, e.g., collisions. Second, their prescriptions are simpler than ours, making them easier to implement.