EVALUATING SYSTEMATIC DEPENDENCIES OF TYPE Ia SUPERNOVAE: THE INFLUENCE OF DEFLAGRATION TO DETONATION DENSITY

Type Ia supernovae (SNe Ia) are bright stellar explosions that are characterized by strong P Cygni features in Si and by a lack of hydrogen in their spectra (see Hillebrandt & Niemeyer 2000; Filippenko 1997 and references therein). Observations of SNe Ia (serving as distance indicators; Phillips 1993; Riess et al. 1996; Albrecht et al. 2006) are at present the most powerful and best proven technique for studying dark energy (Riess et al. 1998; Perlmutter et al. 1999; Kolb et al. 2006; Hicken et al. 2009; Lampeitl et al. 2010), and, accordingly, there are many observational campaigns underway striving to gather information about the systematics of these events and to better measure the expansion history of the universe (see Kirshner 2009 and references therein).

The most widely accepted model for these events is the "single-degenerate" scenario, which is the thermonuclear explosion of a white dwarf (WD) composed principally of ¹²C and ¹⁶O that has accreted mass from a companion (for a review of explosion models, see Hillebrandt & Niemeyer 2000; Livio 2000; Röpke 2006). The peak brightness of the supernova (SN) is set by the radioactive decay of ⁵⁶Ni produced in the explosion to ⁵⁶Co and then to ⁵⁶Fe. The empirical "brighter is broader" relation between the peak brightness of the light curve and the decay time from its maximum is understood to follow from the fact that both the luminosity and opacity are set by the mass of ⁵⁶Ni synthesized in the explosion (Arnett 1982; Pinto & Eastman 2000; Kasen & Woosley 2007). Because of the dependence of the light curve on the amount of radioactive ⁵⁶Ni synthesized during an explosion and the ability to infer the ⁵⁶Ni yield from observations (Mazzali et al. 2007), research into modeling thermonuclear SNe typically focuses on the production and distribution of ⁵⁶Ni as well as other nuclides (such as ²⁸Si) as the measure with which to compare models to observations.

One-dimensional simulations of the single-degenerate case showed that the most successful scenario is an initial deflagration (subsonic reaction front), born in the core of the WD which at some point becomes a (supersonic) detonation, i.e., a deflagration–detonation transition (DDT; Khokhlov 1991; Höflich & Khokhlov 1996). A delayed detonation naturally accounts for the high-velocity Ca features (Kasen & Plewa 2005) and the chemical stratification of the ejecta. While these one-dimensional models are able to reproduce observed features of the light curve and spectra, much of the physics is missing from these models. The presence of fluid instabilities during the deflagration warranted the development of multidimensional models, allowing a physically motivated calculation of the velocity of the burning front and thus removing a free parameter. By relaxing the symmetry constraints on the model, buoyancy instabilities are naturally captured leading to a strong dependence on the initial conditions of the deflagration that often result in too little expansion of the star by the time DDT conditions used in previous one-dimensional studies are met (Niemeyer et al. 1996; Calder et al. 2003, 2004; Livne et al. 2005). Multidimensional models may reach the expected amount of expansion prior to the DDT with the choice of particular ignition conditions and thus retain the desirable features from one-dimensional models (Gamezo et al. 2004; Plewa et al. 2004; Röpke et al. 2006; Jordan et al. 2008).

Due to the strong dependence on ignition conditions, multidimensional simulations of the DDT model are able to produce a wide range of peak luminosities (via the production of a range of ⁵⁶Ni yields) consistent with a common explosion mode suggested from observations (Mazzali et al. 2007). Differences in the mass of synthesized ⁵⁶Ni can follow from properties such as metallicity and central density of the progenitor and/or differences in the details of the explosion mechanism such as the density at which the transition from deflagration to detonation occurs. Timmes et al. (2003) found that metallicity should affect the ⁵⁶Ni yield based on approximate lepton number conservation. The metallicity sets the fractional amount of material in nuclear statistical equilibrium (NSE) that is radioactive ⁵⁶Ni. Bravo et al. (2010) found a stronger dependence on metallicity due to a significant amount of ⁵⁶Ni that is synthesized during incomplete Si burning.

Observational results to date are consistent with a shallow dependence of ⁵⁶Ni mass on metallicity but are unable to rule out a trend entirely (Gallagher et al. 2005, 2008; Neill et al. 2009; Howell et al. 2009). Determining the metallicity dependence is challenging not only because the effect appears to be small, but also due to the difficulty in measuring accurate metallicities for the parent stellar population and problems with strong systematic effects associated with the mass–metallicity relationship within galaxies (Gallagher et al. 2008). This effect is also difficult to decouple from the apparently stronger effect of the age of the parent stellar population on the mean brightness of SNe Ia (Gallagher et al. 2008; Howell et al. 2009; Krueger et al. 2010). Howell et al. (2009) note that the scatter in brightness of this observed relation is unlikely to be explained by the effect of metallicity. In general, the source of scatter can be explained by the development of fluid instabilities during the deflagration phase that contribute to differing rates of expansion between instances of SNe. However, by considering the effect of metallicity on the DDT density, scatter in the metallicity relation may be enhanced beyond its intrinsic value inferred from fluid instabilities.

Townsley et al. (2009) investigated the direct effect of metallicity via initial neutron excess and found it to have a negligible influence on the amount of material synthesized to NSE. However, the neutron excess sets the amount of material in NSE that favors stable Fe-group elements over radioactive ⁵⁶Ni. Therefore, the initial metallicity directly influences the yield of ⁵⁶Ni. In this work, we expand that study to include the potential indirect effect of metallicity in the form of the ²²Ne mass fraction (X₂₂) through its influence on the density at which the DDT takes place. To this end, we explore the effect of varying the transition density, a proxy for varying the microphysics that determine the conditions for a DDT. The conditions under which a DDT occurs are still a subject of debate. Niemeyer & Woosley (1997), Niemeyer & Kerstein (1997), and Khokhlov et al. (1997) proposed that a necessary condition is the transition to a distributed burning regime, in which turbulence disrupts the reaction zone of the flame. More recent numerical studies (Pan et al. 2008; Woosley et al. 2009; Schmidt et al. 2010) that describe the conditions for a DDT include dependencies on the turbulent cascade and the growth of a critical mass of fuel with sufficiently strong turbulence. For this study, we assume the DDT density to be the density at which thermonuclear burning is expected to enter the distributed regime. This choice links the explosion outcome to the dynamical evolution of the progenitor density structure during the deflagration phase. By analyzing the effect of transition density on the NSE yield, we can later analyze how the details of the microphysics affect the DDT density. For the purposes of this study, we will consider only the effect of X₂₂ on the DDT density. In reality, the ¹²C mass fraction will also be important in determining the DDT density, but we choose to leave the exploration of the effect of ¹²C to future work. Many other possible systematic effects exist that are outlined in Townsley et al. (2009), such as the central ignition density (Krueger et al. 2010), which are all held fixed in this study. The interdependence of all of these effects will be considered in the construction of the full theoretical picture in a future study.

We describe the details of our model in Sections 2 and 3, and the properties of our statistical sample in Section 4. We present our findings on the dependence of transition density on NSE yield in Section 5.1. In Sections 5.2 and 5.3, we assume a dependence of the transition density on ²²Ne content and construct the functional dependence of the ⁵⁶Ni yield on metallicity through the ²²Ne content. In Section 6, we discuss our conclusions and future work.

In order to include relevant processes in explosion models, we first estimate the compositional profile of the progenitor WD just prior to the birth of the flame. We begin by estimating the compositional profile resulting from the evolution of the post-main-sequence star that later becomes a WD. Recall that the initial metallicity of the star is, by mass, almost entirely in the form of CNO. As a result of the C–N–O cycle, these all end up in helium layers as ¹⁴N, the target of the slowest proton capture in the cycle. Subsequent reactions during helium burning convert ¹⁴N into ²²Ne; therefore, X₂₂ is proportional to the initial metallicity of the star (Timmes et al. 2003). The composition of the inner portion of the star (≈0.3–0.4 M_☉) is set during core helium burning resulting in a reduced carbon mass fraction with respect to that of the outer layers, which is set by shell burning on the asymptotic giant branch (see Straniero et al. 2003 and references therein).

