SN Ia dependence on {}^{22}Ne

Evaluating Systematic Dependencies of Type Ia Supernovae: The Influence of Progenitor Ne22 Content on Dynamics

Abstract

We present a theoretical framework for formal study of systematic effects in Supernovae Type Ia (SN Ia) that utilizes two-dimensional simulations to implement a form of the deflagration-detonation transition (DDT) explosion scenario. The framework is developed from a randomized initial condition that leads to a sample of simulated SN Ia whose Ni masses have a similar average and range to those observed, and have many other modestly realistic features such the velocity extent of intermediate mass elements. The intended purpose is to enable statistically well-defined studies of both physical and theoretical parameters of the SN Ia explosion simulation. We present here a thorough description of the outcome of the SN Ia explosions produced by our current simulations.

A first application of this framework is utilized to study the dependence of the SN Ia on the content, which is known to be directly influenced by the progenitor stellar population’s metallicity. Our study is very specifically tailored to measure how the content influences the competition between the rise of plumes of burned material and the expansion of the star before these plumes reach DDT conditions. This influence arises from the dependence of the energy release, progenitor structure and laminar flame speed on content. For this study, we explore these three effects for a fixed carbon content and DDT density. By setting the density at which nucleosynthesis takes place during the detonation phase of the explosion, the competition between plume rise and stellar expansion controls the amount of material in nuclear statistical equilibrium (NSE) and therefore Ni produced. Of particular interest is how this influence of content compares to the direct modification of the Ni mass via the inherent neutron excess as discussed by Timmes, Brown & Truran (2003). Although the outcome following from any particular ignition condition can change dramatically with content, with a sample of 20 ignition conditions we find that the systematic change in the expansion of the star prior to detonation is not large enough to compete with the dependence discussed by Timmes, Brown & Truran (2003). In fact, our results show no statistically significant dependence of the pre-detonation expansion on content, pointing to the morphology of the ignition condition as being the dominant dynamical driver of the Ni yield of the explosion. However, variations in the DDT density, which were specifically excluded here, are also expected to be important and to depend systematically on content.

Subject headings:
hydrodynamics — nuclear reactions, nucleosynthesis, abundances — supernovae: general — white dwarfs
10

1. Introduction

Type Ia supernovae (SN Ia) are bright stellar explosions spectroscopically distinguished by strong silicon features near maximum light and a lack of hydrogen features (Filippenko, 1997; Hillebrandt & Niemeyer, 2000). Motivated broadly by the importance of SN Ia for cosmological studies (Phillips, 1993; Riess et al., 1996; Albrecht et al., 2006), contemporary observational campaigns are gathering information about SN Ia at an unprecedented rate (BAOSS, Li et al. 1996, 2001; LOSS, Treffers et al. 1997; Li et al. 2001; SNLS, Astier et al. 2006; CSP, Hamuy et al. 2006; Nearby SN Factory, Copin et al. 2006; Skymapper, Keller et al. 2007; ESSENCE, Wood-Vasey et al. 2007; STRESS, Botticella et al. 2008; SDSS-II, Holtzman et al. 2008; SXDS, Furusawa et al. 2008; Totani et al. 2008; PQ, Djorgovski et al. 2008; CfA, Hicken et al. 2009a; CRTS, Drake et al. 2009; PTF, Kulkarni et al. 2009; Pan-STARRS; the La Silla SN Search) and the future promises even more (the Dark Energy Survey, LSST, JDEM, see Howell et al. 2009a). Interesting systematics have been discovered and continue to be better characterized. All SN Ia appear to burn a similar amount of total material, but can differ widely in the amount of Ni produced, an effect that is closely correlated with their brightness and decline time (Mazzali et al., 2007). Observations indicate that there are two populations of SN Ia that differ in the elapsed time between star formation in the host galaxy and the explosion (Scannapieco & Bildsten, 2005; Mannucci et al., 2006). These two populations also appear to have a different average brightness, with SN Ia occurring in active star-forming galaxies appearing brighter (Howell et al., 2007). Recent observations find correlations suggesting that SNe Ia in galaxies whose populations have a characteristic age greater than 5 Gyr are  1 mag fainter at maximum luminosity than those found in galaxies with younger populations, while progenitor metal abundance has a weaker influence on peak luminosity (Gallagher et al., 2008). The color properties of observed SNe Ia have become an important prospective source of systematic error for cosmic measurements (Hicken et al., 2009b). Alongside the finding that the scatter from the Hubble law depends on host galaxy metallicity (Gallagher et al., 2008), this creates concern that corrections for extinction have become entangled with uncharacterized intrinsic color variation. It is in this context of extensive empirical searches for systematics that our broad but still incomplete theoretical understanding of these events must be pushed as far as possible, to attempt to predict and understand systematic trends.

By invoking a straightforward counting argument based on the fact that the number of protons and neutrons are approximately conserved in the explosion, Timmes et al. (2003) argued that there should be a fairly robust metallicity effect on the average Ni yield of SNe Ia, and therefore, potentially, their brightness. Motivated by the simplicity of this effect and the important implications for cosmological usage of SNe Ia, significant effort has gone into measuring such a metallicity dependence observationally, but a clear effect has proven elusive (Gallagher et al. 2005; Gallagher et al. 2008; Howell et al. 2009b). Constraining metallicity dependence alone is challenging for two reasons beyond the fact that the effect does not appear to be large. It is fairly difficult to measure accurate metallicities for the parent stellar population, and even so there are strong systematic problems with such measurements due to the mass-metallicity relationship within galaxies (Gallagher et al., 2008). Additionally there is an apparently much stronger dependence of mean brightness of SN Ia on the age of the parent stellar population (Gallagher et al., 2008; Howell et al., 2009b).

The most recent observational work leaves the situation still murky. Gallagher et al. (2008) were unsuccessful at finding a metallicity dependence of the average brightness, but did find a such a dependence of the Hubble residual (scatter from the Hubble law) obtained from light-curve fitting. In contrast, Howell et al. (2009b) found a slight dependence of average brightness on metallicity and no similar Hubble residual issue using a different light-curve fitting method. However, the effect is difficult to separate from the dependence on mean stellar age. Given such uncertainty, it is important to continue to evaluate the influence other aspects of the explosion might have on the effect of the simple overall neutron excess outlined by Timmes et al. (2003). The importance is further stressed because it is clear that other effects, including the intrinsic random variation of otherwise similar SNe Ia, can be even larger in magnitude.

The systematic effects of metallicity on the SN Ia outcome have been the topic of a number of theoretical studies with a variety of explosion models. Several of these were accomplished in one dimension (Höflich et al., 1998; Iwamoto et al., 1999; Höflich et al., 2000; Domínguez et al., 2001), and therefore bear revisiting with multi-dimensional models. In addition to changing the overall yield of Ni, the initial neutron excess, as set by the metallicity, is important for the composition of intermediate layers that act as an important opacity source during the photospheric phase of the supernova (Domínguez et al., 2001). The challenges of understanding the SN Ia phenomenon with multi-dimensional numerical models have since somewhat overshadowed this type of study of population systematics. For example, one multi-dimensional study by Röpke et al. (2006) studied systematics using a less parameterized supernova model, but as a result was forced to treat marginally successful deflagration models that do not produce realistic explosions. Additionally, the neutron excess was treated only in post-processing, thus excluding sensitivity to dynamical effects like those studied here.

The number of possible systematic parameters to consider, along with the intrinsic scatter expected from implementation of realistic turbulence characteristics, leads us to develop a theoretical framework for formally evaluating systematic effects and their possible statistical significance. Systematic effects that require study include both physical ones such composition and ignition density and purely theoretical ones such as detonation condition and flame model parameters. In addition to the deflagration-detonation transition (DDT) explosion model itself, the basic component of this framework is an initial condition that defines a theoretical population or ensemble from which we can draw a sample of supernovae and study how the sample as a whole responds to parameter changes, giving a more complete and statistically quantifiable picture of their systematic impact.

The ultimate goal of the study of SN Ia systematics is to understand the dependence of properties of a SN Ia population on characteristics of the parent stellar population. Among other benefits, this understanding would allow known characteristics of the stellar populations of galaxies, such as their cosmic history, to be utilized in understanding systematic trends that may appear in SN Ia, or to better understand the chemical evolution of galaxies. Stellar populations are most basically characterized by their age and metallicity, or more realistically, some mixture of or distribution of these parameters. There are likely to be several other important secondary parameters such as binarity, or environment (e.g. ionizing background) that may have enough influence on a stellar population to be reflected in its SN Ia population. The dominance of stellar population age and metallicity (Scannapieco & Bildsten, 2005; Mannucci et al., 2006; Gallagher et al., 2008; Howell et al., 2009b), or, alternatively, host galaxy morphology or color, in observational investigations to date indicates that such secondary factors are probably small by comparison. However, until the age dependence is fully understood, it is important to be mindful of possible additional contributions.

Our aim in the present paper is far more modest. We begin by considering the dependence of the supernova on the composition of the exploding WD, and refine this exploration to only particular dynamical aspects in order to allow a targeted, conclusive result. The question that we seek to answer is: does the change in the dynamics of the explosion due to a different Ne abundance in the progenitor introduce a significant systematic effect in addition to the neutron excess discussed by Timmes et al. (2003)? This highly targeted scope is motivated by a desire give a clear and extensive description of our framework. This is the first time that we are applying two-dimensional DDT simulations as a method for generating a semi-realistic theoretical sample of supernovae and studying the systematics of that sample.

The DDT model proved to be one of the most successful of the one-dimensional SN Ia models (e.g. Höflich & Khokhlov, 1996). However, it was never satisfying that both the deflagration velocity and the DDT transition density were treated essentially as free parameters. The hope of multi-dimensional models is that burning propagation during the deflagration phase, which was necessarily parameterized in one dimension, can be calculated directly. This would remove another free parameter and lead to more reliable models. Unfortunately the manifestation of buoyancy instabilities in multi-dimensional models became a serious challenge. Even modest asymmetries in the initial conditions of the deflagration led to far too little expansion of the star by the time that a traditional DDT would occur (Niemeyer et al., 1996; Calder et al., 2003, 2004; Livne et al., 2005). This can be ameliorated with particular choices of ignition condition, allowing the main desirable features of the 1-d models to also be obtained from multi-dimensional models (Golombek & Niemeyer, 2005; Gamezo et al., 2005; Röpke & Niemeyer, 2007; Bravo & García-Senz, 2008).

