The Astrophysical Journal Supplement Series, 129:377-398, 2000 July
© 2000. The American Astronomical Society. All rights reserved. Printed in U.S.A.

An Inexpensive Nuclear Energy Generation Network for Stellar Hydrodynamics

F. X. Timmes
Center on Astrophysical Thermonuclear Flashes, University of Chicago, 933 East 56th Street, Chicago, IL 60637
R. D. Hoffman
Lawrence Livermore National Laboratory, Livermore, CA 94550
and
S. E. Woosley
Lawrence Livermore National Laboratory, Livermore, CA 94550, and Department of Astronomy and Astrophysics, UCO/Lick Observatory, University of California at Santa Cruz, Santa Cruz, CA 95064

Received 1999 November 9; accepted 2000 February 18

ABSTRACT

We compare the nuclear energy generation rate and abundance levels given by an -chain nuclear reaction network that contains only seven isotopes with a standard 13 isotope -chain reaction network. The energy generation rate of these two small networks are also compared to the energy generation rate given by a 489 isotope reaction network with weak reactions turned on and off. The comparison between the seven isotope and -chain reaction networks indicate the extent to which one can be replaced by the other, and the comparison with the 489 isotope reaction network roughly indicates under what physical conditions it is safe to use the seven isotope and -chain reaction networks. The seven isotope reaction network reproduces the nuclear energy generation rate of the standard -chain reaction network to within 30%, but often much better, during hydrostatic and explosive helium, carbon, and oxygen burning. It will also provide energy generation rates within 30% of an -chain reaction network for silicon burning at densities less than 10⁷ g cm^-3. Provided there remains an equal number of protons and neutrons (Y_e = 0.5) over the course of the evolution, and that flows through -particle channels dominate, then both of the small reaction networks return energy generation rates that are compatible with the energy generation rate returned by the 489 reaction network. If Y_e is significantly different from 0.5, or if there are substantial flows through neutron and protons channels, then it is not generally safe to employ any -chain based reaction network. The relative accuracy of the seven isotope reaction network, combined with its reduction in the computational cost, suggest that it is a suitable replacement for -chain reaction networks for parameter space surveys of a wide class of multidimensional stellar models.

Subject headings: hydrodynamics; methods: numerical; nuclear reactions, nucleosynthesis, abundances; stars: general

On-line material: source code

1. INTRODUCTION

     Integration of the ordinary differential equations which represent the abundance levels of a set of isotopes serves two functions in models of stellar events. The primary function, as far as the hydrodynamics is concerned, is to provide the magnitude and sign of the thermonuclear energy generation rate. This thermonuclear energy release is usually the largest source of energy in regions conducive to nuclear reactions, which makes its determination vital to accurate hydrodynamic modeling. The second function that the integration serves for the hydrodynamics is to describe the evolution of the abundance levels of the nuclear species. These abundance levels are, of course, fundamental to our understanding of the origin and evolution of the elements.

     Obtaining accurate values for the thermonuclear energy generation rate is, unfortunately, expensive in terms of computer memory and CPU time. The largest block of memory in a stellar hydrodynamic program is usually reserved for storing the isotopic abundances at every grid point. This memory requirement can be quite restrictive for three-dimensional models on present parallel computer architectures. Even with modern methods for solving the associated stiff system of ordinary differential equations, evolving the isotopic abundances begins to dominate the total cost of a multidimensional model when the number of isotopes evolved is 30. To decrease the computer time and memory it takes to calculate a model means making a choice between having fewer isotopes in the reaction network or having less spatial resolution. The general response to this tradeoff has been to evolve a limited number of isotopes, and thus calculate an approximate thermonuclear energy generation rate. The set of 13 nuclei most commonly used for this purpose are ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti, ⁴⁸Cr, ⁵²Fe, and ⁵⁶Ni. This minimal set of nuclei, usually called an -chain network, can reasonably track the abundance levels from helium burning through nuclear statistical equilibrium. More importantly from a hydrodynamics standpoint, and a point that we will demonstrate in this paper, is that an -chain reaction network gives a thermonuclear energy generation rate that is generally, but not always, within 20% of the energy generation rate given by much larger nuclear reaction networks. In essence, one gets most of the energy generated for most thermodynamic conditions at a fraction of the computational cost (memory + CPU). Can the number of isotopes be further reduced, and still give relatively accurate energy generation rates?

     The primary purpose of this paper is to evaluate the accuracy, memory footprint, and execution speed of an -chain reaction network that evolves only seven isotopes; ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, and ⁵⁶Ni. We call this reduced -chain reaction network the iso7 reaction network. The iso7 network is shown to give a good representation of nuclear energy generation rates during helium, carbon, neon, and oxygen burning. It even gives reasonable estimates for the energy generation rates during silicon burning and photodisintegration reactions. Thus, the iso7 reaction network may be useful for exploratory multidimensional calculations where large reaction networks are impractical even on the largest parallel supercomputers.

     The second purpose of this paper is to assess when it is valid to use the iso7, or any -chain, reaction network. The energy generation rate of iso7 and a standard 13 isotope -chain are compared to the energy generation rate given by a 489 isotope reaction network with and without the Fuller, Fowler, & Neumann (1985) weak reactions turned on. The 489 isotope reaction network extends from ¹H to ⁹¹Tc (Timmes 1999). This large reaction network should be sufficient to assess the validity of the energy generation rate returned by the iso7 and -chain reaction networks.

     In § 2 some generalities about nuclear reaction networks, and the details of the iso7 reaction network are discussed. In § 3 the thermonuclear energy generation rate given by the iso7 network is compared to the energy generation rates given by an -chain reaction network for various burning modes in a number of different fuels. A summary discussion is given in § 4, and in Appendix A we list the key FORTRAN routines that implement the iso7 reaction network.