At some point after the WD has formed, it begins to accrete material from its companion. As the mass of the accreting WD approaches the Chandrasekhar mass limit, the core temperature and density increase such that carbon begins to fuse. The energy released by carbon burning drives convection in the core. The convective carbon-burning ("simmering") phase lasts approximately 1000 years before the central temperature is high enough to spark a thermonuclear flame (Woosley et al. 2004). During the simmering phase, carbon is consumed from the convective core. Concurrently, though, the convective zone grows with increasing central temperature, pulling in relatively carbon-rich material from the outer layer. Thus, the net effect is to increase the carbon abundance in the convective region. We show the growth of the convective zone in Figure 1, where dashed lines show the compositional profile prior to the simmering phase and the corresponding solid lines show the compositional profile at the start of our simulations of the explosion.

Zoom In Zoom Out Reset image size

For our WD models, we consider a parameterized three-species compositional structure consisting of ¹²C, ¹⁶O, and ²²Ne, which is sufficient to capture the nuclear burning rates. Since ²²Ne is the only element in our parameterization that has more neutrons than protons, the neutron excess from simmering is accounted for by a parameterized ²²Ne mass fraction

$$\begin{equation} X_{22}^\prime = X_{22}+ \Delta X_{22}(\Delta Y_{e}) {\rm,} \end{equation} \tag{ 1 }$$

where X₂₂ is proportional to the initial metallicity coming from helium burning (Timmes et al. 2003) and ΔX₂₂(ΔY_e) represents the effective enhancement of ²²Ne in the core following from the change in electron fraction during convective carbon burning. The electron fraction (Y_e) is related to X₂₂ by

$$\begin{equation} Y_{e}= \frac{10}{22} X_{22}+ \frac{1}{2}(1-X_{22}) \end{equation} \tag{ 2 }$$

such that a change in X₂₂ constitutes a change in Y_e. The prime on X'₂₂ in Equation (1) indicates inclusion of effects of neutronization during carbon burning and this convention applies to the expressions below. We choose the composition of the WD at the end of the simmering phase to consist of X₁₂ = 0.4 and X'₂₂ = 0.03 in the core and X₁₂ = 0.5 and X'₂₂ = 0.02 in the outer layer. Note that X₂₂ = X'₂₂ in the outer layer because neutronization due to carbon burning only occurs in the convective core.

For comparison, a compositional profile can be estimated for the WD prior to the onset of carbon burning. For simplicity, consider the production of one neutron for every 6 ¹²C burned (i.e., dY_e/dY₁₂ ≈ 1/6, where Y₁₂ is the molar abundance of ¹²C) (Chamulak et al. 2008; Piro & Bildsten 2008). For ΔX₂₂ = 0.01 in our progenitor model, this constitutes burning 0.04 M_☉ of carbon during simmering. Assuming that prior to simmering, the core mass is ≈0.3 M_☉ (Straniero et al. 2003) and the outer layer consists of X₁₂ = 0.5 material, we can account for all ¹²C and conserve the total mass of the WD to estimate X₁₂ ≈ 0.2 in the carbon-reduced core of the WD prior to simmering as shown in Figure 1.

Table 1 shows the composition of the progenitor WD before and after the simmering phase using the parameterized ²²Ne mass fraction as well as the core mass. Note that throughout this study, we choose to neglect any variation of the amount of neutronization during simmering due to the initial Z. Accordingly, ΔX₂₂(ΔY_e) is treated as a constant and, therefore, dX'₂₂ = dX₂₂ and derivatives involving the true ²²Ne mass fraction proportional to metallicity, X₂₂, are equivalent to derivatives involving the parameterized ²²Ne mass fraction, X'₂₂.

Because we consider the enhancement of the laminar flame speed by ²²Ne (Chamulak et al. 2007), we need to consider the effects of our parameterization of the ²²Ne content with care. In actuality, the neutronization during carbon simmering produces ¹³C, ²³Na, and ²⁰Ne and not ²²Ne (Chamulak et al. 2008). A priori, this difference could alter the nuclear burning rates and hence the laminar flame speed and width. Table 2 shows the laminar flame speeds and widths for the same Y_e exchanging ΔX₂₂ for carbon-simmering ash for a composition of 40% ¹²C, 2.0% ²²Ne, 55.5% ¹⁶O, 0.6% ¹³C, 0.9% ²⁰Ne, and 1.0% ²³Na by mass using the same method as described by Chamulak et al. (2007). The flame speeds and widths calculated using the parameterized ²²Ne mass fraction (3% by mass) are denoted with primes. For high densities (≳2.5 × 10⁸ g cm⁻³) in which the flame speed is not dominated by the buoyancy-driven Rayleigh–Taylor (RT) instability, the difference is ≲5%; therefore, our parameterization of the neutronization via X₂₂ accurately captures the corresponding enhancement of the laminar flame speed in this regime.

While the central temperature and central density of the progenitor just prior to the birth of the flame are primarily set by the accretion history of the WD, which varies and is largely unknown, we choose a fiducial central density of ρ_c = 2.2 × 10⁹ g cm⁻³ and central temperature of T_c = 7 × 10⁸ K. We construct isentropic profiles of density and temperature in the (convective) core and isothermal profiles in the (thermally conductive) outer layer while maintaining hydrostatic equilibrium. We choose a fiducial isothermal temperature of the outer layer to be T_iso = 10⁸ K. The total mass of the WD progenitor model is 1.37 M_☉.

Due to the difference in composition between the core and outer layer, requiring a neutral buoyant condition at this interface no longer reduces to a continuous temperature as it does for a compositionally uniform WD (Piro & Chang 2008). In reality, we might expect some convective overshoot and mixing between the core and the outer layer, but because the physics of convective overshoot are complex and not well understood, the width of the transition region is unknown. For simplicity, we assume no mixing region. The composition, density, and thermal profiles of the progenitor used for this study immediately prior to deflagration are the solid lines shown in Figure 1.

We use a customized version of the FLASH Eulerian compressible adaptive-mesh hydrodynamics code (Fryxell et al. 2000; Calder et al. 2002). Modifications to the code include a nuclear burning model, composition-dependent laminar flame speeds, particular mesh refinement criteria, and instructions for determining the conditions for a DDT (for details, see Calder et al. 2007; Townsley et al. 2007, 2009). In particular, the adaptive mesh refinement (AMR) capability of the code was utilized to achieve particular resolutions for burning fronts (4 km) and the initial hydrostatic star (16 km) for which the solution has converged (Townsley et al. 2009).

The burning model is used for both subsonic (deflagration) and supersonic (detonation) burning fronts. The laminar flame width for densities and compositions characteristic of a massive C–O WD is unresolved on the scale of our grid (4 km; Chamulak et al. 2007); therefore, we use an artificially thickened flame represented by the advection–reaction–diffusion equation (Khokhlov 1995; Vladimirova et al. 2006) to which our nuclear energetics scheme is coupled. Appropriate measures have been taken to ensure this coupling is acoustically quiet and stable such that the buoyant instability of the burning front is accurately captured (Townsley et al. 2007). The shock-capturing capabilities of FLASH naturally handle the propagation of detonation fronts given an accurate nuclear energetics scheme (Meakin et al. 2009). We note that 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 an improvement to our burning model and an improvement to the method by which the transition from deflagration to detonation is made. The changes to the burning model reflect recent work to better match steady-state detonation structures that are partially resolved on the grid. This is important for obtaining accurate particle tracks for post-processing nucleosynthetic yields. Additionally, changes to the DDT method were necessary to ensure consistency between individual simulations using different DDT densities and/or different initial configurations of the flame.

3.1. Improved Burning Model

3.2. Improved Detonation Ignition Conditions

In order to study the systematic effects associated with changing the DDT density, we need to minimize any systematics in our method of starting a detonation. Previously, in Townsley et al. (2009), we visualized the simulation data from the deflagration phase and plotted a density contour at the DDT density. When the flame reached ≈64 km away from the contour, we picked a computational cell half-way between the flame front and the DDT density contour to place a detonation ignition point with a radius of 8 km. The placement and size of the detonation point was chosen to be as close to the flame front as possible while still allowing the detonation point to develop into a self-sustained, stable detonation front. If the detonation point is placed too close to the flame front, then the flame will interact with the detonation point before it develops into a self-sustained detonation. Comparisons of simulations from identical initial conditions with 8 and 12 km ignition radii finds the NSE yield differs by ≲0.5% throughout the evolution. This results indicates that the total yield is insensitive to the choice of detonation ignition radius for radii less than the characteristic size of a rising plume (as can be seen in Figure 2). We adopt 12 km for the detonation ignition radius in our study as it produces more robust detonations at low density.

