A&A 460, 565-572 (2006)
DOI: 10.1051/0004-6361:20066129

Population nucleosynthesis in single and binary stars

I. Model[*]

R. G. Izzard1,2,[*] - L. M. Dray2,3 - A. I. Karakas4 - M. Lugaro1,2 - C. A. Tout2

1 - Sterrekundig Instituut Utrecht, Postbus 80000, 3508 TA Utrecht, The Netherlands
2 - Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK
3 - Department of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK
4 - Origins Institute and Department of Physics & Astronomy, McMaster University, 1280 Main St. W., Hamilton, ON L8S 4M1, Canada

Received 28 July 2006 / Accepted 8 September 2006

We present a synthetic algorithm to rapidly calculate nucleosynthetic yields from populations of single and binary stars for use in population synthesis, globular cluster and Galactic chemical evolution simulations. Single star nucleosynthesis is fitted directly to full evolution models and our model includes first, second and third dredge-ups with s-process enhancements, an analytic calculation for hot-bottom burning of CNO, NeNa and MgAl isotopes, surface enhancements due to wind loss in massive stars and core-collapse supernova yields. Even though this algorithm operates about 107 times faster than full evolution and nucleosynthesis calculations, agreement with such models is good. We extend the single star model to include prescriptions of binary star interactions, notably mass loss and gain by stellar winds and Roche-lobe overflow, novae and type Ia supernovae. As examples of the application of our algorithm we present models of some interesting systems containing chemically peculiar stars that may occur in binaries.

Key words: stars: abundances - stars: AGB and post-AGB - stars: binaries: general - stars: chemically peculiar - stars: carbon - stars: Wolf-Rayet

1 Introduction

Most stars spend their lives in some kind of multiple system, and, indeed, perhaps up to half of all stars are in binary systems. Progress in understanding the evolution and nucleosynthesis of binary stars has been hampered by the difficulties posed by the possible interaction of the two components. This depends in turn on the properties of the stars and the orbital parameters of the system. Single stars on the other hand are relatively simple to construct and visualize as one-dimensional balls of gas, which lose matter to the interstellar medium through a combination of winds and explosions. The fate of a star and its contribution to the chemical enrichment of the Galaxy depends mainly on its initial mass and metallicity. Calculation of a time series of structure models of a single star involves following the key elements involved in the main nuclear-energy generation reactions. A complete evolutionary sequence of a star from the zero-age main sequence (ZAMS) to near its death can be computed in less than an hour on a modern desktop personal computer. These computational time estimates can extend to weeks or months when we consider many minor isotopes that are important for nucleosynthesis studies and the structure and nucleosynthesis of the thermally-pulsing asymptotic giant branch (TP-AGB) phase.

The computationally demanding nature of single stellar structure and nucleosynthesis models pose particular problems for population synthesis and chemical evolution studies, where the nucleosynthesis contribution and fate of thousands to millions of stars is required. This problem may be solved by using synthetic stellar evolution algorithms, which are based on analytic fits to detailed single-star models - such as SSE (Hurley et al. 2000) - to rapidly follow the total mass, core mass, luminosity and radius of a star as a function of time. Populations of binary stars pose even more of a problem. Even under the assumption that a binary star code is as fast as a typical single star code, with a mean runtime $\delta t\approx1~{\rm h}$, the size of the parameter space is overwhelming. Instead of variations of  $M_{{\rm ZAMS}}$ and metallicity Z we have two masses, M1and M2, initial separation a and perhaps the orbital eccentricity e (although Hurley et al. 2002 showed that consideration of eis unnecessary). If Z is fixed and all binaries are assumed to be in circular orbits then this leaves M1, M2 and a. For 100 points in each direction this sums to 106 stars, or a minimum of $10^{6}~{\rm h}$ of computing time. Added to this are the many uncertainties in binary star evolutionary theory such as mass transfer by winds, Roche Lobe overflow (RLOF), common-envelope evolution, binary-enhanced mass loss (e.g. the companion reinforced attrition process, CRAP, of Tout & Eggleton 1988), tidal interactions, stellar-wind collisions, stellar mergers, the fate of the accreted matter and thermonuclear explosions. The free parameters and uncertainties from single star evolution (e.g. convection, rotation and magnetic fields) also apply to binaries. One hundred years of computing time is a very conservative estimate to explore all of these uncertainties and free parameters with a detailed binary-star code (such as TWIN of Eggleton & Kiseleva-Eggleton 2002).

In this paper we develop a unique method which follows on from the work of Hurley et al. (2002,2000) and Izzard et al. (2004). The binary star code BSE (Hurley et al. 2002) is used alongside the single-star algorithm SSE (Hurley et al. 2000) to rapidly follow the evolution of single and binary stars and to include mass transfer and orbital dynamics. Nucleosynthesis and details of the TP-AGB phase are modelled by fitting to detailed evolutionary model results and with a synthetic nuclear-burning algorithm (Izzard et al. 2004). The runtime of the synthetic binary code with nucleosynthesis is typically $0.05~{\rm s}$ per system when following 126 isotopes which means we can evolve $10^{6}~{\rm stars}$ in less than 14 h. This is about seven orders of magnitude faster than the computational time required by the Monash stellar structure code plus post-processing nucleosynthesis (e.g. Karakas et al. 2002,2006) whose results for the light-element abundances have been used to prepare the fits for low and intermediate-mass stars. The computational speed-up provided by the synthetic algorithms means that the binary parameter space is no longer inexplorable.

We divide single stars into two groups. Low- and intermediate-mass stars with initial mass $M_{{\rm ZAMS}}$ less than about $8~{M_{\odot}}$, and massive stars with initial masses above about $8~{M_{\odot}}$. Low- and intermediate-mass single stars end their lives as giant branch or TP-AGB stars. Massive single stars are defined as those that ignite carbon in their cores and progress to a core-collapse supernovae (SNe). At solar metallicity the minimum limit for core collapse SNe is $8.25~{M_{\odot}}$, although this limit in detailed models is dependent on many factors including the mass-loss prescription and the treatment of convection, especially during the core H- and He-burning phases. We do not consider the contribution from super-AGB stars (e.g. Iben et al. 1997; Eldridge & Tout 2004; Siess 2006) because of the lack of detailed models available at the present time but work is in progress to resolve this point (Izzard & Poelarends 2006).

2 Low- and intermediate-mass stars

The basis for our synthetic nucleosynthesis model of low- and intermediate-mass stars is detailed in Izzard et al. (2004). Stellar evolution is modelled with the SSE prescription of Hurley et al. (2000). First and second dredge-up are treated as instantaneous merger in surface chemical abundances which are interpolated from the detailed stellar evolution models of Karakas et al. (2002) (hereinafter, the Monash models). TP-AGB evolution and nucleosynthesis are fitted to the Monash models. Third dredge-up is parameterised by the minimum mass for dredge-up  $M_{{\rm c,min}}$ and its efficiency $\lambda$. These are fitted to the detailed model results and are then modified to match observed Magellanic-cloud carbon-star luminosity functions (Izzard et al. 2004). We include hot-bottom burning (HBB) by fitting the temperature  $T_{{\rm bce}}$ and density  $\rho_{{\rm bce}}$ at the base of the convective envelopes and applying an analytic nuclear burning algorithm to obtain isotopic abundances.

Some changes and updates have been made since Izzard et al. (2004). The nucleosynthesis model now includes elements heavier than iron using new fits to theoretical slow-neutron capture intershell abundances from Lugaro et al. (2006), see Sect. 2.1. Better fits have been made to the Z=10-4 Monash models and the prescription for the third dredge-up has been improved: intershell abundances as a function of  $M_{{\rm ZAMS}}$ and pulse number are taken directly from tables of Monash model data (see Appendix A). Our HBB algorithm has been extended to include burning by the NeNa cycle and MgAl chain and it has been re-calibrated with these isotopes (see Sect. A.8).

2.1 The s-process

The production of elements beyond iron is possible in TP-AGB stars by slow neutron capture, known as the s-process. In the He intershell (i.e. the region between the H-burning shell and the He-burning shell) the reactions $^{13}{\rm C}(\alpha,n)^{16}{\rm O}$ and $^{22}{\rm Ne}(\alpha,n)^{25}{\rm Mg}$provide the neutrons. The $^{13}{\rm C}$ reaction activates at temperatures of about 9 $\times$ $10^{7}~{\rm K}$ but the $^{22}{\rm Ne}$ reaction requires T>3 $\times$ $10^{8}~{\rm K}$ and so occurs during thermal pulses only in more massive AGB stars. To create sufficient  $^{13}{\rm C}$, protons must be mixed from the convective envelope into the hydrogen-free intershell region. The physics of this process is still highly uncertain so in our model we assume a  $^{13}{\rm C}$ pocket exists and the abundance of  $^{13}{\rm C}$ (and hence the s-process neutron exposure) is a free parameter. We fit 48 elemental abundances from Ga to Bi to the models of Lugaro et al. (2006), calculated on the basis of stellar structure models produced with the FRANEC stellar evolution code (Straniero et al. 1997; Gallino et al. 1998). The models span the initial mass range $1.5\leq M_{{\rm ZAMS}}/{M_{\odot}}\leq5$and metallicity range $10^{-4}\leq Z\leq0.02$. Details of the fits can be found in Appendix A.7.

