Nucleosynthesis in thermonuclear supernovae

Abstract

We review our understanding of the nucleosynthesis that occurs in thermonuclear supernovae and their contribution to Galactic Chemical evolution. We discuss the prospects to improve the modeling of the nucleosynthesis within simulations of these events.

Graphical abstract

Keywords

supernovae / nuclear reactions / explosive nucleosynthesis / stellar spectra / hydrodynamics / chemical evolution

1  Introduction

The nuclei heavier than lithium observed in stars and galaxies are thought to be the result of thermonuclear processing in previous generations of stars. The most spectacular example of such processing occurs in the explosion of a star as a supernova. Determining the accurate mass distribution of ejecta and the relative elemental abundances of different supernovae types is of the utmost importance in understanding the role each type plays. Achieving a better understanding of the physics of supernovae, by increasing the quality and predictive power of numerical models and the resulting nucleosynthesis calculations, is fundamental to a better understanding of chemical evolution throughout the Universe. Some nuclear abundances can be used as cosmochronometers (to determine the age of the Galaxy), others power light curves, still others appear as anomalous abundances found in tiny meteoritic grains in our Solar System.

The detailed theory of how thermonuclear supernovae (hereafter SN Ia) evolve and explode is still the subject of considerable effort, motivated in part by the role these supernovae have played as cosmological distance indicators. A number of scenarios have been suggested (see Ref. [1], for a review), but there is a broad consensus that Type Ia supernovae are the outcome of the thermonuclear explosion of a carbon-oxygen white dwarf, hereafter CO-WD (see e.g., Refs. [2, 3]). These supernovae show no signs of hydrogen in their spectra but intermediate mass elements such as Si, S, Ca, and Mg are in evidence near maximum light, and many lines of Fe are present at late times. They also show no He. The CO-WD approaches the Chandrasekhar mass (hereafter M_{ch}) by a yet uncertain mechanism, presumably accretion from a companion star, and is completely disrupted by a thermonuclear explosion (see e.g., Ref. [1]). A strong argument in favor of this hypothesis is the ability of models of these explosions to fit the observed light curves and spectra well. Modeling a Type Ia supernova explosion, therefore, requires that we compute the thermonuclear disruption of a white dwarf. However, the evolution of massive white dwarves to explosion is very uncertain, leaving room for some diversity in the allowed set of initial conditions (such as the temperature profile at ignition), and also the physics of thermonuclear burning in degenerate matter is complex and not fully understood.

