A new synthetic model for asymptotic giant branch stars

Abstract

We present a synthetic model for thermally pulsing asymptotic giant branch (TPAGB) evolution constructed by fitting expressions to full evolutionary models in the metallicity range 0.0001 ≤Z≤ 0.02. Our model includes parametrizations of third dredge-up and hot-bottom burning with mass and metallicity. The Large Magellanic Cloud and Small Magellanic Cloud carbon star luminosity functions are used to calibrate third dredge-up. We calculate yields appropriate for galactic chemical evolution models for ¹H, ⁴He, ¹²C, ¹³C, ¹⁴N, ¹⁵N, ¹⁶O and ¹⁷O. The initial–final mass relation is examined for our stars and found to fit to within 0.1 M_⊙ of the observations. We also reproduce well the white dwarf mass function for masses above about 0.58 M_⊙. The new model is to be implemented in a rapid binary star evolution code.

1 INTRODUCTION

Full stellar evolution models of thermally pulsing asymptotic giant branch (TPAGB) stars are difficult and extremely time-consuming to make, and often suffer numerical failure. For this reason, synthetic models based on full stellar evolution models but with the complicated physics replaced by simple expressions are a useful approximation. Such models are well suited to exploring regions of large, multidimensional parameter spaces which would take years to explore with full stellar evolution models. Here we build on the models of Groenewegen & de Jong (1993) and Wagenhuber & Groenewegen (1998) with accurate new parametrizations of third dredge-up and hot-bottom burning (HBB), as well as new fits to stellar structure from our full stellar evolution models. Stellar evolution up to the TPAGB is handled by coupling the new code with the rapid evolution code of Hurley, Pols & Tout (2000) and Hurley, Tout & Pols (2002, H02) which is designed for single, binary and star cluster evolution. The eventual aim of this work is to implement TPAGB evolution and especially nucleosynthesis in the binary star model and, while this is yet some way off, we have made significant progress in constructing a capable single star model. The calculation of binary star yields will result.

Table 1, mostly taken from the review of Henry (2004), summarizes attempts at calculating AGB yields to date. It does not include important works on AGB evolution that do not specifically involve yield calculation, such as the series of works by Boothroyd & Sackmann (e.g. Boothroyd & Sackmann 1988) who investigate Li, Be, B and ¹²C/¹³C ratios in particular, Straniero et al. (1997) who tackle the carbon star formation problem in low-mass (1 ≤M/M_⊙≤ 3) solar-metallicity stars, Lattanzio and collaborators' contribution to HBB (Lattanzio et al. 1997), carbon star formation (Lattanzio 1989) and degenerate pulses (Frost, Lattanzio & Wood 1998), Mowlavi and collaborators' works on ²⁶Al (Mowlavi & Meynet 2000), fluorine (Mowlavi, Jorissen & Arnould 1998) and sodium (Mowlavi 1999), and models by Herwig (2000) detailing convective overshooting and its consequences. There are also countless papers dealing specifically with the s-process in TPAGB stars (e.g. Busso et al. 2001).

Synthetic models have been constructed in the past (e.g. Renzini & Voli 1981; Iben & Renzini 1983) and are still very much in use (Mouhcine & Lançon 2003). Hybrid models which combine aspects of synthetic and full evolution have also been constructed (Marigo, Bressan & Chiosi 1998; Marigo 1999, 2001). While we lean heavily on the work of Groenewegen & de Jong (1993) and Wagenhuber & Groenewegen (1998), our new model has some marked differences in its treatment of dredge-up and HBB. Most previous models assume a constant dredge-up parameter and minimum core mass for dredge-up. We include expressions fitted to our full stellar evolution models (Karakas, Lattanzio & Pols 2002) for these and the stellar structure (luminosity, radius, core mass, etc.). HBB is included in a similar way to Groenewegen & de Jong (1993) but with calibration of the free parameters to our full evolution models.

While the use of a purely synthetic code is inferior in accuracy or detail to full stellar modelling [or the envelope burning technique of Marigo 1999a (M99)], it is the only way to explore a large parameter space such as a full study of binary stars. The single star space could conceivably consist of the mass M, metallicity Z and a few free parameters such as the minimum mass for dredge-up, dredge-up efficiency and perhaps the mass-loss rate. For binaries the problem is worse, there are the two masses, separation and eccentricity, and in addition free parameters associated with uncertain details of common-envelope evolution, mass transfer and accretion and enhanced wind loss owing to the companion. A parameter space takes a time δt×n^N to explore, where δt is the average model time, n is the number of grid-points per free parameter and N is the number of free parameters, so a fast model is desired. Our model has an average execution time (including features not described here, such as a companion star and associated mass transfer, NeNa and MgAl cycles, supernovae and novae) of 0.05 s on a 2.1-GHz AMD Athlon CPU (0.64 s on a Pentium Pro 200-MHz CPU), so for a 10⁶-point grid we have a total execution time of nearly 14 h. An increase of δt to 1 min (which would be an extremely fast full stellar evolution model) increases the total execution time to just less than 2 years.