Zoom In Zoom Out Reset image size

To ensure that we do not introduce unintended systematic effects in this study, we improve our method of detonation ignition point placement over the previous "by hand" method by precisely defining the criteria for a DDT that is used in an algorithm. Parameters in this algorithm are chosen to be consistent with Townsley et al. (2009). Once the flame front reaches the specified DDT density in a simulation, a detonation is ignited 32 km radially outward away from the flame front. Here the reaction–diffusion (RD) front is defined by the variable representing progress of the subsonic burning wave, ϕ_RD, with the leading edge defined as the region between the values 0.001 and 0.01 of this variable. During the deflagration phase of a simulation, ϕ_RD is equivalent to the carbon-burning reaction progress variable, ϕ_fa. At the leading edge of the RD front, very little carbon has burned and the local density is approximately equal to the unburned density. This provides a definition of the DDT density that is much more precise and accurate. We chose these criteria that ignite the detonation ahead of the RD front to avoid any issues with the detonation ignition point overlapping with the artificially thickened flame. If we were to choose criteria that would initiate a detonation inside a thickened flame, the detonation structure would need to be altered in some way to be consistent with the artificial nature of that region.

Our detonation ignition conditions also restrict detonation ignition points to be at least 200 km away from each other. This choice ensures that each rising plume starts 2–3 detonations, which is consistent with Townsley et al. (2009). In the case that multiple points within 200 km meet the detonation ignition conditions, the point furthest from the center of the star is preferred ensuring the ignition of detonations on plume "tops." In reality, the location of the spontaneous detonation points is not well known and is the subject of active research. For instance, Röpke et al. (2007) argue that a spontaneous detonation is triggered by the extreme turbulence found in the roiling fuel underneath the plume caps. Under the assumption that the DDT occurs when the flame reaches the low density for distributed burning, we place ignition points on the tops of the rising plumes. Future studies, however, will explore other physically motivated detonation methods in which the location of the detonation ignition is not necessarily specified relative to a plume, but rather determined by the local turbulence field, composition, density, and temperature (Röpke et al. 2007; Woosley et al. 2009). Figure 2 shows the placement of first one detonation point, and then another more than 200 km away from the first, both at the specified DDT density.

Each detonation ignition point is defined by the profile of the reaction progress variable representing carbon burning, ϕ_fa. Within the detonation ignition radius ϕ_fa = 1, and, because we are igniting in fuel, ϕ_fa = 0 outside this region. The subsequent energy release from this change in ϕ_fa drives a shock strong enough to ignite the surrounding unburned material. This top-hat profile is a simple approach to igniting a detonation and has some drawbacks. At densities below ≈10⁷ g cm⁻³, this detonation ignition method does not produce a sufficiently strong shock to burn the fuel. The improvements to our DDT method place the detonation point radially outward from the specified DDT density; therefore, the density within the detonation point is actually somewhat less than that specified by the DDT density. This differs from the method described in Townsley et al. (2009) in that the detonation point was placed at a density that was somewhat higher than the specified DDT density. Igniting a detonation specifying a DDT density below 10^7.1 g cm⁻³ is impossible using a top-hat profile because the detonation ignition point is actually placed ≈10⁷ g cm⁻³. In future works, we will explore the use of a gradient in ϕ_fa motivated by Seitenzahl et al. (2009a) to describe the detonation ignition point profile which should result in the formation of stronger shocks that will burn the surrounding low-density fuel.

To study the systematic effect of transition density on the ⁵⁶Ni yield synthesized during a simulated explosion, we utilize the statistical framework developed in Townsley et al. (2009). In order to compare results between this study and Townsley et al. (2009), we use the same sample population using the same initial seed. The initial seed defines the starting point to a stream of random numbers used to characterize our sample population of thermonuclear SNe as described in Townsley et al. (2009). In that study, we found that the sample dispersion in the estimated NSE yield does not asymptote until more than 10 realizations. Accordingly, our sample is made up of 30 realizations.

During the early part of the deflagration phase (≈0.1 s), the flame is most affected by convective motion in the core of the WD. Because the velocity field in the core and the number of ignition points are largely uncertain, we choose to characterize the initial flame surface using spherical harmonics, each with a random coefficient picked sequentially from the initial seed. Each realization is defined with a unique perturbation on the initial spherically symmetric flame surface using

$$\begin{equation} r\left(\theta \right) = r_0 + a_0 \sum _{l=l_{\rm min}}^{l_{\rm max}} A_l Y_l^0\left(\theta \right) {\rm,} \end{equation} \tag{ 3 }$$

where r₀ is the radius of flame surface, a₀ is the amplitude of the perturbations, A_l is a randomly chosen coefficient corresponding to the spherical harmonic Y⁰_l, and l_min and l_max set the range of spherical modes used to perturb the flame surface. This method serves to initialize RT-unstable plumes of random relative strengths, which we would expect from various distributions of ignition points and varying strengths of the convective velocity field found in the real population of progenitors. In this study, r₀ = 150 km, a₀ = 30 km, l_min = 12, and l_max = 16. The choices of these parameters are motivated by the resolution of our study and the desire to obtain reasonable ⁵⁶Ni yields inferred from observations. These choices are discussed in further detail in Townsley et al. (2009).

We choose to analyze the dependence on transition density by choosing five different transition densities equidistant in log space at $\mathcal {L}\rho _{\rm DDT}= \lbrace 7.1, 7.2, 7.3, 7.4, 7.5\rbrace$ in cgs units, where $\mathcal {L}\rho _{\rm DDT}= \log _{10}{\rho _{\rm DDT}}$. The range of DDT densities chosen for this study is motivated both by previous work and computational challenges. Below $\mathcal {L}\rho _{\rm DDT}= 7.1$, a more realistic detonation structure is needed to successfully launch a detonation wave as described in the previous section and in Seitenzahl et al. (2009a). While recent studies have suggested a wider range of DDT densities are possible (Pan et al. 2008; Schmidt et al. 2010), we suggest that trends resulting from varying the DDT density are captured with a maximum $\mathcal {L}\rho _{\rm DDT}= 7.5$.

A simulation is performed for each of the 30 realizations in our sample at each DDT density for a total of 150 simulations. We choose not to explore DDT densities below this specified range for computational reasons. Because of the approximate power-law density profile of the WD, DDT densities were chosen in log space because these densities correspond to relatively evenly spaced radial coordinates of the WD. The amount of ⁵⁶Ni synthesized in the explosion principally depends on the amount of expansion during the deflagration phase. Therefore, spatially separated detonation ignition conditions will allow the amount of expansion to vary a non-negligible amount and we can more easily analyze the dependence of the yield on DDT density.

We find that extrapolating our results to a DDT density of $\mathcal {L}\rho _{\rm DDT}= 6.83$ yields an estimated average ⁵⁶Ni yield of ≃0.60 M_☉ with a standard deviation of 0.21 M_☉ (from the fits described below and listed in Table 4) that is consistent with observations (Howell et al. 2009). Because the actual DDT density is unknown and the subject of ongoing research (Pan et al. 2008; Aspden et al. 2008, 2010; Woosley et al. 2009; Schmidt et al. 2010), we are free to choose this value of the DDT density as the fiducial DDT density, $\mathcal {L}\rho _{{\rm DDT},0}= 6.83$. This choice is relevant for analysis and comparison to other works as discussed in Section 5. The distribution of ⁵⁶Ni material synthesized during the explosion is shown in Figure 3 for different transition densities at $\mathcal {L}\rho _{\rm DDT}= 6.83, 7.10$, and 7.30.

Zoom In Zoom Out Reset image size

The NSE mass, M_NSE, is defined as ∫ϕ_qnρdV integrated over the star. We determine the NSE yield by running the simulation until M_NSE is no longer increasing as a function of simulation time. This condition is defined as

$$\begin{equation} \frac{d M_{\rm NSE}}{d t} < 0.01 \frac{M_\odot }{s}{\rm.} \end{equation} \tag{ 4 }$$

Because we estimate the ⁵⁶Ni yield as a fraction of the NSE yield, we consider the ⁵⁶Ni yield to have plateaued when the NSE yield has plateaued. Considering additional ⁵⁶Ni from incomplete Si burning in NSQE material and the efficient capture of excess neutrons onto Fe-group elements changes the final ⁵⁶Ni estimate by ≲1%. Therefore, to good approximation, the ⁵⁶Ni yield is a fraction of the NSE yield. Figure 4 shows the evolution of the NSE yield as a function of simulation time showing the DDT time and the NSE yield plateau time for each DDT density for realization 2. Discernible from this figure, there is a clear dependence of the NSE yield on transition density.