The M_{ch} WD has long been been favored as the progenitor for a majority of SNe Ia (e.g., Refs. [4–7]). The mass of the WD could reach M_{ch} by several evolutionary paths, either by a mass transfer from a giant/main sequence binary companion (single degenerate scenario; e.g., Refs. [8, 9]) or as a result of merging with a degenerate WD binary companion (double degenerate scenario; see e.g., Refs. [10–12]. In addition, sub-M_{ch} explosions have also been investigated (see e.g., Ref. [13]). Typically, in the explosion process the WD material is converted to iron-peak elements (iron, nickel and neighboring elements that form a prominent peak in the solar system abundance profile) and a smaller fraction of intermediate-mass elements (like Si, S, and Ca). However, it is only the radioactive decay of ⁵⁶Ni that powers the observed lightcurve.

In the generally accepted scenario, explosive carbon burning is ignited either at the center of the star or slightly off-center in a couple of ignition spots, depending on the still uncertain details of the previous evolution. Once the thermonuclear flame is ignited, there are two possible modes of the propagation: subsonic deflagration and supersonic detonation. A prompt ignition of the detonation flame is disfavored because the resulting nucleosynthesis yields fail to produce the strong intermediate-mass element features observed in SNe Ia and conflict with Galactic chemical evolution [14]. Thus, the explosion should start as a subsonic deflagration flame. The deflagration stage may last until the end of the explosion (the deflagration model; Ref. [15]). Alternately, it is also possible that the deflagration flame transforms into a detonation wave (the delayed detonation model, or the deflagration-detonation transition (DDT) model; Refs. [16–19]). Numerical models that parameterize the velocity of the burning front have been very successful at reproducing observations, the archetype being the W7 model of Nomoto et al. [15].

Our understanding of nucleosynthesis in these events is based on decades of modeling and nucleosynthesis calculations in spherically symmetric models (for example, Refs. [15, 19, 20]). In recent years, detailed estimates of the nucleosynthetic yields for multidimensional explosion models have become possible. Two dimensional axisymmetric simulations [21–24] and three dimensional simulations (see Refs. [25–27], and references therein) of exploding M_{ch} WD have provided new perspectives on thermonuclear supernova nucleosynthesis calculations. A minimal nuclear reaction network is usually included directly in the hydrodynamic simulations in order to compute the correct energy without drastically increasing the computing time. In direct numerical simulations of localized flame behavior (see e.g., Refs. [28–30]), this often takes the form of an \alpha-network composed of the 13 or 14 lightest nuclear species that have even numbers of protons and equal numbers of neutrons, linked by (\alpha ,\gamma ) and (\gamma ,\alpha ) reactions with a few additional critical heavy ion reactions, like {}_{}{}^{12}\mathrm{C}+{}_{}{}^{12}\mathrm{C} and {}_{}{}^{12}\mathrm{C}+{}_{}{}^{16}\mathrm{O}. In other cases, including most multi-dimensional simulations of the whole star, an even simpler flame scheme which tracks only 2 or 3 compositional stages is used. A much more extended network is therefore considered in a post-processing step using the tracer particles method. The nucleosynthesis experienced by each marker is calculated and the total yield for the supernova is computed as a sum over all the markers, after the decay of unstable isotopes. It is very encouraging that such models predict explosions with energies in the range of observed Type Ia supernovae, and that these models also predict light curves which fit well the observations. In addition, the nucleosynthesis products of these models are in reasonable agreement with expectations. The undeniable influence of SN Ia explosions on, among others, the chemical evolution of galaxies makes the quest for solid theoretical models and nucleosynthetic yields an urgent task.

2 Nucleosynthesis in thermonuclear supernovae: Present understanding

Over the last three decades, one-dimensional spherically symmetric models have been used to study the various channels that may give rise to a successful SN Ia in terms of nucleosynthesis, predicted spectra, and light curves. Much of this work was centered on the M_{ch} scenario wherein a C+O white dwarf accretes H or He from a binary companion [15, W7 model] and ignites explosive carbon burning just before it reaches the critical mass of M_{ch}=1.4M_{\odot }. The subsequent explosion produces enough ⁵⁶Ni (\simeq 0.6 M_{\odot }) and intermediate mass elements to reproduce normal SN Ia lightcurves and spectra, provided that the amount of C+O burned at any given density is suitably chosen. This can be achieved by parameterizing the thermonuclear flame speed and the density at which a transition to supersonic burning (detonation) occurs (see e.g., Ref. [31]). Moreover, some mixing of processed matter must to be assumed in order to fit the observed spectra. On the nucleosynthesis front, these models exhibit considerable production of the α-elements (the elements between Mg and Ca with even proton numbers) and iron-peak elements [32].

Since nuclear networks with hundreds of isotopes are very expensive in terms of CPU time and memory, it is possible to solve fully coupled systems of hydrodynamic and nuclear kinetic equations directly only in one-dimensional stellar codes. Instead in multi-dimensional codes large networks are relegated to a post-processing step. In the past, various authors implemented a method to calculate nucleosynthesis from multi-dimensional simulations (e.g. Ref. [35] for Type II SNe; [36] for Type Ia SNe). This tracer particle method is based on adding a “Lagrangian component” to the Eulerian scheme in the form of tracer (or marker) particles that correspond to fluid mass elements and are passively advected with the flow in course of the Eulerian calculation. To the particles are assigned masses and positions in such a way that a density profile reconstructed from their distribution resembles that of the underlying star. During the hydrodynamical simulation, they are advected by the flow, recording the history of thermodynamic conditions along their path (an example is shown in Fig. 1). The tracers are nevertheless considered massless in the sense that the mass they represent does not couple to the hydrodynamic flow via gravity or inertia. They are simply passively advected by the flow along streamlines and their temperature, density and energy histories are recorded. Given these values, for individual marker particles it is possible to calculate their nucleosynthetic evolution from their initial abundances and compute the total yields. The resolution, i.e., the number of tracer particles filling the volume of the star, has to be high enough in order to get accurate nucleosynthesis results. Seitenzahl et al. [37] performed a resolution study, demonstrating that in 2D SN Ia simulations with 80² tracer particles almost all isotopes up to iron-peak with abundances higher than about 10^-5 are reproduced with an accuracy better than 5%. More recently, Seitenzahl et al. [27] calculated tracer particle nucleosynthesis for different 3D SN Ia models. These authors, extrapolating the tracer resolution and yield convergence study from 2D to 3D, claimed that one million tracer particles (100 per axis) is sufficient to reliably predict the yields for the most abundant nuclides.

A detailed discussion regarding nucleosynthesis results for nuclei up to and including the iron-peak species can be found in Refs. [27, 36, 37]. We note that neutron-rich isotopes like ⁴⁸Ca, ⁵⁰Ti, ⁵⁵Mn, ⁵⁴Fe, ⁵⁸Ni depend only weakly on the particular choice of central density of the WD at the time of ignition. These isotopes are expected to be largely synthesized in nuclear statistical equilibrium condition in SNe Ia. Travaglio et al. [36], with a 3D deflagration model, found that they are strongly underproduced with respect to the standard 1D W7 model presented by Iwamoto et al. [19], but are in good agreement with the Brachwitz et al. [20] predictions. This agreement can be influenced by the differences in the electron capture rates adopted. We note that the most neutron-rich stable Fe-peak isotopes ⁵⁴Cr, ⁵⁸Fe, and ⁶⁴Ni are shielded by ⁵⁴Fe, ⁵⁸Ni, and ⁶⁴Zn from the Z=N line and thus require the highest neutronization for direct production. Consequently, these isotopes are the most sensitive to the central density. Seitenzahl et al. [27] underlined the link between a lower central density and a smaller overproduction of ⁵⁴Fe and ⁵⁸Ni. The very pronounced underproduction of ⁵⁸Fe and ⁶⁴Ni in all models does not pose a problem, as the s-process is the dominant source of these isotopes in galactic chemical evolution. The solar abundances of the remaining isotopes, ⁵⁹Co, ⁶⁰Ni, ⁶¹Ni, and ⁶²Ni, include significant contributions from explosive Si-burning in core collapse SNe, \alpha-rich freeze-out in both types of supernovae, and the s-process in AGB stars. Production factors for SNe Ia, ranging from a few to several tens of percent, are therefore also very reasonable. Seitenzahl et al. [27] interpret the large overproduction factor of ⁵⁴Cr (>3) obtained in their models as an indication that delayed detonation SNe Ia that ignite at high central densities are rare and constitute at most a small fraction of all SN Ia events. The iron-peak isotopes in SN Ia of normal brightness (where ⁵⁶Ni masses are around 0.6 M_{\odot }) are synthesized in the required proportions in Seitenzahl et al. [27], Travaglio et al. [36], indicating that delayed-detonations cannot be ruled out as the dominant SN Ia explosion channel based on solar iron-peak isotopic ratios.

The consequences for nucleosynthesis in SN Ia resulting from progenitor metallicity has been analyzed by several authors [27, 38, 39]. They all found that the ⁵⁶Ni production increases with decreasing initial ²²Ne, following a linear dependence due to the decreasing electron fraction Y_e, with a corresponding reduction in the stable iron-peak isotopes such as ⁵⁸Ni or ⁵⁴Fe. Travaglio et al. [39] found a variation in the ⁵⁶Ni mass of 25% in the metallicity range explored (covering 0.1Z_{\odot } to 3Z_{\odot }). The largest variation in ⁵⁶Ni occurs for metallicities greater than solar. Very little variation occurs in the unburned material ¹²C and ¹⁶O. The largest metallicity effect is seen in the α-elements. This effect must be carefully accounted for in chemical evolution calculations, in order to correctly evaluate the role of SN Ia in creating the Solar System abundance of nuclei.

Using an extended nuclear network with 1024 species from neutron and proton up to ²⁰⁹Bi combined with neutron, proton, and \alpha-induced reactions, recent work by Travaglio et al. [34] has demonstrated that SN Ia may be an important source for p-nuclei. p-nuclei are a class of about 35 neutron-deficient isotopes between ⁷⁴Se and ¹⁹⁶Hg. They are bypassed by the slow (s-process) and rapid (r-process) neutron capture processes, and are typically 10–1000 times less abundant than the corresponding s- and/or r-isotopes in the Solar System. The production mechanism that best fits the creation of these nuclei is photo-disintegration. The p-process nuclei seem to be formed in regions of temperatures around 1.8–3.3 \mathrm{GK}, however the nuclei do not completely photo-disintegrate, implying that some sort of temperature variance, gradient, or something else truncates the photo-disintegrations before reaching ordinary light elements. Few models, up to now, properly predict the abundances of several particular isotopes, specifically the lighter isotopes of molybdenum (92–94) and ruthenium (96–98). Laboratory experiments are being conducted to help determine the nuclear properties of these particular nuclei. One of the few investigations in this direction is that of Travaglio et al. [34]. These authors explored SNe Ia as p-process sources in the framework of 2D SN Ia models and found SNe Ia to be good astrophysical candidates for the production of these mysterious isotopes. But considerable work remains before the role of SN Ia (as well as SNII) in p-process production is fully understood.

The p-process nucleosynthesis occurs in SNe Ia only if there is a prior s-process enrichment. It is therefore essential to determine the source of the s-process enrichment in the exploding WD. In the single-degenerate progenitor model assumed here, there are two sources of s-enrichment: (i) during the AGB phase leading to the formation of the WD, thermal pulses occur during which s-isotopes are produced (TP-AGB phase; see, e.g., Ref. [41]) and (ii) thermal pulses can also occur as matter is accreted onto the WD, enriching the matter accumulating on the WD [42–45]. The first scenario would leave s-enriched matter on the surface of the WD in a layer of 0.1M_{\odot } prior to the accretion phase. During the accretion phase, this layer would be convectively mixed into the WD, diluting the s-seeds so that their abundances are too low to produce significant yields of p-isotopes. In the second scenario, s-process nucleosynthesis occurs in the H-rich matter accreted by the CO-WD, due to recurrent He-flashes [42], with neutrons mainly produced by the ¹³C(\alpha ,n){}^{16}O reaction. The conclusion obtained by these studies is that a flat s-seed distribution directly translates into a flat p-process distribution whose average production factor scales linearly with the adopted level of the s-seeds, giving the indication for a primary origin.

It has been recently demonstrated [40] that SNe Ia could be responsible for at least 50% of the all p-nuclei required for galactic chemical evolution, under the hypothesis that SN Ia are responsible for 2/3 of the solar ⁵⁶Fe, and assumption that delayed detonation model represents the typical SN Ia with a frequency of 70% [46]. Core-collapse supernovae are also expected to give an important contribution to the solar abundances of p-nuclei. Using a simple chemical evolution code and accounting for SNe Ia alone, production factors in comparison to the Solar System composition are shown in Fig. 2. Red labels identify the 35 nuclei classically identified as p-only. As one can see in this Figure, only a few nuclei originally ascribed to the p-only group (¹¹³In, ¹¹⁵Sn, ¹³⁸La, ¹⁵²Gd, ¹⁸⁰Ta) are far below the average production factor of the other p-nuclei. This indicates a different nucleosynthetic stellar origin for these species. As discussed by Dillmann et al. [47] and Nemeth et al. [48], ¹¹³In and ¹¹⁵Sn receive important contributions from delayed r-process decay chains. Woosley et al. [49] demonstrated that neutrino-nucleus interactions can contribute appreciably to the synthesis of ¹³⁸La in the neon shell of core collapse supernovae. ¹⁵²Gd is now thought predominantly to have an s-process origin [50, 51]. Likewise, ¹⁶⁴Er is also of predominantly s-process origin, driven by the β-decay channel at ¹⁶³Dy, which becomes unstable at stellar temperature [52]. ¹⁸⁰Ta, as discussed by Ref. [53], receives an important contribution from the s-process due to the branching at ¹⁷⁹Hf, a stable isotope that becomes unstable at stellar temperatures. Neutrino–Nucleus interactions may also substantially feed ¹⁸⁰Ta. The relative low abundance of the p-only nuclide ¹⁵⁸Dy should be analyzed in the framework of present nuclear uncertainties of the γ-proceses.

The analysis of the yields is fundamental in order to understand the role of SNe Ia in the chemical enrichment of the Galaxy, but how the yields are distributed in mass and velocity space is also important in order to compare the results with observed supernova spectra or light curves. Such analysis has been performed for the multidimensional models quoted above [27, 36, 54]. Travaglio et al. [36], with 3D pure deflagration models, took selected isotope mass fractions at the end of the simulation (at about 1 s after the explosion) and noted that at the lowest velocities (1000 km/s), the dominant component is represented by the unburned material (i.e., ¹²C and ¹⁶O) for a floating bubble model, or by ⁵⁶Fe for a centrally ignited model. In contrast, at the highest velocities (>10 000 km/s) the unburned material dominates in the centrally ignited model, but is comparable to the ⁵⁶Fe component in the floating-bubble high resolution model. Maximum velocities reached are about 12000 km/s in both cases. These models are not yet in homologous expansion, i.e., pressure and gravity still play a role, changing the velocity distribution and possibly also the density, therefore the velocity distribution of the elements has to be taken with care (see Fig. 3). Seitenzahl et al. [27] with delayed detonation 3D models, found that stable iron-peak nuclei are not present at the lowest velocities as predicted by 1D models (e.g., Refs. [15, 16]), but rather at intermediate velocities (3000–10000 km/s). This agrees well with the results found by Maeda et al. [54], with typical velocities of 5000–10000 km/s for the deflagration ash in their two-dimensional DDT model. Stehle et al. [55] find substantial amounts of stable iron out to velocities of about 9000 km/s from analysis of the spectra of SN 2002bo. Seitenzahl et al. [27] also find that the central regions, which were burned by the detonation at high densities to NSE, form homogeneous ⁵⁶Ni clumps.

For the 3D deflagration models presented by Travaglio et al. [36], late time synthetic spectra have been performed by Kozma et al. [56]. They found that the model spectra after 300 to 500 days are in good agreement with the observed Fe II-III features. However, they also show O I and C I lines stronger than those observed in late time spectra. The oxygen and carbon emission originates from the low-velocity unburned material in the central regions of these models. Agreement between the models and observations requires that only a small mass of unburned material may be left in the center after the explosion. This may pose a problem for pure deflagration models. For the same 3D deflagration models of Travaglio et al. [36], synthetic bolometric and broad-band UBVRI light curves have been calculated by Blinnikov et al. [57]. These results show that Chandrasekhar mass models exploded by pure deflagration produce both UBV light curves and photospheric expansion velocities that match well with observed weak to normal SNe Ia. As the majority of the flux during maximum light is emitted within these passbands, this is an encouraging result. Moreover, it is these passbands that are very important in cosmological applications of SNe Ia. It is clear that the models require some improvement in explaining the shapes of the near infrared light curves and in explaining fast spectral features that are observed in many normal events. The bolometric light curves calculated from these deflagration models also evolve slightly slower than what is indicated from observations. These discrepancies hint at the necessity of producing faster moving ejecta and somewhat less mixed chemical composition. The latest deflagration models (see Ref. [58], and references therein) are promising in this respect.

3 Prospects for improvements in the simulation of thermonuclear supernova nucleosynthesis

The utilization of multi-dimensional simulations for thermonuclear supernovae has greatly improved our understanding of the flame propagation in these events. However, in order to offset the increased computational cost of expanded dimensionality, this progress has come at the cost of decreased sophistication in the treatment of the nuclear composition. As discussed in other articles in this volume, multi-dimensional whole star simulations generally employ a three level flame tracking algorithm, dividing the star into regions of carbon & oxygen, magnesium & silicon and iron, nickel & \alpha-particles. In contrast, spherically symmetric models of the entire star have for more than a decade employed much large reaction networks. For example, the models of Höflich (see, e.g., Refs. [59, 60]) employ 218 nuclear species self-consistently within the simulations, while those of Iwamoto et al. [19] and Brachwitz et al. [20] contained 299 nuclear species.

As documented in the previous section, to ameliorate the small networks used in multi-dimensional simulations, the post-processing approach is widely used in supernova nucleosynthesis calculations. In this approach, a simplified approximation is used to calculate the compositional change and resulting rate of nuclear energy generation within the hydrodynamic model. Subsequently, the temperature and density histories for individual mass elements from the supernova model are then used as input for separate, larger nucleosynthesis calculations. Performing post-processing calculations based on one-dimensional models is relatively straightforward, since most one-dimensional simulations are Lagrangian. Thus the needed temporal histories of temperature and density are simply those of the individual Lagrangian mass elements. However, in multi-dimensions, non-smooth fluid motions result in highly tangled Lagrangian grids. As a result, Eulerian hydrodynamics, where the discretization occurs in space rather than mass, is used to perform most multi-dimensional stellar astrophysics simulations. Because Eulerian codes use spatial discretization, the Lagrangian thermodynamic histories which are a natural result in a Lagrangian code are unavailable. Instead passive tracer particles must be employed. Post-Processing nucleosynthesis of tracer particles have been utilized in simulations of core-collapse supernovae since Nagataki et al. [35] and in simulations of thermonuclear supernovae since Travaglio et al. [36].

Despite the wide use of post-processing in supernova nucleosynthesis calculations, this approach has significant drawbacks. The primary limitations of a post-processing approach are i) the accuracy of the energy generation rate provided by the approximation included within the hydrodynamics, ii) limitations in computing the neutron-richness of the matter and iii) an underestimate of the effects of mixing. Calculation of the energetically important stages (carbon, oxygen and silicon burning as well as the freezeout from Nuclear Statistical Equilibrium) is the motivation for spherically symmetric simulations to employ networks of 150–300 nuclei, depending on the desired accuracy and range in neutron-richness. To understand the path forward toward the utilization of similar sized networks in multi-dimensional simulations of thermonuclear supernovae we must analyze the cost of thermonuclear reaction networks before we discuss ways to reduce this cost.

