Exploration of Aspherical Ejecta Properties in Type Ia Supernovae: Progenitor Dependence and Applications to Progenitor Classification

Type Ia supernovae (SNe Ia) have demonstrated a wide diversity (e.g., Taubenberger 2017). To trace the progenitor and the explosion mechanisms, supernova (SN) light curves and spectra are indispensable. (See, e.g., Hillebrand & Niemeyer 2000; Maoz et al. 2014; Nomoto & Leung 2018; and the Appendix for a summary of the progenitors and the explosion mechanism.) Explosion features can be extracted from the light-curve shape (e.g., Blondin et al. 2018), spectral lines such as Ca ii and Ni ii (e.g., Wilk et al. 2018), and the polarization (e.g., Bulla et al. 2016). The element abundance ratio can be the discriminant of current explosion models (e.g., Seitenzahl et al. 2009) by indicating the isotopic abundance ratios (e.g., ⁵⁷Ni/⁵⁶Ni, ⁵⁵Fe/⁵⁶Ni, and ⁴⁴Ti/⁵⁶Fe) and element abundance ratios (e.g., Mn/Fe and Ni/Fe; see Mori et al. 2018). These techniques had been applied to a few well-observed SNe Ia, including SNe 2011fe, 2012cg, 2014J, and 2015F and SN remnant (SNR) 3C 397 (see recent work in, for example, Yamaguchi et al. 2015; Dave et al. 2017; Leung & Nomoto 2017a, 2018; Shen et al. 2018; Zhou et al. 2021). The progenitor mass and metallicity are constrained by this method, which allows further identification of isotopic mass fractions of radioactive Ni isotopes in SN 2012cg (Graur et al. 2016) and SN 2014J (Yang et al. 2018) and stable Ni by the nebular IR spectra (see Dhawan et al. 2018, for the application to SN 2014J).

Radioactive isotopes with a longer half-life, e.g., ⁴⁴Ti with a half-year of 59 yr, can help to distinguish the different subclasses of SNe and their explosion mechanisms by their effects on the decline rate of the late-time light curve and measuring the absolute abundances (see, e.g., Fry et al. 2015).

Based on the tomography and morphology of SNRs, the progenitor constraints are cast on the SNR Tycho (Wang & Li 2014; Lopez et al. 2015). The S and Fe spectral lines obtained from SNRs 0519-69.0, 0509-67.5, and N103B are further examined in the shock-heated nebulae (Seitenzahl et al. 2019). The heated matter can emit X-ray lines that provide direct constraints on the explosion mechanism, in particular the asphericity of SNe Ia.

Analysis of SNRs on larger scales, including the Small and Large Magellanic Clouds, can reveal the SN explosion history (Maggi et al. 2016, 2019). The study of particular elements, e.g., Mn, can point out the major explosion mechanism (McWilliam et al. 2018; de los Reyes et al. 2020; Kobayashi et al. 2020).

The morphology reveals the inherent explosion asymmetry. For individual SNRs, the possible explosion progenitors can be inferred from their X-ray emission, which reveals the chemical abundances (Badenes et al. 2005, 2006, 2008). Differences in the interaction between the SN ejecta and the ambient medium can also point out the progenitor (Martínez-Rodríguez et al. 2018). The history of ejecta interaction with circumstellar matter (CSM) is contained in the late-time light curve as well. ⁶

We have performed nucleosynthesis surveys of SN Ia explosions by the near-Chandrasekhar-mass white dwarf (WD) model using the turbulent deflagration model with the deflagration–detonation transition in Leung & Nomoto (2018, hereafter Paper I), the sub-Chandrasekhar-mass WD model using the double-detonation model in Leung & Nomoto (2020a, hereafter Paper II), and the near-Chandrasekhar-mass model for SNe Iax using the pure turbulent deflagration model in Leung & Nomoto (2020b, hereafter Paper III). In these papers, we conducted an extended parameter survey of the SN Ia models aimed at understanding the implication of its observed diversity to progenitor model diversity, as well as constraints on explosion physics.

In Paper I, we considered models with the masses of M = 1.30–1.38 M_⊙ and metallicity of Z = 0–5 Z_⊙. The supersolar metallicity is suggested to explain the observed SNe Ia, including the SNR 3C 397 ($Z\sim$ 5 Z_⊙) and SN 2012cg ($Z\sim$ 3–5 Z_⊙). In Paper II, we considered WD models with M = 0.9–1.3 M_⊙, Z = 0–5 Z_⊙, and He-envelope mass of M_He = 0.05–0.2 M_⊙. We show that sub-Chandrasekhar-mass models can also explain the isotopic ratio of nearby SNe Ia, as the near-Chandrasekhar-mass model does, but its Mn production cannot explain the nearby [Mn/Fe] trend taken from stars in the solar neighborhood. Explosions of WD models with $M\sim$ 1.2 M_⊙ can provide the key to distinguish the two explosion channels. In Paper III, we specifically model SNe Iax with the pure turbulent deflagration model of the WDs with M = 1.30–1.39 M_⊙ and $Z={Z}_{\odot }$ (Leung & Nomoto 2020b).

In the present Paper IV, we study the asymmetry of the WD models as suggested from SN 2014J, the SNR Tycho, and other features mentioned above. The ejecta geometry is primarily dependent on the explosion progenitor. The near-Chandrasekhar-mass model exploding by the turbulent deflagration model with the deflagration–detonation transition tends to explode spherically (see, e.g., Röpke et al. 2007, for recent three-dimensional realizations). On the contrary, the sub-Chandrasekhar-mass model that explodes by the He-induced double detonations can generate large-scale asymmetry because of the off-center trigger of the explosion (see, e.g., Moll & Woosley 2013; García-Senz et al. 2018; Gronow et al. 2020, for some recent three-dimensional realizations showing aspherical structures).

We thus focus on the sub-Chandrasekhar-mass models, which tend to have a more aspherical structure than the near-Chandrasekhar-mass models. We note that in the literature, there is no extensive work examining how sub-Chandrasekhar-mass models exhibit a large-scale asphericity. In this Paper IV, we study for the first time how the different detonation mechanisms of the sub-Chandrasekhar-mass model can generate the ejecta deviated from the canonical spherical model by multidimensional simulations.

Specifically, we will clarify how the detonation triggered in the He envelope affects the ⁵⁶Ni distribution and ejecta structure in both position and velocity spaces. We also try to understand the underlying principles of the observed irregularities in SNRs.

We choose two-dimensional models so that we may test a larger number of models than three-dimensional models to uncover the global trend of the parameter landscape. Also, quasi-spherical SN 2014J (see Section 4) has encouraged us to explore models with a certain level of symmetry (e.g., rotation and reflection symmetry) assumed in two-dimensional models, instead of arbitrary models without explicit symmetry, as in three-dimensional models.

Ideally, three-dimensional models are naturally desired to match all features self-consistently. However, they are much more computationally expensive. As an exploratory study, we aim at searching the key properties in the detonation setting for the models to contain different features that might ultimately imprint in observational data. We want to understand what kind of shock interaction is necessary for generating the observed features, from which we may obtain hints of the initial detonation pattern. This will guide future three-dimensional simulations in setting up accurate initial models aiming to explain SN 2014J or other SNe.

In Section 2, we describe our methodology and the models to be presented in this article. Then, we briefly review how we compute the explosion models of the sub-Chandrasekhar-mass WD as the progenitor, and we present the stellar parameters, explosion energetics, and essential nucleosynthetic products.

In Section 3, we first present how the near- and sub-Chandrasekhar-mass WDs differ by their large-scale asymmetry. Then, we examine in detail how the explosion ejecta and its chemical composition depend on the viewing (ejecta) angle. We also predict the expected morphology by extracting the representative elements.

In Section 4, we present a detailed case study for SN 2014J to show how the asymmetry and other nucleosynthetic yields can be used for constraining the explosion mechanism and the progenitors of SNe Ia. We cast constraints on the progenitor mass, initial explosion geometry, He-envelope mass, and metallicity.

In Section 5, we further apply our results to some recently observed SNRs reported in the literature to demonstrate how the geometry can provide us the hints.

In Section 6, we discuss how this work is related to other hydrodynamics simulations in the literature.

We follow Leung & Nomoto (2020a) for SNe Ia using the sub-Chandrasekhar-mass WD models. We solve the two-dimensional Euler equations in cylindrical coordinates. We further assume reflection symmetry of the x–y plane so that we model only one quadrant of the star. We use a realistic Helmholtz equation of state (Timmes 1999) for describing the matter with free electrons with an arbitrary relativistic level and degeneracy, nuclei as a classical ideal gas, photons with Planck distribution, and electron–positron pair effects. Different from the near-Chandrasekhar-mass WD models, we include

(1) the notation of He detonation in the simulation,

(2) its energy-generation prescription and timescale, and

(3) its propagation velocity.

We use the same solver for matter in nuclear statistical equilibrium (NSE), as it is independent of the original composition of the matter and only depends on the final density, temperature, and electron fraction. Level-set methods are used for tracing the contour of the He- and C-detonation fronts. In Paper II, we further performed a set of tests to justify the functionality of the code in the Appendix. They include that

(1) the C-detonation trigger is independent of the symmetry boundary we used, ⁷