2.2 Rapid NeNa and 19F burning

HBB occurs in AGB stars more massive than about $4~M_{\odot}$when the base of the convective envelope dips into the top of the H-burning shell. It is a rich source of nucleosynthesis and is a particularly important source of nitrogen produced by the CNO cycle. The NeNa cycle is similar in many regards to the CNO cycle and also operates during HBB. The cycle is a closed loop (if we neglect proton capture on $^{24}{\rm Mg}$as a reasonable approximation) and the isotopic abundances can be solved for by an eigenvalue method (see e.g. Clayton 1983). The net result of NeNa cycling is to destroy $^{21}{\rm Ne}$ and  $^{22}{\rm Ne}$ and create $^{20}{\rm Ne}$ and  $^{23}{\rm Na}$.

Entry via $^{19}{\rm F}$ and exit via $^{23}{\rm Na}(p,\gamma){}^{24}{\rm Mg}$is dealt with by explicitly reducing the abundances on the appropriate timescale. The beta decays have short lifetimes (Tuli 2000) so are assumed to be in equilibrium. The differential equations become

              $\displaystyle %
\frac{{\rm d}}{{\rm d}t}\left[\begin{array}{c}
^{20}{\rm Ne}\\ ...
...rm 21}}{\rm Ne}\\
^{{\rm 22}}{\rm Ne}\\
^{{\rm 23}}{\rm Na}\end{array}\right]$ = $\displaystyle \left[\begin{array}{cccc}
-\frac{1}{\tau_{20}} & 0 & 0 & \frac{1}...
...rm 21}}{\rm Ne}\\
^{{\rm 22}}{\rm Ne}\\
^{{\rm 23}}{\rm Na}\end{array}\right]$  
  = $\displaystyle \frac{{\rm d}}{{\rm d}t}{\bf U} = \left[\Lambda\right]{\bf U},$ (1)

where the isotope symbols represent number densities and the reaction timescales $\tau_{i}$ are calculated from Angulo et al. (1999) with the temperature  $T_{{\rm bce}}$ and density  $\rho_{{\rm bce}}$taken from the fit by Izzard et al. (2004) (with minor improvements, see A.8). This is an eigenvalue problem with solutions of the form

{\bf U}(t)=\sum_{i=0,1,2,3}A_{i}{\rm e}^{\lambda_{i}t}{\bf U}_{i}.
\end{displaymath} (2)

The details of the solution for ${\bf U}$ are given in Appendix B. A detailed comparison of the variation of synthetic and detailed model abundances with time is in Appendix D.

2.3 Rapid MgAl burning

Also active during HBB is the MgAl chain. This chain is similar to the NeNa cycle with neon and sodium replaced by magnesium and aluminium respectively. There are some differences such as the competition between $\beta$-decay and proton capture as the dominant $^{26}{\rm Al}$ destruction method. The $^{27}{\rm Al}(p,\alpha)^{24}{\rm Mg}$ reaction is prohibitively slow at HBB temperatures so the cycle really is a chain, terminating at  $^{27}{\rm Al}$.

We make a number of assumptions in order to solve for the abundances. As mentioned, the chain is terminated at $^{27}{\rm Al}$ because the timescale for $^{27}{\rm Al}(p,\alpha)^{24}{\rm Mg}$ or $^{27}{\rm Al}(p,\gamma)^{28}{\rm Al}(\beta,\nu)^{28}{\rm Si}$ is long compared to the burning time. We do allow for the leak of  $^{27}{\rm Al}$ back to  $^{24}{\rm Mg}$ but this proves to be negligible. The beta-decays of  $^{25}{\rm Al}$ and  $^{27}{\rm Si}$are quick enough to be considered instantaneous. We consider both the ground and metastable states of  $^{26}{\rm Al}$[*] although not including the metastable state has very little influence on the surface abundances of  $^{26}{\rm Al}$.

The rate equations are solved analytically (see Appendix C). Detailed comparison of synthetic to detailed model surface abundances vs. time is given in Appendix D.

3 Massive stars

According to the SSE algorithm, stars in the (initial) mass range $8\la M_{{\rm ZAMS}}/{M_{\odot}}\la24$evolve on to the early-asymptotic giant branch (EAGB) but never make it to the TP-AGB because their cores grows large enough to ignite carbon - such stars explode as type-II supernovae. For masses greater than about $25~M_{\odot}$ (at solar metallicity) the star is luminous enough that line-driven wind mass loss removes the hydrogen envelope to expose the helium core. Such stars are known as Wolf-Rayet (WR) stars (Chiosi & Maeder 1986). As the hydrogen envelope is stripped, deeper and hotter layers of the star are exposed. These layers are rich in $^{14}{\rm N}$ owing to CNO cycling and are the WNL stars. We use the same definition of WR sub-types as in Dray et al. (2003) and Maeder & Meynet (1994), without the effective temperature condition which does not work too well with the SSE code. As the hydrogen itself is removed the star becomes a helium star, this is the WNE phase. The products of helium burning are then also exposed as mass loss continues. This leads to a WC then a WO phase. These stars explode as type Ib/c supernovae if the (degenerate) core mass in any post-helium main sequence stage exceeds $M_{{\rm Ch}}$. The mass-loss rate and hence WR evolution is strongly dependent on metallicity (Kudritzki et al. 1989; Eldridge & Vink 2006).

We fit to the Dray models (see Sects. 3.1 and 3.2). It has proved extremely difficult to implement a simple synthetic burning algorithm because the interplay between burning shells and convective regions in massive stars requires a detailed model. So we simply interpolate from tables to obtain surface abundances of $^{1}{\rm H}$, $^{4}{\rm He}$, $^{12}{\rm C}$, $^{14}{\rm N}$, $^{16}{\rm O}$ and $^{20}{\rm Ne}$. Section 3.3 introduces yields from core-collapse supernovae and Sect. 3.4 describes how we set the supernova remnant mass.

3.1 The Dray models

The Dray models of massive single stars are calculated with the Eggleton STARS stellar evolution code. Details can be found in Dray et al. (2003) and Dray & Tout (2003). The models cover the initial mass range[*] $10\leq M_{{\rm ZAMS}}/M_{\odot}\leq150$ and metallicity range $10^{-4}\leq Z\leq0.03$. The SSE model fits, which we use in our synthetic code to model stellar evolution but not nucleosynthesis, were also made with the Eggleton code so the stellar evolution should be consistent with the Dray models. Subsequent code improvements and the lack of mass-loss in the SSE models means this is not always the case. Nucleosynthesis follows ${\rm H}$, ${\rm He}$, ${\rm CNO}$, ${\rm Ne}$ and Mg. We use two mass-loss rates, referred to as MM (after Maeder & Meynet 1994) and NL (after Nugis & Lamers 2000) details of which can be found in Dray et al. (2003) and Appendix F.

3.2 Surface abundances

For each of the two mass-loss rates we construct a table of surface abundances ( $^{4}{\rm He}$, $^{12}{\rm C}$, $^{14}{\rm N}$, $^{16}{\rm O}$ and ${\rm Ne}$) as a function of metallicity Z, $M_{{\rm ZAMS}}$ and $M/M_{{\rm ZAMS}}$. Linear interpolation is then performed on this table to obtain the surface abundances at any stage of evolution. While the evolutionary timescales of the rapid model do not exactly agree with the Dray models, this technique has the virtue of being simple and gives the same stellar wind yields (assuming the same amount of mass is lost prior to supernova).

Their ${\rm Ne}$ contains both $^{20}{\rm Ne}$ and $^{22}{\rm Ne}$but the degeneracy can be broken by noting that  $^{22}{\rm Ne}$is made from double alpha-capture on $^{14}{\rm N}$, so we calculate $^{22}{\rm Ne}$ from the rate of destruction of  $^{14}{\rm N}$. Once $^{14}{\rm N}$ is exhausted in the WC or WO phase, any extra neon is $^{20}{\rm Ne}$. Surface $^{24}{\rm Mg}$ never changes appreciably in the Dray models. The mass fraction of hydrogen is calculated by subtracting the sum of all the other abundances from one. We approximate the surface  $^{13}{\rm C}$ by assuming it is in equilibrium with  $^{12}{\rm C}$, with an abundance  $^{12}{\rm C}/4$, when surface nitrogen is more than $95\%$ of the total CNO content, that is, when the CN cycle has come to equilibrium.

3.3 Core collapse supernovae

The usual fate of massive single stars is to explode as core-collapse SNe after hydrostatic burning. Neutrinos from electron capture reactions on protons are thought to yield the energy required, up to $10^{53}~{\rm erg}$, to power the ejecta. The remaining matter forms either a neutron star (NS) or black hole (BH). Our supernova yields are interpolated as a function of  $M_{{\rm CO}}$ and Z from either Model sets A - default, B and C Woosley & Weaver (1995), with the envelope removed according to the method of Portinari et al. (1998), or Chieffi & Limongi (2004) as a function of  $M_{{\rm CO}}$, Z and the mass cut $M_{{\rm NS/BH}}$ (see Sect. 3.4 below). To aid future investigations we have included the r-process by yielding a fixed mass per supernova according to the abundances derived in Arlandini et al. (1999) or Simmerer et al. (2004).

