The Intrinsic Stochasticity of the ⁵⁶Ni Distribution of Single-degenerate Near-Chandrasekhar-mass SN Ia

Many classic works have identified white dwarfs (WDs) accreting to near the Chandrasekhar-mass M_ch in binary systems as candidate progenitors of SNe Ia—e.g., Whelan & Iben (1973), and Nomoto et al. (1984). This classic picture of near-Chandrasekhar-mass WDs in the single-degenerate channel was long thought to provide an explanation for the uniformity of brightnesses observed in SNe Ia (Phillips 1993). Although it is possible to achieve both sub-M_Ch and super-M_Ch WD SNe Ia in the SD channel (Maoz & Mannucci 2012), most work has focused upon this classic picture of near-M_Ch SD WDs, which will also be our focus in this paper.

The nature of the dominant production channel for SNe Ia has long been unclear (Branch et al. 1995). More recently, the classic picture of near M_ch WD models has been substantially revised, with single-degenerates now widely believed to be rare in nature. The single-degenerate channel has been shown to be inconsistent with a range of constraints, including the delay-time distribution, the absence of hydrogen, and the absence of companions (Maoz et al. 2014). Single-degenerates are also inconsistent with observational and theoretical rate predictions (Maoz & Mannucci 2012). However, recent observations have provided strong evidence that Chandrasekhar-mass WD SNe Ia do occur in at least some systems in nature. Hard X-ray spectra of the 3C 397 supernova remnant (SNR) are consistent with electron captures, which arise during the nuclear burning at high densities that is typical of Chandrasekhar-mass WDs (Yamaguchi et al. 2014, 2015). Additional X-ray and infrared observations of the Kepler SNR suggest that it was an overluminous single-degenerate supernova (Katsuda et al. 2015). Furthermore, the pre-maximum light shock signature that has been detected in both a subluminous SN Ia 2012cg (Marion et al. 2016) and a normal SN Ia iPTF14atg (Cao et al. 2015) is similar to theoretical predictions of the shock interaction with the companion star (Kasen 2010), although these observations have also been contested (Kromer et al. 2016; Shappee et al. 2018).

A large body of theoretical and computational work has explored possible mechanisms for single-degenerate SNe Ia (Maoz & Mannucci 2012). Many single-degenerate explosion mechanisms begin with a deflagration in the convective core of a near-M_Ch WD (Nomoto et al. 1984). From this common starting point, authors have explored the possibility of pure deflagrations (Röpke et al. 2007a; Jordan et al. 2012a; Kromer et al. 2013), deflagration-to-detonation transitions (DDTs) (Khokhlov 1991; Röpke et al. 2007b, Seitenzahl et al. 2013; Malone et al. 2014; Dave et al. 2017; Martínez-Rodríguez et al. 2017 and many more), and gravitationally confined detonations (GCDs) (Plewa et al. 2004; Röpke et al. 2007b; Townsley et al. 2007; Jordan et al. 2008; Meakin et al. 2009; Seitenzahl et al. 2016). The viability of the proposed explosion mechanisms hinges crucially on the nature of the flame ignition; i.e., the runaway of carbon fusion, during the convective phase. In particular, the GCD mechanism relies upon an offset ignition to buoyantly drive the flame bubble through breakout. Because the vigor of the GCD mechanism relies upon maintaining the WD intact until the ash collides at a point opposite of breakout, its viability is diminished as the ignitions become more centrally concentrated and multipoint. In contrast, a pure deflagration model produces good agreement with observations of the subclass of SNe Iax (Kromer et al. 2013), but requires a vigorous deflagration phase with several simultaneous near-central ignitions. Both the pure deflagration and the GCD mechanism require that the flame surface does not undergo a transition to a detonation prior to breakout, as the DDT model does. Furthermore, it is possible that a detonation does not arise during the initial ash collision subsequent to breakout and that the WD remains gravitationally bound, which leads to a subsequent contraction and a detonation through the pulsationally assisted GCD (PGCD) mechanism (García-Senz & Bravo 2005; Jordan et al. 2012a). Conditions for a classic GCD are harder to achieve in rigidly and differentially rotating near-Mch WDs (García-Senz et al. 2016), but a detonation remains possible during recontraction as a PGCD scenario.

Because the ignition of a flame bubble in the convective core of the WD is inherently stochastic, outcomes ranging from subluminous through overluminous SNe Ia are expected to arise in Chandrasekhar-mass SNe Ia. For example, large initial offsets used in simulations by Jordan et al. (2008) yield overluminous events, while a lopsided multi-bubble setup including ignitions at low offsets as studied by Jordan et al. (2012b) and Kromer et al. (2013) can lead to failed SNe Ia. It is thus important to properly incorporate the ignition distribution to determine likely outcomes in this study.

The ignition arises within a highly turbulent (Reynolds number Re ∼ 10¹⁵) convective flow (Isern et al. 2017), with the detailed outcome critically dependent upon the high-end tail of the temperature distribution. For many years, the distribution of ignition points was poorly constrained by theory and simulation (Garcia-Senz & Woosley 1995; Woosley et al. 2004). Early studies suggested multipoint ignitions as a viable scenario, which has only been revised recently as it became possible to begin to simulate these crucial last minutes of the simmering phase in full three-dimensional simulations. For example, Zingale et al. (2011) and Nonaka et al. (2012) performed a numerical study for a WD with a central density 2.2 × 10⁹ g cm⁻³ and central temperature 6.25 × 10⁸ K to determine the probability distribution of hot spots triggering the deflagration phase. Zingale et al. (2011) demonstrated that most ignitions for the WD model occur at a single point and at radial offsets below 100 km from the center, and most likely at about 50 km. Consequently, these ab initio simulations point toward a low amount of deflagration energy resulting from a small single bubble, buoyancy-driven ignition, which is in contrast to prior simulations that often invoked multiple-bubble ignitions. Thus, as ignitions arise from small single bubbles, we incorporate the implications from this in the discussion of self-consistency of the initial flame setup in the methodology, which has often been neglected in prior publications, where large spherically or perturbed bubbles where initialized inconsistently with the flame's evolution prior to the onset of the simulation. It has been known for some time that such low-deflagration energies generally lead to large amounts of ⁵⁶Ni. Röpke et al. (2007b) ran a series of simulations with off-centered ignitions demonstrating an anti-correlation of the ⁵⁶Ni deflagration yield and ignition offset. However, their initial offsets do not include ignitions below 50 km, thus neglecting roughly half of the ignitions expected from results in Nonaka et al. (2012). Hillebrandt et al. (2007) propose that off-centered, lopsided explosions, such as those following the deflagration phase simulated in Röpke et al. (2007b), might explain overluminous SN Ia events.