(2) the explosion energetics are insensitive to the resolution used,

(3) the shock convergence is less sensitive to the resolution, and

(4) the threshold of the detonation trigger is independent of the resolution.

We further showed in Paper II that our two-dimensional models give results that agree with contemporary one-, two-, and three-dimensional models found in the literature.

In order to keep track of the detailed nucleosynthesis for constructing the isotope distribution, we use the tracer particle scheme (Travaglio et al. 2004; Seitenzahl et al. 2009; Townsley et al. 2016). It makes use of massless Lagrangian tracers. They follow the fluid motion and record their own thermodynamical trajectories. The particles are "massless" in that they do not change the fluid motion. After the hydrodynamical simulations, the tracers are postprocessed with a large nuclear reaction network (a 495-isotope network containing isotopes from ¹H to ⁹¹Tc; Timmes et al. 2000). After postprocessing, the spatial distributions of specific isotopes, such as ¹⁶O, ²⁸Si, and ⁵⁶Ni, are extracted.

In this work, based on the formalism of our previous works, we examine further SN Ia models in the range between 0.95 and 1.0 M_⊙, and we extend systematically to different initial He-detonation structures. We evolve WDs from the onset of He detonation until no significant exothermic reactions take place. This can be because the WD is disrupted completely by both types of detonation, or the He detonation fails to trigger the second detonation and quenches. In Table 1, we list the initial densities and temperatures used in our models. As indicated by observational data, we consider the He-envelope mass M_He from 0.05 to 0.10 M_⊙.

In the table, we also list other related parameters, including the central density ρ_c and interface density ρ_He. Among all of our considered models, the densities range from  ∼10⁷ to 10⁸ g cm⁻³, while the He interface is from  ∼10⁵ to 10⁶ g cm⁻³. Notice that this allows a major part of the star to carry out complete burning from CO matter to ashes in NSE at a density  >5 × 10⁷ g cm⁻³.

In general, three types of events can result. In the table, N/A means that no C detonation occurs; this happens when the He detonation is too weak (without the possible geometric convergence) to heat the CO matter to a sufficient temperature, or collision such that the shock propagates into the core. The results "Cen." and "Off" stand for centered and off-center detonation, respectively. We also list the yielded ⁵⁶Ni and ⁵⁷Ni masses, obtained at the end of the simulations, where most exothermic reactions have ceased.

In Figure 1, we depict the four scenarios used in this work. They include a bubble (B type), a ring (R type), a bubble and a ring (D type), and a sphere (S type). This spans the possible initial He runaway from the lowest to the highest symmetry. In the figure, the cross sections of the WD progenitor are drawn. The B-type detonation corresponds to two bubbles, one at the north pole and one at the south pole.

Zoom In Zoom Out Reset image size

In order to realize the one-bubble event, simulations explicitly modeling the hemisphere are necessary. We recall that the C detonation that is triggered is very similar to the R type, since, in this configuration, the C detonation always starts after the shock converges at the opposite "pole" from where the detonation is initialized. Therefore, we may refer to R-type series to trace how the detonation takes place. We note that, computationally, the detonation trigger is identical, as indicated by Appendix B in Paper II.

In this work, we focus on the asymmetry of the sub-Chandrasekhar-mass WD. We do not focus on the near-Chandrasekhar-mass WD because we observe that typical SN Ia models from the near-Chandrasekhar-mass WD tend to have a more spherical structure. Here we present a comparison to outline the large-scale asymmetry in both models. For the details of these models, we refer interested readers to Papers I, II, and III for more detailed descriptions and implementation.

The two models are chosen from our benchmark models. Both models are chosen to represent the "normal" SN Ia by its ⁵⁶Ni production,  ∼0.6 M_⊙. Furthermore, the near-Chandrasekhar-mass model also produces the necessary amount of ⁵⁵Mn for reproducing the trends of stellar abundance near the solar neighborhood.

In the top panel of Figure 2, we show the temperature in a section through the star of the near-Chandrasekhar-mass model at ∼4 s after the explosion. The model assumes turbulent deflagration with a deflagration–detonation transition. The initial nuclear runaway assumes a centered deflagration wave with angular perturbations. The near-Chandrasekhar-mass WD has a mass of 1.37 M_⊙ and a central density of 3 × 10⁹ g cm⁻³. A centered-flame c3 (with three-finger structure) is used as the initial nuclear runaway. The finger structure is used to enhance hydrodynamical instabilities for the asphericity. The sub-Chandrasekhar-mass model has a total mass of 1.1 M_⊙ with a helium envelope mass of 0.1 M_⊙ and a central density of ∼6 × 10⁷ g cm⁻³. A one-bubble configuration is placed along the rotation axis.

Zoom In Zoom Out Reset image size

For our model, at t ∼ 4–5 s, the global distribution of ejecta structure begins to be frozen. Secondary features, including Rayleigh–Taylor instabilities ("mushrooms"), still continue to grow. But most features are smoothed by the expansion, leaving a surface close to spherical. We can observe that by the external detonation transition, the detonation wave always wraps around the aspherical ash produced during deflagration. Surface asphericity in temperature appears to be significant in the Chandrasekhar-mass model at a radius of ∼30,000 km. However, the matter at the surface is mostly C and O. The area of interest, where iron-peak elements are synthesized, does not show much asphericity. They are produced primarily in the core, which means it is unlikely that large-scale asymmetry features from ⁵⁶Ni and ⁵⁶Co can be generated and exposed. In this work, the core corresponds to the CO-rich matter under the He envelope.

In the bottom panel of Figure 2, we show the temperature profile of a sub-Chandrasekhar-mass CO WD with an He envelope for comparison. The sub-Chandrasekhar-mass model shows a more explicit large-scale asphericity. The detonation is triggered near the "equator." As a result, the He detonation is stronger near the equator, as shown by the hot spot 60,000 km closer to the symmetry axis. Since the C detonation propagates from the place near the equator to the center and then outward, the inner ejecta from the C detonation has a preferred direction along the rotation axis. The high-velocity flow along this direction will be responsible for the later large-scale asphericity.

Having shown that the sub-Chandrasekhar-mass WD is capable of generating large-scale asymmetry, we further examine how the distribution of ⁵⁶Ni deviates from spherical symmetry. In Figure 3, we show the final distributions of the tracer particles, which are the mass fraction contours of ¹⁶O, ²⁸Si, and ⁵⁶Ni. As examples, we consider models 100-100-R, 100-100-B, and 100-100-S.

Zoom In Zoom Out Reset image size

For model 100-100-S, we can see that all of the Ni-rich matter is covered by the thick envelope. We choose this series because all of these models produce comparable ⁵⁶Ni mass as a normal SN Ia.

The S-type model does not show any significant amount of ⁵⁶Ni near the surface. The distribution confirms that the initial spherical symmetry makes the detonation wave propagate only spherically. Except for small-scale perturbation coming from discretization effects from cylindrical to spherical coordinates, the large-scale distribution is, to a good approximation, spherical. The detonation waves do not collide with each other except at the stellar core when the detonation waves converge. There is no detonation wave interaction by shock collision or convergence in the He envelope. The amount of ⁵⁶Ni production becomes lower in the envelope. Hence, this model is less likely to explain the observed early ⁵⁶Ni decay line.

The R- and B-type models show more near-surface ⁵⁶Ni. The ⁵⁶Ni is ejected more along the polar direction for the B-type detonation; meanwhile, the R-type detonation ejects matter more spherically but in a more elongated manner compared with S-type model. The D-type model also behaves similarly to the R-type model, except that the ⁵⁶Ni around the equator is lower in abundance.

We recall that it is the asynchronous burning of helium, coupled with the geometric convergence, that creates the observed asphericity. The asynchronous burning allows He with the same initial density to be burned at different times. When the detonation time is delayed, the density of the He matter decreases, thus making it less likely to generate ⁵⁶Ni in the first place. On the contrary, the geometric convergence creates the hot spot for robustly triggering ⁵⁶Ni synthesis. However, the effect of geometric convergence is more localized. In a similar way, the detonation wave collision can create the high-temperature and high-density zones for synthesizing ⁵⁶Ni.

One of the features in SN 2014J is the early-time ⁵⁶Ni signal and time-dependent velocity for late-time ⁵⁶Co (see Section 4). It is therefore interesting to further study how the initial detonation configuration can give rise to the diversity of the isotope velocity distributions. In Figure 4, we plot the ejecta distribution in the velocity space to analyze how the detonation affects the final ejecta distribution. The distribution is angular-averaged. More samples of SNe Ia show similar early premaximum bumps (see, e.g., Jiang et al. 2017, 2018). The possibility of these SNe Ia forming a subclass shows that SNe Ia with observable He-burning features can have a common evolutionary path. Furthermore, such early gamma-ray flux can be another important sign for future telescopes to capture the early optical evolution of these SNe (Wang et al. 2019).

Zoom In Zoom Out Reset image size