2. NUCLEAR REACTION NETWORKS IN GENERAL

     Let isotope i have Z_i protons, A_i nucleons (protons + neutrons), and a binding energy B_i (in ergs). Let the aggregate total of isotope i have a mass density _i (in g cm^-3 and a number density n_i (in cm^-3) in material with temperature T (in K) and total mass density . Define the dimensionless mass fraction of isotope i as X_i = _i/ = n_iA_i/(N_A), where N_A is Avogadro's number, and define the dimensionless molar abundance of isotope i as Y_i = X_i/A_i = n_i/(N_A).

     The mean number of nucleons per nucleus is defined as = ( X_i/A_i)^-1, the mean charge per nucleus is defined as = Z_iX_i/A_i, and the number of electrons per baryon is defined as Y_e = /. If the composition is pure protons then Y_e = 1, if the composition is pure neutrons then Y_e = 0, and if there are an equal number of protons and neutrons then Y_e = 0.5.

     The general continuity equation for isotope i in a Lagrangian formulation is



In this set of partial differential equations _i is the total reaction rate due to all binary reactions of the form i(j,k)l



where _kj and _jk are the reverse (creation) and forward (destruction) nuclear reaction rates, respectively. Contributions from one-body reactions, such as -decays, and three-body reactions, such as the triple- reaction, are easy to append to equation (2). The mass diffusion velocities _i in equation (1) are obtained from the solution of a multicomponent diffusion equation (Chapman & Cowling 1970; Burgers 1969; Williams 1988), and reflect the fact that mass diffusion processes respond to pressure, temperature, and abundance gradients as well as external gravitational or electrical forces.

     There are two reasons for considering the case of setting the mass diffusion velocities to zero. The first reason is that the physics of the situation may dictate that mass diffusion processes are small over the timescales of interest when compared to other transport process such as thermal diffusion or viscous diffusion (i.e., large Lewis numbers and/or small Prandtl numbers). The second reason is that one may wish to apply the numerical technique of operator splitting. In operator splitting, the various pieces of the physics are decoupled from one another over a given time step. Each piece is solved for independently of the others, with the various pieces being appropriately summed at the end of the time step. This approach is quite common in spherically symmetric stellar evolution programs and in multidimensional hydrodynamic programs. For example, the properties of an adiabatic fluid flow are determined for a given time step. The energy generated and changes in the composition due to thermonuclear reactions over the time step are independently determined. The effects of mass diffusion processes may also be obtained independently. These individual pieces are then summed. Potential disadvantages of the operator splitting technique are that the point at which the decoupling becomes a poor approximation is difficult to estimate a priori, and the ordering of the operators may not be commutative. For either reason, because the physics allows it or because operator splitting is being used, setting the mass diffusion velocity to zero transforms equation (1) into



This set of ordinary differential equations are what constitute a "reaction network."

     The coefficients which appear in the ordinary differential equations of equation (3) span many orders of magnitude because the thermonuclear reaction rates have highly nonlinear dependences on the temperature, incorporate vastly different types of reactions (e.g., neutron capture, charged particle, and weak reactions), and because the abundance levels themselves typically range over several orders of magnitude (e.g., the solar isotopic abundances). As a result, the ordinary differential equations that comprise a nuclear reaction network are "stiff." A set of equations are stiff, in a mathematical context, when the real part of the maximum to minimum eigenvalue ratio of the Jacobian matrix = / is much larger than unity. It is not uncommon for the nuclear reaction networks of astrophysics to have real parts of this eigenvalue ratio that are 10¹⁵. From a physics point of view, this means that at least one of the abundance levels is changing on a much faster timescale than the others. From a practical standpoint, a system of stiff equations usually means that an implicit time integration, hence having the time-dependent Jacobian matrix, is necessary.

     Methods for evolving the stiff system of ordinary differential equations that constitute nuclear reaction networks are surveyed in Timmes (1999). The scaling properties and behavior of three semi-implicit time integration algorithms (a traditional first-order accurate Euler method, a fourth-order accurate Kaps-Rentrop method, and a variable order Bader-Deuflhard method) and eight linear algebra packages (LAPACK, LUDCMP, LEQS, GIFT, MA28, UMFPACK, and Y12M) were investigated by running each of these 24 combinations on seven different nuclear reaction networks (a hardwired 13 isotope network, a hardwired 19 isotope network, a 47 isotope, 76 isotope, 127 isotope, 200 isotope, and 489 isotope reaction networks). The analysis permitted a full assessment of the integration methods and linear algebra packages surveyed, and several pragmatic suggestions to be offered. When the best balance between accuracy, overall efficiency, memory footprint, and ease of use is desirable, Timmes (1999) suggested that the variable order Bader-Deuflhard time integration method coupled with the MA28 sparse matrix package is a good choice. These choices are adopted for this paper, but do not significantly effect the main results.

     Once a nuclear reaction network has been constructed, the instantaneous nuclear energy generation rate _nuc follows from manipulating E = mc² to give



A nuclear reaction network does not need to be evolved in order to obtain the instantaneous energy generation rate, since only the right-hand sides of the ordinary differential equations need be evaluated (e.g., eq. [3]). In an operator split hydrodynamics program it is usually more appropriate to use the average energy generated over a time step t



In this case, a nuclear reaction network does need to be evolved. The average energy generation rate is used in the results below, and similar results are obtained if the instantaneous energy generation rate is used.

2.1. The iso7 Reaction Network

     The iso7 reaction network is based on the nine isotope reaction network of Woosley (1986), the main difference being the omission of the neutron and proton species from the network, principally because neither one contributes to the energy generation rate (both have zero binding energy). This omission does forego the option of including weak reactions on the nucleons, which are important for the late stages of Type II supernova core collapse scenarios, but seems justified for a wider variety of applications. On the positive side this omission reduces the memory required in a multidimensional code, and reduces the execution time to solve the reaction network equations. In addition, the iso7 reaction network improves some of the branching ratio linkages of the heavy ion reactions (in particular, ¹²C + ¹⁶O), clarifies the coding, and splits the formation of the right-hand sides and the Jacobian matrix into two separate routines (useful for using modern methods of solving stiff ordinary differential equations). It also simplifies the conditions for turning on silicon burning. The three FORTRAN routines which implement the iso7 reaction network are listed in Appendix A.

     Note that neither the iso7 reaction network or the nine-isotope reaction network, nor any -chain reaction network, can simulate hydrogen burning (PP chains, CNO cycles, hot CNO cycles, or the rp-process). Such -chain networks are a very poor choice for modeling the hydrogen burning phases of main-sequence stars, classical novae, and some types of X-ray burst. These networks are also of very limited use for neutron-rich environments, such as in the central regions during the advanced stages of stellar evolution in massive stars, the inner regions of core collapse supernovae, and the central regions of Type Ia supernovae at late times. However, the iso7 reaction network offers a reasonable approximation for the energy generation rate of helium, carbon, neon, and oxygen burning under normal stellar evolution, and suitable for parameter space surveys of a wide class of multidimensional stellar models.

     From ⁴He to ²⁴Mg (5 isotopes) the iso7 reaction network is identical to a standard -chain reaction network. If silicon burning or photodisintegration reactions are not occurring, then the iso7 reaction network gives the same nucleosynthesis and energy generation rate as an -chain reaction network.

     The iso7 reaction network is different from an -chain reaction network starting at ²⁸Si. As in an -chain reaction network, the photodisintegration flow of ²⁸Si down to ²⁴Mg is accounted for by a (,) reaction. Unlike an -chain reaction network, there is no flow for moving ²⁸Si into ³²S and beyond by (,) reactions since ³²S is absent from the iso7 reaction network. Instead, the iso7 reaction network moves ²⁸Si plus 7 ⁴He nuclei directly into ⁵⁶Ni, skipping the intermediate isotopes present in an -chain reaction network. The iso7 reaction network makes this direct transition from ²⁸Si to ⁵⁶Ni at a rate given by the product of an effective ⁴⁰Ca abundance and the ⁴⁰Ca(,) reaction rate. Similarly, ⁵⁶Ni photodisintegrates directly into ²⁸Si plus 7 ⁴He nuclei, again skipping six isotopes present in an -chain reaction network, at a rate given by the product of an effective ⁴⁴Ti abundance and the ⁴⁴Ti(,) reaction rate. In effect, the "²⁸Si" abundance in the iso7 reaction network is a silicon group abundance, and the "⁵⁶Si" abundance is an iron group abundance.

     Flows between the silicon group nuclei, iron group nuclei, and -particles are modeled as follows:



where the notation of equations (2) and (3) have been applied. Under the assumption that quasi-equilibrium states exist among and between the silicon and iron group abundances during silicon burning (Bodansky, Clayton, & Fowler 1968; Woosley, Arnett, & Truran 1973; Hix & Thielemann 1996), the ⁴⁰Ca and ⁴⁴Ti abundances are related to the ²⁸Si and ⁵⁶Ni abundances by



where T₉ = T/(10⁹ K). The functions f(T₉) and g(T₉) are branching ratios of the appropriate reaction rates, and approximated in the iso7 reaction network as

     These direct transitions between ²⁸Si and ⁵⁶Ni are the key steps in reducing the network size to seven isotopes. They are also the weak point of the iso7 reaction network. How well these approximations mimic the energy generation rate and nucleosynthesis of an -chain reaction network occupies most of our ensuing discussion about the results. Since the iso7 reaction network is substantially different from an -chain reaction network only for isotopes heavier than ²⁴Mg, the first question to answer is how well the iso7 reaction network models various types of silicon burning.

3. ANALYSIS OF THE iso7 REACTION NETWORK

     Two burning modes are used to analyze the results of iso7 reaction network. The first type of burning occurs at a constant temperature and density, and is referred to below as "hydrostatic burning" because hydrostatic burning in stars occurs under slowly evolving thermodynamic conditions. This mode of burning is very common in operator split hydrodynamic programs; over a given time step the adiabatic properties of the fluid, and the results of a hydrostatic burn are combined. The second type of burning models a region where a shock has heated and compressed the material to some peak temperature T₀ and density ₀, which subsequently expands adiabatically as a radiation-dominated gas (Fowler & Hoyle 1964). This adiabatic expansion is referred to below as "explosive burning." The temperature and density decline in time over a hydrodynamic timescale (a free-fall timescale):

     The results given by the iso7 reaction network are compared to the results given by an -chain reaction network in the analysis below. A definition of what we mean by an -chain reaction network is prudent. A strict -chain reaction network is only composed of (,) and (,) links among the 13 isotopes ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti, ⁴⁸Cr, ⁵²Fe, and ⁵⁶Ni. It is essential, however, to include (,p)(p,) and (,p)(p,) links in order to obtain reasonably accurate energy generation rates and abundance levels when the temperature exceeds 2.5 × 10⁹ K. At these elevated temperatures the flows through the (,p)(p,) sequences are faster than the flows through the (,) channels. An (,p)(p,) sequence is, effectively, an (,) reaction through an intermediate isotope. In our -chain reaction network, we include eight (,p)(p,) sequences plus the corresponding inverse sequences by assuming steady-state proton flows through the intermediate isotopes ²⁷Al, ³¹P, ³⁵Cl, ³⁹K, ⁴³Sc, ⁴⁷V, ⁵¹Mn, and ⁵⁵Co. This strategy permits inclusion of (,p)(p,) sequences without explicitly evolving the proton or intermediate isotope abundances. Thus, our -chain reaction network includes not just (,) and (,) links, but also links through the (,p)(p,) and (,p)(p,) sequences. The iso7 reaction network does not include such sequences (eq. [6]) since it jumps directly from ²⁸Si to ⁵⁶Ni.

3.1. Hydrostatic and Explosive Silicon Burning

     An initially pure silicon composition at a temperature of 6 × 10⁹ K and density of 1 × 10⁷ g cm^-3 was evolved under hydrostatic conditions for 100 s. The thermonuclear energy generation rate as a function of time is shown in Figure 1a. The x-axis gives the elapsed time on a logarithmic scale, and the y-axis gives the energy generation rate on a logarithmic scale. The solid red curve represents the values given by the iso7 reaction network, red circles for the -chain reaction network, gold triangles for the 489 isotope network with Fuller et al. (1985) weak reactions turned off, purple diamonds for the 489 isotope reaction network with Fuller et al. (1985) weak reactions turned on, and green squares show the magnitude of the Itoh et al. (1996) neutrino energy losses. Silicon burning at these constant thermodynamic conditions is an endothermic process because of the relative strength of the (,) photodisintegration reactions, resulting in a negative energy generation rate (an energy sink), labeled as such on the plot. Figure 1a suggests that the iso7 and -chain reaction networks give nuclear energy generation rates that follow each other, at least on energy generation rate scales that span many orders of magnitude.

Fig. 1a Fig. 1b Fig. 1c

Fig. 1   (a) Comparison of the energy generation rate under the hydrostatic conditions of T = 6 × 109 K, = 1 × 107 g cm-3, and an initially pure 28Si composition. The solid red curve is for the iso7 network, red circles for the -chain network, gold triangles for the 489 isotope network with weak reactions turned off, purple diamonds for the 489 isotope network with weak reactions turned on, and the green squares for the magnitude of the neutrino energy-loss rate. For these conditions, silicon burning is endothermic, and the two 489 isotope reaction networks give nearly identical results. (b) Ratio of energy generation rate given by the 489 isotope reaction network, with the Fuller et al. (1985) weak reaction rates turned on, to the -chain (green curve) and iso7 (purple curve) reaction networks. The thermodynamic conditions are the same as in (a). The error between either of these two small networks and the full networks is larger than the error between the iso7 and -chain networks. (c) Evolution of the mass fractions under the hydrostatic conditions of T = 6 × 109 K, = 1 × 107 g cm-3, and an initially pure 28Si composition. Each abundance curve is colored according to the legend at the right of the figure. Solid curves are for the iso7 reaction network, and open circles for the 13 isotope -chain network. For the -chain network, the sum of the 28Si, 32S, 36Ar, 40Ca, and 44Ti mass fractions are plotted as silicon, and the sum of the 48Cr, 52Fe, and 56Ni mass fractions are plotted as nickel.

     Throughout the evolution in Figure 1a there remains an equal number (to four significant figures) of protons and neutrons, Y_e = 0.5, since insufficient time elapses for weak reactions to change Y_e. This is why the curves for the 489 isotope reaction network with weak reactions turned off (gold triangles) and with weak reactions turned on (purple diamonds) overlay each other. This is also why the iso7 and the 13 isotope -chain reaction networks give energy generation rates that are with 30% of the energy generation rates returned by the 489 isotope reaction networks. For the constant thermodynamic conditions being considered, Figure 1a suggests that it is generally safe, for calculation of the energy generation rate, to replace a large reaction network with a smaller -chain reaction network.

     A comparison of the energy generation rates on linear scales is shown in Figure 1b. The x-axis gives the elapsed time on a logarithmic scale, and the y-axis gives the ratio of the energy generation rate given by the 489 isotope reaction network to the energy generation rate given by the -chain (green curve) and the iso7 (purple curve) reaction networks. If the two small reaction networks give the same energy generation rate at the large reaction network, then the ordinate would be zero. If the energy generation rates differed by 10%, then the ordinate would be 0.1. The linear scales of Figure 1b show that the largest source of error is made by choosing one of the small -chain networks instead of a large reaction network. While the energy generation rate given by the iso7 and -chain reaction networks agree to within 30% over the entire evolution, the error associated with not using a large reaction network is larger.

     Evolution of the nuclear mass fractions is shown in Figure 1c. The x-axis gives the elapsed time on a logarithmic scale, and the y-axis gives the mass fraction on a logarithmic scale. The isotopes plotted are labeled and colored according to the legend given to the right of the figure. The solid curves represent the evolution as determined by the iso7 reaction network, and the open circles represent the evolution of the mass fractions as given by the -chain reaction network. The mass fractions of ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, and ⁴⁴Ti in the -chain reaction network are summed, with the total plotted as the "silicon group." Similarly, the mass fractions of the ⁴⁸Cr, ⁵²Fe, and ⁵⁶Ni isotopes in the -chain network are summed and plotted as the "nickel group."

     Figure 1c indicates that the mass fractions as determined by the iso7 reaction network are compatible with the mass fractions determined by the -chain reaction network. There are small differences in the nickel group before 10^-4 s. By 10^-2 s the composition has reached its final distribution, which is a state of nuclear statistical equilibrium. It should be noted that nuclear statistical equilibrium abundances given by a network which only contains -particles and -chain nuclei can differ substantially from the nuclear statistical equilibrium distribution given by a reaction network that includes nuclei with / 0.5 even for compositions with Y_e = 0.5.

     There is about a factor of 3 difference between the final silicon group abundances as given by the iso7 and -chain reaction networks in Figure 1c. It might seem surprising that such a relatively large difference has only small effects on the nuclear energy generation rates (Figs. 1a and 1b). Differences in small abundances have little effect on the nuclear energy released because the nuclear energy generation rate depends linearly on the abundance changes (eq. [5]). It is the largest abundance changes that dominate the energy release (eq. [5]), not the largest abundances.

     An initially pure silicon composition that starts with a temperature of 6 × 10⁹ K and a density of 1 × 10⁷ g cm^-3 was evolved under the adiabatic expansion approximation of equation (9), and the results are shown in Figures 2a, 2b, and 2c. The format, symbols, and color scheme of these plots are the same as those used for the corresponding plots of Figure 1.

Fig. 2a Fig. 2b Fig. 2c

Fig. 2   (a) Comparison of the energy generation rate under adiabatic expansion. The initial conditions are T = 6 × 109 K, = 1 × 107 g cm-3, and a pure 28Si composition. Symbols and colors have the same meaning as in Fig. 1a. Silicon burning under these conditions in endothermic before 10-3 s and exothermic afterward. The transition point is marked with the black vertical line. The two 489 isotope reaction networks, with and without weak reactions, give nearly identical results for these conditions. (b) Ratio of energy generation rate given by the 489 isotope reaction network, with the Fuller et al. (1985) weak reaction rates turned on, to the -chain (green curve) and iso7 (purple curve) reaction networks. The thermodynamic conditions are the same as in Fig. 2a. The error between either of these two small networks and the full networks is larger than the error between the iso7 and -chain networks. (c) Evolution of the mass fractions under adiabatic expansion. The initial conditions are T = 6 × 109 K, = 1 × 107 g cm-3, and a pure 28Si composition. Symbols and colors have the same meaning as in Fig. 1c. Solid curves are for the iso7 reaction network, and open circles for the 13 isotope -chain network. For the -chain network, the sum of the 28Si, 32S, 36Ar, 40Ca, and 44Ti mass fractions are plotted as silicon, and the sum of the 48Cr, 52Fe, and 56Ni mass fractions are plotted as nickel.

     As in the hydrostatic burning case (same initial conditions), the energy generation rate is initially negative because of the strength of the (,) photodisintegration reaction rates. As the temperature falls, however, the strength of the (,) and (,p)(p,) flows begins to dominate. As a result, the energy generation rate changes sign from negative to positive at about 10^-3 s in Figures 2a and 2b. This zero-crossing of the energy generation rate is marked by the vertical black line in Figures 2a and 2b.

     Figure 2a shows that the iso7 and -chain reaction networks give energy generation rates that track with other, at least on global timescales. Each of these reaction networks give the zero-crossing of the energy generation rate at nearly the same point in time. There remains an equal number of electrons per baryon, Y_e = 0.5, throughout the evolution in Figure 2a since not enough time elapses for any weak reactions to significantly change Y_e. This is why the curves for the 489 isotope reaction network with weak reactions turned off (gold triangles) and with weak reactions turned on (purple diamonds) overlay each other. This is also why the iso7 and the 13 isotope -chain reaction networks give energy generation rates that are quite similar to the energy generation rates returned by the 489 isotope reaction networks, except for times larger than about 10^-1 s. At about 10^-1 s the energy generation rate given by the 489 isotope reaction networks is larger than either of the -chain reaction networks because the temperature has fallen enough to impede capture of additional -particles, but not fallen enough to impede flows through various neutron and proton channels, which are not included in the smaller reaction networks.

     For the adiabatic thermodynamic evolutions being considered, Figure 2a suggests that it is generally safe, before 10^-1 s, to calculate the energy generation rate with a small -chain reaction network. It is not generally safe to use an -chain reaction network after 10^-1 s, as the energy generation rate will be in error by an order of magnitude or more.

     The linear scales of Figure 2b indicate that the energy generation rate given by iso7 and -chain reaction networks agree to within 30% over the entire evolution, except at the zero-crossing. This 30% error between the iso7 and -chain networks is only acceptable because the difference between either of the two small networks and the full networks is bigger.

     At times larger than 10^-4 s in Figure 2c, there are large discrepancies between the silicon group abundances as given by the iso7 reaction network (solid curves) and the -chain network (open circles). At 10^-1 s the silicon group abundance of the iso7 reaction network goes below the figure's scale, while the silicon group abundance of the -chain reaction network stays within the figure's scale. These large abundance differences do not cause large differences in the nuclear energy generation rates since differences in small abundances have little effect on the energy release.

     We now consider a case of silicon burning at a larger density where the electron-capture reactions, and flows through the neutron and proton channels, play a more significant role. An initially pure silicon composition at a temperature of 5 × 10⁹ K and density of 1 × 10⁹ g cm^-3 was evolved under hydrostatic conditions for 100 s. The nuclear energy generation rate is shown as a function of time in Figure 3a, the ratio of the -chain to iso7 energy generation rate as a function of time is shown in Figure 3b, and the evolution of the mass fractions are shown in Figure 3c. The format, symbols, and color codes of these plots are the same as those used for the corresponding plots of Figure 1.

Fig. 3a Fig. 3b Fig. 3c

Fig. 3   (a) Comparison of the energy generation rate under hydrostatic conditions of T = 5 × 109 K, = 1 × 109 g cm-3, and an initially pure 28Si composition. Symbols and colors have the same meaning as in Fig. 1a. Silicon burning under these conditions in exothermic. (b) Ratio of energy generation rate given by the 489 isotope reaction network, with the Fuller et al. (1985) weak reaction rates turned on, to the -chain (green curve) and iso7 (purple curve) reaction networks. The thermodynamic conditions are the same as in Fig. 3a. The error between either of these two small networks and the full networks is larger than the error between the iso7 and -chain networks. (c) Evolution of the mass fractions under the hydrostatic conditions of T = 5 × 109 K, = 1 × 109 g cm-3, and an initially pure 28Si composition. Symbols and colors have the same meaning as in Fig. 1c. Solid curves are for the iso7 reaction network, and open circles for the 13 isotope -chain network. For the -chain network, the sum of the 28Si, 32S, 36Ar, 40Ca, and 44Ti mass fractions are plotted as silicon, and the sum of the 48Cr, 52Fe, and 56Ni mass fractions are plotted as nickel.

     The energy generation rate in Figure 3a is initially negative because of the strength of the photodisintegration flows. As the abundance of helium and other -chain nuclei increase due to the photodisintegration of silicon, the strength of the flows which build heavier nuclei increases since these flows are linearly proportional to the helium and -chain nuclei abundances. At about 10^-6 s, these flows toward heavier nuclei dominate, and the energy generation rate changes sign from negative to positive. This zero-crossing of the energy generation rate is marked by the vertical black line in Figures 3a and 3b.

     Figure 3a shows the iso7 and -chain reaction networks do not agree with each other very well, even on scales which span many orders of magnitude. More importantly, neither of these small reaction networks return an energy generation rate that is compatible with the energy generation rate given by the 489 isotope reaction networks. Before 10^-2 s the number of neutrons and protons in the composition are equal, Y_e = 0.5. After 10^-2 s the number of neutrons exceeds the number of protons Y_e < 0.5. This is why the curves for the 489 isotope reaction network with weak reactions turned off (gold triangles) and with weak reactions turned on (purple diamonds) overlay each other before 10^-2 s and follow different tracks afterwards. At the end of the evolution the composition of the 489 isotope reaction network with weak reactions turned on has experienced enough electron capture reactions to become very neutron rich, Y_e = 0.44. The flows through the neutron and proton channels are unavailable to the smaller reaction networks since they have a built-in restriction of Y_e = 0.5.

     The linear scales of Figure 3b magnify the extent of the disagreement. Before about 10^-4 s the energy generation rate given by iso7 and -chain reaction networks can differ by factors of 20. Between the time points where 10% and 90% of the total energy is released, the two small reaction networks give energy generation rates that are within a factor of 3 of each other. More importantly is that neither of the two small reaction networks agree with the more accurate answer given by the full reaction network.

     Figure 3c shows there are large discrepancies not just between the abundances levels, but in their rates of change. This is particularly the case for magnesium, oxygen and nickel. Since the nuclear energy generation rate depends linearly on changes in abundance levels, it is not surprising that the energy generation rates of the two small reaction networks in Figures 3a and 3b disagree substantially.

     An initially pure silicon composition that starts with a temperature of 5 × 10⁹ K and a density of 1 × 10⁹ g cm^-3 was evolved under the adiabatic expansion approximation of equation (9), and the results are shown in Figures 4a, 4b, and 4c. The format, symbols, and color codes of these plots are the same as those used for the corresponding plots of Figure 1.

Fig. 4a Fig. 4b Fig. 4c

Fig. 4   (a) Comparison of the energy generation rate under adiabatic expansion. The initial conditions of T = 5 × 109 K, = 1 × 109 g cm-3, and a pure 28Si composition. Symbols and colors have the same meaning as in Fig. 1a. Silicon burning under these conditions in exothermic before 3 × 10-6 s and exothermic afterward. The transition point is marked with the black vertical line. (b) Ratio of energy generation rate given by the 489 isotope reaction network, with the Fuller et al. (1985) weak reaction rates turned on, to the -chain (green curve) and iso7 (purple curve) reaction networks. The thermodynamic conditions are the same as in Fig. 4a. The error between either of these two small networks and the full networks is larger than the error between the iso7 and -chain networks. (c) Evolution of the mass fractions under adiabatic expansion. The initial conditions of T = 5 × 109 K, = 1 × 109 g cm-3, and a pure 28Si composition. Symbols and colors have the same meaning as in Fig. 1c. Solid curves are for the iso7 reaction network, and open circles for the 13 isotope -chain network. For the -chain network, the sum of the 28Si, 32S, 36Ar, 40Ca, and 44Ti mass fractions are plotted as silicon, and the sum of the 48Cr, 52Fe, and 56Ni mass fractions are plotted as nickel.

     As in the hydrostatic burning case (same initial conditions), the nuclear energy generation rates of the iso7 and -chain reaction networks do not agree. More significantly, neither of these two reaction networks agree with the more accurate answer given by the 489 isotope reaction network. Most of the analysis given for Figure 3 applies to the adiabatic expansion case of Figure 4.

     Figures 3 and 4 indicate that neither the iso7 reaction network, or any -chain network, can accurately track the complexities of silicon burning at densities near 10⁹ g cm^-3. The significant flows through the neutron and proton channels at these thermodynamic conditions cannot be followed by reaction networks with a built-in restriction of Y_e = 0.5. It is not safe to use any -chain reaction network under these conditions, as the energy generation rates will be in error by an order of magnitude or more.

3.2. Hydrostatic and Explosive Carbon-Oxygen Burning

     The flows between ²⁸Si and ⁵⁶Ni in the iso7 reaction network should only be activated when conditions appropriate for silicon burning exist, otherwise these flows should be inactive. Various criteria for the appropriate conditions were experimented with. These included requiring that a fraction of the material exist in heavy (A_i 24) isotopes, that the flows from ²⁴Mg to ²⁸Si were large, and that the time elapsed since any burning began exceeded the time scale for quasi-equilibrium to exist between the silicon and nickel group. The latter conditions adds the memory burden of having to store the time when a zone in a multidimensional model begins to burn, and is particularly unwieldy with adaptive mesh refinement techniques. None of these potential conditions performed as well as turning on silicon burning when the temperature exceeded 2.5 × 10⁹ K and the sum of the carbon and oxygen mass fractions were less than 0.1.

     At the threshold temperature of 2.5 × 10⁹ K, the capture reaction rates are roughly equal to the photodisintegration reaction rates. Slightly smaller or larger threshold temperatures do not significantly change the energy generation rate or nucleosynthesis given by the iso7 reaction network. The precise value for the threshold carbon plus oxygen mass fraction also is not critical, as it serves to shut down premature population of the nickel group abundance. These conditions for turning on silicon burning in the iso7 reaction network are approximations. How well these approximations perform are tested by analyzing the results of hydrostatic and explosive helium, carbon, and oxygen burning.

     An initial composition of half ¹²Chalf ¹⁶O at a temperature of 3 × 10⁹ K and density of 1 × 10⁹ g cm^-3 was evolved under hydrostatic conditions for 10⁹ s. The nuclear energy generation rate is shown as a function of time in Figure 5a, the ratio of the energy generation rate from the -chain and iso7 reaction networks as a function of time is shown in Figure 5b, and the evolution of the mass fractions are shown in Figure 5c. The format, symbols, and color codes of these plots are the same as those used for the corresponding plots of Figure 1.

Fig. 5a Fig. 5b Fig. 5c

Fig. 5   (a) Comparison of the energy generation rate under hydrostatic conditions of T = 3 × 109 K, = 1 × 109 g cm-3, and an initially composition of half 12C and half 16O. Symbols and colors have the same meaning as in Fig. 1a. Carbon-oxygen burning under these conditions is exothermic. (b) Ratio of energy generation rate given by the 489 isotope reaction network, with the Fuller et al. (1985) weak reaction rates turned on, to the -chain (green curve) and iso7 (purple curve) reaction networks. The thermodynamic conditions are the same as in Fig. 5a. The error between either of these two small networks and the full networks is larger than the error between the iso7 and -chain networks. (c) Evolution of the mass fractions under the hydrostatic conditions of T = 3 × 109 K, = 1 × 109 g cm-3, and an initially composition of half 12C and half 16O. Symbols and colors have the same meaning as in Fig. 1c. Solid curves are for the iso7 reaction network, and open circles for the 13 isotope -chain network. For the -chain network, the sum of the 28Si, 32S, 36Ar, 40Ca, and 44Ti mass fractions are plotted as silicon, and the sum of the 48Cr, 52Fe, and 56Ni mass fractions are plotted as nickel.

     Figure 5a shows the iso7 and -chain reaction networks return energy generation rates that track each other, at least on scales that span many orders of magnitude. Before 1 s, both these small reaction networks also give an energy generation rate that is compatible with the energy generation rate given by the 489 isotope reaction networks. In the case of the 489 isotope reaction network with the weak reactions turned off, the flows through various neutron and proton channels are substantial enough at 10² s to give an energy generation rate that is about a factor of 8 larger than the energy generation rate returned by the smaller reaction networks. In the case of the 489 isotope reaction network with the weak reactions turned on, electron capture reactions help drive Y_e below 0.5. Beginning at 1 s, the flows through the larger number of neutron channels are enough to give an energy generation rate that is orders of magnitude larger than the energy generation rate returned by the smaller reaction networks. Again, the built-in restriction of Y_e = 0.5 in -chain reaction networks does not permit an accurate estimate of the energy generation rate if there are substantial flows through neutron and proton channels, or if Y_e 0.5.

     The linear scales of Figure 5b show that the largest source of error in the energy generation rate is made by choosing one of the small -chain networks instead of a large reaction network. While the energy generation rate given by the iso7 and -chain reaction networks agree to within 40% over the entire evolution, the error associated with not using a large reaction network is larger. Neutrino losses (Itoh et al. 1996) are larger than the nuclear energy generation rate after about 10² s, and are indicated as such by the black horizontal line in Figure 5b.

     Figure 5c indicates that the mass fractions as determined by the iso7 reaction network (solid curves) generally track the mass fractions determined by the -chain reaction network (open circles). There are differences in the magnesium abundance between 10^-3 s and 10³ s, and the onset of nickel group abundance around 10³ s. It is the different shapes of the magnesium and nickel group abundance curves that give rise to the differences in the nuclear energy generation rates. Changing the threshold carbon plus oxygen mass fraction from its suggested value of 0.1 to a smaller value (say 0.01) makes the nickel group abundance calculated by the iso7 reaction network rise faster, improving the fit in Figure 5c. An undesirable side-effect of making the suggested composition threshold smaller, however, is to delay the appearance of nickel during silicon burning in Figures 1c4c. Changing the threshold carbon plus oxygen mass fraction to a larger value (say 0.5) has no effect on the silicon burning of Figures 1c-4c, but is produces too much nickel too early during carbon burning in Figure 5c. Our suggested threshold value of 0.1 represents a compromise, and appears to have minimal effects on the silicon or carbon-oxygen burning cases investigated.

     An initial composition of half ¹²Chalf ¹⁶O that starts with a temperature of 3 × 10⁹ K and a density of 1 × 10⁹ g cm^-3 was evolved under the adiabatic expansion approximation of equation (9), and the results are shown in Figures 6a, 6b, and 6c. The format, symbols, and color codes of these plots are the same as those used for the corresponding plots of Figure 1.

Fig. 6a Fig. 6b Fig. 6c

Fig. 6   (a) Comparison of the energy generation rate under adiabatic expansion. The initial conditions are T = 3 × 109 K, = 1 × 109 g cm-3, and a half 12Chalf 16O composition. Symbols and colors have the same meaning as in Fig. 1a. Carbon-oxygen burning under these conditions is exothermic. The two 489 isotope reaction networks, with and without weak reactions, give nearly identical results for these conditions. (b) Ratio of energy generation rate given by the 489 isotope reaction network, with the Fuller et al. (1985) weak reaction rates turned on, to the -chain (green curve) and iso7 (purple curve) reaction networks. The thermodynamic conditions are the same as in Fig. 6a. The error between either of these two small networks and the full networks is larger than the error between the iso7 and -chain networks. (c) Evolution of the mass fractions under adiabatic expansion. The initial conditions of T = 3 × 109 K, = 1 × 109 g cm-3, and a half 12Chalf 16O composition. Symbols and colors have the same meaning as in Fig. 1c. Solid curves are for the iso7 reaction network, and open circles for the 13 isotope -chain network. For the -chain network, the sum of the 28Si, 32S, 36Ar, 40Ca, and 44Ti mass fractions are plotted as silicon, and the sum of the 48Cr, 52Fe, and 56Ni mass fractions are plotted as nickel.

     The nuclear energy generation rate given by the iso7 and -chain reaction networks agree fairly well with each other on global scales of Figure 6a. Both reaction networks also give an energy generation rate that tracks the energy generation rate given by the 489 isotope reaction networks. The differences in the energy generation rate between the large and small reaction networks that arose in the hydrostatic burning case (Fig. 5a) do not occur here because the temperature falls below thermonuclear reaction thresholds at about 10^-1 s. The expansion also slows the electron capture reactions because the density is decreasing. In fact, by about 10^-2 s decays from radioactive nuclei begin to dominate the energy generation rate in the 489 isotope reaction networks.

     The linear scales of Figure 6b indicate that the energy generation rate given by iso7 and -chain reaction networks agree to within 40% over the entire evolution, particularly between the time points where 10% and 90% of the total energy is released. This 40% error between the iso7 and -chain networks is only acceptable because the difference between either of the two small networks and the full networks is bigger. Neutrino losses are larger than the nuclear energy generation rate after about 10^-2 s, and are indicated as such by the black horizontal line in Figure 6b.

     Over the entire evolution in Figure 6c, the iso7 reaction network underproduces the silicon group abundance as compared to the -chain reaction network silicon group abundance, but maintains the same curvatures. Note that the nickel abundance of both the iso7 and -chain reaction networks are below the scale of the plot. Most of the analysis for changes in the silicon burning conditions discussed above for hydrostatic carbon-oxygen burning apply to this case as well.

3.3. Hydrostatic and Explosive Helium Burning

     An initial composition of pure helium at a temperature of 3 × 10⁹ K and density of 1 × 10⁸ g cm^-3 was evolved under hydrostatic conditions for 1 s. The nuclear energy generation rate is shown as a function of time in Figure 7a, the ratio of the -chain and iso7 energy generation rates are shown in Figure 7b, and the evolution of the mass fractions are shown in Figure 7c. The format, symbols, and color codes of these plots are the same as those used for the corresponding plots of Figure 1.

Fig. 7a Fig. 7b Fig. 7c

Fig. 7   (a) Comparison of the energy generation rate under the hydrostatic conditions of T = 3 × 109 K, = 1 × 108 g cm-3, and an initially pure 4He composition. Symbols and colors have the same meaning as in Fig. 1a. Helium burning under these conditions is exothermic. The two 489 isotope reaction networks, with and without weak reactions, give nearly identical results for these conditions. (b) Ratio of energy generation rate given by the 489 isotope reaction network, with the Fuller et al. (1985) weak reaction rates turned on, to the -chain (green curve) and iso7 (purple curve) reaction networks. The thermodynamic conditions are the same as in Fig. 7a. The error between either of these two small networks and the full networks is larger than the error between the iso7 and -chain networks. (c) Evolution of the mass fractions under the hydrostatic conditions of T = 3 × 109 K, = 1 × 108 g cm-3, and an initially pure 4He composition. Symbols and colors have the same meaning as in Fig. 1c. Solid curves are for the iso7 reaction network, and open circles for the 13 isotope -chain network. For the -chain network, the sum of the 28Si, 32S, 36Ar, 40Ca, and 44Ti mass fractions are plotted as silicon, and the sum of the 48Cr, 52Fe, and 56Ni mass fractions are plotted as nickel.

     The energy generation rate returned by the iso7 and -chain reaction networks follow each other, at least on the scales of Figure 7a. Both these small networks also give an energy generation rate that are within 40% of the energy generation rate returned by the 489 isotope reaction network. Over most of the evolution in Figure 7a there remains an equal number of protons and neutrons, Y_e = 0.5. Hence, the curves for the 489 isotope reaction network with weak reactions turned off (gold triangles) and with weak reactions turned on (purple diamonds) overlay each other for most of the evolution. Only for times greater than 1 s, as the compositions settles into its final equilibrium distribution, do weak reactions drive a slight imbalance in the number of neutrons, Y_e = 0.494. These additional flow channels account for the larger energy generation rate returned by the 489 isotope reaction network and with weak reactions turned on.

     The linear scales of Figure 7b show that the energy generation rate given by iso7 and -chain reaction networks agree to within 40% over the entire evolution. This error is acceptable because the difference between either of the two small networks and the full networks is larger.

     Over the entire evolution in Figure 7c, the iso7 reaction network produces silicon group abundances below the scale of the plot. Due to the large helium abundances, silicon group material is transferred immediately into the nickel group. This accounts for the iso7 reaction network giving too large of an energy generation rate before 10^-3 s in Figure 7b, and too small an energy generation rate at later times.

     An initial composition of pure helium that starts with a temperature of 3 × 10⁹ K and a density of 1 × 10⁸ g cm^-3 was evolved under the adiabatic expansion approximation of equation (9), and the results are shown in Figures 8a, 8b, and 8c. The format, symbols, and color codes of these plots are the same as those used for the corresponding plots of Figure 1.

Fig. 8a Fig. 8b Fig. 8c

Fig. 8   (a) Comparison of the energy generation rate under adiabatic expansion. The initial conditions of T = 3 × 109 K, = 1 × 108 g cm-3, and a pure 4He composition. Symbols and colors have the same meaning as in Fig. 1a. Helium burning under these conditions is exothermic. The two 489 isotope reaction networks, with and without weak reactions, give nearly identical results for these conditions. (b) Ratio of energy generation rate given by the 489 isotope reaction network, with the Fuller et al. (1985) weak reaction rates turned on, to the -chain (green curve) and iso7 (purple curve) reaction networks. The thermodynamic conditions are the same as in Fig. 8a. The error between either of these two small networks and the full networks is larger than the error between the iso7 and -chain networks. (c) Evolution of the mass fractions under adiabatic expansion. The initial conditions of T = 3 × 109 K, = 1 × 108 g cm-3, and a pure 4He composition. Symbols and colors have the same meaning as in Fig. 1c. Solid curves are for the iso7 reaction network, and open circles for the 13 isotope -chain network. For the -chain network, the sum of the 28Si, 32S, 36Ar, 40Ca, and 44Ti mass fractions are plotted as silicon, and the sum of the 48Cr, 52Fe, and 56Ni mass fractions are plotted as nickel.

     On the global scales of Figure 8a, the energy generation rate returned by the iso7 and -chain reaction networks are similar, and both of these small reaction networks give an energy generation rate that follows the energy generation rate returned by the 489 isotope reaction networks. The differences in the energy generation rate between the large and small reaction networks that arose in the hydrostatic burning case (Fig. 7a) do not occur in the adiabatic expansion case because the temperature falls below reaction thresholds by about 10^-2 s.

     The linear scales of Figure 8b show that the largest source of error is made by choosing one of the small -chain networks instead of a large reaction network. While the energy generation rate given by the iso7 and -chain reaction networks agree to within 40% over the entire evolution, the error associated with not using a large reaction network is larger.

     In Figure 8c silicon group material is transferred immediately into the nickel group, as in the hydrostatic burning case of Figure 7c. This paucity of silicon group material and the overabundance of nickel group material accounts for the iso7 reaction network giving too large of an energy generation rate before 10^-3 s, and too small an energy generation rate at later times.

4. SUMMARY

     The iso7 reaction network reproduces the nuclear energy generation rate of standard -chain reaction network to within 30% during hydrostatic and explosive helium, carbon, oxygen, and silicon burning whenever Y_e = 0.5 (Figs. 58). If Y_e is outside the rough limits 0.5 ± 0.01, or if there are significant flows through neutron and proton channels, then it is not safe to use any -chain reaction network for the energy generation rate. Under these conditions, the energy generation rate returned by an -chain reaction network can be in error by an order of magnitude or more.

     The iso7 reaction network offers reasonable approximations for the energy generation rate of helium, carbon, neon, and early oxygen burning under normal stellar evolution. For more advanced burning stages the approximations inherent in any -chain reaction network become a poor approximation, the iso7 network included.

     Replacing a 13 isotope -chain reaction network with the 7 isotope iso7 reaction network in a multidimensional hydrodynamics program provides at least three advantages. First, the amount of memory required to store the abundances at each grid point is cut in about half since about half as many isotopes are evolved. Since the largest block of memory in a multidimensional hydrodynamic program is usually reserved for storing the isotopic abundances at every grid point, this memory savings can be quite significant for three-dimensional models on shared memory parallel computers. The memory savings can be reinvested toward covering a larger physical domain in the simulation, or increasing the spatial resolution of the model. Second, the amount of CPU time required to solve the linear system of equations that arise from implicit (or semi-implicit) integration methods is reduced by a factor of 4, assuming standard N² dense matrix solution techniques are used. Third, the amount of CPU time required to evaluate the nuclear reaction rates is reduced by about 75% since the iso7 reaction network evaluates 15 reaction rates while a standard 13 isotope -chain network evaluates about 59 reaction rates (24 rates for the (,) their inverse reactions + 3 heavy ion reaction rates + 32 rates for the (,p)(p,) sequences and their inverses). The total CPU cost of advancing a small nuclear reaction network (less than 30 isotopes) over a time step tends to be dominated by evaluation of the reaction rates and not the linear algebra (Timmes 1999). Thus, reducing the number of reaction rates evaluated is usually more important to improving the overall performance of evolving a reaction network than reducing the linear algebra costs. Another approach to reduce the CPU costs of evaluating the analytic reaction rates is to evaluate the reaction rates by interpolating a one-dimensional table lookup scheme (temperature as the independent variable, the reaction rate as the dependent variable).

     Hix et al. (1998) and Hix & Thielemann (1999) discuss a nuclear reaction network which also has seven isotopes, which they call 7. The main difference between the iso7 and 7 reaction networks is how silicon burning is handled. The iso7 reaction network makes direct transitions between ²⁸Si and ⁵⁶Ni (see § 2), while the 7 reaction network uses the conditions of quasi-equilibrium to solve for the silicon group nuclei (²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti) and the iron group nuclei (⁴⁸Cr, ⁵²Fe, and ⁵⁶Ni). That is, the 7 reaction network must perform a root find for the free -particle abundance at each integration step during silicon burning. In exchange for this small additional computational burden, the 7 reaction network appears to reproduce the abundance levels of a standard -chain network better than the iso7 reaction network.

     An elevated level of attention should be given if the iso7 reaction network, or any other nuclear reaction network, is to be used in Type II core collapse models which employ operator split hydrodynamics. A numerical instability can occur as the nuclei move into their nuclear statistical equilibrium distribution with this common hydrodynamic method. As the temperature rises the photodisintegration reaction rates dominate the nuclear flows. This gives rise to a negative nuclear energy generation rate, which cools the material. Electron capture reaction rates dominate the nuclear flows in the cooled material, which gives a positive energy generation rate, which raises the temperature. This cycle can only be broken by keeping the time steps very short so that the composition does not violently oscillate. A better, if more involved, way to avoid this numerical instability is to solve the burning and hydrodynamics in a fully coupled manner (e.g., Müller 1986).

     The authors thank the referee, William Raphael Hix, for a pragmatic and exhausting review that improved the manuscript. This work was supported by the Department of Energy under grant B341495 to the Center on Astrophysical Thermonuclear Flashes at the University of Chicago, by NASA (NAG5-8128 and MIT SC A292701), NSF (AST 97-31569 and INT-97-31569), and the Department of Energy ASCI Program (W-7405-ENG-48) at the University of California as Santa Cruz, and under the auspices of the US Department of Energy by the University of California Lawrence Livermore National Laboratory under contract W-7405-Eng-48.

APPENDIX A

Iso7 REACTION NETWORK

     Three FORTRAN subroutines that implement the iso7 reaction network are given below. [Formatted source code file also available.] The first subroutine, named iso7, implements the right-hand sides of the seven stiff ordinary differential equations in equation (3). The second subroutine, named iso7jac, returns the 7 × 7 Jacobian matrix in dense matrix format. The third subroutine, named iso7rat, returns the 15 reaction rates required by routines iso7 and iso7jac. In subroutine iso7rat the rates for the binary reactions ¹²C(,)¹⁶O, ¹⁶O(,)²⁰Ne, ²⁰Ne(,)²⁴Mg, ²⁴Mg(,)²⁸Si plus their inverses, the three-body reaction ⁴He(2)¹²C plus its inverse, and the heavy ion reactions ¹²C + ¹²C, ¹²C + ¹⁶O, ¹⁶O + ¹⁶O, are from Caughlan & Fowler (1998), while ⁴⁰Ca(,)⁴⁴Ti plus its inverse is from Woosley et al. (1978). This set of reaction rates was used to generate all the figures in this paper. One may wish to substitute the corresponding reaction rates from the NACRE (Angulo et al. 1999) or Rauscher & Thielemann (1998) compilations, with the realization that the results may then differ from some of the results of this paper. To use the three iso7 reaction network routines one must supply a time integration method and a linear algebra package. Various time integration methods and linear algebra packages when applied to nuclear reaction networks are surveyed in Timmes (1999) and Hix & Thielemann (1999). After a successful integration step, one applies equation (5) to obtain the nuclear energy generation rate.

                           subroutine iso7(y,dydt)

                           implicit none

                           save

     c..

     c..this routine sets up the system of equations for the iso7 network.

     c..

     c..input is y(1:7) the abundance fractions

     c..output is dydt(1:7) the right hand side of the network equations

     c..

     c..isotopes are: he4, c12, o16, ne20, mg24, si28, ni56

     c..

     c..declare

                           integer i

                           double precision y(1),dydt(1),t9,t9i,t932,t9i32,rsi2ni,rni2si

     c..for easy isotope identification

                           integer ihe4,ic12,io16,ine20,img24,isi28,ini56

                           parameter (ihe4=1, ic12=2, io16=3, ine20=4,

                          1                                img24=5, isi28=6, ini56=7)

     c..for easy rate identification

                           integer ircag,iroga,ir3a,irg3a,ir1212,ir1216,ir1616,

                          1                                iroag,irnega,irneag,irmgga,irmgag,irsiga,ircaag,

                          2                                irtiga

                           parameter (ircag=1, iroga=2, ir3a=3, irg3a=4,

                          1                                ir1212=5, ir1216=6, ir1616=7, iroag=8,

                          2                                irnega=9, irneag=10, irmgga=11, irmgag=12,

                          3                                irsiga=13, ircaag=14, irtiga=15)

     c..communicate the temperature, density, and reaction rates

                          integer nrate

                          parameter (nrate=15)

                          double precision rate(nrate),temp,den

                          common /iso7c1/ rate,temp,den

     c..some frequent factors

                          t9 = temp * 1.0d-9

                          t9i = 1.0d0/t9

                          t932 = t9 * sqrt(t9)

                          t9i32 = 1.0d0/t932

     c..rsi2ni is the rate for silicon to nickel

     c..rni2si is the rate for nickel to silicon

                          rsi2ni = 0.0d0

                          rni2si = 0.0d0

                          if (t9 .gt. 2.5 .and. y(ic12)+y(io16) .gt. 0.004) then

                               rsi2ni = (t9i32^**3) * exp(239.42*t9-i;74.741) *

                      1                                den^**3 * y(ihe4)^**3 * rate(ircaag)*y(isi28)

                               rni2si = min(1.0d20,(t932^**3) * exp(-274.12*t9i+74.914)

                      1                                                                            * rate(irtiga)/(den^**3 * y(ihe4)^**3))

                          end if

     c..4he reactions

                          dydt(ihe4) = 3.0 * y(ic12) * rate(irg3a)

                      1                                        - 3.0 * y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a)

                      2                                        + y(io16) * rate(iroga)

                      3                                        - y(ic12) * y(ihe4) * rate(ircag)

                      4                                        + y(ic12) * y(ic12) * rate(ir1212)

                      5                                        + 0.5d0 * y(ic12) * y(io16) * rate(ir1216)

                      6                                        + y(io16) * y(io16) * rate(ir1616)

                      7                                        - y(io16) * y(ihe4) * rate(iroag)

                      8                                        + y(ine20) * rate(irnega)

                      9                                        + y(img24) * rate(irmgga)

                      &                                        - y(ine20) * y(ihe4) * rate(irneag)

                      1                                        + y(isi28) * rate(irsiga)

                      2                                        - y(img24) * y(ihe4) * rate(irmgag)

                      3                                        - 7.0d0 * rsi2ni * y(ihe4)

                      4                                        + 7.0d0 * rni2si * y(ini56)

     c..12c reactions

                          dydt(ic12) = y(ihe4) * y(ihe4) * y(ihe4) * rate(ir3a)

                      1                                        - y(ic12) * rate(irg3a)

                      2                                        + y(io16) * rate(iroga)

                      3                                        - y(ic12) * y(ihe4) * rate(ircag)

                      4                                        - 2.0d0 * y(ic12) * y(ic12) * rate(ir1212)

                      5                                        - y(ic12) * y(io16) * rate(ir1216)

     c..16o reactions

                          dydt(io16) = - y(io16) * rate(iroga)

                      1                                        + y(ic12) * y(ihe4) * rate(ircag)

                      2                                        - y(ic12) * y(io16) * rate(ir1216)

                      3                                        - 2.0d0 * y(io16) * y(io16) * rate(ir1616)

                      4                                        - y(i016) * y(ihe4) * rate(iroag)

                      5                                        + y(ine20) * rate(irnega)

     c..20ne reactions

                          dydt(ine20) = y(ic12) * y(ic12) * rate(ir1212)

                      1                                              + y(io16) * y(ihe4) * rate(iroag)

                      2                                              - y(ine20) * rate(irnega)

                      3                                              + y(img24) * rate(irmgga)

                      4                                              - y(ine20) * y(ihe4) * rate(irneag)

     c..24mg reactions

                          dydt(img24) = 0.5d0 * y(ic12) * y(io16) * rate(ir1216)

                      1                                              - y(img24) * rate(irmgga)

                      2                                              + y(ine20) * y(ihe4) * rate(irneag)

                      3                                              + y(isi28) * rate(irsiga)

                      4                                              - y(img24) * y(ihe4) * rate(irmgag)

     c..28si reactions

                          dydt(isi28) = 0.5d0 * y(ic12) * y(io16) * rate(ir1216)

                      1                                              + y(io16) * y(io16) * rate(ir1616)

                      2                                              - y(isi28) * rate(irsiga)

                      3                                              + y(img24) * y(ihe4) * rate(irmgag)

                      4                                              - rsi2ni * y(ihe4)

                      5                                              + rni2si * y(ini56)

     c..ni56 reactions

                          dydt(ini56) = rsi2ni * y(ihe4)

                      1                                              -rni2si*y(ini56)

                          return

                          end

                          subroutine iso7jac(y,dfdy)

                          implicit none

                          save

     c..

     c..this routine sets up the jacobian of the iso7 network in dense matrix form.

     c..

     c..input is y(1:7) the abundance fractions

     c..output is dfdy(1:7,1:7) the jacobian of the network equations

     c..

     c..declare

                           integer i,j

                           double precision y(1),dfdy(7,7),t9,t9i,t932,t9i32,

                          1                                               rsi2ni,rni2si,rsi2nida,rsi2nidsi,rni2sida

     c..for easy isotope identification

                           integer ihe4,ic12,io16,ine20,img24,isi28,ini56

                           parameter (ihe4=1, ic12=2, io16=3, ine20=4,

                          1                                img24=5, isi28=6, ini56=7)

     c..for easy rate identification

                           integer ircag,iroga,ir3a,irg3a,ir1212,ir1216,ir1616,

                          1                                iroag,irnega,irneag,irmgga,irmgag,irsiga,ircaag,

                          2                                irtiga

                           parameter (ircag=1, iroga=2, ir3a=3, irg3a=4,

                          1                                ir1212=5, ir1216=6, ir1616=7, iroag=8,

                          2                                irnega=9, irneag=10, irmgga=11, irmgag=12,

                          3                                irsiga=13, ircaag=14, irtiga=15)

     c..communicate the temperature, density, and reaction rates

                           integer nrate

                           parameter (nrate=15)

                           double precision rate(nrate),temp,den

                           common /iso7c1/ rate,temp,den

     c..zero the jacobian

                           do j=1,7

                               do i=1,7

                                   dfdy(i,j) = 0.0d0

                               enddo

                           enddo

     c..some frequent factors

                           t9 = temp * 1.0d-9

                           t9i = 1.0d0/t9

                           t932 = t9 * sqrt(t9)

                           t9i32 = 1.0d0/t932

     c..rsi2ni is the rate for silicon to nickel

     c..rni2si is the rate for nickel to silicon

                           rsi2ni = 0.0d0

                           rsi2nida = 0.0d0

                           rsi2nidsi = 0.0d0

                           rni2si = 0.0d0

                           rni2sida = 0.0d0

                           if (t9 .gt. 2.5 .and. y(ic12)+y(io16) .gt. 0.004) then

                               rsi2ni = (t9i32^**3) * exp(239.42*t9i-74.741) *

                          1                                    den^**3 * y(ihe4)^**3 * rate(ircaag)*y(isi28)

                               rsi2nida = 3.0d0 * rsi2ni/y(ihe4)

                               rsi2nidsi = rsi2ni/y(isi28)

                               rni2si = min(1.0d20,(t932^**3) * exp(-274.12*t9i+74.914)

                          1                                                                        * rate(irtiga)/(den^**3*y(ihe4)^**3))

                               rni2sida = -3.0d0 * rni2si/y(ihe4)

                               if (rni2si .eq. 1.0d20) rni2sida = 0.0d0

                           end if

     c..4he jacobian elements

                           dfdy(ihe4,ihe4) = -6.0d0 * y(ihe4) * y(ihe4) * rate(ir3a)

                          1                                                              -y(ic12) * rate(ircag) - y(io16) * rate(iroag)

                          2                                                              -y(ine20)*rate(irneag) - y(img24)*rate(irmgag)

                          3                                                              -7.0d0 * rsi2ni - 7.0d0 * rsi2nida * y(ihe4)

                          4                                                              + 7.0d0 * rni2sida * y(ini56)

                           dfdy(ihe4,ic12) = 3.0d0 * rate(irg3a) - y(ihe4) * rate(ircag)

                          1                                                              + 2.0d0 * y(ic12) * rate(ir1212)

                          2                                                              + 0.5d0 * y(io16) * rate(ir1216)

                           dfdy(ihe4,io16) = rate(iroga) + 0.5d0 * y(ic12) * rate(ir1216)

                          1                                                              + 2.0d0 * y(io16) * rate(ir1616)

                          2                                                              - y(ihe4) * rate(iroag)

                           dfdy(ihe4,ine20) = rate(irnega) - y(ihe4) * rate(irneag)

                           dfdy(ihe4,img24) = rate(irmgga) - y(ihe4) * rate(irmgag)

                           dfdy(ihe4,isi28) = rate(irsiga) - 7.0d0 * rsi2nidsi * y(ihe4)

                           dfdy(ihe4,ini56) = 7.0d0 * rni2si

     c..12c jacobian elements

                           dfdy(ic12,ihe4) = 3.0d0 * y(ihe4) * y(ihe4) * rate(ir3a)

                          1                                                              - y(ic12) * rate(ircag)

                           dfdy(ic12,ic12) = -rate(irg3a) - y(ihe4) * rate(ircag)

                          1                                                              - 4.0d0 * y(ic12) * rate(ir1212)

                          2                                                              - y(io16) * rate(ir1216)

                           dfdy(ic12,io16) = rate(iroga) - y(ic12) * rate(ir1216)

     c..16o jacobian elements

                           dfdy(io16,ihe4) = y(ic12) * rate(ircag) - y(io16) * rate(iroag)

                           dfdy(io16,ic12) = y(ihe4) * rate(ircag) - y(io16) * rate(ir1216)

                           dfdy(io16,io16) = -rate(iroga) - y(ic12) * rate(ir1216)

                          2                                                              - 4.0d0 * y(io16) * rate(ir1616)

                          3                                                              - y(ihe4) * rate(iroag)

                           dfdy(io16,ine20) = rate(irnega)

     c..20ne jacobian elements

                           dfdy(ine20,ihe4) = y(io16)*rate(iroag) - y(ine20)*rate(irneag)

                           dfdy(ine20,ic12) = 2.0d0 * y(ic12) * rate(ir1212)

                           dfdy(ine20,io16) = y(ihe4) * rate(iroag)

                           dfdy(ine20,ine20) = -rate(irnega) - y(ihe4) * rate(irneag)

                           dfdy(ine20,img24) = -rate(irmgga)

     c..24mg jacobian elements

                           dfdy(img24,ihe4) = y(ine20)*rate(irneag) - y(img24)*rate(irmgag)

                           dfdy(img24,ic12) = 0.5d0 * y(io16) * rate(ir1216)

                           dfdy(img24,io16) = 0.5d0 * y(ic12) * rate(ir1216)

                           dfdy(img24,ine20) = y(ihe4) * rate(irneag)

                           dfdy(img24,img24) = -rate(irmgga) - y(ihe4) * rate(irmgag)

                           dfdy(img24,isi28) = rate(irsiga)

     c..28si jacobian elements

                           dfdy(isi28,ihe4) = y(img24) * rate(irmgag) - rsi2ni

                          1                                                              - rsi2nida * y(ihe4) + rni2sida * y(ini56)

                           dfdy(isi28,ic12) = 0.5d0 * y(io16) * rate(ir1216)

                           dfdy(isi28,io16) = 2.0d0 * y(io16) * rate(ir1616)

                          1                                                              + 0.5d0 * y(ic12) * rate(ir1216)

                           dfdy(isi28,img24) = y(ihe4) * rate(irmgag)

                           dfdy(isi28,isi28) = -rate(irsiga) - rsi2nidsi * y(ihe4)

                           dfdy(isi28,ini56) = rni2si

     c..ni56 jacobian elements

                           dfdy(ini56,ihe4) = rsi2ni + rsi2nida * y(ihe4)

                          1                                                              - rni2sida * y(ini56)

                           dfdy(ini56,isi28) = rsi2nidsi * y(ihe4)

                           dfdy(ini56,ini56) = -rni2si

                           return

                           end

                           subroutine iso7rat

                           implicit none

                           save

     c..

     c..this routine generates 15 nuclear reaction rates for the iso7 network:

     c.. r3a = triple-alfa rg3a = c12 to triple alfa

     c.. rcag = c12(ag) roga = o16(ga)

     c.. roag = o16(ag) rnega = ne20(ga)

     c.. rneag = ne20(ag) rmgga = mg24(ga)

     c.. rmgag = mg24(ag) rsiga = si28(ga)

     c.. rcaag = ca40(ag) rtiga = ti44(ga)

     c.. r1212 = c12+c12 r1216 = o16+c12 r1515 = o16+o16

     c..

     c..all the rates below are from caughlan and fowler 1988, except ca40(ag) and

     c..its inverse, which is from woosley et al. 1978. one may wish to substitute

     c..the corresponding reaction rates from the nacre (angulo et al. 1999) or

     c..rauscher \& thielemann (1998) libraries.

     c..

     c..declare

                           integer i

                           double precision tt9,t9r,t9,t912,t913,t923,t943,t953,t932,

                          1                                                   t92,t93,t972,t9r32,t9i,t9i13,t9i23,t9i32,

                          2                                                   t9i12,t9ri,term,term1,term2,term3,rev,

                          3                                                   r2abe,t9a,t9a13,t9a56,t9a23,gt9h,

                          4                                                   rbeac,oneth,fivsix

                           parameter (oneth = 1.0d0/3.0d0, fivsix=5.0d0/6.0d0)

     c..for easy isotope identification

                           integer ihe4,ic12,io16,ine20,img24,isi28,ini56

                           parameter (ihe4=1, ic12=2, io16=3, ine20=4,

                          1                                img24=5, isi28=6, ini56=7)

     c..for easy rate identification

                           integer ircag,iroga,ir3a,irg3a,ir1212,ir1216,ir1616,

                          1                                iroag,irnega,irneag,irmgga,irmgag,irsiga,ircaag,

                          2                                irtiga

                           parameter  (ircag=1, iroga=2, ir3a=3, irg3a=4,

                          1                                ir1212=5, ir1216=6, ir1616=7, iroag=8,

                          2                                irnega=9, irneag=10, irmgga=11, irmgag=12,

                          3                                irsiga=13, ircaag=14, irtiga=15)

     c..communicate the temperature, density, and reaction rates

                           integer nrate

                           parameter (nrate=15)

                           double precision rate(nrate),temp,den

                           common /iso7c1/ rate,temp,den

     c..zero the rates

                           do i=1,nrate

                               rate(i) = 0.0d0

                           enddo

     c..some temperature factors

     c..limit t9 to 10, except for the reverse ratios

                           tt9 = temp * 1.0e-9

                           t9r  = tt9

                           t9 = min(tt9,10.0d0)

                           t912 = sqrt(t9)

                           t913 = t9^**oneth

                           t923 = t913*t913

                           t943 = t9*t913

                           t953 = t9*t923

                           t932  = t9*t912

                           t92 = t9*t9

                           t93 = t9*t92

                           t972 = t912*t93

                           t9r32 = t9r*sqrt(t9r)

                           t9i = 1.0d0/t9

                           t9i13 = 1.0d0/t913

                           t9i23 = 1.0d0/t923

                           t9i32 = 1.0d0/t932

                           t9i12 = 1.0d0/t912

                           t9ri = 1.0d0/t9r

     c..triple alpha to c12

                           r2abe = (7.40d+05 * t9i32) * exp(-1.0663*t9i) +

                          1                          4.164d+09 * t9i23 * exp(-13.49*t9i13 - t92/0.009604) *

                          2                          (1.0d0 + 0.031*t913 + 8.009*t923 + 1.732*t9 +

                          3                          49.883*t943 + 27.426*t953)

                           rbeac = (130.*t9i32) * exp(-3.3364*t9i) +

                          1                          2.510d+07 * t9i23 * exp(-23.57*t9i13 - t92/0.055225)*

                          2                          (1.0d0 + 0.018*t913 + 5.249*t923 + 0.650*t9 +

                          3                          19.176*t943 + 6.034*t953)

                           if (t9.gt.0.08) then

                               term = 2.90d-16 * (r2abe*rbeac) +

                          1                          0.1 * 1.35d-07 * t9i32 * exp(-24.811*t9i)

                           else

                               term = 2.90d-16*(r2abe*rbeac) *

                          1                      (0.01 + 0.2*(1.0d0 + 4.0d0*exp(-(0.025*t9i)^**3.263))/

                          2                      (1.0d0 + 4.0d0*exp(-(t9/0.025)^**9.227))) +

                          3                      0.1 * 1.35d-07 * t9i32 * exp(-24.811*t9i)

                           end if

                           rate(ir3a) = term * (den*den)/6.0d0

                           rev = 2.00d+20*exp(-84.424*t9ri)

                           rate(irg3a) = rev*(t9r*t9r*t9r) * term

     c..c12 + c12 reaction from caughlan and fowler 1988

                           t9a = t9/(1.0d0 + 0.0396*t9)

                           t9a13 = t9a^**oneth

                           t9a56 = t9a^**fivsix

                           term = 4.27d+26 * t9a56/t932 *

                          1                                                       exp(-84.165/t9a13-2.12d-03*t9*t9*t9)

                           rate(ir1212) = 0.5d0 * den * term

     c..c12 + o16 reaction

                           if (t9.ge.0.5) then

                               t9a = t9/(1.+0.055*t9)

                                t9a13 = t9a^**oneth

                                t9a23 = t9a13*t9a13

                                t9a56 = t9a^**fivsix

                                term = 1.72d+31 * t9a56 * t9i32 * exp(-106.594/t9a13) /

                          1                              (exp(-0.18*t9a*t9a) + 1.06d-03*exp(2.562*t9a23))

                                rate(ir1216) = den * term

                           else

                                rate(ir1216) = 0.0d0

                           endif

     c..16o+16o rate

                           term = 7.10d36 * t9i23 *

                          1                                                        exp(-135.93 * t9i13 - 0.629*t923 -

                          2                                                        0.445*t943 + 0.0103*t9*t9)

                           rate(ir1616) = 0.5d0 * den * term

     c..12c(ag)16o and inverse

                           term = 1.04d8/(t92*(1.0d0 + 0.0489d0*t9i23)^**2) *

                          1                                                        exp(-32.120d0*t9i13 - t92/12.222016d0)

                          2                                                   + 1.76d8/(t92*(1.0d0 + 0.2654d0*t9i23)^**2) *

                          3                                                        exp(-32.120d0*t9i13)

                          4                                                   + 1.25d3*t9i32*exp(-27.499*t9i)

                          5                                                   + 1.43d-2*t92*t93*exp(-15.541*t9i)

                           term = 1.7d0 * term

                           rate(ircag) = term * den

                           rev = 5.13d10 * t9r32 * exp(-83.111*t9ri)

                           rate(iroga) = rev * term

     c..16o(ag)20ne + inverse

                           term1 = 9.37d9 * t9i23 *

                          1                                                        exp(-39.757*t9i13 - t92/2.515396)

                           term2 = 62.1 * t9i32 * exp(-10.297*t9i) +

                          1                                                        538.0d0 * t9i32 * exp(-12.226*t9i) +

                          2                                                        * t92 * exp(-20.093*t9i)

                           term = term1 + term2

                           rate(iroag) = den * term

                           rev = 5.65d+10*t9r32*exp(-54.937*t9ri)

                           rate(irnega) = rev * term

     c..20ne(ag)24mg + inverse

                           term1 = 4.11d+11 * t9i23 *

                          1                                                        exp(-46.766*t9i13 - t92/4.923961) *

                          2                                                        (1.0d0 + 0.009*t913 + 0.882*t923 + 0.055*t9 +

                          3                                                        0.749*t943 + 0.119*t953)

                           term2 = 5.27d+03 * t9i32 * exp(-15.869*t9i) +

                          1                                                        6.51d+03 * t912 * exp(-16.223*t9i)

                           term3 = 0.1d0 * (42.1 * t9i32 * exp(-9.115*t9i) +

                          1                                                        32.0 * t9i23 * exp(-9.383*t9i))

                           term = (term1+term2+term3)/

                          1                                                        (1.0d0 + 5.0d0*exp(-18.960*t9i))

                           rate(irneag) = den * term

                           rev = 6.01d+10 * t9r32 * exp(-108.059*t9ri)

                           rate(irmgga) = term * rev

     c..24mg(ag)28si + inverse

                           term1 = (1.0d0 + 5.0d0*exp(-15.882*t9i))

                           term = (4.78d+01 * t9i32 * exp(-13.506*t9i) +

                          1                                                        2.38d+03 * t9i32 * exp(-15.218*t9i) +

                          2                                                        2.47d+02 * t932 * exp(-15.147*t9i) +

                          3                                                        0.1 * (1.72d-09 * t9i32 * exp(-5.028*t9i) +

                          4                                                        1.25d-03 * t9i32 * exp(-7.929*t9i) +

                          5                                                        2.43d+01 * t9i * exp(-11.523*t9i)))/term1

                           rate(irmgag) = den * term

                           rev = 6.27d+10 * t9r32 * exp(-115.862*t9ri)

                           rate(irsiga) = rev * term

     c..40ca(ag)44ti + inverse

                           term = 4.66d+24 * t9i23 * exp(-76.435 * t9i13 *

                          1                                                        (1.0d0 + 1.650d-02*t9 + 5.973d-03*t92 -

                          2                                                             3.889d-04*t93))

                           rev = 6.843d+10 * t9r32 * exp(-59.510*t9ri)

                           rate(ircaag) = den*term

                           rate(irtiga) = rev*term

     c..the user may want to insert screening factors for the relevant rates here

                           return

                           end

