SN Ia Nucleosynthesis

A Tracer Method for Computing Type Ia Supernova Yields:
Burning Model Calibration,
Reconstruction of Thickened Flames,
and Verification for Planar Detonations

Dean M. Townsley11affiliation: Department of Physics & Astronomy, University of Alabama, Tuscaloosa, AL; , Broxton J. Miles11affiliation: Department of Physics & Astronomy, University of Alabama, Tuscaloosa, AL; , F. X. Timmes22affiliation: School of Earth and Space Exploration, Arizona State University, Tempe, AZ 33affiliation: The Joint Institute for Nuclear Astrophysics , Alan C. Calder44affiliation: Department of Physics & Astronomy, Stony Brook University, Stony Brook, NY 55affiliation: Institute for Advanced Computational Sciences, Stony Brook University, Stony Brook, NY , and Edward F. Brown66affiliation: Department of Physics & Astronomy, Michigan State University, East Lansing, MI 33affiliation: The Joint Institute for Nuclear Astrophysics

We refine our previously introduced parameterized model for explosive carbon-oxygen fusion during thermonuclear supernovae (SN Ia) by adding corrections to post-processing of recorded Lagrangian fluid element histories to obtain more accurate isotopic yields. Deflagration and detonation products are verified for propagation in a uniform density medium. A new method is introduced for reconstructing the temperature-density history within the artificially thick model deflagration front. We obtain better than 5% consistency between the electron capture computed by the burning model and yields from post-processing. For detonations, we compare to a benchmark calculation of the structure of driven steady-state planar detonations performed with a large nuclear reaction network and error-controlled integration. We verify that, for steady-state planar detonations down to a density of  g cm, our post processing matches the major abundances in the benchmark solution typically to better than 10% for times greater than 0.01 s after the shock front passage. As a test case to demonstrate the method, presented here with post-processing for the first time, we perform a two dimensional simulation of a SN Ia in the Chandrasekhar-mass deflagration-detonation transition (DDT) scenario. We find that reconstruction of deflagration tracks leads to slightly more complete silicon burning than without reconstruction. The resulting abundance structure of the ejecta is consistent with inferences from spectroscopic studies of observed SNe Ia. We confirm the absence of a central region of stable Fe-group material for the multi-dimensional DDT scenario. Detailed isotopic yields are tabulated and only change modestly when using deflagration reconstruction.

Subject headings:
supernovae: general – nuclear reactions, nucleosynthesis, abundances – methods: numerical

1. Introduction

Type Ia Supernovae (SNe Ia) are at once a pillar of modern cosmology and one of the persistent puzzles of stellar physics. These bright stellar transients are characterized by strong P Cygni features in Si and a lack of hydrogen or helium in their spectra. It has generally been accepted that these events follow from the thermonuclear incineration of a white dwarf (WD) star producing between 0.3 and 0.9 of radioactive Ni, the decay of which powers the light curve (see Filippenko, 1997; Hillebrandt & Niemeyer, 2000; Röpke, 2006; Calder et al., 2013, and references therein). The light curves of SNe Ia have the property that the brightness of an event is correlated with its duration (Phillips, 1993). This relation is the basis for light curve calibration that allows use of these events as distance indicators for cosmological studies (see Conley et al. 2011 for a contemporary example). However, their exact stellar origin remains unclear, even in the face of extensive observational and theoretical study. Recent early-time observations of the nearby SN Ia 2011fe are challenging for a variety of common progenitor scenarios, both single and double degenerate (Nugent et al., 2011; Li et al., 2011; Bloom et al., 2012; Chomiuk et al., 2012). Fitting of a wide variety of light curves with a simplified ejecta model appears to require ejecta masses both at and below the Chandrasekhar mass (Scalzo et al., 2014), however, indicating a variety of progenitors may be present.

In thermonuclear supernovae, explosive nuclear combustion of a degenerate carbon oxygen mixture proceeds in one or both of the deflagration and detonation combustion modes. In a deflagration, or flame, the reaction front propagates by thermal conduction (Timmes & Woosley, 1992; Chamulak et al., 2007), and is therefore subsonic. In a detonation, the reaction front propagates via a shock that moves supersonically with respect to the fuel (Khokhlov, 1989; Sharpe, 1999). These two combustion modes have been used to construct a variety of possible explosion scenarios, either in combination, as in the deflagration-detonation transition (DDT) scenario (Khokhlov, 1991), an example of which is presented in this work, or singly as in the double-detonation model (Livne & Arnett, 1995; Fink et al., 2010) or the pure-deflagration model (Fink et al., 2014).

Figure 1.— Length and time scales of burning processes as a function of fuel density for a mixture of C:O:Ne in the proportion 50:48:2 compared to typical simulation resolution (4km, solid black lines) and scales of the star (dashed black lines). Stellar scales are taken to be the radius of the initial star and the dynamical time. The top panel shows the size scale of various reaction front features while the bottom shows the self-crossing time of these features at the propagation speed of the reaction front. The behavior of the C flame is shown by the solid blue curve that extends to high density (Chamulak et al., 2007). Scales for planar steady-state detonations are shown in red. We show the distance from the shock to three points in the detonation structure: the peak of Si abundance (dashed), which is also the end of O consumption; the sonic point (dash-dot), also called the pathological point, which is also the size of the detonation driving region; and the attainment of the fully burned NSE state (solid), which is the completion of consumption of Si. The detonation driving region transitions from being just resolved to being larger than the radius of the star between densities of about  g cm and  g cm.

A major challenge in simulations of SNe Ia is capturing these burning processes with confidence and accuracy. The carbon-oxygen reaction fronts transition from being unresolved by many orders of magnitude, to being partially resolved, to finally being larger than the time and length scales of the star. Figure 1 shows length and time scales for detonations (red) and deflagrations (blue, Chamulak et al. 2007) at various densities. For the stellar scales we take the initial WD radius,  cm, and the dynamical time  s where is the WD mass. A representative simulation resolution of 4 km is shown, along with the corresponding timestep of approximately the sound crossing time of a cell. This resolution was found by Townsley et al. (2009) to be sufficient to give convergence in 1D with the thickened flame reaction front. As a result our multi-dimensional (multi-D) simulations are commonly performed between 4 and 1 km to study resolution dependence in a regime in which convergence is demonstrated in 1D. Several different stages within a steady-state planar detonation front are indicated, with distances measured from the shock that initiates the reactions and propagates the front. The shortest length scale shown (dashed line) is that when the Si abundance peaks in time, which also corresponds to the completion of O consumption. The next length scale (dot-dashed line) is the size of the detonation driving region, which is the distance to the sonic point. Finally the solid line shows the distance to completion of burning, reaching the Fe-group element (IGE) dominated nuclear statistical equilibrium (NSE) state.

The total yields of the explosion are determined by how and when the reaction fronts stop propagating as well as what portion of the burning occurs within the reaction front as opposed to what occurs after the reaction front itself has passed. The latter can then be influenced by the expansion of the star. The fairly thin range of densities, , being density in units of  g cm, in which the detonation driving region transitions from being unresolved to being larger than the radius of the star is a manifestation of the difficulty of capturing the reaction dynamics appropriately. As the driving time and length scales get large, the detonation may not be able to attain the planar steady-state structure. Curvature of the front on scales comparable to the driving length, which will occur due to the structure of the star, reduce the detonation speed and the completeness of the burning (Sharpe, 2001; Dunkley et al., 2013). The long reaction times also mean that an ignited detonation may not reach steady state before the star expands (Townsley et al., 2012).

Many recent results on multi-D simulations of SNe Ia have computed nucleosynthetic yields by post-processing the density and temperature recorded by a Lagrangian fluid history during the simulation (e.g. Travaglio et al., 2004). A large nuclear reaction network is used to integrate a set of species subject to this , history. The burning model used in the simulation is therefore critical, as it determines these histories. Recent multi-D work (Maeda et al., 2010; Seitenzahl et al., 2010; Ciaraldi-Schoolmann et al., 2013; Seitenzahl et al., 2013) has used the method described in the appendix of Fink et al. (2010) to set the energetics of the burning model used in the hydrodynamics. In this technique, the results of the post-processing are used to revise the output of the model of burning and the process is iterated until the yields no longer change.