The drawback of a fast code is a loss of accuracy and, while we try to fit to our full evolution models as well as we can, it is impossible to fit perfectly. We have to interpolate and sometimes extrapolate into regions of parameter space where we cannot be sure that we get the correct answer. Since we aim to investigate binary stars with this code, we put up with these limitations and keep in mind that detailed models of binary stars may differ. There are some things that a synthetic model avoids, such as numerical breakdown which can occur with a full evolution model. The synthetic model is no worse than our detailed model for star-to-star analyses.

Table 2 lists some common variables used in this paper. Section 2 describes our full stellar evolution models. Section 3 contains the (gory) details of our synthetic model. Detailed calibration and analysis of the variation in the free parameters introduced in the HBB model are made in Section 4. Calibration of third dredge-up using carbon star luminosity functions is described in Section 5.1. The initial–final mass relation and white dwarf mass functions derived from the model are compared with observations in Sections 5.4 and 5.5. Yields from single stars are calculated in Section 6 with comparison to our full stellar evolution model yields and the yields of van den Hoek & Groenewegen (1997) and Marigo (2001). The appendix (in the online version of the article only) contains the coefficients for the fitting formulae, yield tables and synthetic-detailed model composition comparisons.

2 FULL EVOLUTIONARY MODELS

Our full stellar evolution models used are those described in Karakas et al. (2002) (hereafter K02). They were constructed with the Monash version of the Mt Stromlo Stellar Evolution Code (Wood & Zarro 1981; Frost 1997) updated to use the OPAL opacity tables of Iglesias & Rogers (1996). The thermally pulsing phase of the AGB is covered by the models until mass loss makes convergence impossible. Mass loss is parametrized on the red giant branch using the Kudritzki & Reimers (1978) formula with η= 0.4 and on the AGB using the prescription of (Vassiliadis & Wood 1993, VW93). The mixing length parameter α is set to 1.75 and convective overshooting is not included in the models.

We have high-resolution evolution (taken every 100 models from the stellar structure code) and nucleosynthesis model data for Z= 0.02, M_i= 3, 4, 5, 6, 6.5 M_⊙, Z= 0.008, M_i= 4, 5, 6 M_⊙, Z= 0.004, M_i= 4, 5, 6 M_⊙ and Z= 0.0001, M_i= 1.25, 2, 2.25 M_⊙, and lower resolution data (every 1000 stellar structure models) for Z= 0.02, M_i= 1, 1.25, 1.9, 2.5, 3.5 M_⊙, Z= 0.008, M_i= 1, 1.9, 2.5, 3, 3.5 M_⊙, Z= 0.004, M_i= 1, 2, 2.5, 3, 3.5 M_⊙ and Z= 0.0001, M_i= 1.75 M_⊙, where M_i is the initial (zero-age main-sequence, ZAMS) mass of the star. There are typically a few thousand evolutionary models per interpulse period.

3 OUR SYNTHETIC MODEL

Stellar evolution from the ZAMS up to the thermally pulsing AGB is already dealt with in the rapid evolution code (Hurley, Pols & Tout 2000). The main sequence, giant branch evolution and early AGB (EAGB) abundance changes can be represented by simple formulae dealing with first and second dredge-up. All abundances are mass fractions. Coefficients for fits are in Appendix A unless otherwise stated. The fits are made using a Levenberg–Marquardt gradient descent iterative ��²-minimization code (Press et al. 1992).

The state of the star at the beginning of the TPAGB is known from the fits given in Sections 3.2 and 3.3. The star is evolved forward in time pulse by pulse (see Fig. 1). Between pulses the star loses mass at a rate from the envelope in a wind. The core grows owing to hydrogen burning and the envelope material may experience HBB. At every time-step (usually coincident with a thermal pulse) the HBB algorithm is activated and, if the time since the previous pulse exceeds the interpulse period, is immediately followed by third dredge-up. The change in core mass due to nuclear burning and dredge-up combined with the effect of wind loss (see Section 3.6) determines the time evolution of the star.

The rapid stellar evolution code defines a time-step δt such that certain variables (e.g. radius or angular momentum) may not change by a large amount during that time-step. We add an additional constraint that the time-step must be at most the interpulse period. Because we expect the time-step to be the same as the interpulse period most of the time (especially for detached binaries), luminosity variations owing to the pulse cycle are averaged out, although luminosity changes owing to an increase in core mass are followed if the time-step is small enough.

3.1 From the ZAMS to EAGB

Stellar evolution on the main sequence does not affect the surface abundances, except in some rare cases not considered here (e.g. when the lifetime of the star is long enough that diffusion becomes an important transport mechanism). Only during hydrogen shell burning, when the star is a red giant, does the convective envelope reach down into burned material and mix the products of nuclear processing. This mixing event changes the surface abundances and is known as the first dredge-up. A similar process takes place during the early AGB when second dredge-up follows the transition to helium shell burning.

3.1.1 Initial abundances

To some extent the synthetic model is independent of the initial abundances because HBB is dealt with by solution of the appropriate differential equations. Abundance changes at first, second and third dredge-up are independent of modest changes in the initial abundances because, to first order, the stellar structure does not depend on the abundance mix at a given Z.

3.1.2 First dredge-up

3.1.3 Second dredge-up

3.2 The start of the TPAGB

