CARBON DEFLAGRATION IN TYPE Ia SUPERNOVA. I. CENTRALLY IGNITED MODELS

Common Type Ia supernovae (SNe Ia) are believed to be the thermonuclear explosions of carbon–oxygen white dwarfs in binary systems (Hoyle & Fowler 1960; Hillebrandt & Niemeyer 2000). Beyond this generality though lies room for considerable diversity. The white dwarf can grow to criticality by accreting mass slowly from a companion star, the so called "single-degenerate" model, or by merging with another white dwarf, the "double degenerate model." Depending upon the accretion rate and angular momentum transfer, the mass of the white dwarf when the nuclear runaway begins can vary upward (Howell et al. 2006) or downward (e.g., Woosley & Kasen 2011) from the Chandrasekhar value (1.38 M_☉ for ignition). Ignition itself may happen in the center, slightly off-center, or near the surface of the white dwarf. Historically, greatest attention has been given to the centrally ignited Chandrasekhar mass model (the MCh model). This focus partially reflected the predominance, until recently, of one-dimensional (1D; hence spherically symmetric) codes and models, but also the success of the model itself. Given the liberty of varying uncertain flame physics, including the possibility of a delayed detonation, the MCh model has demonstrated striking agreement with the observed energetics, light curves, spectra, and nucleosynthesis of common SN Ia events (Nomoto et al. 1984; Branch et al. 1985; Höflich & Khokhlov 1996; Mazzali et al. 2007).

More recently, though, it has been realized that other models are probably needed, especially if one is to explain unusual events (e.g., Howell et al. 2006; Stritzinger et al. 2006; Perets et al. 2010; Pakmor et al. 2011). Simulations with increased complexity and realism have also shown that the MCh model itself probably does not ignite exactly in its center (Garcia-Senz & Woosley 1995; Höflich & Stein 2002; Kuhlen et al. 2006; Zingale et al. 2011; Nonaka et al. 2012). This symmetry breaking may have important observational consequences for the explosion mechanism (Plewa et al. 2004; Röpke et al. 2007b), as well as the spectrum and light curve (Kasen et al. 2009). Still it remains interesting to explore the traditional model. Despite X-ray limits on the progenitor system (Gilfanov & Bogdán 2010), the single-degenerate model remains viable (Voss & Nelemans 2008; Hachisu et al. 2010), and even in double degenerate systems, the mergers of white dwarf pairs may produce MCh models as well as super-MCh models (Yoon et al. 2007). Studies of off-center ignition also show a small but non-negligible probability of ignition occurring so close to the center that it probably burns through to the other side (Zingale & Dursi 2007; Nonaka et al. 2012).

In the last decade, the symmetric MCh model has also attracted considerable attention from the three-dimensional (3D) simulators. Despite some initial optimism (Reinecke et al. 2002; Röpke et al. 2006), these calculations (Gamezo et al. 2003; García-Senz & Bravo 2005; Schmidt & Niemeyer 2006; Röpke et al. 2007a) have shown that even though pure deflagration, i.e., deflagration without a subsequent detonation, released enough energy to unbind the star, it is unable to provide the larger kinetic energies characteristic of common SN Ia. Equally problematic, the deflagration-only models produce too little ⁵⁶Ni, leave behind significant unburned carbon and oxygen with low speed near the center of the ejecta, and lack intermediate mass elements (IMEs) at high velocity. These aspects make it difficult for pure deflagration alone to produce a common SN Ia consistent with observations, although see Kromer et al. (2013) for pure deflagration models of sub-luminous SNe Ia.

Here we revisit this classic problem of centrally ignited deflagration in MCh models, not with the intention to prove that they work, but to better understand why they do not work. Our study has several motivations. First, different groups have assumed different initial conditions, codes, grid geometries, and resolutions. With the power of adaptive mesh refinement (AMR) and larger computers, we are able to define the initial conditions more precisely and track the burning more realistically. Previous studies also disagree concerning how much carbon is left at the center after the explosion is over and on the amount of IMEs produced. We want to understand this and give an alternative set of results. We are also implementing a novel, inexpensive approach to flame propagation that uses tables of yields generated off-line assuming isobaric burning and want to test the sensitivity of results to uncertain parameters in a well-studied case. Finally, in order to study delayed detonation in the MCh model as we plan to do in a subsequent paper, it is necessary to better understand the conditions that exist at the transition.

In Section 2 we briefly describe the fully compressible hydrodynamics code, CASTRO used in these studies (see also Almgren et al. 2010), and describe our treatment of the flame physics. Details of the calculations are presented in Section 3, followed by the results of our simulations in Section 4. A discussion and comparison to previous work are presented in Section 5.

2.1. CASTRO

2.2. Mapping and Initialization

The presupernova model employed here is a 1.38 M_☉ white dwarf with an initial composition of 50% carbon and 50% oxygen, generated using KEPLER, a 1D implicit Lagrangian code (Weaver et al. 1978; Woosley et al. 2002). This white dwarf model has been previously used many times as a starting point for studies of SN Ia (e.g., Niemeyer & Woosley 1997; Woosley et al. 2007) and is similar to the initial model of Nomoto's model "W7" (Nomoto et al. 1984). The central density of the star at the time of the final runaway was 2.6 × 10⁹ g cm⁻³, the central temperature, 6 × 10⁸ K, and the initial radius, about 1700 km. To ensure hydrostatic equilibrium on the new Cartesian grid, the distribution of physical parameters was adjusted from their values in KEPLER to satisfy the differencing formulae, 0.5(ρ_i + ρ_{i − 1})g = (P_i − P_{i − 1})/dr. Constant entropy was assumed in the convection zone which includes most of the mass of the presupernova star. Here g = GM(r)/r² was defined on the interface between zones, where M(r) is the enclosed mass at the radius of r. For the multi-dimensional calculations, a very fine grid was first generated in hydrostatic equilibrium in the 1D version of CASTRO using spherical coordinates. This grid was then interpolated onto the Cartesian grid AMR hierarchy to give the starting model in 3D. The Helmholtz stellar EOS of Frank Timmes was used in all of our calculations (Timmes & Swesty 2000). Because the Eulerian grid requires non-zero values in all zones, the white dwarf was surrounded by an artificial, but inert, low-density medium with temperature 5 × 10⁶ K and density 10⁻⁴ g cm⁻³. This low density medium had no dynamic effect on the supernova explosion calculation. A numerical "sponge" was implemented to damp the velocity in the surrounding medium where densities were between 10⁻⁴ and 1000 g cm⁻³. This was done by multiplying the three velocity components by a factor that varied smoothly from nearly 0 to 1 in this density range (Zingale et al. 2009). The sole purpose of the sponge was to keep the time step in the low density, artificial circumstellar medium from affecting the calculation as would happen if the medium began to accrete onto the white dwarf.