Recent theoretical work has explored the physics of stochastic ignition close to the center of the WD using semi-analytic methods and have demonstrated that single-bubble ignitions are generally buoyancy-dominated, leading to a weak deflagration phase (Fisher & Jumper 2015). Consequently, as Fisher & Jumper (2015) argued, single-bubble ignitions tend to lead to the production of a relatively large amount of ⁵⁶Ni and hence an overluminous SN Ia. This theoretical work was soon given observational support when spectral modeling of the nebular phase of SNe Ia revealed that the canonical bright event SN 1991T had an inferred ejecta mass of 1.4 M_⊙ (Childress et al. 2015). Most recently, Jiang et al. (2018) have examined the early-phase light curves of 40 SNe Ia in the optical, UV, and NUV, and demonstrated that all six luminous 91T- and 99aa-like events in their sample are associated with an early-excess that is consistent with a ⁵⁶Ni-abundant outer layer, as expected in the GCD scenario. Subsequent three-dimensional simulations of a buoyantly driven single-bubble ignition confirmed a large amount of ⁵⁶Ni consistent with SN 1991T (Seitenzahl et al. 2016). However, because the stable iron group elements (IGEs—specifically, elements from Ti–Cu) tend to be buoyantly driven in the GCD model, Seitenzahl et al. (2016) found that the observed stable IGEs at low velocities in their model could only be reproduced along a line of sight centered around the detonation region. While there are systemic differences in how the LEAFS code used by Seitenzahl et al. (2016) treats subgrid scale turbulent nuclear burning in comparison to FLASH (e.g., see Jordan et al. 2008), it is possible that the bulk of this inconsistency could be rectified by a DDT as opposed to a GCD model. In particular, Fisher & Jumper (2015) noted that buoyantly driven ignitions will lead to a large amount of ⁵⁶Ni and an overluminous SNe Ia in both the DDT and GCD models.

This recent observational and theoretical progress motivates the current study, in which we explore the inherent stochasticity of near-Chandrasekhar-mass WDs in the single-degenerate channel, from ignition through to detonation. In Section 2, we briefly summarize the simulation setup and the assumed initial hot spot distribution. In Section 3, we describe the WD evolution from ignition to its possible detonation depending on the ignition's offset to the center of mass and we link our findings to the initial hot spot distribution. In Section 4, we discuss possible uncertainties in our modeling before summarizing our findings in Section 5.

Our simulations were performed with the three-dimensional Eulerian adaptive mesh refinement (AMR) code FLASH 4.3 (Fryxell et al. 2000) and we solved the hydrodynamic equation with the directionally split piecewise-parabolic method (PPM). We use a tabular Helmholtz equation of state and take into account radiation, nuclei, electrons, positrons, and corrections for Coulomb effects, which remains valid in the electron-degenerate relativistic regime (Timmes & Swesty 2000). The flame physics are modeled by an advection-diffusion-reaction equation. Nuclear energy generation is incorporated using a simplified treatment of the flame energetics (Townsley et al. 2007, 2009, 2016). Self-gravity is accounted for by a multipole solver (Couch et al. 2013) up to order l = 6 with isolated boundary conditions.

The WD model used assumes a mass of 1.38 M_⊙ and a uniform 50/50 carbon/oxygen (C/O) composition. See Section 4 for a discussion of the impact of non-zero WD model metallicity. The WD has a central temperature stratification, including an adiabatic core with central density 2.2 × 10⁹ g cm⁻³ and temperature 7 × 10⁸ K, pressure-matched onto an isothermal envelope with temperature 10⁷ K (Jackson et al. 2010; Krueger et al. 2012). Furthermore, the central density of our WD model is a standard value commonly considered in the literature because higher central density WD models produce anomalously high abundances of neutron-rich isotopes, including ⁴⁸Ca, ⁵⁴Cr, and ⁶⁶Zn (Meyer et al. 1996; Nomoto et al. 1997; Woosley 1997; Brachwitz et al. 2000; Dave et al. 2017; Mori et al. 2018).

A very low density region surrounding the WD, which is sometimes referred to in the literature as "fluff," is required by Eulerian grid-based simulations, which cannot treat empty space without some matter density. The fluff is chosen to have an initial density of 10⁻³ g cm⁻³ and temperature of 3 × 10⁷ K, and is dynamically unimportant for the duration of the models presented here.

Since the deflagration energy release and the nucleosynthetic yield of ⁵⁶Ni hinges critically on the bubble initial conditions, we investigate the earliest phases of the bubble evolution in Section 3.1 in two dimensions, while performing full three-dimensional runs when evolving runs through deflagration and detonation. This choice was largely influenced by the fundamentally different behavior of the turbulent cascade in two dimensions and three dimensions. In particular, in two dimensions, the turbulent cascade is inverse, proceeding from smaller to larger scales (Kraichnan 1967). In contrast, in three dimensions, the turbulent cascade proceeds directly, from larger to smaller scales, where the energy is dissipated at the smallest scales due to viscosity. In the early-time simulations, the bubble remains laminar, and can be simulated in two dimensions. The fundamental distinctions between turbulence in two dimensions and three dimensions have major ramifications for the study of longer timescales, on which the flame becomes fully turbulent, since physically motivated flame-turbulence interaction subgrid models can only be realized in three dimensions. Consequently, all longer-time simulations in which the bubble enters a turbulent state have been run in three dimensions with Cartesian geometry with a turbulence-flame interaction model to capture the enhanced burning on subgrid scales. All three-dimensional simulations were performed both with and without the turbulence-flame interaction (TFI, also: 'enhanced burning') model (which will be described later), thereby spanning a range of possible outcomes on the flame propagation resulting from unresolved turbulence. Test two-dimensional simulations in cylindrical coordinates led to unphysical behavior in the turbulent phase, including spurious surface protuberances burning in the radial direction and thus significantly altering the simulation outcomes in comparison to the three-dimensional models. Artificial outcomes were particularly significant for runs with ignition points close to the center of mass of the WD, where unphysical burning in the radial direction has the largest impact.