The degree to which such a symmetric deflagration is physical is still a matter of some debate. The distribution of flame ignition regions in time and space remains fairly uncertain (Woosley et al., 2004; Röpke et al., 2006) as well as the degree to which turbulence-flame interaction after ignition can influence the subsequent spread of the flame (Röpke et al., 2007; Jordan et al., 2008). When placed alongside the physical uncertainty as to whether or not a transition to detonation can actually occur (Niemeyer, 1999), in a certain light the DDT scenario may simply be too contrived. However, so far the evidence is not sufficient to disprove it, and it continues to provide one of the best prospects for matching the observations of the typical SN Ia.

We begin in Section 2 by discussing the variety of physical effects through which composition can lead to systematic variation of SN Ia properties. This forms the context for how and why we limit the scope of this first study to dynamic effects only. Following this, in Section 3, some improvements to the numerical model of the explosion are discussed and the extensions necessary to treat detonation-flame interaction and the presence of are presented. The burning model, flame speed, and mesh refinement are each discussed. In Section 4, we present the ignition condition on which our ensemble of SNe Ia is built and use it to construct a framework for evaluating the significance of systematic effects. The varieties of yield arising in the theoretical SN Ia population are discussed along with how the initial condition appears to control the outcome. A metric for measuring the expansion prior to DDT is introduced and calibrated to reflect the mass of nuclear statistical equilibrium (NSE) material synthesized. Section 5 presents a detailed description of the outcome of a representative two-dimensional simulation and how it generally compares with observations. Some points of the implementation of the explosion model that necessitate referring to the details of the explosion are also described here. Finally, in Section 6 the framework is applied to measure the effects of content on the expansion of the star at the time of the DDT. We then summarize our conclusions and discuss future work in Section 7.

2. Systematic Influences of Composition

We begin with an overview of the physical effects via which composition can influence the process and outcome of a DDT explosion. The most important material constituents of the WD are C, O, and Ne, and we will frame our discussion in terms of these. Composition can influence the explosion through changes in several physical processes and properties involved in the explosion. These include ignition density, DDT density, energy release, flame speed, WD structure, and neutron excess. Detailed treatment of the first two of these, ignition density and DDT density, will be deferred to future work. The others will all be treated here, with the last, neutron excess, taken as the baseline effect to which others should be compared. Several of these effects involve significant uncertainties, and it is useful to highlight each in turn.

The compositional structure of the carbon-oxygen WD that explodes is determined principally by the post-main sequence evolution of the star of which it is a remnant (Domínguez et al., 2001). The inner several tenths of a are formed during the convective core helium burning phase, and the layers outside this in shell burning on the asymptotic giant branch (Straniero et al., 2003, and references therein). During helium burning, the C, N, and O which was present in the initial star is transformed into Ne (Timmes et al., 2003), leading to a direct dependence on metallicity of the parent stellar population. In addition, Ne is formed in the pre-explosion convective carbon burning core at a comparable abundance (Piro & Bildsten, 2008; Chamulak et al., 2007).

Ignition density – The ignition density characterizes when, as the result of accretion, the WD core begins runaway heating due to carbon fusion outpacing neutrino losses (Nomoto, 1982; Woosley & Weaver, 1986). The central density at flame ignition, when convection is insufficient to spread the heating throughout the core, is slightly less (see e.g. Piro & Bildsten 2008), and is also often called the ignition density, the meaning usually being discernible from context. The flame ignition density is generally near  g cm for successful models of SN Ia (e.g. Nomoto et al., 1984; Höflich & Khokhlov, 1996). As we do here, Höflich et al. (1998, 2000) adopted a single value for the ignition density for their studies of composition dependence. Although we exclude it here, the variation of ignition density is expected to be a significant effect. Along with the energetic variations due to the carbon content discussed below, it is likely an important contribution to the observed dependence on stellar population age, which has proven to be stronger than any metallicity dependence (Gallagher et al., 2008; Howell et al., 2009b).

The precise value of the ignition density is sensitive to several factors, each of which has remarkable uncertainties. One of the reasons we will simply fix the ignition density for this study is the variety of uncertainties which must be considered if it is varied. The energy generation rate depends on the C abundance in the core, as set by the evolution of the star that formed the WD. Authoritative calculation of this abundance remains elusive due to its dependence on the details of convection during late core helium burning (e.g. Straniero et al., 2003). The core temperature of the WD is also important, and so the thermal history of the core, notably the accretion rate and possibly properties of the helium flashes (Nomoto et al., 1984). Either carbon composition or thermal state could lead to metallicity dependencies that are closely involved with both the evolution of the parent stellar population and the still very poorly understood (Branch et al., 1995) process of progenitor system formation. There are also uncertainties in the screening enhancement of nuclear reactions at these high densities (Gasques et al., 2005; Yakovlev et al., 2006), and in the reaction rates themselves. In particular, the reaction cross-section is not experimentally determined down to the low center- of-mass energies relevant for ignition in white dwarfs. There is evidence, from heavy-ion fusion reactions, of “hinderance”—a suppression of the astrophysical S-factor—at sub-barrier energies (Jiang et al., 2007; Gasques et al., 2007). In the case of , however, resonances are predicted in the energy range of interest (Michaud & Vogt, 1972; Perez-Torres et al., 2006), which could significantly increase the cross section by as much as two orders of magnitude (Cooper et al., 2009) at temperatures .

DDT density – We make the presumption that there is a unique characteristic density at which there is a transition from deflagration to detonation. In future work, this will be extended to comparison of flame width and turbulent state (Niemeyer & Woosley, 1997; Khokhlov et al., 1997; Golombek & Niemeyer, 2005), giving a more physical transition point. While the details of this transition remain difficult to fully quantify (Aspden et al., 2008; Woosley et al., 2008), it is very likely that it will have direct dependencies on composition because both the reactivity and the flame width depend on both the and abundances. Note that in this case, the important composition is that in the outer portions of the star. In addition to having a higher abundance than the center (e.g. Höflich et al., 2000) because it is the product of shell burning, this region will have a lower abundance because it will not have participated in the core convection during which is enhanced in the convection zone (Piro & Bildsten, 2008; Chamulak et al., 2008). Given the outcome of this study, showing that dynamical factors have little impact on the NSE yield, the influence on the DDT density is expected to be the principal avenue, beyond inherited neutron excess, for dependence on metallicity. The laminar flame studies of Chamulak et al. (2008) predict a 10% reduction in the DDT density for a increase in the abundance of 0.02. From results found here, we estimate that this might correspond to a 3% decrease in NSE mass. This is about half as strong as the dependence from inherited neutron excess.

Energy release – The amount of energy release per unit mass can affect both the global dynamics of the expansion and eventual disruption of the WD, as well as the more local dynamics related to the buoyancy that drives the acceleration of the burning front during the deflagration phase. The composition is essential for both these effects in two ways. First, the gross nuclear energy available per unit mass is mostly sensitive to the abundance of . Due to mixing during the smoldering phase, the abundance over a broader region of the WD is involved in this case, as opposed to precisely at the center as is the case for the ignition density above. Thus, the products of both the helium core and shell burning are important. Secondarily, the abundance, by changing the balance of protons and neutrons, can influence the NSE state to which material will flash. More neutron rich material favors more tightly bound material and therefore releases more energy (Calder et al., 2007).

Flame speed – The propagation speed of the subsonic burning front propagated by thermal diffusion, the flame, is sensitive to both the and abundances due to their effects on the energy release and the speed of the early stages of the nuclear fusion (Timmes & Woosley, 1992; Chamulak et al., 2007). This laminar propagation speed is likely to be much less important after the strong turbulence that results from the buoyancy instability develops. Turbulence has the ability to add a nearly arbitrary amount of local area to the flame surface, therefore making the burning effectively independent of the laminar propagations speed (Niemeyer & Hillebrandt, 1995). However, there is a period of time early in the deflagration phase of the supernova where the interaction of the laminar flame speed with the lower strength turbulent velocity field from the core convection will be important for setting the morphology of the burned region at the point when strong buoyancy takes over (Zingale & Dursi, 2007). In this work, we find an important sensitivity to the assumed outcome of these earliest stages.

WD structure – The WD is supported by pressure of degenerate electrons. The neutron excess, and therefore the abundance, sets the number of baryons worth of weight each electron must support. Thus a WD of the same mass with a higher neutron excess will be more compact because it will have fewer total electrons. (Conversely at the same central density, a star with more , and thus more neutrons, will be slightly less massive.) This has a concomitant effect on the density distribution throughout the star. This is a fairly small effect (Höflich et al., 1998), but important to treat appropriately due to the marginally bound nature of a near-Chandrasekhar mass WD.

Neutron excess – In addition to the indirect effects mentioned in relation to energy release above, the neutron excess, as set principally by the content, has a direct influence over the final nucleosynthetic products, particularly the amount of Ni. Timmes et al. (2003) showed that the decrease in the Ni produced in the explosion, absent other factors, should be linearly proportional to the content and therefore the original metallicity of the stellar population. The distribution of will again be important, with that in the inner bulk of the star being important for the gross yield, and that in surface layers for opacity sources in those regions of the ejecta.

In order to bring the scope of our study within tractable limits, it is necessary to either neglect or exclude some of these effects with assumptions. Most notably we would like to avoid for now the possibly quite complex dependence of ignition density on composition. This can naturally be done by considering only a single value of , and assuming that we are only studying the variation for progenitors that resulted from the same formation history. Note that since metallicity changes main sequence lifetimes, for example, this simplification may be a rather unnatural one in contrast to comparing SN Ia at the same stellar population age (Höflich et al., 1998). Our second simplifying assumption will be to neglect the dependence of the DDT density on composition. Exploring the effect of DDT density will be undertaken in immediate future work along with the compositional inhomogeneity with which it is intertwined. These two sets of assumptions leave us to study the dependence arising from how changes in abundance modify the energy release, flame speed, WD structure and neutron excess. Because Timmes et al. (2003) have provided an excellent discussion of the direct influence of neutron excess, we focus mainly on how these other effects may modify the robust conclusion they reached.

3. Improvement and Extension of Numerical Model

The essential components of the numerical code used in this work were presented in Townsley et al. (2007) and work referenced therein, but we will give a brief overview. The overarching code is formed by the FLASH Eulerian compressible hydrodynamics code (Fryxell et al., 2000; Calder et al., 2002) with modules added for the nuclear burning which occurs in SNe Ia. FLASH uses a high-order shock-capturing compressible hydrodynamics method, the piecewise parabolic method (PPM, Colella & Woodward, 1984), adapted to treat a general equation of state (EOS, Colella & Glaz, 1985). We use FLASH’s tabulated fully ionized electron-ion plasma EOS (Timmes & Swesty, 2000; Fryxell et al., 2000). FLASH applies this hydrodynamics method on an adaptively refined, tree-structured, non-moving Eulerian grid. We make extensive usage of this adaptive mesh refinement (AMR) capability, using different resolutions for burning fronts (4 km), the initial hydrostatic star (16 km), and the region of negligible density initially outside the star (as coarse as thousands of km). Extensive detail on our refinement scheme and tests of resolution are given in section 3.3 below.