Zoom In Zoom Out Reset image size

Within the DDT paradigm for thermonuclear SNe, the conditions under which a transition occurs are still not completely understood (Aspden et al. 2008, 2010). Generally, though, a hypothesized transition occurs when the flame enters the distributed burning regime, which occurs when the flame speed equals the turbulent velocity at the scale of the flame width assuming a Kolmogorov turbulent cascade (Niemeyer & Woosley 1997). Recent work has placed more stringent requirements on the DDT (Pan et al. 2008; Woosley et al. 2009; Schmidt et al. 2010) and fundamental questions concerning deflagration in the limit of disruptive turbulence remain (Poludnenko & Oran 2010). Regardless of the actual DDT mechanism, it is likely that the DDT conditions will depend on composition because both the width and burning rate of the flame depend on the abundances ¹²C and ²²Ne. For the current study, we assume that a DDT will occur at a unique density given a particular composition, regardless of the microphysics involved. This assumption is reasonable given that the characteristics of the flame depend strongly on density. Using this assumption, we can delay the analysis of the particular microphysics that lead to a specific transition density and analyze the dependence of the amount of material synthesized to NSE during the explosion on the transition density. Therefore, we can think of each transition density as a proxy for changing the composition that determines the conditions for a DDT via the appropriate microphysics.

As discussed in detail in Townsley et al. (2009), many systematic effects exist that influence the outcome of an SN. In that study, the direct effect of ²²Ne was explored and found to have a negligible influence on the NSE yield. Therefore, we do not vary the initial ²²Ne mass fraction, but rather study the effect of varying the DDT density with the expectation that the NSE yield will be influenced indirectly by ²²Ne and ¹²C abundances through the DDT density. We focus on the indirect effect of the ²²Ne abundance on the NSE yield, following up previous work in Townsley et al. (2009) and neglect the effect of varying the carbon abundance. Additionally, we neglect effects due to the central ignition density, compositional and thermodynamic WD structure, and the total WD mass. The effects due to these variables will be studied in turn in future works with the goal of addressing interdependencies once individual effects are better understood.

5.1. Dependence on Transition Density

The evolution of the amount of material above a density threshold ($M_{\rho > \rho _{\rm thres}} = \int _{\rho _{\rm thres}} \rho dV$, where we take ρ_thres = 2 × 10⁷ g cm⁻³) principally determines the dependence of the NSE yield on DDT density. This is due largely to the linear relationship between the NSE yield and $M_{\rho _7>2}(t=t_{\rm DDT})$ shown in Figure 5, where t_DDT is defined as the time the flame first reaches the specified DDT density. This definition is consistent with our assumption that DDT conditions are met on the tops of rising plumes. As mentioned above, other DDT locations are possible, such as the highly turbulent region underneath plume "caps," which would result in a higher DDT density correlated to our present definition via the density structure local to the rising plume. As shown in Figure 4, the deflagration phase burns only a relatively small fraction of the WD and the majority of material is burned during the detonation. Once a detonation has started, the propagation speed of the burning wave is much greater than the rate of expansion; therefore, the NSE yield is essentially independent of the evolution of $M_{\rho _7>2}$ for t>t_DDT. The number of detonation points and their corresponding distribution in time and location contributes to the variance of the relation between the NSE yield and $M_{\rho _7>2}$. The relation between the NSE yield and DDT density via $M_{\rho _7>2}$ also depends on the acceleration of the RT-unstable plumes not being too great near t_DDT. For a constant plume rise rate near t_DDT, there is a linear relationship between the DDT density and t_DDT. This relationship allows the evolution of $M_{\rho _7>2}$ to translate directly into a dependence of the NSE yield on DDT density. Figure 6 shows the evolution of $M_{\rho _7>2}$ and ρ_min as a function of simulation time for three different realizations where ρ_min is the minimum density at which the flame is burning material. The evolution of ρ_min shows the linear relationship between the log of the DDT density and t_DDT where the filled circles represent the t_DDT times for each of the five transition densities and the filled squares show the fiducial DDT density. For the times corresponding to the DDT conditions, the evolution of $M_{\rho _7>2}$ is in a region where the rate of expansion of the star becomes significant and $M_{\rho _7>2}$ begins to drop off relatively quickly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find, for each individual realization, the NSE yield depends quadratically on transition density. Table 3 lists the NSE yields for each realization at each of five transition densities equidistant in log space. Table 4 shows the coefficients used to fit the quadratic dependence of NSE yield on transition density for

$$\begin{equation} M_{\rm NSE} (M_\odot) = a \left(\mathcal {L}\rho _{\rm DDT}\right)^2 + b \mathcal {L}\rho _{\rm DDT}+ c {\rm.} \end{equation} \tag{ 5 }$$

Figure 7 shows fits to the NSE yield for each realization and the average yields at each transition density. The NSE yields for individual realizations are plotted as blue crosses and the individual curves showing the dependence of NSE yield on DDT density are shown in gray. The average NSE yields at each DDT density are shown in red. At a transition density of 10^7.5 g cm⁻³, notice that the dependence on DDT density has flattened out for almost all realizations and the variance of the yields among all realizations is small. This behavior is due to the fact that the progenitor WD star has a finite mass of 1.37 M_☉. At the highest transition density, the expansion rate of the star is relatively small, and the evolution of $M_{\rho _7>2}$ is relatively slowly declining as compared to lower DDT densities (see Figure 6).

Zoom In Zoom Out Reset image size

Table 3. NSE Yields in M_☉ for each Realization at each Transition Density, $\mathcal {L}\rho _{\rm DDT}$

Rel.	$\mathcal {L}\rho _{\rm DDT}$
#No.	7.1	7.2	7.3	7.4	7.5
1	0.989	1.082	1.152	1.207	1.236
2	0.970	1.056	1.131	1.184	1.215
3	0.948	1.043	1.119	1.174	1.209
4	0.951	1.062	1.144	1.197	1.226
5	1.027	1.101	1.168	1.214	1.242
6	1.074	1.138	1.185	1.219	1.248
7	1.076	1.132	1.179	1.215	1.243
8	1.110	1.167	1.210	1.239	1.260
9	0.917	1.022	1.105	1.163	1.199
10	0.791	0.929	1.038	1.113	1.166
11	1.092	1.147	1.189	1.229	1.255
12	1.240	1.258	1.273	1.289	1.301
13	0.882	0.987	1.083	1.145	1.183
14	1.070	1.125	1.185	1.223	1.250
15	1.119	1.176	1.218	1.254	1.274
16	0.883	0.994	1.081	1.142	1.183
17	1.026	1.115	1.171	1.213	1.241
18	1.181	1.223	1.249	1.270	1.290
19	0.967	1.057	1.132	1.185	1.220
20	1.005	1.089	1.138	1.180	1.218
21	1.222	1.248	1.268	1.284	1.297
22	1.045	1.104	1.170	1.216	1.244
23	1.001	1.085	1.150	1.187	1.223
24	0.900	1.012	1.091	1.130	1.176
25	0.989	1.057	1.139	1.193	1.226
26	1.092	1.150	1.192	1.228	1.254
27	1.192	1.222	1.251	1.274	1.290
28	1.145	1.183	1.220	1.253	1.273
29	1.063	1.131	1.175	1.210	1.235
30	1.179	1.225	1.249	1.273	1.288
$\bar{M}$	1.038	1.111	1.169	1.210	1.239
σ	0.110	0.082	0.059	0.046	0.036
$\sigma _{\bar{M}}$	0.020	0.015	0.011	0.008	0.007

Download table as:  ASCIITypeset image

Table 4. Coefficients used to Fit the Quadratic Dependence of Transition Density on NSE Yield Using Equation (5)

