Three-dimensional Core-collapse Supernova Simulations with 160 Isotopic Species Evolved to Shock Breakout

Both two-dimensional and three-dimensional simulations were run in spherical geometry. The spherical geometry is natural for self-gravitating objects and allowed us to easily "remove" the region of space containing the PNS. Following the studies of Wongwathanarat et al. (2015) and Stockinger et al. (2020), whose similar 3D simulations used parameterized models from the Prometheus-HOTB code (Scheck et al. 2006, 2008), we use an inner radial boundary of 500 km to excise the PNS. Because of the high sound speed and fine zoning in that region, the excision helps alleviate the CFL time step constraint. A point mass was placed at the origin to replace the mass of the excised PNS.

These simulations are intended to accurately capture the explosion throughout the entire star, approximately 10⁸ km in radius. An efficient way to accomplish this in spherical coordinates is to use a logarithmically spaced radial grid, as described in Fernández (2012) and shown in Wongwathanarat et al. (2015). This type of grid provides the ability to more easily maintain "square" zones with constant Δθ ≈ Δr/r and can more accurately track the near power-law density structure of stars. Though adaptive mesh refinement (AMR) provides an excellent way to resolve specific regions of the explosion, while efficiently ignoring others, FLASH's AMR is incompatible with log spacing. Consequently, we have implemented a log-spaced version of FLASH's uniform grid using logarithmically spaced blocks and uniformly spaced cells within each block along the radial dimension.

As outlined in Fernández (2012), we similarly define the domain between ${r}_{\min }$ and ${r}_{\max }$ such that consecutive block sizes have a ratio Δr_i+1/Δr_i = ζ > 1, where i is the block number, which increases with increasing radius. Logarithmic block spacing is achieved by setting

$$\begin{eqnarray}&&\zeta ={({r}_{\max }/{r}_{\min })}^{1/{N}_{r}},\end{eqnarray} \tag{ 1 }$$

where N_r is the number of radial blocks. The grid is then created over 0 ≤ q ≤ N_r by defining the inner edge of each block as

$$\begin{eqnarray}&&{r}_{q}={\zeta }^{q}\ {r}_{\min },\end{eqnarray} \tag{ 2 }$$

where ${r}_{0}={r}_{\min }$ and ${r}_{{N}_{r}}={r}_{\max }$. Each logarithmically spaced block contains 16 uniformly spaced cells in the radial direction.

The inner and outer radial grid boundaries are diode and outflow, respectively, the polar grid boundaries are reflecting, and the azimuthal boundaries are periodic. The diode boundary condition is similar to outflow but only allows matter to flow out of the domain, as opposed to letting matter freely enter the domain as well. The inner boundary is fixed until the first R-T instabilities begin to develop (∼2–3 s) and then is shifted to larger radii, following the progress of the shock. This is accomplished by removing the innermost radial block whenever the inner boundary becomes smaller than 1% of the minimum shock radius. This removes the region where the PNS, absent from our model, may have influence in the form of a PNS wind, and it makes the simulation computationally cheaper, progressively reducing the number of radial zones, which, in turn, relaxes the CFL time step constraint. The mass loss caused by moving the inner boundary is small but not negligible—it is consistently ∼10⁻⁵ M_☉ in all of our 3D simulations. As we will discuss in Section 3.3, this accounts for only ∼1.5% of the total mass lost, while the rest is due to fallback (matter passing through the inner boundary).

In addition, for all simulations, we define angle-averaged grid quantities as

$$\begin{eqnarray}&&\langle X(r)\rangle =\displaystyle \frac{\int X(r)\,d{\rm{\Omega }}}{\int d{\rm{\Omega }}},\end{eqnarray} \tag{ 3 }$$

where dΩ is the differential of the solid angle.

For each 3D model, the grid initially consists of 2304 × 192 × 384 total cells in r, θ, ϕ, respectively. The radial section of the grid extends logarithmically from 500 km ≤ r ≤ R_⋆, where R_⋆ is the stellar radius of the progenitor (see Table 1). As with Hammer et al. (2010), cones were excised along the poles to help further relax the CFL condition in 3D—in this case having a half-opening angle of 5°. The grid therefore covers 0.0278π ≤θ ≤ 0.972π at δ θ = 0.885° and azimuthal angles 0 ≤ ϕ ≤ 2π at δ ϕ = 0.938°. In two dimensions, the grid covers the same radial extent and polar angles 0 ≤ θ ≤ π with 2304 × 204 cells, respectively.

This leads to a radial resolution of nearly constant Δr/r of 5.7 × 10⁻³ and 6.1 × 10⁻³ for models D9.6 and D10, respectively. This can be compared to 6.9 × 10⁻³ of Müller et al. (2018), 8.9 × 10⁻³ of Stockinger et al. (2020), and 1 × 10⁻² for Hammer et al. (2010) and Wongwathanarat et al. (2015). All of our 3D runs have angular resolutions <1°, compared to 1° for Hammer et al. (2010), 1.6° for Müller et al. (2018), and 2° for both Wongwathanarat et al. (2015) and Stockinger et al. (2020).

The removal of the PNS from the grid allows us to neglect the high densities and temperatures present there and use FLASH's implementation of the Helmholtz equation of state (EOS; Timmes & Swesty 2000), which displays perfect thermodynamic consistency and includes contributions to internal energy from ions, electrons, positrons, and radiation. Because the Helmholtz EOS assumes full ionization, we halt each simulation when the shock front reaches the region of the progenitor where this criterion is no longer true (T ≲ 10,000 K), which happens only a few zones before shock breakout for our models.

FLASH allows us to track large numbers of species and utilize a multispecies network for nuclear burning, in this case the FLASH implementation of XNet (Hix & Thielemann 1999). Nuclear burning does not occur during these extended FLASH runs unless the Chimera runs serving as our initial conditions were stopped before nuclear burning was complete, which is the case with the D9.6–3d3d model, where XNet was on for the first 279 ms. Regardless, the inclusion of XNet gives us the unprecedented ability to track the composition of 160 nuclear species throughout the evolution of the explosion, which leads to a more accurate analysis of the ejecta seen at shock breakout. The species list of the sn160 network from XNet is given in Table 2. Of particular note are ⁶⁴Ge and ⁶⁶Zn, where proton-rich and neutron-rich flows that would progress to higher atomic number in nature stagnate in this network.

Each of our FLASH simulations are initialized from the final step of a 2D or 3D Chimera neutrino radiation hydrodynamics simulation of the supernova mechanism that has reached an asymptotic explosion energy. The Chimera simulations are initialized from two 1D progenitors calculated with the Kepler stellar evolution code (Weaver et al. 1978) up to the moment of Fe core collapse. These simulations are part of Chimera's "D-series" (and so prefixed) and are substantially similar in input physics to prior Chimera simulations except for the inclusion of the larger 160-species nuclear network to better handle the formation of neutron-rich ejecta. The explosion dynamics of the input Chimera simulations will be described in forthcoming publications. As the Chimera simulations do not include parts of the outer core or the envelope of the progenitor, the missing portions of the Kepler progenitor are reattached outside the limit of each Chimera simulation when mapping to FLASH. The Chimera simulations are ended when the computational intensity required to simulate the neutrino mechanism is no longer needed.