3.4 SN remnants and their masses

The SSE model defines the remnant neutron star or black hole baryonic mass (the mass cut) as a linear function of the core mass

M_{{\rm NS/ BH}}=1.17+0.09M_{{\rm CO}}
\end{displaymath} (3)

with BH formation for $M_{{\rm CO}}>7~{M_{\odot}}$. This gives a minimum BH baryonic mass, and corresponding maximum NS baryonic mass, of $1.8~{M_{\odot}}$ and a BH baryonic mass of $3~{M_{\odot}}$for a  $20~{M_{\odot}}$ carbon-oxygen (CO) core. Observations of BHs in X-ray binaries (Orosz 2003) suggest far higher masses, perhaps up to $18~{M_{\odot}}$ (Fryer & Kalogera 2001). We have added the BH mass prescription of Belczynski et al. (2002) to the SSE code (see Appendix E) to provide an option for high-mass single star BHs. Note that for $M_{{\rm CO}}\geq7~{M_{\odot}}$ according to this prescription the entire stellar mass is swallowed into a BH, corresponding to an upper initial mass for ejection of mass in a SN of $70~{M_{\odot}}$ for the SSE wind prescription (see Appendix F.1), $85~{M_{\odot}}$ for MM and $45~{M_{\odot}}$ for NL.

4 Binary stars

The BSE algorithm already deals with the dynamics of binary evolution, including orbital motion, tidal interaction, wind loss and accretion, common envelopes, mergers, novae and type Ia supernovae. It remains to link these processes with our synthetic nucleosynthesis model.

4.1 Wind accretion

As in Hurley et al. (2002) wind accretion is assumed to occur by the Bondi-Hoyle mechanism (Bondi & Hoyle 1944) averaged over a binary period. If both stars have stellar winds the BSE model ignores any interaction between them but, for nucleosynthetic purposes, we must follow the composition of the accreted matter. We model the collision of the two winds by balancing momentum fluxes in a similar way to Huang & Weigert (1982). The momentum flux at a distance r from each star is $J=\dot{M}v_{{\rm w}}/4\pi r^{2}$where  $v_{{\rm w}}$ is the velocity of the wind and $\dot{M}$the mass-loss rate. We assume the momentum flux of both stars is equal at the point where they interact and shock J1=J2 such that

\end{displaymath} (4)

where the $\dot{M}v$ terms are factored into fv. The stars are a distance a apart (the orbit is assumed circular) so r1+r2=awhich gives r1/a=fv/(1+fv) and r2/a=1/(1+fv). For each star ri is compared to the stellar radius Ri. If ri>Ri the shock is outside the star and none of the donor material reaches the surface. If $r_{i}\leq R_{i}$ the shock is inside, or rather on the surface of, the star and all the accreted material is from the donor. In reality some mixing will occur at the interface of the two winds (e.g. Ruffert 1994; Foglizzo & Ruffert 1999; Walder & Folini 2000b,a) but the complexities are too poorly understood to do more than the simple prescription described above.

4.2 Thermohaline mixing and the accretion process

Material which accretes onto a star either remains on the surface or sinks and mixes due to thermohaline instability (Kippenhahn et al. 1980; Chen & Han 2004; Proffitt 1989) or convective instability in the stellar envelope. We assume that accreted material forms a separate accretion layer on the surface of the star or, if one already exists, the accreted material is mixed into it. In this way the number of sets of abundances that we keep track of is only two. While this is required to keep the code from being too slow it is at the expense of the consideration of abundance gradients in the accreted material. Such consideration is rarely required. The accretion layer usually has a higher molecular weight than the envelope because it comes from a star with significant wind loss which occurs during the later stages of stellar evolution when heavy isotopes are mixed to the surface. In this case, or when the envelope is convective, we mix the accretion layer instantaneously with the stellar envelope to form a new, polluted envelope at each timestep.

4.3 RLOF and common envelopes

The prescription used in the BSE code for RLOF is described in detail in Hurley et al. (2002). There are several consequences of RLOF for nucleosynthesis. The first is the truncation of phases of stellar evolution when the radius is large as in giant branch (GB) and AGB stars. This means there are relatively fewer giants in binaries and the nucleosynthesis associated with dredge-up in giants does not take place, to a greater or lesser extent depending on the initial binary distributions. The other effect of RLOF is due to accretion, or the lack of it, onto the secondary star. RLOF need not be conservative so matter can be lost to the interstellar medium (ISM) directly. Rapid accretion can also lead to a common envelope (CE) forming around the stars (Paczynski 1976). While their cores spiral toward one another the CE may be expelled into the ISM. Alternatively, the cores of the stars may merge in the envelope to form a new star.

We assume that CE evolution occurs quickly relative to the nuclear burning timescales at the surface of the in-spiralling cores so there is no nuclear processing. This is partly for simplicity and partly because the CE process is already uncertain. Detailed and reliable hydrodynamical models of CEs are lacking. It is also assumed that no matter accretes from the common envelope on to either star (Hjellming & Taam 1991).

If there is no CE, the stars may completely merge (Sect. 4.4) or relatively slow accretion can lead to explosive events such as type Ia supernovae (Sect. 4.5), novae (Sect. 4.6) or the acquisition of a new stellar envelope and rejuvenation of the star.

A first dredge-up is forced for Hertzsprung gap (HG) and GB stars which undergo a common envelope phase. We expect that either dynamical mixing effects owing to the spiral-in process completely mix the envelope (in an analogous way to convection at the dredge-up) or the envelope is completely ejected. In the latter case we must also eject previously burnt material from inner layers of the envelope: so we force a dredge-up to mix burnt material into the envelope and then eject it.

4.4 Stellar mergers

The treatment of stellar mergers in Hurley et al. (2002) deals with each of the 15 $\times$ 15 possibilities involved, there being 15 stellar types (Hurley et al. 2002 or Izzard et al. 2004 for definitions). Symmetry reduces this by almost a factor of two but still there are more than 100 possibilities. Fortunately, for nucleosynthesis, all that is required is knowledge of where the material goes. We assume there is no additional nucleosynthesis during the merger process but burning rates could increase if protons are mixed into hot material (Ivanova et al. 2002). This is reasonable given the long nuclear timescales of stars compared to their dynamical timescales. The new stellar envelope is formed by one of the following routes.

Star 2 is compact compared to the envelope of star 1 so star 2 merges with the core of star 1 to form a new core. The new envelope composition is the same as that of star 1 prior to the merger. This happens when a helium star or white dwarf (WD) merges with a giant. It is also possible to form a new envelope from the whole of star 1 if star 1 is a main sequence (MS) star accreting on to e.g. a helium star.

As above with 1 and 2 interchanged.

Both cores are compact compared to their envelopes so the new envelope is a mixture of the existing envelopes and both cores sink.

A Thorne-Zytkow object is formed by the merger of a neutron star (NS) or black hole (BH) with a pre-compact object star. The envelope is removed instantaneously without any extra nucleosynthesis (Cannon 1993). WDs and other NS/BHs collide with a NS/BH to form a heavier NS/BH (BH if $M>1.8~{M_{\odot}}$) again without any nucleosynthesis.
The route taken by each type of merger is shown in Appendix G, where the stellar type numbers correspond to those of Hurley et al. (2002). Note that collision of two helium white dwarfs (HeWDs) may lead to an explosion (see below).

4.5 Type Ia supernovae and AICs

Three types of type Ia supernova (SNIa) are considered here. Edge-lit detonations are thought to occur when $0.15~{M_{\odot}}$of helium-rich matter is accreted on a dynamical timescale on to a sub- $M_{{\rm Ch}}$ carbon-oxygen white dwarf (COWD). This was modelled in two dimensions by Livne & Arnett (1995) who made eight models with CO core mass $0.55\leq M_{{\rm CO}}/{M_{\odot}}\leq0.9$and helium layer mass  $M_{{\rm He}}$ between 0.1 and $0.2~{M_{\odot}}$. Their yields are fitted to functions of $M_{{\rm CO}}$ and  $M_{{\rm He}}$(see Appendix H).

A SNIa caused by a COWD which reaches the Chandrasekhar mass by steady accretion of hydrogen or helium-rich material, accretion from another COWD or merger with another COWD, is modelled with the explosive yields of the Iwamoto et al. (1999) DD2 model. They claim this model best fits observed spectra and lightcurves. We do not include the disc wind of Hachisu et al. (1996) that enables large numbers of SNeIa to be formed by accretion from sub giants and early giants. We can, however, get the right rate.

Helium white dwarfs which accrete helium-rich matter until their total mass exceeds $0.7~{M_{\odot}}$ explode with the yields of Woosley et al. (1986), scaled to the ejected mass $\Delta M$by a factor $\Delta M/0.664~{M_{\odot}}$. Strictly these yields are applicable only for the accretion of helium on to the helium WD but are used in the absence of other models for the merger of two HeWDs.