In this work we pursue a different route toward construction of our burning model and post-processing methods which, instead of an iterative bootstrap, is based on comparison to separate resolved calculations of the deflagration and detonation modes. The burning model and post-processing method are then constructed with the goal that the post-processed results reproduce the results of resolved calculations of the steady-state structure of the reaction front even though the actual reaction front is unresolved. The resolved calculations to which we want to compare are standard methods (e.g. Fickett & Davis, 1979) for the computation of reaction front structure that can be performed with fairly complete nuclear reaction networks and using error-controlled time integration methods to eliminate most computational uncertainty. Here we succeed in matching steady-state yields for detonations at high densities and in planar geometry. Further development of benchmarks and methods for lower densities and other geometries in future work will enable confident higher-accuracy yields for an even larger fraction of the ejecta.

The burning model presented here is the successor to that initially presented by Calder et al. (2007) and Townsley et al. (2007), with tabulations presented by Seitenzahl et al. (2009b), that has been used in a number of studies using large multi-D simulations of SN Ia (Jordan et al., 2008; Meakin et al., 2009; Jordan et al., 2012b, a; Kim et al., 2013; Long et al., 2014). The capability to treat neutron-enriched fuel was added by Townsley et al. (2009) in order to study how neutron enrichment in the progenitor might influence the explosion. The model presented here includes a change in dynamics to better match iron-group production in detonations and extends the treatment of initial composition to spatially non-uniform abundances, allowing more realistic WD progenitors. This has been used in work exploring systematic effects of progenitor WD composition and central density in the DDT scenario, (Krueger et al., 2010; Jackson et al., 2010; Krueger et al., 2012), as well as a study of the turbulence-flame interaction during the deflagration phase (Jackson et al., 2014), and consideration of hybrid C-O-Ne progenitor WDs (Willcox et al., 2016). Those studies, however, did not proceed to nucleosynthetic post-processing, which is discussed in detail here for the first time for our burning model. The first work utilizing the post-processing for astrophysical study is an investigation of spectral indicators of progenitor metallicity (Miles et al., 2015).

We present below the structure of our burning model and post-processing methods, along with particular assumptions currently in use in our SN Ia simulations, as well as tests performed so far comparing to calculations of steady-state deflagrations and detonations. Our burning model is based on tabulation of physical quantities and fits of parameters based on resolved steady-state calculations. To improve accuracy in post-processing, we explore supplementing the Lagrangian - history recorded during the hydrodynamic simulation with a reconstruction of unresolved processes based on conditions near the reaction front when the fluid element is burned.

In section 2 we present the structure of our model for carbon-oxygen burning including the basic variables and the form of their dynamics. Following this, we discuss our post processing treatment for tracks (fluid elements) burned by the deflagration front in section 3. This section is fairly brief since the application of a burning model like that presented here to deflagrations was a major topic of previous work detailed by Calder et al. (2007) and Townsley et al. (2007). Detonations are discussed in two sections. Section 4 develops the error-controlled computation of steady-state detonation structure that we use as a benchmark, calibrates the timescales in the burning model dynamics based on this, and compares the resulting dynamics of the burning model in hydrodynamic tests to the benchmark calculations. The full method including track post-processing is then outlined and tested in section 5. Finally in section 6 we detail how results from full-star simulations are post-processed, and in section 7 we show the results of applying these methods to compute the yields of a 2D simulation of the DDT model of SN Ia, including a consideration of what we can infer about current uncertainties. We summarize conclusions in section 8.

2. Improved Parameterized Model for Explosive Carbon-Oxygen Fusion

We present here our current parameterized model for the thermonuclear burning of carbon and oxygen fuel. The model is intended to capture the dynamics of burning for densities relevant to SNe Ia for either the deflagration or detonation mode of combustion. Conversion of protons to neutrons (neutronization or deleptonization) is included. The initial abundance of carbon and neutron-rich elements (e.g., Ne) are allowed to vary with position in the WD. The model is constructed to use a small number of scalars to track the reaction state and products in order to improve computational efficiency. Accurate final-state energy release and electron capture rates are obtained by tabulation. Abundances of intermediate burning stages are approximated and the formalism can be further refined by adjusting these if necessary.

The process of explosive carbon-oxygen fusion can be roughly divided into 3 stages – C consumption, O consumption, and conversion of Si-group to Fe-group material (Khokhlov, 1989, 2000; Calder et al., 2007). The main processes involved in each of these stages are: C destroyed to produce additional O, Si, Ne, and Mg; O destroyed to produce Si, S, Ar, and Ca, generally in nuclear quasi statistical equilibrium (QSE, sometimes called NSQE); particles liberated by photodisintegration are then captured until this material is converted into Fe-group, eventually reaching full nuclear statistical equilibrium (NSE). Due to differences in the rates of the nuclear processes involved, at densities of interest these stages are well-separated, in logarithmic time, and sequential. This structure makes it possibly to greatly simplify the complex reaction state and dynamics to the behavior of a model containing just a few reaction progress variables.

Individual cells are allowed to contain both unburned and fully burned material in order to allow modeling of reaction fronts that are much thinner than the grid scale. This conceptual structure is shown in Figure 2, where curves are shown to represent two distinct processes in the overall burn which occur on different timescales. In this example we will use O consumption as the shorter-timescale process and Si consumption as the longer-timescale one. The curves indicate the contour on which the O abundance reaches half of its value in the fuel (solid) and where the Ni abundance reaches half its final value (dashed). Intermediate Ni abundances are represented by dotted lines, which would not be distinct at high densities (the separation between stages is exaggerated at high density). At high densities all reaction stages are localized on scales much smaller than the grid, as indicated by the reaction length scales shown in Figure 1, leading to cells which are volumetrically divided into fuel and ash. At lower densities, some burning stages become resolvable, while others remain thin compared to the grid. For resolvable stages, the actual abundance structure more closely resembles a spatial interpolation of the coarse grid values. This is demonstrated in the lower panel of Figure 2. A structure like this is present for both detonation and deflagration combustion modes, though in the turbulent deflagration phase the thin reaction front structure can me much more irregular than shown in this diagram.

Figure 2.— Diagram of structure of thin, multistage reaction fronts embedded in a coarse computational grid of control volume cells. At high densities the reaction front is completely localized – spatially thin – such that a mixed computational cell contains physically separated regions of fuel and ash material. At lower densities, some reaction stages remain localized, while others extend over multiple cells, so that the reactant and product abundances are smoothly varying on subgrid scales.

Our burning model is currently implemented in the Flash code, an adaptive-mesh reactive hydrodynamics code with additional physics for astrophysical applications developed at the University of Chicago (Fryxell et al., 2000; Dubey et al., 2009). The model is readily adaptable for use with other similar reactive hydrodynamics software and the source code is available as add-on Units for the Flash, distributed separately to allow a more liberal license.

2.1. Definition of Stages and Relation to Fluid Properties

The first step in abstraction of the fusion processes is defining the relation of our progress variables to the actual physical properties of the fluid. The transformations taking place via nuclear reactions act most fundamentally on the abundances in the fluid. Since we will reduce the burning processes to just a few stages, we must define first how these stages are related to the actual abundances. After this we will proceed to develop reaction kinetics that will reproduce the effects that the nuclear reactions have on the actual abundances and how that manifests in the corresponding abstracted stages.

Our basic physics will be phrased in terms of baryon fractions, denoted by the symbol . These are the fraction of the total number of baryons which are in the form of the nuclide indicated by the label . This is very similar to the traditional definition of mass fractions, but avoids the ambiguity that rest mass is not conserved as nuclear reactions take place due to energy release. Since Baryon number is a conserved quantity, in the absence of sources the baryon number density, , satisfies the continuity equation:


where is the fluid velocity. For reasons of convenience in a non-relativistic fluid code, we will make the definition