Our Cartesian domain extends from −6.5536 × 10⁵ to +6.5536 × 10⁵ km in each direction, with a maximal refinement down to Δ = 4 km in three dimensions. We employ several refinement criteria, which are designed to follow the nuclear burning of the models at high resolution, while also minimizing the resolution in the very low density regions outside the WD itself. Our simulations seek to maintain the highest resolution in the burning region behind the flame surface and we employ a standard density gradient criterion to refine when the density gradient parameter exceeds 0.1 and derefine when it is beneath 0.0375. Further refinement criteria seek to derefine in the fluff and in regions outside of active burning, derefining one level if the energy generation rate is lower than 5 × 10¹⁷ erg g⁻¹ s⁻¹, and completely to level one if the density is below 10³ g cm⁻³. Except for their resolution and threshold, these criteria are the same as in Townsley et al. (2009). Furthermore, because the ejected ash continues to expand over time, the computational cost of following the ejected ash grows without bound. Consequently, we impose an additional derefinement outside a radius of 4000 km to Δ = 128 km, which only impacts the ejected mass. We increased this derefinement radius to 6000 km for offsets r₀ ≲ 20 km, where the pre-expansion can reach similar radii.

The initial size of the single flame bubble is limited by the hydrodynamic resolution of our simulations. At a resolution of Δ = 4 km for the flame front, we assume an initial spherical shape with a radius of R₀ = 16 km. For this to be a reasonable assumption, a self-consistent evolution since appearance of the hot spot should yield a negligible velocity profile and a self-similar evolution preserving sphericality, as discussed in Vladimirova (2007).

This consistency of assuming a spherical ignition point can be assessed with simple physical arguments. The flame polishing scale ${\lambda }_{\mathrm{fp}}=4\pi {S}_{l}^{2}/(\mathrm{Ag})$, below which perturbations on the surface are polished out (Timmes & Woosley 1992), implies that even if the initial ignition was non-spherical, then it would still become spherical soon afterwards. An explicit numerical test confirming this was performed by Malone et al. (2014). Later perturbations to the sphericality can arise from the background flow and the buoyant rise. We conceptually split the background flow into two contributions: possible large scale convective plumes and small-scale turbulence. The first, the convective plume velocity, while similar or exceeding the laminar flame velocity of ∼100 km s⁻¹ (Zingale et al. 2011; Piersanti et al. 2017), primarily induces a change in the reference frame rather than a change in the flame bubble's early dynamics. Its impact is shown in Section 4. The latter, the convective turbulence, is small in comparison to the laminar flame velocity (Nonaka et al. 2012), so that sphericality should initially be sustained. In addition to the initial background velocity field, the change of the velocity field due to the buoyant rise of the bubble itself prior to onset of the simulation needs to be considered for self-consistency. As a characteristic scale of a negligible built-up velocity for the bubble's evolution, we suggest stretching scale l_fl, defined as bubble size at which the buoyant velocity by the gravitational acceleration g reaches the order of the flame front's laminar velocity. This should approximately corresponds to ${l}_{\mathrm{fl}}=2{S}_{l}^{2}/(\mathrm{Ag})$ (Malone et al. 2014), which primarily depends on the offset near the WDs center. A is the Atwood number of fuel and ash density. Note that this criterion is stricter by a factor of 2π compared to the criterion for sphericality due to perturbations: the flame bubble is expected to stretch radially before the wrinkles in the flame front are no longer polished out. As an alternative to this estimate for l_fl, we integrate the flame's evolution based on Fisher & Jumper (2015) to determine when the bubble's velocity reaches the laminar flame speed, which is shown in Table 1. Only the small radial bubble offsets r₀ from the center, which we are particularly interested in, fulfill this condition. This length scales with ${r}_{0}^{-1}$ and therefore the condition allows large bubbles at low offsets. As Fisher & Jumper (2015) argue, there is a critical offset at which the deflagration will burn through the core, vastly changing the overall deflagration yield and thus its possible detonation. Therefore, a completely self-consistent evolution would be desirable at these offsets. However, we are effectively limited by the required computational resources and resolution. As an alternative for such a self-consistent treatment, we evolve two-dimensional models for the linear phase, where deviations to three-dimensional outcomes should be negligible, which we can resolve sufficiently well.

We incorporate a turbulence-flame interaction model presented in Jackson et al. (2014) to implement a specific model of power-law wrinkling based on that proposed by Charlette et al. (2002). The reaction front is modeled by a reaction-diffusion front, which propagates with a speed based upon the estimated physical features of the wrinkled physical flame whose width is, for most of the interior of the WD, many orders of magnitude smaller than the computational grid scale. Due to the interaction of turbulence with the flame, the location of the reaction front, as coarsened to a filter scale Δ consisting of a few grid cells, is approximated to propagate at a turbulent flame speed s_t = Ξs_l, where s_l is the physical laminar flame speed and Ξ is called the wrinkling factor. The wrinkling is given by

$$\begin{eqnarray}&&{\rm{\Xi }}={\left(1+\displaystyle \frac{{\rm{\Delta }}}{{\eta }_{c}}\right)}^{1/3},\end{eqnarray} \tag{ 1 }$$

where η_c is the cutoff scale for wrinkling, and is dependent upon local properties of both the turbulence and the physical flame. In this model, η_c is the inverse of the mean curvature of the flame surface, and is determined by assuming equilibrium between subgrid flame surface creation due to wrinkling by the turbulence and flame surface destruction by flame surface propagation and diffusion. This turbulence-flame interaction model leads to a turbulent flame speed s_t that is approximately equal to the characteristic speed of turbulent fluctuations on the filter scale, ${u}_{{\rm{\Delta }}}^{{\prime} }$ at intermediate densities, 10⁸–10⁹ g cm⁻³, as can be seen in Figure 4 in Jackson et al. (2014). The turbulent flame speed falls off to the laminar flame speed at lower densities, where the flame is too thick and slow to support wrinkling. At high densities where the flame is effectively polished by the high laminar speeds, the turbulent flame speed also approaches to the laminar flame speed value. Performing the calculation of the cutoff scale for wrinkling, η_c, requires a measurement of the turbulence on the filter scale, ${u}_{{\rm{\Delta }}}^{{\prime} }$, and makes the physical assumption that the subgrid turbulence is homogeneous, isotropic, and follows Kolmogorov's theory on the filter scale. As shown by Zingale et al. (2005), buoyancy-driven turbulence becomes increasingly homogeneous and isotropic on small scales, implying that the last assumption is valid provided that the filter scale is sufficiently small.