The nuclear burning processes, beginning with carbon fusion, and extending to electron capture in material in NSE are implemented with a nuclear energetics model described by Townsley et al. (2007), and here in section 3.1, and calibrated to reproduce the features of nuclear processes that occur in SN Ia as calculated using hundreds of nuclides (Calder et al., 2007). This burning model is used for both subsonic (deflagration) and supersonic (detonation) burning fronts. Because the flame physics that governs the propagation speed of the subsonic burning front is unresolved at 4 km (Chamulak et al., 2008), we propagate this front by coupling our nuclear energetics to an artificial, resolved reaction front given by the advection-reaction-diffusion (ARD or ADR) equation (Khokhlov, 1995; Vladimirova et al., 2006). A variety of special measures are necessary to ensure this coupling is acoustically quiet and stable, and therefore appropriate for simulating the buoyant instability of the burning front (Townsley et al., 2007). Detonation fronts are handled somewhat naturally by the shock-capturing features of the hydrodynamics scheme (Meakin et al., 2008). The only common components between our code and that of Plewa (2007) are those publicly available as components of FLASH, which excludes all components treating the nuclear burning; differences are discussed in Townsley et al. (2007).

In this section we discuss various changes made to improve the overall consistency of the numerical modeling as well as the extensions necessary for the burning model and flame speed treatments to include an arbitrary Ne content. Changes to the burning model, flame speed and refinement are discussed in turn, with separate subsections for improvements and extensions. The improvements include a backwards-differenced time integration method for the nuclear energetics model, more well-defined methods for defining the density used to compute the flame speed and the Atwood number used to calculate the buoyancy-compensated (”turbulent”) flame speed, as well as simplified and better tested mesh refinement prescriptions. Extensions include modifications to the coupling of the RD front and nuclear energetics in order to better account for the interaction of the detonation front with the artificially thickened sub-sonic burning front, and inclusion of the dependence of the laminar flame speed on composition.

3.1. Burning Model

Although the model of nuclear energy release used here is functionally identical to that of Townsley et al. (2007), it is useful to restate it in a way that makes adjusting the abundance of the fuel easier and that is more amenable to backwards differencing or sub-stepping in time. This will also provide the context in which we can present our prescriptions for how the detonation interacts with the artificially broad flame front.

Simplified Dynamical Equations

The first step is to abstract the properties of the fuel and the ashes of carbon burning into adjustable parameters. The properties of interest are the number of electrons per baryon, , the number of fluid ions per baryon, , and the average nuclear binding energy per baryon, . These can be constructed from mass fractions per

(1)
(2)
(3)

where is the nuclear charge, is the atomic mass number (number of baryons), and is the nuclear binding energy, where and , , and are the rest masses of the proton, neutron and nucleus of nuclide . Note that because the quantities , and are ”per baryon”, it is more clear to think of the density on the grid as representing the baryon number density. The advantage being that, unlike mass, this is actually a conserved quantity, making it possible to calculate the overall change in rest mass energy in a well-defined way. Baryon number divided by Avogadro’s number makes an extremely good approximation for the mass in grams, and this correspondence will be used so extensively that in situations where the difference is of no consequence, we will almost always refer to the mass density instead. In this interpretation the are actually baryon number fractions and also become conserved quantities in the absence of transformations.

We begin by defining the properties of ”pure” fuel or carbon-burning ashes as and respectively. As before, we also define three progress variables, which measure the progress of the burning through various stages. These are , for processing of fuel to carbon-burning ashes, , for processing of these ashes to Nuclear Statistical Quasi-Equilibrium (NSQE), and finally for relaxation of NSQE to full nuclear statistical equilibrium (NSE). Finally we will define the properties of the NSQE+NSE material, in such a way that the bulk properties are

(4)
(5)
(6)

The usage of, for example, rather than , as was done in Townsley et al. (2007), is to avoid issues of how the evolution of the NSE material should be treated when is very small.

As discussed in Townsley et al. (2007) it is possible to calculate a ”final” NSE state by calculating the endpoint of an isobaric or isochoric burning based on the instantaneous local conditions. If we denote this state by ”NSE” and also use the NSQE and NSE timescales parameterized based on the expected temperature of this final state, and (Calder et al., 2007), we can posit Lagrangian source terms.

(7)
(8)
(9)
(10)
(11)
(12)

Here is the contribution from the reaction-diffusion model for the burning front propagation (Townsley et al., 2007),

(13)

where is the progress variable for the (A)RD front and the coefficients and are determined from the prescribed propagation speed of the front, , and its desired width. , with tunable parameters , and , provides a stable and acoustically quiet reaction front propagation as discussed in Townsley et al. (2007). Thermal fusion is included via the reaction term

(14)

where is the carbon-carbon fusion rate from Caughlan & Fowler (1988), is the local density, and is the carbon mass fraction in pure fuel. While this rate is subject to several important uncertainties as described in section 2 above, the features of the detonation are fairly insensitive to its precise value, depending mainly on the overall energy release. Some additional prescriptions for treating the temperature used to calculate this rate are required in regions where the RD front is active, , which will be discussed below. In order to accurately propagate the detonation, is set to zero inside shocks (Fryxell et al., 1989). Finally, is a very rough estimate of the ion abundance for NSQE material given by

(15)

Conservation of mass energy yields the evolution of the mass-specific energy due to these nuclear processes,

(16)

where is obtained using (6) and the above source terms, , , and are the proton, neutron and electron masses respectively, is the speed of light, and is the energy loss rate to neutrinos due to weak processes in the predicted NSE state (Calder et al., 2007; Seitenzahl et al., 2009).

All the dynamical quantities in equations (7)-(12) are mass-specific, or more appropriately baryon-specific, quantities. Their full hydrodynamic evolution is thus given, for a representative quantity, , by

(17)

where is the velocity field. The portions indicated as and are treated in an operator split fashion, each acting consecutively, with being the conservative hydrodynamic operator implemented with Piecewise Parabolic Method (PPM) (Colella & Woodward, 1984) and the burning scheme described above. Because the operator also does not change , equations (7)-(12) reduce to simple time derivatives after operator splitting. For the application of the operator, we further assume that the temperature used in the calculation of as well as the timescales and properties of the NSE final state are constant, so that equations (7)-(12) can be backwards-differenced and solved algebraically. This is intended to increase the stability of the treatment, particularly when is similar to the size of a hydrodynamic time step. While the assumption that is constant is not a good one for the thermal term, this is likely to have fairly little impact as this term is largely only active in propagating the detonation, for which the burning length is unresolved. The NSE final state prediction should only vary slowly with time, even in burning fronts, since it depends on the combination or for the isochoric and isobaric predictions respectively. The former combination only varies on the hydrodynamic timescale and the timescale of weak processes, and the latter should vary fairly slowly due to the subsonic propagation of the RD front.

Detonation-Flame Interaction

The RD front that we utilize to model the propagation of the subsonic burning front is several computational zones thick in order to enable its stable propagation. This presents several complications when it is desirable to also treat thermally activated burning that is not related to the flame front. At high densities, the physical flame is very thin, and as a result, there is very little partially burned material. Thus, in this density regime, partially burned material on the grid represents regions where burned and unburned material are well separated, but the interface is not resolved. In addition, this means that the temperature of the zone is not representative of the temperature in the unburned material, which should be used in the calculation of . In previous work (Meakin et al., 2008), the thermally activated term, in equation (7) was suppressed (set to zero) within the RD front, where exceeded some small threshold. This is not too bad an approximation in the case of the gravitationally confined detonation (GCD) scenario because only a small portion of the star is burned via the deflagration, so that an even smaller portion is left in an unrealistic partially burned state. However, such a prescription would leave a large amount of partially burned material in a simulation in the DDT scenario.

The principal intent of allowing a thermal contribution within the RD front is to allow the detonation to consume nearly all of the partially burned material as it encounters the RD front, as it would if encountering a thin flame. In order to accomplish this, two measures are applied. First, is suppressed where unless within one flame width (4 zone widths) of a given cell. This prevents spurious thermal burning. Second, where is enabled within the RD front, an estimate of the temperature of the unburned material is used in place of the local temperature in the calculation of . This temperature estimate is obtained by estimating the properties the local material would have in the absence of the RD-front based fuel consumption. The fraction burned due to thermal processes is . The properties in the absence of the RD front would then be

(18)
(19)

where we estimate

(20)
(21)

Then the temperature estimate for the thermally activated burning, , is found by evaluating the equation of state for the local density, , and mass-specific internal energy given by . With a sufficiently accurate and robust estimate of the temperature, the first condition, suppressing in the RD front except in proximity to other burning, might prove unnecessary. This was not the case in tests we have run so far. The above prescriptions appear sufficient for the purposes of this study.

3.2. Flame Speed

The improvements to the flame speed treatment are intended only to improve numerical consistency by obtaining a more well-defined and stable value for the front propagation speed. The improvements are to estimate the unburned density and to use a tabulation of the Atwood number instead of a crude estimate. In doing this, we also develop methodology for testing the accuracy of tabulations to ensure that the tabulation grid is sufficiently fine. Following this, we present details of our extension of the computation of the laminar flame speed to account for dependence on fuel composition. This is one of the physical effects that our study of Ne systematics is evaluating. The tabulation utilized for the laminar flame speed is also evaluated for accuracy.

Consistency Improvements of Flame Speed Treatment

To be consistent with the burning model, the input flame front speed is a function of the local chemical composition of fuel and the density of that fuel. We define the flame front speed to be the speed of the carbon burning front with respect to the fuel. We parameterize the chemical composition of the fuel by and mass fractions, with the remaining material being O. In order to prevent the flame from being torn apart by turbulence induced by the Rayleigh-Taylor (RT) instability, we introduce a minimum flame speed,

(22)

where is the Atwood number for carbon burning, is the gravitational acceleration, is the grid resolution, and is a calibrated parameter determined to be  (Townsley et al., 2008). In this work, we use . The Atwood number and gravitational acceleration describe the strength of the RT-induced turbulence. The resulting speed used to set the propagation of the ADR front is

(23)