where is the atomic mass unit and is Avogadro’s number. Here is used to denote ”is numerically equivalent to in cgs units.” Our intention is to place distinction between mass and (binding) energy in the gravitational treatment; we could calculate the mass-energy density if necessary, but computation of gravity in our simulation will just use to approximate it. The baryon fractions then identify directly the number of baryons in species and therefore, in the absence of reactions, also follow a conservation equation of the form


Where may be or one of the progress variable defined below that will be constructed as linear combinations of the . Taken together Equation (3) and Equation (2) mean that any linear combination of baryon fractions can be treated as “mass scalars” by the advection infrastructure in conservative fluid dynamics software (e.g. Flash; Fryxell et al. 2000; Dubey et al. 2009).

In order to start from quantities that satisfy Equation (1), for purposes of tracking 3 stages of burning (that is three transitions), we conceive of having four sets of all nuclides, each of which represents a certain ”type” of material:

These denote, respectively, the baryon fractions of individual species comprising fuel, (intermediate) ash (product of carbon consumption), a quasi-equilibrium (QSE) group, and a terminal (NSE) group. These stages and the various symbols used here are laid out in the diagram in Figure 3. This means that any given baryon has two labels: the type of nucleus in which it resides (e.g. silicon), and whether we call that material part of, for example, the ash or the QSE material. It is convenient to define 4 “superspecies” by

Dominant species:material in stage:Amount ofStage:Abundances:Properties:Progress variable:
Figure 3.— Diagram of burning stages for C+O burning and associated symbols used here. The progress of burning of a fluid element is from left to right, with the horizontal axis indicating distance with respect to the foremost part of the reaction front or time as a fluid element passes through the reaction front. See section 4.3 for comparison to an actual reaction front structure.

Since each of the burning stages follows in sequence from the earlier ones – a feature unlike general nuclear species in a reaction network – it is convenient to define progress variables such that




By virtue of the property and thus , we see that . Also, since the are simply linear combinations of the they also satisfy continuity, Equation (1), in the absence of sources.

We define a set of specific abundances:


so that . This is useful because we will, at times, want to specify by specifying . This subtle distinction was left unaddressed in our previous revisions of this burning model (Calder et al., 2007; Townsley et al., 2007). Note that since the are quotients of the , they are no longer linear combinations. Nonlinear terms are any that contain products or quotients of fields that are position-dependent. Linear combinations are required in order for the numerical scheme to be explicitly conservative. While any general algebraic combination, including a nonlinear one like a product or power, of quantities satisfying Equation (1) still satisfies Equation (1), once the fields , , and are discretized into values averaged over control volumes, i.e., mesh cells, conservation of ’s no longer implies conservation of ’s due to nonlinearities in the advection scheme. This results from the property that the average of a quotient is not the quotient of the averages. Since we will not compute Equation (1) for both the ’s and ’s, we must make a choice of which will satisfy explicit conservation to numerical accuracy, as performed by a conservative advection scheme like that in Flash. Since overall energy release is important, we choose to compute conservative evolution for quantities that are linear combinations of the ’s.

In order to evaluate fluid properties and follow nuclear energy release into the fluid, we must be able to obtain several bulk quantities. We will express these in units per baryon, such that obtaining units per cm is trivial using the Baryon density . The two fluid quantities necessary are:


where, as customary, is the number of protons and is the number of protons plus neutrons in nuclide . Here we have assumed charge neutrality between the number of protons and the number of non-thermal electrons, and defined to include only the net non-thermal electrons. Some - pairs are created thermally at high temperatures and these are accounted for in the EOS (Timmes & Arnett, 1999), and therefore are, in effect, advected with the energy field instead of as fluid electrons included in our definition of .

For energetic purposes we need to be able to track the rest-mass of our material rather than the approximation mentioned above. This is accomplished by tracking the nuclear binding energy per baryon:


where and are the masses of the (free) proton and neutron respectively and is the mass of 1 nucleus of nuclide . Note that is not the atomic mass, which is often given in mass tables and includes electrons and their binding energy. The average mass of a baryon in the fluid is


Thus the actual rest mass density is . Note that because the nuclear binding energy, , is defined with respect to free protons and neutrons in the same proportion as the material, calculation of the average rest mass requires both and , with the latter specifying the overall relative numbers of protons and neutrons in the material.

We may now define the group-specific quantities


so that


It is again important to note that the group-specific quantities such as are not linear combinations of the because the are quotients of linear combinations of the . In order to maintain machine-precision advection of the discretized field , for example, we will need to perform a conservative advection scheme on the product instead of separately on and , since their product will not evolve conservatively to machine precision. A similar statement holds for and .

We will derive the quantities , , , and , , from the initial state. Our simulation begins with fuel of known abundances,


which may vary in space as indicated. These initial abundance will satisfy Equation (1) throughout our simulation; they will have no sources. This allows us to, throughout the burning process, know how much of the local baryons were in what form initially. From these we define the properties of the fuel,


Additionally the ashes of the first stage of burning are assumed to be only a function of the initial composition. Thus

and is therefore also position dependent. Then