From Nonaka et al. (2012) we obtain the probability density function (PDF) P(r) of hot spots forming at a certain distance from the center of mass. Using the raw data and the same methodology, we create a histogram for such distribution and a fit to this distribution, see Figure 1. The PDF per unit length is shown. Under the assumption that the probability density per volume only mildly changes near the center of mass, this implies a P(r)dr ∝ r²dr scaling at low offsets due to the shrinking volume available for hot spots to occur. While not exactly fulfilling this consideration, we obtain a reasonable fit using a β-distribution, obtaining an expectation value of $\langle {r}_{0}\rangle =48$ km and a probability of 2.2% for hot spots forming at r₀ < 16 km, the critical ignition radius determined by Fisher & Jumper (2015).

Zoom In Zoom Out Reset image size

We simulate the outcomes of varying ignition offsets for the given WD model according to the this probability distribution and choose a representative range of initial offsets with an initial bubble radius of R₀ = 16 km shown in Table 2. We utilize both two-dimensional and three-dimensional simulations. The size of the initial bubble is naturally limited by the simulation resolution. As demonstrated analytically in Fisher & Jumper (2015), the flame bubble's dynamics vastly change at low initial offsets. We employ two-dimensional simulations to investigate the initial stages of the bubble dynamics at very high resolution. Moving to two dimensions is a reasonable strategy in this case because we are only interested in the initial, linear phase. In two dimensions, we employ a maximal resolution of 0.25 km and an initial bubble radius of 2 km, at initial offsets ranging from 4 to 50 km as listed in Table 1.

The three-dimensional simulations are evolved until they undergo a detonation as GCD. The precise conditions under which a DDT may arise are still a matter of active investigation, although recent three-dimensional simulations may shed further light on this issue (Poludnenko et al. 2011; Fisher et al. 2019). We adopt conservative criteria for detonation initiation based upon studies of the the Zel'dovich gradient mechanism (Seitenzahl et al. 2009), which demonstrate that the critical length above temperatures ≃2 × 10⁹ K at a density 10⁷ g cm⁻³ becomes of order 1 km, and a detonation is deemed likely once the ¹²C–¹²C reaction runs away. Furthermore, in this paper, we evolve all three-dimensional models within the context of the GCD scenario. As discussed in the introduction, current ab initio calculations point toward offset single-point ignitions, which favor both the GCD and DDT scenario over pure deflagrations. Because the GCD model involves further evolution post-bubble breakout, it generally predicts a greater deflagration energy release than the DDT model, for an otherwise identical WD model and flame bubble ignition model. Consequently, consideration of the GCD model yields a lower limit for the mass of ⁵⁶Ni due to a lower central density ρ_c at the time of detonation, in comparison to DDT models.

Artificial detonations can occur due to temperature oscillations arising as numerical artifacts from degenerate stellar equation of state coupled with hydrodynamics close to discontinuities (Zingale & Katz 2015). These oscillations are particularly prominent during the flame bubble's buoyant rise. To prevent detonations arising from these artifacts, we restrict detonations to occur in the southern hemisphere (z < 0 km).

We first discuss the early phase of bubble evolution. The early linear phase of evolution is similar in both two dimensions and three dimensions, so we investigate the early linear evolution in high-resolution two-dimensional models, and compare these against semi-analytic predictions. We then move on to examine the subsequent nonlinear evolution through breakout and detonation in full three dimensions, which we also evolve starting with the linear phase but at lower resolution.

3.1. Early Linear Evolution in Two Dimensions

Figure 2 shows the slices of the flame bubble's evolution in its laminar phase for r₀ = 10 km. The burned material grows spherically as long as the buoyant velocity is small and the bubble's size stays below the flame polishing scale. As the bubble grows, the acceleration for material at the northern and southern flame front start to differ and the bubble becomes elongated along the initial offset's direction until a plume forms at the northern front. Interestingly, the southern flame front's laminar speed seems to be countered by the background flow from the buoyant rise at the northern front.

Zoom In Zoom Out Reset image size

In Figure 3, we show the evolution of the flame bubble's radius R(r) as a function of the bubble offset r(t). For the simulation data, the volume-equivalent spherical radius ($R=\sqrt[3]{3V/4\pi }$) deduced from the burned volume V is shown. We compare our results from the initial phase of linear growth with the analytic model presented in Fisher & Jumper (2015) and find them to be in good agreement for the first tenths of seconds, particularly for larger initial offsets. However, the analytic description starts to fail as the velocity from the buoyant rise becomes inhomogeneous across the growing bubble, effectively stretching the bubble due to a lower speed at the southern flame front. Figure 4 shows the position of the southern flame front's position for the analytic and numerical evolution, which start to differ because the analytic solution does not incorporate an inhomogeneous velocity or acceleration field. The resulting elongation gives rise to a stem that is left behind the rising plume at the northern front. Even when the southern flame front crosses the center of mass, it will not buoyantly rise toward the opposite pole but is confined close to the center of mass on the relevant timescale due to the background flow caused by the buoyant rise of the ash on the northern hemisphere.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Nonlinear Evolution and Detonation

With the formation of a rising plume, the evolution becomes nonlinear and depends on the imposed flame model as presented in Section 2. To capture the flame's turbulent rise, we evolve three-dimensional models from ignition to detonation for the parameters listed in Table 2. Given that we find the three-dimensional runs to remain mostly symmetric, we show slices in the z–x plane restricted by x ≥ 0 km and y = 0 km. Figures 5 and 6 show the evolution of the WD for 20 and 100 km offset. Figure 7 also shows the evolution of a 20 km offset model, but without the enhanced burning that the prior two models use. Because the evolutionary timescales for each run depend on the initial conditions chosen, the slices for each run are chosen with respect to the state of the flame and not in absolute time. In particular, in each plot, the first frame shows the breakout of the flame at the star's surface. The next frame depicts the post-breakout flame crossing the z = 0 equator on the star's surface. The last frame shows the model just prior to detonation across from the point of breakout.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While low offsets also have a slightly larger distance to the surface of the WD, the evolution of shown models demonstrate the delay to breakout in comparison to larger offsets due to the smaller buoyant force near the center of mass, increasing the breakout time by roughly 0.3 s for r₀ = 20 km over r₀ = 100 km. Smaller initial offsets lead to a slightly increased plume size in both radial and tangential directions with respect to the center of mass across different initial offsets.