3.1 Solving thermonuclear reaction networks

The details of nuclear reaction network calculations have been extensively reviewed elsewhere (see, e.g., Refs. [61–65]), so we will include here only a brief outline, primarily to define terms to aid the following discussion. The nuclear population in each cell of a simulation is frequently described in terms of the nuclear abundance of each species, Y=n/\rho N_A, where n is the number density of the species, ρ is the total mass density, and N_A is Avagadro’s number. For a nucleus with atomic mass number A, AY\equiv X is commonly referred to as the mass fraction of this nucleus and the sum of all mass fractions (\sum A_i{Y}_i) is by construction equal to 1. The expression \sum A_i{Y}_i=1 is in reality an expression of nucleon number conservation rather than mass conservation, since the nuclear binding energy can approach 1% of the nuclear mass. In fact the conversion of mass into energy is one of the primary reasons to include nuclear reactions in an astrophysical model. The equation of charge conservation becomes \sum Z_i{Y}_i=Y_e, where Y_e(=n_e/\rho N_A) is the electron abundance, which provides a measure of the neutron-richness of the matter.

If we operator-split spatial advection from nuclear transmutations, the evolution of the nuclear abundances in a thermonuclear reaction network takes the form of a set of ordinary differential equations for the time derivatives \dot{Y} which depend only on nuclear reaction rates and local thermodynamic conditions.

\begin{array}{}\text{(1)}& \begin{array}{cc}{\dot{Y}}_i=\hfill & {\displaystyle \sum _j}{\mathcal{N}}_j^i{\lambda }_j{Y}_j+{\displaystyle \sum _{j,k}}{\mathcal{N}}_{j,k}^i\rho N_A{⟨\sigma v⟩}_{j,k}{Y}_j{Y}_k\hfill \\ \hfill & +{\displaystyle \sum _{j,k,l}}{\mathcal{N}}_{j,k,l}^i{\rho }^2{N}_A^2{⟨\sigma v⟩}_{j,k,l}{Y}_j{Y}_k{Y}_l\hfill \end{array}\end{array}

where the decay constants (\lambda ) and velocity-integrated reaction cross sections (⟨\sigma v⟩) contain the essential nuclear data. The \mathcal{N}s provide for proper accounting of numbers of nuclei with {\mathcal{N}}_j^i=N_i, {\mathcal{N}}_{j,k}^i=N_i/{\prod }_{m=1}^{n_m}|N_{j_m}|!, and {\mathcal{N}}_{j,k,l}^i=N_i/{\prod }_{m=1}^{n_m}|N_{j_m}|!, where the N_i^{\prime }s are positive (or negative) numbers that specify how many particles of species i are created (or destroyed) in a reaction, while the denominators, including factorials, run over the n_m different species destroyed in the reaction and avoid double counting of the number of reactions when identical particles react with each other.

In principle, the initial value problem presented by the nuclear reaction network can be solved by any of a large number of methods discussed in the literature. However the physical nature of the problem, reflected in the decay constants and velocity-integrated reaction cross sections, greatly restricts the optimal choice. The large number of reactions display a wide range of reaction timescales, with the mean lifetime of particle j against destruction by reaction with particle k given by

\begin{array}{}\text{(2)}& {\tau }_k(j)={\displaystyle \frac{1}{{⟨\sigma v⟩}_{j,k}{Y}_k}}\text{or}{\tau }_j(j)={\displaystyle \frac{1}{{\lambda }_j}}\end{array}

for a decay. Systems whose solutions depend on a wide range of timescales are termed stiff. Gear [66] demonstrated that even a single equation can be stiff if it has both rapidly and slowly varying components. Practically, stiffness occurs when the limitation of the timestep size is due to numerical stability rather than accuracy. A more rigorous definition [67] is that a system of equations is stiff if the ratio of the real part of the largest eigenvalue {\lambda }_{max} of the Jacobian matrix of the system of equations is much larger than the real part of the smallest eigenvalue {\lambda }_{min}. Because of the wide range in timescales between strong, electromagnetic and weak reactions, {\lambda }_{max}/{\lambda }_{min}>{10}^{15} is not uncommon in astrophysics, making thermonuclear reaction networks extraordinarily stiff.

Nucleosynthesis calculations belong to the more general field of reactive flows, and therefore share some characteristics with related terrestrial fields. In particular, chemical kinetics, the study of the evolution of chemical abundances, is an important part of atmospheric and combustion physics and produces sets of equations much like Eq. (1) (see Ref. [68] for a good introduction). These chemical kinetics systems are known for their stiffness and a great deal of effort has been expended on developing methods to solve these equations. Many of the considerations for the choice of solution method for chemical kinetics also apply to nucleosynthesis calculations. In both cases, temporal integration of the reaction rate equations is broken up into short intervals because of the need to update the hydrodynamics variables. This favors one step, self starting algorithms. Because abundances must be tracked for a large number of computational cells (hundreds to thousands for one dimensional models, millions to perhaps billions for the coming generation of three dimensional models), memory storage concerns favor low order methods since they do not require the storage of as much data from prior steps. In any event, both the errors in fluid dynamics and in the reaction rates are typically a few percent or more, so the greater precision of these higher order methods often does not result in greater accuracy.

For a set of nuclear abundances \overrightarrow{Y}, one can calculate the time derivatives of the abundances, \dot{\overrightarrow{Y}} using Eq. (1). The desired solution is the abundance at a future time, \overrightarrow{Y}(t+\mathrm{\Delta }t), where \mathrm{\Delta }t is the network timestep. Since coupling with hydrodynamics favors low order, one step methods, many nucleosynthesis calculations use the simple finite difference prescription

\begin{array}{}\text{(3)}& {\displaystyle \frac{\overrightarrow{Y}(t+\mathrm{\Delta }t)-\overrightarrow{Y}(t)}{\mathrm{\Delta }t}}=(1-\Theta )\dot{\overrightarrow{Y}}(t+\mathrm{\Delta }t)+\Theta \dot{\overrightarrow{Y}}(t)\end{array}

With \Theta =1, Eq. (3) becomes the explicit Euler method while for \Theta =0 it is the implicit backward Euler method, both of which are first order accurate. For \Theta =1/2, Eq. (3) is the semi-implicit trapezoidal method, which is second order accurate. For the stiff set of non-linear differential equations which form most nuclear reaction networks, the implicit backward Euler method, introduced to astrophysics by Arnett and Truran [69], is generally most successful, though the trapezoidal method has been used in Big Bang nucleosynthesis calculations [70], where coupling to hydrodynamics is less important. Solving the fully implicit version of Eq. (3) is equivalent to finding the zeros of the set of equations

\begin{array}{}\text{(4)}& \overrightarrow{\mathcal{Z}}(t+\mathrm{\Delta }t)\equiv {\displaystyle \frac{\overrightarrow{Y}(t+\mathrm{\Delta }t)-\overrightarrow{Y}(t)}{\mathrm{\Delta }t}}-\dot{\overrightarrow{Y}}(t+\mathrm{\Delta }t)=0\end{array}

This is done using the Newton–Raphson method (see, e.g., Ref. [71]), which is based on the Taylor series expansion of \overrightarrow{\mathcal{Z}}(t+\mathrm{\Delta }t), with the trial change in abundances given by

\begin{array}{}\text{(5)}& \mathrm{\Delta }\overrightarrow{Y}={\left({\displaystyle \frac{\partial \overrightarrow{\mathcal{Z}}(t+\mathrm{\Delta }t)}{\partial \overrightarrow{Y}(t+\mathrm{\Delta }t)}}\right)}^{-1}\overrightarrow{\mathcal{Z}}\end{array}

where \partial \overrightarrow{\mathcal{Z}}/\partial \overrightarrow{Y} is the Jacobian of \overrightarrow{\mathcal{Z}}. Iteration continues until \overrightarrow{Y}(t+\mathrm{\Delta }t) converges. Timmes [72] has advanced higher order semi-implicit methods (e.g. Bader–Deuflhard) for network solution. These also involve a matrix solution, thus they share the same computational bottleneck as the backward Euler method.

3.1.1 The structure of the Jacobian matrix