For model 100-100-B, the one-bubble configuration allows ejecta to be concentrated with ⁵⁶Ni-rich material (see Figure 3) near the equator. As a result, there is a multilayered distribution of ⁵⁶Ni. Below 6000 km s⁻¹, the ejecta are filled with ^56-58Ni. From 6000 to 8000 km s⁻¹, ²⁸Si and ³²S are the major isotopes. For 8000–11,000 km s⁻¹, unburned oxygen is the major element, and outside He is in the main element in the ejecta. We remark that the ⁵⁶Ni distribution is not monotonically decreasing, as compared to the classical spherical model. It first drops around 10,000 km s⁻¹, showing that the detonation reaches the low-density region for C burning. After that, the mass fraction of ⁵⁶Ni rises again and quickly drops in its abundance. This shows that the shock collision in the He envelope allows formation of ⁵⁶Ni directly. But the shock is not strong enough to channel the outburst of ⁵⁶Ni, as seen by the covering layer of ⁴He.

For model 100-100-D, the strong collision away from the axis allows an outburst of ⁵⁶Ni during ⁴He burning at early time. This is also reflected in the ejecta distribution. Ejecta with a velocity below 10,000 km s⁻¹ are again filled with Ni isotopes. Here ²⁸Si and ³²S are the major isotopes in the velocity range 10,000–11,000 km s⁻¹. Products of incomplete C burning, such as unburned ¹⁶O, can be found most abundantly up to 12,000 km s⁻¹. Outside that, ⁵⁶Ni and ¹⁶O are the major isotopes. From this, it can be seen that multiple ignitions allow ⁴He to be burned quickly along the α-chain. Such outermost ⁵⁶Ni can be readily ejected and seen through its decay.

For model 100-100-R, the geometric convergence takes place at the pole that is strong enough to create similar pinching to the He envelope. The ejecta are covered with intermediate mass elements (IMEs) from 10,000 to 11,000 km s⁻¹ and ¹⁶O from 11,000 to 12,000 km s⁻¹. The outermost layer is mixed with ⁵⁶Ni and ⁴He with hints of ²⁸Si and ³²S. Similar to model 100-100-D, there is a thin layer of IMEs between 10,000 and 11,000 km s⁻¹ and then incomplete C-burning products between 11,000 and 12,000 km s⁻¹. The outermost ejecta are a mixture of ⁵⁶Ni, ⁴He, ²⁸Si, ³²S, and ¹⁶O.

For model 100–100-S, the spherical detonation allows a stratified structure in the ejecta. Ejecta with a velocity below 12,000 km s⁻¹ are dominated by ⁵⁶Ni, ⁵⁷Ni, and ⁵⁸Ni. The IMEs, including ²⁸Si and ³²S, are mostly found between 12,000 and 14,000 km s⁻¹. Ejecta with a velocity  >14,000 km s⁻¹ are occupied by ⁴He.

Comparing these four models demonstrates the possibility of mixing ⁵⁶Ni at high velocity, namely at the outermost ejecta for the D and B types. The oblique shock and geometric convergence of detonation inside the He envelope can provide the necessary thrust for channeling the ⁵⁶Ni produced in the He detonation to the outermost part of the ejecta. When the ejecta quickly expands, the gamma rays produced by the decay of ⁵⁶Ni into ⁵⁶Co may be directly seen after the surface matter becomes optically thin.

In the previous section, we examine the angular average ejecta velocity composition, and we show that the ejecta composition depends on the explosion geometry. Here we further analyze the ejecta composition by choosing specific angular slices. In particular, we choose the angular slices at 0°–9° from the rotation axis and from the symmetry axis (i.e., 81°–90° from the rotation axis) to contrast the ejecta composition. This demonstrates the difference between metal production with and without geometric convergence.

3.4.1. Bubble-type Explosion

The one-bubble configuration has been shown to be weaker than its one-ring counterpart due to the lack of shock convergence. The explosion is weaker with lower ⁵⁶Ni production. In the top left and right panels of Figure 5, we plot the ejecta distribution in the velocity space of model 100-100-B for two angular slices near the rotation axis (left) and symmetry axis (right).

Zoom In Zoom Out Reset image size

The two slices show very different ejecta structures. The ejecta along the rotation axis is slightly faster than the rotation axis. This is because the initial detonation triggered along the poles. However, we remark that this does not mean there is more substance along the axis because the velocity space does not directly follow the mass coordinate. The iron-peak elements in the core also differ a lot. Along the poles, the representatives of iron-peak elements, including ^56,57,58Ni, occupy the innermost ∼7000 km s⁻¹, and then the abundance quickly drops off. On the other hand, the iron-peak element-rich core extends up to 12,000 km s⁻¹ along the equator. The large difference comes from the C detonation. It is triggered near the equator. As a result, the C detonation first burns the material along the equator and reaches the center, and then it burns the matter along the equator. There is more time for the matter to move outward and expand, thus yielding a weaker heating effect.

The IMEs form the middle layer, from 7000 to 17,000 km s⁻¹ along the poles and from 12,000 to 15,000 km s ⁻¹ along the equator. Along the poles, as the velocity increases, which corresponds to lower-density matter, some unburned ¹²C can be seen. However, there is no such trace along the equator. Again, this demonstrates that the detonation along the poles is weaker than that along the equator due to the time lapse during expansion.

From 17,000 to 21,000 km s⁻¹ along the poles, ⁵⁶Ni dominates the ejecta again. These are the products of the He detonation, as the composition is mostly ⁴He. They share similar mass fractions up to the surface. On the other hand, along the equator, there is more ⁵⁶Ni from 15,000 to 19,000 km s⁻¹ and more ⁴He from 19,000 km s⁻¹ onward.

3.4.2. Ring-type Explosion

3.4.3. Model Comparison

The explosion of a sub-Chandrasekhar-mass WD is known to be sensitive to the initial WD mass, M. This is because the WD is degenerate, so the central density varies largely from 10⁷ g cm⁻³ (for M ∼ 0.95 M_⊙) to 10⁹ g cm⁻³ (for M ∼ 1.20 M_⊙). The difference in the central density corresponds to a difference in the average density of WD matter, where the CO matter can be completely burned to ⁵⁶Ni when it has a typical density  of ∼5 × 10⁷ g cm⁻³. As a result, the amount of final ⁵⁶Ni drastically varies by a factor of 2–6 when M increases from 0.9 to 1.2 M_⊙.

In Figure 7, we show the final ⁵⁶Ni mass against M for the four different detonation morphologies described in the previous subsections. Here the magenta box represents the range of ⁵⁶Ni derived from the gamma-ray data from SN 2014J. See Section 4 for the application to SN 2014J. Models in which the C detonation cannot be triggered are omitted, i.e., 095-050-B, 095-050-B, 100-050-B, 105-050-B, 110-050-B, and 095-050-D. All of these show a typical He nova event without the C detonation (Kippenhahn & Thomas 1978; Piro & Bildsten 2004), which is inconsistent with the supernova observation that the whole star is disrupted after the explosion. From both figures, we can see that the S series is the strongest, and then the D and R series. The B series is the weakest for the same initial M.

The possibility of triggering the C detonation in the low M_He limit relies on the detonation symmetry. The CO WD models with a lower mass from M = 0.95 to 1.10 M_⊙ show that with the lowest symmetry (B type), there is no geometric convergence. The only shock collision occurs when the detonation reaches the equator of the WD. The assumed boundary condition (reflection symmetry) allows the arriving shock waves to collide in a laminar way. As indicated in Paper II, when M increases, the minimum M_He required for triggering the second detonation decreases. According to similar work without assuming reflection symmetry (e.g., Fink et al. 2014), the minimum necessary M_He drops from 0.126 to 0.0035 M_⊙ when the CO WD mass increases from 0.810 to 1.385 M_⊙. Both works give us an insight that a higher M_He is necessary to provide sufficient shock strength in triggering the C detonation.

By only considering the WD mass, a higher M means that the transition from the CO core to the He envelopes takes place at a higher density. This increases the typical reaction rate and hence the energy production. The postshock temperature in the He envelope is, therefore, higher for more massive WDs, where the burned matter can reach the threshold temperature easier, independent of additional geometrical convergence. Thus, a higher-M model favors the trigger of the second detonation.

We recall that in simulations using a hemisphere of a WD, the dependence on the M_He is stronger than in simulations using a quadrant. This is because the detonation starts from one pole, and then the detonation wave wraps over the He envelope and converges at the other pole. In this situation, a shock convergence similar to the R-type detonation always happens. As shown in the table for the R-type detonation, the corresponding minimum M_He for the second detonation is lower.

For D, R, and S, they have a higher symmetry where there is geometric convergence by means of oblique shock, two-dimensional shock convergence (from a ring to a point), and three-dimensional shock convergence (from a sphere to a point), respectively. The resultant temperature in the CO core can be much enhanced by the converged shock strength. Therefore, the minimum M_He required to trigger the second detonation is more relaxed.

The difference in the final ⁵⁶Ni mass for the same M at different He detonation is related to the propagation of the C-detonation direction. The S model is always the strongest because the C detonation begins at the center and propagates outward, so that most of the star remains approximately static before the detonation wave arrives. This ensures that the matter remains less expanded and hence maintains a higher density, which results in a stronger explosion. On the contrary, in the B, D, and R models, the off-center C detonation means that the C detonation has to overcome the density gradient in order to reach the high-density matter in the core. This implies that the relative explosion strength is weaker because of the density gradient.