Burning was initiated by a nearly spherical bubble of hot iron ash surrounding the center of the star with a radius of 100 km. The temperature in the ash was set to 8.5 billion K, and binding energy per nucleon, BE, and mean atomic weight, $\bar{A}$, were set using the nuclear burning tables (Tables 1 and 2) to 8.17 MeV nucleon⁻¹ and 11.6, respectively. The ash density was derived from the EOS using this temperature and the same pressure as the surrounding fuel. Since the initial ash was in nuclear statistical equilibrium (NSE), its abundances and Y_e started to evolve immediately. The surface of this sphere was perturbed using spherical harmonic functions assuming a superposition of all l-wave modes from 14 to 20 and a maximum perturbation amplitude of about 25 km. The choice of 100 km for the radius was influenced by studies of ignition (e.g., Nonaka et al. 2012) which showed that prompt ignition outside this radius was very unlikely. This is a smaller initial radius than has been employed by some groups (Röpke et al. 2006, 2007a). The initial zoning and AMR structure employed are described in Section 3.1. Pressure was kept unchanged across the ash fuel interface, though the ash had lower density and higher temperature. The initial perturbed sphere was therefore Rayleigh–Taylor (RT) unstable, though it took some time for these initial instabilities to grow. This use of a single large deformed ignition region seems more physical to us than ignition at many numerous distinct points (Röpke et al. 2006), especially since recent studies show that the flame only ignites once (Nonaka et al. 2012).

Table 1. Final Composition—Deflagration Ashes

Log Density	He	C	O	Ne	Mg	Si	Fe
(g cm−3)
6.40	...	0.056	0.341	0.411	0.177	0.015	...
6.50	...	0.050	0.348	0.395	0.189	0.019	...
6.60	...	0.037	0.358	0.373	0.208	0.024	...
6.70	...	0.014	0.399	0.285	0.255	0.046	...
6.80	...	0.002	0.549	0.012	0.244	0.194	...
6.90	...	...	0.560	...	0.213	0.227	...
7.00	...	...	0.561	...	0.194	0.245	...
7.10	...	...	0.548	...	0.167	0.285	...
7.20	...	...	0.382	...	0.031	0.587	...
7.30	...	...	...	...	...	0.992	0.007
7.40	0.001	...	...	...	...	0.971	0.028
7.50	0.003	...	...	...	...	0.825	0.171
7.60	0.013	...	...	...	...	0.370	0.617
7.70	0.030	...	...	...	...	0.079	0.891
7.80	0.048	...	...	...	...	0.033	0.918
7.90	0.064	...	...	...	...	0.037	0.899
8.00	0.081	...	...	...	...	0.040	0.879
8.10	0.097	...	...	...	...	0.044	0.860
8.20	0.112	...	...	...	...	0.047	0.842
8.30	0.125	...	...	...	...	0.050	0.825
8.40	0.137	...	...	...	...	0.054	0.809
8.50	0.148	...	...	...	...	0.057	0.795
8.60	0.157	...	...	...	...	0.061	0.781
8.70	0.166	...	...	...	...	0.066	0.769
8.80	0.173	...	...	...	...	0.070	0.757
8.90	0.179	...	...	...	...	0.075	0.746
9.00	0.185	...	...	...	...	0.080	0.735
9.10	0.190	...	...	...	...	0.085	0.725
9.20	0.194	...	...	...	...	0.091	0.715
9.30	0.198	...	...	...	0.001	0.097	0.705
9.40	0.201	...	...	...	0.001	0.103	0.695
9.50	0.203	...	...	...	0.001	0.110	0.685
9.60	0.205	...	...	...	0.001	0.118	0.675
9.70	0.207	...	0.001	...	0.001	0.126	0.665
9.80	0.208	...	0.001	...	0.001	0.135	0.655
9.90	0.208	0.001	0.001	...	0.002	0.144	0.644
10.0	0.207	0.001	0.001	...	0.002	0.156	0.639

Download table as:  ASCIITypeset image

Table 2. Bulk Properties—Deflagration Ashes

Log ρ	T9f	ρ7f	dρ7/dX12	BE/A	$\bar{A}$
(g cm−3)	(109 K)	(107 g cm−3)	(107 g cm−3)	(MeV nucleon−1)
6.40	1.790	0.122	0.258	8.040	18.340
6.50	1.903	0.158	0.316	8.046	18.430
6.60	2.036	0.203	0.390	8.057	18.630
6.70	2.216	0.255	0.493	8.082	18.980
6.80	2.471	0.305	0.651	8.137	19.190
6.90	2.631	0.399	0.792	8.144	19.240
7.00	2.791	0.522	0.955	8.147	19.280
7.10	2.978	0.677	1.163	8.158	19.460
7.20	3.382	0.794	1.582	8.262	21.680
7.30	3.899	0.878	2.235	8.464	29.270
7.40	4.154	1.158	2.708	8.467	29.440
7.50	4.444	1.507	3.310	8.488	30.480
7.60	4.781	1.932	4.099	8.527	31.730
7.70	5.071	2.533	4.959	8.518	28.210
7.80	5.311	3.375	5.869	8.479	23.970
7.90	5.542	4.489	6.909	8.442	21.130
8.00	5.761	5.948	8.104	8.406	19.000
8.10	5.969	7.848	9.482	8.373	17.410
8.20	6.169	10.31	11.08	8.343	16.190
8.30	6.362	13.49	12.93	8.316	15.240
8.40	6.550	17.58	15.08	8.293	14.500
8.50	6.736	22.83	17.60	8.272	13.910
8.60	6.921	29.54	20.54	8.254	13.440
8.70	7.105	38.13	23.98	8.238	13.050
8.80	7.291	49.09	28.02	8.224	12.720
8.90	7.479	63.05	32.76	8.211	12.450
9.00	7.670	80.84	38.33	8.200	12.220
9.10	7.866	103.50	44.87	8.190	12.020
9.20	8.065	132.20	52.57	8.181	11.860
9.30	8.270	168.70	61.63	8.173	11.720
9.40	8.480	215.00	72.30	8.166	11.600
9.50	8.697	273.80	84.88	8.160	11.500
9.60	8.922	348.30	99.71	8.154	11.420
9.70	9.153	442.60	117.20	8.149	11.360
9.80	9.393	562.00	137.90	8.145	11.310
9.90	9.641	713.20	162.30	8.141	11.270
10.0	10.090	902.10	195.80	8.194	11.260

Download table as:  ASCIITypeset image

2.3. Nuclear Physics of the Burning

The nuclear physics of high temperature carbon deflagration can be quite complex. The products of the burning depend upon the density and may or may not be in NSE. After the deflagration front has passed, the composition continues to evolve through electron capture and the reassembly of partially photodisintegrated matter. At low density, nuclear burning can continue. Electron capture changes the electron mole number, Y_e, to which the pressure is quite sensitive. Reassembly of alpha-particles into iron-group elements provides a substantial energy boost as the composition cools by expansion. Nuclear fusion can appreciably alter the final composition, especially at lower temperatures where NSE is not achieved. These effects can only be followed accurately with a large nuclear reaction network. Yet, at the same time, coupling a large network directly to a 3D simulation would be very computationally expensive. Worse still, the nuclear time step can be very small, especially in NSE, and the solution to the stiff system of Ordinary Differential Equations (ODEs) is prone to numerical instability.