where is the speed of the laminar flame front. is then rolled off to zero for densities below  g cm. In reality, the RT-induced turbulence wrinkles the flame, increasing the surface area of the burning front and accordingly increasing the burning rate. The minimum flame speed for our model thus serves to compensate for this boosted burning rate (Khokhlov, 1995). By construction, this prescription does not properly account for the interaction of isolated turbulence with the flame surface. Thus while the integrity of a given plume is maintained, the interaction of turbulence created by one plume with others is not necessarily captured accurately. We perform resolution studies and vary the parameter to check for precisely these issues (Section 5.3) and find that they do not appear to adversely affect the principal metric measured in this study, the expansion prior to DDT. Inclusion of more direct models of flame-turbulence interaction (Colin et al., 2000; Schmidt et al., 2006) is planned for future work.

Both and are best characterized by a dependency on the composition and density of the unburned fuel. This presents a difficulty when must be constructed in a partially burned cell, where the density and composition has changed from that of the unburned material due to burning and energy release. For the composition, we simply store separate mass scalars that represent the initial composition: the mass fractions of material that was initially in the form of and . While the density varies across the subsonic burning front, the pressure is approximately constant. We therefore use the local pressure, , to estimate the density of the unburned state. This is greatly simplified by the fact that the unburned material in the WD is at a high degree of degeneracy, so that the temperature of the unburned state is not too important and simple forms for the pressure-density relation of a degenerate gas can be used (Hansen & Kawaler, 1994). Our estimate for the unburned density is

(24)

where is the pressure in cgs units (erg cm), is the electron fraction, and 0.92 is adjusted to provide a good fit to the pressure-density relation in the initial star. The difference between this estimated density based on the local pressure and the actual density is no more than 6% in the initial star for  g cm and varies smoothly from an overestimate at low pressures to an underestimate at high pressures.

The Atwood number used in the calculation of the minimum flame speed is determined based on the local fuel composition and the resulting energy release from carbon burning. Given these parameters, we estimate the density of the ash using the Rankine-Hugoniot jump condition across the flame front (Vladimirova et al., 2006). For computational efficiency, we calculate the Atwood number at the beginning of the simulation for twenty-one equidistant log densities between 6 and 9.6 and ten equidistant carbon mass fractions between 0.3 and 1.0 using a representative neon mass fraction appropriate for the simulation. Because the Atwood number changes less than over the relevant range of , we reduce the dimensionality of the interpolation and memory requirements by introducing a characteristic neon mass fraction, which for this study is equal to the global neon mass fraction. During a simulation, the code performs a bi-linear interpolation of the Atwood number given the local initial mass fraction of and the estimated unburned density from Eq. (24).

In order to test the accuracy of this interpolation procedure, we estimate the uncertainty for the Atwood number throughout the table. Our method of uncertainty estimation solves for the second-order term in a polynomial interpolation of the table in each direction. Because the Atwood number is a smooth and slowly varying function of and as shown in Figure 1, we can assume that higher order terms in the Taylor expansion are negligible. In this case, the second-order term serves as a correction to the linear interpolation in a particular direction. We use this correction to define the uncertainty estimate

(25)

where the Atwood number, , is calculated using a 1st and 2nd order interpolation. The second-order interpolation, , is performed in only one variable-dimension, . By subtracting the first-order from the second-order interpolation, we are left with a term proportional to the second partial derivative with respect to . This allows us to assess the curvature of our table in each direction and whether we have enough points in our table in each dimension such that linearly interpolating the data provides an accurate estimate of the Atwood number. For the Atwood table, we are able to calculate an exact value to compare against the interpolated estimate. However, for the flame speed table discussed later, we are not able to calculate an exact value to test against. We use the Atwood table to verify our method of uncertainty estimation used later to analyze our flame speed table. The comparison of the uncertainty estimate, , to the actual error, , in Table 1 verifies our method of uncertainty estimation.

Figure 1.— The Atwood number is shown as a function of log density for several carbon mass fractions. This displays the sufficient smoothness of the Atwood number as a function of log density and carbon mass fraction such that our method of uncertainty estimation is valid.

The results of performing these tests at the grid points and midpoints of the Atwood table are shown in Table 1. We show the minimum, maximum, absolute average, and average uncertainty estimate, , in each direction of the Atwood table. In columns 1 and 2, we evaluate the uncertainty estimate at a grid-edge (a midpoint along holding the other variable at tabulated values). In column 3, we evaluate the uncertainty estimate at a cell-centered point (where both variables are at midpoints). In this case, we sum the uncertainty estimates in each direction, . We do not perform a bi-quadratic interpolation. While the cross-term is potentially important in this case, we do not include it in the uncertainty estimate of a cell-centered point due to the dependence on the order of interpolation in multi-dimensional polynomial interpolations. In columns 4 and 5, we calculate the actual error, , for any and , respectively. In this study, we use models with . The error is calculated at all midpoints and grid points (the same values used for the uncertainty estimates). Due to the way we initialize the Atwood table, the test points are equivalent to forty-two equidistant log-densities between 6.0 and 9.6 and twenty equidistant carbon mass fractions between 0.3 and 1.0 with .

Because the table has positive curvature in and negative curvature in , we expect to find the maximum magnitude in the uncertainty estimate on a grid-edge as opposed to a cell-centered point. Our method of uncertainty estimation provides a maximum expected uncertainty of -10.99% which is verified by the true maximum error of -9.99%. Both of these values were evaluated at and . This shows that our method of uncertainty estimation works for the Atwood table and can be applied to the flame speed table that will be discussed next. For the purposes of this study, the relevant maximum error is obtained from comparing the interpolation against calculated Atwood numbers at midpoint log-densities all at .

(%) (%) (%) (%) (%)
-0.98 3.30 0.00 0.03 -0.15
-10.99 6.88 -5.19 -9.99 -7.60
4.82 5.52 1.93 4.11 3.69
-4.82 5.52 0.82 0.95 -3.10
Table 1Error estimates for tabulating and linearly interpolating the Atwood number

Composition Dependence of the Laminar Flame Speed

For the laminar flame speed, we give preference to values calculated by Chamulak et al. (2007) using a 430-nuclide reaction network for a variety of initial carbon and neon mass fractions and a range of densities. Similarly to the Atwood numbers, at runtime the code performs a linear interpolation to obtain the laminar flame speed from a table of previously calculated results. The method used by Chamulak et al. (2007) is not well suited to solve for the laminar flame speed at low density; therefore, we use the results from Timmes & Woosley (1992) to supplement this table. Because the flame speeds for each case were calculated using different initial carbon mass fractions, we merged the data by linearly interpolating Timmes & Woosley (1992) values onto the Chamulak et al. (2007) grid. The results of this merger are more clearly shown in Figure 2.

Figure 2.— Shown are laminar flame speeds as a function of log density from Chamulak et al. (2007), Timmes & Woosley (1992), and our resulting flame speeds from merging these datasets all at zero Ne mass fraction. Each panel compares these datasets at different carbon mass fractions, . (Note that Timmes & Woosley (1992) performed their study at . We expect slightly higher laminar flame speeds at .)

Some discussion of the accuracy of this method of obtaining flame speeds at low densities is warranted. Timmes & Woosley (1992) do not track Ne dependence with their method of determining the laminar flame speed. For identical points in parameter space, the two methods produce laminar flame speeds that differ on average by . This difference is of the same order as the effect of adding Ne where we see a speed-up depending on density. For the models considered in this study, low densities occur near the surface of the white dwarf star. In these regions, the input flame speed (23) is dominated by the Rayleigh-Taylor driven turbulence such that . In fact, this transition occurs at  g cm in the initial star. Therefore, this tabular method of estimating the laminar flame speed is sufficient at low densities for this study.

The tri-linear interpolation occurs within the three-dimensional parameter space of carbon and neon mass fraction and log-density to calculate the laminar flame speed. We cannot easily evaluate the uncertainty in our interpolation method by comparing with a direct calculation of the laminar flame speed due to the computational cost. Therefore, we apply our method of uncertainty estimation discussed for the Atwood number to the laminar flame speed table, Eq. (25). We calculate an uncertainty estimate, , for our method of interpolation at the quarter-points and grid-points for each parameter in our table. We limit our analysis to densities above  g cm. The laminar flame speed becomes unimportant below  g cm because the buoyancy-compensated flame speed takes over in Equation (23) at roughly this density. The results of these calculations are given in Table 2 and Table 3. Table 3 shows results for , which is more relevant for this study, while Table 2 shows the behavior of our flame speed table in general. As was the case for the Atwood number, the uncertainty estimate along a single variable dimension, , is evaluated along a grid-edge (at the quarter-points of holding the other variables at tabulated values). This is the case for columns 1-3. In column 4, we evaluate cell-centered points (at the quarter-points of and holding at tabulated values). To evaluate the uncertainty estimate on cell-centered points, we sum the individual uncertainty estimates in each direction. In column 5, we calculate the uncertainty estimate at the quarter-points in all directions such that , , and are not at tabulated values.

For the flame speed table, both the and variable dimensions have negative curvature, while has mostly positive curvature. Because the uncertainties along each grid-edge are comparable, we expect the maximum uncertainty in our table to occur either along a grid-edge or on a cell-centered point. We do not expect a maximum uncertainty estimate involving and any of the other two variables due to opposing curvatures.

(%) (%) (%) (%) (%)
0.00 0.00 0.24 0.45 0.39
8.74 -10.96 5.81 8.63 8.27
1.14 0.25 1.76 2.73 2.68
1.01 -0.22 1.76 2.73 2.68
Table 2Error estimates for tabulating and linearly interpolating the laminar flame speed (above )
(%) (%) (%)
0.00 0.70 0.72
0.62 3.35 3.34
0.06 1.82 1.83
0.06 1.82 1.83
Table 3Error estimates for tabulating and linearly interpolating the laminar flame speed (above and )

We determined that the laminar flame speed table has sufficient resolution at densities above with the magnitude of estimated uncertainties as shown in Table 2. The estimated uncertainties relevant for this study at are as shown in Table 3.

3.3. Mesh Refinement

The goals of the design of our refinement scheme are to be as simple as possible while both capturing interesting or physically important features and doing so with good efficiency. These three goals are somewhat at odds, and therefore provide a wide latitude for choosing refinement prescriptions. We proceed by defining three regions of the physical domain:

  1. fluff: regions with

  2. star: non-fluff regions,

  3. energy generation:
    regions with or

where is the progress variable in the reaction-diffusion model of the flame front. These are indicated with a ””, ”” or ”eg” subscript respectively. We assign to each of these a consecutively increasing maximum refinement level. For simplicity, we will here use the minimum cell size rather than the refinement level, .