The early ⁵⁶Ni in the ejecta near the surface provides distinctive hints of the asphericity of the detonation model. Furthermore, the distribution of iron-peak elements, which are shock-heated by the backward shock in the ejecta, can reach reach a temperature of  ∼10⁶ K, which is sufficiently hot for X-ray emission for further diagnosis. This will give further constraints on the explosion models. For example, the X-ray spectra of the SNR 3C 397 have been used as a diagnosis of a Chandrasekhar-mass WD progenitor (Yamaguchi et al. 2015; Dave et al. 2017; Leung & Nomoto 2018).

In Figure 6, we plot the element distribution of Cr, Mn, Fe, and Ni for models 100-100-B, 100-100-R, 100-100-D, and 100-100-S. We assume that after the star reaches homologous expansion, the distribution of ejecta remains mostly unchanged. However, we also remark that the radiation energy during the decay of radioactive isotopes can still trigger inner motion and affect the element distribution. But the effect is secondary compared to the initial distribution during explosion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The distributions of the iron-peak elements show more diversity than the major elements, as shown in the previous section. From the contour shape, we observe two features. First, Cr and Mn almost follow each other, and Fe and Ni follow each other. Second, Cr and Mn tend to have a more spherical distribution, while Fe and Ni follow the detonation geometry more closely.

Model 100-100-R shows the largest deviation from a spherical structure among all four elements. Model 100-100-B has spherical Cr and Mn but aspherical Fe and Ni. Similar characteristics appear for model 100-100-D. All elements are spherical in model 100-100-S, as anticipated by the initial spherical symmetry. To distinguish between models 100-100-B and 100-100-D, we notice that the distribution of Cr is more irregular in model 100-100-D compared to the quasi-spherical distribution in model 100-100-B. However, the difference is subtle.

The closest SN Ia in the last four decades, SN 2014J is a special example of SNe Ia exploded in the nearby galaxy M82, just 3.3 Mpc away. Its closeness to the Milky Way has provided the chance for detailed multiband observations, including the radio (Pérez-Torres et al. 2014), infrared, optical (Goobar et al. 2014; Kawabata et al. 2014), UV (Foley et al. 2014), X-ray (Terada et al. 2016), and gamma-ray (Diehl et al. 2014; Churazov et al. 2015; Diehl et al. 2015; Diehl 2015; Siegert & Diehl 2015; Isern et al. 2016) bands, with its spectra at early and late times (Ashall et al. 2014; Jack et al. 2015; Dhawan et al. 2018).

Measurements of SN 2014J have been made in different works in the literature. In Churazov et al. (2014), the estimated ${M}_{{}^{56}\mathrm{Ni}}=0.56{\pm }_{0.06}^{0.14}\,{M}_{\odot }$, and the estimated ejected mass is $1.2{\pm }_{0.5}^{1.9}\,{M}_{\odot }$. A similar measurement is found in Diehl et al. (2015), who gave ${M}_{{}^{56}\mathrm{Ni}}=0.50\pm 0.12\,{M}_{\odot }$. For ⁵⁷Ni, Yang et al. (2018) reported that ⁵⁷Ni/⁵⁶Ni has a mass ratio of $0.065{\pm }_{0.004}^{0.005}$ based on the B-band maximum light and $0.066{\pm }_{0.008}^{0.009}$ based on the pseudobolometric light curve. Stable Ni mass is constrained at 0.053 ± 0.018 M_⊙ (Dhawan et al. 2018).

It shows gamma-ray features that agree well with the classical spherical pure deflagration model W7 (Churazov et al. 2014; Diehl et al. 2015). In photometry, SN 2014J appears to be a normal SN Ia (Isern et al. 2016) that shows a comparable structure with the W7 model (Nomoto et al. 1984) and a similar abundance profile as the heavy elements (e.g., ⁵⁶Ni) in the core and lighter elements (e.g., Si and S) in the envelope (Ashall et al. 2014).

However, detailed examinations of the observational data of SN 2014J reveal some differences from ordinary SNe Ia. For example, the rise of the UVOIR light curve with time (e.g., Nugent et al. 2011) shows its delay in SN 2014J. The early light curve of SN 2014J suggests a "shoulder" only a few days after the first light (Goobar et al. 2014). The late-time evolution (beyond a few hundred days) shows derivations from classical SNe Ia such as SN 2011fe, where the slower decline rate suggests interactions with CSM (Foley et al. 2014; Yang et al. 2018). The ultraviolet data of SN 2014J show large extinction (Brown et al. 2015). Such extinction and CSM can be in the dusk disk structure (Nagao et al. 2017). Images around SN 2014J do not show an observable companion star, thus making its companion as a red giant unlikely (Kelly et al. 2014; Margutti et al. 2014; Pérez-Torres et al. 2014).

The progenitor of SN 2014J is still a question of debate (e.g., Margutti et al. 2014; Pérez-Torres et al. 2014; Dragulin & Hoeflich 2016; Graur & Woods 2019). We note that these constraints on CSM cannot be applied to the pre-SN environment of a uniformly rotating WD with a slightly super-Chandrasekhar mass in the single-degenerate scenario, as calculated by Benvenuto et al. (2015). The gamma-ray signal cannot distinguish with high significance which class SN 2014J belongs to (Terada et al. 2016).

According to further examinations, the observational data of SN 2014J show features that deviate from the spherical approximation as discussed below.

First, the early gamma-ray observations with INTEGRAL discovered the lines at 158 and 812 keV (T_1/2 = 6.6 days) that are characteristic for the ⁵⁶Ni decay around 17.5 days after the inferred explosion date (Diehl et al. 2014; Isern et al. 2016). The model fits of Isern et al. (2016) appeared more consistent with a redshifted and broadened ⁵⁶Ni emission, so they suggested ejection of ⁵⁶Ni-rich material in a blob moving away from the observer. The analysis by Diehl et al. (2014) was performed in finer energy bins and without any model bias; their sampling of possible spectral solutions suggests a narrow emission line from ⁵⁶Ni at the laboratory energy value, with indications of blue- as well as redshifted satellite. This led them to suggest a model with ⁵⁶Ni ejected perpendicular to the observer's line of sight. But the total significance of this surface ⁵⁶Ni line emission is only 3σ–4σ; hence, both interpretations remain possible.

In any case, the existence of ⁵⁶Ni near the surface is witnessed by the INTEGRAL data and is the main motivation for our study. The missing C and O absorption lines are possibly connected to the exposed ⁵⁶Ni (Goobar et al. 2014). Also, expected Ni lines are however not seen in the infrared at this time.

The second aspherical feature appears at later times. When the energy output is dominated by the decay of ⁵⁶Co (T_1/2 = 77.1 days), the measured Doppler shifts of the ⁵⁶Co decay lines show the expected behavior, plus an additional structure of (at least) three blobs of distinctive velocities (i.e., showing an early blueshift when the redshifted part is opaque in the line of sight and then becoming symmetric with no Doppler shift later, as pointed out in Diehl et al. 2015). Instead, at least three distinctive centroid energies could be identified. This may correspond to fluid parcels containing ⁵⁶Co being ejected with different velocities with respect to the observer's frame. This time-dependent variation of the Co decay line frequency, denoted as "flickering," is observed in SNe Ia for the first time. This suggests the possibility that "blobs," or large-scale asymmetries, developed during the explosion.

To summarize, SN 2014J has shown conflicts with the classical spherical model. These include (1) the observed Doppler-broadened ⁵⁶Co lines, which have an irregular appearance with time (Diehl et al. 2015); (2) hints of the ⁵⁶Ni decay lines on the surface (Diehl et al. 2014); and (3) an enhanced ionization on the outer part of the star, as revealed by the early atomic line spectra.

Among the above mentioned features of SN 2014J, the existence of ⁵⁶Ni near the surface suggested by the INTEGRAL data is the main motivation for the present study.

Before discussing the sub-Chandrasekhar mass models in the following subsections, we note that near-Chandrasekhar mass models have variations depending mainly on the mass accretion rate, and some models produce ⁵⁶Ni near the surface of WDs. In the classical picture of the near-Chandrasekhar mass model, ⁵⁶Ni is concentrated in the inner core (The & Burrows 2014). However, ⁵⁶Ni is produced near the surface in the following models.

Nomoto (1982b) showed several models where the WD mass increases to the near-Chandrasekhar mass with slow accretion of He. If the accretion rate is low enough, the accreted He is too cold to be ignited, thus being just accumulated near the surface. Eventually, the WD becomes massive enough to ignite the central C-deflagration.

In the late detonation model by Yamaoka (1992), the deflagration-detonation transition can occur in the outer layer of the near-Chandrasekhar mass WD and burns He to produce ⁵⁶Ni near the surface (see model W7DHE in Figure 3 of Yamaoka 1992). Such a deflagration-detonation transition might occur at the very steep density gradient near the WD surface even with a small amount of He. The deflagration-detonation transition near the surface would be likely to occur in an aspherical manner which might be interesting for further study.

Although the production of ⁵⁶Ni near the surface is possible in the late detonation of the near-Chandrasekhar-mass model (W7DHE), here we focus on the sub-Chandrasekhar-mass models in order to apply the results of this Paper IV to the "asphericity'' of SN 2014J.