Our solution was to use tables to describe the energy generation and composition changes. The tables were generated off-line using a large nuclear reaction network whose accuracy has been tested many times in the Kepler stellar evolution code. Different tables were used to describe the burning in the flame Section 2.3.1 and the volumetric burning in NSE Section 2.3.2. The implementation of this burning in a flame model is discussed in Section 2.4.

2.3.1. Burning in the Flame

The first stage of the burning occurs as the flame passes through the unburned carbon and oxygen. At the high density of the deep white dwarf interior, this flame is very thin, though highly convoluted, essentially a moving unresolved discontinuity that turns carbon and oxygen into iron-group elements, alpha-particles, and nucleons in NSE. At lower densities, the heat capacity of the matter rises and the burning results in a lower temperature. Carbon, oxygen, neon and silicon burning regions in the flame grow thicker and eventually burning no longer culminates in NSE, but in the production of intermediate mass elements, silicon through calcium. At still lower density, the oxygen does not burn to completion and one has neon, oxygen and magnesium. Eventually the burning goes out.

A characteristic shared at all densities though is that the burning in the flame transpires at constant pressure. Detonation does not occur in the models being studied, and because the flame speed during the deflagration stage is such a small fraction of the sound speed (∼1%), pressure equilibrates across the flame. Since the gas is initially very degenerate and has Y_e = 0.50, the pressure in the fuel, and hence in the ash, is completely determined by the fuel's density and is not very sensitive to its (low) temperature. Assuming a reasonable limit to the timescale for the burning to continue once the energy generation saturates at a maximum, the change in composition at a given density can be uniquely determined from the initial composition and the isobaric condition. This allows the energy yield and ash composition to be computed off-line using a reaction network of arbitrary size.

Beginning at a low temperature, taken here to be 6 × 10⁸ K, and for a grid of initial densities from 10^6.4 to 10^10.0 g cm⁻³ in steps of log ρ of 0.05, the results of nuclear burning were calculated using a nuclear reaction network that contained 188 isotopes, including neutron, proton, ⁴He, ^12–14C, ^14–15N, ^16–18O, ^18–19F, ^19–22Ne, ^22–23Na, ^23–28Mg, ^26–29Al, ^27–32Si, ^31–36P, ^30–33S, ^34–39Cl, ^35–42Ar, ^38–42K, ^39–48Ca, ^41–49Sc, ^43–52Ti, ^45–53V, ^47–56Cr, ^49–57Mn, ^51–60Fe, ^53–61Co, ^55–64Ni, ^57–65Cu, ^59–68Zn, ^61–69Ga and ^63–72Ge. The initial composition was 50% ¹²C and 50% ¹⁶O by mass. Burning proceeded with an imposed isobaric condition. That is, as time passed the temperature rose and density fell in response to the energy generation, but the changes in each were such as to keep the pressure constant for the appropriate EOS. The "Helmholtz" EOS was used to obtain the pressure, its derivatives, and the heat capacity at constant pressure during this iteration (Timmes & Swesty 2000).

For each point in the ρ–T grid, after a lengthy time, which would be much shorter in the star due to conduction, the temperature gradually began to rise and the rise accelerated along with the nuclear energy generation. For each density, there came a point where the nuclear energy generation reached a maximum and then began to decline due to the depletion of the most incendiary parts of the composition. Once the energy generation had declined to 80% of that maximum value, a counter was set and 0.1 s later the calculation was halted and the characteristics of the ash sampled. The 0.1 s is arbitrary, but the composition and temperature were not changing much during this time. The value chosen is a lower bound to the hydrodynamic expansion time of the star and an upper bound to the time for burning to cross the typical five zone thickness of the (thickened) flame.

This procedure is an approximation to what would actually happen in a multi-zone calculation of a steady state laminar flame. In a flame, conduction, not gradual nuclear burning, would heat the material to the point that it flashed on a timescale sufficiently short to sustain a self-similar front with a thickness given by the Landau condition (Landau & Lifshitz 1959), which says that the nuclear time balances the energy diffusion time. Later in the burning, a small amount of heat would diffuse out of the ash into the fuel in order to ignite it and keep the flame going. However, because of the extreme temperature sensitivity of the burning rates and the small heat capacity of the degenerate matter, the amount of heat that must diffuse to propagate the flame is small. The major part of the burning thus occurs in thermal isolation at constant pressure and can be approximated by a single zone as we have done here. It is necessary only to resolve the energy generation spike that is most sensitive to the density of the fuel, and the small amount of burning that occurs after that spike. If the flame were very much broader or slower moving and the star did not expand so that the 0.1 s timescale was orders of magnitude longer, the ashes would be different, but the necessary change in timescale to have a major effect is very large. In the future, if flames need to be studied in very different circumstances, it would not be difficult to couple the off-line reaction network to a multi-zone (1D) flame model (Woosley et al. 2009).

The validity of the approximation is demonstrated by the good agreement with previous fine-zoned studies of laminar flames under the same conditions. Table 3 shows the comparison with other studies that also used either a 50% C, 50% O composition (Timmes & Woosley 1992; Woosley 2007) or a 40% C, 60% O composition (indicated with an asterisk; Woosley et al. 2011). Overall the agreement between ash temperature and density changes found in these studies is excellent. For the limited compositional information given in these references, the agreement was also very good.

Table 3. Comparison with Published Flame Models

Log Density	Tash	Published Tash	Δρ/ρ	Published Δρ/ρ	Reference
(g cm−3)	(109 K)	(109 K)
6.78	2.42	2.45	0.514	0.527	Woosley (2007)
7.0	2.79	2.79	0.478	0.477	Woosley (2007)
7.398*	4.11	4.1	0.529	...	Woosley et al. (2011)
7.477*	4.33	4.3	0.515	...	Woosley et al. (2011)
8.0	5.76	...	0.405	0.426	Timmes & Woosley (1992)
9.0	7.67	...	0.192	0.192	Timmes & Woosley (1992)
10.0	10.1	...	0.098	0.094	Timmes & Woosley (1992)

Download table as:  ASCIITypeset image

The resulting composition from our network simulations was packed into the seven abundance variables carried by the 3D code depending on the atomic weight of the product nucleus. Species lighter than A = 12 were packed into helium; species from 12 through 15 into carbon; 16 through 19 into oxygen; 20 through 23 into neon; 24 through 27 into magnesium; and 28 through 47 into silicon. Hence all IMEs, silicon through calcium are contained in the "silicon" variable. Species heavier than A = 47 were generically considered to be "iron" (Table 1). Additionally the final temperature, density, binding energy per nucleon, and mean atomic weight, $\bar{A}$, were determined for the ash (Table 2).

Once the change in energy and composition across the flame is known as a function of the initial density, the burning can be implemented without any recourse to temperature-dependent reaction rates or network solution. The flame is essentially a thin phase transition moving with a desired speed. The interface can be treated either using a level-set approach, or as was done here, as a thickened flame. The progress variable is the carbon mass fraction. As its mass fraction varies from 0.5 to zero, the binding energy per nucleon changes linearly from its value in the fuel, 7.828 MeV per nucleon, to its value given for carbon–free ash in the tables. The seven abundance variables are interpolated linearly from their values in the fuel (zero for all but C and O) and in the ash. The nuclear energy generation rate was