On the TPAGB a star consists of a hydrogen-rich convective envelope and a hydrogen-deficient core. The core is defined here as the mass inside which the hydrogen abundance is less than half the surface abundance [≲35 per cent by mass, as in Wagenhuber & Groenewegen (1998).] The envelope is assumed to be fully convective.

At the beginning of the TPAGB we calculate the core mass (Section 3.2.1) as a function of initial TPAGB mass M_1TP and metallicity Z. The luminosity, radius and initial interpulse period are then calculated.

3.2.1 Mass and core mass

Our initial core masses are compared with the models of the Padova group (Girardi et al. 2000, hereafter G00, as used by Marigo 2001) in Fig. 2. The G00 convective overshooting models (which are evolved without mass loss) all have a dip in M_c,1TP at around 2 M_⊙ which is not so pronounced in any of the non-overshooting models. The G00 models have a lower M_c,1TP by up to a few hundredths of a solar mass for M_1TP < 3 M_⊙.

3.3 Evolution on the TPAGB as a function of M, M_c and Z

3.3.1 Interpulse period

3.3.2 Luminosity

3.3.3 Radius

3.3.4 Temperature at the base of the convective envelope

The temperature at the base of the convective envelope T_bce is the critical factor which governs the rate of HBB (see also Section 4.3.1). If the temperature is sufficiently high (T_bce > T_HBBmin≈ 10^7.5 K) it is possible that hydrogen burning occurs in the convective envelope of the star, altering the surface abundances. From our full stellar evolution models we see that the temperature rises quickly in these hot-bottom envelopes and then stays at a roughly constant value until the envelope mass becomes small. In order to model the HBB, a fit to the temperature in the burning zone is needed.

3.3.5 Density at the base of the convective envelope

3.4 Third dredge-up

The dredged-up material is instantaneously mixed with the convective envelope of the star. We note that there is a possibility of a degenerate thermal pulse in some stars (Frost et al. 1998); however, the effect is to increase the amount of ¹²C dredged up by a factor of 4, making one degenerate pulse equivalent to four normal pulses. Frost et al. (1998) report degenerate thermal pulses in a M= 5 M_⊙, Z= 0.004 star which would also undergo many dozens of non-degenerate third dredge-up events, so the effect of one or two degenerate pulses is small compared with the effect of non-degenerate pulses. We neglect the phenomenon.

3.4.1 Lambda parametrization and minimum mass for dredge-up

In Fig. 5 (left-hand panel) we show the temperature at the base of the convective envelope from our full stellar evolution models. Note that use of the Marigo, Girardi & Bressan (1999a) prescription for dredge-up above log₁₀T_bce= 6.4 would lead to dredge-up in all our stars. Also in Fig. 5 is a comparison of our M^min_c with Marigo et al. (1999a)'s prescription (with log₁₀T_bce= 6.4) for solar metallicity. Even after calibration of our synthetic model by comparison to carbon star luminosity functions (see Section 5.1), which leads to a reduction in M^min_c, our synthetic models still have a slightly higher M^min_c than Marigo et al. (1999a).

3.4.2 Intershell abundances

Our full stellar evolution models do not obtain high values of intershell ¹⁶O such as the 2 per cent reported by Boothroyd & Sackmann (1988). There is some debate on the exact composition in the intershell region. The inclusion of diffusive overshooting (Herwig 2000) increases the abundance of ¹²C and ¹⁶O at the expense of ⁴He.

3.4.3 Dredge-up of the hydrogen-burning shell

Others (notably GdJ93, equation 35) include nuclear burning of third dredge-up material as well as envelope burning (Section 3.7 below). The reason for this is that material brought up by third dredge-up is preferentially exposed to high temperatures at the base of the convective envelope. Our full stellar evolution models do not show this phenomenon although, in our low-metallicity models, there is dredge-up of ¹³C and ¹⁴N that leads to a similar effect which cannot possibly result from helium burning. These isotopes are enhanced in the envelope by dredge-up of material previously in the hydrogen-burning shell (burnt during the interpulse period) but not mixed into the intershell convective zone.

We account for this by burning a fraction f_DUP of ΔM_dredge for a fraction of the interpulse period f_burn,DUP and at the temperature and density at the base of the envelope (extrapolated from equations 37 and 41). The hydrogen abundance of the material to be burned is set to the envelope hydrogen abundance (even though it may be somewhat lower due to interpulse hydrogen shell burning). Because f_DUP and f_burn,DUP are fitted to the full evolution models, any problems are circumvented by the calibration. The burning algorithm is described in Section 3.7.1 and is the same algorithm as used for HBB. The hydrogen-burned material is immediately mixed with both the helium-burned intershell material and the whole convective envelope to give the post-third dredge-up envelope abundances.

Note that when normal HBB occurs it is the dominant burning mechanism. At metallicity greater than 0.004 the change of ¹³C and ¹⁴N in the envelope due to the dredge-up of the H-burning shell is negligible compared with the abundance of ¹³C and ¹⁴N already in the envelope. The model used here is only approximate (and does not reflect the actual hydrogen-burning process), and we hope to improve upon it soon when more low-Z full evolution models become available.

3.5 Core growth during the interpulse period

3.6 Wind loss prescription