As an aside, some concrete examples are useful. In the Townsley et al. (2007) burning model, the initial abundances were , constant in space, and with other abundances zero. Also the ashes of carbon consumption were specified by and , with others again zero. In Townsley et al. (2009) the abundances of the fuel and carbon-consumption ash stages were effectively modified to add a small amount of Ne, whose abundance was still uniform in space, so that the initial abundances were , constant in space, and carbon-consumption ash abundances were {. In the model at hand we will use two parameters and that vary in space to define the initial state, and the are defined as previously. More detailed fuel abundances, or those containing additional major constituents such as Ne or Mg also fit naturally into this scheme.

The fluid properties of the fuel and ashes of just the carbon burning step depend almost entirely on the initial abundances. For the equilibrium groups (QSE and NSE), however, all of these properties change dynamically as the nuclear processing continues at high temperatures. The broad rearrangements of abundances which lead to the variation of properties like and in the more processed ashes are precisely the dynamics that we would like to abstract down to a few parameters for the sake of computational efficiency. To this end, we will treat gross properties of the quasi-equilibrium and equilibrium groups together. For convenience we define another superabundance representing the total amount of material in either QSE or NSE, . This allows us to collect the properties of the equilibrium groups by defining:


The prepended to the quantities here helps to indicate the somewhat odd units involved. For example, is the number of QSE+NSE ions (nuclei) per fluid Baryon, whereas itself is the number of QSE+NSE ions per QSE+NSE Baryon. This unit convention is the most awkward for . To restate why this is desirable: If we had chosen instead to treat directly, that would cause the total nuclear energy , which is a nonlinear combination of and , to not be explicitly conserved by the conservative hydrodynamics scheme. Applying the conservative hydrodynamics to maintains conservation of the total nuclear energy. This form also makes it straightforward to derive appropriate dynamics, which we do below.

Using the progress variables, intermediate state definitions, and QSE+NSE material definitions, the bulk fluid properties can be restated as


This defines the relationship of our burning model variables to the physical fluid properties.

2.2. Posited source terms

The previous subsection developed a framework in which the properties of the parameterized burning stages can be expressed in a way that can be advected in the absence of sources. This leaves us to define dynamical equations (source terms) for the themselves and the properties of the equilibrium materials (the prefixed quantities). By specifying these source terms here, we complete the form of the burning model.

First a brief note on the form of the source terms which we will posit. Typically we will write down source terms by specifying the Lagrangian time derivative


In Eulerian form this gives


Thus we are specifying the term in the evolution of the conserved quantity that is due to transformations rather than just advection.

The evolution of the first stage of burning, , can be set from a flame-tracking scheme or via a thermal reaction rate. This is done just as it is in Townsley et al. (2009),


where is the reaction due to the reaction-diffusion (RD) flame propagation calculation, and is thermally activated carbon-carbon fusion (Townsley et al., 2009). The other progress variables then obey


Here the is the timescale previously determined in Calder et al. (2007) for oxygen consumption. reaches completion at the peak Si abundance, when all oxygen is consumed. The time and length scales for completion of this stage were given as the dashed lines in Figure 1. As can be seen there, this stage is mostly unresolved in our simulations for steady-state detonations, including at all densities important for Fe-group production, where the scales for completion of burning are less than the stellar scales.

In contrast, as shown by the solid lines in Figure 1, the completion of processing of Si- to Fe-group, the phase, can occur on resolved scales for  g cm. Additionally, this stage can be left incomplete for  g cm by the limited length scales in the star and time of expansion of the star. Therefore, the dynamics of this phase are very important for accurate total Fe-group yields. The dynamics we are now using for , Equation (27), differs from that used previously (Townsley et al., 2009; Jordan et al., 2008; Meakin et al., 2009; Jordan et al., 2012b, a; Kim et al., 2013; Long et al., 2014), . In the process of performing the comparisons to benchmark detonation abundance structures presented in section 4.3, it was found that the dynamics used previously led to an approximately exponential relaxation of that did not match the time dependence of the consumption of Si as well as was hoped. In order to improve accuracy of our recorded Lagrangian histories, the dynamics applied to was altered to that of Equation (27). This also necessitates recalibration of the parameter , which will be performed below in section 4.2.

Both of the parameterized timescales above, and , depend on temperature, , and some of the values used below also depend on density. However, there will be significant regions in the artificial flame reaction front – where is not near 0 or 1 – that have a cell-averaged temperature and density that is not a good representation of the temperature of most of the fluid in the cell. These are regions where, as shown in Figure 2, a cell at the reaction front in reality consists partly of unburned fuel and partly of fully burned material separated by a thin front. In this region our use of an artificially thickened reaction front gives a temperature intermediate between that of the fuel and ash. The evaluation of the timescales also needs to be stable as changes to obtain reasonable burning dynamics. By assuming the rest of the burning will occur at either constant density or constant pressure, the final burned state, , , and abundances can be determined based on the current local abundances and thermal state (Calder et al., 2007). The constant pressure prediction provides a reasonable approximation for the final burning state that will be reached by the flame, and so the and of this final state are used to evaluate , , , , and (see below) in regions were . Otherwise, in regions away from the artificial flame the local temperature is used to evaluate and the temperature predicted for an isochoric evolution is used to evaluate , , , and .

Evolution of due to electron capture occurs mainly by conversion of Fe-group material, that is material that has at some point fully relaxed to NSE. At the densities relevant to our SNIa computation, the timescale for relaxation to NSE and the timescale for electron capture are well enough separated that electron capture in material that is only partially fully relaxed to NSE is not an issue. However, due to the artificially thickened reaction front in our SN Ia simulations, a single cell at high densities will consist of an artificial mixture of unburned fuel and fully relaxed NSE ash undergoing electron capture. To constrain electron capture evolution to relaxed NSE material, we separate the components of further into QSE and NSE portions:


For all but the NSE material, . This simplifies Equation (20) to


Applying the chain rule to gives


The right hand side terms each have a distinct physical interpretation. The first is the change due to newly produced material, while the second is due to the adjustment of the pre-existing material. New NSE material is created with and old NSE material evolves according to the tabulated , which naturally gets scaled by the fraction of material currently fully relaxed to NSE, , so that


Next we consider . This represents the average binding energy of all material involved in incomplete Si burning, whether currently in QSE or having progressed fully to NSE. Using the chain rule on Equation (19) splits this into two contributions,


In earlier versions (Townsley et al. 2009 and prior) of this burning model, we posited dynamics in which the binding energy relaxed to the NSE value on the shorter relaxation timescale, . However, in verification comparisons to detonation structures it was found that at low densities this released energy too quickly and led to under-prediction of the temperature just behind the unresolved portion of the detonation front. In order to improve this behavior, we here introduce a that changes as Si-group material is converted to Fe-group, as measured by the progress variable ,


Relaxation toward this value is assumed to occur via capture or photodisintegration, and thus take place on the shorter timescale . To capture these two timescales we posit the following dynamics,


where the evolution on the timescale is contained in . Here represents the of the material at the completion of O consumption, that is the initial QSE state. This state is less easily quantified at high densities, as it may contain a significant, and density-dependent, fraction of particles, but it will only be important at low densities when the Si burning is resolved in the simulation. While, therefore, the most appropriate value for is likely to be density- and composition-dependent, for simplicity we will use , which appears mostly sufficient in the verification tests performed.

The evolution of , or equivalently the ion mean molecular weight, , poses a similar challenge to . Each of the QSE and NSE materials will relax the balance between heavies and /protons/neutrons on approximately the NSQE relaxation time, whereas the conversion between QSE and NSE occurs more slowly. We resolve this by using the scalar that tracks relaxation toward NSE, , to appropriately mix approximations of for the QSE and NSE states and then set our dynamics to move toward this value on the NSQE timescale. Working in a way similar to the construction of Equation (34),

It is left to obtain a suitable estimate of for the QSE state. We found that at densities  g cm that the simple estimate provides a well-behaved approximation that matches produced by benchmark detonation calculations within 10% (see section 4.3). A somewhat complex approximation was proposed in Townsley et al. (2009), but it did not yield a better match to in testing.

The dynamics of our parameterized model for CO burning are contained in Equations (25), (26), (27), (31), (34), and (2.2). The energy release is computed based on conservation of energy, giving the energy release rate per mass,


where , , and are the masses of the proton, electron, and neutron respectively, and is the energy loss to emission of neutrinos based on the local predicted NSE abundances. While the burning dynamics has been stated analytically, the resulting differential equations must now be implemented in a way that is numerically efficient. It is possible to exploit some aspects of the separation between timescales and the strict ordering of the burning stages to make the integration of these dynamical equations extremely efficient. This is discussed in Appendix E.

2.3. Calculation of Nucleosynthesis Using Post-Processed Lagrangian particle Histories

The burning model presented here is intended to give approximately the right energy release, as determined by direct computation of steady-state reaction front structure with large, complete nuclear networks and error-controlled numerical methods, but with a relatively low computational cost. In order to recover detailed abundances, Lagrangian fluid histories are recorded from the hydrodynamic simulation and post-processed. Our post-processing is described in later sections.

The Flash code includes the capability to produce Lagrangian fluid histories through the use of ”tracer” particles (Dubey et al., 2012). These are particles whose position is calculated as


Where the time-dependent velocity field is simply that determined by the hydrodynamic evolution. Generally the number of particles followed and the distribution of initial positions are chosen to provide a sampling that is useful for nucleosynthesis (Seitenzahl et al., 2010), though here we use a simple weighting in which each tracer represents an equal mass and initial positions are chosen randomly to follow the mass distribution. This random distribution is achieved as follows: The domain is decomposed into blocks of cells, where is the dimension, 2 in this case, and we are using blocks that are 8 cells on a side. The mesh structure in Flash provides an ordering for these blocks, called the Morton ordering (Fryxell et al., 2000). We split the mass of the star into segments based on how much mass is contained in each block, using the same order as the Morton ordering. A random number between zero and the total mass is then generated for each particle, and the segment in which it falls determines the block in which that particle is initially placed. A similar procedure is repeated at the block level, using the mass of material in each cell. Once the cell in which the particle will be placed is chosen, each coordinate of the location of the particle within the cell is chosen randomly and uniformly across each dimension of the cell. The impact of the finite sampling represented by this distribution on the uncertainty of our results is discussed in Appendix C.

The Lagrangian tracks are then computed at the same time as the hydrodynamics. The method used to perform the integration of the particle positions is essentially a second-order Runge-Kutta scheme with the velocity field sampled at the end of each combined hydrodynamics and energy source step and linearly interpolated to the particle position. Note that in the directionally-split hydrodynamics solver, which is used here, each hydrodynamics step consists of multiple sweeps of the 1D PPM method to allow for multi-D problems (Fryxell et al., 2000).

3. Deflagration Fronts

Particles representing fluid burned by a deflagration front must be treated differently from those undergoing detonation because the true burning structure differs from the effective one used in the simulation. In some ways treatment of the particles undergoing deflagration is more straightforward because the combustion in the hydrodynamical calculation has been made into a spatially resolved process by coupling it to the RD front as given in Equation (25). The parameterized dynamics used for the RD front are the same as those discussed in Townsley et al. (2009), basically causing the 4-zone wide reaction front to propagate at a specified speed. However, since the flame is generally quite subsonic, with Mach number , it will typically take many time-steps, approximately , for a fluid element tracked by a Lagrangian tracer particle to pass fully through the RD front. In our simulations this is several tenths of a second, as can be seen by the progress variable and temperature histories shown by the solid black lines in Figure 4. During this time, by construction (Townsley et al., 2009, and section 2.2 above), the local temperature is not physical, but a mixture between burned and unburned states in approximate pressure equilibrium. This makes it essential to perform a reconstruction of the portion of the particle’s thermodynamic history during which it is still inside the RD front, before the fully burned state is reached.

Figure 4.— Thermal and burning progress histories for fluid burned by a deflagration front. Shown is the temperature (upper panel thin solid), reaction-diffusion front progress variable (, lower panel dashed) and QSE to NSE progress variable (, lower panel thin solid) recorded at the position of the Lagrangian tracer particle embedded in material ejected at approximately 5 000 km s. and are identical up to about 1.2 s, at which time the fluid element reaches low enough temperature and density that the separation between burning stages begins to become spatially resolvable. Also shown (thick red solid lines) are the reconstructed temperature history used in the post-processing calculation of nucleosynthetic yields and an analog of constructed from the detailed abundances, , computed during post processing (see Equation 38).

The black line in the upper panel of Figure 4 shows a typical temperature history for a tracer particle embedded in material ejected in a DDT SNIa at approximately  km s. The bottom panel shows the evolution of the progress variable representing relaxation toward Fe-group, (solid black line). As can be seen, the transition from unburned to nearly fully burned covers times from about 0.6 s to 1.2 s, and the slow rise in temperature seen in the upper panel covers a similar time range. During this interval, the density and temperature are not representative of a physical burning process, but are instead the average of the burned and unburned states based on the fraction of the cell burned as indicated by the artificially thickened reaction front (see Figure 2). This makes calculation of, for example, the electron-capture history of this fluid element based on a direct post-processing of the , history inappropriate.

We attempt to reconstruct a reasonable approximation to the temperature-density history that a fluid element would have undergone passing through a flame of realistic thickness. The reconstruction of the portion of the fluid history that elapses while the particle is within the artificially broad reaction region is obtained by assuming that the pressure jump across the flame is small, (Vladimirova et al., 2006; Calder et al., 2007). This will be true as long as the Mach number of the flame propagation is low, as is the case for our simulations. Under this assumption, although the local density and temperature are not representative of the actual values, the local pressure should be similar to that near the actual thin flame front to within approximately the Mach number. In order to use this feature, we perform self-heating calculations with a pressure history specified from the fluid histories extracted from the hydrodynamic simulation. This novel mode of specified-pressure-history self-heating was added to the TORCH nuclear reaction network (Timmes, 1999)222Original sources available from Our modifications are available from The set of 225 nuclides used includes all those indicated in the discussion of weak reactions in Calder et al. (2007), which includes an extension to neutron-rich nuclides near the Fe group compared to the standard 200 nuclide set used in TORCH. Weak cross sections were taken from Fuller et al. (1985), Oda et al. (1994), and Langanke & Martínez-Pinedo (2001), with newer rates superseding earlier ones.

Assuming that the fluid element actually crosses the flame front when the progress variable passes through , the reconstructed temperature history is shown by the red curve in the upper panel of Figure 4. It is notable that the temperature peak is much higher and occurs about 0.2 seconds sooner. The initial condition for the trajectory is found by performing a short computation at constant pressure that was raised high enough for the C to begin burning ( K), continuing until the C abundance is 0.1. The specified-pressure self-heating follows this. Once the fluid element exits the artificial reaction front, post-processing can proceed from there using the recorded temperature-density history. We take this point to be when in the recorded history, or when  erg cm, whichever comes first. This corresponds roughly to when burning of heavier elements will cease, when  g cm and , and it is more convenient to impose the condition in than in or directly. In Figure 4, this transition occurs just after  s.

The red line in the bottom panel of Figure 4 shows a progress variable constructed from the full set of species treated in the post-processing,




This effective progress variable measures the process that is intended to track, the conversion of Si-group, or generally intermediate mass elements (IME) to IGE. In NSE, there can also be a significant fraction of light elements (LE, protons, neutrons, ’s) that will be present throughout the transition, but will eventually be captured to form more IGE as the temperature falls. Here the completeness of processing from IME to IGE is comparable between the parameterized burning performed in the hydrodynamic simulation and the post-processed values, with the post-processing giving complete conversion to IGE and indicating more than 95% converted to IGE. The reduction in at late times, starting at approximately 1.8 s, is due to mixing with surrounding zones in the hydrodynamic simulation as the grid is coarsened from 4 to 16 or 32 km cells in order to accommodate the expanding ejecta.

Figure 5.— The time history of for two fluid elements burned by the artificial flame starting at a density of  g cm. The time history on the grid computed from the simplified burning model is shown in black, whereas the history computed in post processing is shown in red (reconstructed portion) and blue (direct post-processed portion). Two fluid elements are treated, one that begins the simulation in the burned state (lower black), and one that the flame passes through after 0.5 s (upper black).

In order to verify that the neutronization is captured well by this method, we turn to 1D simulations in a spatially uniform density medium. For a low Mach number flame in these conditions, a constant-pressure self-heating calculation is a good approximation to the correct fluid history (Vladimirova et al., 2006; Calder et al., 2007; Chamulak et al., 2007). A comparison of the history obtained from the hydrodynamics and the fluid element history post-processed as described above is shown in Figure 5. The histories of two fluid elements are taken from a simulation in which an artificial flame is propagated from a hard wall into 50:50 CO fuel at a uniform density of  g cm with a flame speed of  cm s. The first fluid element begins in the burned region. Its initial state is determined by the Rankine-Hugoniot jump conditions satisfied across the flame front as used to set the initial condition of the simulation. The second fluid element is taken from a position that the flame passes at about 0.5 s. The reconstructed portion of the post-processing is shown in Figure 5 by the solid red portion of the curves and the direct density-temperature history post-processed portion is shown by the blue dashed lines. The black curves show the according to the burning model, Equation (20), at the fluid element position in the hydrodynamic simulation. The agreement is fairly good, with the change in from 0.5 matching within a few percent for both the initially burned case and the case passing through the reconstructed portion at 3 s, a few times longer than expected exposure in an explosion simulation. This provides confirmation that scaling with in the burning model provides a reasonable behavior even with a thickened reaction front. The difference between the time history given by the simulation and the post-processing appears consistent with the use of a larger nuclide set to compute the neutronization rate tables used in the burning model (Seitenzahl et al., 2009b). Using a larger nuclear network for post-processing would improve this difference at some cost to efficiency.

Ideally the histories of the two fluid elements would just be shifted by a time delay based on when their burning began. However, the flame propagation in physical space is slowing somewhat due to the loss of pressure due to neutronization of the earlier burned material. This causes the later burned fluid element to be at a slightly higher density at a given time interval after burning began. The first several tenths of a second of evolution match well in both cases, demonstrating that the post-flame state is consistent with the Rankine-Hugoniot calculation as expected.

4. Detonation Hydrodynamics

In this section we demonstrate the detonation structure we wish to reproduce and we test the burning model in hydrodynamic simulations in comparison to this benchmark. Although it was developed initially for deflagrations in carbon-oxygen mixtures, the reaction structure of detonations is similar enough (Khokhlov, 1983, 1989) that the 3-stage model can also be applied to them. In the simplest form, this just involves identifying the first stage, C consumption, with the rate of the actual reaction, and then following the later burning stages. This was done, for example, in Meakin et al. (2009), and we will do something similar here, with some adjustments for improved accuracy.

As can be inferred from the length scales shown in Figure 1, the actual burning structure is not resolved in full-star simulations. Therefore, somewhat like in the case of the deflagration, the dynamics which lead to the reaction front propagation are not the same in the simulation as in reality. The physics is similar; the energy release determines the strength, and therefore speed, of the detonation shock. However, the acoustic structure in the simulation is not the same as the physical detonation structure. Reactions must be suppressed in the numerically unresolved shock in order to prevent numerical diffusion from dominating the propagation of the reaction front (Fryxell et al., 1989). This creates an artificial separation of a few zones between the shock and the reaction zone. In addition, the reactions may run to near completion within the single zone in which reactions are re-enabled downstream of the shock. We show in Appendix A that the widely-used technique of disabling reactions in the zones adjacent to the shock reproduces the steady state detonation speed and the resolved portions of the reaction structure.

Here we present the error-controlled calculation of the 1D structure of planar detonations that we will use as our benchmark for both the burning model in hydrodynamics and the Lagrangian post-processing. After introducing this benchmark, the remainder of this section will focus on how, in comparison, the burning model acts in hydrodynamic simulations. Post-processing will be discussed in section 5. As already mentioned in the presentation of the burning model in section 2, simply treating according to the C reaction rate and then proceeding as discussed in Townsley et al. (2007) turned out in testing to not reproduce partially-resolved detonation temperature and abundance structures at intermediate densities, - g cm. The successful comparison to benchmarks shown in this section is the result of making the required adjustments to the burning model timescales discussed in section 4.2.

4.1. Verification Benchmark: The ZND structure

In order to evaluate the realism of our simplified model of burning, it is necessary to define an authoritative reference with which it will be compared. Since, as one might expect, no direct experimental validation of nuclear detonations in stellar matter are available, we instead turn to a hierarchical approach to validation (Calder et al., 2002). Following this practice, our interest is in verifying that burning characteristics of our models are similar enough to those computed with methods in which we have more confidence. A typical benchmark in a hierarchical verification like this would be a direct numerical simulation (DNS) of similar phenomenon with more detailed, and typically separately verified, treatments of physical processes. Another source of benchmarks is particular configurations or steady states that can be computed more easily, for example in lower dimension, or in more detail and with better numerical error control.

As one of the two combustion modes in SN Ia explosions, the predicted outcome of C-O detonations have been discussed in some detail previously in the astrophysical literature. Khokhlov (1989) presented an overview of the microscopic structure of steady-state planar C-O and He detonations at a variety of densities. Further work by Sharpe (1999) extended calculations of the structure of the planar steady-state structure and products beyond the sonic point in the detonation wave, allowing the completion of burning to be computed at a wider range of densities. Sharpe (2001) followed this up with computations of detonation speeds and structure for non-planar, i.e. curved, detonation fronts in steady state, still in one dimension. Gamezo et al. (1999) and Timmes et al. (2000) investigated the multi-D structure of C-fueled detonations with high resolution reactive hydrodynamics for cases important for SNe Ia. Recently, Domínguez & Khokhlov (2011) performed a high-resolution investigation into the stability of C-fueled detonations in 1 spatial dimension at low densities,  g cm.

We are interested here in an inherently transient phenomena as the detonation traverses different densities within the star. As a result, the ideal benchmark is simulations of the reactive Euler equations which include all relevant nuclides (and therefore all relevant reactions) and in which all important length scales are resolved. The component models of such a DNS have been separately validated in many contexts, and their limitations are fairly well understood. Unfortunately, a DNS is challenging for the nuclear processes under consideration here. In order to fully capture the reaction kinetics, it is necessary to include hundreds of species. The more severe limitation, however, as demonstrated in Figure 1, is the large separation of time and length scales between the final reaction stages – those which process Si-group to Fe-group elements or perform electron captures – and the reactions which drive the burning front forward, fusion of carbon. At the densities of most interest, where the nucleosynthetic processing to Fe-group is incomplete due to the finite size of the star, a few  g cm, these length scales are  cm and  cm respectively.

In this work, we will compare our results with those obtained from the well-known Zel’dovich, von Neumann and During (ZND) model of detonations (Zel’dovich, 1940; von Neumann, 1942, 1963; Döring, 1943; Fickett & Davis, 1979). This model predicts both the detonation velocities and final products as well as the detailed 1D thermal and compositional structure in space for steady state denotations. It can also be computed with error-controlled methods with a large reaction network including all relevant reactions. Matching these detailed structures during burning is crucial for our application. The yield of the supernova will be determined by the burning processes that lead to these structures. Therefore, if our burning model, including particle post processing steps, can accurately reproduce the abundance profiles predicted by the ZND model, it increases our confidence in the yields that it predicts in more general cases.

The ZND equations describe the detonation structure between the detonation shock front and the sonic point. Beyond the sonic point, where the following flow is moving away from the detonation front at the local sound speed, disturbances cannot move upstream to change the detonation flow. The portion of a propagating steady-state detonation between the shock and the sonic point is a static (i.e. time invariant) structure which propagates in space at the detonation speed. The flow beyond the sonic point is typically not static, and its form depends on the boundary condition of the following flow. Sharpe (1999) computes the flow beyond the sonic point for asymptotically free propagation, but we do not undertake that here.

Before moving further in the discussion, it is useful to state the ZND equations in the form in which we will use them, for plane-parallel, steady-state detonations (Fickett & Davis, 1979; Khokhlov, 1989). In the frame of the detonation front,


Here dot indicates an ordinary time derivative, , is the flow velocity (with respect to the detonation front), is the unburned density, is the detonation speed, is the frozen (evaluated with constant ) adiabatic sound speed, is the pressure. To be consistent with the above conventions, is the number of nuclei of nuclide per fluid baryon. Thus , where is the mass number of nuclide . The are given by the nuclear reactions. The energy release function, in the absence of weak interactions, is


where is the internal energy. The integration of these equations is begun just behind the leading shock, whose properties are related to those of the fresh fuel by the detonation speed and the usual shock conservation equations.

A diagram of the form of typical solutions are shown in Figure 6. Equation (42) is singular at the sonic point, where , unless is also zero there. There is a large class of solutions for which is high enough that the entire following flow is subsonic. That is, , and therefore , changes sign from negative to positive before increases to , thus avoiding an encounter with this singularity. This type of solution has a higher pressure in the final state than in the reaction zone, and is called ”overdriven” or ”supported” since it is effectively being pushed from behind by an overpressure. In this case the full flow, including itself, has an inherent dependence on this boundary condition. As the pressure in the final state, or at the ”piston” following the detonation, is decreased, also drops, and eventually a sonic point will appear. For detonations with lower pressures in the following flow, the steady portion of the detonation flow then becomes an eigenvalue problem such that at the sonic point.

In simplified reaction systems, at the sonic point because that is the point at which fuel consumption completes. This is called a Chapman-Jouget detonation (Fickett & Davis, 1979), and its speed can be computed from just the energy release and the EOS, without a need for the full ZND equations (Khokhlov, 1989; Gamezo et al., 1999). For reaction systems with complex or reversible reactions or changes in mean molecular weight, the heat release function may not reach or cross zero at a unique level of progress toward the fully burned state. That is, may be attained before burning is ”complete” and a static final state reached. In this case, the sonic point, and thus the end of the static portion of the detonation profile, also occurs before a stable final state is reached. Such a detonation is termed ”pathological” or ”eigenvalue” and the sonic point, where the singularity appears in the ZND equations, and where, therefore , is called the pathological point. This is, in fact, the more common case, and eigenvalue detonation structures in this case represented a major advancement manifest by the ZND model (Fickett & Davis, 1979).

Figure 6.— Detonation pressure structure in space and time for a steady-state detonation in 1 dimension. The lower solid and dashed lines that pass through the pathological point represent possible solutions for the cases in which the sonic point is reached before completion of burning. These detonations are termed ”pathological” or ”eigenvalue” detonations.

The ZND integration can be continued after passing through the pathological point, but there is more than one way to exit this point (Sharpe, 1999). Figure 6 shows a diagrammatic representation of the relation of the pressure profile in overdriven and self-sustained, or unsupported, detonation. The lowest overdriven detonation which can be fully integrated using just the ZND equations without traversing a singularity is that which passes just above the pathological point. While it is possible with special methods to traverse the pathological point and obtain the self-sustained solution (Sharpe, 1999; Moore et al., 2013), we do not undertake this here due to our large set of species and complex reactions. This seems prudent because even some of the profiles obtained by Sharpe (1999) using this method show clear indications of having further zero-crossings of beyond the pathological point. How these would manifest in the detonation structure is unclear from this level of analysis.

It is now possible to choose a well-defined verification benchmark problem whose solution can be calculated with both the ZND model with a fairly complete reaction set and a 1D hydrodynamic simulation with our simplified burning model. We choose our benchmark to be the slightly-overdriven state found by tuning to be a small amount above the eigenvalue that leads to the pathological point. This configuration can be replicated in a 1D hydrodynamic calculation by manipulation of the boundary conditions in the following flow to have the appropriate pressure in the fully burned state. The static portion of the benchmark structure, between the shock and the sonic point, can be also used as a reference solution for self-sustained detonations once they reach steady-state.

4.2. Calibration of Timescale for Si Consumption

In order to make use of the simplified dynamics for the transition from the QSE to the NSE state, given by Equation (27), we must calibrate the timescale . In Calder et al. (2007), was calibrated by computing the consumption timescale in isochoric self-heating as a function of temperature and then using a fit to that timescale for . Here we will compare the time evolution of Si group element abundances for our benchmark detonation, computed using the ZND equations, directly to the behavior posited in our burning model by Equation (27).

Figure 7.— Effective progress variable for the conversion of Si- to Fe-group material, or relaxation toward full NSE (see text for definition). Driven detonations computed by a ZND integration are shown at two pre-shock densities of  g cm (black, shorter timescale), and  g cm (red, longer timescale). Also shown are fits to the dynamics in the current burning model (dashed lines) and exponential relaxation (dot-dashed lines).

In order to make comparisons we use based on Equation (38), where the are the abundances from the ZND calculation computed with a large network for a driven solution. In Figure 7, is shown for two densities spanning the range of interest, 0.5 and  g cm. From Figure 1 we see that at these densities the synthesis of IGE from IME will occur as a partially or mostly resolved process on the grid during the explosion of the star, and will largely determine the IGE yield of the explosion. Expansion times for the star are in the range of a few tenths of a second and the hydrodynamic timestep is around  s for typical simulation resolutions of a few km.

As will be shown in section 4.3, the early rise to in both curves is due to IGE+LE produced during the oxygen consumption stage. Therefore we will proceed by fitting the latter part of the curve only to get a better characterization of the transition timescale. If necessary, the inclusion of some IGE in the intermediate state, in Figure 3, could be introduced in converting the progress variables to abundances. However since we use post-processed yields for our final abundances this is not necessary.

In the C-O burning process, the stages are well enough separated in time that oxygen consumption, which is complete about the same time the Si abundance peaks, completes before the transition from Si- to Fe-group proceeds very far. This can be seen clearly in Figure 1 as the 5 orders of magnitude separating the time of maximum Si abundance (dashed red line) and the completion of burning (solid red line). We therefore assume and for a characteristic value of we may analytically integrate Equation (27) to obtain


Here and are taken from at the end of oxygen consumption and in the final state respectively. In this case they are about 0.15 and 0.99. might be different if we performed this calibration with different initial abundances. This form can now be fit to the curves shown in Figure 7 using a non-linear least squares fit. We use a fitting region , to capture the major portion of the evolution. The resulting fits are shown by the dashed lines. The fit timescale is not sensitive to choice of ; a 5% variation in only changes the fit by 1%. The maximum error in the fits occurs when , and is about 0.06 and 0.03 for the higher and lower density shown in Figure 7 respectively. We will discuss below in section 4.3 how well the resulting burning model performance in hydrodynamics compares to the detonation benchmark, and extend this comparison to abundances in post-processing, compared to those in the benchmark, in section 5.

This fitting procedure has been repeated at several densities, between 0.3 and 10  g cm. At each of these densities the conversion of Si- to Fe-group takes place at a declining temperature. The decline during this burning stage is, however, much less than the variation from one density to another. By evaluating the temperature when the relaxation is approximately half complete, we can construct and fit a relation between and . We obtain


The timescale found here is not directly comparable to previous work because we have used different burning dynamics. However, a similar fit can be performed with the exponential decay form that results from the simpler dynamics previously posited, (Townsley et al., 2007). This is shown by the dot-dashed lines in Figure 7 when fit to the same region indicated above. This form does not appear to provide a good reproduction of the late-time behavior of . Also the timescales obtained for the exponential fit are approximately a factor of 10 to 20 shorter than those given for in Calder et al. (2007). This is understandable because the definition used in that work measured a timescale to reach a fairly complete burning stage, whereas we have fit an exponential form directly.

4.3. Comparison of parameterized burning hydrodynamics against 200-nuclide ZND Structure

The verification that we are attempting to perform involves demonstrating that the abundance structure produced by post-processing particle tracks from the hydrodynamics which utilizes the parameterized burning matches the ZND structure for a steady state detonation. That comparison will be done in section 5, but first it is useful to compare the intermediate result obtained from the parameterized burning model in the hydrodynamics simulation alone. This will provide a check on the realism of spatial thermodynamic structure without the added complication of the integration of the Lagrangian tracks, and also give some diagnostics concerning whether the parameters within the burning model are behaving as expected.

Figure 8.— Detonation thermal and composition structure displayed by our parameterized model for C-O burning in 1D hydrodynamic simulations compared to our detonation benchmark of the steady-state ZND solution of the equivalent detonation. This case is at a density of  g cm and an initial composition of 50% C, 48% O and 2% Ne. Left: Thermal structure of simulations at spatial resolutions of 4 km (red) and 1 km (blue) compared to the benchmark steady-state ZND solution (black). Right: Top: The composition structure of the benchmark ZND calculation computed with a 200-nuclide network. Middle: Effective progress variables derived from the abundances in the benchmark (red, Equation (47); blue, Equation (48); black, Equation (38)), compared to the progress variable obtained in the hydrodynamic simulations at 4 km (green) and 1 km (magenta) resolution. Bottom: Average number of nucleons per nucleus, , derived from the full abundances in the benchmark (black) and obtained from the progress variables, Equation (21), in 1D hydrodynamic simulations at 4 km (green) and 1 km (magenta).
Figure 9.— Similar to Figure 8 but for  g cm.

Our benchmark is, as described in section 4.1, the ZND solution for a steady-state, planar, slightly overdriven detonation in 1 dimension. This solution is shown as the reference curves in Figures 8 and 9, with thermodynamic quantities, , , , in the left panel (black), and abundances in the top right panel. The initial condition for the 1D hydrodynamic simulations is material at spatially constant density and temperature away from the ignition point. We consider cases here with this background temperature set to  K. The domain extends from to  km in order to allow the detonation to approach steady state. Two resolutions, 4 km and 1 km, similar to the resolution of production supernova simulations (Townsley et al., 2009), are used to confirm insensitivity to resolution. We will refrain from using the term convergence here, reserving it for circumstances in which gradients are resolved. The boundary condition on the opposite end of the domain from the ignition is reflecting, but has no impact on the simulation due to the supersonic nature of the detonation and since the simulation is stopped before the front reaches it. The left boundary, at , is a zero-gradient boundary. The initial perturbation is made in both temperature and velocity. Along with inflow from the zero-gradient boundary, the latter will serve to support the detonation from behind. Both temperature and velocity are placed as linear gradients decreasing from a maximum at to the background values of  K and velocity of zero over a size we will call the size of the ignition region. The velocity is tuned by hand until the pressure far behind the detonation front and near the boundary matches the late-time pressure found for the slightly overdriven ZND solution. Ignition region sizes were 1024 km and 128 km for and  g cm respectively.

As above, we will focus on densities at which the transition from Si-group burning products to Fe-group burning products is fully or partially resolved on the spatial grid. At a density of  g cm, as indicated by Figure 1, nearly the entire Si- to Fe-group transition is resolved at 4 km resolution for the steady state detonation. The spatial structure obtained from the ZND calculation and the hydrodynamics, which uses the parameterized burning, is shown in Figure 8. The thermodynamic quantities, , , and , are shown in the left panel. The hydrodynamic result is at an evolution time of 5.45 seconds, when nearly the entire domain has been consumed. The zero point for the distance behind the shock in the hydrodynamic simulations is taken as the last zone in which the shock detection considers the cell inside a shock, thereby suppressing the reactions in that zone. See Appendix A for more on this suppression. In steady state, the shock region in which the reactions are suppressed is a well-localized region of approximately 4-5 zones. As a result of this, the first point from the hydrodynamic simulations, indicated with stair-stepped lines, is at 4 km and 1 km for simulations of those respective resolutions. The top right panel shows the spatial abundance structure of a selection of nuclides for the steady-state detonation from the 200-nuclide ZND calculation. From the ZND abundance and thermal structures shown in Figure 8 we see that both the C and O consumption stages are entirely unresolved because they take place on length scales of approximately 1 cm and several  cm respectively. In the span of less than a single zone, the burning reaches the Si-rich QSE.

The values of , , and at the point chosen as zero distance behind the shock in the simulation are not quite the same as the post-shock values expected based on the detonation speed. This is presumably the result of numerical mixing in the vicinity of the under-resolved shock and burning front. The post-shock density is about 35% lower than the peak value predicted by the ZND calculation and the pressure, rather than peaking at the shock, peaks in the first zone in which burning is allowed at a value about 10% lower than expected. The peak, which also occurs in the first zone in which reactions are allowed, is about 3% higher than the peak in the benchmark. This transient is also larger in time and space than the true burning structure due to the resolution, but the thermal state appears to relax back toward a good approximation of the QSE state very quickly, within 2 zones. After this and a small undershoot, the hydrodynamic solution is a very good match, with 3% in and , and within about 1% in , all the way out to the pathological point. There is noise on a similar level, but more so in and than . An artifact of the initial ignition is evident at the end of the hydrodynamic curves for and . We also find very good consistency between resolutions after the first few zones behind the shock, matching within a percent, with noise in each case slightly larger than that. The hydrodynamic result is probably not completely relaxed to the steady-state overdriven solution, since there is no pressure minimum. However, the pathological point occurs quite close to the end of the domain even for this large domain and the pressure minimum is expected to be fairly shallow.

A comparison of some of the parameters in the burning model are shown in the lower right two panels in Figure 8. In order to make a comparison of the progress variables we have defined some effective progress variables for the 200-nuclide set. In addition to Equation (38) above, we define


The spatial structure of both and are unresolved at this density and these resolutions. Thus they are both 1 in the first zone behind the shock-detection suppression of burning because our data dumps always follow a reaction sub-step in our operator split time evolution. For this reason the and from the hydrodynamic simulations are not shown.

We find a good match between the evolution of and the effective equivalent defined for the 200-nuclide set. The largest discrepancy is due to the production of some Fe-group material with Si-group in the benchmark. After the discrepancy is less than , and after it is less than . This indicates that our temperature-dependent fits of the timescales for this evolution, described in section 4.2, are acting satisfactorily. The bottom right panel of Figure 8 shows how the mean ion molecular weight compares to the equivalent quantity from the parameterized burning . This quantity is systematically about 4% low, probably due to our choice of as an estimate of the of the QSE state. The QSE state is not pure Si, and therefore this estimate is slightly off and creates a systematic offset in the consecutive evolution toward . The difference observed in the test may also be magnified by the hydrodynamic simulation having not quite reached the steady overdriven state. In either case the discrepancy in only leads to less than 1% discrepancy in , as found above, so this level of agreement appears sufficient for producing accurate thermodynamic histories for particle post-processing.

As a second case, shown in Figure 9, we perform a similar calculation at an ambient density of  g cm. At this density, more than half of the transition from the Si-group dominated QSE to the Fe-group dominated NSE is unresolved on a 4 km grid. This is according to the profile of predicted by the 200-nuclide ZND calculation, shown in the middle right panel of Figure 9. We see a region, similar to that in the first case, of about 2 zones in which is about 3% higher than the expected peak and and are intermediate between the expected post-shock values and the QSE values, after which all of these relax to within 3% of the benchmark values. The largest source of discrepancy is due to the lack of the expected minimum near the pathological point at a distance of  cm behind the shock. Instead, the hydrodynamic solution monotonically relaxes to the state given by the boundary condition pressure. However, even with this discrepancy the maximum difference between the benchmark and hydrodynamic result is about 5% in and and less than 2% in . As before the two resolutions match very well, within 0.5%.

In terms of progress variables, during the partially-resolved transition from Si- to Fe-group, the progress variable for this process, , is about 0.1 higher than the benchmark predicts at a given distance behind the shock. This seems like a reasonable indication of the uncertainty of the progress variable’s reproduction of the real process for partially-resolved cases like this one. The determined in the final state by the burning model in the hydrodynamics is only about 2% lower than the benchmark. However, as for the thermal profiles, the non-monotonic behavior near the pathological point is not captured.

The main difference from the benchmark in this case is due to the lack of a clear pathological point in the hydrodynamic result. It is unclear if this is due to the limited resolution, a deficiency in the burning model, or insufficient time to relax to the steady state. In any case, the discrepancy in the thermal quantities used for post processing is, at maximum, a fairly modest 5% in and 2% in . We will accept this as the approximate uncertainty in the thermal histories produced by the burning model, and proceed to investigate the abundances produced in post-processing directly below. At higher densities than about  g cm, as can be seen from the length scale for completion of Si- to Fe-group conversion, the conversion will be nearly complete on scales smaller than the resolution. The burning model shows good reproduction of the final state, within a few percent, so that denser cases should also have similar good accuracy.

5. Verification of Lagrangian particle nucleosynthesis against ZND solution

While it is important that the progress variables provide a good reproduction of the detonation structure, in the end the yields will be computed by post-processing Lagrangian tracer particle histories. In this section we compare computed Lagrangian track yields to the steady-state ZND solutions that we are using as a benchmark. Detonation yields are computed by a direct integration of the , history recorded by the Lagrangian tracer particle from the hydrodynamic simulation, using them to set the reaction rates in the nuclear reaction network.

Figure 10.— Abundance histories computed from post-processed Lagrangian histories from hydrodynamic simulations (dashed lines) compared to benchmark steady-state detonation structures computed from the ZND equations (solid lines). Shown as mass fractions are the major abundances after C consumption, O (blue), Si (green), and Ni (black), as well as Ca (orange) and Fe (brown). Two densities are shown,  g cm  (left column) and  g cm (right column). Each of these is shown from a simulation with 1 km (top row) and 4 km (middle row) spatial resolutions with histories recorded at full time resolution, and at 4 km spatial resolution with history recorded at a reduced time resolution of  s (bottom row).

The results of the integration of the reactions over the Lagrangian history are compared with benchmark calculations in Figure 10 for the same two densities, and  g cm (left and right columns), and two spatial resolutions, 1 km and 4 km (top and middle row), as used in Section 4.3. We also consider a case with a reduced time resolution for the recording of the Lagrangian history (bottom row). Comparison can now be made directly with actual abundances. We show the major abundances for stages beginning at oxygen consumption, O, Si, and Ni, as well as the major neutron-rich nuclide produced before other Fe-group material, Fe (Bravo et al., 2010), and the spectroscopically important Ca. Each plot shows two curves for each nuclide: the benchmark solution (solid lines) computed using the ZND equations and the post-processing of the , history (dashed lines).

In order to compare structures we must choose a zero time during the Lagrangian history. Zero time for the benchmark ZND integration corresponds to the downstream side of the shock. We have chosen the zero time for the Lagrangian history to be at the first timestep that reaches 1% above the ambient temperature. This makes the entire reaction region visible on these plots because the C and O consumption timescales in the benchmark are shorter than the timestep in all cases. The abundances in the first part of the reaction region are unrealistic, as expected. Notably at the