After breakout, the flame front travels around the WD and eventually reaches the point opposite breakout. The envelope material is pushed in front of this flame front into this opposing point. At larger offsets than 40 km, such as the shown 100 km run, the ram pressure building up suffices to trigger a detonation before the ash reaches the opposing point initiated by the ¹²C–¹²C reaction running away. For offsets smaller than roughly 40 km, such as the runs with an initial offset of 20 km, the ram pressure is insufficient to trigger a detonation upon the flame reaching the opposite pole and a detonation only occurs after a subsequent partial recontraction of the WD delaying the detonation. For some offsets lower than the shown 20 km, the WD might not detonate upon recontraction.

The difference of the enhanced burning model seems to be only moderate for the shown slices at offset r₀ = 20 km. The evolutionary phases represented by the slices coincide between flame models. Meanwhile, without enhanced burning, a larger fraction of the material ejected from the WD seems to be burned.

Figure 8 shows the estimated ⁵⁶Ni yield over time relative to the time of ignition. The ⁵⁶Ni yield in each model is obtained from the electron mass fraction Y_e, and by assuming the neutronization of IGE occurs in equal parts by mass of ⁵⁴Fe and ⁵⁸Ni for all Y_e, which holds within 2% in tabulated yields from previous models (Meakin et al. 2009; Townsley et al. 2009). The rapid increase of M_Ni56 occurring at t ≳ 2.0 s indicates the onset of detonation, except for offset r₀ = 0 km, where a large growth sets in at t = 1.0 s due to turbulent deflagration.

Zoom In Zoom Out Reset image size

The vast majority of ignitions producing large amounts of ⁵⁶Ni, as here, produce ∼0.1 M_⊙ of IMEs (Meakin et al. 2009). This relatively low yield of IMEs is well-known to be inconsistent with normal Ia. A full comparison against observation is non-trivial because the spectral signature of the observations depends not solely upon the IME yield but also on its distribution. Given the high ⁵⁶Ni yields and low IME yields, the resulting model synthetic light curves and spectra are, however, unlikely to correspond to normal SNe Ia.

We classify the WD model's evolution into three different classes: failed, GCD and PGCD. The PGCD scenario introduced by Jordan et al. (2012a) shows a strong recontraction phase ($2\cdot {t}_{\det ,\mathrm{GCD}}\lesssim {t}_{\det ,\mathrm{PGCD}}$) due a significantly increased deflagration yield from a many-bubble-ignition setup. Given the smaller deflagration yields in our simulations due to the single ignition hot spot, PGCDs here only show a mild recontraction (as e.g., indicated by the evolution of the central density), making the transition between PGCDs and GCDs gradual. Sometimes, there is no detonation upon buildup of ram pressure but only when the fuel-ash mixture reaches the southern pole, even if no clear recontraction is present. Therefore, we do not impose a binary criterion between those scenarios.

As Figure 9 and Table 2 show, the ⁵⁶Ni yield seems to converge toward overluminous events with roughly 1.21 M_⊙ at large offsets r ⪆ 100 km with the WD models undergoing the GCD scenario. As suspected by Fisher & Jumper (2015), the ⁵⁶Ni yield decreases with lower initial offset as visualized in Figure 9. For runs ≲40 km a transition toward PGCD-like scenarios takes place. The transition is not monotonic but has a stochastic component, whether or not the ram pressure from the initial deflagration will suffice for an imminent detonation. For example, in our simulations a 32 km offset suffices for a GCD, despite the onset of a PGCD already at r₀ = 40 km.

Zoom In Zoom Out Reset image size

At some point we might expect the PGCD scenario to fail. However, we lack the spatial resolution at very small radial offsets to determine the location of this transition. Thus, we are left with the artificial case of central ignition for which the ⁵⁶Ni yield from deflagration increases to 0.56 M_⊙/0.35 M_⊙ (TFI/no TFI), depending on the flame model. However, there is no detonation, so that the total yield drops to these values as compared to higher initial offsets. Due to the numerical expense of resolving the energy generating regions at maximal resolution with AMR, we had to derefine to a maximal resolution of Δ = 8 km at t = 1.82/2.00 s for the TFI/no TFI scenarios with r₀ = 0 km. Based on our parameter sampling, the transition from PGCD to failed must occur at 0 km < r₀ < 16 km for the chosen WD model.

At large offsets, the yields with and without enhanced burning model vary only marginally. However, at lower offsets the enhanced burning significantly adds to the ⁵⁶Ni yield, not only for the failed detonations where differences add up to 60% but also for PGCD events in the order of up to 10%.

We note that in the context of the current paper, a "failed" model is simply one that does not detonate—e.g., it is a pure deflagration. Prior work in the literature has made a strong connection between such pure deflagration models and the observed class of subluminous Iax events (Kromer et al. 2013).

3.3. Likelihood of ⁵⁶Ni Yields

We next compute the probability distribution of M_Ni56 outcomes for the presented GCD SD channel models. The transformation from the hot spot probability distribution P(r₀) to the probability distribution $P({M}_{\mathrm{Ni}56})$ of ⁵⁶Ni outcomes is given as

$$\begin{eqnarray}&&P({M}_{\mathrm{Ni}56})=\displaystyle \sum _{{r}_{0}\in {g}^{-1}({M}_{\mathrm{Ni}56})}\displaystyle \frac{P({r}_{0})}{\left|g^{\prime} ({r}_{0})\right|}.\end{eqnarray} \tag{ 2 }$$

Here g (r₀) ≡ M_Ni56(r₀) is the amount of ⁵⁶Ni produced as a function of offset radius r₀. P(r₀) is the hot spot distribution found in Zingale et al. (2011) shown in Figure 1. This relationship may be derived from Bayes' Theorem with minimal assumptions. We start with