Accretion-induced collapse (AIC) to a NS owing to accretion of material on to an oxygen-neon white dwarf (ONeWD), which is not really a SNIa and only occurs in binaries, produces zero yield according to Nomoto & Kondo (1991). This has recently been challenged by Qian & Wasserburg (2003) who speculate that there may be some r-processing in a wind leaving the nascent NS. The situation is unclear and in any case there are no published yields to include in our synthetic model so we assume there is zero yield from an AIC.

4.6 Novae

Accretion of hydrogen-rich material on to a WD at a rate $\dot{M}<1.03\times10^{-7}~{M_{\odot}}~{\rm yr}^{-1}$ (Warner 1995) leads to unstable nuclear burning in explosive novae. During the explosion hydrogen is converted to helium and the temperature is high enough to activate the CNO, NeNa and MgAl cycles. Novae are thought to contribute a significant fraction of the Galactic content of $^{13}{\rm C}$, $^{15}{\rm N}$ and  $^{17}{\rm O}$. The most complete set of yields is that of José & Hernanz (1998) who evolve 14 sequences spanning a CO/ONeWD mass range of $0.8{-}1.35~{M_{\odot}}$. Mixing of accreted material with the surface layer of the WD is essential to the explosion. While mixing fractions of $25{-}75\%$ are considered by José & Hernanz (1998), we adopt the $50\%$ mixed models (models CO 2, 3 and 5, ONe 1, 3, 5 and 6) and linearly interpolate a table of yields. The fraction of accreted matter retained after the explosion, $\epsilon_{{\rm nov}}$, is set to 10-3.

4.7 Other binary-specific nucleosynthesis

The impact of a supernova on the companion star is not included in BSE, other than its effect on the orbital dynamics. Matter from the SN explosion can either be accreted by a companion or can strip the companion of matter. The latter is more likely (Wheeler et al. 1975; Marietta et al. 2000) and the companion star probably survives the stripping (Taam & Fryxell 1984). Accretion may occur from a weak supernova such as an AIC and such a process may explain stars which are simultaneously r-process and s-process rich (Qian & Wasserburg 2003; Ivans et al. 2005). Given these uncertainties mass accretion and stripping from SNe are not currently included in our synthetic model.

Our massive-star nucleosynthesis model (Sect. 3) necessarily uses $M_{{\rm ZAMS}}$ as a parameter but this is a problem during binary evolution because a significant amount of mass can be accreted from a companion. The accreted material is dealt with in Sect. 4.2 above and it probably mixes by a thermohaline instability with the stellar envelope. When the star evolves further its surface abundances change according to a prescription based on $M_{{\rm ZAMS}}$ rather than the now larger M so $M_{{\rm ZAMS}}$ is set to M when matter is accreted. Surface CNO abundances are calculated by adding the change in abundance from the ZAMS value due to evolution and mass loss (which would also occur in an equivalent mass single star) to the change due to accretion. The helium abundance is taken to be the maximum of the envelope abundance after accretion and of the equivalent abundance of the Dray model with the same mass, fraction of mass stripped and metallicity.

There is also the problem of stellar phases not modelled by any of the full stellar evolution codes. A good example is a low-mass helium star formed by stripping of a red giant envelope before core helium burning has finished. We simply convert hydrogen to helium, and all CNO to  $^{14}{\rm N}$ in these stars, under the assumption that there is no change in the heavier isotopes. This is reasonable because they are cooler than their more massive equivalents, such as those in the Dray models, which also do not burn, for example, neon or magnesium.

5 Results: example systems

In this section we present four systems which demonstrate the unusual nucleosynthesis that can occur in binary systems and the ability of our synthetic model to reproduce these situations. Our synthetic code has been and is being applied to other problems in both single and binary stars, such as an investigation of dim carbon stars in the Magellanic clouds (Izzard & Tout 2004), calibration of the s-process efficiency in giants and post-AGB stars (Bonacic Marinovic et al. 2006; Bonacic Marinovic et al. in preparation), models of super-AGB stars (Izzard & Poelarends 2006), a study of the effect of nuclear reaction rate uncertainties in hot-bottom burning AGB stars (Izzard et al. 2006, in prep., and Izzard et al. 2006), an investigation into the origin of the R-stars (Izzard et al, in prep.). It has also been coupled with a galactic chemical evolution code to directly investigate the effect of binaries on chemical evolution. The real power of our synthetic model is not only displayed by the few examples shown here, but also in the population studies mentioned above.

\par\mbox{\includegraphics[width=8cm,clip]{aa6129f1a.eps}\includegraphics[width=8.35cm,clip]{aa6129f1b.eps} }
\end{figure*} Figure 1: Evolution in the surface abundance-mass plane for a binary with initial conditions $M_{1}=7~{M_{\odot}},\; M_{2}=5~{M_{\odot}},\; P=40~{\rm days},\; Z=0.02$, the primary is black, the secondary is grey. Z marks the zero-age main sequence, while the numbers show mass-transfer by RLOF at 1, 3 and 4 and common envelope evolution at 2 and 5. The arrows show the direction of mass-transfer. A detailed description of the evolution is given in the text.

\par\mbox{\includegraphics[width=7.7cm,clip]{aa6129f2a.eps}\includegraphics[width=8.2cm,clip]{aa6129f2b.eps} }
\end{figure*} Figure 2: Evolution in the surface abundance-mass plane for a binary with initial conditions $M_{1}=3~{M_{\odot}},\; M_{2}=1~{M_{\odot}},\; P=3700~{\rm days},\; Z=10^{-4}$. Symbols are as in Fig. 1, with RLOF at point 1. The secondary accretes enough $^{12}{\rm C}$-rich material to become an extrinsic carbon star, with a very high $^{12}{\rm C}/^{13}{\rm C}$ ratio.

\includegraphics[width=8.1cm,clip]{aa6129f3b.eps} }
\end{figure*} Figure 3: Evolution in the surface abundance-mass plane for a binary with initial conditions $M_{1}=6~{M_{\odot}},\; M_{2}=5~{M_{\odot}},\; P=5000~{\rm days},\; Z=10^{-4}$. Symbols are as in Fig. 1, with a common envelope phase at point 1. Due to the high mass of the primary, hydrogen-burning products (notably $^{13}{\rm C}$ and  $^{14}{\rm N}$) are transferred to the secondary.

\includegraphics[width=8.1cm,clip]{aa6129f4b.eps} }
\end{figure*} Figure 4: Evolution in the abundance-mass plane for a binary with initial conditions $M_{1}=40~{M_{\odot}},\; M_{2}=20~{M_{\odot}},\; P=2000~{\rm days},\; Z=0.03$. Symbols are as in Fig. 1, with RLOF at point 1.

Here we present four example systems which were calculated with a default set of parameters and the initial masses, period and metallicity are varied. Notable parameters are eccentricity e=0, the mass-loss prescription of Hurley et al. (2002), common envelope parameter $\alpha=3.0$, a supernova kick velocity dispersion of $190~{\rm km~ s^{-1}}$and no extra CRAP.

6 Conclusions

Our synthetic nucleosynthesis model for low- and intermediate-mass single stars has been extended from that of Izzard et al. (2004). The fits for abundance changes at first and second dredge-up have been improved. Isotopes up to iron and the s-process elements are now included in the intershell material dredged up on the TP-AGB. Our HBB algorithm now burns the CNO and NeNa cycles and the MgAl chain and can reasonably well fit the Monash detailed model yields even though it runs more than seven orders of magnitude faster.

We have constructed a synthetic nucleosynthesis model of stars in the mass range $8\leq M/{M_{\odot}}\leq100$ and metallicity range $10^{-4}\leq Z\leq0.03$ with two mass-loss prescriptions by fitting surface abundances to detailed stellar evolution models. Hydrogen, helium, the CNO isotopes, neon and magnesium are fitted for the pre-SN evolutionary phases. The supernova yields of Woosley & Weaver (1995) are fitted for all available isotopes with pre-SN envelope removal according to the prescription of Portinari et al. (1998). The synthetic model yields reasonably well approximate the detailed model yields.

We have presented a model for synthetic nucleosynthesis in binary stars. The details of mass transfer, accretion, mixing and possible explosions are considered. Some example systems are presented to demonstrate the flexibility of our model and its ability to produce exotic stars.

We are now in a position to use our model: in Paper II we shall give results from a parameter space study and Paper III will compare yields from single and binary stars.