The first progenitor is a 9.6 M_☉ zero-metallicity star, provided by A. Heger (private communication). That choice is motivated by 3D simulations of this same progenitor by other groups (Melson et al. 2015b; Müller et al. 2019; Stockinger et al. 2020) and is used to explore progenitors relating to low-mass cores and to demonstrate the 160-species network. The diagnostic explosion energies at the end of the Chimera runs are 1.91 × 10⁵⁰ erg for Chimera model "D9.6-sn160-2D" and 1.68 × 10⁵⁰ erg for Chimera model "D9.6-sn160-3D," where the difference is due, in part, to the 2D run having been evolved nearly 200 ms further than the 3D model (E. J. Lentz et al., in preparation). We refer to these models collectively as D9.6 in this work. The explosion energy in 3D is ∼95% higher than the energy reported for the same progenitor in Stockinger et al. (2020) at the time of their mapping from Vertex-Prometheus to Prometheus-HOTB. Like low-mass oxygen-neon supernovae, the shock did not stall after bounce and the initial ejecta includes neutron-rich material drawn from the vicinity of the PNS that is not exposed to neutrino radiation in protracted pre-explosion convective heating. As a result, the outer ejecta of this explosion is enhanced with neutron-rich species like ⁴⁸Ca, ⁶⁰Ni, and ⁶⁶Zn that are not seen in typical iron core CCSNe and less of isotopes like ⁵⁶Ni and ⁴⁴Ti that are more common in supernovae that take longer to explode. This neutron-rich ejecta is noticeable at t_map in Figure 1 (left) and can be seen in the ⁶⁰Ni peak that rivals typical CCSN ejecta like ⁵⁶Ni.

Zoom In Zoom Out Reset image size

At the time it is mapped into FLASH, the D9.6 model has a relatively featureless density profile that gradually decreases ahead of the shock front until the edge of the star (Figure 2, left). The mean shock position at t_map resides in the He shell, which extends from 6.95 × 10³ km to 1.40 × 10⁷ km and accounts for 0.33 M_☉ of the total mass. An extensive H envelope spans from the edge of the He shell to the edge of the star at 1.50 × 10⁸ km and accounts for 7.85 M_☉ of the total mass. The compactness parameters, as described in O'Connor & Ott (2011), are ξ_2.5 = 7.65 × 10⁻⁵ and ξ_1.5 = 2.34 × 10⁻⁴, which are smaller overall compactness than our second progenitor.

Zoom In Zoom Out Reset image size

The second progenitor, a 10 M_☉ solar-metallicity star, was presented in Sukhbold et al. (2016) as a part of their study of 200 pre-supernova models. The 2D Chimera model ("D10-sn160-SEWBJ16"; J. A. Harris et al., in preparation) has a diagnostic explosion energy of 3.075 × 10⁵⁰ erg, which is almost double the energy of the D9.6 3D at its t_map and ∼50% higher than the D9.6 2D energy.

The D10 Chimera model is a traditional CCSN model with the shock stalling shortly after bounce and significant accretion onto the PNS occurring. This explains the lack of neutron-rich material when comparing composition profiles of the two models in Figure 1—note the lack of a ⁶⁰Ni peak. Combined with a significantly higher presence of ¹²C and ¹⁶O, this leads to a different profile in the ejecta lying behind the shock. More fluctuations are also noticeable within the composition profile owing to the more prolate shock front in the D10 compared to the relatively spherical shock front of the D9.6, leading to a less uniform angular distribution of ejecta.

In further contrast to the D9.6, this progenitor has a rather erratic density profile, especially noticeable in ρ r³, with a dramatic change in density gradient at the He/H interface (Figure 2, right). The mean shock position at t_map resides in the former He-burning shell, which extends from 2.02 × 10⁴ km to 3.20 × 10⁵ km. This shell is a key feature in the ρ r³ profile, for it is the source of a dramatic acceleration that the shock experiences when transitioning to the inert He layer residing above. The He-burning shell contributes 0.44 M_☉, while the remaining He layer, which ends at 4.32 × 10⁶ km, provides a comparable 0.43 M_☉ for a total mass of 0.87 M_☉ of the entire He shell. The similar masses for each section of the He shell spread across widely different spatial extents explain the change in density gradient at the transition point. A hydrogen envelope spans from the edge of the He shell to the edge of the star at 3.57 × 10⁸ km and accounts for 7.2 M_☉ of the total mass. An additional density feature can be seen near the edge of the hydrogen shell located at 1.49 × 10⁸ km. The compactness parameters are ξ_2.5 = 2.04 × 10⁻⁴ and ξ_1.5 = 4.32 × 10⁻¹. The large difference between ξ_2.5 and ξ_1.5 is the result of the (C+O)/He interface lying at 1.61 M_☉. Details of key interfaces for both progenitors, as well as mapping times from Chimera to FLASH, are given in Table 1.

Model D9.6–3d3d was mapped into FLASH at a time t_map = 466.6 ms after the bounce that marks the formation of the PNS, having been simulated to that point with Chimera. At this point in the explosion, the mean shock radius is at ∼1.0 × 10⁴ km, just across the (C+O)/He interface. As noted by Stockinger et al. (2020), this progenitor is in the process of a second dredge-up of the He shell that has created a section at the base of the shell that contains minimal hydrogen (in contrast to the rest of the He shell). The shock encountering changes to ρ r³ in this region explains the slight deviation in the trend of the shock velocity at ∼1.7 × 10⁴ km and ∼6.0 × 10⁴ km (see Figure 3).

Zoom In Zoom Out Reset image size

At the time of mapping from Chimera to FLASH, the metal-rich shell lying behind the shock is mainly composed of ¹²C and ¹⁶O. It is this shell, which is quasi-spherical, that begins to deform and starts to develop the initial R-T instabilities. This shell, located at the green-to-yellow transition in Figure 4(a) at ∼8.5 × 10³ km, is also the location of the reverse shock created upon crossing the (C+O)/He interface. The departure from sphericity is imprinted on the reverse shock at its creation. Although the location of the mass shell is the position of the reverse shock in this scenario, this is not always the case as we will see with the He/H interface discussed below. Because the main shock has only just crossed the (C+O)/He interface, there are still portions of it that are still traveling down the density cliff; thus, overall the shock is still accelerating at this point—represented by the velocity spike shown in Figure 3. There is also a wind termination shock (also in Stockinger et al. 2020) that resides close to the inner boundary and will eventually collapse inward owing to the absence of the PNS and its wind from the simulation. From this point forward, the explosion propagates through the He core, and the deformed metal-rich shell starts to mix with that material.

Zoom In Zoom Out Reset image size