$$\begin{equation} \epsilon _{\rm nuc} = 9.64 \times 10^{17} (\Delta {\rm BE}/0.5)(d X(^{12}{\rm C})/dt) \ \ {\rm erg \ g^{-1} \ s^{-1}}. \end{equation} \tag{ 1 }$$

Here 0.5 is the initial mass fraction of carbon and ΔBE is the change in binding energy per nucleon following the burning of a composition of 50% C and 50% O to ash (Table 2). Though the oxygen usually burns as well, carbon was the progress variable. The change in the carbon abundance itself was calculated using a thickened flame model (Section 2.4). Several 1D studies were carried out with the CASTRO code that showed this prescription for energy generation and composition change gave the same resultant temperatures and densities for a conductive flame as in the one zone network simulation.

For fuel densities below 10^6.5 g cm⁻³ it was assumed that no nuclear burning took place and no burning was allowed on the grid until conduction had raised the local temperature in a carbon-rich composition above 2 × 10⁹ K.

2.3.2. Volumetric Burning in Nuclear Statistical Equilibrium

In NSE, all strong and electromagnetic nuclear reactions have come into a state of near balance with their reverse reactions, e.g., (n,γ), (p,γ), and (α, γ) reactions are balanced by (γ,n), (γ,p), and (γ, α) reactions. The distribution of abundances is then completely specified by just three parameters: the temperature, density and the neutron excess. The latter is usually measured by the electron mole number,

$$\begin{equation} Y_e= \Sigma Z_i Y_i, \end{equation} \tag{ 2 }$$

where Y_i is the mass fraction of the species "i" divided by its atomic weight and Z_i is its nuclear charge. For the initial composition and for the ash immediately behind the flame, there is neutron–proton equality so that Y_e = 0.50.

The NSE composition can evolve, however, both as a consequence of the changing temperature and density and the decrease of Y_e due to electron capture on nuclei at high densities ρ ≳ 10⁸ g cm⁻³. Of particular importance is the fact that, at the high temperatures appropriate to the central regions of exploding white dwarfs, NSE favors the partial photodisintegration of iron group nuclei into alpha-particles. This tends to reduce the initial increase in nuclear binding energy and decrease the average atomic mass of the ions. The temperature of the ash is not as high as one would obtain for fusion just to iron. However, as the matter expands and cools the alpha-particles reassemble, providing late time energy input to the explosion. Also important is the reduction of Y_e by electron capture. Since relativistic electrons are providing the pressure and since that pressure depends on the combination (ρY_e)^4/3, this can affect the dynamics. At constant pressure, a decrease in Y_e means an increase in the density and a reduction in buoyancy so plumes rise less rapidly. The value of Y_e also affects the nucleosynthesis and light curve of the supernova. For values close to 0.50, the principle nucleosynthetic product is radioactive ⁵⁶Ni, but for Y_e ≲ 0.48, ⁵⁴Fe and ⁵⁸Ni are more abundant.

The evolution of the NSE composition can only be followed using a large network of nuclei since many nuclei participate in capturing electrons at different densities and many reactions are involved in maintaining the equilibrium. We used the 127 isotope equilibrium network from the Kepler stellar evolution code (Weaver et al. 1978) with weak interaction rates taken from a variety of sources as discussed in Heger et al. (2001). Particularly important are the revisions to earlier rates by Langanke & Martinez-Pinedo (2000) and the inclusion of accurate rates for beta-decay and positron emission as well as electron capture. Temperature-dependent partition functions were used in calculating the abundances. The isotopes included in the NSE network were neutrons, ^1–3H, ^3–5He, ⁵Li, ⁹Be, ¹²C, ¹⁶O, ²⁰Ne, ²³Na, ^24–26Mg, ^26–28Al, ^28–30Si, ^30–33P, ^31–34S, ^35–37Cl, ^36–40Ar, ^39–43K, ^40–48Ca, ^44–49Sc, ^44–52Ti, ^47–54V, ^48–56Cr, ^51–58Mn, ^52–62Fe, ^54–64Co, ^56–66Ni, ⁵⁹Cu, and ⁶⁰Zn.

The table was constructed to include, as a function of the three independent variables, ρ, T, and Y_e, the binding energy per nucleon of the NSE distribution, its average atomic weight, $\bar{A}$, the helium mass fraction, and the time rate of change of Y_e resulting from all the weak interactions. Thirty one values of density from 10⁷ to 10¹⁰ g cm⁻³, 71 temperatures from 10⁹ to 10^10.4 K, and 21 values of Y_e from 0.50 to 0.40 were included in the table (46221 grid points). The temperature and density grid was in equal steps of the base 10 log. The values of ΔY_e were also equal so table look up time was minimized. Linear interpolation in log T, log ρ and Y_e thus gave the binding energy, mean atomic mass, and time rate of change of Y_e for any point of interest in the calculation. NSE was only assumed to exist in a zone if its density was over 10⁸ g cm⁻³, its temperature above 3 × 10⁹ K, its carbon mass fraction was below 1%, and the combined mass fractions of iron and helium were over 88%.

Due to the stiff dependence of abundances and binding energy on temperature and density, NSE calculations are subject to a well known instability. If the temperature goes up at constant density, for example, photodisintegration can give a large negative energy generation rate. That will tend to reduce the temperature for the next call which results in recombination and a large positive energy generation. Since the timescale for this adjustment is instantaneous, the code may take very small times steps and still have convergence problems. This difficulty was circumvented here by the addition of numerical inertia to the NSE calculation. The change in mean binding energy implied by a temperature–density change was not applied fully in a single time step but spread over three time steps. In practice, this kept the calculation stable and asymptotic values for Y_e and the binding energy were attained much more quickly than the timescale for the temperature or density to vary appreciably.

A few sample points from the NSE table are given in Table 4. In the actual table, six figure precision was maintained. The NSE network was not large enough, nor reliable weak interaction information available for the accurate calculation of abundances and $\dot{Y}_e$ for Y_e below 0.43. Fortunately, the smallest value of Y_e in any of our simulations here was above 0.45.

Table 4. Sample Characteristics of Nuclear Statistical Equilibrium