A “fluff” region outside the star is necessary because the hydrodynamics method in FLASH has no explicit mechanism for treating empty (zero-density) cells or free surfaces. To ameliorate this, would-be empty cells are filled with material of extremely low density which will not effect the dynamics of the more dense material of interest in the simulation (see also Zingale et al., 2002). Here we set it to  g cm. Because this material is of no physical interest and has negligible contribution to the dynamics of other material, is taken as large as possible, generally being the total domain size divided by the block size. Recall that the smallest independently refinable unit in PARAMESH is a ”block”, which in all our simulations is a cell region.

The finest resolution will be treated as the resolution of the simulation, though in fact before the burning spreads, most of the star is only refinable to . Some care must be taken in defining the thresholds such that regions near the thresholds do not cycle between refinement levels. Cycling is avoided by checking both parent and child blocks (i.e. blocks at the finest and next to finest refinement levels) reconciling the results and using some hysteresis in defining the boundaries. For the simulations presented here  g cm throughout the simulation. If a simulation is run far enough into expansion, this restriction would need to be relaxed to lower density so that stellar material remains refined as it expands. The energy generating thresholds are placed at  erg g s and  s for the deflagration and detonation phases, but are relaxed to  erg g s and infinity, respectively, for calculations of the expansion phase.

Beyond the definition of these regions of maximum refinement, refinement is triggered by gradients in the physical quantities. Here we use the density and the first progress variable in the burning model, , to trigger refinement. This choice guarantees that all burning fronts are refined to the degree allowable. In order to detect gradients, we utilize the built-in tests included with PARAMESH and FLASH. These measure a ratio between the second and first derivatives of the fields being checked. The reader is referred to the source code for the exact generalization to multiple dimensions. A threshold for triggering refinement is then defined for each quantity checked, the default being 0.8. In order for the initial star to be fully refined, it was necessary to drop this threshold to 0.1 for the test of the density field. With this threshold, de-refinement set in during expansion in one-dimensional simulations, but not appreciably so in two dimensions.

In order to verify that sufficient resolution was being utilized and that the adaptive refinement was not adversely affecting the accuracy of the solution, a one-dimensional convergence test was performed. Convergence is expected in one dimension because instabilities are suppressed. The absence of the instabilities that accelerate the burning in a multidimensional simulation necessitates a significant artificial enhancement of the burning rate in order to unbind the star. Because this is only a numerical test, we simply increased the parameter (Equation [22]) until an explosion was obtained. Then the combination , where is the finest resolution, generally , is held fixed as the resolution is varied. We used  km. We are not interested in the order of convergence here, only that reasonable convergence is obtained. Therefore, we simply compare the density distribution of the outgoing ejecta to a uniformly refined case with additional levels of refinement,  km. This is the solid line in Figure 3, which shows density distributions for simulations at various resolutions at a time of 5.6 seconds after ignition. The fluff was not refined in any of these simulations, but a test was performed with the entire domain refined to confirm that this had no effect in the density range shown here. The same initial WD was used in all cases, which was created on an 8 km grid and mapped onto the new grid by averaging density and temperature without a reconstruction (or equivalently a piecewise constant one).

Figure 3.— Density distribution of outgoing ejecta in radius for one-dimensional test of resolution. Convergence is obtained when energy-generating regions are refined to 4 km and the rest of the star is initially refined to 16 km (dashed line). The reference case in which the entire star is uniformly refined to 1 km is shown by the solid line, while cases in which each of the refinement regions are relaxed by one level from the coarsest converged values are shown by the dot-dashed and dot-dot-dashed lines.

We found that a resolution of  km and  km was necessary to satisfactorily match the reference solution. Also it was necessary for the initial star to be fully refined at , and the trigger threshold for the density gradient test was adjusted to achieve this as described above. Subsequent de-refinement as the star expanded did not appear to have any adverse effect on the outcome of the simulation. The de-refinement threshold was set to 0.0375. The case with these parameters matches the reference case extremely well down to a density of about  g cm. Most likely the initial resolution is insufficient to resolve the hydrostatic equilibrium in these low density layers. For comparison, we show in Figure 3 the cases with a factor of 2 too little resolution separately for the energy generation region and the star. We find the outcome is most sensitive to the refinement level of the star. If this is too low, the hydrodynamics of the stellar expansion are not properly captured.

This convergence study also provides an important check on the model of nuclear energy release, particularly demonstrating that sub-stepping between hydrodynamic steps does not appear necessary for such a simple reaction model, as it likely would for even a highly reduced nuclear network. It is possible that lower resolutions of the energy generation region could be made viable by introducing a sub-stepping mechanism.

4. Framework for Evaluating Systematic Dependencies

As a starting point for our study of systematic effects, we must consider a model of the supernova explosion that reproduces many observed characteristics. In light of its success in one-dimensional work (e.g. Höflich & Khokhlov, 1996), we take the deflagration-detonation transition as this starting point. The defining feature of this scenario is that the flame front will undergo a transition to detonation when turbulent mixing on the scale of the flame front becomes more vigorous compared to the propagation speed of the flame front (Niemeyer & Woosley, 1997; Khokhlov et al., 1997). One-dimensional work such as Höflich & Khokhlov (1996) utilized the density at which the DDT occurs, , as a non-unique parameter. Variation of then led to the observed variety of SN Ia outcomes. Instead of this, we make the assumption that the conditions that lead to the DDT, while dependent on local properties (e.g., composition and turbulence strength) are otherwise unique, though not currently known with precision. Variety in the outcome for the same initial WD relies on the ignition configuration being non-unique due to the turbulence present in the convective core at the time of deflagration ignition. In this section we construct the essential component of our framework for evaluating systematic dependencies: a theoretical population of SNe Ia obtained by sampling a defined ensemble of ignition conditions. Some basic characteristics of this population are discussed. The next section, 5, describes a typical outcome in detail along with additional technical details of the implementation of the explosion.

The ability of isolated plumes to rise to low densities from a central ignition before the star expands significantly, as demonstrated by a number of studies (Calder et al., 2003; Plewa et al., 2004; Calder et al., 2004; Livne et al., 2005; Plewa, 2007; Townsley et al., 2007; Röpke et al., 2007; Jordan et al., 2008), presents a significant challenge to producing a realistic explosion with a DDT. If the transition to detonation takes place as early as these simulations imply, only a very thin layer of Si-group elements, the hallmark of a Ia spectrum, are produced. However, Röpke et al. (2006) found that a multi-spot deflagration ignition led to a much more symmetric deflagration. Building upon this, somewhat like other authors (Röpke & Niemeyer, 2007; Bravo & García-Senz, 2008), we here are able to obtain realistic DDT explosions by considering deflagration ignition conditions chosen to have no low-order power.

4.1. Constructing a Theoretical Population of SN Ia

In order to evaluate systematic dependencies in the SN Ia population, we require a theoretical sample of supernovae that mimic the properties of supernovae as a population. SN Ia possess an intrinsically large scatter in the Ni mass synthesized during the explosion, ranging typically between about 0.3 and 0.8. (See Howell et al. 2009b for one example sample distribution.) With the DDT scenario, we have found that this degree of scatter can be obtained by introducing a certain degree of randomness into the initial condition of the deflagration phase. We now motivate and describe this initial condition and then describe the samples that it produces.

Our initial condition is motivated as a possible abstraction from prospective three-dimensional studies of the very early deflagration phase. During the first  s of the deflagration, the flame surface will be heavily influenced by the pre-existing convection field in the core. In order to develop a randomized sample of a variety of ignition conditions, we choose to parameterize the possible outcomes of this early evolution as harmonic structure in a flame surface when it reaches a modest distance from the core. We place the position of the initial flame surface before perturbation at a radius of  km. To this base radial position, perturbations with definite harmonic content are added in the form of spherical harmonics. In the case of axisymmetry, only spherical harmonics can be included, reducing to Legendre polynomials. However, this technique should extend in a straightforward way to three dimensions. The initial position of the flame surface is defined to be

(26)

The normalization convention of Jackson (1999) is used. The amplitudes are chosen from a normal distribution. The maximum harmonic content is chosen so that the perturbations in the flame surface are resolved. For the 4 km resolution simulation performed here, we choose . As demonstrated by a number of studies of localized, off-center ignition in the absence of a convection field (Calder et al., 2003; Plewa et al., 2004; Calder et al., 2004; Plewa, 2007; Röpke et al., 2007; Townsley et al., 2007; Jordan et al., 2008), low-order perturbations on such a flame surface would lead to an early DDT with too little expansion of the WD to give realistic yields. To prevent this effect, low order modes are left out by choice of . Here we have used , which appears to give a reasonable sample (see below). An extensive study of the sample that arises from different values was not undertaken, so the sensitivity of the sample to this choice is uncertain. One additional restriction was imposed due to the axisymmetry. The perturbation, the second additive term in equation (26), is restricted to be negative on the symmetry axis, , . This suppresses the formation of slender ”jets” of burned material, which are likely artifacts of the axis singularity. Such jets are much smaller than the rising bubbles observed by Townsley et al. (2007) and similar simulations, and bear more resemblance to the ”tails” seen in those cases.

Within the framework defined by equation (26), besides and , it in necessary to specify the implementation of the random choices of the .11 We will refer to a set of these as a realization of the initial condition. Rather than choosing completely unrelated random seed for each realization of the initial condition, a more well-defined and reproducible sample is obtained by drawing random numbers for each consecutive realization from a single random number stream. Thus the entire sample is represented by the initial seed of the random number stream, along with algorithmic details. Also this gives a definite ordering to the realizations, arising from the order of the random number stream. We do not require a large number of random numbers, and therefore opt for a simple linear congruential generator (LCG) pseudo-random number generator with a 31-bit seed/output value. The ending seed from one realization is simply used as the starting seed for the next. Candidate realizations that do not satisfy the property of having a negative perturbation on the symmetry axis are dropped and another is generated. These dropped candidates are not included in the numbering of the realizations used in the next section. We use the LCG discussed in section 16.1.3 of Newman & Barkema (1999). The initial seed for the set of realizations presented here was obtained from the standard Linux kernel random number source (/dev/random) and is 1866936915.

These sets of initial conditions provide a framework within which we may study a wide variety of possible systematic effects. This includes both physical systematics such as those explored in this study, or those arising from physical uncertainties in the numerical model. As study of the central ignition mechanism for SNe Ia advances, with improved flame models and DDT conditions, for example, it may be necessary to adjust the overall parameters of this framework, notably , in order to maintain a realistic sample. Also, as discussed previously, these choices can be compared to simulations of the early deflagration phase in order to both inform future choices and retrospectively understand the context of previously performed studies. One advantage of such a controlled initial condition is that it provides a slightly stronger probe of systematic effects than would generally be available observationally. This arises from the fact that it is possible to compare the same ignition sample rather than independent samples, though the latter can also be performed if desirable.