For networks of the size we desire to include in multi-dimensional simulations of thermonuclear supernovae, the implicit and semi-implicit methods require solution of moderate-sized (N=100–300) matrix equations. Since general solution of a dense matrix scales as O(N^3), this can make these large networks progressively much more expensive. While in principal, every species reacts with each of the hundreds of others, resulting in a dense Jacobian matrix, in practice it is possible to neglect most of these reactions. Because of the Z_i{Z}_j dependence of the repulsive Coulomb term in the nuclear potential, captures of free neutrons and isotopes of H and He on heavy nuclei occur much faster than fusions of heavier nuclei. Furthermore, reactions involving secondary isotopes of H (deuterium and tritium) and He are neglectable in a thermonuclear supernova. Likewise, photodisintegrations tend to eject free nucleons or \alpha-particles. Thus, with a few important exceptions, for each nucleus we need only consider twelve reactions linking it to its nuclear neighbors by the capture of an n,p,\alpha or γ and release a different one of these four. The exceptions to this rule are the few heavy ion reactions important for burning stages like carbon and oxygen burning where the dearth of light nuclei cause the heavy ion collisions to dominate.

Figure 5 demonstrates the sparseness of the resulting Jacobian matrix, for a reaction network with a desirable 300 nuclei. Of the 90 000 matrix elements, less than 5000 are non-zero. In terms of the standard forms for sparse matrices, this Jacobian is best described as doubly bordered, band diagonal. With a border width, \mathrm{\Delta }B, of 45 necessary to include the heavy ion reactions among {}_{}{}^{12}\mathrm{C}, {}_{}{}^{16}\mathrm{O} and {}_{}{}^{20}\mathrm{Ne} along with the free neutrons, protons and \alpha-particles and a band diagonal width, \mathrm{\Delta }D, of 54, even this sparse form includes almost 50 000 elements. With solution of the matrix equation consuming 80% or more of the computational time, there is a clear advantage for solvers that take advantage of the sparseness of the Jacobian. To date best results for small () matrices are obtained with machine optimized dense solvers (e.g. LAPACK) or matrix specific solvers generated by symbolic processing [73, 74]. For large matrices, generalized sparse solvers, both custom built and from software libraries, perform best (see, e.g., Ref. [72]), though the benefit is small for systems of the size one would like to include with a thermonuclear supernova model.

A potential numerical problem with the solution of Eq. (4) is the singularity of the Jacobian matrix, \partial \overrightarrow{\mathcal{Z}}(t+\mathrm{\Delta }t)/\partial \overrightarrow{Y}(t+\mathrm{\Delta }t). From Eq. (4), the individual matrix elements of the Jacobian have the form

\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial {\mathcal{Z}}_i}{\partial Y_j}}={\displaystyle \frac{{\delta }_{ij}}{\mathrm{\Delta }t}}-{\displaystyle \frac{\partial {\dot{Y}}_i}{\partial Y_j}}={\displaystyle \frac{{\delta }_{ij}}{\mathrm{\Delta }t}}-{\displaystyle \sum \frac{1}{{\tau }_j(i)}}\end{array}

where {\delta }_{ij} is the Kronecker delta, and {\tau }_j(i) is defined in Eq. (2). The sum accounts for the fact that there may be more than one reaction by which nucleus j is involved in the creation or destruction of nucleus i. Along the diagonal of the Jacobian, there are two competing terms, 1/\mathrm{\Delta }t and \sum 1/{\tau }_i(i). This sum is over all reactions which destroy nucleus i, and is dominated by the fastest reactions. As a result, \sum 1/{\tau }_i(i) can be orders of magnitude larger than the reciprocal of the desired timestep, 1/\mathrm{\Delta }t. This is especially a problem near equilibrium, where both destruction and the balancing production timescales are very short in comparison to the preferred timestep size, resulting in differences close to the numerical accuracy (i.e., 14 or more orders of magnitude). In such cases, the term 1/\mathrm{\Delta }t is numerically neglected, leading to numerically singular matrices.

3.1.2 Equilibria in nuclear astrophysics

While the approach to equilibrium can play havoc with the numerics of the matrix solution, it also offers a way around this problem. As is the case in many disciplines, equilibrium expressions can be employed to simplify nuclear abundance calculations. The simplest application of equilibrium is the expression for the relative abundances of two species when the reactions linking them are equilibrated. For reactions of the form i+j\leftrightarrow k+\mathrm{\ell }, it is four abundances that are interrelated

\begin{array}{}\text{(7)}& {\displaystyle \frac{Y_k{Y}_{\mathrm{\ell }}}{Y_i{Y}_j}}={\displaystyle \frac{{⟨\sigma v⟩}_{i,j}}{{⟨\sigma v⟩}_{k,\mathrm{\ell }}}}={\displaystyle \frac{G_k{G}_{\mathrm{\ell }}}{G_i{G}_j}}{\displaystyle \frac{1+{\delta }_{ij}}{1+{\delta }_{k\mathrm{\ell }}}}{\left({\displaystyle \frac{{\mu }_{k\mathrm{\ell }}}{{\mu }_{ij}}}\right)}^{\frac{3}{2}}\exp \left({\displaystyle \frac{Q}{kT}}\right)\end{array}