The object SN 2014J was first proposed to be the ignition of an "He belt" accumulated in the orbital plane of the binary (Diehl et al. 2014). The He-belt model was motivated by the early ⁵⁶Ni decay line detected with a close-to-zero line Doppler shift. This picture can explain the origin of the early-time ⁵⁶Ni decay line. However, on top of that, in Diehl et al. (2014), such a decay line has a very small redshift. The He-belt model ejects Ni synthesized in the He envelope along the equator direction. ⁸ As a result, it provides the source of Ni with a small Doppler shift. But, it remains unclear if the He belt can be formed and stably maintained during the binary accretion (Kippenhahn & Thomas 1978; Piro & Bildsten 2004).

In the following discussion, we search for a qualitative model that may resemble the observed characteristics of SN 2014J. We note that the model asphericity should not be too strong, as otherwise it would violate the similarity of SN 2014J with the classical spherical model. Thus we set the He envelope in spherical form. To generate the required surface Ni, different He-detonation configurations are studied.

Our attempt is to elucidate the physical conditions that are able to qualitatively reproduce the distinctive features of SN 2014J. Owing to the complexity of multiple data constraints and the subtlety in the interpretations of different observations, we avoid scrutinizing an exact or complete model that can explain all of the quantitative features for SN 2014J. Our focus is, therefore, on the early ⁵⁶Ni line emission, the ⁵⁶Ni mass as derived from the peak (V-band) luminosity, large-scale asymmetries as implied by the flickering ⁵⁶Co decay line profiles, and the ⁵⁷Ni/⁵⁶Ni mass fraction ratio from the late-time light curve. Using ⁵⁷Ni/⁵⁶Ni to constrain the SN explosion models as proposed in Seitenzahl et al. (2009) has been done in other SNe Ia, for example, SN 2012cg.

To recapitulate, based on the two-dimensional models, we aim at searching for (1) WD parameters that correspond to the general features of SN 2014J and (2) appropriate He-detonation triggers that produce the global asymmetries and qualitative features of SN 2014J, including large-scale asymmetry and/or near-surface production of ⁵⁶Ni.

The ⁵⁶Ni mass observed in SN 2014J provides the important constraint on the progenitor WD mass. In Diehl (2015), ⁹ the ⁵⁶Ni mass is estimated to be 0.49 ± 0.10 M_⊙. From our simulation results in Table 1, we identify the possible mass range to be M = 1.0–1.1 M_⊙.

To account for the observed ⁵⁶Ni mass in SN 2014J shown by the magenta box in Figure 7, we require M = 0.95–1.00 M_⊙ for the "S"-series and M = 1.00–1.10 M_⊙ for the "B"-series (M_He = 0.10 M_⊙). "D"- and "R"-series require M = 1.05–1.10 (1.00–1.05) M_⊙ for M_He = 0.05 (0.10) M_⊙.

Another aspect by which to constrain SN 2014J is the explosion geometry. The gamma-ray signal of SN 2014J has suggested the (near-)surface ⁵⁶Ni. We examine which initial detonation geometry allows the formation of ⁵⁶Ni at these regions. To extract the final distribution of ⁵⁶Ni, we use the tracer particle data when the ejecta reaches homologous expansion and the structure of the ejecta is frozen out.

Future observations of SN 2014J, which will be able to disentangle the morphology, may thus provide clues to the initial configuration from measurements of different elements. For example, in Grefenstette et al. (2014, 2017), such a technique is pioneered in showing ⁴⁴Ti in Cas A to demonstrate how to disentangle and map the gamma rays.

In Figure 3, we show the distributions of some representative elements for models with different detonation geometry. For example, for models starting with an He-detonation bubble, the Ni and Si ejecta are in a cocoon shape, compared to the spherical shape in the spherical detonation. The B type is also highly distinctive due to the thick layer of O-rich ejecta, compared to the spherical counterpart. On the other hand, the D and R types are different from the other two by the quasi-spherical O- and Si-rich ejecta, while the Ni ejecta maintain an observable ellipticity. However, between the D- and R-type models, Si and Ni show a similar distribution.

We now examine the dependence of the large-scale asymmetry on the He-envelope mass. In Paper II, we presented a parameter survey on the nucleosynthesis yield of SNe Ia using the sub-Chandrasekhar-mass WD as the initial progenitor. It is shown that the He-envelope mass can strongly enhance the production of some Fe-peak isotopes, including ⁴⁸Ti, ^50,51V, and ⁵²Cr.

In Figure 8, we plot the ⁵⁶Ni distribution for models 100-050-R and 100-100-R. They differ by the mass of the He envelope from 0.05 to 0.1 M_⊙. The distribution of ⁵⁶Ni after the explosion can show the existence of large-scale asymmetry. In model 100-100-R, with a more massive He envelope, the He detonation is strong enough to drive the inward propagation of the detonation to the core; this results in complete disruption, where most of the matter in the CO core is spherical. Therefore, this scenario will be less likely to exhibit a flickering, as seen in SN 2014J. However, its ⁵⁶Ni also distributes close to the surface, which is an important feature of SN 2014J. As a result, the ejected ⁵⁶Ni around 10° from the rotation axis has a higher radial velocity. The lower He-envelope mass makes the production of ⁵⁶Ni lower, and the synthesized ⁵⁶Ni is covered by the original He envelope, which can block the gamma rays emitted by the radioactive decay. This means that a lower-M_He model has more difficulty reproducing the early gamma rays observed from SN 2014J.

Zoom In Zoom Out Reset image size

Following the expansion of the ejecta, when the matter becomes optically thin to gamma rays, the inner structure of the ⁵⁶Ni and ⁵⁶Co distribution will be exposed. The velocity fluctuation, as seen from Figure 8, depending on the ejecta angle and the exact time to become optically thin, may coincide with the flickering feature, as seen in SN 2014J.

In Figure 9, we plot the ⁵⁷Ni and ⁵⁶Ni yields (prior to decay) of our simulated SN Ia models. Shown are models for M_He = 0.05 (solid lines) and 0.10 (dashed lines) at solar metallicity. We also show the SN 2014J data (crosses; Yang et al. 2018). Apparently, the models listed there are not sufficient to explain the high ⁵⁷Ni value of this SN.

Zoom In Zoom Out Reset image size

Similarly, to demonstrate the effects of metallicity, we plot in the figure a red dashed–double-dotted line for model 110-050-B at 0, 1, 3, and 5 Z_⊙ (corresponding to the data points from the bottom to the top). To be consistent with the SN 2014J data, a model $Z=\,\sim$ 5 Z_⊙ is required. All of the trend lines for other models behave similarly. In Paper I, we have showed that the metallicity has a minor impact on the global energetic explosion. Instead, it shows its influence on the relative abundance ratio, in particular the high-Y_e isotopes including ⁵⁴Fe, ⁵⁸Ni, and ⁵⁵Mn.

From this figure, we observe that in order to explain the abundance pattern of SN 2014J, a high-metallicity model ≳4 Z_⊙ is necessary to be consistent with the high ⁵⁷Ni mass relative to ⁵⁶Ni.

The morphology of the SN ejecta and the shape of the SNR, such as Tycho (Ferrand et al. 2019), may directly link to the initial explosion configuration. By tracing the line emission of shock-heated ejecta, the abundance of the measured elements can be revealed. In Seitenzahl et al. (2019), the tomography of three very young SNRs, 0519-69.0, 0509-67.5, and N103B, in the Large Magellanic Cloud (Hughes et al. 1995) are studied for the first time. By examining the S xii, Fe ix, and Fe xv lines, they reconstructed the large-scale distributions of these elements in these remnants. These objects are sufficiently young such that the shocked-heated matter remains clearly visible and the shock front has not completely swept through the remnant.

The SNR 0519-69.0 has a more spherical shape but with small-scale perturbations on the surface, as depicted by the X-rays. From our simulations, such features are possible when we consider the Rayleigh–Taylor instabilities that freeze out during the expansion of the ejecta. As discussed in Seitenzahl et al. (2019), this remnant can be fitted by a Chandrasekhar-mass WD of mass 1.4 M_⊙. Their total Fe mass is about 0.4 M_⊙. This belongs to the lower side of SN Ia production (Leung & Nomoto 2018). Such features can be obtained for a higher-mass progenitor of ∼1.4 M_⊙ with an initial central density >3 × 10⁹ g cm⁻³.

The SNR 0509-67.5 also has a spherical shape where the ejecta, including the X-ray emitting part, demonstrates a close-to-spherical structure. The estimated mass of this remnant is ∼1.0 M_⊙ with 0.15 M_⊙ He. To achieve such a spherical shape, a spherical He detonation is necessary avoiding any large-scale asymmetry created during the C detonation. In fact, when such a heavy He envelope is included, spherical detonation is preferred because the expected nuclear runaway time can easily be shorter than the convection timescale, thus triggering simultaneous burning in the spherical layer near the He–CO interface. Also, as discussed in Paper II, the thick He envelope >0.1 M_⊙ can lead to a very severe excess in light iron-peaked elements such as Ti, V, and Cr. Future detection of these elements may provide further confirmation of this explosion picture. At last, its total Fe mass,  ∼0.5 M_⊙, can be mapped consistently with our models, such as 100-100-S.