4.2. The Theoretical Population

The above development provides a clear definition of a population of ignition conditions from which we may draw a sample. Lacking further information about the nature of the harmonic content of the initial condition, we will assume cases drawn from this sample will have equal weight (likelihood). The population of supernovae that results is a purely theoretical construct. Any observed population will have a variety of progenitors that will have a distribution of intrinsic properties (composition structure, accretion history, etc). Therefore, while this population of ignition conditions will be very useful for studying systematic effects, caution should be exercised when drawing conclusions about observational ramifications. Future studies will address a more observationally motivated sample.

In order to understand the diversity in our sample, and the qualitative changes arising from changes in Ne fraction, the first 5 realizations from the sample sequence were run through the end of the detonation phase. The basic outcome properties of these ten cases are listed in table 4. The DDT time, , is defined as the time at which any part of the flame surface first passes a density of  g cm. This is the time at which the first detonation front is launched. The mass at high density at the DDT time, in this case that above a density of  g cm, will be discussed more below. The NSE mass, , is defined as integrated over the star. The Ni mass is determined by using eq. (29) below to estimate the local mass fraction from the and again integrating over the whole ejecta. The Si-group and O-Si masses are defined as and respectively and are discussed more in section 5, where the total nuclear energy released is also defined. All of these yields are determined after the detonation phase is complete, and the gross production of burning products is complete (see section 5). The amplitudes of the additive components that make up the initial conditions are given in Table 5, though this does not bear out any particular pattern in the results under discussion.

realization [s] energy released [ erg]
1 1.17 1.00 0.91 0.79 0.34 0.11 1.93
2 1.19 0.96 0.90 0.78 0.34 0.11 1.92
3 1.36 0.78 0.75 0.64 0.45 0.14 1.88
4 1.36 0.77 0.66 0.55 0.53 0.15 1.85
5 1.16 0.96 0.89 0.77 0.36 0.11 1.92
1 1.21 0.87 0.88 0.72 0.35 0.11 1.93
2 1.15 0.97 0.84 0.68 0.38 0.13 1.91
3 1.43 0.54 0.60 0.46 0.57 0.17 1.83
4 1.29 0.81 0.73 0.58 0.48 0.14 1.86
5 1.17 0.92 0.83 0.67 0.40 0.13 1.92

Table 4Outcomes for subset of theoretical sample
realization 13 14 15 16
1 -1.20 -1.16 0.49 0.21 -1.16
2 -0.73 1.17 0.95 -1.18 -0.62
3 -0.94 -0.64 0.69 0.14 -0.61
4 -1.15 -1.00 0.03 1.11 0.60
5 0.70 -1.69 -1.17 0.29 -1.05
Table 5 Initial condition amplitudes

Even this small sample from the theoretical population produces a diverse set of supernovae. The estimated yield of Ni spans a range between 0.45 and 0.8. The average is in the upper portion of this range. For the average is , with a sample standard deviation () of 0.11. For the average is , also with a sample standard deviation of 0.11. The decrease of the Ni yield with increasing Ne content is also directly reflected in the fraction of the NSE material that becomes Ni. This fraction is 0.86 with for , but is 0.80 with for . This is a direct result of the differences in initial neutron fraction due to the presence of Ne.

The main determining factor in the gross amounts of products synthesized in the explosions is the degree of expansion when the detonation ignites. This dependence has been discussed by previous authors (e.g. Röpke & Niemeyer, 2007; Townsley et al., 2007; Meakin et al., 2008), but the details of such a relation will differ among different delayed detonation scenarios. It is useful to characterize the degree of expansion by the mass at high density (above some density threshold), which forms an indicator of how much material will be processed to NSE. A part of this will become Ni and determine the overall brightness of the supernova. Several density thresholds for defining the mass at high density were considered. The mass above a density of  g cm, , appears the most appropriate. This conjecture is demonstrated by the open symbols in figure 4, which lie near a 1-1 relation (dashed line).

Figure 4.— Mass of all Fe-group (NSE, open symbols) elements and estimated mass of Ni (solid symbols) produced in the explosion plotted against the mass of material above a density of  g cm when the first detonation begins, an indicator of the degree of expansion of the star prior to the DDT. Two initial abundances of Ne are included, (squares) and (circles). Both of these yields are well-correlated with the mass at high density at the DDT transition, making the latter a useful intermediate indicator of the explosion outcome.

The variation in expansion at the DDT appears to arise from a competition between the rise of the highest plume and the expansion of the star in response to energy input. Both these processes have inherently similar timescales, the star’s dynamical time. Figure 5 compares the state of the star approximately 0.1 seconds after the launch of the first detonation for the 5 realizations at two Ne fractions run through the detonation phase. The top row displays the and the bottom , with each column being a separate realization of the ignition condition. The coloring indicate the nucleosynthetic yield, with black indicating NSE material, green Si-group material, and red O-Si mixture. It is immediately evident that the lowest Ni-yield case, realization 3 for , is also that which is the most expanded at the DDT transition.

Figure 5.— Comparison of burning products approximately seconds after first detonation is launched for different realization of the initial flame surface (columns) and for two abundances of Ne (rows). Simulations are performed in axisymmetry. Fuel and burning products are indicated by color: unburned C, O, Ne (yellow), O-Si (red), Si-group (green), Fe-group (NSE, black). Density in g cm is indicated by contours logarithmically spaced at integer powers of 10 starting from  g cm at the outside. One extra contour (red) is added at a density of  g cm.

The outcome of the deflagration appears to arise to some degree from the morphology of the deflagration, and therefore presumably from (randomly determined) characteristics of the ignition condition. Cases for which the highest plume is near the axis, realization 1, 2 and 5, expanded less before the detonation began than more equatorial lead plumes as in realization 3. Realization 4 appears to form an intermediate case, which in turns makes its outcome the most sensitive to the inclusion of Ne among these trials. This effect can be reasonably understood in terms of the interaction of the Rayleigh-Taylor (R-T) instability with the axisymmetric geometry. Plumes near the equator are dynamically more like rings instead of columns, therefore being more two-dimensional structures, and therefore are driven somewhat more weakly by the R-T instability. This effect may account for their slower rise in these simulations, and the resultant longer delay before detonation ignition. Additionally, equatorial plumes will generally burn material more quickly due to a larger integrated surface and volume in axisymmetry, and therefore could also enhance the expansion of the star directly.

This interaction with the geometry is important for extending these results to three dimensions. While the direct differential suppression of R-T will no longer be important, the presence or absence of localized ”spikes” in the initial condition is likely to become the most important determining factor in the competition between stellar expansion and plume rise, and therefore the explosion yield. Thus, assumed geometry is no longer important, but geometric features of the initial condition serve a similar role. This also creates the possibility that the specific morphology imparted to the flame surface by the convection field in the early deflagration phase is an essential aspect of the outcome of the DDT scenario. Filamentary or sheet-like structures might act quite differently from lumps of burned material once the strong R-T growth phase takes over.

5. A Deflagration Detonation Transition Supernova in two dimensions

In this section we describe the outcome of a typical case from the ensemble that will be used to study systematic effects, the case labeled ”realization 2” or ”r2” for . This case is considered typical because the mass at high density at the time of transition to detonation (see Section 4.2) and therefore the total mass of Ni synthesized, , is similar to the most probable values for both the ensemble of section 4, and, for Ni mass, the observed distribution of SNe Ia (e.g., Howell et al., 2009b).

5.1. The Explosion

The explosion can be roughly divided into 3 phases, the deflagration phase, the detonation phase, and the expansion phase. The deflagration phase lasts for the first 1.2 to 1.4 seconds of the simulation, depending on the particular ignition condition. During this phase a subsonic burning front, accelerated by turbulence and buoyancy, burns material in the inner portion of the star and expands it appreciably. As plumes of burning material rise from the core, the deflagration front eventually reaches low enough densities to transition to a detonation. The first of these to initiate a detonation will begin the detonation phase, which will burn the entire star within a few tenths of a second. Once all the material is burned, the star will continue to expand, beginning the expansion phase. The outgoing ejecta will eventually reach a state of free expansion in which the radial velocity of material is a linear function of radius and therefore the density becomes a simple function of time. We will now discuss each of these stages in turn for a typical case from our ensemble, including some details about how the simulation is accomplished.

Deflagration

Figure 6.— Time evolution of star and burning products during the explosion and establishment of the expansion. This sequence is for realization 2 with from the sample shown in Figure 5, which is a typical case based on exploration of our sample in section 4.2. By the end of this sequence ( s) all burning has ceased and the star is nearly in free expansion. Fuel and burning products are indicated by color: unburned C, O, Ne (yellow), O-Si (red), Si-group (green), Fe-group (NSE, black). Density in g cm is indicated by logarithmically spaced at integer powers of 10 starting from  g cm at the outside. One extra contour (red) is added at a density of  g cm. The final panel also includes contours of radial velocity at 5, 10, 15, and 20 km s (thick gray). The adaptive mesh is indicated in the early panels by outlines of blocks of cells, the smallest unit of contiguous refinement. The first detonation points are initiated at the moment depicted in the 5th panel at  s.

The initial burned region and the progression of the deflagration phase for realization 2 is shown in the top row of Figure 6. The burned material is colored black throughout this phase, because burning in the flame at these densities always results in Fe-group (NSE) material within the width of the artificial flame. Black lines show density contours on a logarithmic scale in integer powers of 10 beginning at 1 in the outer regions and reaching to 9 in the first panel. The initial WD has a central density of  g cm. The adaptive grid is shown by outline of the ”blocks.” These are logical mesh regions of cells that represent the smallest region that can be independently refined to the next level. Note also that PARAMESH restricts neighboring refinement regions to differ by at most a single refinement level (a factor of 2 in cell size). As discussed in section 3.3, the background star is refined to 16 km resolution, while energy-generating regions are refined to 4 km. The ”fluff” outside the star is not refined except as necessary to accommodate the interior grid regions.

Detonation and Unbinding

In the simulations, the transition of the burning to detonation must be initiated artificially because thermally activated burning is explicitly suppressed within the RD front (the artificially broadened flame). While section 3.1 describes a method for allowing thermal burning within the RD front for the purpose of allowing the detonation to propagate into material that has been ”partially burned,” the necessary estimate of temperature has so far proven to be of limited utility beyond detonation propagation. Spurious detonations typically occur if the thermal burning proximity detection is not utilized. Detonations are initiated by setting , the progress variable that represents the consumption of C, to 1 in one or several neighboring zones in one time step. Note that this involves no addition of energy to raise the temperature. The increased pressure resulting from the energy release following from the change in composition can then initiate a detonation that propagates away from the ignition point. It was found that for detonations at the density used here ( g cm), lighting a single  km cell does not always successfully launch a detonation. The outcome of small ignition points appeared to depend on the local flow characteristics, with fast-rising plumes being one of the more commonly problematic regions. This problem could be ameliorated by simply increasing the size of the lighted region. We found good success using a region with a linear radius of 8 km, resulting in the simultaneous ignition of a cell region, still small compared to the overall flow characteristics and scale for changes in density.