Note that the radius used in the above prescription is that of equation (35) without the small-envelope correction of H02. This is because the corrected radius drops as the envelope mass becomes small, leading to a small wind-loss rate – the use of the fitted radius (which diverges as the envelope becomes small) ensures that the envelope is lost rapidly at the end of the AGB phase.

3.7 Hot-bottom burning

If the hydrogen envelope of an AGB star is sufficiently massive, the hydrogen-burning shell can extend into the convective region, a process known as hot-bottom burning (HBB). We deal with HBB by burning a fraction of the convective envelope f_HBB for a fraction of the time-step period f_burn at the temperature and density as fitted in Sections 3.3.4 and 3.3.5. The burned fraction is mixed with the rest of the convective envelope at the end of the time-step (usually coincident with the end of the interpulse period). The validity of this approach is discussed in Section 3.7.2. We burn only CNO, but HBB affects other elements via the NeNa and MgAl chains (see Karakas & Lattanzio 2003) as well as lithium via pp-burning. These are difficult problems which are being worked on but are outside the scope of this paper. HBB also depends on the mixing-length parameter because the base of the convective envelope is linked to superadiabatic layers near the surface by an adiabat; however, we ignore this and use a fixed value α= 1.75.

3.7.1 Clayton's CNO cycle

We burn the CNO elements according to Clayton's CNO bi-cycle (Clayton 1983). He claims this is accurate to 1 per cent for the temperature range we are considering (log₁₀T/K < 8), and this was confirmed by GdJ93 who used a similar approach. We calibrate out any errors (see Section 4). It is much faster than solving the differential equations of a complete nuclear reaction network.

The CNO cycle can be simplified from the full set of differential equations if ¹³N, ¹⁵N, ¹⁵O and ¹⁷F are in equilibrium. The cycle then splits into two, the CN cycle and the ON cycle, with branching ratios α_CN= 1 −γ and γ≃ 7 × 10⁻⁴ respectively (Angulo et al. 1999). The small value of γ reflects the fact that the time-scales in the ON cycle are many thousands of times those required to bring the CN cycle into equilibrium so we can treat the cycles separately.

The method of solution for equation (60) is found in Clayton (1983) (but beware the typographical error!). The time-step δt is substituted for t in equation (61) to calculate the abundances at the end of the current time-step.

The ON cycle equations are identical to the CN cycle, with ¹²C, ¹³C and ¹⁴N replaced by ¹⁴N, ¹⁶O and ¹⁷O and with appropriate τ_i (see Table 4 and Clayton 1983). The minor species ¹⁵N and ¹⁷O are assumed to be in equilibrium. Any ¹⁷F produced is assumed to decay immediately to ¹⁷O (τ_1/2≈ 65 s).

For short burning times (δt < τ₁₂) only the CN part of the cycle is necessary. For longer times, the CN cycle is burned to equilibrium before the ON cycle is activated. Even in the most massive AGB stars undergoing vigorous HBB, X_O16 does not change much so the ON cycle never approaches equilibrium. Nuclear reaction rates are taken from the formulae in the NACRE compilation (Angulo et al. 1999) except the beta decay constants which are taken from the compilation of Tuli (2000). The rates compare well to table 5.3 (p. 393) of Clayton (1983).

3.7.2 Thin shell burning versus whole envelope burning approximation

Between pulses the convective envelope of an AGB star turns over many thousands of times. It is impossible to model this using little CPU time because the burn → mix → burn → mix… process is computationally expensive, especially if the code is to be extended to more isotopes than just ¹H, ⁴He and CNO. Given the uncertainties involved in convective mixing and local mixing at the base of the envelope, it is simpler and preferable to approximate the burning (many times) of a thin HBB layer at the base of the convective envelope with a single burning of a larger portion of the envelope.

This can be justified by considering the size of the HBB region. For a 5-M_⊙ star d log₁₀(T/K)/dm at the base of the envelope is typically 3 × 10³M⁻¹_⊙. HBB ceases at log₁₀(T/K) ≈ 7.6 and the temperature at the base of the envelope is typically log₁₀(T/K) ≈ 8 for most of the TPAGB. So ΔM_HBB≈ 10⁻⁴M_⊙. This is much smaller than the size of the convective envelope (about 4 M_⊙ for a 5-M_⊙ star), so the HBB shell can be considered as thin.

When the thin HBB shell is burned and then mixed into the envelope the abundances in the envelope are essentially unchanged. Only once a significant number of mixings (of the order M_env/ΔM_HBB) have occurred will the envelope abundances change noticeably, so in our approximation we burn a fraction of the envelope, f_HBB, for a fraction of the interpulse period time f_burn and fit f_HBB and f_burn to our full stellar evolution models. In reality, some parts of the envelope burn more than once but this is absorbed into the calibration of f_HBB. The burned shell and the rest of the envelope are mixed at the end of the time-step.

To calibrate the model we allow mixing 10 times per interpulse period so the shape of the abundances versus time profiles between pulses can be observed. The code is designed so that the result is identical to that obtained if there is only one mixing per interpulse in the way that the code is expected to be used in population synthesis runs.

Note that this technique differs from that of GdJ93 where the average abundance rather than the final abundance between time 0 and f_burnτ_ip is mixed into the envelope.