Log T	Log Density	Ye	X(He)	X(Si–Ca)	X(Fe group)	$\bar{A}$	BE/A	dYe/dt
9.80	9.30	0.500	1.601(−2)	5.174(−2)	0.9322	35.29	8.55635	2.617(−1)
9.80	9.30	0.495	1.205(−2)	4.097(−2)	0.9469	40.79	8.60340	1.920(−1)
9.80	9.30	0.490	9.798(−3)	2.832(−2)	0.9618	44.58	8.63952	1.313(−1)
9.80	9.30	0.485	8.243(−3)	1.645(−2)	0.9752	47.26	8.66989	8.148(−2)
9.80	9.30	0.480	7.000(−3)	7.951(−3)	0.9850	49.21	8.69495	4.668(−2)
9.80	9.30	0.475	5.975(−3)	3.575(−3)	0.9904	50.57	8.71502	2.645(−2)
9.80	9.30	0.470	5.042(−3)	1.601(−3)	0.9933	51.64	8.73201	1.449(−2)
9.80	9.30	0.465	4.072(−3)	7.147(−4)	0.9952	52.58	8.74541	7.197(−3)
9.80	9.30	0.460	3.159(−3)	3.524(−4)	0.9964	53.37	8.75220	3.470(−3)
9.80	9.30	0.455	2.439(−3)	1.989(−4)	0.9973	53.96	8.75367	1.745(−3)
9.80	9.30	0.450	1.841(−3)	1.115(−4)	0.9980	54.40	8.75201	8.264(−4)
9.80	9.30	0.445	1.339(−3)	5.506(−5)	0.9986	54.74	8.74627	3.423(−4)
9.80	9.30	0.440	1.001(−3)	2.451(−5)	0.9989	55.16	8.73595	1.219(−4)
9.80	9.30	0.435	8.112(−4)	1.029(−5)	0.9991	55.70	8.72337	1.182(−5)
9.80	9.30	0.430	7.294(−4)	3.667(−6)	0.9992	56.02	8.70984	−5.946(−5)
9.80	9.30	0.425	8.355(−4)	8.179(−7)	0.9991	55.10	8.69408	−1.119(−4)
9.80	9.30	0.420	1.933(−3)	6.037(−8)	0.9980	49.51	8.66726	−1.310(−4)

Download table as:  ASCIITypeset image

2.4. Flame Model

Burning during the deflagration stage of a SN Ia is quite complex and has challenged model builders for years. The most interesting burning during this stage happens at sufficiently low Karlovitz number, Ka ≲ 1, that the concept of a well defined laminar flame remains valid. At higher values of Ka, which happen for densities below about 10⁷ g cm⁻³, turbulence may lead to mixing, at the microscopic level, of cold fuel and hot ash and a possible detonation (Niemeyer & Woosley 1997; Khokhlov et al. 1997; Woosley et al. 2009). That possibility is ignored here. Thus, the picture of the flame is of a surface deformed by floatation, various instabilities, and possibly turbulence. Large scale motions caused by buoyancy are captured on the grid, but the smaller scale wrinkling is not. The manner in which this unresolved burning is treated is usually called the "subgrid model" for burning. Following Damköhler (1940), it is often assumed that if the large scale motions governing the spread of the flame are sufficiently resolved, e.g., the integral scale of the turbulence in the vicinity of the flame is well resolved, small scale burning will adapt so as to "digest" any entrained material at the necessary rate.

Adapting tools developed by the chemical combustion community, two approaches have been used by astrophysicists. One approach uses a "level set" (Osher & Sethian 1988; Reinecke et al. 1999) and a subgrid model for turbulence (Schmidt et al. 2006a, 2006b). The ash–fuel interface, found by solving a "G-equation," is moved at a speed given by sampling the velocity field on the grid and estimating the turbulent intensity. Typical speeds are ∼100 km s⁻¹, though a wide range of values is found (Röpke 2007). An alternative approach is the "thickened flame model" (Davies & Taylor 1950; Sharp 1984), where an artificial conductivity and energy generation are employed to spread the flame over several zones and force it to move at a prescribed speed. Both the flame speed and width can be much greater than the corresponding laminar values. The width is determined by the need to resolve the flame while the speed is sometimes taken to be the value to "burn out" structure at the grid level. It is well known (e.g., Timmes & Woosley 1992), that flame perturbations smaller than

$$\begin{equation} \lambda _{\rm min} \approx 4 \pi V^2 \left(g \frac{\delta \rho }{\rho }\right)^{-1} {} \end{equation} \tag{ 3 }$$

will be burned before they can grow. This is sometimes called the "fire-polishing length." Conversely, if one sets λ_min equal to the resolution, say a few zones, one can solve for a characteristic speed, V. Thus, $V \sim f \sqrt{g_{\rm eff} \Delta x}$, where g_eff = g(Δρ)/ρ, Δx is the grid resolution, and f < 1. This approach has been adopted by most groups studying the problem in the U.S. (Khokhlov 1995, 1996; Gamezo et al. 2003; Townsley et al. 2007). For a resolution of several km and g_eff ∼ 10⁹ cm s⁻², this prescription also gives a characteristic speed ∼100 km s⁻¹.

While we see great merit in the level set approach and are working on developing our own, the current study used a thickened flame. In particular, the nuclear burning timescale and diffusion timescale are defined here as,

$$\begin{equation} \tau _{{\rm nuc}} = \tau _{{\rm diff}} = \frac{n\Delta x}{V}, \end{equation} \tag{ 4 }$$

where n is the desired number of grid zones across the flame surface and Δx is the grid spacing, so that nΔx is the desired flame width. V is the desired flame speed. The nuclear energy generation rate and the conductivity are then derived from prescribed flame width and velocity.

In the flame, where the carbon abundance is larger than 1% and temperature is beyond 2 billion degrees, the carbon consumption rate is defined as,

$$\begin{equation} \frac{d X(^{12}{\rm C})}{dt} = {-} \frac{0.5}{\tau _{{\rm nuc}}}, \end{equation} \tag{ 5 }$$

where 0.5 is the initial carbon mass fraction. Additionally,

$$\begin{equation} \frac{d({\rm BE})}{d X(^{12}{\rm C})} = {-} \frac{\Delta {\rm BE}}{0.5}, \end{equation} \tag{ 6 }$$

so that nuclear burning rate in Equation (1) can be derived.

For partially burned cells inside the flame, the current BE and $\bar{A}$ are interpolated using the equations:

$$\begin{equation} \frac{{\rm BE} - {\rm BE}_{f}}{X(^{12}{\rm C})-0.5} = -\frac{{\rm BE}_f-{\rm BE}_i}{0.5}, \end{equation} \tag{ 7 }$$

$$\begin{equation} \frac{\bar{A} - \bar{A}_{f}}{X(^{12}{\rm C})-0.5} = -\frac{\bar{A}_f-\bar{A}_i}{0.5}. \end{equation} \tag{ 8 }$$

The opacity is also defined using τ_diff from the thickened flame model as

$$\begin{equation} \kappa = \frac{c \tau _{\rm diff}}{3 \rho \Delta ^2}, \end{equation} \tag{ 9 }$$

and conductive coefficient is,

$$\begin{equation} D = 4acT^3/3 \rho \kappa , \end{equation} \tag{ 10 }$$

where c is light speed, Δ is the flame thickness, and a = 7.5646 × 10⁻¹⁵ erg cm⁻³ K⁻⁴ is radiation density constant. The conductivity is only taken to be non-zero in zones whose temperature is above 1 billion degree. This criterion allows us to distinguish the flame and post-flame region from material heated diffusively in front of the flame. Suppressing spurious burning in the diffusion front ahead of burning front is necessary to prevent an artificial "diamond" shape from appearing in the propagation of the flame.