For this study we have chosen to characterize the DDT point by a simple density criteria,  g cm. Whenever the flame front, as represented in the simulation by the RD front, reaches this density, we light a detonation. Because the top of the deflagration is generally characterized by a few dominant plumes, the positioning of the detonation ignition points is relatively unambiguous. While quantitatively defines our detonation point, it is further necessary specify the method in which this density is used for placing the ignition point or points. In order to increase the chance of obtaining a robust detonation ignition, the ignition point is placed slightly outside of the rising burned region. Points are ignited when the RD front is approximately 64 km below the  g cm contour, at a point halfway between the RD front and the contour. Typically two points are ignited per plume, but sometimes three if the plume is wide, as is the case for the first ignited plume in realization 2 shown in Figure 6. This method places the detonation points at a density of between 1.05 and 1.1g cm. A density contour of g cm is shown in Figure 6 in order to show that the ambiguity in detonation time and place, introduced by the slight difference between and the actual point of detonation ignition, should be a fairly small factor in the outcome of the overall explosion, assuming it is applied consistently.

It should be emphasized that this treatment of detonation ignition is only a first, simplest option. The correct placement of the detonation ignition point is currently unknown, but at minimum a more realistic condition would take into account the local flow characteristics in an estimate of the Gibson scale and compare it to the flame width. It is even possible that the detonation ignition is qualitatively different. For example, we have ignited the detonation on the ”top” of the rising plumes, although the most pronounced mixing and tumble is occurring on the underside of plume caps. It may be that the detonation does not occur until this region reaches some density threshold, as mixing on the plume tops is less vigorous. This condition would lead to a systematic delay of the detonation ignition from the methodology implemented here. We leave evaluation of various alternatives for ignition of the detonation to future work.

Once the detonation is ignited, most of the remainder of the star is burned in a few tenths of a second, between 1.12 s and 1.6 s in realization 2. In the sequence shown in Figure 6, additional ignition points launch from where plumes in the lower half-plane reach . The refinement tracks the detonation at lower densities, but much of the higher density material remains refined longer because of the additional energy release as NSE material expands and alpha particles are recaptured. The display of the grid has been eliminated after the sixth panel due to it becoming too dense in this visualization. Nearly all of the material in the previously deflagrated region is burned, mostly to Fe-group material, demonstrating the success of the method used to allow the detonation to propagate into material within the RD front. Some material in heavily mixed regions is not fully burned, but given the low resolution achievable in the simulations, it is unclear how realistic this is. This incomplete burning does not appear to be a significant issue for gross yields from the explosion, but detailed nucleosynthesis, which will be undertaken in future work, will need to provide a better treatment of how deflagration ashes interact with the detonation front. This is discussed more with the final yields below.

Figure 7.— Dynamics of energy content during explosion. Shown are the total binding energy, , including the gravitational binding energy, the internal energy and the kinetic energy, the kinetic energy as a separate component, (dashed), the accumulated energy lost in the form of neutrinos, (dotted), the change in the total rest mass energy of the star due to nuclear processes, and the change in the total of all these energetic components . The term has been multiplied by 10 for display purposes.

The energetic history of the explosion is shown in Figure 7, where we find that, as expected, most of the energy is deposited during the detonation phase, and that the star is, in fact, bound up to this point. Shown is the total binding energy, , which includes the gravitational binding energy, the internal energy, and the kinetic energy. The exclusive source of energy is the nuclear energy release. The source of this energy is actually the change in rest mass energy of material due to nuclear processes. We quantify this by measuring the difference between the rest mass of all the material on the grid and the same number of baryons of symmetric matter () in a completely unbound state of free protons, neutrons and electrons:

(27)

where is the rest mass energy per baryons (approximately 1 g of baryons), , , and are the rest masses of the proton, neutron and electron respectively, and is the average nuclear binding energy as defined in section 3.1. The total rest mass energy is then , where the integral is over the computational domain. As discussed in section 3.1, we are using the convention that density as represented on the grid is actually the baryon density divided by , and is only approximately the density in g cm. Due to the binding energy of the unburned material, is in fact a large negative number at the beginning of the simulation. Therefore, Figure 7 only displays the change during the simulation, . It is immediately evident that this is the principal energy source, as its dependence is the direct converse of the energy deposition.

The final contribution to the total energy balance is the energy lost to neutrinos during the neutronization process. We define the cumulative neutrino loss by

(28)

where is the neutrino loss for the local conditions as tabulated by Seitenzahl et al. (2009) (see section 3.1). The energy loss due to neutrinos is also shown in Figure 7, where it has been multiplied by a factor of ten in order to be visible on this scale. We observe that there is no significant neutrino losses during the detonation phase, indicating that the neutronization takes place exclusively during the deflagration.

Adding all the energy terms together provides a useful check on the energy conservation of our code. The total is formed by . This sum again has a large offset due to the binding energy of the unburned C-O-Ne material, so we only display the change in total energy during the simulation. Very good energy conservation is observed, though there appears to be a small loss of energy (barely perceptible in Figure 7) during heavy refinement and derefinement, apparently related to . This loss may be a unexpected interaction of interpolation with the representation of the conserved quantities in the burning model. The loss is small enough to not be a concern, but will be addressed in the ongoing verification of our code components.

Transition to Expansion

In order to obtain an accurate evaluation of the distribution of ejected material in velocity, which is critical for spectral properties of the explosion, it is necessary to continue the simulation until the ejecta has reached approximate free expansion. The detonation has propagated throughout the interior of the star by just after  s as can be seen from Figure 6. Some nucleosynthesis takes place later as the NSE freezes out, but the need for fine refinement of the energy generating regions ends at this time since there are no propagating burning fronts. As a result, at 1.8 s in this simulation, the refinement of energy-generating regions, including regions in which the propagation of RD front is still taking place, is disabled. Also the refinement of non-energy-generating regions is coarsened from 16 to 32 km. Ideally, de-refinement with further expansion could take place automatically with appropriate detection of the expansion of flow features. This somewhat delicate balance was not attempted in these calculations, and nearly the entire ejecta remains refined at 32 km through the end of the calculation.

Figure 8.— Run of radius and mass with radial velocity at s, demonstrating that free expansion has been approximately attained.

The calculation was halted as material began to flow off of the grid around  s. The domain was chosen to be of such a size that this time was late enough for the ejecta to be in approximate free expansion. We take free expansion to be indicated by a linear relation between the radial velocity and the radius. Profiles of the radius against the radial velocity of the outgoing ejecta are shown in Figure 8. Shown are both line-out profiles at , 90, and 135, as well as the averaged profile (offset by  Mm s for clarity). The average here is the mass weighted average radius for material moving at a radial velocity within a bin of 0.2 km s. While the average profile has a fairly linear relation, there is some departure from linear in the outer regions for certain directions in the ejecta. Additionally, the radial velocity component completely dominates any motion at this time. Radial velocity contours are also shown on the final panel of Figure 6 and are fairly symmetric, though there is a noticeable degree of asymmetry in the abundance profiles.

5.2. Explosion Yield and Ejecta Structure

Figure 9.— Masses of synthesized material in time during explosion, see text for definitions. Note that these are estimates of final outcome, and therefore the relaxation of the NSE material dominantly Fe-group constituents is not evident in this plot.

The time history of the material produced in the incineration of the star is shown in Figure 9. In the interest of simplicity, we leave tracer particle post-processing and detailed nucleosynthesis to future work and instead just make use of the progress variables defined in the burning model (see Section 3.1). The ashes of the first stage of burning are predominantly O and Si, with some other intermediate mass elements (Ne, Mg), and its local mass fraction is given by . Further burning produces a mixture of various Si-group nuclides, predominantly Si, which begins in NSQE. We are concerned with the form that material will take in the outgoing ejecta, after NSQE and NSE material have completely relaxed. As such, Figure 9 shows the final yield masses as simply Si- and Fe-group. The mass fraction of the Si-group material is given by . Material in NSE is given by , and will consist, in the ejecta, of almost entirely Fe-group elements. These definitions are also used below to characterize the abundance profile of the ejecta and allow the colors displayed in figure 6 to be interpreted roughly in terms of the nucleosynthetic products.

By tracking , we also have a local measure of the neutronization that has occurred, and from this we can estimate how much of the NSE material will be in the form of stable, neutron-rich Fe-group nuclides instead of Ni. For , the ejected material is nearly pure Ni (Meakin et al., 2008). In order to estimate the local mass fraction of Ni of a neutronized fluid element being ejected, we assume that the next most abundant nuclides are an equal admixture of Fe and Ni, as is observed in the tabulation of nucleosynthetic outcome with by Meakin et al. (2008). The local mass fraction of Ni is then estimated by

(29)

This estimated mass fraction is also used below in discussing abundance profiles.

From Figure 9, nearly all of the intermediate mass elements are produced during the detonation phase, with a small amount of Si group material produced in the deflagration. This latter is likely to be greater for more expanded explosions with lower NSE and Ni yields. The isolation of the neutronization to the deflagration phase is also evident. While only half of the NSE material produced in the deflagration phase will be ejected as Ni, nearly all of that produced in the detonation phase will end up as Ni. Note that since for this explosion, the maximum even for non-neutronized material is 0.95. Although it appears all the burning is complete by approximately 1.6 seconds, this is only true in terms of gross yields. At that time much of the Fe-group material will still be in NSE, and will relax as the star expands. The small glitch in the Ni curve at 1.8 seconds is due to the cell-mixing upon zone de-refinement and the resulting change in our estimate of the final Ni yield via equation (29).

Figure 10.— The electron molar fraction at  s after free expansion has been essentially attained. Contours indicated radial velocities of 5, 10, 15 and  km s. View is the same as the last panel of Figure 6.

The elemental ejecta structure is shown by the combination of the last panel of Figure 6, which shows radial velocity contours overlayed on the color-coded gross yields as defined above, along with Figure 10 that gives the spatial distribution of at the same moment with the same radial velocity contours shown. Material below is expected to have fairly little Ni, favoring more neutron-rich nuclides instead.

Figure 11.— Run of burning products with radial velocity at s.