The SNR N103B also has an aspherical ejecta shape. However, not much analysis of this object is reported in Seitenzahl et al. (2019). Despite that, from the morphology of the ejecta, an arrow or cone shape of Fe xiv distribution can be observed. Such an aspherical shape with a pointing effect can indicate the focused shock in the CO core and its consecutive breakout, as hinted from Figure 2. In our results, due to the reflection symmetry, the cone shape is always paired on both sides. We expect that, if we allow a single bubble without assuming reflection symmetry, a one-sided feature can emerge, as in more general simulations such as Tanikawa et al. (2019) and Gronow et al. (2020).

Another possibility for tracing the remnant progenitor is to examine the abundance patterns. As described in previous sections, the X-ray spectra from SNRs can provide essential clues to the relative amount of elements, especially iron-peaked elements including Cr, Mn, Fe, and Ni (Yamaguchi et al. 2015; Martínez-Rodríguez et al. 2017). They can directly constrain the type of WD progenitor and the explosion physics. In their work, the SNR spectra, including Tycho, Kepler, 3C 397, and G337.2-0.7 in the Milky Way and N103B in the Large Magellanic Cloud, are analyzed. It is suggested that, to explain the variety of remnant chemical abundance patterns, a factor beyond neutronization (i.e., metallicity tracer by ²²Ne) is necessary. Here we examine how the metallicity affects the elemental distribution and ⁵⁶Ni.

In Figure 10, we plot the mass and Z dependence of ⁵⁶Ni production for the sub-Chandrasekhar-mass WD models with M = 0.9–1.2 M_⊙ and Z = 0–0.04. The ⁵⁶Ni production is strongly M-dependent and monotonically increasing with M. The transition from complete to incomplete burning shows the strongest effects between 0.9 and 1.0 M_⊙ by a change of almost a factor of 10 within this mass range. The effects of metallicity can be seen but are small compared to the effects of M. In general, an  ∼10% difference can be seen among the Z range examined here.

Zoom In Zoom Out Reset image size

In Figure 11, we plot the mass fraction ratio of [Mn/Fe], [Ni/Fe], [Cr/Fe], and [Ca/S] for our sub-Chandrasekhar-mass WD models with a mass of 1.0 M_⊙ for metallicity Z = 0, 0.002, 0.01, 0.02, and 0.04. The metallicity effect is much larger. Major iron-peaked elements such as [Mn/Fe] and [Ni/Fe] differ by ∼0.6 dex between models with Z = 0 and 0.04. Both of them increase when Z increases. On the other hand, minor iron-peaked elements such as [Cr/Fe] and IMEs such as [Ca/S] show much smaller variations by only 0.1 dex in the Z range examined. They decrease when Z increases.

Zoom In Zoom Out Reset image size

However, we remark that the exact values of [Ca/S] and [Cr/Fe] are more prone to systematic uncertainties. This is because the light iron-peaked elements, including Ca and Cr, can be produced in both NSE and incomplete Si burning. For Mn, Fe, and Ni, they are robustly produced in NSE matter, where complete and instantaneous energy release can be assumed. On the other hand, for Ca and Cr, the numerical scheme of low-density burning can affect the abundance of these elements. For the same amount of energy release, a longer energy deposition time can slow down the thermal expansion. As a result, the matter has more time to carry out a slow nuclear reaction in the α-rich freeze-out regime before the ejecta expands and becomes too cold for any significant nuclear reactions. The case for low-density matter contains complications because the actual reaction depends on the detailed chemical composition, which is not well traced in multidimensional composition. The simplified chemical composition (seven-isotope network) may not provide an accurate estimation of how fast those reactions and the associated energy production take place.

Our work has suggested that He detonation is the key feature for explaining some well-observed SNe Ia, such as SN 2014J. However, it also remains to be understood in detail how the ignition of surface He can be coupled to the current explosion mechanism that centrally ignites carbon in the WD. Three-dimensional low Mach number hydrodynamical simulations of the He-burning envelope are necessary for realizing how the initial hydrostatic He burning develops into nuclear runaway. The detonation size is typically assumed to be as large as the pressure scale height (Bildsten et al. 2007; Shen & Bildsten 2009; Dan et al. 2014). On the other hand, temperature fluctuations tend to favor the ignition at a single spot, which has a much smaller size. A stringent limit on the temperature fluctuation appears because of its small size (Holcomb et al. 2013). Contamination of C-rich matter from the CO core can be an alternative to decreasing the necessary size of nuclear burning (Shen & Moore 2014).

The prerunaway phase of the double-detonation model was unclear until recent large-scale works on clarifying the possibility of such a proposal (Jacobs et al. 2016). Depending on the convection flow, different He-runaway pattern scenarios could occur from the most nonspherical extreme, i.e., a bubble, to the most spherically symmetric case. The modeling of such processes is typically much longer than the hydrodynamics timescale in order to directly capture the first nuclear runaway from nuclear reactions. To resolve the first runaway, the simulation requires the hydrodynamics timescale t_hyd to be smaller than the timescale of nuclear burning t_burn and convection t_conv, i.e., t_hyd < t_burn < t_conv (Glasner et al. 2018). Only recently have there been a few pioneering models using hydrodynamics simulations in the low Mach number regime to follow how the convection develops into nuclear runaway self-consistently (Zingale et al. 2013; Jacobs et al. 2016). Therefore, to trace the possible origin and explain the observations of SN 2014J, we considered different possibilities for how the He detonation can be triggered.

On the other hand, once a detonation spot forms, the second detonation is generally inevitable. Recent three-dimensional large-scale hydrodynamics simulations of one- or multiple-spot He detonation has been found to be robust in triggering the off-center detonation by geometric convergence in a quiet He envelope (Moll & Woosley 2013). The He-envelope mass required to trigger the second detonation can drop significantly from  ∼0.1 to  ∼10⁻³ M_⊙ for a WD mass increasing from 0.8 to 1.3 M_⊙ (Fink et al. 2010), while observationally, an He envelope below 0.05 M_⊙ is favored due to the discrepancy with the theoretical light curve when a massive He envelope is applied (Woosley & Kasen 2011).

Such diversities in the He-detonation trigger and detonation wave interactions have provided the flexibility to account for the diversity of SNe Ia. In Diehl et al. (2015), a multiple-plume structure is proposed to illustrate the apparent flickering of the ⁵⁶Co decay line. Based on our models, it is possible that such a feature can be realized by multiple spots in the He envelope. For a quantitative comparison, three-dimensional models are required for a one-to-one matching of the observables and predicted signatures. Nevertheless, using the virtue that geometric convergence and laminar wave shock collision do not differentiate between two- and three-dimensional simulations, our models can shed light on what kind of shock interaction, and hence what kind of detonation pattern, is necessary for reproducing features taken from SN 2014J.

In this work, we explore the possible triggering and ignition mechanisms that might lead to the asymmetric properties demonstrated by the observational features of SN 2014J. Among all models, the closest model we obtain is model 110-050-R, with a metallicity  ≳4 Z_⊙.

The total (WD) mass required is 1.0–1.1 M_⊙ at the intermediate-to-high-mass end of a sub-Chandrasekhar-mass WD. Stellar evolution theory suggests a progenitor mass constrained in a range likely to be ∼6–7 M_⊙ (see, e.g., Catalán et al. 2008; Doherty et al. 2015 for the progenitor final mass relation). ¹⁰ For a sub-Chandrasekhar-mass WD to produce 0.6 M_⊙ ⁵⁶Ni as a normal SN Ia, theory suggests an initial mass of  ∼1.0 M_⊙ in one-dimensional models (Shigeyama et al. 1992; Nomoto & Leung 2018; Shen et al. 2018). The ⁵⁶Ni mass production is known to be sensitive to the progenitor mass because of the density dependence of ⁵⁶Ni production, with a minimum  >5 × 10⁷ g cm⁻³. Hence, a higher WD mass ensures not only a higher central density but also a higher energy release at the center, which favors the propagation of the detonation wave. In the one-dimensional models, the 1.1 M_⊙ case gives rise to a bright SN Ia for its  ∼0.8 M_⊙ production of ⁵⁶Ni. However, in multidimensional models, there are variations in the ⁵⁶Ni production based on the initial He-detonation structure. The aspherical He detonation tends to give a lower ⁵⁶Ni due to off-center ignition of C + O detonation because of the density gradient, as discussed in previous sections.

The He-envelope mass required by our calculations is ∼0.05 M_⊙. This is a marginal value for the He detonation to be observed (Woosley & Kasen 2011), where the optical observational features remain compatible with the data of normal SNe Ia. A high He-envelope mass likely overproduces some iron-peak elements, including Cr and V, near the surface. This changes the typical isotope and element abundance distributions compared to normal SNe Ia, where such elements are produced in a deeper layer, e.g., the Chandrasekhar-mass WD with a deflagration–detonation transition. Such high-opacity material in the envelope may make the explosion appear redder in optical spectra (Polin et al. 2019). However, their results are based on a one-dimensional model where the He-rich matter is always burned from high to low density. Therefore, the high-density matter always has a longer time to carry out nuclear reactions, which favors the production of such iron-peak elements. However, for aspherical detonation, this is not always true. The Cr and V production depends on how the He detonation spreads around the He envelope.