4 HBB CALIBRATION AND COMPARISON OF SYNTHETIC MODELS WITH FULL EVOLUTIONARY MODELS

The free parameters,

• f_HBB – the fraction of the envelope of the star that is burned in the HBB shell,

• f_burn – the fraction of the interpulse period for which the HBB shell burns,

• f_DUP – the fraction of the dredged-up material that is hydrogen-burned before being mixed into the envelope to simulate dredge-up of the hydrogen shell in low-Z stars,

• f_DUP,burn – the fraction of the interpulse period for which the dredged-up material is burned, and

• N_rise – the factor used to define how quickly the HBB temperature reaches T_max,

are to be calibrated to our full stellar evolution models.

Previous authors (e.g. GdJ93) have used constant values for f_HBB, f_burn, f_DUP and f_DUP,burn, with a different prescription for the temperature (not requiring N_rise). Here we assume that these values differ for each star, so we attempt to parameterize them in terms of M_1TP and Z. Often we quote 10⁶f_burn instead of f_burn because f_burn≲ 10⁻⁶.

4.1 Calibration method

A Monte Carlo (MC) method is used to test the above free parameters with ranges 0.0 < f_HBB < 1.0, 0.0 < 10⁶f_burn < 10.0 and 0 < N_Trise < 20. f_DUP and f_DUP,burn are chosen to be zero and are only increased when necessary.

A weighted sum of squares is constructed from our full stellar evolution model nucleosynthesis data versus the corresponding synthetic model nucleosynthesis results to enable comparison between MC model runs. A is defined such that higher numbers mean a better fit where the weights are w_i= (w_C12, w_C13, w_N14, w_O16, w_C/O, w_C12/C13) = (1, 10, 1, 1, 5, 5) and s_i is the sum of squares difference between our full stellar evolution and synthetic models for the isotope (or ratio) i. The ratios X_C12/X_C13 and X_C/X_O are weighted preferentially because these are important observed nucleosynthetic constraints on AGB stars. ¹³C is also boosted because its abundance is small. 1D slices and 2D projections of the resulting parameter space are then examined and compared with the best fit obtained by this method. Human intervention comes last but proves invaluable when trying to fit the details. Appendix B contains the details of the calibration results.

4.2 Free parameter Heaven or Hell

The results of the MC runs for each star are shown in Tables 5, 6 and 7. Ranges are given when the MC technique cannot distinguish a unique solution. In such cases we choose a value that aids the fit of the free parameter to M_1TP and Z or such that 10⁶f_burn≈f_HBB. The chosen value is shown under the range. If the value is ‘–’ then there is no HBB so f_DUP=f_burn= 0.0. Where no value is given there is no full stellar evolution model with which to compare. It is not possible to use the best MC values for every star because there is too much non-systematic scatter.

For Z < 0.004 we include some immediate burning of dredge-up material.

4.3 Sensitivity to f_HBB, f_burn and N_rise

With the model described above and an appropriate choice of f_HBB, f_burn and N_rise, it is possible to match our full stellar evolution models to our synthetic models quite accurately. Problems occur with the fitted values of f_HBB, f_burn and N_rise because minor deviations in the fit of the free parameters produce very different output (the sensitivity does not help us to pin down unique values for f_burn and f_HBB because of their inherent degeneracy). For example, if N_rise is too small then ¹³C rises and falls too early in the M_i= 6 M_⊙ stars. If f_burn is even slightly too small then the ON cycle does not get switched on. If f_burn is slightly too large then more ¹⁴N is produced at the expense of ¹²C and ¹³C. A slight rise in f_HBB causes a large rise in HBB products, especially ¹⁴N. The sensitivity to f_HBB and f_burn is compounded when both are erroneous in the same direction.

A change to f_HBB affects the ¹²C surface abundance evolution for M_i= 5 M_⊙, Z= 0.02 (see Fig. 6). An increased f_HBB better fits the drop in ¹²C which occurs when HBB sets in, but by the end of the evolution the surface abundances are not very different. At most X_C12 is wrong by a factor of 1.4 at any point over the entire evolution, while overall it changes by a factor of 5. Final ¹³C and ¹⁴N have a similar scatter in log₁₀X of about 0.15. These effects are hardly visible in the case of the M_i= 6 M_⊙, Z= 0.02 star and, because no burning occurs at all at M_i/M_⊙= 4, it is only in the M_i/M_⊙≈ 4.5–5 transition zone that we have to worry about this.

Alteration of the burning time, f_burn, has essentially the same effect as a similar change in f_HBB with the exception of oxygen which is burned in the ON cycle when f_burn is long enough. The amount of ¹⁶O burned for M_i= 6 M_⊙, Z= 0.02 is very small in our full stellar evolution models (Δ log₁₀X_O16=−0.04). This is about twice the size of the spread with Δf_burn=±0.13 so we should not read too much into this. Significant oxygen burning occurs for M_i= 6 M_⊙, Z= 0.004 but the standard model deals quite well with this (see Fig. 7). The carbon and nitrogen abundances are weakly affected at M_i= 6 M_⊙ but at the transition mass (M_i= 5 M_⊙ for Z= 0.02, M_i= 4 M_⊙ for Z= 0.004) the surface abundance is sensitive to f_burn.