We thank Roberto Gallino for sharing unpublished results. RGI wishes to thank PPARC for a studentship while at the IoA, CIQuA ( for computing and monetary support, the NWO for his current fellowship in Utrecht and Axel Bonacic, Ross Church, John Eldridge, Evert Glebbeek, Simon Jeffery, John Lattanzio, Carolina Ödman, Onno Pols, Sean Ryan and Richard Stancliffe for useful discussion and advice. LMD gratefully acknowledges funding from the Leicester PPARC rolling grant for theoretical astrophysics and PPARC for a studentship while at the IoA. CAT thanks Churchill college for a fellowship. ML gratefully acknowledges the support of NWO through the VENI grant.



7 Online Material

Appendix A: New fits for low and intermediate mass stars

Please note that in this section of the Appendix M refers to the initial stellar mass (otherwise known as $M_{{\rm ZAMS}}$) unless stated otherwise.

A.1 TP-AGB luminosity

We have updated the luminosity of our TP-AGB stars function to better fit the Monash models, especially at low metallicity. We treat the total luminosity L as a combination of that from the radiative region of hydrogen burning, $L_{{\rm c}}$, and from hot-bottom burning  $L_{{\rm env}}$.

$\displaystyle %
L_{{\rm c}} = -39557+62747M_{{\rm c,nodup}}^{0.59552},$     (A.1)

where $M_{{\rm c,nodup}}$ is the core mass in the absence of third dredge up (see Izzard et al. 2004), and if $M_{{\rm c,nodup}}>0.6$we apply
$\displaystyle %
L_{{\rm c}} = \max\Big(L_{{\rm c}},\left[16853+4959M_{{\rm c,nodup}}\right]
\left[M_{{\rm c,nodup}}-0.5576\right]^{0.32988}\Big).$     (A.2)

We calculate the envelope luminosity from
$\displaystyle %
L_{{\rm env}} = L_{{\rm env,peak}}\left(\frac{M_{{\rm env}}/{M_{\odot}}-1}{M_{{\rm env,1TP}}/{M_{\odot}}-1}\right),$     (A.3)

$\displaystyle %
L_{{\rm env,peak}}/{L_{\odot}} = \max\left(0,~(6113.6-2071.9M_{{\rm 1TP}}\right)
\left(1+1.3864\ln Z\right))$     (A.4)

and set $L_{{\rm env}}$ to zero if $M_{{\rm env}}<1~{M_{\odot}}$. The final luminosity L is calculated by adding $L_{{\rm c}}$and  $L_{{\rm env}}$ with a turn-on modulation factor which is fitted for each Monash model and interpolated from Table A.1,
$\displaystyle %
L/{L_{\odot}} = \max\left(1,f_{{\rm turnon}}\left[L_{{\rm c}}+L_{{\rm env}}\right]\right),$     (A.5)

$\displaystyle %
f_{{\rm turnon}} = \max(f_{{\rm turnon,min}},1-(1-f_{{\rm turnon,min}})
\exp\left(-\max\left[1,N_{{\rm TP}}/N_{{\rm TO}}\right]\right),$     (A.6)

and $\alpha$ is adjusted by fitting the minimum luminosity at the first thermal pulse and  $N_{{\rm TO}}$ is adjusts the rise rate during the first few pulses (see Table A.1).

Table A.1: Luminosity parameters.
$M/{M}_\odot$ Z $N_{{\rm TO}}$ $f_{\rm turnon,min}$
1 0.0001 4 0.8
1 0.004 3.5 0.25
1 0.008 3.5 0.25
1 0.02 11 0.4
2 0.0001 9 0.3
2 0.004 4.0 0.25
2 0.008 2.0 0.25
2 0.02 12 0.35
3 0.0001 5 0.5
3 0.004 30 0.55
3 0.008 29 0.55
3 0.02 15 0.35
4 0.0001 18 0.25
4 0.004 28 0.35
4 0.008 22.5 0.4
4 0.02 15 0.4
5 0.0001 45 0.25
5 0.004 45 0.35
5 0.008 52 0.4
5 0.02 24 0.4
6 0.0001 40 0.25
6 0.004 35 0.35
6 0.008 52 0.4
6 0.02 17 0.4

A.2 TP-AGB radius

The TP-AGB radius follows the luminosity closely, so we fit

$\displaystyle %
\log_{10}R = \alpha+\beta\log_{10}L-\frac{1}{3}\log_{10}\left(\frac{M}{M_{\rm 1TP}}\right),$     (A.7)

where $\alpha$ and $\beta$ are fitted to the Monash models and then linearly interpolated from a table as a function of mass and metallicity (see Table A.2).

Table A.2: Radius fitting parameters.
Z $M/{M}_\odot$ $\alpha$ $\beta$
0.02 1 -0.5047 0.8163
0.02 2 -0.44956 0.75637
0.02 3 -0.85366 0.84665
0.02 4 -0.97215 0.85531
0.02 5 -0.80457 0.80302
0.02 6 -0.083275 0.62613
0.008 1 -0.43725 0.76396
0.008 2 -0.558475 0.766392
0.008 3 -1.21167 0.919587
0.008 4 -1.345 0.92926
0.008 5 -0.47064 0.70514
0.008 6 -0.62494 0.72929
0.004 1 -0.28398 0.69964
0.004 2 -1.0385 0.89392
0.004 3 -1.3532 0.94287
0.004 4 -1.0266 0.8409
0.004 5 -0.35658 0.66565
0.004 6 -0.39832 0.66837
0.0001 1 -0.45898 0.69964
0.0001 2 -1.2135 0.89392
0.0001 3 -1.5282 0.94287
0.0001 4 -1.25 0.8409
0.0001 5 -0.47158 0.66565
0.0001 6 -0.47158 0.655665

A.3 Core masses and dredge-up parameter ${\lambda _{{max}}}$

The core mass at the first thermal pulse, $M_{{\rm c,1TP}}$, is fitted as in Karakas et al. (2002) with an extra set of coefficients for the Z=10-4 models: $p_{1\dots7}=-0.0202032$, 3.33072, 0.847026, 0.0376153, 0.493466, 2.52581, -0.318271. Similarly, the coefficients for $M_{{\rm c,min}}$ for Z=10-4are $a_{1\dots4}=0.788873$, -0.437194, 0.241785, -0.0311767 and for  $\lambda_{{\rm max}}$ are $b_{1\dots4}=-0.823671$, $\min(1,1.25{-}0.1~M/{M_{\odot}})\times0.855046$, 0.0759697, 0.0945916 (both formulae are given in Karakas et al. 2002).

A.4 First dredge up

Abundance changes at first dredge up are interpolated as a function of mass and metallicity from a table of Monash model results.

A.5 Second dredge up

Second dredge-up occurs in sufficiently massive stars ( $M_{{\rm c},{\rm bagb}}\geq0.8~{M_{\odot}}$where  $M_{{\rm c,bagb}}$ is the core mass of the star at the start of the Early-AGB) at the end of the Early-AGB when twin shell burning begins. Following Renzini & Voli (1981) and Groenewegen & de Jong (1993), with alterations to better fit the Monash models, the hydrogen-rich fraction of the envelope is defined as

a=\frac{M-M_{{\rm c,bagb}}}{M-M_{{\rm c}}^{A}}
\end{displaymath} (A.8)

while the fraction which is hydrogen-burned is

b =\frac{M_{{\rm c,bagb}}-M_{{\rm c}}^{A}}{M-M_{{\rm c}}^{A}},
\end{displaymath} (A.9)

where $M_{{\rm c}}^{A}$ is the core mass just after second dredge-up, which is assumed to be equal to the core mass at the first thermal pulse  $M_{{\rm c}}^{{\rm 1TP}}$. Then the surface abundances after second dredge-up X'i are given by

\begin{eqnarray*}X'_{i} = aX_{i,{\rm env}}+bX_{i,{\rm burnt}}

where $X_{i,{\rm env}}$ is the envelope abundance just prior to second dredge up, and  $X_{i,{\rm burnt}}$ is modified to take into account hydrogen burning by converted $^{1}{\rm H}$ to  $^{4}{\rm He}$and all CNO to  $^{14}{\rm N}$, and setting $X_{^{26}{\rm Al},{\rm burnt}}=2.5$ $\times$ 10-5Z.

A.6 Intershell abundances

We calculate intershell abundances by linearly interpolating the Monash model results as a function of core mass, thermal pulse number and metallicity. S-process isotopes are dealt with below.

A.7 The s-process

Our models have masses M=1.5, 3 and $5~{M_{\odot}}$ and metallicities Z=0.02, 0.006, 0.002, $5\times10^{-4}$ and 10-4. The $^{13}{\rm C}$pocket has mass $9.3\times10^{-4}~{M_{\odot}}$ and is assumed to contain no $^{14}{\rm N}$ (which would act as a neutron poison). We define $\xi$ to be proportional to the abundance $X_{{\rm C13}}^{{\rm p}}$of  $^{13}{\rm C}$ in the pocket

\xi=X_{{\rm C13}}^{{\rm p}}/0.00382.
\end{displaymath} (A.10)

We interpolate directly from the detailed model results to obtain intershell abundances of Ag, As, Au, Ba, Bi, Br, Cd, Ce, Co, Cs, Cu, Dy, Er, Eu, Fe, Ga, Gd, Ge, Hf, Hg, Ho, I, In, Ir, Kr, La, Lu, Mo, Nb, Nd, Ni, Os, Pb, Pd, Pm, Po, Pr, Pt, Rb, Re, Rh, Ru, Sb, Se, Sm, Sn, Sr, Ta, Tb, Tc, Te, Tl, Tm, W, Xe, Y, Yb, Zn and Zr as a function of total mass, pulse number and $\xi$. We also have the option of using data directly from the detailed models for all the isotopes available, but this would result in us following more than 400 isotopes in our code which is both computationally intensive and not required for our present work.

A.8 HBB Recalibration

Izzard et al. (2004) calibrated the fraction of the convective envelope exposed to HBB, $f_{{\rm HBB}}$, and the time for which this is burned as a fraction of the interpulse time, $f_{{\rm burn}}$, as a function of M and Z by comparing the surface abundances as a function of time of the synthetic model to the Monash models. As discussed in Izzard et al. (2004) there is some degeneracy between the free parameters. It was worthwhile to repeat the process here because we have more isotopes, some of which are more sensitive to  $T_{{\rm bce}}$, $f_{{\rm HBB}}$ or $f_{{\rm burn}}$, so the parameters can be better constrained. We use the same Levenberg-Marquardt code as Izzard et al. (2004), which is based on the code of Numerical Recipes in C (Press et al. 1992).

We slightly reduce the peak temperature of HBB, $\log_{10}T_{{\rm bce}}$, by $0.18\%$ prior to the shifts given below, compared to the Izzard et al. (2004) fit, to obtain better agreement between the synthetic and Monash models for the Ne, Na, Mg and Al isotopes. This is not entirely unphysical, as $T_{{\rm bce}}$ is taken from the base of the envelope and there will always be HBB at a lower temperature than this just above the base. In the case of the CNO isotopes this could be accounted for by changing $f_{{\rm burn}}$ and/or $f_{{\rm HBB}}$, but for NeNa and MgAl this proved impossible, probably because the appropriate reaction rates are more steeply dependent on temperature. We also refit $N_{{\rm rise}}$, the temperature turn-on factor (see Izzard et al. 2004).

Table A.3 shows the new parameters. We linearly interpolate on this table for a given M and Z (the initial TP-AGB mass is used for M). The final column is $f_{{\rm T}}$, a multiplicative factor on the log temperature after it is reduced by $0.18\%$ - it is reassuring that this factor is usually very close to 1.0.

Table A.3: Our refit of the Izzard et al. (2004) HBB calibration parameters.
$M/{M}_\odot$ Z $N_{{\rm rise}}$ $f_{{\rm burn}}$ $f_{{\rm HBB}}$ $f_{{\rm T}}$
4.0 0.0001 1.68 0.94 0.62 1
4.0 0.004 3.18 0.34 0.585 0.998
4.0 0.008 2.98 0.09 0.09 1
4.0 0.02 3 0.0 0.0 1
5.0 0.0001 0.76 0.13 0.99405 1.0055
5.0 0.004 2.52 1.5 0.64 1.00048
5.0 0.008 2.743 1.33 0.75 1
5.0 0.02 2.52 1.21 0.144 1
6.0 0.0001 1.0619 0.892308 0.716098 1.00059
6.0 0.004 0.16 0.92 0.758 1
6.0 0.008 1.58 1.0413 0.6913 0.99975
6.0 0.02 1.1 1.5 0.54 1
6.5 0.0001 1.19292 0.91236 0.689636 1.00043
6.5 0.004 1.05016 0.947269 0.696531 1.00019
6.5 0.008 1.40533 0.992912 0.66681 1.00007
6.5 0.02 1.55 0.96 0.55 1

Appendix B: NeNa-cycle solution method

The general problem is

$\displaystyle %
\frac{{\rm d}}{{\rm d}t}\mathbf{x}(t) = \mathbf{Ux}(t),$     (B.1)

where the abundances form the vector x(t),
$\displaystyle %
\mathbf{x}(t) = \left(\begin{array}{c}
^{20}{\rm Ne}\\
x_{4}\end{array}\right)$     (B.2)

$\displaystyle %
\mathbf{U} = \left[\begin{array}{cccc}
\frac{-1}{\tau_{20}} & 0...
...}} & 0\\
0 & 0 & \frac{1}{\tau_{22}} & \frac{-1}{\tau_{23}}\end{array}\right].$     (B.3)