$$\begin{eqnarray}&&P({M}_{\mathrm{Ni}56}| {r}_{0})P({r}_{0})=P({r}_{0}| {M}_{\mathrm{Ni}56})P({M}_{\mathrm{Ni}56}),\end{eqnarray} \tag{ 3 }$$

and assume a simplification that the ⁵⁶Ni yield is solely determined by the offset position, leaving out possible uncertainties from velocity flow and early bifurcations arising in the turbulent phase, one finds

$$\begin{eqnarray*}\begin{array}{rcl}P({M}_{\mathrm{Ni}56}| {r}_{0}) & = & \delta ({g}_{{M}_{\mathrm{Ni}56}}({r}_{0})-{M}_{\mathrm{Ni}56}),\\ P({r}_{0}| {M}_{\mathrm{Ni}56}) & = & \delta ({r}_{0}-{g}^{-1}({M}_{\mathrm{Ni}56})).\end{array}\end{eqnarray*}$$

Here δ(x) is the Dirac delta distribution. Finally, using the identity

$$\begin{eqnarray}&&\delta (f(x))=\displaystyle \sum _{i}\displaystyle \frac{\delta (x-{x}_{i})}{\left|{f}^{{\prime} }({x}_{i})\right|},\end{eqnarray} \tag{ 4 }$$

where we sum over the roots x_i of f(x), we can rewrite Equation (3) as (2).

Even for a very similar stellar structure, we expect bifurcations arising from the turbulent nature of the flame bubble's buoyant rise affecting the final ⁵⁶Ni yield. In our computations, we see a similar phenomenon from slight offset changes and perturbations of the initial flame bubble. Therefore, the yields that we have obtained do not need to follow a monotonous relationship because only one possible realization is drawn at a given offset. Given the numerical expense, we can neither evaluate multiple runs at the same offset with slightly modified WD structure or flame bubble nor can we afford to run more models at different offsets. However, if deemed significant, additional parameters such as varied background velocity field, could easily be incorporated by marginalizing over such parameters. With our limited data sample, we nevertheless try to obtain insights into the resulting spread of ⁵⁶Ni yields given the stochastic nature of the initial ignition offset. To do so, we impose a strictly monotonous fit function for M_Ni56(r₀) with an asymptotic yield at high offsets r₀ for which we use

$$\begin{eqnarray}&&y({r}_{0})=\displaystyle \frac{{y}_{\max }+{\rm{\Delta }}y}{2}+\displaystyle \frac{{y}_{\max }-{\rm{\Delta }}y}{2}\cdot \tanh \left(s\cdot ({r}_{0}-{r}_{s})\right),\end{eqnarray} \tag{ 5 }$$

where y_max is the asymptotic yield at high offsets, Δy the spread between the two asymptotic branches, r_s the position of the turning point and s characterizes the sensitivity of the ⁵⁶Ni yield with respect to the initial offset r₀. We fix the asymptotic yield to the approximate value y_max found earlier.

The resulting distributions P(M_Ni56) for our WD model with and without enhanced burning are shown in Figure 10. We show the probability distributions for events with offsets larger than 16 km, accounting for 97.8% of the ignitions. The distribution shows a slightly larger spread for the non-TFI models due to the lower ⁵⁶Ni yield at low ignition offsets. Nevertheless, we find that the majority of ignitions result in a very confined ⁵⁶Ni yield.

Zoom In Zoom Out Reset image size

For a given hot spot distribution P(r₀) and the outcomes M_Ni56(r₀) of the simulated WD model, we get a stochastic spread P(M_Ni56) that strongly favors overluminous events with a ⁵⁶Ni yield of ∼1.2 M_⊙ with a standard deviation of σ ∼ 0.03. However, the total spread in outcomes of detonating models $\delta =\max ({M}_{\mathrm{Ni}56}({r}_{0}))-\min ({M}_{\mathrm{Ni}56}({r}_{0}))$ due to the stochasticity of hot spots forming is significantly larger. We are limited by the hydrodynamic resolution, but find that δ ≳ 0.2 M_⊙.

We find σ ≪ δ for both TFI and no TFI models as the ⁵⁶Ni yield is already close to asymptotic value of 1.21 M_⊙ at radii r₀ ∼ 50 km, which is the most likely point of ignition for the assumed hot spot distribution. For a hot spot distribution peaking closer to the WD's center of mass, the standard deviation could be significantly higher.

The exact shape of the left tail of the distribution is highly uncertain because it depends on the chosen fit function given our sparse sampling. Similarly, we expect modeling uncertainties, such as due to the lack of a velocity field, to propagate most severely into the ⁵⁶Ni yield and the resulting probability distribution at low offset radii as the buoyant evolution can be strongly enhanced or delayed. Other stochastic parameters, such as the state of the WD's velocity field, were not considered but would have to be marginalized to obtain the probability distribution in a more elaborate study.

Given that we only study a single ignition bubble setup, the result might be altered by taking into account the possibility of a multi-ignition setup. First, the analysis of the hot spot distribution by Nonaka et al. (2012), which we use here, disfavors multiple ignition points. Thus, as a matter of self-consistency, we do not consider these. Second, even if considering the option of a multi-ignition setup, given the limited analysis that has been possible for the data set provided in Nonaka et al. (2012) on this matter, we expect the impact of additional ignition bubbles to be limited because the background velocity flow induced by the buoyant rise of the first bubble quickly carries secondary bubbles out of the center of the WD into regions where the stem of the original bubble resides. The biggest impact of multiple ignition points is expected when this background velocity field is not yet established (${\rm{\Delta }}t\lesssim {l}_{\mathrm{fl}}$), so that the deflagration yield might increase substantially.

We have shown how the ignition offset probability distribution directly links to a range of SNe Ia outcomes, parameterized by the ⁵⁶Ni yield. This range of SNe Ia outcomes is intrinsically connected to the turbulent convective velocity field in the near-M_Ch WD, which causes the ignition of the SD channel to be inherently stochastic and unpredictable. The physics of the SD channel is complex and is subject to numerous modeling uncertainties: the pre-WD stellar evolution, accretion from the companion, possibly impacting the WD initial composition and structure, the physics of the simmering phase leading up to ignition, and the physics of turbulent nuclear burning and detonation. In the following, we discuss this range of modeling uncertainties and to what extent each of these effects may impact our conclusions.