By ∼2 s, the shock has crossed fully into the He layer with the initial R-T plumes appearing as ripples in the fragmenting metal-rich shell. The instabilities begin to develop their typical mushroom state at 10 s and are still mainly composed of the species from the metal-rich shell (see ripples at ∼2.3 × 10⁵ km in Figures 4(b), 5(b)). Beginning at approximately 30 s, the inner regions of the ejecta (the "hot bubble") are injected into the rear of the instabilities, including the key isotopes of nickel like ⁵⁶Ni and ⁶⁰Ni. To track the bullets, we have combined the mass fractions of ⁵⁶Ni and neutron-rich iron group nuclei (${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$) and have taken a 3% isosurface of the result, which enables a direct comparison to the tracking of ${X}_{{}^{56}\mathrm{Ni}+\mathrm{Tr}}$ bullets in Stockinger et al. (2020). To approximate the crude tracer nucleus of Stockinger et al. (2020), we define "neutron-rich iron group" as all species in our network falling in the range of ⁴⁹Cr–⁶⁴Ni, excluding ⁵²Fe and ⁵⁶Ni.

Zoom In Zoom Out Reset image size

Large-scale features start to form at ∼60 s, where the radial shock position is ∼1.0 × 10⁶ km, which can be seen in Figure 3 as the ${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$ maximum velocity (dashed blue) curve crosses the shock's (solid blue) curve. Fluctuations in the velocity after the crossing point are due to plume interactions with the shock. As the fastest-moving bullet penetrates the shock, that bullet slows, and the maximum velocity shifts to the next-fastest bullet. This continues until all of the fast-moving clumps eventually interact with the shock front, which then results in a steady decline of the maximum velocity of these R-T plumes. These features can explicitly be seen penetrating the shock in Figure 4(d).

By 150 s, there is no semblance left of the metal-rich shell, as the inner ejecta from the hot bubble has completely engulfed it. The R-T fingers have grown significantly by this point and have reached the back of the shock (see large mushroom features in the middle panel of Figure 6). As the shock continues to progress through the He core, the R-T fingers progress with it, remaining near the rear of the shock (see elongated fingers penetrating the shock in Figures 4(e) and 5(e)). Whether the R-T fingers penetrate the shock is key to the morphology of the remnant. The shock experiences a gradual deceleration in this region of the progenitor due to the increasing ρ r³ and the extent of the He shell. Additionally, the reverse shock created at the first density interface has continued to propagate inward in mass and starts to shred the inner regions (the blue region at ∼2.5 × 10⁶ km in Figure 4(e)). This reverse shock is not spherical as a consequence of the asphericity and timing of the main shock's interactions with the prior composition interface.

Zoom In Zoom Out Reset image size

At ∼1000 s (Figures 4(f) and 5(f)), the shock crosses the He/H composition interface located at 1.4 × 10⁷ km and creates a weak pressure wave due to the minimal change in ρ r³ (see Figure 3) that propagates inward in mass and radius before eventually steepening into a second reverse shock. This delay ensures that the second reverse shock location is quickly decoupled from position of the mass shell at the He/H interface, whereas the (C+O)/He mass shell and first reverse shock positions coincided. Although slight, the deceleration gives the closest R-T plumes to the shock front the opportunity to interact with the rear of the shock. We only see this interaction happen in D9.6–3d3d and D9.6–2D3D _Tilted, which end up having higher overall velocities compared to the other simulations (discussed in Sections 3.2 and 3.3). This interaction not only seeds new instabilities but also further develops the most dominant R-T fingers into even larger mushroom-shaped plumes that are able to penetrate and reshape the shock. These will be the fastest bullets at shock breakout, though the shock must still propagate through most of the envelope before they reach that point.

By 2500 s, the inner regions of the explosion are completely shredded by the first reverse shock. The mixing effects can be seen when transitioning from Figure 5(f) to Figure 5(g) (note that the distribution of ejecta is much less uniform in the inner regions after the transition). Having been born quasi-spherical and propagated through inhomogeneous regions of the star leaves the reverse shock aspherical. As a result, although some portions of the collapsing reverse shock pass through the inner radial boundary (located at r = 1.8 × 10⁵ km at this point in the simulation), most of the reverse shock bypasses the boundary altogether and collides with itself off-center rather than at the PNS. This sets up the creation of a secondary forward-propagating blast wave that is reminiscent of the implosions discussed in supernova remnant theory (Truelove & McKee 1999; Cioffi et al. 1988). The blast wave can be seen at ∼5.0 × 10⁶ km in Figure 4(g). This causes significantly more mixing, as the inner regions also bounce off the reflecting boundaries of the grid, and the R-T plumes grow to be quite abundant.

What once were primarily metal-rich mushrooms are now heavily coated in helium, for the propagation through the He core has filled the gaps between the R-T fingers and has shaped them further. However, as noted earlier, the original inner regions of the ejecta still form the "bulk" of the inner anatomy of a single finger (see Figure 7) owing to the injection through the metal-rich shell. Most notably represented in the main anatomy of an instability are the Ni isotopes, as expected, with the most abundant isotope occupying the bullets being ⁶⁰Ni, from the early, neutron-rich portion of the hot bubble.

Zoom In Zoom Out Reset image size

As the shock continues to expand, the pressure wave created at the He/H interface reaches the center of the grid at ∼20,000 s, after steepening into a reverse shock at ∼16,000 s. As the first reverse shock collided with itself, so does the second, but it does not rebound as hard as the former, because this second reverse shock has been weak since its launch as a result of the slight deceleration of the main shock noted above. Nevertheless, reaching the center still creates another forward-propagating blast wave that further influences the inner ejecta.

The shock continues to propagate through the H envelope until hitting the edge of the star, and grid (1.5 × 10⁸ km), at ∼70,000 s (19.4 hr). We "rewind" and declare the end of our simulations at ∼62,000 s, as this is the time where the shock enters the region of the progenitor where the Helmholtz EOS's assumption of fully ionized hydrogen is no longer valid (T ≲ 10,000 K). As the shock has only been backtracked to ∼1.4 × 10⁸ km from 1.5 × 10⁸ km, the changes to ejecta morphology, yields, and speeds are negligible. Carrying the models further would require accounting for the circumstellar environment and radiation hydrodynamics of shock breakout.

The synchronous conversion between kinetic and internal energy through the entire evolution of the explosion can be seen in Figure 8. Both quantities respond to changes in v_Shock as the shock front moves through the density structure of the star. At the end of our simulations, the internal energy is still in the process of converting to kinetic energy, which starts to converge toward the total energy of the system.

Zoom In Zoom Out Reset image size

Though the D9.6–3d3d model clearly has significant extended plume features present at shock breakout (see Figures 6 and 7), the majority of trailing R-T bullets have only made it to approximately 1.0 × 10⁸ km and are therefore well short of the surface of the star. The early development of features within the He shell, in combination with a smaller relative velocity gap between the shock and fastest-moving Ni bullets, enables the further spawning of large R-T mushrooms at the He/H interface. This model demonstrates how the early-time asymmetries can impact the late-time evolution of a CCSN. Asymmetries at the time t_map seed the initial instabilities from the (C+O)/He interface that spawn further R-T plumes upon reaching the He/H interface—provided that they are moving fast enough relative to the shock. This, in turn, affects the efficiency of radial mixing in the outer envelope. As seen in the top panels of Figure 9, the instabilities spawned at the He/H interface in this model drive the mixing of metal-rich ejecta beyond the inner 4 M_☉ to the edge of the star. The bulk of the bullets end up with a peak centered around 500 km s⁻¹ with the yields for ¹⁶O, ²⁸Si, ⁵⁶Ni, and ⁶⁰Ni in that region all in the range between 1 × 10⁻³ M_☉ and 1 × 10⁻⁴ M_☉. However, the extent of radial mixing is quite apparent, with a high-velocity tail reaching to ∼1750 km s⁻¹, including yields between 1 × 10⁻⁵ and 1 × 10⁻⁷ M_☉ for these same species.

Zoom In Zoom Out Reset image size

These bullets are heavily coated in ⁴He, with the maximum bins of Figure 9 exceeding 1 × 10⁻¹ M_☉ around 1000 km s⁻¹ and roughly 1 × 10⁻³ M_☉ around 1500 km s⁻¹. The most unusual aspect is the internal anatomy of these metal-rich clumps. The typical isotope associated with these types of bullets in the literature has consistently been ⁵⁶Ni; however, ⁶⁰Ni seems to fill that role in this star. We suspect that this is due to the nature of very low mass CCSNe, but additional simulations with realistic networks are required (to be discussed in E. J. Lentz et al., in preparation). Across mass and velocity spaces, ⁶⁰Ni is the most abundant of our isotopes in the iron group, and it occupies more of the large-scale features, whereas the ⁵⁶Ni resides more in the microstructure (Figure 7, right). Although surprising, the distribution of these isotopes at the time of t_map from Chimera (Figure 1, left) makes this the most logical outcome. The explosion is surrounded by a shell of ¹²C, ¹⁶O, and ²⁸Si, but the two relevant Ni isotopes are distributed in such a way that the ⁵⁶Ni occupies the innermost ejecta, whereas the ⁶⁰Ni is more extended (Figure 7, left). The extended ⁶⁰Ni features present at the time of t_map grow into further extended structures as the explosion progresses, thus mixing more effectively in mass and velocity in the explosion.

Figure 10 provides a more detailed look at the compositional structure of the bullets. In the left panel of Figure 10, we plot the angular distribution of the composition of R-T plumes residing in the x-y plane (θ = 90°) at a constant radius of 4.0 × 10⁷ km—marked as the dashed circle slicing the plumes at this radius in Figure 4(h). The right panel of Figure 10 displays the composition versus radius of a specific plume residing at ϕ = 18° in the x-y plane. We choose this time (10,000 s) for this inspection because the plumes are much more distinct in Figure 4(h) compared to their later appearance and the R-T plumes are simply expanding beyond this point in their evolution. In this plane, there exists six extended R-T plumes residing at 18°, 43°, 86°, 115°, 155°, and 171°. Two additional, less extended R-T plumes can also be seen at 9° and 350°. All of the extended bullets are consistent in composition, having a ⁴He coating that surrounds a metal-rich interior dominated by ⁵⁶Ni+IG. We find that in addition to the significant amount of ⁵⁶Ni+IG present in a single bullet (X ≃ 0.07), there also exists a substantial amount of ¹⁶O (X ≃ 0.02) and lesser amounts of ⁶⁶Zn (X ≃ 0.008) and ²⁸Si present (X ≃ 0.001). The relatively large presence of ⁶⁶Zn in this model is representative of the enhanced α-rich, neutron-rich ejecta seen at the end of the Chimera run. (⁶⁶Zn is the neutron-rich upper limit of the sn160 network.) A closer inspection of the 18° bullet can be seen in the inset in the left panel of Figure 10. The vertical dashed lines indicate the angular positions of the edge of the 40% ⁴He coating, and the intersection with the horizontal dashed line indicates what value this represents in ${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$ (0.001 or 0.1%), demonstrating the correspondence between these isosurfaces. We will track the cocoon of ⁴He that encases the heavy-element bullets in Section 3.3 by creating an isosurface at 0.1% ${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$ instead of 40% ⁴He, which allows us to analyze the evolution of the external coating without having additional noise at early times from the He shell.

Zoom In Zoom Out Reset image size

The D9.6–2D3D model exhibits similar behavior when it comes to the general progression of the shock front, yet differences can be seen when analyzing the leading R-T bullets. Although the Ni bullets are able to catch up to the rear of the shock, they are never able to fully interact with it in this model owing to a sufficiently large gap in the relative velocity between v_Shock and v_bullets. This can be seen explicitly in Figure 3 in the shock and ${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$ velocity curves (green lines). The bullets are closest to the shock when the shock front hits the He/H interface at ∼1.0 × 10⁷ km and ∼1000 s. The maximum velocity of the bullets in this model never rises above the average shock velocity in Figure 3, explaining why the plumes only minimally interact with the shock. Furthermore, in contrast to the D9.6–3d3d and D9.6–2D3D _Tilted models, the maximum radial position of the ${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$ isosurface (Figure 11, upper edge of the green shaded region) always stays just below the curve representing the average position of the shock (green solid line), which highlights the absence of extended features. This minimal interaction leads to the scarcity of large-scale structures and asymmetries in D9.6–2D3D. Despite that D9.6–2D3D is mapped roughly 200 ms later than D9.6–3d3d and the explosion energy is ∼13% higher in D9.6–2D3D, the enhanced growth rate of the R-T plumes enabled by the 3D initial state allows D9.6–3d3d to retain higher-velocity bullets.

Zoom In Zoom Out Reset image size

The lack of macrostructure is apparent in the yields of key isotopes at shock breakout (Figure 9). The velocity distribution of the ejecta extends only to ∼1225 km s⁻¹ in D9.6–2D3D, much less than the typical velocities associated with SN 1987A, and ∼30% lower in maximum velocity than D9.6–3d3d (1750 km s⁻¹), where the plumes interact with the shock. Because the shock is plowing through ⁴He and ¹H, the shock can be seen as the "hump" in the ⁴He curve centered at ∼1000 km s⁻¹, whereas the bulk of the bullets can be seen as the metal-rich hump further behind, peaking at ∼500 km s⁻¹. The gap between the humps shows how large the relative velocity between the shock front and metal-rich clumps is in D9.6–2D3D. The distribution gap in velocity further explains the inefficiency of mixing in mass space as well, with most of the metal-rich ejecta only extending to just within 4 M_☉. Large-scale mushrooms are never spawned from the interaction of the shock with the He/H interface; thus, the metal-rich ejecta stays trapped behind the wall of He and is unable to extend its radial mixing.

The shock in D9.6–2D3D remains roughly spherical, even late in the evolution of the supernova. The D9.6–2D3D model is nearly identical to our D9.6–2D2D simulations, and this run can be viewed as, in essence, a 2D simulation existing in 3D. Without transverse velocities in the 2D Chimera model, due to the initially assumed axisymmetry, true 3D behavior never develops (note the unbroken axisymmetry in Figure 12, left panel). However, the absence of structure in the shape of the shock is more than made up for by the amount of microstructure present in the inner regions containing the bulk of the R-T instabilities. This, in general, is similar to the results of Kifonidis et al. (2003) and Wongwathanarat et al. (2015), though with different progenitors. From a yields perspective, in both mass and velocity spaces, D9.6–2D3D is the most similar to the "3d3d" model presented in Stockinger et al. (2020), which uses the same progenitor, referred to as "z9.6" (discussed further in Section 3.4).

Zoom In Zoom Out Reset image size

The failure of the D9.6–2D3D model to depart significantly from its 2D origin is a concrete reminder that the development of fluid flows in response to instabilities is determined by both the strength of the instabilities and the perturbations that seed the fluid flows. In the D9.6–2D3D case, even though mapping to 3D allowed instabilities in the longitudinal direction to develop, the much stronger latitudinal velocity perturbations result in the growth of instabilities that maintain the axisymmetric character of the model. To test this assertion, we constructed the D9.6–2D3D _Tilted model based on a simple coordinate transform. Tilting the original 2D Chimera model by 90° on its axis through a coordinate transform introduces longitudinal velocities similar in scale to the previously 100% latitudinal velocity field. Other choices of rotation angle would presumably have similar effects, as long as the misalignment of the old and new axes produces significant longitudinal velocities. The presence of both longitudinal and latitudinal velocities seeds the development of features in both coordinate directions. Similar to D9.6–3d3d, the D9.6–2D3D _Tilted model establishes extended features in its explosion, which allows us to further investigate the morphology of this system. Features develop in this version of the "2D" model owing to the evolution of a spherical-bubble structure rather than a pure toroidal structure imposed by axisymmetry. Hence, rotation of the angular velocities enables the explosion to deviate from the initial toroidal structure and start developing bubble-type features when forming the R-T bullets. Although only demonstrating a slight deviation initially, the bubble structures are able to retain higher velocities owing to experiencing a lower drag-to-buoyant-force ratio and deviate further from axisymmetry as the explosion progresses. Echoes of the axisymmetric origin persist, but they do not dominate the entire model to the extent seen in the D9.6–2D3D model. Side-by-side comparisons of all D9.6 models shown in Figure 12 demonstrate the true impact that tilting the model has on the resulting morphology of the explosion. Clearly, the D9.6–2D3D _Tilted model, though retaining a grossly axisymmetric form, looks more like its true 3D counterpart, while the D9.6–2D3D model remains almost purely toroidal.

The similarities between D9.6–2D3D _Tilted and D9.6–3d3d are also apparent in the distribution of radial velocity across the entire grid (Figure 13). Unlike D9.6–2D3D, both D9.6–2D3D _Tilted and D9.6–3d3d have a significant number of grid cells occupying high-velocity space beyond a radius of 7.0 × 10⁷ km. Additionally, the overall shape of the D9.6–2D3D _Tilted distribution at larger radii looks similar to the D9.6–3d3d distribution, with a high peak of grid cells before the shock front resulting from the extended features produced in those models. Although obscured at lower radii by the D9.6–2D3D data, D9.6–2D3D _Tilted and D9.6–3d3d are still consistent, where D9.6–2D3D is an outlier.

Zoom In Zoom Out Reset image size

Analyzing the 3D surface area of the ${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$ isosurface shown in Figure 14 further illustrates the divergence between the D9.6–2D3D _Tilted (orange lines) and D9.6–2D3D (green lines) models. The surface area representing the inner anatomy of the bullets (the 3% isosurface) is nearly identical for D9.6–3d3d and D9.6–2D3D _Tilted, while D9.6–2D3D quickly falls behind in the development of surface area. The divergence starts when the shock front encounters the He/H interface, because, once encountering this region, the D9.6–2D3D model does not have extended features penetrating the interface, which would significantly contribute to the surface area, while the other two models do.

Zoom In Zoom Out Reset image size

The external coating of the bullets (represented by the 0.1% isosurface) is visualized in Figure 12 and displayed more quantitatively by the respective surface area curves in Figure 14 (dashed lines). Despite the formation of extended structures in D9.6–2D3D _Tilted, those features do not occupy as much overall surface area as in the D9.6–3d3d model, which can be seen in the isosurface plot. The biggest plumes in D9.6–2D3D _Tilted do not grow as large as the biggest plumes in D9.6–3d3d, which affects the rate of change of the surface area. The largest contribution to the surface area occurs at the peak in the bottom panel of Figure 14, which represents the time that the bullets and shock hit the He/H interface. The addition of the He/H mass shell to the coating of the bullets provides this boost owing to the significant amount of extra volume it adds to the bullets. As the bullets expand, their surface area grows and reaches a point near 10,000 s where the slopes converge toward the contribution provided by shock expansion. The D9.6–2D3D model converges much faster, as it has no large plumes contributing to its evolution, whereas the other two models are able to stay above the shock expansion curve for longer. Although overall converging toward the contributions from the shock, the D9.6–2D3D _Tilted model is able to achieve shock breakout while the total surface area resides above the curve represented by the shock (top panel of Figure 14). This does not occur in D9.6–2D3D, as the total surface area is dominated by the shock starting at ∼40,000 s while almost all of the bullets are trapped behind the He/H mass shell and constantly outpaced by the shock.

Although the external coating isosurface tracks the larger structures of each model, its surface area contribution is not as large as the 3% isosurface. This is due to the fact that while the 0.1% tracks larger structures that produce overall greater individual contributions, the amount of smaller individual contributions from the 3% isosurface is more numerous and builds up to occupy more of the volume, thus representing a larger total surface area.

The overall radial progression of the shock for all three models is nearly identical whether viewed as average shock radius relative to velocity (Figure 3) or time (Figure 11). Similar to the D9.6–3d3d model, D9.6–2D3D _Tilted develops its larger-scale features in the middle of the He shell near 1.0 × 10⁶ km. This can be seen explicitly in Figure 3, as the maximum velocity of the bullets surpasses the average shock velocity. The same types of variations in the velocity profile after 1.0 × 10⁶ km occur in the D9.6–2D3D _Tilted model as they did in D9.6–3d3d before the velocities of the bullets steadily decline until shock breakout is achieved. The plume penetration into the shock is seen further in Figure 11, as the bullets' maximum radial extent in the D9.6–2D3D _Tilted model (orange shaded region) surpasses the average shock radius and even exceeds that of the equivalent highlighted range in the D9.6–3d3d model. The average radius of the ${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$ clumps in the D9.6–2D3D _Tilted model is nearly identical to that of D9.6–3d3d, while the D9.6–2D3D model deviates around 60 s—the time when D9.6–2D3D _Tilted and D9.6–3d3d form their large-scale structures.

Due to the formation of extended features in D9.6–2D3D _Tilted, this model bridges the gap between the D9.6–2D3D and D9.6–3d3d yields of metal-rich ejecta in both mass and velocity spaces (see Figure 9, bottom panels). Not only is the extent of radial mixing similar to that of D9.6–3d3d, but the maximum velocity of the high-velocity tail is 200 km s⁻¹ larger than the D9.6–3d3d model (1950 km s⁻¹). Once again, the bulk of metal-rich bullets peak at ∼500 km s⁻¹ and 1.0 × 10⁻⁴ M_☉. The dominant isotope of the iron group in D9.6–2D3D _Tilted is ⁶⁰Ni, matching the D9.6–3d3d model. The total yields are relatively consistent across all models (Table 3), with the largest differences arising owing to evolution within Chimera for the 2D initial condition, while also having a different t_map than the 3D initial condition. Additionally, we see relatively low mass loss across all models, ∼7 × 10⁻⁴ M_☉ lost, where ∼1.5% of this is due to the moving inner boundary (removal of innermost grid cells) and the remaining ejecta lost is due to fallback (matter falling through the inner boundary), most of which is ⁴He (∼3 × 10⁻⁴ M_☉).

Previous works have also studied supernovae from the same progenitor as our D9.6 models, but only one has studied the long-time evolution. Stockinger et al. (2020) started from the Melson et al. (2015b) simulation of the neutrino heating phase using similar microphysics in 3D in their version of the 9.6 M_☉ progenitor, but with a smaller nuclear network (15 species + n + p) and lower resolution than our initial state (see Section 2.2). They report metal-rich clumps centered around ∼300 km s⁻¹ and extending to a maximum of ∼500 km s⁻¹ in velocity space, which is clearly slower than our high-velocity tails extending to 1225, 1750, and 1950 km s⁻¹ for D9.6–2D3D, D9.6–3d3d, and D9.6–2D3D _Tilted, respectively. The lower clump velocities in their run also lead to less efficient radial mixing, with the metal-rich ejecta only falling within the inner 2 M_☉. These results are starkly different from the results of our respective D9.6–3d3d model, which shows mixing to the surface, and are less well mixed than even our simulations initiated from 2D Chimera models.

We believe that the notable differences in Stockinger et al. (2020) derive from the lower overall diagnostic explosion energy reported at their t_map, for our initial energy is ∼95% larger in the 3D3D case and even larger in the 2D3D cases. Consequently, our shock achieves breakout nearly 12 hr sooner than their reported breakout of ∼31 hr, which is an approximately 60% difference compared to our shock escape before rewinding (∼19.4 hr). The weaker overall explosion helps explain why their model does not produce large structures during its evolution, despite it being a 3D model starting with 3D initial conditions. The relative velocity gap is too large between the shock front and the leading metal-rich bullets, which enables the R-T plumes to get trapped behind the He/H mass shell as opposed to spawning large features from it (similar to our D9.6–2D3D model, but to a greater extent). This is seen explicitly in Figure 13 from Stockinger et al. (2020) (equivalent to our Figure 3), where their isosurface of ${X}_{{}^{56}\mathrm{Ni}+\mathrm{Tr}}$ never reaches maximum velocities that are larger than their average shock velocity. Further comparison can be seen in Figure 20 from Stockinger et al. (2020) (equivalent to our Figure 6), as our bullets look distinctly more elongated while propagating through the He shell. The total yields at shock breakout for D9.6–3d3d are shown in Table 3. They differ slightly from the yields at the end of the Chimera D9.6–3D model owing to the 279 ms of additional burning that occurred. The largest relative changes are in trace abundances of CNO-cycle products (¹⁴O, ¹⁷F, ¹⁸Ne). More important are enhancements as large as factors of several in species formed during proton-rich, α-rich freezeout (e.g., ^53,54Co, ^57,59Cu), including the precursors to the ν p-process (^63,64Ga, ⁶⁴Ge). Similar enhancements occur during the same interval in the Chimera D9.6–2D model. The yields are also broadly comparable to those listed in Stockinger et al. (2020). However, differences can be seen in the form of ²⁸Si, ⁵⁶Ni, and the iron group tracer material. The amount of ²⁸Si present at the end of our simulation is 128% greater than that reported by Stockinger et al. (2020). They report approximately 68% more ⁵⁶Ni than our total, which may result from the inclusion of the mass of all iron group species not included in their network increasing the ⁵⁶Ni yield, as they discuss. Overall, the amount of neutron-rich iron group material across all of our models is larger by an order of magnitude (∼900%). Considering the amount of ⁶⁰Ni present in our model, and that it has essentially "replaced" ⁵⁶Ni as the traditional bullet material in this simulation, we see comparable or greater total nickel and iron group yields.

Neutron-rich iron peak isotopes are an area where the sn160 network we employed has significant advantage over the smaller network of Stockinger et al. (2020), even with their tracer species. Since we see little mass loss of the ejecta during our extended FLASH runs and no significant formation of iron group nuclei during our short period of nuclear burning, the main cause of the discrepancy in the yields seen at shock breakout between D9.6–3d3d and Stockinger et al. (2020) must be how the species were evolved in the Chimera and Vertex-Prometheus portions of the runs. As Stockinger et al. (2020) also initiated their shock breakout run from an early-time CCSN simulation (Melson et al. 2015b), and because there is no discussion of notable mass loss during their late-time evolution, we stress how critical the initial conditions are in this yields comparison.

To a lesser extent, we believe that our higher resolution also impacts the morphology of the system. As we will discuss further in Section 4.1, resolution directly affects the number of R-T plumes spawned when a mass shell fragments. The fragmenting phenomenon determines how the shock front is able to be reshaped by the metal-rich bullets. A greater number of extant R-T plumes allow for more shock interaction across the entire domain, directly affecting the development of large-scale features. However, a more extensive resolution study is required to support this supposition.

For a more general comparison of our D9.6–3d3d model, we look to the morphology analysis in Wongwathanarat et al. (2015), who categorized the late-time metal-rich ejecta into three types: (1) small clusters of R-T bullets having the fingerprint of early-time asymmetries as in their 15 M_☉ red supergiants (RSGs), (2) fragmented and squished round features as in their 20 M_☉ blue supergiant (BSG), and (3) long extended fingers as in their two 15 M_☉ BSGs. Our RSG D9.6 simulations do not seem to fall completely into one of these regimes, but the reasoning outlined by Wongwathanarat et al. (2015) does explain why our models look the way they do. In their 15 M_☉ BSG models, the steep rise of ρ r³ inside the He layer causes a steady deceleration of the shock front as it propagates, and the acceleration/deceleration at the He/H interface is nearly nonexistent. This allows the bullets to stay close behind the shock and avoid interaction with any of the reverse shocks. This is the type of density profile found in D9.6–3d3d, as the metal-rich bullets are able to catch up to the rear of the shock in the middle of the He layer owing to higher maximum velocities than the average shock velocity.

Despite the presence of extended features in our simulation, they are not as extreme and distinct as those produced by the 15 M_☉ BSG simulations of Wongwathanarat et al. (2015). Our features look like slightly more extended versions of their clustered RSG fingers. These clustered structures are associated with early-time asymmetries and are clearly visible in Figure 7; however, the journey for our clumps is different. In Wongwathanarat et al. (2015), the large gap between the shock and the trailing bullets allows for more momentum to build before they collide with the reverse shock produced by the dramatic deceleration at the He/H interface, which does not occur in our simulations owing to the smoother density profile of the D9.6 progenitor. The development of the bullet shape is strongly impacted by the interaction of the fingers with reverse shocks, which squash the clumps. Because the majority of our R-T plumes spawned ahead of the first reverse shock and out of the first mass shell, they completely avoid any reverse shock interaction. (The second reverse shock also forms behind the bullets).

As important as the dynamics of the reverse shock are, neither Stockinger et al. (2020) nor Wongwathanarat et al. (2015) discuss the phenomenon of the first reverse shock setting up a point-like rebound blast wave as it approaches the inner boundary as seen in our D9.6 models. We suspect that this event is missing owing to how they moved their inner boundary. As discussed in Section 2.1, we mimicked Wongwathanarat et al. (2015) in the handling of our inner boundary of the grid but used a 1% of shock radius criterion as opposed to their 2%. This means that we waited longer to move our boundary, thus allowing more accurate interactions near the center of the grid. If the collapsing reverse shock encounters the inner boundary when its radial excision is too large, then the reverse shock does not have the opportunity to set up a point-like blast and instead exits the grid. Regardless, the impact of this event on the morphology of the system, as well as its interaction with the PNS wind, needs to be explored further.

Model D10–2D3D was mapped into FLASH at a much later t_map (∼1.76 s), as nuclear burning ceased much later in this explosion compared to the D9.6 models. In addition, the shock front in the D10 is significantly more aspherical from the start (Figures 15(a) and 16(a)) compared to the D9.6 (Figures 4(a) and 5(a)). This is a common feature in 2D models of iron CCSNe, because the development of the explosion depends on the development of large plumes that can deliver energy from the regions of most intense neutrino heating to the shock. The three big convective plumes have already shaped the shock front at t_map and will eventually create shear Kelvin–Helmholtz (K-H) instabilities at the interface between the dominant plumes later in its evolution. At t_map, the mean shock radius is ∼2.5 × 10⁴ km, just across the (C+O)/He interface, and is now in the former He-burning shell. Unlike the D9.6 simulations that have essentially two phases to the shock progression (pre- and post-He/H interface), the D10 progenitor's density profile gives rise to four distinct phases of evolution. Fragmentation of the (C+O)/He shell (phase one) is followed by the acceleration of the shock away from carbon encompassed ejecta (phase two), where this material gets injected into the rear of the fragmenting He/H shell after the shock's deceleration (phase three) before slowly expanding to shock breakout (phase four).

Zoom In Zoom Out Reset image size

Starting with phase one, the shock is still briefly accelerating after crossing the (C+O)/He interface, which creates a large reverse shock from the subsequent deceleration once fully into the He-burning shell. This reverse shock is coupled with the location of the mass shell that once marked the (C+O)/He interface, analogous to D9.6's first reverse shock. This promptly shreds and shapes the inner ejecta, as it starts to propagate inward in mass and soon in radius. The aspherical shock hits the (C+O)/He interface at slightly different times, leaving a fingerprint in the form of nonuniform fragmentation in its wake. Four main R-T plumes are quickly spawned (ignoring the poles), which are reminiscent of the dominant plumes that caused the shock to hit this composition interface, unlike the many small R-T plumes seen in D9.6. The main R-T plumes are located at approximately 20°, 35°, 80°, and 130° from the right pole as seen in Figure 15(b) and more distinctly in Figure 15(c). Of this "phase one" material, three out of the four plumes are significantly metal-rich (red features in Figure 16(c) at 20°, 35°, and 130°), while the last bullet is rich in ¹²C (the darker-blue plume at 80°). These features are mirrored owing to the axisymmetry of the initial state, and we will omit the mirror features from further discussion because they exhibit the same behavior as their original counterparts.

Zoom In Zoom Out Reset image size

By 30 s, the shock reaches the density interface at the transition from the He-burning shell to the rest of the He layer (Figures 15(c) and 16(c)). Due to the dramatic change in ρ r³, the shock encountering this shell starts an acceleration that continues until the shock has fully entered the hydrogen envelope (see change in v_Shock starting at ∼3.5 × 10⁵ km in Figure 17). The He-burning shell can be seen throughout the remaining evolution of the explosion as it is propelled forward by the shock. The enhanced ¹²C from partial He burning appears as a dark-blue carbon "bubble" surrounding the inner ejecta in Figure 16(c) at ∼3.0 × 10⁵ km while it interacts with the unburned He shell, similar to the helium bubble surrounding the inner ejecta in the D9.6 simulations. The unburned He shell in Figure 16 appears white (lower than the color map limit of 10⁻³), as it has converted its ¹²C to ¹⁴N via the CNO cycle. At this point we have entered phase two.

Zoom In Zoom Out Reset image size

During phase two, until hitting the He/H interface, the shock front significantly outpaces the inner ejecta. By 60 s, the four main R-T fingers stretch with extremely thin stems at the base while the reverse shock collapses the material behind them. Previously the carbon bubble kept the shape of the shock front, but now the bubble starts to shear at one of the points on its perimeter where the uneven spherical arcs of the shock front hit it prior (at the shock triple point seen earlier at ∼50° in Figures 15(c) and 16(c)). This starts to split the bubble and drive a physical wedge between the inner ejecta, which eventually develops into the dramatic dip seen at much later times in the northeastern quadrant of the dark-blue bubble in Figures 16(d)–(e) and (f).

At ∼250 s, the entirety of the first reverse shock has reached the inner boundary of the grid (now at r = 3.4 × 10⁴ km) as we near the end of phase two. By this point, the inner ejecta has been completely collapsed by the reverse shock, with all shape and distribution being either huddled close to the inner boundary or pushed into the four main R-T fingers. The first reverse shock approaches the inner boundary significantly more centered about the origin than in the D9.6 models; thus, most of it is carried off the grid as opposed to colliding with itself and creating a point-like explosion as seen in the D9.6 models. Because of the irregular shape of the first reverse shock, sections of it reach the inner boundary at different times (with the earliest portion reaching the center at ∼100 s), which further enables the opportunity to evade collision. The only collision that occurs is the portion of the reverse shock produced along the poles that are able to avoid the inner boundary and start to impede the collapsing pressure waves on the opposite side.

Phase three begins at ∼400 s (Figures 15(d) and 16(d)) when the main shock encounters the He/H interface and launches a strong reverse shock owing to the significant shock deceleration (see in Figure 17 the sharp change in v_Shock at ∼4.5 × 10⁶ km). The second reverse shock is coupled to the location of the mass shell of the He/H interface, unlike the decoupled second reverse shock in the D9.6 models. This shell starts to fragment quickly, with an R-T instability forming promptly at the point where the main shock hit unevenly. Since the eastern side of the shock encounters the interface first, this region of the explosion develops its "phase three" R-T plumes the quickest, which can be seen explicitly as the two white instabilities in the northeastern quadrant in Figure 16(e) at a radius of ∼9.0 × 10⁶ km. Additional fragmentation occurs later, forming singular, but dominant, R-T instabilities in succession. In our higher-resolution 2D tests, we find the development of R-T instabilities to be much more abundant and the fragmentation to be much more uniform.

The deceleration of the shock front allows the trailing phase two and phase one material (the carbon bubble and R-T plumes within this bubble, respectively) to eventually get injected into the rear of the fragmenting He/H shell. The carbon bubble achieves this first, as a portion of it first reaches the fragmenting He/H shell and its reverse shock at 3000 s. By 10,000 s (Figure 16(f)), some of the He/H R-T plumes have penetrated the rear of the shock front, and the second reverse shock continues to propagate inward in mass, which allows the remaining regions of the phase two carbon bubble to catch up to it. (The dark-blue bubble in Figure 16(f) catches up to the white.) Note that the two metal-rich phase one R-T plumes in the northeast have at this point merged and burrowed through both the carbon bubble and the second reverse shock (see red R-T plume at ∼40° in Figure 16(f)). At 30,000 s, the remaining phase one R-T instabilities reach and interact with this shell as well. (The 80° and 130° R-T plumes reach the front edges of the dark blue and white at ∼10⁸ km in Figure 16(g).) Additionally, the fragmenting shell, which was once only composed of helium and hydrogen, is now enriched in the phase two carbon. (The blue bubble now occupies the inner anatomy of the previously white R-T plumes in Figure 16(g).)

Phase four is the simplest of all our phases, as most features within the explosion are solely expanding radially. At about 40,000 s, the shock crosses a sudden density spike in the middle of the H shell (ρ r³ spike at ∼1.5 × 10⁸ km in Figure 17). This does not produce a third reverse shock, but it does spawn a noticeable pressure wave that starts propagating inward in mass (and eventually in radius), as the shock experiences a jolt seen as fluctuations in its velocity starting at this point (see v_Shock, solid lines, in Figure 17). Although some of the He/H R-T plumes penetrated the rear of the main shock earlier, they have lost momentum trying to dig their way through the shock and are now being outpaced by it. By 60,000 s, this model has partial shock breakout at the poles and the shock exits the grid along the pole. As these polar flows are artifacts of the assumed symmetry in Chimera, we continue the simulation to determine when the remainder of the shock front would achieve shock breakout. From this point forward, we provide analysis on the wedge of data that exclude the polar regions. (The wedge considers polar angles 30° ≤ θ ≤ 150° across all ϕ.)

The second reverse shock further collapses the phase two carbon bubble, as well as the stems of the phase one R-T plumes within it, as it starts to progress inward in radius at ∼70,000 s. This continues until full shock breakout is achieved when the (nonpole) shock leaves the grid (3.57 × 10⁸ km) at ∼140,000 s (38.8 hr). We rewind the end of our simulation to ∼110,000 s, when the aforementioned "wedge" of the shock enters the region of the progenitor that is partially ionized. By this time, the majority of trailing R-T bullets are at ∼2.0 × 10⁸ km, approximately 12 hr behind the shock front.

The D10–2D3D model keeps its toroidal shape through its entire evolution, like the D9.6–2D3D model. The average velocity of the metal-rich clumps is significantly lower than the average velocity of the shock (see the consistent gap between the green curves in Figure 17). The velocity gap between the two increases when the shock front starts to accelerate down the density gradient as it approaches the He/H interface, which enlarges the relative velocity gap to a difference of ∼7000 km s⁻¹. Although this does not allow for any interaction with the main shock, it does allow for the main R-T clumps to grow rather elongated before encountering the He/H mass shell and reverse shock.

Burrowing through the He/H mass shell is what establishes the final morphology of the CCSN, as this greatly shapes the ejecta and has the ability to spawn further R-T plumes. However, the fragmentation of this shell is quite minimal, and the perturbation from the trailing R-T clumps only seems to add to its bulk at the point of collision. Although some R-T plumes are seeded from this event, the development of the extended structures echoes only the previously trailing asymmetries, rather than having a fully fragmented shell across all angles. Figure 18 shows how different the environment is between D10–2D3D and a high-resolution D10–2D2D model. The D10–2D3D model has three main He/H R-T features forming out of the fragmenting shell as the trailing instabilities catch up to it, while the high-resolution D10–2D2D simulation has numerous R-T plumes developing at the equivalent time.

Zoom In Zoom Out Reset image size

Naturally, the greater number of R-T plumes is not surprising given a much higher resolution, but we provide it here as an example of how the morphology can evolve much differently if the trailing R-T plumes encounter a fragmenting shell equivalent to that of the D10–2D2D high-resolution model. While the bullets in the high-resolution D10–2D2D model still have a fingerprint of the clumps that collided with the second reverse shock, there is a much more complex angular distribution of ejecta with much more mixing close to the rear of the shock front. This complex environment does not occur in the D10–2D3D model (or in D10–2D3D _Tilted, as we will discuss in Section 4.2), which shows a morphological environment that echoes the asymmetries of the past. The D10–2D3D model is eerily similar to its D10–2D2D counterpart of the same resolution (compare middle to left panel of Figure 18). As was apparent with the D9.6–2D3D model, a basic 2D3D mapping does not provide much benefit over running a 2D simulation with similar resolution, due to the absence of longitudinal velocities.

The D10–2D3D model does not seem to accurately portray the long-term evolution of a strongly axisymmetric explosion, due to the lack of initial longitudinal velocities and exaggerated polar flows from an unfortunate interaction between the 2D Chimera model's polar flow and the excised cone in FLASH. The D10–2D3D _Tilted model alleviates the interaction with the excised cone, though the polar flow itself is still present, as it has been tilted fully onto the FLASH grid.

At first glance, the first column of the yields in Figure 19 does not show much change in the ejecta distribution in velocity space between D10–2D3D and D10–2D3D _Tilted. The dominance of the poles in both models drowns out contributions from the rest of the ejecta to the higher-velocity matter and hides the microstructure in the first column of Figure 19. Because the poles in the D10 model are so dramatic, this provides a counterexample to the argument that the yields distribution in the D9.6–2D3D _Tilted model could potentially be misleading owing to more of the polar flow being present on the grid compared to its respective D9.6–2D3D model. If that were the case, then we would see a more dramatic difference in the distribution of the ejecta when comparing the top and bottom panels of the first column in Figure 19. Clearly, we do not.

Zoom In Zoom Out Reset image size

To reveal microstructure obscured by the poles, we further analyze the yields by considering a wedge of the models that excludes contributions on the grid from the polar flows (second and third columns of Figure 19). The wedge for D10–2D3D _Tilted is the same wedge discussed in Section 4.1 for the D10–2D3D model (30° ≤ θ ≤ 150° across all ϕ) but is applied after a 90° coordinate transform (i.e., after "undoing" the tilt). Through this, we actually see more of an effect that tilting the model has provided, as D10–2D3D _Tilted has an apparent higher-velocity tail (∼2500 km s⁻¹) compared to its D10–2D3D counterpart (∼1900 km s⁻¹) when comparing models in the second column of Figure 19. Comparing models in mass space (third column of Figure 19) shows higher yields for D10–2D3D _Tilted in the outer regions. This can clearly be seen in the extent of ⁵⁶Ni and ⁴⁴Ti, which both drop significantly in the D10–2D3D model at 7.5 M_☉ (top row, third column of Figure 19). In contrast, for the D10–2D3D _Tilted model both ⁵⁶Ni and ⁴⁴Ti extend to 8.5 M_☉, joining the lighter elements in the ejecta (bottom row, third column of Figure 19). The total yields (Table 4) further reveal this difference, with roughly 6% and 13% greater ⁴⁴Ti and ⁵⁶Ni yields, respectively, in the D10–2D3D _Tilted model. Due to more of the polar flows, which originate from the hot bubble, being included on the grid, these isotopes (plus ⁵²Fe) are some of the key differences relative to the D10–2D3D model, while the rest of the yields are relatively consistent between D10 models. Although the poles are excised for both models in Table 4, the relevant species are more heavily mixed into surrounding areas during the explosion in the D10–2D3D _Tilted model; thus, these species are more abundant than for D10–2D3D.

The consistency of the ejecta for the two D10 models is matched by the consistency in shock progression (colored solid curves in Figure 20). Even the average radii of the ${X}_{{}^{56}\mathrm{Ni}+\mathrm{IG}}$ isosurfaces (colored dashed curves) are quite similar. Despite this, the D10–2D3D _Tilted model develops more spherical-bubble structures during its evolution (Figure 21, right), due to the initial longitudinal and latitudinal velocities. This is consistent with what happened in the D9.6–2D3D _Tilted model. D10–2D3D _Tilted is slightly less axisymmetric than D10–2D3D in Figure 21 and has more structure in its central and outer regions. Therefore, these metal-rich clumps in the D10–2D3D _Tilted model retain slightly higher velocities (dashed orange curve in Figure 17) over its D10–2D3D counterpart (dashed green curve in Figure 17) until the He/H interface when the shock starts to decelerate and the second reverse shock forms. The second reverse shock dictates the subsequent velocity profile of the clumps, limiting their velocities as they try to burrow through it, bringing the average clump velocities back together as both dashed curves decrease until shock breakout. Although the average velocities of the clumps for both models obtain similar values near shock breakout, the overall velocity distribution across the analysis wedge domain (Figure 22) shows that the D10–2D3D _Tilted model still retains higher velocities in the outer envelope.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dynamics of the small features are further demonstrated by the growth in the isosurface areas (Figure 23). The total area for both the external coating (0.1% isosurface) and inner anatomy (3% isosurface) of the ⁵⁶Ni+IG-rich plumes start to diverge early during the dramatic acceleration of the shock. After encountering the reverse shock at ∼10,000 s, the total surface area represented by the external coating (0.1% isosurface) of the bullets diverges further, as the bullets in D10–2D3D _Tilted are able to burrow through it more efficiently owing to the somewhat higher velocities that result from the spherical-like structures created upon the deviation from axisymmetry. The second divergence between models is not present in the 3% isosurface (inner anatomy) curves. This is not surprising, due to the relatively similar distribution of metal-rich ejecta in both simulations, with the key differences occurring at larger mass coordinates and higher velocities that are inherently captured by the external coating isosurface instead. As with D9.6, the contributions to the total surface area converge back toward those provided by the expansion of the shock and more dramatically for D10–2D3D, which stays more axisymmetric and lacks the spherical-bubble structures that retain higher velocities and prolong the convergence to the shock-driven area increase. In contrast to the D9.6 models, the 0.1% isosurface in the D10 models has a larger total surface area than the 3% isosurface owing to the considerable amount of fallback caused by the reverse shocks combined with a more condensed angular distribution of the metal-rich ejecta due to fewer R-T plumes spanning the whole volume.

Zoom In Zoom Out Reset image size

The D10–2D3D _Tilted model achieves greater velocities compared to D10–2D3D, although not as striking as D9.6–2D3D _Tilted. We believe that this is less dramatic in this simulation owing to the resolution-limited spawning of only a few clumps at the He/H fragmentation, whereas D9.6 has a wider range of bullets developing from its (C+O)/He fragmentation. We would expect a larger deviation from the D10–2D3D toroid shape if the fragmentation environment was more similar to the high-resolution test of Figure 18. Most importantly, despite all this, the end result of the D10–2D3D _Tilted model no longer looks like a rotated 2D model and provides a more faithful 3D model of a polar-dominated explosion than the untilted D10–2D3D model.