The final equality in Eq. (7) is the result of the detailed balance relation for a time-reversible processes (see, for example, Refs. [75–77]. From this relation, the relative abundances of the nuclear species depend only on their properties (G_i is the level density (2J_i+1), {\mu }_{ij} is the reduced mass, and Q is the energy released by the forward reaction, m_k+m_{\mathrm{\ell }}-m_i+m_j, where m_i the nuclear mass of species i).

In the highest temperature regions of a thermonuclear supernova, the fast strong and electromagnetic reactions can reach equilibrium while those involving the weak nuclear force do not. Since the weak reactions are not equilibrated, the resulting Nuclear Statistical Equilibrium (NSE) requires monitoring of weak reaction activity. The expression for NSE is commonly derived using either chemical potentials or detailed balance (see, e.g., Refs. [75, 78, 79]). For a nucleus {}_{}{}^{A}Z, composed of Z protons and N=(A-Z) neutrons, in equilibrium with these free nucleons, the chemical potential of {}_{}{}^{A}Z can be expressed in terms of the chemical potentials of the free nucleons

\begin{array}{}\text{(8)}& {\mu }_{Z,A}=Z{\mu }_p+N{\mu }_n\end{array}

For a collection of particles obeying Boltzmann statistics, the chemical potential, including rest mass, of each species is given by

\begin{array}{}\text{(9)}& {\mu }_i=m_i{c}^2+k_{B}T\ln \left[\rho N_A{\displaystyle \frac{Y_i}{G_i}}{\left({\displaystyle \frac{2\pi {\mathrm{\hslash }}^2}{m_i{k}_{\mathrm{B}}T}}\right)}^{\frac{3}{2}}\right]\end{array}

Substituting Eq. (9) into Eq. (8) allows derivation of an expression for the abundance of every nuclear species in terms of the abundances of the free protons (Y_p) and neutrons (Y_n),

\begin{array}{}\text{(10)}& \begin{array}{cc}Y{(}^{A}Z)=\hfill & {\displaystyle \frac{G{(}^{A}Z)}{2^A}}{\left({\displaystyle \frac{\rho N_A}{\Theta }}\right)}^{A-1}{A}^{\frac{3}{2}}\exp \left({\displaystyle \frac{B{(}^{A}Z)}{k_{\mathrm{B}}T}}\right)\hfill \\ \hfill & \cdot Y_n^{}{}_{}{}^{N}Y_p^{}{}_{}{}^Z\equiv C{(}^{A}Z)Y_n^{}{}_{}{}^{N}Y_p^{}{}_{}{}^Z\hfill \end{array}\end{array}

where G{(}^{A}Z) and B{(}^{A}Z) are the partition function and binding energy of the nucleus {}_{}{}^{A}Z, N_A is Avagadro’s number, k_B is Boltzmann’s constant, ρ and T are the density and temperature of the plasma, and \Theta is given by

\begin{array}{}\text{(11)}& \Theta ={\left({\displaystyle \frac{m_u{k}_{\mathrm{B}}T}{2\pi {\mathrm{\hslash }}^2}}\right)}^{3/2}\end{array}

Thus abundances of all nuclear species can be expressed as functions of two. Two constraints are required to uniquely define these abundances. Strong and electromagnetic reactions do not alter the total abundances of protons and neutrons, thus the total proton abundance, \sum ZY and the total neutron abundance, \sum NY, provide the needed constraints. More frequently, total nucleon number conservation, (\sum (N+Z)Y=\sum AY=1) is used as one constraint. This is colloquially referred to as mass conservation, since A, the number of nucleons in the nucleus, is also referred to as the atomic mass number. The second constraint then records the amount of weak reaction activity and is often expressed in terms of the total proton abundance, \sum ZY, which charge conservation requires equal the electron abundance, Y_e. Thus the nuclear abundances are uniquely determined for a given (T,\rho ,Y_e).

Where applicable, NSE offers significant advantages, since hundreds of abundances are uniquely defined by the thermodynamic conditions and a single measure of the weak interaction history or the degree of neutronization. Computationally, this reduction in the number of independent variables greatly reduces the cost of nuclear abundance evolution. Because there are fewer variables to follow within a hydrodynamic model, the memory footprint of the nuclear abundances is also reduced, an issue of importance in modern multi-dimensional models of supernovae. These advantages allow Seitenzahl et al. [80] to include the effects of electron capture on 443 species within multi-dimensional whole star simulations, at least in regions where NSE applies.

As with any equilibrium distribution, there are limitations on the applicability of NSE. The first requirement for NSE to provide a good estimate of the nuclear abundances is that the temperature be sufficient for the endoergic reaction of each reaction pair to occur. Since for all particle-stable nuclei between the proton and neutron drip lines (with the exception of nuclei unstable against alpha decay), the photodisintegrations are endoergic, with typical Q-values among (β) stable nuclei of 8–12 \mathrm{MeV}, this requirement reduces to T>3\mathrm{GK}. While this requirement is necessary, it is not sufficient. In the rapidly changing conditions that occur within supernovae, even when this condition is met, hydrodynamic flows can cause thermodynamic conditions to change on timescales more rapid than NSE can be maintained, breaking NSE down first between {}_{}{}^4\mathrm{He} and {}_{}{}^{12}\mathrm{C}, at T\sim 6\mathrm{GK} [81] and later between the species near silicon and the iron-peak nuclei, at T\sim 4\mathrm{GK} [82]. The breakdown of NSE leaves large groups of nuclei that remain in mutual equilibrium, termed Quasi-Equilibrium [QSE; 83, 84], linked by series of reactions through out-of-equilibrium species. As the temperature continues to decline, these QSE groups fragment until eventually no reactions are in equilibrium. Note, QSE is also referred to as Nuclear Statistical Quasi-Equilibrium (NSQE) by some authors [85]. The relative abundance of a nucleus, {}_{}{}^{A}Z, in one of these QSE groups, compared to some focal species {}_{}{}^{A^{\prime }}Z_{}^{\prime } is given by

\begin{array}{}\text{(12)}& Y_{QSE{,}^{A^{\prime }}{Z}^{\prime }}{(}^{A}Z)={\displaystyle \frac{C{(}^{A}Z)}{C{(}^{A^{\prime }}{Z}^{\prime })}}Y{(}^{A^{\prime }}{Z}^{\prime })Y_p^{Z-Z^{\prime }}{Y}_n^{N-N^{\prime }}\end{array}

where C{(}^{A}Z) is defined in Eq. (10). Through the QSE relation, the abundances of large groups of nuclei can be expressed in terms of those of a single group member and the free nucleons.

3.2 The quasi-equilibrium-reduced reaction network

If one can exploit this reduction in the number of species that must be tracked, the size and therefore cost of the nuclear reaction network can be reduced without reducing the network’s accuracy. In place of tracking the abundances of each group member, we track the abundance of the group as a whole,

\begin{array}{}\text{(13)}& Y_{group}{(}^{A^{\prime }}{Z}^{\prime })={\displaystyle \sum _{{}_{}{}^{A}Z{\in }^{A^{\prime }}{Z}^{\prime }group}}Y{(}^{A}Z)\end{array}

Under conditions where large QSE groups exist, the rapid series of reactions within the group which establishes QSE also determines the free proton and neutron abundances, causing a QSE group involving the light nuclei (up to He) to form. The abundances in this light QSE group obey the NSE relation [Eq. (10)], except that the free nucleon abundances are not determined by global constraints on the total neutron and proton abundance, but instead by constraints for these populations within the QSE groups. This leads to additional abundance equations that must be evolved,

\begin{array}{}\text{(14)}& \begin{array}{cc}{Y}_{NG}=\hfill & {\displaystyle \sum _{{}_{}{}^{A}Z\in Ltgroup}}NY+{\displaystyle \sum _{{}_{}{}^{A}Z{\in }^{A^{\prime }}{Z}^{\prime }group}}(N-N^{\prime })Y+\mathrm{\cdots }\hfill \\ Y_{ZG}=\hfill & {\displaystyle \sum _{{}_{}{}^{A}Z\in Ltgroup}}ZY+{\displaystyle \sum _{{}_{}{}^{A}Z{\in }^{A^{\prime }}{Z}^{\prime }group}}(Z-Z^{\prime })Y+\mathrm{\cdots }\hfill \end{array}\end{array}

with the sums carried out over each QSE group. In the special case of the \alpha-network, free neutrons and protons are absent leaving \alpha-particles as the only light reactants. Therefore in the QSE-reduced \alpha-network, Y_{NG} and Y_{ZG} are replaced by a single equation

\begin{array}{}\text{(15)}& Y_{\alpha G}=Y_{\alpha }+{\displaystyle \sum _{{}_{}{}^{A}Z{\in }^{A^{\prime }}{Z}^{\prime }group}}{\displaystyle \frac{A-A^{\prime }}{4}}Y+\mathrm{\cdots }\end{array}

A system of abundance equations including Y_{NG}, Y_{ZG} and Y_{group}{(}^{A^{\prime }}{Z}^{\prime }) for each QSE group, along with the abundances of the species which are not members of any QSE groups, can be evolved in the same fashion as the conventional network [see Eq. (3)]. One complication is that the \dot{Y} terms from Eq. (1) that appear on the righthand side of Eq. (3) depend on the abundances of individual species Y{(}^{A}Z) rather than these group abundances. This requires evolution of the relative abundances with the QSE groups, or equivalently the abundances of the free nucleons. In the simpler case of the \alpha-network, the iso7 network of Timmes et al. [86] utilizes a direct solution of the free \alpha-particle abundance, the only free particle abundance in that simplified system. This proves sufficient for carbon and oxygen burning, but is less successful for silicon burning, where the reaction flow between QSE groups centered on silicon and the iron-peak species is important. The \alpha7 network of Hix et al. [87] takes a two step approach wherein the QSE group abundances (Y_{\alpha G}, Y_{group}({}_{}{}^{28}\mathrm{Si}), Y_{group}({}_{}{}^{56}\mathrm{Ni}),\mathrm{\cdots } ) are evolved, then the individual abundances within the groups are determined in each iteration. Y_{\alpha G} and one group abundance per QSE group provide sufficient constraints to iteratively solve for Y_{\alpha }, and Y{(}^{A^{\prime }}{Z}^{\prime }) for each group, thereby determining the abundances of all QSE group species. This approach proves more successful during silicon burning at reproducing the abundances calculated from a complete \alpha-network.

Hix et al. [88] extended this approach to the general case, wherein Y_{NG}, Y_{ZG} are evolved in addition to group abundances for silicon and iron-peak QSE groups, providing sufficient constraints to solve for Y_n, Y_p, Y({}_{}{}^{28}\mathrm{Si}) and Y({}_{}{}^{56}\mathrm{Ni}). This again allows the calculation of abundances for all species within the light, silicon and iron-peak groups and accurate determination of reaction rates that are needed to evolve the next iteration. Figure 6 illustrates the accuracy with which the QSE-reduced network replicates the results of a complete reaction network, showing the evolution of key species and the energy generation rate for a case of incomplete silicon burning as might occur in the outer layers of a thermonuclear supernova. The energy generation figure includes results from two QSE-reduced network calculations, one (LQSE) assuming large QSE groups and, as a result, relatively few non-QSE species, and the other (SQSE) assuming small QSE groups, leaving more species outside of the QSE groups. There is a trade-off of greater speed in the LQSE case for moderately greater accuracy in the SQSE case.

It is important to note that unlike a conventional network, where abundance changes naturally cease as the temperature declines, the QSE relations continue to alter the relative abundances of the species in the QSE groups even after reaction flows between the QSE groups have ceased. As Fig. 6 illustrates, once the temperature drops below 3.5 \mathrm{GK} (the temperature at which QSE physically breaks down within the QSE groups), one must either continue the evolution in a conventional network or freeze the abundances. Thus the QSE-reduced network does not represent global replacement for a conventional network, as it is inaccurate for following the abundance evolution at low temperatures. However, by significantly reducing the computational cost of the most expensive parts of a supernova simulation, the QSE-reduced network will be a valuable tool in accelerating the calculation of multi-dimensional models which track the desired hundreds of species. Hix et al. [88], using 3 fixed QSE groups, demonstrated that the QSE-reduced network can be 5–10 times faster for these expensive nucleosynthesis stages (depending on the size of the network and assumed QSE groups, as well as the processor and matrix solver packages used). Parete-Koon et al. [89, 90] have shown that a further factor of 2 can be realized through dynamic changes in the QSE group assignments, primarily by taking advantage of the extended period where the composition is determined by 2 QSE groups separated by a few species between helium and magnesium.

3.3 Explicit methods to solve thermonuclear reaction networks

The success of the QSE-reduced network in reducing the cost of thermonuclear evolution by reducing the cost of building and solving the Jacobian matrix is simultaneously suggestive of a general approach to such speed improvements and highly specialized to the conditions under which QSE applies. A more generally applicable method to evolve nuclear abundances without building and solving a matrix is highly desirable. Such methods are widely used in other disciplines to solve stiff systems of equations that bear algebraic resemblance to the equations of a nuclear reaction network. Two examples which have been applied with some success to astrophysical nuclear reaction networks are the asymptotic method [91, 92] and the quasi-steady state (QSS) method [93, 94]. Here we will discuss the simpler asymptotic method, though the QSS is generally moderately faster.

3.3.1 The asymptotic method

Asymptotic methods have been widely studied in other fields (see, e.g., Refs. [95, 96]). They exploit the fact that exponential terms which impose stiffness often decay rapidly, limiting the impact errors made in their approximation may have. The primary properties we note here are that the asymptotic method is fully explicit and never requires the creation or solution of a Jacobian matrix, they are A-stable, and explicit asymptotic methods such as described here scale linearly with network size. These properties are shared by the more recently developed quasi-steady state (QSS) method [93, 94].

The system of differential equations shown in Eq. (1) to represent a nuclear reaction network can be written as

\begin{array}{}\text{(16)}& \dot{Y_i}=F_i(\overrightarrow{Y},t)={\displaystyle \sum _j}{F}_{ij}(t)\end{array}

where F_i(\overrightarrow{Y},t) is the net flux contributing to population i, t is the time, and F_{ij} is the individual flux between population i and population j. For the purposes of the asymptotic approximation, it is convenient to group these fluxes into two new variables: one which is the sum of all fluxes which add to the population of i (F_i^{+}) and one which is the sum of all fluxes which decrease population i (F_i^{-}). These variables are chosen with a sign convention such that they are both non-negative. Adopting this notation, Eq. (16) can be expressed as

\begin{array}{}\text{(17)}& \dot{Y_i}=F_i^{+}(t)-F_i^{-}(t)\end{array}

For astrophysical networks the flux out of species i is proportional to Y_i, such that

\begin{array}{}\text{(18)}& F_i^{-}(t)=[k_1^i(t)+k_2^i(t)+\mathrm{\cdots }+k_m^i(t)]Y_i(t)\end{array}

and equivalently

\begin{array}{}\text{(19)}& F_i^{-}(t)=k^i(t)Y_i(t)\end{array}

where the k_n^i are rate parameters (either {\lambda }_n or {⟨\sigma v⟩}_{in}{Y}_n, in units of inverse time) for each of the m reactions that can remove population from Y_i, and kⁱ is an effective rate parameter summing all of the other rate parameters for species i. Note that kⁱ may retain a dependence on populations Y_n for reactions with 2 or more reactants. Equation (19) allows Eq. (17) to be written as

\begin{array}{}\text{(20)}& Y_i(t)={\displaystyle \frac{F_i^{+}(t)-\dot{Y_i}}{k^i(t)}}\end{array}

Taking a finite-difference approximation of Eq. (20) at timestep t_n produces

\begin{array}{}\text{(21)}& Y_i(t_n)=F^{+}(t_n)k^i(t_n)-{\displaystyle \frac{\dot{Y_i}(t_n)}{k^i(t_n)}}\end{array}

Define an asymptotic limit for species i such that F_i^{+}(t)\simeq F_i^{-}(t), which implies \dot{Y_i}(t)\simeq 0. Therefore Eq. (21) can be rewritten as

\begin{array}{}\text{(22)}& Y_i(t_n)={\displaystyle \frac{F_i^{+}(t_n)}{k^i(t_n)}}\end{array}

which is a first approximation to the asymptotic abundance Y_i(t) at time t=t_n.

This approximation can be improved by the addition of a correction term. We can approximate \dot{Y_i}(t_n) by the backward difference

\begin{array}{}\text{(23)}& \dot{Y_i}(t_n)\simeq {\displaystyle \frac{Y_i(t_n)-Y_i(t_{n-1})}{\mathrm{\Delta }t}}\end{array}

where \mathrm{\Delta }t is the timestep (t_n-t_{n-1}), Y_i(t_n) is the approximated abundance evaluated at timestep n, and Y_i(t_{n-1}) is the same evaluated at timestep n-1. Substituting Eq. (22) into Eq. (23) gives

\begin{array}{}\text{(24)}& \dot{Y_i}(t_n)={\displaystyle \frac{1}{\mathrm{\Delta }t}}\left({\displaystyle \frac{F_i^{+}(t_n)}{k^i(t_n)}}-{\displaystyle \frac{F_i^{+}(t_{n-1})}{k^i(t_{n-1})}}\right)\end{array}

An improved approximation is gained by substituting Eq. (24) back into Eq. (21),

\begin{array}{}\text{(25)}& Y_i(t_n)\simeq {\displaystyle \frac{F_i^{+}(t_n)}{k^i(t_n)}}-{\displaystyle \frac{1}{k^i(t_n)\mathrm{\Delta }t}}\left({\displaystyle \frac{F_i^{+}(t_n)}{k^i(t_n)}}-{\displaystyle \frac{F_i^{+}(t_{n-1})}{k^i(t_{n-1})}}\right)\end{array}

This approximation is most valid when the right-hand term, which came from the approximated derivative, is small. This can be accomplished if the fluxes into population i are small, but is more likely to occur because the product k\mathrm{\Delta }t is large. If k_i\mathrm{\Delta }t⩾1, then the flux out of species i is greater than the population of species i, and in an ordinary forward Euler timestep, the population of species i would become over-depleted and have a negative value. This would lead to the solution of the system of equations becoming unstable. If however, we choose this over-depletion point as the criteria for using the asymptotic formalism, the over-depletion is prevented and the system will not become unstable. In any given timestep, the asymptotic approximation of Eq. (25) is then employed only for those populations that meet the criteria of k_i\mathrm{\Delta }t⩾1. The remaining species are treated with a different numerical technique, for example, forward Euler, but in principle higher order methods can be used to provide greater accuracy if needed.

3.3.2 Applying explicit methods to thermonuclear supernovae

To examine the utility of the asymptotic method to thermonuclear supernovae, we present a test case drawn from Ref. [97], a thermodynamic profile extracted from a deflagrating white dwarf model from Brown et al. [98]. It begins at an initial temperature of 0.033 \mathrm{GK} with a density of 9.4\times {10}^8\mathrm{g}\cdot {\mathrm{cm}}^{-3}, reaches a peak temperature of 4.7 \mathrm{GK} after 1.57 seconds, by which point the density has dropped to 7.5\times {10}^7\mathrm{g}\cdot {\mathrm{cm}}^{-3}. The calculation concludes at 2.0 seconds, with the temperature of 2.1 \mathrm{GK} and density of 4.5\times {10}^6\mathrm{g}\cdot {\mathrm{cm}}^{-3}. Two nuclear reaction networks are used, a 14 species \alpha-network and more complete network with 150 species. Table 1 compares the final mass fractions for the \alpha-network calculated by the asymptotic solver to those computed by a standard backward Euler solver.

Figure 7 shows a comparison of the timestep behavior of the asymptotic solver to that of the backward Euler solver for this example. The upper panel illustrates the behavior of the \alpha-network. Prior to 1.5 seconds, the timestep behavior of the asymptotic solver is quite comparable to that of the backward Euler, indicating that the asymptotic solver is adept at Carbon burning. However once the temperature exceeds 2 \mathrm{GK} the asymptotic solver begins to trail the backward Euler solver dramatically. At the peak temperature, the asymptotic solver takes more than 10 000 timesteps for each backward Euler timestep. While the asymptotic solver’s timestep grows closer to that of the backward Euler solver as the temperature declines from the peak, it does not achieve equality until 1.8 seconds into the calculation, at which point the temperature has dropped below 3 \mathrm{GK} and substantial nuclear processing has effectively ceased with this reaction network. As a result of the much smaller timesteps near the peak temperature, the total number of timesteps taken by the asymptotic solver is much larger than the backward Euler solver requires. In the case of the alpha network, the calculation using the asymptotic solver took 4.8 million timesteps compared to the implicit solver taking only 636 timesteps. With matrix operations in the alpha network requiring 66% (or less) of the calculation time in the backward Euler case, we can estimate the relative speed of the asymptotic method as (1/(1-66\%))\sim 3 that of the implicit solver, implying that a mature asymptotic solver would compute each timestep at best three times faster than the implicit method. Even with this per timestep advantage, the asymptotic calculation is completely inadequate from a performance perspective, except during carbon burning where it is very competitive. This mirrors previous results from Mott et al. [93] for a similar problem.

Asymptotic calculations with larger networks performed even more poorly at high temperatures. The lower panel of Fig. 7 show the timestep comparisons for the 150 isotope network. As in the \alpha-network case, the early timestep behavior during carbon burning is comparable between asymptotic and backward Euler solvers, but diverge once the temperature exceeds 2 \mathrm{GK}. At peak both of the larger networks require the asymptotic solver to take nearly 1 000 000 timesteps for each timestep taken by the backward Euler. While the disparity between the asymptotic solver and backward Euler solver diminishes with time, it does not do so as sharply in the large network cases as it did in the \alpha-network case. Even at the end of the calculation, 2 seconds from its start with a temperature of 2.1 \mathrm{GK} and density of 4.5\times {10}^6\mathrm{g}\cdot {\mathrm{cm}}^{-3}, the asymptotic solver is taking several hundred timesteps for each backward Euler timestep. This difference occurs because the larger networks contain free neutrons and protons, which continue to react, driving nuclear abundance evolution, after declining temperatures prevent \alpha captures from occurring at a noticeable rate. Hence, freezeout occurs much earlier in the \alpha-network because of its omission of free nucleons. The result of this behavior is that the entire model requires more than 10 million timesteps for the asymptotic solver to complete these large network models, while the backward Euler solver requires less than 2600. Even taking into account the reduced cost per timestep, estimated to be \sim 1/(1-0.8)=5 for the 150-isotope network where the matrix build and solve takes 80% of the implicit method’s cost, the asymptotic method is not competitive except during carbon burning.

It is very revealing that the asymptotic solver is competitive with implicit methods during carbon burning (and also during hydrogen and helium burning, Refs. [92, 97]), but not for hotter burning stages. While it is tempting to ascribe the asymptotic solvers difficulties near the peak temperature of this model simply to the temperature being in excess of 2 \mathrm{GK}, the actual cause is more subtle. For temperatures in excess of 2 \mathrm{GK}, many pairs of forward and inverse reaction reach equilibrium. Figure 8 plots the fraction of equilibrated reaction pairs as well as the fraction of nuclear species for which the asymptotic algorithm is used to update the abundances, those for which k_i\mathrm{\Delta }t⩾1. While a significant number of reactions are defined as “equilibrated” at early times, these fluxes are in fact nearly zero. The rise of the fraction of equilibrated reaction pairs after 1.5 seconds, peaking above 75% for much of the burning, is the establishment of detailed balance and corresponds to the rise in temperature above 2 \mathrm{GK}. At this same time, the majority of species also satisfy the asymptotic condition. Unfortunately, updating the abundance using Eq. (25) performs poorly in cases where the reaction pairs involved in the solution are equilibrated because the approximation of Eq. (22) is not well justified, forcing a smaller timestep. In contrast, the modest level of equilibration shown for hydrogen, helium and carbon burning does not affect species being evolved via the asymptotic approximation [92, 97]. This absence of partial equilibrium allows the asymptotic solver to perform competitively in these stages.

3.3.3 Application of partial equilibria to explicit integration methods

Just as was the case for the implicit network, were the QSE-reduced network turned the numerical difficulty posed by equilibrium into a virtue, equilibrium can also be harnessed to speed up the solution of the explicit methods. In general, the existence of equilibrium abundance relations like Eq. (7) can be used in two ways. First, they can be used to reduce the number of equations to be solved by grouping abundances together. Second, they can be used to remove the equilibrated reactions from the system of equations by providing algebraic ways to update the equilibrated abundances. Since these equilibrated reactions are in general among the fastest in the network, their removal has the effect of reducing the stiffness of the system. In the case of the implicit network solution, the first approach is most efficacious, because the implicit timestep is limited by accuracy not stability and the need to solve a matrix equation places a high premium on the number of equations, and hence the size of the matrix, to be solved. For the explicit network, the cost of additional differential equations is relatively small, but as was shown in the preceding section, instability can place a severe constraint on the timestep. Thus for explicit systems, the second approach to utilizing equilibrium is most useful.

To explore the use of partial equilibrium in explicit reaction network integrators, we begin by re-writing the total fluxes, F in Eq. (17) in terms of their component fluxes, f,

\begin{array}{}\text{(26)}& \dot{Y_i}={\displaystyle \sum _m}{f}_i^{m+}(t)-{\displaystyle \sum _m}{f}_i^{m-}(t)={\displaystyle \sum _m}(f_i^{m+}(t)-f_i^{m-}(t))\end{array}

The second equality reflects a re-grouping of the terms into individual pairs of forward and reverse reactions. Mott [91, 99] introduced the idea of using partial equilibrium to set f_i^{m+}(t)-f_i^{m-}(t)=0 for equilibrated reaction pairs, thereby removing their short timescales from the integration. This idea has recently been extended and improved for astrophysical problems by Guidry et al. [100]. The removal of equilibrated reactions results in a two stage abundance update. First, the equilibrated rates pairs are set to 0 and the isotopic populations are evolved by the system of differential equations including the remaining non-equilibrated reactions. Second, the effects of the equilibrated reactions are included via equilibrium constraint equations derived from Eq. (7), restoring equilibrium. Guidry et al. [100] contains a complete list of these constraint equations, which differ depending on the number of products and reactants involved in a reaction.

One complication in this approach for utilizing partial equilibrium within explicit integration is the evolution of abundances which participate in more than one equilibrated reaction pair. A common example is the free nucleons and \alpha-particles, which are among the reactants or products in the vast majority of reactions at work in an astrophysical nuclear reaction network, but even heavier species may fall into this category. Indeed, under conditions where QSE applies, most heavy nuclei belong to this category since every strong or electromagnetic reaction among QSE group members is equilibrated. Since each equilibrium constraint suggests a (slightly) different equilibrium abundance value for the shared species, some method must be used to find a common abundance that satisfies all equilibrium constraints. The most general approach is to find all of these equilibrium abundances simultaneously by iteration across the entire network. This solution is directly analogous to the iterative determination of the abundances of the focal species from the group abundances in the QSE-reduced network described in Section 3.2. In that case, the iteration, which requires a matrix solution if the Newton–Raphson method is used, is over a small set of group abundance equations and is therefore of small computational cost. In the case of partial equilibrium applied to explicit methods, the iteration is over a larger number of equilibrium constraints for individual reaction pairs, making the computational cost sizable. This is particularly distasteful for a method whose chief virtue is the avoidance of a matrix solution during the network integration in time. Fortunately, Guidry et al. [100] have found that a simpler approach, obtaining each equilibrium abundance by averaging over the values computed by the different relevant equilibrium constraints, is sufficiently accurate, at least for the small networks tested thus far.

Figure 9 shows an example of one such test, integration of an \alpha-network over a simple thermodynamic profile (shown in the lower panel of Fig. 9) similar to the results of the passage of a thermonuclear flame through the white dwarf. As in the case documented in Fig. 7, the timestep used by the asymptotic solver trails that of the implicit solver by almost 5 orders of magnitude. However, the use of partial equilibrium within the asymptotic solver (labeled PE in Fig. 9) improves the timestep by more than a factor of 1000, making the explicit solver with partial equilibrium much more competitive with the implicit solver in terms of time. Further refinement is needed in this new solver (for example, the PE solver does not employ the same accuracy limit on the timestep, likely the reason the PE solver’s time steps exceed those of the implicit solver in the low temperature portion of Fig. 9). Nonetheless it is clear that this explicit integrator, including partial equilibrium where appropriate, is much more effective and efficient for conditions encountered during thermonuclear supernovae than previous results found.

4 Conclusion

Our understanding of the mechanism that produces Type Ia supernova has grown over time with models of increasing sophistication. Initially, spherically symmetric (1D) models with small nuclear networks revealed the need for the flame to begin as a deflagration. Increasingly sophisticated 1D models, using larger nuclear reaction networks, now provide excellent reproduction of the observed spectra and light curves, albeit using proscribed flame speeds and mixing. These models predictions of the nucleosynthesis have revealed the contribution that thermonuclear supernovae make to galactic chemical evolution, especially for the iron-peak species.

The transition over the past decade to two dimensional, axisymmetric models and now to three dimensional models has removed the need for prescriptions on flame speed and mixing, but at a cost of simplified thermonuclear evolution. This has been addressed in part by post-processes calculations relying on tracer particles, for example such studies have revealed the potential for the p-process to occur in these supernovae. However, the next step is to perform multi-dimensional simulations with larger nuclear reaction networks, matching the capabilities of 1D models. This next step however comes at a considerable computational cost, as 10 times more nuclear species requires at least 100 times more computation. It seems clear that quasi-equibrium will be critical to managing this computational cost, either in the context of the fully implicit solutions to the network equations (the QSE-reduced network) or as an component in explicit solvers. With the savings that these methods afford, the next step in simulating thermonuclear supernovae moves closer to hand.

References

[1] W. Hillebrandt and J. C. Niemeyer, Ann. Rev. Astron. Astrophys., 2000, 38: 191, arXiv: astro-ph/0006305

[2] K. Nomoto, K. Iwamoto, and N. Kishimoto, Science, 1997, 276(5317): 1378, arXiv: astro-ph/9706007

[3] D. Branch, Annu. Rev. Astron. Astrophys., 1998, 36(1): 17, arXiv: astro-ph/9801065

[4] P. Höflich and A. Khokhlov, Astrophys. J., 1996, 457: 500

[5] P. Nugent, E. Baron, D. Branch, A. Fisher, and P. H. Hauschildt, Astrophys. J., 1997, 485(2): 812, arXiv: astro-ph/9612044

[6] M. Fink, W. Hillebrandt, and F. K. Röpke, Astron. Astrophys., 2007, 476(3): 1133, arXiv: 0710.5486

[7] P. A. Mazzali, F. K. Röpke, S. Benetti, and W. Hillebrandt, Science, 2007, 315(5813): 825, arXiv: astro-ph/0702351

[8] J. Whelan and I. Iben, Jr., Astrophys. J., 1973, 186: 1007

[9] K. Nomoto, Astrophys. J., 1982, 253: 798

[10] I. Iben, Jr., and A. V. Tutukov, Astrophys. J. Suppl., 1984, 54: 335

[11] R. F. Webbink, Astrophys. J., 1984, 277: 355

[12] F. K. Röpke, M. Kromer, I. R. Seitenzahl, R. Pakmor, et al., Astrophys. J., 2012, 750(1): L19, arXiv: 1203.4839/astro-ph.SR

[13] D. Arnett and E. Livne, Astrophys. J., 1994, 427: 330

[14] W. D. Arnett, Astrophys. J., 1969, 157: 1369

[15] K. Nomoto, F. K. Thielemann, and K. Yokoi, Astrophys. J., 1984, 286: 644

[16] A. M. Khokhlov, Astron. Astrophys., 1991, 246: 383

[17] H. Yamaoka, K. Nomoto, T. Shigeyama, and F. K. Thielemann, Astrophys. J., 1992, 393: L55

[18] S. E. Woosley and T. A. Weaver, in: Supernovae, NATO Advanced Science Institutes (ASI) Series C: Mathematical and Physical Sciences, edited by S. A. Bludman, R. Mochkovitch, and J. Zinn-Justin, Amsterdam: Elsevier, 1994: 63

[19] K. Iwamoto, F. Brachwitz, K. Nomoto, N. Kishimoto, W. Hix, and F. K. Thielemann, Astrophys. J. Suppl., 1999, 125(2): 439

[20] F. Brachwitz, D. Dean, W. Hix, K. Iwamoto, K. Langanke, G. Martinez-Pinedo, K. Nomoto, M. R. Strayer, F. K. Thielemann, and H. Umeda, Astrophys. J., 2000, 536(2): 934

[21] E. Livne, Astrophys. J., 1993, 412: 634

[22] M. Reinecke, W. Hillebrandt, J. C. Niemeyer, R. Klein, and A. Gröbl, Astron. Astrophys., 1999, 347: 724, arXiv: astro-ph/9812119

[23] A. M. Lisewski, W. Hillebrandt, and S. E. Woosley, Astrophys. J., 2000, 538(2): 831, arXiv: astro-ph/9910056

[24] F. K. Röpke and W. Schmidt, in: Interdisciplinary Aspects of Turbulence, Lecture Notes in Physics, Vol. 756, edited by W. Hillebrandt and F. Kupka, Berlin: Springer-Verlag, 2009: 255

[25] M. Reinecke, W. Hillebrandt, and J. C. Niemeyer, Astron. Astrophys., 2002, 391(3): 1167, arXiv: astro-ph/0206459

[26] V. N. Gamezo, A. M. Khokhlov, and E. S. Oran, Astrophys. J., 2005, 623(1): 337, arXiv: astro-ph/0409598

[27] I. R. Seitenzahl, F. Ciaraldi-Schoolmann, F. K. Röpke, M. Fink, W. Hillebrandt, M. Kromer, R. Pakmor, A. J. Ruiter, S. A. Sim, and S. Taubenberger, Mon. Not. R. Astron. Soc., 2013, 429(2): 1156, arXiv: 1211.3015/astro-ph.SR

[28] J. R. Boisseau, J. C. Wheeler, E. S. Oran, and A. M. Khokhlov, Astrophys. J., 1996, 471(2): L99

[29] F. X. Timmes, M. Zingale, K. Olson, B. Fryxell, P. Ricker, A. C. Calder, L. J. Dursi, H. Tufo, P. MacNeice, J. W. Truran, and R. Rosner, Astrophys. J., 2000, 543(2): 938

[30] A. C. Calder, D. M. Townsley, I. R. Seitenzahl, F. Peng, O. E. B. Messer, N. Vladimirova, E. F. Brown, J. W. Truran, and D. Q. Lamb, Astrophys. J., 2007, 656(1): 313, arXiv: astro-ph/0611009

[31] J. C. Niemeyer, Astrophys. J., 1999, 523(1): L57, arXiv:astro-ph/9906142

[32] F. X. Timmes, S. E. Woosley, and T. A. Weaver, Astrophys. J. Suppl., 1995, 98: 617, arXiv: astro-ph/9411003

[33] F. K. Röpke, M. Gieseler, M. Reinecke, C. Travaglio, and W. Hillebrandt, Astron. Astrophys., 2006, 453(1): 203, arXiv: astro-ph/0506107

[34] C. Travaglio, F. K. Röpke, R. Gallino, and W. Hillebrandt, Astrophys. J., 2011, 739(2): 93, arXiv: 1106.0582/astro-ph.SR

[35] S. Nagataki, M. Hashimoto, K. Sato, and S. Yamada, Astrophys. J., 1997, 486(2): 1026

[36] C. Travaglio, W. Hillebrandt, M. Reinecke, and F. K. Thielemann, Astron. Astrophys., 2004, 425(3): 1029, arXiv: astro-ph/0406281

[37] I. R. Seitenzahl, F. K. Röpke, M. Fink, and R. Pakmor, Mon. Not. R. Astron. Soc., 2010, 407(4): 2297, arXiv: 1005.5071/astro-ph.SR

[38] F. X. Timmes, E. F. Brown, and J. W. Truran, Astrophys. J., 2003, 590(2): L83, arXiv: astro-ph/0305114

[39] C. Travaglio, W. Hillebrandt, and M. Reinecke, Astron. Astrophys., 2005, 443(3): 1007, arXiv: astro-ph/0507510

[40] C. Travaglio, R. Gallino, and F. K. Röpke, in: Proceedings of Nuclei in the Cosmos XII, SISSA Proceedings of Science, 2012: 45

[41] I. Domínguez, P. Höflich, and O. Straniero, Astrophys. J., 2001, 557(1): 279, arXiv: astro-ph/0104257

[42] I. Iben, Jr., Astrophys. J., 1981, 243: 987

[43] I. Iben, Jr., and A. V. Tutukov, Astrophys. J., 1991, 370: 615

[44] W. M. Howard and B. S. Meyer, in: Nuclei in the Cosmos 2, edited by F. Kaeppeler and K. Wisshak, Bristol: IOP Publishing, 1993: 575

[45] M. Kusakabe, N. Iwamoto, and K. Nomoto, Astrophys. J., 2011, 726(1): 25, arXiv: 1001.0145/astro-ph.SR

[46] W. Li, J. Leaman, R. Chornock, A. V. Filippenko, D. Poznanski, M. Ganeshalingam, X. Wang, M. Modjaz, S. Jha, R. J. Foley, and N. Smith, Mon. Not. R. Astron. Soc., 2011, 412(3): 1441, arXiv: 1006.4612/astro-ph.SR

[47] I. Dillmann, T. Rauscher, M. Heil, F. Käppeler, W. Rapp, and F. K. Thielemann, J. Phys. G, 2008, 35(1): 014029, arXiv: 0805.4756

[48] Z. Nemeth, F. Kaeppeler, C. Theis, T. Belgya, and S. W. Yates, Astrophys. J., 1994, 426: 357

[49] S. E. Woosley, D. H. Hartmann, R. D. Hoffman, and W. C. Haxton, Astrophys. J., 1990, 356: 272

[50] C. Arlandini, F. Käppeler, K. Wisshak, R. Gallino, M. Lugaro, M. Busso, and O. Straniero, Astrophys. J., 1999, 525(2): 886, arXiv: astro-ph/9906266

[51] F. Käppeler, R. Gallino, S. Bisterzo, and W. Aoki, Rev. Mod. Phys., 2011, 83(1): 157

[52] K. Takahashi and K. Yokoi, At. Data Nuc. Data Tab., 1987, 36(3): 375

[53] P. Mohr, F. Käppeler, and R. Gallino, Phys. Rev. C, 2007, 75(1): 012802, arXiv: astro-ph/0612427

[54] K. Maeda, F. K. Röpke, M. Fink, W. Hillebrandt, C. Travaglio, and F. K. Thielemann, Astrophys. J., 2010, 712(1): 624, arXiv: 1002.2153/astro-ph.SR

[55] M. Stehle, P. A. Mazzali, S. Benetti, and W. Hillebrandt, Mon. Not. R. Astron. Soc., 2005, 360(4): 1231, arXiv: astro-ph/0409342

[56] C. Kozma, C. Fransson, W. Hillebrandt, C. Travaglio, J. Sollerman, M. Reinecke, F. K. Röpke, and J. Spyromilio, Astron. Astrophys., 2005, 437(3): 983, arXiv: astro-ph/0504317

[57] S. I. Blinnikov, F. K. Röpke, E. I. Sorokina, M. Gieseler, M. Reinecke, C. Travaglio, W. Hillebrandt, and M. Stritzinger, Astron. Astrophys., 2006, 453(1): 229, arXiv: astro-ph/0603036

[58] M. Kromer, M. Fink, V. Stanishev, S. Taubenberger, F. Ciaraldi-Schoolman, R. Pakmor, F. K. Röpke, A. J. Ruiter, I. R. Seitenzahl, S. A. Sim, G. Blanc, N. Elias-Rosa, and W. Hillebrandt, Mon. Not. R. Astron. Soc., 2013, 429(3): 2287, arXiv: 1210.5243 [astro-ph.HE]

[59] E. Baron, P. Höflich, K. Krisciunas, I. Dominguez, A. M. Khokhlov, M. M. Phillips, N. Suntzeff, and L. Wang, Astrophys. J., 2012, 753(2): 105, arXiv: 1205.0814/astroph.SR

[60] P. Höflich, J. C. Wheeler, and F. K. Thielemann, Astrophys. J., 1998, 495(2): 617

[61] S. E. Woosley, in: 16th Saas Fee Advanced Course, Nucleosynthesis and Chemical Evolution, edited by B. Houck, A. Maeder, and G. Meynet, Sauverny: Geneva Obs., 1986: 1

[62] F. K. Thielemann, K. Nomoto, and M. Hashimoto, in: Supernovae, Les Houches, Session LIV, edited by S. Bludman, R. Mochkovitch, and J. Zinn-Justin, Amsterdam: Elsevier, 1994: 629

[63] W. D. Arnett, Supernovae and nucleosynthesis: An investigation of the history of matter, from the big bang to the present, Princeton: Princeton University Press, 1996

[64] W. R. Hix and F. Thielemann, J. Comput. Appl. Math., 1999, 109(1–2): 321, arXiv: astro-ph/9906478

[65] W. R. Hix and B. S. Meyer, Nucl. Phys. A, 2006, 777: 188, arXiv: astro-ph/0509698

[66] C. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Englewood Cliffs, NJ: Prentice-Hall, 1971

[67] J. Lambert, in: Computational Techniques for Ordinary Differential Equations, edited by I. Gladwell and D. Sayars, New York: Academic Press, 1980: 19

[68] E. Oran and J. Boris, Numerical Simulation of Reactive Flow, New York: Elsevier, 1987

[69] W. D. Arnett and J. W. Truran, Astrophys. J., 1969, 157: 339

[70] R. V. Wagoner, Astrophys. J., 1973, 179: 343

[71] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes, 2nd Ed., Cambridge: Cambridge University Press, 1992

[72] F. X. Timmes, Astrophys. J. Suppl., 1999, 124(1): 241

[73] E. Müller, Astron. Astrophys., 1986, 162: 103

[74] E. Müller, in: Saas-Fee Advanced Course 27: Computational Methods for Astrophysical Fluid Flow, 1998: 343

[75] D. D. Clayton, Principles of Stellar Evolution and Nucleosynthesis, Chicago: University of Chicago Press, 1983

[76] C. Rolfs and W. Rodney, Cauldrons in the Cosmos, Chicago: University of Chicago Press, 1988

[77] C. Iliadis, Nuclear Physics of Stars, by Christian Il-iadis, Wenheim: Wiley-VCH, 2007

[78] F. E. Clifford and R. J. Tayler, Mem. RAS, 1965, 69: 21

[79] D. Hartmann, S. E. Woosley, and M. F. El Eid, Astrophys. J., 1985, 297: 837

[80] I. R. Seitenzahl, D. M. Townsley, F. Peng, and J. W. Truran, At. Data Nuc. Data Tab., 2009, 95(1): 96

[81] B. S. Meyer, T. D. Krishnan, and D. D. Clayton, Astrophys. J., 1998, 498(2): 808

[82] W. R. Hix and F. K. Thielemann, Astrophys. J., 1999, 511(2): 862

[83] D. Bodansky, D. D. Clayton, and W. A. Fowler, Astrophys. J. Suppl., 1968, 16: 299

[84] S. E. Woosley, W. D. Arnett, and D. D. Clayton, Astrophys. J. Suppl., 1973, 26: 231

[85] A. M. Khokhlov, Sov. Astron. Lett., 1981, 7: 410

[86] F. X. Timmes, R. D. Hoffman, and S. E. Woosley, Astrophys. J. Suppl., 2000, 129(1): 377

[87] W. R. Hix, A. M. Khokhlov, J. C. Wheeler, and F. K. Thielemann, Astrophys. J., 1998, 503(1): 332

[88] W. R. Hix, S. T. Parete-Koon, C. Freiburghaus, and F. K. Thielemann, Astrophys. J., 2007, 667(1): 476

[89] S. T. Parete-Koon, W. R. Hix, and F. K. Thielemann, in: Proceedings of Nuclei in the Cosmos X, edited by H. Schatz et al., SISSA Proceedings of Science, 2008: 157

[90] S. T. Parete-Koon, W. R. Hix, and F. K. Thielemann, Astrophys. J., 2013 (in preparation)

[91] D. Mott, New Quasi-Steady-State and Partial-Equilibrium Methods for Integrating Chemically Reacting Systems, Ph.D. thesis, University of Michigan, Ann Arbor, 1999

[92] M. W. Guidry, R. Budiardja, E. Feger, J. J. Billings, W. R. Hix, O. E. B. Messer, K. J. Roche, E. McMahon, and M. He, Computational Science and Discovery, 2013, 6(1): 015001, arXiv: 1112.4716 [astro-ph.SR]

[93] D. Mott, E. Oran, and B. van Leer, J. Comput. Phys., 2000, 164(2): 407

[94] M. W. Guidry and J. A. Harris, Comp. Sci. Disc., 2013, 6(1): 015002, arXiv: 1112.4750 [astro-ph.SR]

[95] T. R. Young and J. P. Boris, J. Phys. Chem., 1977, 81(25): 2424

[96] E. Oran, and J. Boris, Numerical Simulation of Reactive Flow, 2nd Ed., Cambridge: Cambridge University Press, 2001

[97] E. D. Feger, M. W. Guidry, and W. R. Hix, Astrophys. J. Suppl., 2013 (submitted)

[98] E. F. Brown, A. C. Calder, T. Plewa, P. M. Ricker, K. Robinson, and J. B. Gallagher, Nucl. Phys. A, 2005, 758: 451, arXiv: astro-ph/0505417

[99] D. Mott, E. Oran, and B. van Leer, in: AIAA-2003-667, 41st Aerospace Sciences Meeting and Exhibit, Reno, Nevada, Jan. 6–9, 2003: 389

[100] M. W. Guidry, J. J. Billings, and W. R. Hix, Comp. Sci. Disc., 2013, 6(1): 015003, arXiv: 1112.4738 [astro-ph.SR]

Rights & permissions

Higher Education Press and Springer-Verlag Berlin Heidelberg 2013