In this work, we have incorporated the statistical distribution of ignition points drawn from actual three-dimensional simulations of the convective simmering phase of near-M_Ch WDs leading up to ignition. While this approach has clear advantages over the majority of prior work, which typically adopted ignition points in an arbitrary fashion, it is nonetheless still limited by the fact that there has only been one high-quality three-dimensional simmering phase simulation in a single WD model completed to date. The simulation has been performed at increasing resolution and the distribution of hot spot offsets appears to be converged (Zingale et al. 2011; Nonaka et al. 2012; Malone et al. 2014). However, it is conceivable that the distribution of hot spots could be more centrally condensed in WD models with higher central density. Additionally, the physical and chemical structure of the simulated WD was assumed ad hoc, similar to the WD used here, and was not a self-consistent result that traces the whole simmering phase of 10³–10⁴ yr (Martínez-Rodríguez et al. 2017; Piersanti et al. 2017), which adds to the shortcomings of the hot spot study and the uncertainties of the hot spot probability distribution used here.

The mechanism underlying the initiation of detonation plays an important role in SNe Ia theory and much effort has focused upon whether the detonation mechanism in a near-M_Ch SD scenario is a DDT or GCD (Röpke et al. 2007b; Seitenzahl et al. 2016; Dave et al. 2017). For example, because the DDT detonates prior to bubble breakout, the stratification of the ⁵⁶Ni and IGEs is generally more centrally condensed, which is in broader agreement with observations of high M_Ni56, overluminous SNe Ia like 91T (Seitenzahl et al. 2016). In the current work, we have focused upon the GCD mechanism in our simulations when inferring the intrinsic variation of the ⁵⁶Ni production resulting from stochastic ignition. However, if we were to have instead adopted a DDT criterion for detonation initiation, the ⁵⁶Ni distribution would be even more heavily left-skewed because, given identically the same WD model and ignition point, the DDT detonates prior to breakout, and consequently always results in a less pre-expanded WD model than a GCD (Dave et al. 2017). Therefore, the conclusion that the stochastic variance in ⁵⁶Ni yields is small, and the mean ⁵⁶Ni yield is large, is not expected to be qualitatively modified under the DDT scenario for our simulation setup.

Another crucial piece of the physics underlying both the simmering phase and the nuclear burning within the SNe Ia is the rate for C12 + C12 fusion. Recent experiments have measured this reaction rate for the first time for center-of-mass energies in the range of 0.8–2.5 MeV and have demonstrated an enhancement in the cross sections by as much as a factor of 25 in a key temperature range relevant to SNe Ia using the Trojan horse method (Tumino et al. 2018a). While this work has been contested by other authors (Mukhamedzhanov & Pang 2018; and subsequently rebutted by Tumino et al. 2018b) and will ultimately await additional confirmation, it is important to recognize the possible impact that the uncertainties in this key reaction rate may have upon the SNe Ia modeling. A higher reaction rate would increase the flame speed and might particularly change the M_Ni56 outcomes at low offsets where the buoyant evolution is most sensitive to changes of our fiducial model.

Furthermore, we began with a quiescent WD, although the ignition arises in the WD interior, which is itself convective and, as a consequence of the transport of angular momentum from the accretion stream, may itself be rotating. Indeed, recent work has shown that the effect of rotation may be significant enough to weaken the convergence in the detonation region of a classical GCD, although a detonation can still occur during recontraction (García-Senz et al. 2016) as a PGCD (Jordan et al. 2012a). Furthermore, at low ignition offsets, the magnitude of the initial convective velocity field may have an impact on the early flame bubble's evolution, and thus the ⁵⁶Ni yield. Because we start our simulations with zero velocity, this adds an additional uncertainty in the resulting M_Ni56 distribution. On the one hand, there is small-scale turbulence, which distorts the flame front early on. The expected velocities for this are small (∼10 km s⁻¹) with regard to the laminar flame speed (∼100 km s⁻¹), so that the flame bubble's sphericality is still mostly unaffected until broken by its buoyant rise. Therefore, minor shifts in a possible failed-to-detonated transition radius might be expected. On the other hand, there is a possibility of the ignition point occurring in a larger convective flow, similar or exceeding the laminar flame speed (Zingale et al. 2011; Nonaka et al. 2012). If these hot spots form in convectively outward moving regions as found by Nonaka et al. (2012), then this will further decrease the probability for ignitions that burn through the WD's center. Malone et al. (2014) ran a series of numerical simulations similar to our setup for the deflagration phase but also included a comparison of a setup with self-consistent convective velocity field and one without any velocity field. In these simulations, the authors find that the influence of the initial flow field increases as the initial ignition point is set closer to the center of mass, particularly for an exactly centered ignition.

Stellar composition influences the final nucleosynthetic yield of a SD SNe Ia through a variety of effects. The CNO metals of the WD stellar progenitor ultimately yield ²²Ne during He burning. Umeda et al. (1999) suggested that a variation in the carbon abundance within the WD model in the single-degenerate channel would impact the production of ⁵⁶Ni. In particular, Umeda et al. (1999) conjectured that WDs with a richer C/O ratio would lead to a more turbulent flame, an earlier transition from deflagration to detonation at higher densities, and hence a greater production of ⁵⁶Ni. The neutron excess carried by ²²Ne results in a decrease in the M_Ni56 of the SN Ia event (Timmes et al. 2003; Bravo et al. 2010), which linearly declines with the abundance of ²²Ne. Chamulak et al. (2007) studied the speedup of the laminar flame with increasing ²²Ne, which was further expanded by Townsley et al. (2009), considering additional compositional effects influencing the final nucleosynthetic yields, including the ignition density, the energy release, the WD structure, and the density at which a possible deflagration-to-detonation transition arises. The simulations with ²²Ne mass fractions increasing from 0 to 0.02, which were run long enough to determine a final ⁵⁶Ni yield, demonstrate that the combination of these effects result in a roughly 10% decrease in M_Ni56. Similarly, we expect a slight decrease in M_Ni56 based on complementary work by Jackson et al. (2010) investigating the impact of the ²²Ne content on the DDT density and the resulting ⁵⁶Ni mass. Furthermore, computational simulations of single-degenerate SNe Ia have subsequently explored the influence of varying the C/O ratio within the WD model in the context of the DDT model (Höflich et al. 1998; Domínguez et al. 2001; Krueger et al. 2010; Ohlmann et al. 2014). These investigations have demonstrated that higher C/O ratios yield more energetic and more luminous SNe Ia.