No.	a (M☉)	b (M☉)	c (M☉)	$M_{\rm NSE}(\mathcal {L}\rho _{{\rm DDT},0})$
1	−1.024	15.56	−57.9	0.637
2	−0.933	14.24	−53.1	0.633
3	−1.010	15.40	−57.5	0.589
4	−1.388	20.95	−77.8	0.515
5	−0.815	12.43	−46.2	0.732
6	−0.595	9.12	−33.7	0.852
7	−0.480	7.42	−27.4	0.876
8	−0.619	9.41	−34.5	0.898
9	−1.164	17.70	−66.1	0.515
10	−1.454	22.17	−83.3	0.279
11	−0.442	6.87	−25.4	0.902
12	−0.083	1.36	−4.3	1.183
13	−1.212	18.45	−69.0	0.456
14	−0.557	8.58	−31.8	0.844
15	−0.575	8.78	−32.2	0.910
16	−1.173	17.88	−66.9	0.470
17	−0.963	14.58	−54.0	0.712
18	−0.348	5.34	−19.2	1.048
19	−0.926	14.15	−52.8	0.628
20	−0.708	10.86	−40.4	0.740
21	−0.210	3.25	−11.3	1.134
22	−0.583	9.03	−33.6	0.798
23	−1.093	16.47	−60.9	0.696
24	−1.219	18.46	−68.7	0.502
25	−0.782	12.04	−45.1	0.671
26	−0.490	7.56	−27.9	0.896
27	−0.248	3.86	−13.7	1.080
28	−0.276	4.36	−15.9	1.006
29	−0.685	10.42	−38.4	0.827
30	−0.429	6.54	−23.6	1.032

Note. The NSE mass in units of M_☉ is evaluated at $\mathcal {L}\rho _{{\rm DDT},0}$ using the coefficients.

Download table as:  ASCIITypeset image

The curvature of the NSE yield dependence on DDT density, a, is well correlated with the NSE yield at a given DDT density. Figure 8 shows this correlation for $\mathcal {L}\rho _{{\rm DDT},0}=6.83$; however, a correlation exists for all DDT densities as may be discernible from Figure 7 noting that most black lines representing individual realizations do not cross. The lower the NSE yield for a given realization, the stronger the dependence on transition density. This result is likely due to realizations with lower yields having multiple competing plumes that release more energy allowing the star to expand more rapidly leading to a stronger dependence on the transition density.

Zoom In Zoom Out Reset image size

After fitting the dependence on DDT density for each realization, we find an interesting correlation between the fitting parameters from Equation (5). Figure 9 shows the correlation between the fitting parameters a, b, and c with c on the horizontal axis. The tight correlation between these parameters indicates that a single parameter describes the dependence on DDT density for a given realization. The fitting parameters a and b can be expressed as a function of c:

$$\begin{eqnarray} a &=& \alpha c + \beta, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} b &=& \delta c + \gamma {\rm,} \end{eqnarray} \tag{ 7 }$$

where α = (1.754 ± 0.003) × 10⁻², β = −0.005 ± 0.002, δ = −0.2646 ±  0.0005, and γ = 0.21 ±  0.02. These parameters were calculated using a least-squares method and the associated errors were calculated from the corresponding covariance matrix.

Zoom In Zoom Out Reset image size

Each realization has a different random initial perturbation of the central ignition condition that sets the plume morphology. Dominant single plumes tend to allow the star to expand less prior to reaching the conditions for a detonation, while multiple competing plumes tend to release more energy during the deflagration phase, allowing more expansion prior to reaching the conditions for a DDT. Shown in Figure 10 are the RT unstable plumes for the realizations with the highest and lowest yields a few tenths of seconds into the deflagration. Unfortunately, no strong correlation exists between the properties of the initial flame surface, such as mass enclosed, surface to volume ratio, or amount of power in the perturbation, for a given realization and the single fitting parameter c. This lack of a correlation implies there is no way to tell whether a particular initial condition will seed a single dominant plume or multiple competing plumes.

Zoom In Zoom Out Reset image size

Physically, we might expect two competing effects that influence the NSE yield: plume morphology (rise time) and the rate of expansion. In our two-dimensional simulations, these two effects appear to be correlated as seen in Figure 6. Our results indicate that plumes near the symmetry axis tend to rise faster than plumes rising near the equator and we attribute this result to the fact that our simulations are two-dimensional. A plume that develops near the symmetry axis naturally becomes dominant in two dimensions and determines the DDT time and thus the NSE yield. In addition, a plume near the symmetry axis represents less volume than a plume of similar size near the equator. This indicates that energy is being deposited into a smaller volume allowing for a lesser rate of expansion. This explains why the overall rate of expansion of the star is correlated to the rise time of the first plume to reach DDT conditions.

In order to evaluate the dependence on DDT density without this unphysical correlation, we need to perform three-dimensional simulations. A suite of three-dimensional simulations will likely create a two-parameter family (plume rise and expansion rates) of solutions which describe the yield as a function DDT density. Additionally, we expect three-dimensional plumes to behave similarly to two-dimensional on-axis plumes implying a faster plume rise time.

As a result of the tight correlation for a and b and the determination of Equations (6) and (7), the average relation is characterized by just the average c—or equivalently, the average NSE yield at a specified ρ_DDT. We evaluate the average and standard deviation of the fitting parameter c as well as the statistical properties of the NSE yield at the fiducial transition density of 10^6.83 g cm⁻³ shown in Table 5.

Table 5. Statistical Properties of the Fitting Parameter c and the NSE Yield at ρ_DDT = 6.76 × 10⁶ g cm⁻³

Parameter	Mean	Std. Dev.	Std. Dev. of Mean
c (M☉)	−42	21	4
MNSE,0 (M☉)	0.77	0.22	0.04

Download table as:  ASCIITypeset image

5.2. NSE Yield Dependence on ²²Ne

Recall that for the purposes of this study, we are neglecting the effects of varying the core carbon abundance, the central ignition density, the WD structure, etc., except for the indirect effect of X₂₂, the neutron excess produced during simmering, and the DDT density. In order to derive the NSE yield dependence on X₂₂, we need to expand the derivative of M_NSE with respect to X₂₂ which involves only a couple of terms given our assumptions

$$\begin{equation} \frac{d M_{\rm NSE}}{d X_{22}} = \frac{\partial M_{\rm NSE}}{\partial \mathcal {L}\rho _{\rm DDT}} \frac{d \mathcal {L}\rho _{\rm DDT}}{d X_{22}} + \frac{\partial M_{\rm NSE}}{\partial X_{22}}{\rm.} \end{equation} \tag{ 8 }$$

Townsley et al. (2009) showed that $\frac{\partial M_{\rm NSE}}{\partial X_{22}} = 0$. In that study, they employed an estimate of the mass burned to NSE by measuring the amount of mass above a density of 2 × 10⁷ g cm⁻³. That correlation is confirmed by our current study as shown in Figure 5 where we plot the NSE yield of all realizations at each transition density. The correlation between NSE yield and mass above 2 × 10⁷ g cm⁻³ at t_DDT for a transition density of 1.26 × 10⁷ g cm⁻³, closest to the transition density used for that study, is highlighted. Given our assumptions and the result that the NSE yield does not depend directly on X'₂₂ (Townsley et al. 2009), the NSE yield only depends directly on the DDT density. The effect of X'₂₂ on the yield enters through the DDT density. Therefore, we will construct the functional dependence of $\mathcal {L}\rho _{\rm DDT}$ on X'₂₂.

The DDT density depends on X₂₂ via the microphysics involved in determining the conditions under which the flame transitions from a deflagration to a detonation, which, as noted above, is incompletely understood. We assume the transition occurs at a particular density at which the flame enters the distributed burning regime, but more stringent conditions for the DDT include dependences on the turbulent cascade and the growth of a critical mass of fuel with sufficiently strong turbulence (Pan et al. 2008; Woosley et al. 2009; Schmidt et al. 2010). The dependence on composition of these models may be explored and applied to the trends with DDT density presented in this study. Under our assumptions, the flame enters the distributed burning regime when the laminar flame width becomes of the order of the Gibson length (l_G), where

$$\begin{equation} l_{\rm G} = L \left(\frac{s_l}{u^{\prime }} \right)^3 {\rm.} \end{equation} \tag{ 9 }$$