REFERENCES

• Angulo, C., et al. 1999, Nucl. Phys., A, 656, 3 First citation in article | CrossRef Link | NASA ADS | Order from Infotrieve
• Bodansky, D., Clayton, D. D., & Fowler, W. A. 1968, ApJS, 16, 299 First citation in article | CrossRef Link | NASA ADS | Order from Infotrieve
• Burgers, J. M. 1969, Flow Equations for Composite Gases (New York: Academic) First citation in article | NASA ADS
• Caughlan, G. R., & Fowler, W. A. 1988, Atomic Data and Nuclear Data Tables, 40, 283 First citation in article | CrossRef Link | NASA ADS | Order from Infotrieve
• Chapman, S., & Cowling, T. G. 1970, The Mathematical Theory of Non-uniform Gases (Cambridge: Cambridge Univ. Press) First citation in article | NASA ADS
• Fowler, R. A., & Hoyle, F. 1964, ApJ, 1, 1 First citation in article | Order from Infotrieve
• Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1 First citation in article | CrossRef Link | NASA ADS | Order from Infotrieve
• Hix, W. R., Khokhlov, A. M., Wheeler, J. G., & Thielemann, F.-K. 1998, ApJ, 503, 332 First citation in article | IOP Article
• Hix, W. R., & Thielemann, F.-K. 1996, ApJ, 406, 869 First citation in article | CrossRef Link | NASA ADS | Order from Infotrieve
• . 1999, JCAM, 109, 321 First citation in article | CrossRef Link | NASA ADS | Inspec Abstract | Order from Infotrieve
• Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411 First citation in article | CrossRef Link | NASA ADS | Order from Infotrieve
• Müller, E. 1986, A&A, 162, 103 First citation in article | NASA ADS | Inspec Abstract | Order from Infotrieve
• Rauscher, T., & Thielemann, F.-K. 1998, in Stellar Evolution, Stellar Explosions, and Galactic Chemical Evolution, ed. A. Mezzacappa (Philadelphia: Inst. Phys.), 134 First citation in article
• Timmes, F. X. 1999, ApJS, 124, 241 First citation in article | IOP Article
• Woosley, S. E. 1986, in Nucleosynthesis and Chemical Evolution, 16th Adv. Course, Swiss Soc. Astrophys. Astron., ed. B. Hauck, A. Maeder, & G. Meynet (Geneva: Geneva Obs.), 1 First citation in article | NASA ADS
• Woosley, S. E., Arnett, W. D., & Clayton, D. D. 1973, ApJS, 26, 231 First citation in article | CrossRef Link | NASA ADS | Order from Infotrieve
• Woosley, S. E., Fowler, W. A., Holmes, J. A., & Zimmerman, B. A. 1978, At. Data Nucl. Data Tables, 22, 71 First citation in article | CrossRef Link | Order from Infotrieve
• Williams, F. A. 1988, Combustion Theory (Menlo Park; Benjamin-Cummings) First citation in article