Taken together, this body of work on SD SNe Ia generally supports the view that stellar progenitor C/O ratio and metallicity play a role in determining the brightness of a SN Ia event. However, at the same time, these models have demonstrated that additional free parameters, including both the number and distribution of ignition points, as well as the DDT transition density, have a combined effect on the explosion energy comparable to that of the C/O ratio and stellar progenitor metallicity. Moreover, based upon this body of work, the combined influence of both a decrease in the C/O ratio and an increase in the stellar progenitor metallicity from the values assumed here (50/50 and 0, respectively) would result in a 10%–20% decrease in the ⁵⁶Ni yields, which would quantitatively impact our predicted M_Ni56 distribution but would not alone yield a distribution more closely resembling normal SNe Ia.

Most simulation models of near-M_Ch WDs adopt a central density ρ_c ≃ 2 × 10⁹ g cm⁻³, as we have in this paper. Because the electron capture rates are highly sensitive to the density, higher central density WDs generally produce greater amounts of stable IGE, and a lower ⁵⁶Ni yield. Higher central density WDs significantly overproduce (relative to solar) a range of neutron-rich isotopes, including ⁵⁰Ti, ⁵⁴Cr, ⁵⁸Fe, and ⁶²Ni, and were consequently generally excluded from consideration as near-M_Ch WD models (Meyer et al. 1996; Nomoto et al. 1997; Woosley 1997; Brachwitz et al. 2000). However, if SD near-M_Ch WDs constitute a small fraction of all SNe Ia, then such high-central density WDs may not be rare occurrences. If the central density of the near-M_Ch WD is indeed higher than ρ_c ≃ 2 × 10⁹ g cm⁻³, then the flame speed and the consequent deflagration energy release can be greater than considered here. This can in turn lead to greater pre-expansion and a reduced amount of ⁵⁶Ni as shown in two-dimensional simulations (Krueger et al. 2012; Dave et al. 2017), possibly consistent with a normal or even a failed SNe Ia. However, the qualitative outcome of an increased higher central density can vary, as shown by Seitenzahl et al. (2011). In their three-dimensional numerical study of the DDT scenario, the authors of the latter study show that the central density is only a secondary parameter. However, their study assumed multipoint ignition over a wide range of ignition kernels. When single-point ignitions are adopted, increased electron capture rates at higher central densities lead to higher abundances of neutron-rich IGEs at the expense of ⁵⁶Ni (Dave et al. 2017).

In this paper, for the first time, a small single ignition point as predicted by ab initio models of the convective phase leading up to the ignition is placed in a fiducial 50/50 C/O WD with a central density of 2.2 × 10⁹ g cm⁻³ and an adiabatic temperature profile leading up to SNe Ia in the GCD scenario.

We showed that a transition to failed SNe Ia (i.e., those events lacking a GCD) occurs as r₀ falls below some offset below 16 km. Even for those WDs detonating, the ⁵⁶Ni yield spawns a range of outcomes changing by 10%–20% with a decreasing yield as r₀ approaches the radius where no detonation is triggered.

We summarize our key conclusions as follows:

• 1.  
  Stochastic range of outcomes. For the chosen WD model, this corresponds to a spread of δ ≳ 0.2 M_⊙ for detonating models, even though the M_Ni56 distribution is strongly left-skewed so that low M_Ni56 are unlikely for the given probability distribution. This range of outcomes is stochastic and will add to other variations from the different WD models' stellar structure and evolution.
• 2.  
  For non-centered ignitions, all ignitions lead up to an overluminous SNe Ia. While some prior simulations also demonstrate low pre-expansion as a result of off-centered single ignition setups (Townsley et al. 2007; Jordan et al. 2008) and possibly an overluminous SNe Ia for a successful GCD scenario, we confirm this conclusion for a self-consistent modeling of the flame bubble's evolution at even very small offsets. We do not find a viable scenario from a single-bubble ignition leading to a normal Type Ia for the WD model adopted here, which is also commonly referenced in literature. This disfavors near-Chandrasekhar-mass WD explosions as a contributing channel to failed and normal type Ia SNe. If this channel was to contribute to failed and normal type Ia supernovae, then this would require readjustment and better understanding of the stellar structure and evolution, and also the flame dynamics.
• 3.  
  (Quasi-)symmetric deflagrations around the center of mass, as commonly used in numerical studies, are most likely artificial constructs: ignitions very close to the center are rare, as shown by Nonaka et al. (2012), and even if such events occur, then a strong asymmetry with a one sided plume is most likely to arise. We demonstrated that this happens as a background flow induced by the buoyant rise of the outermost flame front counteracts the burning into other directions. However, future in-depth studies on the likelihood of multi-ignition occurrences and their correlations with the turbulent velocity field might leave room for rare occurrences of symmetric deflagrations.

The authors thank Mike Zingale for generously providing the simulation data from his group's convective simmering phase models, which was used in Figure 1. The authors thank Pranav Dave and Rahul Kashyap for insightful conversations. R.T.F. also thanks the Institute for Theory and Computation at the Harvard-Smithsonian Center for Astrophysics for visiting support, during which a portion of this work was undertaken. C.B. acknowledges support from the Deutscher Akademischer Austauschdienst (DAAD). R.T.F. acknowledges support from NASA ATP award 80NSSC18K1013. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) Stampede 2 supercomputer at the University of Texas at Austin's Texas Advanced Computing Center through allocation TG-AST100038. XSEDE is supported by National Science Foundation grant number ACI-1548562 (Towns et al. 2014).

We use a modified version of the FLASH code 4.3 (Fryxell et al. 2000), which was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago, for our simulations including the SN Ia modeling presented in Townsley et al. (2016). The authors have made use of Frank Timmes' hot white dwarf model (http://cococubed.asu.edu/code_pages/adiabatic_white_dwarf.shtml). For our analysis, we acknowledge use of the Python programming language (Van Rossum & de Boer 1991), the use of the Numpy (van der Walt et al. 2011), IPython (Perez & Granger 2007), and Matplotlib (Hunter 2007) packages. Our analysis and plots strongly benefited from the use of the yt package (Turk et al. 2011).