From a survey of 1D flames using the above prescription, we determined the multiplier, f_w and f_v (both of the order of unity), of the flame width and speed to get prescribed values at different densities (see Table 5). In the calculation, the left side is set to be hot ashes, and right side is fresh fuel (Figure 1). The boundary condition is impermeable at the right side, free flow at the left side. Flame speeds are measured with respect to the fuel (Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. Flame Adjustment Parameters, f_w,f_v, at Different Densities

log(ρ)	fw	fv
6.90	2.40	2.90
7.00	2.00	2.20
7.15	1.80	1.95
7.30	1.40	1.55
7.50	1.00	1.20
7.70	1.00	1.05
7.85	1.20	1.10
8.00	1.20	1.15
8.30	1.60	1.30
8.50	1.80	1.40
8.70	2.00	1.60
8.85	2.20	1.70
9.00	2.40	1.75
9.20	2.80	2.00
9.30	3.00	2.20
9.40	3.40	2.30
9.45	3.40	2.30

Download table as:  ASCIITypeset image

2.5. Choice of Flame Speed Parameter

It thus remains only to pick the desired flame speed. As discussed in Section 2.4, previous studies using a subgrid model to describe the turbulent flow (Reinecke et al. 1999, 2002; Röpke & Hillebrandt 2005) give a speed on the order of 100 km s⁻¹. In particular, Röpke (2007) studied the probability density function (PDF) of turbulent velocity at a density interval of 1–3 × 10⁷ g cm⁻³ in the deflagration phase. His plot of the PDF shows a turbulent velocity peaking approximately 100–300 km s⁻¹. Thickened flame models for typical values of resolution give a comparable, though somewhat smaller, value for V. The laminar speed, which sets a lower bound to the overall flame speed, is almost this large at the high density where burning begins. At a density of 2 × 10⁹ g cm⁻³, the laminar speed is 76 km s⁻¹ (Timmes & Woosley 1992), though it declines rapidly for densities much less than 10⁹ g cm⁻³. By the time the density has declined, though, floatation is dominating with typical speeds of thousands of km s⁻¹ and the turbulent estimate or the fire polishing estimate is valid.

This all suggests that a constant flame speed of 100 km s⁻¹ is a good place to begin and this is the value used in our standard calculation. However, we also carried out other otherwise identical calculations using fixed V = 50 and 200 km s⁻¹. The lower value here is actually a little less than the laminar speed for a composition of 50% carbon, 50% oxygen at a density of 2.4 × 10⁹ g cm⁻³ where the laminar speed has a value of 89 km s⁻¹ (Timmes & Woosley 1992). However, this high density and laminar speed are only maintained a short time and, in the interest of maintaining a constant value of V throughout the run, the lower value of 50 km s⁻¹ was used. Interestingly, had the carbon mass fraction been 20%, 50 km s⁻¹ would have been the correct value even at this high density.

To test the sensitivity of the results to the local turbulence, we also performed a simulation where the flame speed depended quite simply on the local turbulent intensity. In particular, for a zone within the flame (as described by the nuclear burning regime discussed in Section 2.3.1) the turbulent velocity fluctuations, v_turb, are calculated in the usual way:

$$\begin{equation} v_{\rm turb} = (\overline{v^2} - \overline{v}^2)^{1/2}, \end{equation} \tag{ 11 }$$

where the overbar denotes spatial averaging over a cube of 3³ zones centered about the zone in question. The value of v_turb from Equation (11) is then rescaled to the grid scale assuming Kolmogorov scaling. The local flame speed, V, is then set to the maximum between the laminar speed and the local turbulent intensity. As we shall see, the answer is gratifyingly insensitive to the choice of V—constant or set to the local turbulent intensity—within the limited range explored.

3.1. Calculations

3.2. Parsing Iron

Once the composition had ceased to evolve, the electron mole number, Y_e, and helium mass fraction were used to distinguish the dominant forms of iron (^{54, 56}Fe, ^{56, 58}Ni) produced in the explosion. Three cases are considered based on neutron excess η = 1–2 Ye (Hartmann et al. 1985),

$$\begin{equation} \eta =\sum \frac{(N_{i}-Z_{i}) \ X_{i}}{A_{i}}. \end{equation} \tag{ 12 }$$

For the conditions studied here, most of the neutrons will be bound in ⁵⁴Fe, ⁵⁶Fe, and ⁵⁸Ni with the relative proportions depending upon the values of η and the helium mass fraction, X(He). For near proton–neutron equality, i.e., η less than 1/27, and X(He) less than 0.20, it is assumed that most of the neutrons reside in ⁵⁴Fe with a small fraction in ⁵⁸Ni. Specifically η = 1/27 X(⁵⁴Fe) + 1/29 X(⁵⁸Ni) and X(⁵⁸Ni) = f X(⁵⁴Fe) with f = 0.35 (Hartmann et al. 1985). If X(He) is larger than 20%, however, the so-called α-rich freeze out is assumed to produce mostly ⁵⁸Ni so that f = 10. For η less than 1/14, but greater than 1/27, X(⁵⁶Ni) = 0 and the sum of the ⁵⁴Fe, ⁵⁶Fe, and ⁵⁸Ni mass fractions is 1−X(He). Neutron conservation additionally gives η = 1/27 X(⁵⁴Fe)+1/29 X(⁵⁸Ni)+1/14 X(⁵⁶Fe), and the condition X(⁵⁸Ni) = f X(⁵⁴Fe), as before, can be solved to give the three isotopic abundances separately. For η greater than 1/14, which was only true for a small fraction of the matter where burning happened at very high density, ⁵⁶Fe was assumed to be the only isotope for iron group elements.

Figure 3 shows the propagation of the flame front at different star times for our standard run, model A100. The flame is wrinkled and distorted by the RT instability, with buoyant, hot ash bubbles floating toward the surface and spreading while cold, denser matter moves outward more slowly and is overtaken. Viewed in Lagrangian coordinates, the unburned matter appears to migrate inward. At very early times, one sees finger-like structures resulting from the initial perturbations, but these quickly mushroom and become nonlinear. New plumes form and the structure becomes complex.

Zoom In Zoom Out Reset image size

During the early growth of the reactive RT instability, burning effectively smooths out perturbations below the fire polishing length (Equation (3)). As the RT modes begin to merge, interact, and break down, they become wrinkled on ever smaller scales (e.g., Zingale et al. 2005). Indeed, the flame experiences the effects of turbulence on scales down to the Gibson scale—the scale at which the laminar speed equals the turbulent speed. Furthermore, the turbulent scaling study of Ciaraldi-Schoolmann et al. (2009) shows that the nature of the turbulence shifts from being dominated by RT scaling (Equation (3)) to Kolmogorov scaling for isotropic eddies on scales smaller than ∼15 km, or about the initial λ_min. It should be noted that with our choice of grid resolution in the present studies, we are not resolving the turbulent flame on scales smaller than λ_min.

Several effects lead to the evolution of λ_min. As the flame moves out, the greater interior mass enclosed tends to increase the gravitational acceleration, but at the same time the star expands, causing it to decrease. The density contrast, Δρ/ρ also gets larger at lower density. The combined effect is to decrease g_eff slowly, though its value remains roughly constant at about 10⁹ cm s⁻² throughout much of the burning. Since the flame speed is also a constant, 100 km s⁻¹ in this case, λ_min does not vary greatly from about 10 km. While there are always at least several zones within λ_min, this does not necessarily mean that structures of this scale are well resolved. Since the thickened flame itself is typically a few to five zones wide and any structure must be bounded on all sides where there is fuel by flames, λ_min is only marginally resolved. The smallest structures in Figure 6 are roughly 40 km across. The length scale λ_min is more fully resolved for the model A200 run, but less resolved for model A50. Ideally, greater resolution should be employed, but it would be prohibitively expensive to run with even twice the current resolution since run timescales as the fourth power of the resolution, three for the three dimensions plus one for the shorter Courant time step. Fortunately, the most interesting burning happens on the finest grid scale, which is 1 to 2 km during the nucleosynthesis epoch. By 1.2 s when the finest grid becomes 4 km, g ∼ 8 × 10⁸ cm s⁻² and Δρ/ρ ∼ 0.3, so λ_min has already increased to 50 km. Visually, one also sees more structure in the model A50 run (Figure 6), suggesting that the structure is sensitive to V and we are actually resolving λ_min.

The flame first breaks out of the surface in the model A100 run 1.21 s after ignition, by which time the star has already expanded by a factor of about 2.5 and the central density has declined to 5 × 10⁷ g cm⁻³. This decrease in density significantly affects the products of the burning and allows the production of significant IMEs. The largest speed up to this point is ∼10,000 km s⁻¹ in the outermost layers of the star, but the flow is still subsonic in most of the star. Nuclear burning is over at this point and, as the internal energy converts into kinetic energy, the maximum speed rises to 25,000 km s⁻¹ at 2.76 s.

Figure 4 shows how kinetic energy, binding energy and internal energy evolve with time. At the beginning of the explosion, the net binding energy, i.e., the gravitational binding energy minus the internal energy, is about 5 × 10⁵⁰ erg. The final kinetic energy at infinity is about 6 × 10⁵⁰ erg, so the total energy released thermonuclear burning energy is 1.1 × 10⁵¹ erg. The supernova ejecta is mainly composed of 0.63 M_☉ iron group elements and 0.16 M_☉ IME. The mass of ⁵⁶Ni and the final kinetic energy suggests a relatively weak event in SN Ia.

Zoom In Zoom Out Reset image size

A sequence of plots in Figure 5 shows the evolution of the angle-averaged supernova composition. At early time, iron group elements are the major products, but after the flame spreads to the low density region (less than 5 × 10⁷ g cm⁻³) at 1 s, more IME are made. As the iron-group ashes continue to move out at higher speeds than the unburned fuel, boosted along by the recombination of helium, there is a gradual accumulation of carbon and oxygen near the center of the explosion. At the end, about 20% unburned carbon and oxygen are left in the central region.

Zoom In Zoom Out Reset image size

The calculation was repeated using two different flame speed parameters, models A50 and A200, and followed until 1.8 s when the nuclear burning was over (but not into the final coasting phase). Comparisons of the flame front morphology using the different flame speeds are shown in Figure 6. Slower flame speeds result in increased structure as noted, but, as Figure 7 shows, the total amount of burned material is roughly the same for flame speeds of 100 km s⁻¹ and 200 km s⁻¹. The slower burning speed is compensated for by a larger surface area. For a flame speed of 50 km s⁻¹, the calculation produces 0.73 M_☉ iron, 0.17 M_☉ IME, and released 1.3 × 10⁵¹ erg. This is significantly more burning than either of the other two series A runs but still within 10% of their average. We do not know if this limited variation is due to the small number of runs or reflects real physical dependence on the parameter. There is a tendency, however, for more burning to occur in an explosion that starts slower because it remains denser longer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The B series simulations used a different random number generator seed than the A series simulations; the random number generator applied spherical perturbations of random size and amplitude to the initial flame surface. Figure 8 compares the X(¹²C) = 0.49 isosurface of the constant flame speed model B100 (left) to that of the turbulent flame speed model BT (right) at 0.3 s (top), 0.5 s (middle), and 0.7 s (bottom). Morphologically, models B100 and BT are nearly identical. A comparison of the left column of Figure 8 (B100) with the center column of Figure 6 (A100) at the same time elucidates the difference in initial model conditions between the two series. The A series perturbations are all nearly of the same magnitude, whereas the B series initially has a few slightly dominant perturbations that grow more rapidly and develop into large-scale RT plumes; this feature is similar to published 2D studies of central ignition with perturbed spheres (see Jackson et al. 2010; Krueger et al. 2012, for example). The extraordinarily large appearance of the plume in the lower edge of the plots is a result of that plume's coming toward the viewer. Neglecting the dominant plumes, at any given time the size of the ash bubble is about the same in all runs.

Zoom In Zoom Out Reset image size

Figure 9 shows the total mass of iron group elements (red) and IME (blue) produced for the B series of models. The dotted lines are for model B100, and the dashed lines are model BT. The two thin, solid lines show the percent difference between the two models where positive (negative) values indicate model B100 producing more (less) material than model BT. These two different flame models agree remarkably well. Initially, the constant flame speed for model B100 is greater than the laminar flame speed, which is the speed at which model BT burns in the absence of turbulence. This difference causes model B100 to burn more material more rapidly than model BT. By about 0.2 s, some zones within the flame experience turbulent velocities that are greater than the local laminar flame speed and/or the constant 100 km s⁻¹ flame speed of model B100; the flame in model BT begins to increase in speed, which burns more material that then decreases the difference between the two flame prescriptions. As the flame grew in size, the number of fine-resolution zones increased. After ∼0.31 s of evolution, there were over one billion zones in the finest level and we de-refined the simulation for computational efficiency. This increase by a factor of two in zone size (in each direction) caused a corresponding increase in the turbulent intensity on the grid scale by a factor of 2^1/3 according to the Kolmogorov scaling we assumed in our turbulent flame speed model. This is evident in Figure 9 as a small but sudden increase in the percent difference between elemental yields with the constant flame speed model and the turbulent flame speed model. Another de-refinement occurred around 0.61 s, but the change in percent difference between the models is lost in the noise.

Zoom In Zoom Out Reset image size

Around 0.4 s, the turbulent flame in model BT has "caught up" to the constant flame speed model B100 in terms of total mass produced in each elemental group. The PDFs for the turbulent speed at various times in Figure 10 show that indeed the most likely turbulent flame speed used around this time in model BT is ∼100 km s⁻¹ as is used in model B100. The distribution does have an extended tail, which shows some subset of the zones within the flame experience turbulent intensities greater than the constant flame speed model B100. This causes the turbulent flame speed model BT to overproduce iron and IME relative to model B100 (negative percent differences in Figure 9). It is around this point in the evolution that simple interpretation of the differences in iron and IME yields between the two flame speed prescriptions becomes difficult. There is a constant battle between (1) the small differences in flame structure being amplified by the (dominant) acceleration due to buoyancy, (2) the dominant plumes reaching a region where the density gradient becomes steep, and (3) the shear induced turbulence that tends to redistribute the turbulent flame speed to ever higher velocities. Even though the total mass of iron and IME is increasing in both models, the differences between model B100 and BT fluctuate between positive and negative due to the interplay of the above mentioned forces.

Zoom In Zoom Out Reset image size

Throughout the B series runs, the choice of flame speed model (constant versus turbulent) made less than a 7% difference in the iron and IME group total yields. Comparing, however, slightly different initial conditions but the same (constant) flame speed model—models A100 (Figure 7) and B100 (Figure 10)—at, say, 0.6 s shows nearly a 50% difference in total iron group yield. This seems to suggest that the simulation results are much more sensitive to the initial conditions of the hot ash than they are to the flame speed prescription, at least in the case of a constant flame speed model and a model where the flame speed depends simply on the turbulent intensity on the grid within the flame.

The elemental yields of our models and the nuclear energy released by each are given in Table 6. It is important to note that the B series models were not run to completion so their tabulated elemental yields should not be compared with the other models in the table, but only with each other to assess the role of a simplified turbulent flame model. Results from the A series studies with V = 100 and 200 km s⁻¹ resemble each other closely and are also similar to those of Röpke et al. (2006). They do show, however, significantly less production of IMEs than Röpke et al. (2007a). We do not know the specific cause of the large production of IMEs in the Röpke study, but note that they are significantly larger than in any other study. Perhaps this reflects the initial conditions and igniting a very large number of ignition points or an unusually high flame speed from the turbulent subgrid model at low densities (see below). García-Senz & Bravo (2005) found around 0.45 M_☉ ⁵⁶Ni was produced in the explosion with ∼15 initial igniting bubbles near the center, which is similar to our series A results, but IME production in their simulation was much smaller. Gamezo et al. (2005) found only 0.47 M_☉ iron group elements were produced during the deflagration assuming an initial ignition region with a radius of only 30 km.

The actual mass of ⁵⁶Ni implied by observations of SN Ia is highly variable ranging from 0.1 to 0.9 M_☉ (Mazzali et al. 2007), but a typical value is about 0.6 M_☉ (Branch & Khokhlov 1995). A lower limit to the total burned mass derived from observations of the velocity of silicon, applying the density–velocity structure in W7 model, is around 1.05 M_☉ (Mazzali et al. 2007). Our series A calculations found that approximately 0.9–1.0 M_☉ of the 1.38 M_☉ white dwarf was burned to heavier elements regardless of the assumed flame speed, not too much smaller than the optimal value. This burning certainly released enough energy to unbind the star, but is well short of what is needed to account for the brighter SN Ia.

Perhaps this is not too surprising. If typical SN Ia ignite off center and later detonate, then one expects a correlation between the displacement of the ignition and the brightness of the supernova. Explosions igniting off-center burn less before detonating and are thus denser when the detonation occurs. They thus make more ⁵⁶Ni and are brighter. Central ignition is a limiting case where the burning in the deflagration stage is maximal. A subsequent detonation makes less ⁵⁶Ni and the total synthesis may be less, hence a fainter supernova. The puzzle then is not why deflagrations are faint, but the origin of some SN Ia that are fainter still.

The main difficulty faced by pure deflagrations, though, is not their brightness, but their spectrum. They have too much low velocity unburned fuel near their middles (Marion et al. 2003). More than 70% of the central ejecta in our standard model A100 was carbon and oxygen. Even though the central regions initially burned to completion, unburned fuel is later found near the center of the star because the ashes have floated above the fuel, leaving it behind as they move outward with greater speed than their cold surroundings (see Figure 11). Iron group elements and IME in the ejecta continue mixing radially even after burning has gone out; the layered structure of different chemical products present at the end of burning but before much expansion is hardly seen. This is contrary to the observed spectral lines of carbon and silicon in the supernova ejecta, which indicates velocities mostly beyond 10,000 km s⁻¹ (Filippenko 1997; Branch et al. 2003; Benetti et al. 2005; Parrent et al. 2011). Previous numerical studies have demonstrated the same shortcomings out of the centrally ignited pure deflagration model. For example, Gamezo et al. (2003) assumed a centrally burned region of radius 30 km, and carried out the following deflagration simulation with grid-length and effective gravity-dependent flame speeds in the thick-flame model. Eighty to ninety percent of the low velocity materials near the center are composed of carbon and oxygen at 1.9 s in their study. Other numerical works initialized a number of hot bubbles around the center instead of the whole central region, including 2D and 3D studies using a different flame tracking method and subgrid model (Reinecke et al. 2002; Röpke et al. 2006), and 3D calculations using a smoothed particle hydrodynamics code (García-Senz & Bravo 2005). Central carbon concentrations were found in all these studies. However, parameter studies have shown that larger numbers (more than 150) of initially Gaussian distributed hot bubbles with a maximum ignition radius of 180 km leave much less unburned fuel in central regions, resulting a spectrum better consistent with the observations (Röpke et al. 2006, 2007a). This specific ignition pattern is unlikely to happen at least in mildly turbulent convective core of the white dwarf, though (Kuhlen et al. 2006; Zingale et al. 2009, 2011).

Zoom In Zoom Out Reset image size

This weakness of the pure deflagration model might be alleviated by a later transition from deflagration to detonation (DDT) in the center or at a certain radius (Gamezo et al. 2005). The supersonic detonation wave is expected to burn off the inhomogeneous unburned materials, release enough energy to match classic SN Ia, and give stratified shells of product elements. The synthetic light curves and spectra of the DDT explosion agreed well with the observations of a normal SN Ia (Kasen et al. 2009). It has been suspected that a DDT may take place when the flame enters the distributed burning regime at the density of 1 ∼ 3 × 10⁷ g cm⁻³, where turbulence starts to deform the internal flame structure and mix the hot ashes and fuel (Niemeyer & Woosley 1997). The existence of strong turbulence, with velocity fluctuations larger than one-fifth of local sound speed, is crucial for a spontaneous detonation (Woosley et al. 2009), and Röpke (2007) suggested the possibility of these high turbulent velocities at the corresponding density range from their deflagration simulations. We plan to study the outcome of a later detonation following the current deflagration calculations with two free parameters, the location and time of the detonation wave initiation.

We thank Andy Aspden and Mike Zingale for helpful discussions regarding the turbulent flame physics and the set up of problems in hydrostatic equilibrium using the CASTRO code. The work at LBNL was supported by the SciDAC Program of the DOE Office of High Energy Physics and by the Applied Mathematics Program of the DOE Office of Advance Scientific Computing Research under U.S. Department of Energy under contract No. DE-AC02-05CH11231. The work at UCSC was supported by the DOE SciDAC program, under grant No. DE-FC02-06ER41438 and by the NASA Theory Program (NNX09AK36G). Computer time for the calculations in this paper was provided through a DOE INCITE award at the Oak Ridge Leadership Computational Facility (OLCF) at Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC05-00OR22725. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC02-05CH11231. Early developmental work necessary to this study was carried out on the Pleiades supercluster at UCSC. Pleiades was purchased using a grant from the MRI Program of the NSF (AST-0521566). Visualizations were performed using the VisIt package. We thank Gunther Weber and Hank Childs for their assistance with VisIt.