Another theoretical uncertainty is the exact He mass when the first nuclear runaway starts. The exact He-envelope mass depends on the mass accretion rate and the type of binary system (single- or double-degenerate). A higher M_He is more likely from the double-degenerate scenario, while a lower M_He is more likely from the single-degenerate scenario. Such a calculation was done by Kawai et al. (1988) from the stellar evolutionary perspective. The steady-state accretion of He on C + O and O + Ne + Mg WDs is investigated in the single-degenerate scenario. It is shown that the He-envelope mass drops sharply with the C + O core mass, with  ∼10⁻² M_⊙ for a 0.7 M_⊙ C + O core down to  ∼10⁻⁶ M_⊙ for a 1.36 M_⊙ core. This shows that the steady-state accretion in the single-degenerate scenario may not provide a robust way for accumulating an He envelope beyond 10⁻² M_⊙ in a WD of mass 1.00 M_⊙ or above.

The detonation required at the beginning is from a ring around the equator. Such a configuration is shown to produce more aspherical features in the Ni distribution, which would be compatible with the multiple redshifted ⁵⁶Co decay lines measured in SN 2014J. But we have also found that other types of initial detonation, such as the D and B types, may also produce similar characteristics, although less pronounced. How the He detonation is initialized is a matter of debate.

For a one-dimensional model, an entire mass shell is ignited simultaneously because of the assumed symmetry. However, it is unclear whether such symmetry can be maintained prior to the ignition. For example, in the single-degenerate scenario, the accretion of matter from the companion star through Roche lobe overflow generally has a high angular momentum. Such rapidly rotating matter, when accumulated on the stellar surface, may create strong dragging, which disturbs the material near the surface. Also, the He burning near the core–envelope interface may trigger convective motion (see, e.g., Jacobs et al. 2016). This creates a highly turbulent background due to the shear between the quasi-static C + O core and the rapidly rotating He-rich matter.

Knowing that the runaway of He is highly temperature-sensitive, it is conceivable that the ignition may occur at random locations and spherical symmetry is broken. In the most extreme case, only one spot can be ignited, which corresponds to the B-type explosion. If the rotation symmetry may be preserved, then an R-type explosion is one of the possibilities. However, the exact configuration will be best estimated from the detailed multidimensional hydrodynamics simulations for the last minutes before the runaway to capture how all of these processes interfere with each other (see, e.g., Zingale et al. 2011; Malone et al. 2014).

The early ⁵⁶Ni signal can act as a tracer of the explosion mechanism. In Figure 5, we show that how the He detonation is initiated can strongly influence the surface ejecta composition and its angular dependence. In particular, the early low-redshift ⁵⁶Ni line implies the possibility that we are observing SN 2014J close to where the shock is initialized. If we observe the shock converging point, the abundant elements of ²⁸Si and ⁴He may easily block the gamma rays. The later ⁵⁶Co line can have an origin in multiple shock convergence on the He envelope. This can be triggered by, for example, multiple rings or bubbles with nonuniform orientation. However, the exact details may require future study because the multiple-plume feature (Diehl et al. 2015) indicates the history of multiple shock locations when the carbon detonation is triggered.

Our explorations favor a high metallicity for the WD compared to the solar metallicity. In fact, this feature is common to the observed SNe Ia whose chemical abundance is extracted from their light curves and spectra. Such high metallicity appears to be common in recently observed SNe Ia, e.g., SN 2012cg (Graur et al. 2016; Leung & Nomoto 2018, 2020a) and SNR 3C 397 (Yamaguchi et al. 2015; Leung & Nomoto 2018). These works demonstrated that a supersolar metallicity is paramount to boost certain isotope or element ratios, especially the ⁵⁷Ni/⁵⁶Ni or Mn/Ni ratios. Such effects cannot be completely replaced by tuning other major parameters, such as the WD mass, M, or the nuclear runaway structure. These examples demonstrated that the metallicity of exploding WDs can be higher than solar metallicity. A detailed evolutionary path for such a high-metallicity WD progenitor would be interesting future work.

We note that there exist controversies regarding two-dimensional modeling containing symmetries, which might not necessarily be realized in reality. Ideally, three-dimensional models are required to provide a comprehensive and self-consistent explanation to match the explanation in a one-to-one correspondence. Here we briefly recapitulate how we use two-dimensional models and why this can still provide reliable estimates.

First, two-dimensional models allow a more time-effective search of appropriate models. As indicated in our previous works (Papers I and II), the parameter space that is suitable for SNe Ia is large. The running time for one n-dimensional hydrodynamical model scales with Nⁿ for a grid mesh of Nⁿ . Typical resolution requires ∼500 grids for one direction. This means that three-dimensional models are, at minimum, ∼500 times more computationally expensive. It generally takes 3–5 days for our two-dimensional model to complete its hydrodynamical simulation and nucleosynthesis. This simple scaling implies that a three-dimensional model requires, at minimum, months for a single model. This is beyond the computational time we can afford for a practical model investigation.

Second, the large-scale aspherical effect can be well captured by two-dimensional models. We recall that two-dimensional models are capable of naturally producing three-dimensional aspherical structure, such as a bubble or ring, as well as one-dimensional spherical structure. Three-dimensional simulations can produce more complex structures in the form of multiple bubbles, for example. In fact, the processes determining how the second detonation starts depend on the wave collision details. These can already be captured by one of the two-dimensional scenarios. Furthermore, two-dimensional models provide the minimum perturbation from the spherical symmetry. Notice that SN 2014J has features that can be explained by the classical W7 model in a broader picture, even though it has aspherical features. Therefore, it might be considered a generalization to start from models that behave almost spherically and extend them to lower symmetries.

Third, the symmetry is conserved in the simulation. As demonstrated in Moll & Woosley (2013), where one-, two-, and three-dimensional simulations of the sub-Chandrasekhar-mass models are carried out, one of their explicit three-dimensional models with a two-dimensional counterpart gives results that agree with each other. This provides some support that in the explosion phase, symmetry does not break during its propagation. Similar three-dimensional models in Gronow et al. (2020) also demonstrated similar features, in that the detonation propagates like a two-dimensional front. Furthermore, in this work, we further show that the two-dimensional spherical models explode spherically, as in the one-dimensional case. Both results support that our axisymmetric model remains axisymmetric throughout the simulations, as long as turbulent motion is unimportant. This is true during the explosion phase (∼1 s). When time is sufficiently long, i.e., in the nebular phase, we expect that the Rayleigh–Taylor instabilities play a role and perturb the morphology. Then, initial seeds break the rotation symmetry. However, such effects are secondary compared to the large-scale asymmetry and require more time to grow. Also, similar to models in the literature (see, e.g., Fink et al. 2007, 2010; Sim et al. 2010), we set up WDs in hydrostatic equilibrium as the initial condition. The quiet environment suggests that the turbulent motion is suppressed.

Finally, in this work, we focus on the common features that exist in both two- and three-dimensional models. In particular, we investigate how shock waves superpose, interact with each other, or grow by themselves through geometric convergence and consequently generate the structure that breaks the spherical symmetry. As indicated in Paper II, how the wave interacts is independent of the boundary condition. Thus, our two-dimensional model can offer the necessary starting point to explore which kind of detonation structure is necessary to generate the corresponding interaction for creating the large-scale asymmetry.

In this paper, we explored the parameter space in the classical double-detonation model that can produce observables indicating deviations from spherical symmetry. We studied how the initial detonation geometry affects the final explosion morphology by examining the ejecta composition in the spatial distribution and velocity space. We studied how the spherical symmetry can be broken for creating large-scale asymmetry. The sub-Chandrasekhar-mass WD progenitor tends to produce more pronounced asymmetry than the Chandrasekhar-mass WD progenitor. The surface He detonation can be the origin of the early ⁵⁶Ni gamma-ray line of some SNe Ia, e.g., SN 2014J, and the recently observed early bumps in the observed light curves of some SNe Ia.

We have examined how the initial mass and He-detonation geometry affect the final explosion results, in particular the ejecta geometry and element distribution in both spatial and velocity phase space. We observe that starting the He detonation as a bubble (with the lowest symmetry), then to a ring, and then a sphere (with the highest symmetry) may give observable differences in the ejecta morphology and velocity for the characteristic elements, including He, O, Si, S, Fe, and Ni, and their directional dependence.

We have provided a detailed case study on searching for models that may resemble the qualitative features observed in SN 2014J based on the gamma-ray line detections and the late-time photometry of the optical band. Four key aspects of SN Ia explosion are

(1) the total mass M of the WD determines the total ⁵⁶Ni production,

(2) the He-envelope mass M_He determines the large-scale asymmetry in the radial distribution of ⁵⁶Ni,

(3) the metallicity determines the required ⁵⁷Ni/⁵⁶Ni mass, and

(4) the initial He-runaway geometry determines the surface ⁵⁶Ni distribution.

From our explorations simulating a set of key scenarios, we conclude that the SN 2014J progenitor should have the following properties:

(1) an initial He detonation in the orbital plane set by the binary companion,

(2) a WD mass in the range from 1.00 to 1.10 M_⊙,

(3) a WD metallicity in the range from 3 to 5 Z_⊙, and 

(4) an He-envelope mass of ∼0.05–0.10 M_⊙.