In components
$\displaystyle %
\frac{{\rm d}x_{i}(t)}{{\rm d}t}$ = Uijxj(t), (B.4)

where the sum over j is implicit. The matrix  U has 4 eigenvalues  $\lambda_{k}$ and corresponding eigenvectors  U(k)where k=1,2,3,4 is a label over the eigenvectors. Solutions are of the form
$\displaystyle %
x_{i}(t) = \sum_{k=1,2,3,4}A_{k}{\rm e}^{\lambda_{k}t}U_{i}^{(k)}.$     (B.5)

The eigenvalues $\lambda_{k}$ are solutions of
$\displaystyle %
...1}{\tau_{23}}+\lambda\right) -\frac{1}{\tau_{20}\tau_{21}\tau_{22}\tau_{23}}=0,$     (B.6)

one is zero and the others can be calculated with the cubic formula. Substitution gives the eigenvectors

\end{displaymath} (B.7)

where $\alpha=1/\tau_{20}$, $\beta=1/\tau_{21}$, $\gamma=1/\tau_{22}$and $\delta=1/\tau_{23}$. The constants Ak are calculated from the initial abundance xi(t=0) i.e.
$\displaystyle %
x_{i}(t=0) = \sum_{k=1,2,3,4}A_{k}U_{i}^{(k)},$     (B.8)

which can be rewritten as
x(t=0) = UA,     (B.9)

$\displaystyle %
\mathbf{A} = \left(\begin{array}{c}
A_{4}\end{array}\right).$     (B.10)

Pre-multiplying by the inverse of U gives A=U-1xas required, and hence the solution for  x(t).

Appendix C: MgAl-chain solution method

The following is the general solution to the MgAl chain when $^{27}{\rm Al}$acts as a sink. The differential equation set to be solved is then

\frac{{\rm d}^{24}{\rm Mg}}{{\rm d}t}=-\frac{^{24}{\rm Mg}}{\tau_{24}},
\end{displaymath} (C.1)