In summary, stars in the zone of transition between non-HBB and HBB (M_i= 4 M_⊙ for Z= 0.004, M_i= 5 M_⊙ for Z= 0.02) are the most troublesome when it comes to errors in the fit to f_HBB and f_burn. However, this transition is quite sharp, so not too many stars in a population would have the wrong surface abundances.

4.3.1 Temperature sensitivity

If the fit to the temperature of the HBB layer (Section 3.3.4) is allowed to vary even by a tiny amount, while the other free parameters are kept constant, CNO element production varies significantly. To show this, log₁₀T_max is allowed to vary from the fitted value by a factor of 0.98 < f_T < 1.02, just 2 per cent variations (5 per cent in T_max), and the abundance versus time profiles are compared for the case M_i= 6 M_⊙, Z= 0.02 which would ordinarily undergo large amounts of HBB (see Fig. 8). We do not always expect f_T= 1.00 to be the best fit because in reality the HBB layer has a temperature profile while in our synthetic model it does not. Note that the log₁₀T_max limit of 7.95 has been disabled for these tests, leading to numerical problems due to imaginary eigenvalues at f_T= 1.02.

• The final ¹²C is not greatly affected by temperature changes but f_T= 0.98 effectively switches off HBB. Paradoxically f_T= 1.02 burns less ¹²C than f_T= 1.01 during most of the evolution. This is because f_T= 1.02 puts the temperature above the log₁₀(T/K) = 7.95 limit of applicability of the burning code. The best fit is for f_T= 1.0.

• ¹³C is affected in a more subtle way. At low temperature more ¹³C is produced by incomplete CNO cycling. At the higher temperatures this ¹³C is converted to ¹⁴N. The final abundances for f_T > 1.0 are all similar because CN equilibrium is achieved, while for f_T < 1.00 there has not been enough conversion of ¹³C to ¹⁴N. Again f_T= 1.00 is the best fit.

• The log of the final surface abundance of nitrogen varies from −2.55 at f_T= 0.98 (the same as the abundance at the start of the TPAGB) to −1.85 at f_T= 1.02. The best fit is f_T= 1.01 although f_T= 1.00 is not too bad. For f_T≥ 1.01 nitrogen abundances are high because of excessive ON cycling.

• Out of all the CNO elements the surface oxygen abundance varies the most with temperature. For f_T≤ 1.00 there is little change in oxygen abundance, just as in our full stellar evolution models. An increase in the value of f_T to just 1.02 causes the surface oxygen abundance to drop by a factor of 10. This is not seen in our full stellar evolution models, so a value of f_T= 1.00 is certainly justified while f_T= 1.01 also gives too large a drop.

4.3.2 Dangerous interpolations and extrapolations

The TPAGB code is designed to be inserted into the rapid stellar evolution code, which deals with both single and binary stars for 0.1 ≤M_i/M_⊙≤ 100.0 and 0.0001 ≤Z≤ 0.03. This means that the expressions in the TPAGB code, developed here for 1 ≲M_i/M_⊙≲ 6.5 and 0.004 ≤Z≤ 0.02, or M_i≤ 2.25 M_⊙, Z= 0.0001, must be interpolated over a factor of 40 in Z, or extrapolated beyond Z= 0.004 for M≳ 2.5 M_⊙, into regimes where they have not been fitted to full stellar evolution models. Our expressions are designed to give sensible results when extrapolated or interpolated to low metallicity, but we have no way to tell if the results are correct.

5 DREDGE-UP CALIBRATION AND COMPARISON WITH OBSERVATIONS

5.1 Carbon stars

The carbon star luminosity function (CSLF) is defined as the number of carbon stars per unit bolometric magnitude for a particular population, i.e. it is a probability density function. We model a population of stars in the mass range 0.5 ≤M_i/M_⊙≤ 8.0 where the probability for each star is taken from the initial mass function of Kroupa, Tout & Gilmore (1993) (KTG93, see also Appendix A7) and we assume a constant star formation rate. Results are compared with the Large Magellanic Cloud (LMC) (Z= 0.008) and Small Magellanic Cloud (SMC) (Z= 0.004) data taken from Groenewegen (2002) (2002; see also Groenewegen 1999). The theoretical distributions are binned identically to the observed data in 0.25-mag bins. All distributions are normalized such that the integrated probability is 1.0 (so are independent of the star formation rate if we assume it is constant). Because the bin widths are fixed, the probability density is directly proportional to the probability per bin (i.e. the number of stars per bin) so is directly compared to the (suitably normalized) observations.

5.2 Modification of dredge-up parameters

Our best-fitting single-star models have ΔM^min_c=−0.07 and λ_min= 0.5 for the LMC and ΔM^min_c=−0.07 and λ_min= 0.65 for the SMC with the VW93 wind with a superwind at a Mira period of 500 d. These values are similar to those of Marigo (2001), noting that ΔM^min_c=−0.07 gives M^min_c≈ 0.59 M_⊙ at M_i≈ 1.5 M_⊙ (a typical mass for C-stars: Wallerstein & Knapp 1998).

Comparison with Fig. 5 shows that a value of ΔM^min_c nearer to −0.09 M_⊙ is necessary to match Marigo's prescription although the functional form of the prescription is otherwise similar.