Here, L is the length scale on which the strength of the turbulent velocity (u') is evaluated, and s_l is the laminar flame speed (see, e.g., Peters 2000). By using the laminar flame speeds and widths which depend on X₂₂ from Chamulak et al. (2007), we choose the conditions for a detonation are met at X₂₂ = 0.02 at a particular transition density ($\mathcal {L}\rho _{\rm DDT}$) when the Gibson length is equal to the laminar flame width. The strength of the turbulent velocity field is calculated assuming a Kolmogorov turbulence cascade:

$$\begin{equation} u^{\prime } = s_l(\rho,X_{22}=0.02) \left(\frac{\delta _l(\rho,X_{22}=0.02)}{L} \right)^{1/3} {\rm.} \end{equation} \tag{ 10 }$$

Using u', we solve for the change in transition density, $\Delta \mathcal {L}\rho _{\rm DDT}$, by changing ²²Ne to X₂₂ = 0.06 and again setting the Gibson length equal to the laminar flame width where L cancels out of the equation. This evaluation yields a change in transition density over a change in X₂₂, or $\frac{d \mathcal {L}\rho _{\rm DDT}}{d X_{22}}$.

We use log–log fits to the flame speed and flame width as a function of density using the table generated by Chamulak et al. (2007). We use only densities above 10⁸ g cm⁻³ for X₁₂ = 0.3 and above 5 × 10⁸ for X₁₂ = 0.5. Within this parameter space, a power-law dependence of the flame speeds and widths on the log of the density is well defined and the algorithm used to solve the flame characteristics is more stable at higher densities. Our resulting expressions for the density dependence on flame speed, s, and flame width, δ_l, at the carbon mass fraction used in this study (X₁₂ = 0.4) for X₂₂ = 0.02, 0.06 are given by

$$\begin{equation} \ln {s} = a \ln {\rho } + b \end{equation} \tag{ 11 }$$

with coefficients given in Table 6.

Table 6. Coefficients for log–log Fits to Equation (11)

Parameter	X22 = 0.02	X22 = 0.06
as	0.7942	0.7745
bs	−12.735	−12.121
aδ	−1.3550	−1.3507
bδ	19.17	18.88

Note. The subscript denotes whether fitting for the laminar flame speed (s) or flame width (δ).

Download table as:  ASCIITypeset image

Our derived expression for the derivative of transition density as a function of X₂₂ is given by

$$\begin{eqnarray} \frac{d \mathcal {L}\rho _{\rm DDT}}{d X_{22}} &=& \frac{ b_{\delta,6} - b_{\delta,2} + 3\, (b_{s,2} - b_{s,6}) + \ln {\rho _{\rm DDT}}\, (3 a_{s,2} - a_{\delta,2}) }{ \ln {10}\, (3 a_{s,6} - a_{\delta,6}) \Delta X_{22}} \nonumber \\ &&- \frac{\mathcal {L}\rho _{\rm DDT}}{\Delta X_{22}} \nonumber \\ &=& u \mathcal {L}\rho _{\rm DDT}+ v {\rm,} \end{eqnarray} \tag{ 12 }$$

where u = 0.4315 and v = −6.301. We solve this first-order differential equation and express the DDT density as a function of X'₂₂:

$$\begin{equation} \mathcal {L}\rho _{\rm DDT}(X_{22}^\prime) = \frac{v}{u} \big(e^{u (X_{22}^\prime - X_{22,0}^\prime)} - 1 \big) + \mathcal {L}\rho _{{\rm DDT},0}e^{ u (X_{22}^\prime - X_{22,0}^\prime)} {\rm,} \end{equation} \tag{ 13 }$$

where X'_22,0 is the parameterized ²²Ne mass fraction chosen to be associated with the fiducial transition density, $\mathcal {L}\rho _{{\rm DDT},0}$. For all the transition densities considered in this study, the interface between the core composition and the composition of the outer layer is at a lower density. Therefore, the relevant ²²Ne content to consider is that of the core. Plugging in Equation (13) into Equation (5), we obtain the functional dependence of M_NSE on X'₂₂, such that

$$\begin{equation} M_{\rm NSE} = M_{\rm NSE}(\mathcal {L}\rho _{\rm DDT}(X_{22}^\prime)) {\rm.} \end{equation} \tag{ 14 }$$

Equation (14) is evaluated and plotted in Figure 11 for the fiducial transition density at the ²²Ne mass fraction used in this study,

$$\begin{equation} \mathcal {L}\rho _{\rm DDT}(X_{22,0}^\prime =0.03) = \mathcal {L}\rho _{{\rm DDT},0}= 6.83 {\rm.} \end{equation} \tag{ 15 }$$

We propagate the standard deviation of the mean evaluated at $\mathcal {L}\rho _{{\rm DDT},0}$ for a range of X'₂₂. This is calculated by considering the relation between the standard deviation of the mean of the NSE mass and the standard deviation of the mean of the fitting parameter, c, given by

$$\begin{equation} \sigma _{\rm NSE}=(\alpha [\mathcal {L}\rho _{\rm DDT}(X_{22}^\prime)]^2 + \delta \mathcal {L}\rho _{\rm DDT}(X_{22}^\prime) + 1) \sigma _c {\rm,} \end{equation} \tag{ 16 }$$

where σ_NSE is the standard deviation of the mean of the NSE mass and σ_c is the standard deviation of the mean in the fitting parameter c. We evaluate σ_c by inverting Equation (16) and solving for X'₂₂ = X'_22,0. Plugging in this solution, the standard deviation of the mean of the NSE mass as a function of X'₂₂ becomes

$$\begin{equation} \sigma _{\rm NSE} (X_{22}^\prime) = \frac{\alpha [\mathcal {L}\rho _{\rm DDT}(X_{22}^\prime)]^2 + \delta \mathcal {L}\rho _{\rm DDT}(X_{22}^\prime) + 1}{\alpha [\mathcal {L}\rho _{{\rm DDT},0}]^2 + \delta \mathcal {L}\rho _{{\rm DDT},0}+ 1} \sigma _{{\rm NSE},0} {\rm.} \end{equation} \tag{ 17 }$$

Zoom In Zoom Out Reset image size

5.3. ⁵⁶Ni Yield Dependence on Metallicity

Now that we have constructed the functional dependence of NSE yield on ²²Ne we need to consider the dependence of the amount of ⁵⁶Ni synthesized in the explosion on metallicity through the ²²Ne content. A fractional amount of the NSE material is radioactive ⁵⁶Ni, which powers the SN light curve. This fraction depends on Y_e (and ²²Ne) as described in Timmes et al. (2003):

$$\begin{equation} M_{56} = f_{56} M_{\rm NSE} {\rm.} \end{equation} \tag{ 18 }$$

In order to determine the dependence of M₅₆ on metallicity through the DDT density, we must construct f₅₆ and explore any dependencies on DDT density and metallicity. From our simulations, we calculate M_NSE directly and estimate M₅₆ from Y_e choosing the non-⁵⁶Ni NSE material to be 50/50 ⁵⁴Fe and ⁵⁸Ni by mass. Recall that X'₂₂ from Equation (1) contains a term, ΔX₂₂, which is a parameterization of the change in Y_e due to neutronization during the carbon-simmering phase. However, not all of the neutron excess evaluated at the M_NSE plateau time comes from ²²Ne or the carbon-simmering products. Some change in Y_e is due to weak reactions occurring during the explosion that are included in our burning model (Calder et al. 2007; Townsley et al. 2007).

First, we consider the amount of non-⁵⁶Ni NSE material determined by X'₂₂ by equating the initial Y_e to the Y_e of material in NSE. Using baryon and lepton conservation for NSE material, we describe the electron fraction by contributions from ⁵⁴Fe, ⁵⁸Ni, and ⁵⁶Ni

$$\begin{equation} Y_{e}= \frac{26}{54}\left(\frac{1}{2}f_{{\rm non}-56}\right) + \frac{28}{58}\left(\frac{1}{2}f_{{\rm non}-56}\right) + \frac{28}{56}\left(1-f_{{\rm non}-56}\right) {\rm,} \end{equation} \tag{ 19 }$$

where f_non−56 is the mass fraction of non-⁵⁶Ni NSE material. For the following evaluation, we approximate the composition to be that of the core because most of the NSE material is within the core. We write the initial Y_e as

$$\begin{equation} Y_{e,i} = \frac{10}{22}X_{22}^\prime + \frac{1}{2}(1-X_{22}^\prime) {\rm.} \end{equation} \tag{ 20 }$$

We can then solve for f_non−56 by equating Equations (19) and (20) and add a term X_n to represent the additional contribution due to neutronization from weak reactions occurring during the explosion

$$\begin{equation} f_{{\rm non}-56} = \frac{783}{308}X_{22}^\prime + X_n {\rm.} \end{equation} \tag{ 21 }$$

Then, the ⁵⁶Ni fraction of material in NSE is

$$\begin{equation} f_{56} = 1 - \frac{783}{308}X_{22}^\prime - X_n {\rm.} \end{equation} \tag{ 22 }$$

We note that the rate of weak reactions occurring during the explosion may depend on the initial composition as well. Accordingly, X_n may have a dependence on X'₂₂. We also expect X_n to vary as a function of the transition density because a higher transition density will have less time for weak reactions to occur.

We construct a statistical sample of f₅₆ using the ratio of the ⁵⁶Ni and NSE yields produced in the simulations. We calculate the dependence of X_n on transition density using a least-squares method and the result is shown in Figure 12. Evaluating the partial derivative of X_n at the fiducial transition density, we find a shallow dependence:

$$\begin{equation} \left.\frac{\partial X_n}{\partial \mathcal {L}\rho _{\rm DDT}}\right|_{\mathcal {L}\rho _{\rm DDT}\,=\,\mathcal {L}\rho _{{\rm DDT},0}} = -0.096 {\rm.} \end{equation} \tag{ 23 }$$

Zoom In Zoom Out Reset image size

We want to know the dependence of M₅₆ on X₂₂. Since f₅₆ depends on the neutronization from weak reactions during the explosion, we must consider any dependencies of this neutronization on DDT density or X'₂₂. We show that the dependence of f₅₆ on X'₂₂ via the effect of X'₂₂ on the in situ neutronization is much weaker than the direct dependence due to lepton number conservation. Similar to Equation (18), we write

$$\begin{equation} \frac{d M_{56}}{d X_{22}} = \frac{\partial M_{56}}{\partial M_{\rm NSE}} \frac{d M_{\rm NSE}}{d X_{22}} + \frac{\partial M_{56}}{\partial X_{22}} {\rm.} \end{equation} \tag{ 24 }$$

Using Equation (18) and expanding the partial derivative of M₅₆ on X₂₂, we obtain

$$\begin{equation} \frac{d M_{56}}{d X_{22}} = f_{56} \frac{d M_{\rm NSE}}{d X_{22}} + M_{\rm NSE}\frac{d f_{56}}{d X_{22}} {\rm.} \end{equation} \tag{ 25 }$$

Taking the derivative of Equation (22) and expanding on X₂₂, we obtain

$$\begin{equation} \frac{d f_{56}}{d X_{22}} = \frac{\partial f_{56}}{\partial X_{22}} - \frac{d X_n}{d X_{22}} {\rm.} \end{equation} \tag{ 26 }$$

Now we wish to evaluate whether $\frac{d X_n}{d X_{22}}$ is an important contribution to the overall evaluation of the ⁵⁶Ni mass. Expanding this term yields

$$\begin{equation} \frac{d X_n}{d X_{22}} = \frac{\partial X_n}{\partial \mathcal {L}\rho _{\rm DDT}} \frac{d \mathcal {L}\rho _{\rm DDT}}{d X_{22}} + \frac{\partial X_n}{\partial X_{22}} {\rm.} \end{equation} \tag{ 27 }$$

Referring to Townsley et al. (2009), we can estimate $\frac{\partial X_n}{\partial X_{22}}$ by calculating the average ratio of M₅₆ to M_NSE at X₂₂ = 0, 0.02 for the first five realizations whose detonation phases were simulated. The result is $\frac{\partial X_n}{\partial X_{22}} \sim -0.2$. The first term in Equation (27) can be evaluated from multiplying Equation (23) and Equation (12) for $\mathcal {L}\rho _{\rm DDT}$ in the range 7.0–7.5. For the fiducial transition density ($\mathcal {L}\rho _{{\rm DDT},0}$), we find

$$\begin{equation} \frac{d X_n}{d X_{22}} = \left(-0.126 \right) \left(-3.35 \right) - 0.2 \sim 0.2 {\rm.} \end{equation} \tag{ 28 }$$

Comparing to $\frac{\partial f_{56}}{\partial X_{22}} \simeq -2.5$, we find that the magnitude of $\frac{d X_n}{d X_{22}}$ is much smaller and is unimportant for our evaluation of the dependence of the ⁵⁶Ni yield on X₂₂. Therefore, we can ignore this term in the expansion of $\frac{d f_{56}}{d X_{22}}$ in Equation (26) and let $\frac{d X_n}{d X_{22}} \sim 0$. The full derivative of f₅₆ with respect to X₂₂ can now be written as a partial derivative such that

$$\begin{equation} \frac{d M_{56}}{d X_{22}} = f_{56} \frac{d M_{\rm NSE}}{d X_{22}} + M_{\rm NSE} \frac{\partial f_{56}}{\partial X_{22}} {\rm.} \end{equation} \tag{ 29 }$$

While we approximate X_n as constant, we choose to evaluate it at $\mathcal {L}\rho _{{\rm DDT},0}$ using the best-fit curve from Figure 12 obtaining X_n = 0.133. Now we can relate the metallicity to ²²Ne since X₂₂ traces metallicity (Timmes et al. 2003). Substituting X₂₂ = 0.014(Z/Z_☉) in Equation (1), we obtain

$$\begin{equation} X_{22}^\prime = 0.014 \left(\frac{Z}{Z_\odot } \right) + \Delta X_{22}(\Delta Y_{e}) {\rm,} \end{equation} \tag{ 30 }$$

where ΔX₂₂ = 0.01 for this study and X'₂₂ is our parameterization of the actual ²²Ne mass fraction, X₂₂. The ⁵⁶Ni yield as a function of X'₂₂ is calculated by multiplying Equation (22) by Equation (14) and using X_n = 0.133. The ⁵⁶Ni yield is plotted in Figure 13 using Equation (30) to relate X'₂₂ to Z/Z_☉.

Zoom In Zoom Out Reset image size

We have analyzed the influence of the DDT density on the total ⁵⁶Ni synthesized during thermonuclear SNe. We determined that the dependence on DDT density is quadratic in nature, but this dependence can be described by a single parameter. We estimated the dependence of the SN brightness (⁵⁶Ni yield) on metallicity by assuming a DDT occurs when the flame enters the distributed burning regime and extrapolating the laminar flame speeds and widths down to low densities. We find

$$\begin{equation} \left.\frac{d M_{56}}{d (Z/Z_\odot)}\right|_{Z\,=\,Z_\odot } = -0.067 \pm 0.004 \,M_\odot {,} \end{equation} \tag{ 31 }$$

which is slightly steeper than Timmes et al. (2003) as seen in Figure 13. The uncertainty was calculated using the standard deviation of the mean of the fitting parameter c—or equivalently, the standard deviation of the mean of the NSE yield at a particular DDT density. The uncertainty in the assumptions about the normalization of the transition density as a function of X'₂₂ was not considered for the purpose of evaluating the uncertainty in the derivative. We find the effect of metallicity on DDT density influences the production of NSE material; however, the ratio of the ⁵⁶Ni yield to the overall NSE yield does not change as significantly and remains similar to the relation estimated from approximate lepton number conservation. We also find that the scatter in SN brightness increases with decreasing transition density.

The very recent work of Bravo et al. (2010) on metallicity as a source of dispersion in the luminosity–width relationship of bolometric light curves addresses many of the same issues as our study and warrants discussion. In particular, Bravo et al. (2010) derive a metallicity dependence on the DDT density similarly to our study. However, the DDT density in their one-dimensional simulations is very different from the DDT density in ours. The principal difference between the work described in this manuscript and that of Bravo et al. (2010) follows from the use of multidimensional simulations. While three-dimensional simulations are required to correctly capture the effects of fluid dynamics and the RT instability, two-dimensional simulations incorporate these effects and relax the assumption of symmetry. Breaking the symmetry and assuming that a DDT initiates as a rising plume approaches ρ_DDT produce an expansion history very different from what would be observed in one-dimensional simulations. In fact, using one-dimensional simulations implies that ρ_DDT plays an unphysical role. The real physical description of the DDT in an SN Ia will depend on flame–turbulence interactions that will themselves depend on multidimensional effects in the flow. One-dimensional models may be able to capture these effects, but such models must be motivated by more physical multidimensional studies.

In addition, our statistical framework with realizations from randomized initial conditions allows calculation of the scatter inherent in the multidimensional models. The standard deviations calculated for our models averaged over a set of realizations demonstrate that there can be considerable variation following from the randomized initial conditions. These variations follow from the different amounts of expansion occurring during the deflagration phase that follow from the different rising plume morphologies. In one-dimensional models, even with the progenitor metallicity determining ρ_DDT, there is a one-to-one correspondence between expansion and ρ_DDT for a given progenitor.

Bravo et al. (2010) explored two scenarios, a linear and a nonlinear dependence of the ⁵⁶Ni yield on Z, and report excellent agreement between the nonlinear scenario and observations reported by Gallagher et al. (2008) with the caveat that they used one-dimensional models to arrive at their conclusion. Without comparing results involving the metallicity dependence on ρ_DDT, Bravo et al. (2010) find a stronger dependence of the ⁵⁶Ni synthesized from incomplete Si burning on Z than ⁵⁶Ni coming from NSE material. We find a shallower dependence of the ⁵⁶Ni yield on Z than Bravo et al. (2010) in part due to our assumption that ⁵⁶Ni is synthesized as a fraction of NSE material; however, comparing their stratified models to the Timmes et al. (2003) relation results in only a ∼20% difference in the dependence of ⁵⁶Ni yield on Z. This effect is compounded in their "nonlinear" scenario, but again, we must emphasize that ρ_DDT plays a completely unphysical role in their simulations and the degree to which this effect is actually enhanced is unknown.

The light curve width calculations and subsequent population synthesis performed by Bravo et al. (2010) after finding a dependence of the ⁵⁶Ni yield on metallicity are useful in estimating a metallicity dependence in the Hubble residual. For this purpose, we compare the results of the nonlinear scenario of Bravo et al. (2010) to our study along with the expected dependence due to lepton number conservation (Timmes et al. 2003) in Figure 13. Bravo et al. (2010) find an ∼30% steeper dependence of the ⁵⁶Ni yield on Z/Z_☉ evaluated at Z = 1Z_☉ than our result in Equation (31). In addition, they find that for high Z the ⁵⁶Ni yield tends to flatten out becoming more similar to our results. For a particular subrange of metallicites, our results are very similar indicating that by performing the same analysis as Bravo et al. (2010) our results should also agree with the metallicity dependence found by Gallagher et al. (2008).

This similarity depends on the choice of mean ⁵⁶Ni yield as both Timmes et al. (2003) and Bravo et al. (2010) relations are proportional to this quantity. This implies that the steepness of the relations is affected by the choice of mean ⁵⁶Ni yield. Our results are not sensitive to this choice. Choosing a higher fiducial DDT density of 10⁷ g cm⁻³ that results in a higher mean ⁵⁶Ni yield of 0.77 M_☉ increases Equation (31) by only 0.005 M_☉. We see less dependence on the mean ⁵⁶Ni yield because, although we have a similar dependence on X₂₂ due to lepton number conservation, the dependence of the yield on $\mathcal {L}\rho _{\rm DDT}$ is stronger at lower fiducial $\mathcal {L}\rho _{{\rm DDT},0}$, as shown in previous sections. This effect is not captured in a simple proportionality relation like that quoted by Bravo et al. (2010), even if this effect is present in their calculations.

It has been discussed whether the source of scatter in peak brightness itself can be attributed to metallicity (Timmes et al. 2003; Howell et al. 2009; Bravo et al. 2010). We submit that multidimensional effects following from fluid instabilities during the deflagration phase leading to varying amounts of expansion provide scatter consistent with observations. The influence of metallicity on the DDT density affects the duration of the deflagration which will secondarily influence the magnitude of the scatter. Our results show that the primary parameter is the degree of (pre-)expansion before DDT—which determines the amount of mass at high density. This is, in turn, controlled by both the expansion rate and the plume rise time. These are also expected to be the basic ingredients in reality.

While our study is consistent with the expected theoretical brightness trend with metallicity (Timmes et al. 2003), observations to date have not been able to confirm this prediction (Gallagher et al. 2005, 2008; Howell et al. 2009). The progenitor age is difficult to decouple from metallicity given the mass–metallicity relationship within galaxy types (Gallagher et al. 2008). Additionally, the dependence of mean brightness of SNe Ia on the age of the parent stellar population appears to be much stronger than any dependence on metallicity (Gallagher et al. 2008; Howell et al. 2009). Seemingly, the only way to observe a dependence on metallicity is to constrain the mean stellar age by selecting galaxies of the same type. This approach was used by Gallagher et al. (2008), but the difficulties in accurately measuring metallicities for the parent stellar population have only constrained the magnitude of the effect on mean brightness, but thus far have not proved that a metallicity effect exists.

In order to determine the dependence of SN brightness (⁵⁶Ni yield) on metallicity, a better understanding of how the transition density is affected by changes in the ²²Ne content is needed. Currently, we have extrapolated data from Chamulak et al. (2007) down to the range of expected transition densities using flame data from densities above 10⁸ g cm⁻³. The trend observed by Chamulak et al. (2007) that increasing X₂₂ increases the laminar flame speed is valid for densities above 10⁸ g cm⁻³; however, it is unclear if this trend will hold for densities below 10⁷ g cm⁻³. Direct numerical simulations of the flame are needed for these lower densities to determine whether this trend still holds. However, in this parameter space, the laminar flame is extremely slow requiring a low-Mach-number treatment to model the flow. Regardless of the choice of model for the mechanism that produces a spontaneous detonation, a necessary condition is thought to be distributed burning. Properties of the laminar flame are necessary to estimate when burning becomes distributed. Future studies exploring other DDT mechanisms must also determine their compositional dependencies. So far, we have only explored the effect of varying the ²²Ne content directly (Townsley et al. 2009) and indirectly through the DDT density. Other effects exist that we have yet to study such as the carbon composition and flame ignition density (set by the average progenitor age). These properties may be influenced by metallicity such that the net effect on mean brightness is negligible. Additionally, if our assumptions about the dependence of DDT density on metallicity are incorrect and further studies indicate the opposite trend, the metallicity effect on DDT density could negate the effect due to lepton number conservation.

This study stresses the importance of multidimensional effects, but has provided evidence of the limitations of two-dimensional simulations of SNe. In the immediate future, we plan to perform three-dimensional simulations for better realism. In reality, we might expect a two-parameter family in which the morphology of the dominant RT plume is independent of the rate of expansion; however, for our two-dimensional simulations, these two effects appear to be correlated by the choice of cylindrical geometry. A plume developing along the symmetry axis is like three-dimensional RT compared to one near the equator, which is more like two-dimensional RT, and the latter is weaker (slower to develop). In addition to developing faster, plumes near the symmetry axis represent less volume and, therefore, expand the star less. A fast rising plume combined with a slow rate of expansion indicates there will be a shallow dependence on DDT density (since the plume will move through the various DDT densities faster and with less expansion). This result may also depend on our choice that DDT conditions are met at the tops of rising plumes. For these reasons, three-dimensional studies are necessary to ascertain whether there is a physically motivated correlation between the dominant plume morphology and the rate of expansion and whether this correlation depends on the choice of location of the DDT. In transitioning to three-dimensional simulations, we expect the growth of many more plumes; however, the rate of expansion will likely not exceed that found in two dimensions. Because RT develops faster in three dimensions, the conditions for DDT will be met sooner leading to less overall expansion of the star and more mass at high density at the first DDT time. Therefore, we expect three-dimensional results similar to the higher yielding two-dimensional realizations with a shallower dependence on DDT density. In any case, this study indicates that it is possible to determine the susceptibility of a particular model to a change of DDT density (and, hence, a change in the composition which alters the microphysics that set the DDT density). Once the plume morphology has been established several tenths of seconds into the simulated explosion, it should be possible to estimate the total ⁵⁶Ni yield to fairly high accuracy given DDT conditions. However, the treatment of turbulent flame properties may be important.

This work was supported by the Department of Energy through grants DE-FG02-07ER41516, DE-FG02-08ER41570, and DE-FG02-08ER41565, and by NASA through grant NNX09AD19G. A.C.C. also acknowledges support from the Department of Energy under grant DE-FG02-87ER40317. D.M.T. received support from the Bart J. Bok fellowship at the University of Arizona for part of this work. The authors acknowledge the hospitality of the Kavli Institute for Theoretical Physics, which is supported by the NSF under grant PHY05-51164, during the programs "Accretion and Explosion: the Astrophysics of Degenerate Stars" and "Stellar Death and Supernovae." The software used in this work was in part developed by the DOE-supported ASC/Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago. We thank Nathan Hearn for making his QuickFlash analysis tools publicly available at http://quickflash.sourceforge.net. We also thank the anonymous referee for a careful reading of the manuscript and constructive comments that improved this work. This work was supported in part by the US Department of Energy, Office of Nuclear Physics, under contract DE-AC02-06CH11357 and utilized resources at the New York Center for Computational Sciences at Stony Brook University/Brookhaven National Laboratory which is supported by the U.S. Department of Energy under Contract No. DE-AC02-98CH10886 and by the State of New York.