\frac{{\rm d}^{25}{\rm Mg}}{{\rm d}t}=-\frac{^{25}{\rm {Mg}...
...c{^{25}{\rm {Mg}}}{\tau_{25}}+\frac{^{24}{\rm Mg}}{\tau_{24}},
\end{displaymath} (C.2)

where the reactions $^{25}{\rm Mg}(p,\gamma)^{26}{\rm Al}^{g}$and $^{25}{\rm Mg}(p,\gamma)^{26}{\rm Al}^{m}$ are combined into one timescale, $\tau_{25m}^{-1}+\tau_{25g}^{-1}=\tau_{25}^{-1}$,

\frac{{\rm d}^{26}{\rm Al}}{{\rm d}t}=\frac{^{25}{\rm Mg}}{...
...c{^{25}{\rm Mg}}{\tau_{25g}}-\frac{^{26}{\rm Al}}{\tau'_{26}},
\end{displaymath} (C.3)

where $\tau_{26'}^{-1}+\tau_{\beta26}^{-1}=\tau'_{26}\!^{-1}$,

\frac{{\rm d}^{26}{\rm {Mg}}}{{\rm d}t}=\frac{^{26}{\rm Al}...
...rac{^{25}{\rm Mg}}{\tau_{25m}}-\frac{^{26}{\rm Mg}}{\tau_{26}}
\end{displaymath} (C.4)


\frac{{\rm d}^{27}{\rm Al}}{{\rm d}t}=\frac{^{26}{\rm Mg}}{\tau_{26}}+\frac{^{26}{\rm Al}}{\tau_{26'}}>0,
\end{displaymath} (C.5)

where the isotopes again represent their number densities.

A useful result is the solution to the differential equation

\dot{y}+\lambda y=\alpha+\sum_{i}\beta_{i}{\rm e}^{-\gamma_{i}t}
\end{displaymath} (C.6)

which is

...rac{\beta_{i}}{\lambda-\gamma_{i}}\right){\rm e}^{-\lambda t}.
\end{displaymath} (C.7)

The first equation of the chain has the solution

^{24}{\rm Mg}=a_{24}+b_{24}{\rm e}^{-t/\tau_{24}},
\end{displaymath} (C.8)

where a24=0 (if there is no input to the chain) and b24= $^{24}{\rm Mg}_{0}-a_{24}$. The $^{25}{\rm Mg}$ equation can be arranged as
        $\displaystyle %
\frac{{\rm d}{}^{25}{\rm Mg}}{{\rm d}t}+\frac{^{25}{\rm Mg}}{\tau_{25}}$ = $\displaystyle \frac{^{24}{\rm Mg}}{\tau_{24}}=\frac{a_{24}}{\tau_{24}}+\frac{b_{24}}{\tau_{24}}{\rm e}^{-t/\tau_{24}}$  
  = $\displaystyle a'_{24}+b'_{24}{\rm e}^{-t/\tau_{24}},$ (C.9)

which has the solution

^{25}{\rm Mg}=a_{25}+b_{25}{\rm e}^{-t/\tau_{24}}+c_{25}{\rm e}^{-t/\tau_{25}},
\end{displaymath} (C.10)

where $a_{25}=\tau_{25}a'_{24}$, $b_{25}=b'_{24}/(\frac{1}{\tau_{25}}-\frac{1}{\tau_{24}}$) and $c_{25}=^{25}{\rm Mg}_{0}-b_{25}-a_{25}$. The $^{26}{\rm Al}$ equation is
          $\displaystyle %
\frac{{\rm d}{}^{26}{\rm Al}}{{\rm d}t}+\frac{^{26}{\rm Al}}{\tau_{26}'}$ = $\displaystyle \frac{^{25}{\rm Mg}}{\tau_{25}}$  
  = $\displaystyle a'_{25}+b'_{25}{\rm e}^{-t/\tau_{24}}+c'_{25}{\rm e}^{-t/\tau_{25}},$ (C.11)


\end{displaymath} (C.12)

and $a'_{25}=a_{25}/\tau_{25g}$, $b'_{25}=b_{25}/\tau_{25g}$ and $c'_{25}=c_{25}/\tau_{25g}$. The solution is

^{26}{\rm Al}=a_{26}+b_{26}{\rm e}^{-t/\tau_{24}}+c_{26}{\rm e}^{-t/\tau_{25}}+d_{26}e^{-t/\tau'_{26}},
\end{displaymath} (C.13)

where $a_{26}=\tau'_{26}a'_{25}$, $b_{26}=b'_{25}/(\frac{1}{\tau'_{26}}-\frac{1}{\tau_{24}})$, $c_{26}=c'_{25}/(\frac{1}{\tau'}-\frac{1}{\tau_{25}})$ and $d_{26}=^{26}{\rm Al}_{0}-a_{26}-b_{26}-c_{26}$.

Substitution into the $^{26}{\rm Mg}$ equation gives

$\displaystyle %
\frac{{\rm d}{}^{26}{\rm Mg}}{{\rm d}t}+\frac{^{26}{\rm Mg}}{\t...
+c'_{26}{\rm e}^{-t/\tau_{25}}+d'_{26}{\rm e}^{-t/\tau'_{26}},$     (C.14)

where $a'_{26}=a_{26}/\tau_{\beta26}+a_{25}/\tau_{25m}$, $b'_{26}=b_{26}/\tau_{\beta26}+b_{25}/\tau_{25m}$, $c'=c_{26}/\tau_{\beta26}+c_{25}/\tau_{25m}$ and $d'_{26}=d_{26}/\tau_{\beta26}$. So

^{26}{\rm Mg}\!=\!f_{26}+g_{26}{\rm e}^{-t/\tau_{24}}\!+\!h...
..._{26}{\rm e}^{-t/\tau'_{26}}\!+\!j_{26}{\rm e}^{-t/\tau_{26}},
\end{displaymath} (C.15)

where $f_{26}=\tau_{26}a'_{26}$, $g_{26}=b'_{26}/(\frac{1}{\tau_{26}}-\frac{1}{\tau_{24}})$, $h_{26}=c'_{26}/(\frac{1}{\tau_{26}}-\frac{1}{\tau_{25}})$, $i_{26}=d'_{26}(\frac{1}{\tau_{26}}-\frac{1}{\tau'_{26}})$and $j_{26}=^{26}{\rm Mg}_{0}-f_{26}-g_{26}-h_{26}-i_{26}$. Finally for $^{27}{\rm Al}$

\frac{{\rm d}^{27}{\rm Al}}{{\rm d}t}=\frac{^{26}{\rm Mg}}{\tau_{26}}+\frac{^{26}{\rm Al}}{\tau'_{26}}>0,
\end{displaymath} (C.16)

which is of a slightly different form to the other equations. It can be rewritten as
                      $\displaystyle %
\frac{{\rm d}{}^{27}{\rm Al}}{{\rm d}t}$ = $\displaystyle f'_{26}+g'_{26}{\rm e}^{-t/\tau_{24}}+h'_{26}{\rm e}^{-t/\tau_{25}}+i'_{26}{\rm e}^{-t/\tau'}+j'_{26}{\rm e}^{-t/\tau_{26}}$  
    $\displaystyle +a_{26}^{*}+b_{26}^{*}{\rm e}^{-t/\tau_{24}}+c_{26}^{*}{\rm e}^{-t/\tau_{25}}+d_{26}^{*}{\rm e}^{-t/\tau'_{26}},$ (C.17)

where $f'_{26}=f_{26}/\tau_{26}$, $g'_{26}=g_{26}/\tau_{26}$, $h'_{26}=h_{26}/\tau_{26}$, $i'_{26}=i_{26}/\tau_{26}$, $j'_{26}=j_{26}/\tau_{26}$, $a_{26}^{*}=a_{26}/\tau'_{26}$, $b_{26}^{*}=b_{26}/\tau'_{26}$, $c_{26}^{*}=c_{26}/\tau'_{26}$ and $d_{26}^{*}=d_{26}/\tau'_{26}$. This can be further simplified to

\frac{{\rm d}{}^{27}{\rm Al}}{{\rm d}t}\!=\!a_{27}+b_{27}{\...
..._{27}{\rm e}^{-t/\tau_{26}}\!+\!f_{27}{\rm e}^{-t/\tau'_{26}},
\end{displaymath} (C.18)

where a27=f'26+a26*, b27=g'26+b26*, c27=h'26+c26*, d27=j'26 and f27=i'26+d26*. Integration gives
                      $\displaystyle %
^{27}{\rm Al}$ = $\displaystyle a_{27}t-(b_{27}\tau_{24}{\rm e}^{-t/\tau_{24}}+c_{27}\tau_{25}{\rm e}^{-t/\tau_{25}}$  
    $\displaystyle +d_{27}\tau_{26}{\rm e}^{-t/\tau_{26}}+f_{27}\tau'{\rm e}^{-t/\tau'_{26}}){\rm ~+~ const.},$ (C.19)

where the constant of integration is found from the initial value  $^{27}{\rm Al}_{0}$. The full solution is
                      $\displaystyle %
^{27}{\rm Al}$ = $\displaystyle {}^{27}{\rm Al}_{0}\!+\!a_{27}t-(b_{27}\tau_{24}{\rm e}^{-t/\tau_...
...c_{27}\tau_{25}{\rm e}^{-t/\tau_{25}}\!+\!d_{27}\tau_{26}{\rm e}^{-t/\tau_{26}}$  
    $\displaystyle +f_{27}\tau'{\rm e}^{-t/\tau'_{26}}) \!+\!(b_{27}\tau_{24}\!+\!c_{27}\tau_{25}\!+\!d_{27}\tau_{26}\!+\!f_{27}\tau').$ (C.20)

After the calculation is complete, we allow for an explicit exponential decay of $^{27}{\rm Al}$ to  $^{24}{\rm Mg}$, but the rate is negligible.

Appendix D: Comparison of the synthetic and Monash models

In Figs. D.1-D.3 we compare the detailed Monash models for hot-bottom burning stars ( $M/{\rm ~ M_{\odot}}=4,~5$ and 6) to our synthetic models.

\par\includegraphics[width=14cm,clip]{6129figD1}\end{figure} Figure D.1: Comparison of surface abundances (log mass fraction) vs. time (in units of $10^{5}~{\rm years}$) for $^{12}{\rm C}$to  $^{20}{\rm Ne}$ for our synthetic models (solid lines) and the Monash detailed models (dashed lines) with Z=0.004, 0.008and 0.02. The three columns represent M=4, 5 and $6~{M_{\odot}}$respectively from left to right.

\par\includegraphics[width=14cm,clip]{6129figD2}\end{figure} Figure D.2: As Fig. D.1 for $^{21}{\rm Ne}$to $^{26}{\rm Mg}$.

\par\includegraphics[width=14cm,clip]{6129figD3}\end{figure} Figure D.3: As Fig. D.1 for $^{26}{\rm Al}$and $^{27}{\rm Al}$.

Appendix E: Belczynski NS/BH mass prescription

The Belczynski et al. (2002) NS/BH mass is coded according to a prescription provided by Jarrod Hurley (private communication).

First, $M_{{\rm cx}}$ is set

$\displaystyle %
M_{{\rm cx}} = \left\{ \begin{array}{cc}
0.161767M_{{\rm c}}+1....
...314154M_{{\rm c}}+0.686088, & M_{{\rm c}}\geq2.5~{M_{\odot}},\end{array}\right.$     (E.1)

then the remnant mass is given by
$\displaystyle %
M_{{\rm NS/BH}} = \left\{ \begin{array}{cc}
M_{{\rm cx}}, & M_{...
...eq M_{{\rm c}}/{M_{\odot}<7.6},\\
M, & M\geq7.6~{M_{\odot}}.\end{array}\right.$     (E.2)

Appendix F: Mass-loss prescriptions

The following mass-loss prescriptions are included here for completeness. No justification to the terms is given, for such details see Hurley et al. (2002), Dray et al. (2003) and references included below.

F.1 SSE (Hurley et al. 2002)

F.1.1 GB and post-GB stars

The formula of Kudritzki & Reimers (1978) is used:

\dot{M}_{{\rm R}}=\eta~4\times10^{-13}\frac{(L/{L_{\odot}})(R/{R_{\odot}})}{(M/{M_{\odot}})}~{M_{\odot}}~{\rm yr}^{-1},
\end{displaymath} (F.1)

with $\eta=0.5$.

F.1.2 AGB stars

The Vassiliadis & Wood (1993) rate

\log\dot{M}_{{\rm VW}}\!=\!-11.4\!+\!0.0125\left[P_{0}\!-\!100\max\left(M/{M_{\odot}}\!-\!2.5,0.0\right)\right]
\end{displaymath} (F.2)

is applied where P0 is the Mira period pulsation given by
$\displaystyle %
\log(P_{0}/{\rm d})=\min\left[3.3,-2.07-0.9\log(M/{M_{\odot}})
+1.94\log(R/R_{\odot})\right].$     (F.3)

A cap (the $\min\left[3.3,...\right]$ term) is applied at $\dot{M}_{{\rm VW}}=1.36\times10^{-9}(L/{L_{\odot}}){~ M_{\odot}~{\rm yr}^{-1}}$to model the superwind phase.

F.1.3 Massive stars

The rates of Nieuwenhuijzen & de Jager (1990) are applied for $L>4000~{L_{\odot}}$,

$\displaystyle %
\dot{M}_{{\rm NJ}}=\left(\frac{Z}{Z_{\odot}}\right)^{0.5}9.6\ti...
...0.81} ~ (L/L_{\odot})^{1.24}(M/{M_{\odot}})^{0.16}{\; M_{\odot}~{\rm yr}}^{-1},$     (F.4)

modified by the factor Z0.5.

F.1.4 Small envelopes


\mu=\frac{M_{{\rm env}}}{M}\min\left(5.0,\max\left[1.2,\left\{ \frac{L}{L_{0}}\right\} ^{\kappa}\right]\right),
\end{displaymath} (F.5)

where $L_{0}=7.0\times10^{4}$ and $\kappa=-0.5$ for hydrogen-rich envelopes. For $\mu<1.0$ a Wolf-Rayet-like mass loss is included in the form

\dot{M}_{{\rm WR}}=10^{-13}L^{1.5}(1.0-\mu){~ M_{\odot}}~{\rm yr}^{-1}.
\end{displaymath} (F.6)

F.1.5 Finally...

The total mass-loss rate is the dominant rate from the above choices

\dot{M}=\max(\dot{M}_{{\rm R}},\dot{M}_{{\rm VW}},\dot{M}_{{\rm NJ}},\dot{M}_{{\rm WR}}).
\end{displaymath} (F.7)

A Luminous-Blue-Variable-like mass loss is added for stars beyond the Humphreys-Davidson limit (Humphreys & Davidson 1994),
$\displaystyle %
\dot{M}_{{\rm LBV}}=0.1(10^{-5}(R/{R_{\odot}})(L/{L_{\odot}})^{...
...\left(\frac{L}{6\times10^{5}~{L_{\odot}}}-1.0\right){~ M_{\odot}~{\rm yr}^{-1}}$     (F.8)

which is added to $\dot{M}$ if $L>6\times10^{5}~{L_{\odot}}$and $10^{-5}(R/{R_{\odot}})(L/{L_{\odot}})^{0.5}>1.0$. For naked helium stars the WR mass-loss rate is used with $\mu=0$to give

\dot{M}=\max\left(\dot{M}_{{\rm R}},\dot{M}_{{\rm WR}}(\mu=0)\right).
\end{displaymath} (F.9)

F.2 MM

The MM rates are similar to the enhanced mass-loss rates of Maeder & Meynet (1994).

F.2.1 Pre-WR evolution

The empirical mass-loss rate of de Jager et al. (1988) is used, but note the extra factor of 2,

$\displaystyle %
-\log(-\dot{M}/{M_{\odot}~{\rm yr}^{-1}})=\sum_{n=0}^{N}\sum_{1...
...m K})-4.05}{0.75}\right)
T_{j}\left(\frac{\log(L/{L_{\odot}})-4.6}{2.1}\right),$     (F.10)


T_{j}(x)=2\cos(j~\arccos x)
\end{displaymath} (F.11)

and aij are tabulated in de Jager et al. (1988).

F.2.2 WNL phase

A constant rate of $8\times10^{-5}~{M_{\odot}~{\rm yr}^{-1}}$is used (Conti et al. 1988).

F.2.3 WNE, WC and WO phases

The theoretical mass-loss rate of Langer (1989),

\dot{M}_{{\rm WR}}=0.6\times10^{-7}\left(\frac{M_{{\rm WR}}}{{M_{\odot}}}\right)^{2.5}~{M_{\odot}}~{\rm yr}^{-1},
\end{displaymath} (F.12)

is used (with a pre-multiplier of 0.6).

F.3 NL

F.3.1 Pre-WR evolution

As with the MM mass-loss rate the de Jager et al. (1988) rates are used, although without the factor of 2 in Eq. (F.11).

F.3.2 WNL and WNE phases

The mass-loss rate of Nugis & Lamers (2000) is used:

\log(\dot{M}_{{\rm WN}}/{M_{\odot}~{\rm yr}^{^{-1}})\!=\!-13.6+1.63\log(L/}L_{\odot})\!+\!2.22\log Y
\end{displaymath} (F.13)

where Y is the surface helium abundance (by mass fraction).

F.3.3 WC phase

The Nugis & Lamers (2000) mass-loss rate

$\displaystyle %
\log(\dot{M}_{{\rm WC}}/{M_{\odot}}~{\rm yr}^{-1})=-8.3+0.84\log(L/{L_{\odot}})
+2.04\log Y+1.04\log Z$     (F.14)

is used.

F.3.4 WO phase

A constant rate of $1.9\times10^{-5}~{M_{\odot}}~{\rm yr}^{-1}$is used.

Appendix G: Stellar mergers

Table G.1 shows the outcome of stellar mergers according to the prescription given in Sect. 4.4.

Table G.1: Collision matrix for the product of a stellar merger, see Hurley et al. 2002 for a definition of the stellar types, and Sect. 4.4 for the meaning of the resulting numbers.
  Primary stellar type
  0-6 7-9 10 11, 12 13, 14  
  0-6 3 2 2 2 4
Secondary stellar type 7-10 1 3 3 2 4
  11, 12 1 1 1 3 4
  13, 14 4 4 4 4 4

Appendix H: Binary star explosions

The coefficients for the fits to the alpha elements from Livne & Arnett (1995) are in Table H.1. The yields are fitted to functions of the form $\Delta M(a+b M_{{\rm CO}}+c M_{{\rm CO}}^{2})\times(1+d M_{{\rm He}})$where $\Delta M$ is the amount of mass ejected in the explosion and the coefficients are given in Table H.1. COWDs that accrete hydrogen-rich matter are treated in the same way, with the hydrogen steadily burnt to helium then CO on the COWD surface prior to the explosion.

Table H.1: Coefficients to the fits to the SNIa yields of Livne & Arnett (1995).
Isotope a b c d
$^{4}{\rm He}$ $1.885\times10^{-2}$ $1.919\times10^{-1}$ $-1.90040\times10^{-1}$ $7.398\times10^{-1}$
$^{12}{\rm C}$ $3.2448\times10^{-1}$ $-7.5668\times10^{-1}$ $4.4693\times10^{-1}$ -2.5692
$^{16}{\rm O}$ $7.9324\times10^{-1}$ -1.4184 $6.7159\times10^{-1}$ -2.2394
$^{20}{\rm Ne}$ $5.6865\times10^{-2}$ $-1.1574\times10^{-1}$ $6.0577\times10^{-2}$ -2.4772
$^{24}{\rm Mg}$ $6.8268\times10^{-2}$ $-1.8269\times10^{-2}$ $-5.339\times10^{-2}$ -1.9527
$^{28}{\rm Si}$ $1.9537\times10^{-2}$ $4.9225\times10^{-1}$ $-4.2238\times10^{-1}$ -1.229
$^{32}{\rm S}$ $3.3483\times10^{-1}$ $-7.4649\times10^{-1}$ $5.2431\times10^{-1}$ -1.0288
$^{36}{\rm Ar}$ $1.553\times10^{-1}$ $-4.2591\times10^{-1}$ $3.0771\times10^{-1}$ $5.0104\times10^{-1}$
$^{40}{\rm Ca}$ $4.7899\times10^{-2}$ $-1.3997\times10^{-1}$ $1.0286\times10^{-1}$ $3.8396\times10^{1}$
$^{44}{\rm Ti}$ $1.5964\times10^{-1}$ $-3.9855\times10^{-1}$ $2.5757\times10^{-1}$ -2.4052
$^{48}{\rm Cr}$ $1.2483\times10^{-1}$ $-3.1778\times10^{-1}$ $2.0946\times10^{-1}$ -2.2352
$^{52}{\rm Fe}$ $3.5674\times10^{-1}$ $-9.9596\times10^{-1}$ $7.0243\times10^{-1}$ $-2.2407\times10^{-1}$
$^{56}{\rm Fe}$ $-7.9924\times10^{-1}$ 2.2283 -1.2529 $1.3922\times10^{1}$

Copyright ESO 2006