5.3 Number ratios

We briefly consider the C/late-M number counts using the observations compiled by Groenewegen (2002) and the spectral type table of Jaschek & Jaschek (1995). We assume the CSLF dredge-up calibration above as our standard model and a constant star formation rate. It turns out to be very difficult to match the number ratios with our models unless we again change our free parameters. A good LMC fit (ΔM^min_c=−0.09, λ_min= 0.7, see Table 8) is obtained by increasing the amount of dredge-up, but then the CSLF peak is too dim. The SMC behaves in the opposite direction, with a decrease in the amount of dredge-up (ΔM^min_c=−0.06, λ_min= 0.4, see Table 9) which leads to a good fit but then the CSLF peak is too bright. We also show the effect of choosing a different wind prescription, with an η= 3 Reimers rate (Kudritzki & Reimers 1978) giving an even poorer fit than the VW93 wind (and too many bright stars in the CSLF). There are, however, many caveats to this simple approach. First, we have neglected the effect of a varying star formation rate which can be important when comparing number ratios (see Mouhcine & Lançon 2003). Secondly, the ratios are highly dependent on the spectral type which in turn depends on the stellar radius. It proved to be very difficult to fit the radius so any slight error on the fit leads to a resulting error in C/M. Thirdly, we neglect binary stars which produce giant branch and EAGB carbon stars (Izzard & Tout 2004). Fourthly, we do not take into account any observational selection effect, which is effectively the same as the assumption that the M-star surveys of the LMC and SMC are complete. This is unlikely to be true, while the C-star surveys (and resulting luminosity functions) are thought to be complete, although this helps us with the LMC and not the SMC. To simulate a magnitude-limited survey is beyond the scope of this paper. Since there is no systematic effect we will continue, in our ignorance, to use the dredge-up prescription calibrated by the CSLF.

5.4 Initial–final mass relation

The initial–final mass (M_i–M_f) relation is another check on the consistency of our models. Our synthetic results, including the dredge-up calibration, are plotted against the relation derived by Weidemann (2000) in Fig. 11. Weidemann 2000's results are partly based on evolution models and partly on observation, so it is difficult to draw any firm constraints from this comparison. The agreement is excellent for our Z= 0.02 models while the Z= 0.004 models are systematically high but always within 0.1 M_⊙. The effect of the choice of mass-loss prescription at intermediate metallicity (Z= 0.008) is shown in Fig. 12, where we use the rates of Kudritzki & Reimers (1978) and Blöcker (1995). No single mass-loss relation gives the best agreement, but a Reimers η= 1 rate is discrepant by more than 0.1 M_⊙ for M_i≳ 3 M_⊙, perhaps ruling this out.

5.5 White dwarf mass distribution

The observed white dwarf mass distribution provides an additional constraint on our models. The inherent uncertainty in the initial–final mass relation, which is due to the use of evolutionary models to calculate an initial mass, is removed. The comparison of our synthetic white dwarf mass distribution with observations is shown in Fig. 13. The observations are compiled from Bergeron, Saffer & Liebert (1992), Bergeron, Liebert & Fulbright (1995), Bragaglia, Renzini & Bergeron (1995), Dreizler & Werner (1996), Bergeron, Ruiz & Leggett (1997), Finley, Koester & Basri (1997), Marsh et al. (1997), Vennes et al. (1997), Napiwotzki, Green & Saffer (1999), Vennes (1999), Bergeron, Leggett & Ruiz (2001), Claver et al. (2001) and Silvestri et al. (2001), with no attempt to take into account any observational selection effects such as dimming of old white dwarfs or systematic biases such as the effect of a binary companion, metallicity or errors inherent in different white dwarf mass determination techniques (although for duplicate stars the newer data are believed over the old). The data are plotted as a histogram in 10 bins, with an equal number of stars in each bin. Our synthetic model curves for Z= 0.02, 0.008 and 0.004 are calculated with 10⁵stars for each metallicity between 0.1 and 8.0 M_⊙ and are binned in an identical way. Our dredge-up calibration of Section 5.1 is used. The initial mass function of KTG93 is used to calculate the contribution of each star to the histogram, and no attempt is made to mimic selection effects. The peak of the observations and theoretical curves are normalized to 1.0.

The peak position of the mass distribution is matched well by our synthetic models. We have a dearth of white dwarfs above 0.7 M_⊙ and below about 0.55 M_⊙. These are probably due to the omission of binary stars from our analysis. Binary-induced mass loss leads to low-mass helium white dwarfs and mergers to high mass (possibly oxygen–neon) white dwarfs.

6 SINGLE STAR YIELDS

In this section we calculate stellar yields appropriate for use in galactic chemical evolution models. We compare these yields with our full stellar evolution model yields and previously published yields of van den Hoek & Groenewegen (1997) and Marigo (2001).

6.1 Calculation of yields

The yield in equations (70) and (71) can be negative if an isotope is consumed (e.g. the hydrogen yield is always negative).

The contributions to the yields are calculated from equation (71) at each time-step. The initial mass is in the range 0.1 ≤M_i/M_⊙≤ 8.0. The stars are evolved for 16 Gyr and the yields (without supernova yields) shown are for the entire stellar lifetime. Full tables of yields are in Appendix D. We show the yields from our full stellar evolution models without attempting to compensate for mass loss that would occur after model breakdown (see Karakas & Lattanzio 2003).