For further detailed scrutiny, we have extracted the abundance pattern along radial lines at several latitudes and the averaged abundance profile for 0.2 km s bins, all shown in Figure 11 in radial velocity. In this figure the Ni mass fraction is estimated from equation (29). In the case of the averaged profile this is done before averaging. The distribution of the mass in radial velocity is shown in Figure 8, showing that very little mass is outside of approximately 15,000 km s.

Figure 12.— Distribution of burning products with mass coordinate in ejecta.

Finally the distribution of material in mass is shown in Figure 12, where the mass coordinate is defined as the mass enclosed by consecutive radial velocity shells.

5.3. Resolution Dependence

Given the small scale of the expected flame surface structure with respect to the grid scale in these simulations (4 km), and the very basic treatment of turbulent acceleration of the burning, some dependence on resolution is expected. The important question is whether or not such resolution effects are strong enough to preclude the use of suites of 4 km simulations from being used to evaluate systematic effects. We find that the apparent level of uncertainty from the low resolution is acceptable when considering particular metrics, but that there are important dependencies on resolution that we hope can be resolved as treatments of subgrid flame structure and their usage becomes more advanced.

The metric used in our study of systematics in section 6 below is the mass at high density at the DDT time. This time is defined by the moment at which the lowest density of fresh fuel being encountered by the flame, , falls below , which we are assuming here is  g cm. It was found in section 4.2 that a good density threshold to use to define the mass at high density is  g cm. Thus it is useful to look at how the evolution of as the star expands, and the flame rises, depends on resolution. The top panel of Figure 13 displays the dependence of on . Evolution in this figure moves from right to left as the flame rises through the star. Two different realizations from the sample developed in section 4.2 are shown. One is that described in detail in this section (realization 2), and the other is a more extreme case that produces much less NSE material (realization 3).

Figure 13.— Investigation of the dependence of results on resolution. The top panel shows the mass above the density threshold  g cm, found in section 4.2 to correlate well with total NSE yield. The curves are plotted against the minimum density of fuel being currently burned by the flame, , such that the evolution follows the curve from right to left as the flame rises and the star expands. The bottom panel shows the total mass burned, essentially all of which is NSE material. While the value of when appears fairly robust with resolution, providing some confidence in its usage to study systematic effects in section 4.1, the total burned mass shows more resolution dependence.

In addition to several resolutions, a case was run with the value of the parameter that appears in the buoyancy-compensation prescription for the flame speed (see section 3.2) doubled. By simultaneously doubling this parameter and halving the resolution, the actual value of the prescribed flame speed stays the same.

Our metric for measuring the expansion of the star, , shows a fairly low sensitivity to resolution for these cases, especially for realization 2. There is also fairly little sensitivity to the value of the parameter . This insensitivity is likely due to the fact that the two competing processes at work, the pulsational expansion of the star and the rise of burned plumes through the star, are not overly sensitive to the poorly resolved flame surface. The small-scale flame structure is not important for the selection of the dominant plume because this is largely determined by the initial condition (see also discussion in section 4.2). Also, due to the delay in manifestation of the turbulence from the Rayleigh-Taylor instability of burned plumes, the additional flame surface area and concomitant energy deposition may come after the critical phase for the launch of the pulsational expansion of the star. This delay in energy release might lessen the impact of uncertainty in the burning rate during the intermediate stages of the deflagration, the time when it is most difficult to model accurately. However, this situation highlights the necessity for a careful understanding of the very early portion of the deflagration phase and the interaction of the flame surface with the turbulent convection field arising from the simmering phase.

There is a stronger and more systematic dependence of the total burned mass on resolution, as shown by the bottom panel of Figure 13. However, the scatter in the consecutive resolution cases performed makes them difficult to interpret clearly. In an observational sense, the resolution dependence might be closely related to the weakness of the neutron-enriched core observed in the yields from realization 2, as discussed below. The additional enhancement of burning due to turbulence should be most effective in the inner regions of the star where the intersecting wakes of rising plumes will lead to strong turbulent shearing of the flame surface. This boost may burn core material early enough in the expansion to enable this material to undergo significant electron captures. The two-dimensional simulations performed here (in which there is no physically realistic turbulence cascade) and the lack of an explicit treatment of turbulence-flame interaction make it inappropriate for this to be addressed in the current study. We will simply leave the degree to which turbulence influences the neutronization of material an open question, which is not expected to depend on the Ne content under study here due to the predominance of turbulence in setting the fuel consumption rate under these conditions.

5.4. Brief Comparison with Observed Properties of SN Ia

As seen from Figures 6, 10, 11, and 12, our explosion reproduces fairly well the observed radial velocity structure of the intermediate mass elements ejected in the SNe Ia (Filippenko, 1997). We will leave a fully detailed comparison until full abundances can be obtained from post-processing tracer particles, but many important features warrant mentioning here. There is a prominence of Si in abundance for the velocity range of roughly 10 to 15 Mm s, accounting for nearly of the final ejecta, with some Si-group material extending down as far as 8 Mm s. This is broadly consistent with the spectral evolution of SNe Ia, though possibly 1-2 Mm s high in velocity compared to the spectral evolution of the most standard cases like 1994D (Patat et al., 1996). Since the mass distribution (Figure 8) is weighted toward lower velocities, conclusive comparison will require radiative transfer calculations. Also, this is only a single case of many possible, which may have Si at lower velocities. The ejecta is fairly symmetric, with a slight asymmetry arising from the location and timing of the detonation sites. While the degree of asymmetry is likely not too unrealistic, the pattern imprinted may be a bit different in a three-dimensional treatment. It does appear that the velocity extent of Si features near the photosphere would be viewing-angle-dependent. Note that the jet-like columns of Fe-group material on the symmetry axes are unrealistic, but appear to have a fairly small contribution to the averaged profile shown in the last panel of Figure 11. Detailed radiative transfer calculations will be required to determine if the patchiness of the border between the Si and Fe-group material provides a good match to polarimetric observations (e.g. Wang et al., 2003).

The distribution of neutron-rich material demonstrated in Figures 10 and 11 compares less well with the distribution inferred from observations (Stehle et al., 2005; Mazzali et al., 2008). On the positive side, the impurity created by the distribution of deflagration material at velocities around 3-8 Mm s dilutes the pure Ni to a degree that provides a reasonable agreement with observations. However, notably absent is the predominance of stable Fe at low velocities,  Mm s, which is inferred from observations of 2002bo (Stehle et al., 2005) and 2004eo (Mazzali et al., 2008). While the treatment of the ejecta evolution leading to this inferred abundance is approximate, it seems unlikely that such a stark difference could be reconciled without a change in the model. The most apparent reason for this shortcoming of the model is the flame treatment being utilized here. The buoyancy-compensation form of Khokhlov (1995), which we are using, does not account for interaction of the turbulent wakes with other burning fronts. Thus as plumes rise out of the core, the turbulence they leave behind does not lead to the enhancement of flame spreading that it should. This shortcoming will be addressed in future work, with more sophisticated treatments of the turbulence-flame interaction (Colin et al., 2000; Schmidt et al., 2006). Investigation of this effect should be undertaken with attention to the detailed nuclear products, as extensive central burning early in the deflagration phase can lead to ”over”-neutronization (e.g. Brachwitz et al., 2000).

Many SNe Ia have been observed with high-velocity Ca features (Mazzali et al., 2005). This material is generally produced during the transition from NSQE to NSE (e.g. Nomoto et al., 1984) and thus should appear at the border between Si-group (green) and NSE (black) regions in Figure 6. In this simulation, these yields are restricted to velocities  km s, where the main component also lies in observed supernovae. Some features, however, do indicate that this material might be produced at higher velocities by intersecting detonation fronts. The case of the symmetry axes, while an artifact here, suggests that points where 3 detonations spreading through the star meet might be good candidate sites for the formation of such high-velocity features. This is reinforced by features seen in some detonation interactions in other realizations (see the comparison of ignition conditions in Figure 5). Realistic tests will await DDT calculations in three dimensions. Additionally, we have assumed that detonations are ignited when the first flame front passes through the DDT density. This may be too stringent a condition, and some number of deflagration plumes may extend well above this transition region (which also corresponds closely to the division between Si- and NSE-rich ashes) before the detonation takes place.

6. Systematic Effects of Ne Arising from Dynamics

We have implemented the current study as an attempt to determine whether the presence of Ne has a significant systematic effect on the expansion state of the star at the time of first detonation. In light of the relation between the expansion and final yields discussed in Section 4.2, this is a likely candidate for an effect of Ne that can compete with that expected due to the simple overall neutron excess highlighted by Timmes et al. (2003). Our focus on the degree of expansion, quantified by , in such specificity is motivated by a desire to separate the contribution of Ne to the overall nucleosynthesis into a dynamical component distinguishable from the gross neutron excess. The use of the outcome of the deflagration phase alone also allows many more cases to be run at the same computational cost. This is very important given the high variance in the outcome of the supernovae. Averaging over many realizations is necessary to obtain uncertainties that are small compared to the expected effect of neutron excess. Additionally, targeting our study so specifically gives more power for tertiary conclusions about the impact of other effects that might modify the propagation speed of the burning front.

Variation due to Ne was explored by simulating the deflagration phase for 20 realizations with , 0.02, and 0.04. As discussed in section 2, the carbon abundance was kept fixed at by mass as was the central density of the initial WD. This results in the WD being slightly less massive with a higher Ne content. Our WD, as mapped on the initial grid, is , or  g for , 0.02, and 0.04 respectively. Samples with different Ne content are compared in Figure 14, where for the Ne enhanced cases is plotted against the same outcome for no Ne. The dashed line indicates a 1-1 relation. Systematic effects should arise as a departure from this relation on average. The two nearly coincident points with error bars indicate the averages of the samples with different Ne abundances. The outer error bars indicate the standard deviation and the inner the standard deviation of the mean. The sample with has a slightly higher average of 0.851 compared to for , and both of these are somewhat lower than the for no . Standard deviations of the mean are 0.042 for both of the nonzero cases and 0.046 for the zero .

Figure 14.— Relation of degree of pre-expansion for cases with (circles) and (squares), to a case with no Ne. Degree of pre-expansion is indicated by the mass above a density of  g cm when the first detonation initiates, which is an indicator of how much Fe-group material will be ejected in the supernova. The point with vertical and horizontal error bars indicates the average, sample standard deviation (outer error bars) and standard deviation of the mean (inner error bars). The sample with has a slightly higher average. There is not a significant difference between the average expansion in the three cases.

The overlap of the inner error bars with the 1-1 relation indicates that there is no statistically significant impact of Ne content on the expansion prior to DDT over this range of mass fractions. On the basis of neutron excess (Timmes et al., 2003), the decrease in the