We also derived the detailed velocity distributions of some major isotopes, for example ¹⁶O, ²⁸Si, ⁵⁴Fe, and ^56-58Ni, and the spatial distributions of major iron-peaked elements, including Cr, Mn, Fe, and Ni. Future observations of the ejecta morphology by specific elements (e.g., Seitenzahl et al. 2019) can provide a strong constraint on the models presented in this work. Large-scale features in these objects might reveal how the detonation has interacted during its propagation, thus shedding light on its initial detonation pattern.

At last, we discuss the recent application of SN tomography as presented in Seitenzahl et al. (2019) for the SNRs 0519-69.0, 0509-67.5, and N103B. From how aspherical the SN ejecta in the reverse-shock-heated region and their corresponding Fe masses are, we can deduce the fundamental properties of the progenitor, including whether it is a Chandrasekhar or sub-Chandrasekhar WD, the expected initial mass, and the detonation geometry. We also summarize the M and Z dependence of the major element ratios typically found in the spectra of SNRs.

When ¹²C burning is ignited in the center of the near-Chandrasekhar-mass WD, the central density is as high as 2−3 × 10⁹ g cm⁻³. Electrons are strongly degenerate, so the gas pressure is not so sensitive to the temperature. The temperature rise becomes unregulated, as the carbon-burning rate is strongly temperature-sensitive (∼T³³). This plants the first seed of nuclear runaway. Simultaneously, the rapid temperature rise does not trigger a shock because the pressure growth is small compared to the temperature rise due to the strong degeneracy, even when the central temperature becomes as high as  ∼10¹⁰ K (Nomoto et al. 1976; Nomoto & Sugimoto 1977). As a result, the temperature jump becomes a localized event, where the temperature discontinuity propagates by microscopic processes instead of macroscopic shock compression. Thus, the initial runaway is likely to be in the form of a deflagration wave (Nomoto et al. 1976, 1984). The short mean free path of electron conduction in such a density implies a very thin flame front (∼10⁻³ cm) compared to the size of a WD (∼10³ km; Timmes & Woosley 1992).

Despite the turbulent motion emerging down to the Kolmogorov scale (10⁻³ cm, assuming a typical Reynolds number of  ∼10¹⁴), the Gibson scale decreases with density and is generally larger than the Kolmogorov scale. Flame structure with a size below the Gibson scale (∼10⁻¹ km at the center) is smoothed (Niemeyer & Hillebrandt 1995; Röpke et al. 2003, 2004a), albeit eddies can appear below this scale.

If the propagation speed of the subsonic deflagration is fast enough, the deflagration efficiently releases the necessary energy for unbinding a WD and creates a successful explosion like the W7 model (Nomoto et al. 1984), as does the detonation wave (Arnett 1969).

The explosive nuclear burning at high densities synthesizes iron-peak elements (Thielemann et al. 1986; Iwamoto et al. 1999). However, the observed intermediate-mass elements should be synthesized at lower densities, which suggests that the explosion consists of subsonic burning, i.e., deflagration (Nomoto et al. 1976; Nomoto & Sugimoto 1977), which decreases the densities at the flame front.

If the propagation of the deflagration wave is slow, it may not unbind the star (but in some cases may cause pulsation) because the stellar expansion makes the deflagration wave quench (Nomoto et al. 1976). A subsequent transition from deflagration to detonation is vital for explaining a successful SN Ia explosion, known as the deflagration–detonation transition (Nomoto et al. 1984; Khokhlov 1991; Arnett & Livne 1994; Iwamoto et al. 1999; Gamezo et al. 2003, 2004; Röpke et al. 2007). However, the deflagration–detonation transition requires a turbulence strength that is less likely to be reached (see, e.g., Khokhlov et al. 1997; Niemeyer et al. 1999; Lisewski et al. 2000; Gamezo et al. 2005; Woosley et al. 2009). In order to realize this effect, ab initio numerical experiments with very fine resolutions (∼0.1 km) are necessary (Kushnir et al. 2012) that are 1–2 orders of magnitude below the affordable resolution. Recent direct experiments, both numerical and laboratory ones using a methane–air mixture to explicitly resolve turbulent motion, have demonstrated that turbulent acceleration can be an important key factor (Poludenko et al. 2011; Poludnenko et al. 2019).

The motion of the deflagration wave can be convoluted. The subsonic propagation (∼1% of the speed of sound) implies that the deflagration wave structure is coupled with the underlying fluid motion, which means that the flame structure is also susceptible to various hydrodynamics instabilities, such as Rayleigh–Taylor (Bell et al. 2004a; Zingale et al. 2005; Hicks 2015, 2019), Kelvin–Helmholtz, Landau–Derrieus (Bell et al. 2004b; Röpke et al. 2004b), and pulsational (Glazyrin 2013, 2014; Glazyrin et al. 2013; Poludnenko 2015) instabilities. On the contrary, the supersonic detonation is less sensitive to fluid motion. However, direct numerical simulations of the small-scale detonation show spontaneous cellular structure formation behind the detonation wave front (Gamezo et al. 1999).

To model the explosion, following how the deflagration propagates reveals how the energy is released. The subgrid scale of the reaction front indicates that on-site modeling is inaccurate, but a subgrid-scale model is necessary to describe partial cell burning and the irregular wave front inside the cell. This relies on the subgrid-scale turbulence model (Clement 1993; Niemeyer & Hillebrandt 1995; Shih et al. 1995a, 1995b) and the flame-tracking scheme.

The subgrid-scale turbulence model assumes that eddy motion below the resolved scale can be well described by statistical models. This gives an accurate approximation given the large difference between the resolved scale and the much smaller Kolmogorov scale. Scaling relation has been studied explicitly in direct simulations (see, e.g., Fisher et al. 2019). The model tracks the generation and dissipation of eddy motions by channels including shear stress, compression, Rayleigh–Taylor instabilities, and so on (see, e.g., Niemeyer & Hillebrandt 1995; Shih et al. 1995a, 1995b; Schmidt et al. 2006).

Flame-tracking schemes are algorithms designed for resolving subgrid-scale features. There are multiple representations, including (1) the advective–diffusive–reactive equations (Khokhlov 1995; Vladimirova et al. 2006; Townsley et al. 2007), (2) level-set methods (Osher & Sethian 1988; Sethian 2001), and (3) point-set methods (Glimm et al. 1999, 2000; Zhang 2009; Leung et al. 2015a). The main idea is to introduce additional variables with model parameters that represent how much the grid is partially burned; from that, the actual flame-front geometry is reconstructed.

In a sub-Chandrasekhar-mass WD, the less degenerate matter with a lower density means that the detonation is more likely (Nomoto 1982a). Such a scenario is viable when the surface energy production is faster than its heat loss by convection or expansion (Jacobs et al. 2016). The initial nuclear runaway can be triggered by accretion from its companion star in the single-degenerate scenario or through a violent merger in the double-degenerate scenario (see, e.g., Tanikawa et al. 2015, 2019). Pure CO matter has a high ignition threshold for its high temperature (∼1–2 × 10⁹ K; e.g., Sato et al. 2015), which is shown to be difficult to trigger and sensitive to the way of contact (Dan et al. 2012). A WD merger with a helium envelope can suppress this ignition condition, but its required hot-spot size can be unrealizable in a thin WD envelope (see, e.g., Shen & Bildsten 2009). Mixing with C/O-rich matter through, for example, turbulence can resolve this difficulty (Holcomb et al. 2013; Piro 2015). The violent merger of two CO WDs (Pakmor et al. 2012) is therefore challenging for a robust ignition, as the collision can fail to generate spots that are sufficiently hot for the first runaway (Dan et al. 2012, 2014). Mixture with helium provided by its He envelope or the companion star lowers the ignition temperature such that the detonation trigger is less sensitive to the merging dynamics (Shen & Moore 2014). A thin layer of He (∼0.01 M_⊙) can already trigger the second explosion more robustly (Pakmor et al. 2013). When the He detonation fails to trigger the second detonation, the star develops like a nova and explodes as a so-called SN Ia (Bildsten et al. 2007; Shen et al. 2010; Waldman et al. 2011).

Even when the C-detonation trigger becomes robust with the aid of He, the exact position and timing of the detonation are unclear because they depend on the dynamics of the He atmosphere. These require multidimensional low Mach number simulations of the atmosphere for multiple eddy turnover times. Multiple possibilities exist, including direct He ignition, C ignition, or ignition after the merger process when a Chandrasekhar-mass WD is formed (Dan et al. 2011; Shen et al. 2018; Tanikawa et al. 2019). Geometric convergence in a low-mass WD is more difficult to achieve (Shen & Moore 2014). High-resolution simulations using the ab initio approach are necessary to trace when and where the first hot spot appears (Fenn & Plewa 2017). The asymmetry in a three-dimensional simulation tends to suppress the prompt detonation as the geometrical convergence breaks down in the violent merger scenario (Fenn et al. 2016).

In a low-mass WD, the detonation front has a size comparable with the resolved scale. This allows direct coupling of the hydrodynamics with a nuclear reaction network (Shen et al. 2018; Polin et al. 2019). In a more massive CO WD (central density  ≥10⁸ g cm⁻³), the detonation width can be much smaller than the resolved grid size (∼10 km). Subgrid-scale methods or adaptive mesh refinement are often used in the literature.