6.2 Results

Most of the mass from each star that is expelled in the stellar wind is hydrogen or helium, and most of this is expelled in the TPAGB stage of the evolution of the star. Stars with M_i≳ 7–8 M_⊙ do not have a TPAGB stage, but explode as supernovae first, so the stellar wind yield from these stars is negligible. Stars with M_i≲ 0.8 M_⊙ do not evolve to the TPAGB in 16 Gyr so also have negligible yield. Comparison with the yields from (van den Hoek & Groenewegen 1997, (HG97), Marigo (2001) M01) and our full stellar evolution models is made. The yields of HG97 were calculated using a Kudritzki & Reimers (1978) mass-loss relation, with η= 4 (and η= 2 for Z= 0.004). The yields from M01 were calculated using the wind of VW93 (presumably with the M_i > 2.5 M_⊙ correction) for Z= 0.019 (but we compare them with our Z= 0.02 models) and mixing length parameter α_MLT= 1.68 rather than 1.75.

6.2.1 Comparison with our full stellar evolution models

Figs 14 and 15 compare our synthetic model yields (equation 71) with our full stellar evolution model yields using the HBB calibration of Section 4, without the dredge-up calibration of Section 5.1 and with the KTG93 initial mass function weighting. Our synthetic model does an excellent job of reproducing our full stellar evolution models for all isotopes considered except ²²Ne. Slight overproduction of ¹²C, ¹³C and ¹⁴N for M_i≳ 5 M_⊙ is due to small differences between our synthetic and full evolution models. The spike in production of ¹³C at around M_i≈ 4–5 M_⊙ is not an anomaly – rather it is not resolved on our full stellar evolution model mass grid unless it happens to lie on an integer multiple of M_⊙ (as is the case for Z= 0.02). We do not include ²²Ne destruction reactions so our synthetic yields are overestimates for M_i≳ 4 M_⊙.

6.2.2 Comparison with other models

We use the HBB and dredge-up calibrations made above (Sections 4 and 5.1) to compare our synthetic model yields with those of M01 and HG97. The yields are weighted with the KTG93 initial mass function – see Figs 16 and 17 and Appendix A7.

For most isotopes our yields lie between the values given by M01 and HG97. Our synthetic models experience less dredge-up than M01's in the range 1 ≲M_i/M_⊙≲ 3, leading to a relative underproduction of carbon and oxygen, though slightly more than HG97. The ¹³C spike position differs between our synthetic models and M01's by 0.5 M_⊙ at Z= 0.004, suggesting that we invoke HBB at a slightly higher mass, and also by a factor of about 2.5 (perhaps due to excess ¹²C seed). At lower metallicities (Z≤ 0.008) we overproduce nitrogen compared with the other yield sets, but M01 does not have AGB stars above 5 M_⊙ and our results are not dissimilar to HG97's Z= 0.004, η= 2 case (they publish no η= 2 yields for Z= 0.008). Our ¹⁵N agrees well with M01. The huge ¹⁶O production at low mass in HG97's data set for Z= 0.02 cannot be reproduced, although otherwise our yields agree well except for the lack of oxygen from dredge-up compared with M01 (also our intershell abundance is on average about 1 per cent compared with Marigo's 2 per cent). We produce negligible ¹⁷O compared with M01 although the amount is still rather small. The initial mass function-weighted ²²Ne yields show that the region of HBB (M_i≳ 4–5 M_⊙) is a relatively unimportant contributor to the total yield. It is also possible that the lower initial core masses from the Padova models have an effect on the evolution of M01's models, with lower core mass leading to lower luminosity and so a longer TPAGB phase with more dredge-up and enhanced CO yields.

7 CONCLUSIONS

We have presented a fast yet accurate synthetic model for TPAGB evolution based on state-of-the-art full evolution models. Yields calculated by the synthetic model and the full evolution model agree closely, with continuation of the synthetic model beyond the breakdown point of the detailed model for high-mass AGB stars. Calibration of third dredge-up to the Magellanic Cloud carbon star luminosity functions forces us to reduce our theoretically derived minimum mass for core mass by −0.07 M_⊙ and enforce a value of at least 0.8 − 37.5Z on the maximum dredge-up efficiency λ_max. The calibrated yields in general lie between the previously published yields of van den Hoek & Groenewegen (1997) and Marigo (2001), although our models produce more nitrogen and less carbon at low metallicity. Our synthetic model fits the initial–final mass relation to within 0.1 M_⊙ and reasonably fits the observed white dwarf mass distribution (we require binary stars to fit it properly). The next step is to include this synthetic model in the Hurley et al. (2002) binary code, and much progress has already been made in this direction (Izzard & Tout 2004,2003).

Acknowledgments

We thank John Lattanzio, Maria Lugaro, Richard Stancliffe, Carolina Ödman and especially the anonymous referee for useful discussion, encouragement and suggestions. RGI thanks PPARC for a scholarship, Monash Mathematics Department for putting up with him, and his own pocket for decent computing facilities. CAT thanks Churchill College for a fellowship.

REFERENCES

Footnotes