Three-dimensional Hydrodynamics Simulations of Precollapse Shell Burning in the Si- and O-rich Layers

As already mentioned above, we employ two 1D progenitors from Yoshida et al. (2019): one is a 22 M_⊙ star with an extended Si/O-rich layer (hereafter model 22L), and the other is a 27 M_⊙ star with a more extended O–Si layer (hereafter model 27L_A) than model 22L. Note that "L" in the model name represents the choice of the overshooting parameter, which is calibrated to explain early B-type stars in the Large Magellanic Cloud (Brott et al. 2011), and the subscript "A" denotes the inclusion of a convective overshoot during the advanced stages.⁶ Both models are nonrotating (see Yoshida et al. 2019; Takahashi et al. 2016, 2018, for more details about the 1D stellar evolution code "HOSHI" ).

Figure 1 shows the initial structure for the 3D simulations of model 22L (left panel) and model 27L_A (right panel). These initial structures correspond to the 1D structures at ∼130 s before the central temperature is 10¹⁰ K. In four rows from the top to the bottom, we show the mass fraction profiles, the radial velocity profiles, the profiles of Brunt–Väisälä frequency and turnover time of the convection, and the temperature and entropy profiles. Brunt–Väisälä frequency in convective layers is defined as

$$\begin{eqnarray}&&{\omega }_{\mathrm{BV}}=\sqrt{-\displaystyle \frac{{{GM}}_{r}}{{r}^{2}}\displaystyle \frac{\delta }{{H}_{P}}({{\rm{\nabla }}}_{\mathrm{ad}}-{\rm{\nabla }}+\displaystyle \frac{\varphi }{\delta }{{\rm{\nabla }}}_{\mu })},\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, M_r is the mass coordinate, r is radius, $\delta \equiv -{(\partial \mathrm{ln}\rho /\partial \mathrm{ln}T)}_{P,\mu }$, H_P is pressure scale height, ${{\rm{\nabla }}}_{{\rm{ad}}}\equiv {({\rm{\partial }}{\rm{ln}}T/{\rm{\partial }}{\rm{ln}}P)}_{{\rm{ad}}}$ is the adiabatic temperature gradient, ${\rm{\nabla }}\equiv \partial \mathrm{ln}T/\partial \mathrm{ln}P$ is the temperature gradient, $\varphi \,\equiv {(\partial \mathrm{ln}\rho /\partial \mathrm{ln}\mu )}_{P,T}$, and ${{\rm{\nabla }}}_{\mu }\equiv \partial \mathrm{ln}\mu /\partial \mathrm{ln}P$. Note in this definition that ${\omega }_{\mathrm{BV}}^{2}\gt 0$ corresponds to convective instability. The turnover time of convection in a convective layer is calculated using the equation t_turn = r_cv/v_cv, where r_cv is the width of the convection layer and v_cv is the convection velocity calculated using MLT.

Zoom In Zoom Out Reset image size

The left panel shows the initial structure of model 22L. The top row indicates the mass fraction profiles. The green hatched region shows the location of the Si/O layer (3 × 10⁸ cm ≤ r ≤ 1.2 × 10⁹ cm and 1.77 M_⊙ ≤ M_r ≤ 2.71 M_⊙). The width of the Si/O layer is 9 × 10⁸ cm. This layer is convective owing to O shell burning. In fact, the temperature of the base of the layer (at r ∼ 3 × 10⁸ cm, the left end of the hatched region) is about 3 × 10⁹ K (see the black line in the bottom left panel), which exceeds the ignition temperature of O burning. The entropy profile (red line in the bottom left panel) is nearly flat in the Si/O layer (shaded region), which is a result of convective energy transport treated in the 1D stellar evolution calculation. Note that the regions inside and outside the shaded region correspond to the Fe/Si core and O/Ne layer, respectively. Both of the regions are convectively inert in the 3D simulation. The former is due to a positive entropy gradient in the core; the latter is due to too short computational timescale for our 3D simulation to see the convective activity for the outer region (i.e., longer convective timescale; see the third row, left panel).

Note that the composition structure of model 22L is similar to model 25M, of which 2D and 3D hydrodynamics simulation was investigated in Yoshida et al. (2019). In the previous study, the 2D dynamical simulation has also been done for model 22L, which gives a similar result to the 2D simulation of model 25M. The maximum turbulent Mach number reaches 0.108, and the spectrum analysis showed a peak of the radial turbulent velocity at a low mode of ℓ = 2. For model 22L, we start the 3D calculation mapping from the 1D data when the temperature and density at the center are 6.40 × 10⁹ K and 1.32 ×10⁹ g cm⁻³, respectively. We follow the 3D evolution for 65.5 s until the onset of CC.

The right panel is for model 27L_A. It has an extended oxygen- and silicon-rich (O/Si-rich) layer at 5 × 10⁸ cm ≤ r ≤ 60 × 10⁸ cm (2.53 M_⊙ ≤ M_r ≤ 7.99 M_⊙). The inner neon-depleted region (r < 5.8 × 10⁸ cm) is formed by oxygen shell burning after oxygen core burning. The outer region (r > 5.8 × 10⁸ cm) with a small amount of neon is the former O/Ne layer. Neon in this region (5.8 × 10⁸ cm ≤ r ≤ 60 ×10⁸ cm) is burned through Ne shell burning after the Si core burning phase but still remains after the Fe core formation. We call this layer the O/Si/Ne layer, which is indicated by the green hatched region. The O/Si/Ne layer is convective and has constant entropy profile (red line). The shell convection is powered by Ne shell burning, as the base temperature exceeds the ignition temperature of Ne burning of 2 × 10⁹ K (black line).

The 2D model of 27L_A showed the maximum turbulent Mach number 0.179 with a peak mode of the radial Mach number ℓ = 2 in Yoshida et al. (2019). The convective region extends to the whole O/Si/Ne layer. Note that the region outside and inside of this layer is convectively inert in the 3D simulation owing to the same reason as for model 22L as mentioned above. For this model, we start the 3D calculation by mapping from the 1D data when the temperature and density at the center are 6.87 × 10⁹ K and 9.44 × 10⁸ g cm⁻³, respectively. We follow the 3D evolution for 218.51 s until the onset of CC.

The numerical schemes in this work are the same as those in Yoshida et al. (2019). We employ the 3DnSEV code, which solves Newtonian hydrodynamics equations. The 3D models are computed on a spherical polar coordinate grid with a resolution of n_r × n_θ × n_ϕ = 512 × 64 × 128 zones. The radial grid is logarithmically spaced and covers from the center up to the outer boundary of 10¹⁰ cm. For the polar and azimuthal angle, the grid covers all 4π sr. For boundary conditions, reflective boundary conditions are adopted for the inner boundary. Fixed-boundary conditions are used for the outer boundary, except for the gravitational potential that is inversely proportional to the radius at outer ghost cells. To focus on the convective activity of the shaded regions in Figure 1, the inner 1000 km is solved in 1D. Gravity is included by assuming a 1D, monopole, gravitational potential. Such a treatment is indispensable for reducing the computational time; the nonlinear coupling between the core and the surrounding shells (e.g., Fuller et al. 2015) is beyond the scope of this study.

A piecewise linear method with geometrical correction of the spherical coordinates is used to reconstruct variables at the cell interface, where a modified van Leer limiter is employed to satisfy the condition of total variation diminishing (Mignone 2014). The numerical flux is calculated by the HLLC solver (Toro et al. 1994). The consistent multifluid advection (CMA) method of Plewa and Müller (1999) is used for evaluating the numerical flux of isotopes.

We use the "Helmholtz" equation of state (Timmes & Swesty 2000). Neutrino cooling is taken into account (Itoh et al. 1996) as a sink term in the energy equation. As for nuclear reaction, a reaction network of 21-isotopes (aprox21; Paxton et al. 2011) is applied when the temperature is lower than 5 × 10⁹ K. This network includes ⁵⁴Fe, ⁵⁶Fe, and ⁵⁶Cr, which are crucial to treat a low electron fraction Y_e ≳ 0.43 in the pre-SN stage. Also, it is as large as that of Couch et al. (2015) and a little larger than 19-isotopes of Müller et al. (2016) and Yadav et al. (2020). When the temperature is higher than 5 × 10⁹ K, the chemical composition is calculated assuming nuclear statistical equilibrium (NSE) instead of solving the nuclear reaction network.

Figure 2 shows the electron fraction profile as a function of the mass coordinate in M_r < 2.0 M_⊙ at the last step when the central temperature reaches 10¹⁰ K for models 22L and 27L_A. For both models neutronization proceeds more slowly in the 3D simulations compared with the 1D simulations. To correctly handle the neutronization of heavy elements from Si to the iron group and the gradual shift of the nuclear abundances, one needs to use a sufficient number of isotopes (∼100; Arnett & Meakin 2011), which is currently computationally and technically very challenging. Since the NSE region appears mainly in the Fe core, this treatment may not significantly affect our results in which we focus on convection in the outer layers (e.g., Si/O-rich and O/Si/Ne layers).

Zoom In Zoom Out Reset image size

When we start the 3D runs, seed perturbations to trigger nonspherical motions are imposed on the 1D data by introducing random perturbations of 1% in density on the whole computational grid. We terminate our 3D runs when the central temperature exceeds 10¹⁰ K, because the core is dynamically collapsing at this time.

Figure 3 shows the spatiotemporal evolution of the angle-averaged turbulent Mach number (left panels, 〈Ma²〉^1/2(r)) and the turbulent velocity in units of 10⁷ cm s⁻¹ (right panels, v_turb,7) for models 22L (top panels) and 27L_A (bottom panels), respectively. Here the definition of 〈Ma²〉^1/2(r) is the same as in Equation (8) in Yoshida et al. (2019),

$$\begin{eqnarray}&&{\left\langle {\mathrm{Ma}}^{2}\right\rangle }^{1/2}(r)={\left[\displaystyle \frac{\int \rho \{{\left({v}_{r}-\langle {v}_{r}\rangle \right)}^{2}+{v}_{\theta }^{2}+{v}_{\phi }^{2}\}d{\rm{\Omega }}}{\int \rho {c}_{s}^{2}d{\rm{\Omega }}}\right]}^{1/2},\end{eqnarray} \tag{ 2 }$$

where ρ is the density; v_r, v_θ, and v_ϕ are radial, tangential, and azimuthal velocities, respectively; 〈v_r〉 is the angle-averaged radial component velocity obtained by Reynolds averaging; Ω is the solid angle; and c_s is the sound speed. The definition of the turbulent velocity v_turb is defined as

$$\begin{eqnarray}&&{v}_{\mathrm{turb}}={\left[\displaystyle \frac{\int \rho \{{\left({v}_{r}-\langle {v}_{r}\rangle \right)}^{2}+{v}_{\theta }^{2}+{v}_{\phi }^{2}\}d{\rm{\Omega }}}{\int \rho d{\rm{\Omega }}}\right]}^{1/2}.\end{eqnarray} \tag{ 3 }$$

Figure 4 is the time evolution of the turbulent kinetic energy E_kin,turb in the Si- and O-rich region defined as

$$\begin{eqnarray}&&{E}_{\mathrm{kin},\mathrm{turb}}=\int \displaystyle \frac{1}{2}\rho \{{({v}_{r}-\langle {v}_{r}\rangle )}^{2}+{v}_{\theta }^{2}+{v}_{\phi }^{2}\}{r}^{2}{drd}{\rm{\Omega }}.\end{eqnarray} \tag{ 4 }$$

The inner and outer boundaries of the Si- and O-rich region are determined by the conditions of the O mass fraction larger than 0.01 and the Si mass fraction larger than 0.1, respectively. Note that we use the term turbulent velocity as a velocity fluctuation in this paper even when the fluctuation has not developed to turbulence.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The top panels of Figure 3 show that in model 22L turbulence starts to develop from t = 10 s at the base of the Si/O layer (r ∼ 3 × 10⁸ cm), which is driven by O shell burning. Note that t = 0 is defined as the epoch of the start of the 3D simulation. As seen, turbulent flows with the Mach number of ∼0.1 (greenish region) spread over the entire Si/O layer (up to the upper white dashed line) before the final simulation time of this model, t = 65.5 s (before CC onset). At t ∼ 50 s, a strong turbulence of the Mach number ∼0.16 and the turbulent velocity ∼10⁸ cm s⁻¹ develops transiently from the base of the Si/O layer (again triggered by O shell burning). Figure 4 shows that the kinetic energy of the turbulent motion increases after ∼50 s. At this time, the nuclear energy generation rate at the bottom of this convective layer reaches ∼6 × 10¹⁶ erg g⁻¹ s⁻¹, whereas it is less than 1 × 10¹⁶ erg g⁻¹ s⁻¹ before that time. This kind of episodic burning causes a slight expansion and the subsequent contraction of the above layer, which was already found in 3D models of 25M (Yoshida et al. 2019) and also of the 18.8 M_⊙ star by Yadav et al. (2020). In Figure 19 in Appendix B, we show 2D slices on the x-z plane of the radial turbulent velocity distributions at 10, 40, 50, and 65.5 s in model 22L.

Note the difference of the radial scale between the top and bottom panels of Figure 3. The turbulent flows develop in the more extended region for model 27L_A, whereas the maximum Mach number (∼0.1) is smaller than that of model 22L (compare the left and right panels). In model 27L_A, the temperature at the bottom of the convective layer is lower and the neon shell burning occurs. The nuclear generation rate by the neon shell burning (model 27L_A) is much smaller than the oxygen shell burning (model 22L). Given the initial composition profile (i.e., the broader shaded region; see Figure 1), this may not be very surprising. However, such a large-scale mixing has been first obtained in our model series. Hence, we pay particular attention to this model. We show 2D slices on the x-z plane for the radial turbulent velocity distributions at 10, 90, 140, and 218.5 s for model 27L_A in Figure 19 in Appendix B.

For better understanding, we divide the entire evolution of model 27L_A into the three phases, namely, phase I (0 s ≤ t ≤ 107 s), phase II (107 s ≤ t ≤ 127 s), and phase III (127 s ≤ t ≤ 218 s), respectively. In phase I, turbulence starts from the base of the O/Si/Ne layer (r ∼ 6 × 10⁸ cm at t = 10 s), which gradually extends outward. In phases II and III, the more enhanced turbulence stems, the maximum angle-averaged turbulent velocity is 4.2 × 10⁷ cm s⁻¹, from stronger (Ne) shell burning. We explain more in detail in the next section.

We move on to explore how the turbulent activity shown in Figure 3 is triggered in association with the change in the element composition in the burning shells.

3.2.1. Model 22L

First, let us focus on the mass fraction distribution for model 22L. Figure 5 shows the spatiotemporal evolution of the angle-averaged O mass fraction for model 22L. One sees that shortly after the start of the calculation (t = 0), the oxygen-depleted materials are mixed into the Si/O layer (0 s ≤ t ≤ 20 s, shown as a region below the cyan line). Subsequently, the turbulent mixing encompasses the entire Si/O layer (3 × 10⁸ cm ≤ r ≤ 12 × 10⁸ cm, corresponding to the shaded region in the left panel of Figure 1) up to t ∼ 35 s with the diminishing turbulent activity (see the top panels of Figure 3). After t ∼ 50 s, the O mass fraction decreases at the bottom of the Si/O layer owing to (the episodic) O shell burning, and the turbulent flows radially growing outward result in the reduction of the O mass fraction in the outer layer (see the bluish region for 50 s ≤t ≤ 60 s and Figure 21).

Zoom In Zoom Out Reset image size

3.2.2. Model 27L_A

Second, we will show the mass fraction distributions of O, Ne, and Si in model 27L_A. We discuss shell burning processes in the Si/O-rich and O/Si/Ne layers and the development of the turbulent motion in these layers.

Ne mass fraction.—Figure 6 shows the spatiotemporal evolution of the Ne mass fraction. The bottom panel is focusing on the region close to the lower boundary of the O/Si/Ne layer. Comparing with the bottom right panel of Figure 3, a cautious reader would notice that the reduction of the Ne mass fraction is associated with the high turbulent velocity in this layer. We will describe details below.

Zoom In Zoom Out Reset image size

In phase I (0 s ≤ t ≤ 107 s), the Ne burning near the base of the layer close to the lower horizontal line results in the reduction of the Ne mass fraction in the region r ≲ 1.5 × 10⁹ cm. Looking at the same region to the right up to t ∼ 50–100 s, one can see that the radial gradient of the Ne mass fraction becomes less steep owing to turbulent matter mixing, by which the Ne-rich material is mixed downward via downflows into the Ne-poor region. We also see the above trend in the top panel of Figure 7, which indicates the radial profiles of the angle-averaged Ne mass fraction at definite times with intervals of 30 s. The growth of the turbulent region from r ∼ 7 × 10⁸ cm at t ∼ 30 s can also be depicted in the spacetime diagram of the turbulent velocity, which continues into phase III up to the radius of the outer boundary of the layer at r ∼ 60 × 10⁸ cm. In phase II and also in phase III, the growth of high turbulent activity is observed as a red region above the shell (the lower horizontal white dashed line). This is again predominantly because of the Ne shell burning there. The reduction of the Ne mass fraction is seen in the regions r ≲ 1 × 10⁹ cm and r ∼ 4 × 10⁹ cm around t ∼ 150 s in the top panel. Besides, one can see the episodic burning below the base of the O-rich layer (see the two red regions near t ∼ 100 s and 127 s ≤ t ≲ 150 s below the lower horizontal white dashed line). This plays a key role for enhancing the turbulent mixing as we will explain in what follows.

Zoom In Zoom Out Reset image size

Silicon and oxygen mass fractions in the region close to the inner boundary of the O/Si/Ne layer.—Figure 8 is for the Si (left panel) and O (right panel) mass fraction of model 27L_A, focusing on the region close to the inner boundary of the O/Si/Ne layer (note the location of the horizontal white dashed line at 5.8 × 10⁸ cm). By these zoom-in plots, we focus on why the episodic expansion occurs, which is seen as the three bumps in the red regions in phases I, II, and III in the left panel.

Zoom In Zoom Out Reset image size

From the left panel of Figure 8, one can see that the gradual decrease of the Si mass fraction takes place at 2 × 10⁸ cm ≲ r ≲ 3 × 10⁸ cm, M_r ∼ 2 M_⊙, at 0 s ≤ t ≲ 107 s in phase I. Figure 9 shows that the Si mass in this region (1.8 M_⊙ ≲M_r ≲ 2.2 M_⊙) gradually decreases for the former ∼50 s. This gradual Si shell burning causes the first bump, which can be seen as an upward lift of the Si-rich material up to M_r ≲ 2.6 M_⊙. Looking at the same region to the right (namely, phase II), one can see the abrupt decrease of the Si mass fraction (the color changes from orange to blue/green), which is due to the episodic Si shell burning. This also reflects the decrease in the total Si mass in this region (see Figure 9). The energy release results in the expansion and the following contraction above the Si-rich layer with M_r ≳ 2.1 M_⊙. This expansion and contraction can be also visible in the O mass fraction in phase II (right panel of Figure 8).

Zoom In Zoom Out Reset image size

In phase III (127 s ≤ t ≤ 218 s), a dorm-like transient expansion and contraction can be seen in Figure 8. This is predominantly triggered by O shell burning in the range of M_r = 2–2.6 M_⊙, which corresponds to the bluish low O mass fraction regions in phase III. Figure 10 shows the time evolution of the total O mass in the range of 2.0–2.6 M_⊙, where the time variation of the O mass fraction is mainly seen. We also see the decrease in the oxygen mass in this region in 130–180 s. On the other hand, the Si mass fraction increases in the region M_r ∼ 2.2 M_⊙ by the O shell burning. The Si mass in the region 1.8–2.2 M_⊙ also increases between ∼130 and 150 s (see Figure 9). Because of the O shell burning and inward turbulent mixing, the Si mass fraction increases in the dorm-like region concurrently (see the red region in the top right panel).

Zoom In Zoom Out Reset image size

Silicon mass fraction in the O/Si/Ne layer.—Let us move to the viewpoint of the Si mass fraction in the O/Si/Ne layer. The top left panel of Figure 11 is for the Si mass fraction for the entire O/Si/Ne layer. Si is an ash of the Ne and O burning. By inspecting the Si mass fraction distribution, hereafter we try to understand the complicated interplay of the Ne burning (top right panel) and the O burning (middle right panel), which take place in a different manner.

Zoom In Zoom Out Reset image size

In phase I (0 s ≤ t ≤ 107 s), the increase of the Si mass fraction, which takes place from the near base region (10 ×10⁸ cm) to the outer region (20 × 10⁸ cm), is clearly visible. This Si is produced through the Ne shell burning at the bottom of this layer. By comparing with the extension of the turbulent activity shown by the bottom panels of Figure 3, the correspondence between the two regions is confirmed. This naturally supports that the turbulent mixing globally takes place in the extended O/Ne/Si layer.

We note that the mass shells from 2 to 8 M_⊙ move in a synchronized fashion in each phase. This expansion and contraction propagates from the bottom of the Si/O-rich layer inward and outward with the sound velocity of (5–7) × 10⁸ and (2–5) × 10⁸ cm s⁻¹, respectively. When an episodic burning occurs, the temperature and the pressure gradient at the burning region increase. The regions inside and outside the burning shell expand to achieve a hydrostatic equilibrium. In the region inside the Si/O-rich layer, the structural change propagates via the sound wave propagation within 1 s.

2D slices and 3D evolution.—Figure 12 depicts the evolution of the 2D slice on the x-z plane for the Si mass fraction for four time snapshots from 120 to 180 s in phases II and III. The 3D contours of the Si mass fraction at the last step are also shown in Figure 13. The size of the Si-rich region grows up through the O shell burning and the expansion as shown in Figure 11. In phase III, the size of the inhomogeneities of the Si mass fraction in an inner region of the O/Si/Ne layer is shown to become bigger with time. This is because of the stronger turbulence mixing encompassing the whole layer (e.g., the bottom panels of Figure 3 for the region between the two horizontal white dashed lines in phase III) predominantly triggered by the Ne shell burning (e.g., bottom right panel of Figure 11).

Download figure:

Video Standard image High-resolution image

Download figure:

Video Standard image High-resolution image

Similar to Figure 12, Figure 14 shows 2D slices for the Ne mass fraction. The 3D contours of the Ne mass fraction at the last step are also shown in Figure 15. In the relatively early phase of phase III (at ∼140 s), one can see yellow to red regions in the northeast and southwest directions, which shows the penetration of the high Ne fraction matter downward into the low Ne fraction region (greenish and bluish regions). This convection mixing develops more strongly with time (see the right panel). In fact, the Ne mass fraction is homogenized to a large extent in the O/Ne/Si layer as already shown in Figure 7 (see the radial composition profile at ≳210 s).

Download figure:

Video Standard image High-resolution image

Download figure:

Video Standard image High-resolution image

Radial composition profile and angular dispersion.—Finally, in Figure 16, we compare the composition distributions of representative α-elements at the initial (dashed lines) and final time steps (solid lines) of the 3D simulations of model 27L_A. In the corresponding animations, we also present the accompanying two thin lines, which denote the angular dispersion relative to the angle-average one (thick line). We note that in the 1D evolution the radial profiles of the O, Ne, and Si mass fractions scarcely change during the final ∼100 s up to the central temperature of 10¹⁰ K.

Download figure:

Video Standard image High-resolution image

Regarding the O mass fraction (blue lines), the angular dispersion is less than 2% in the O/Si/Ne layer (e.g., the shaded region of 7 × 10⁸ cm ≲ r ≲ 6 × 10⁹ cm and 2.7 M_⊙ ≲M_r ≲ 8.4 M_⊙). Since the O mass fraction is large and almost constant in this layer at the initial step of the simulation, the effect of turbulent mixing is not seen in the O mass fraction profile except for the region close to the base of the layer (as seen from the deviation from the thin dashed line and the thick solid line).

As is expected from the global mixing in the O/Si/Ne layer (Figures 14 and 15), the initial Ne composition that gradually increases with the radius in the layer (see the dashed magenta line) is flattened owing to the mixing at the final simulation time (see the solid magenta line). This means that the 3D model shows smaller radial dependence than the 1D model for the composition profile in the O/Si/Ne layer at the last step. Regarding the angular dispersion, it becomes larger compared to those of O (blue lines) and Si (green lines). It is about 10%–20% in the inner layer of r ∼ (7–20) × 10⁸ cm and M_r ∼ 2.7–4.4 M_⊙, reaching 20%–30% near the upper boundary of the layer. Regarding the Si mass fraction (green lines), the initial distribution (thin dashed green line) is slightly shifted outward, and the initial composition gradient in the layer is also flattened owing to the mixing (compare the thin with thick green lines at 6 × 10⁸ cm ≲ r ≲ 6 × 10⁹ cm and 2.6 M_⊙ ≲ M_r ≲ 8.4 M_⊙). The angular dispersion is less than 10%–20% in the layer.

Irrespective of the composition differences, it can be commonly observed that the angular dispersion becomes bigger at the layer interface, where the nuclear burning occurs vigorously and the multidimensional effects become significant. In fact, the largest dispersion of 20%–30% of the Ne (magenta line) is obtained near the boundary between the O/Ne/Si layer ((6–7) × 10⁸ cm) and the O/C layer (∼6 × 10⁹ cm), and the dispersion of the O mass fraction distribution (blue line) becomes bigger at the interface between the O/Ne/Si layer and the Si layer inside.

The effect of shell burning on the turbulent mixing can be found from the time variation of the mass fraction profile in the animated version of Figure 16. When a shell burning occurs, the radial distribution of angle-averaged mass fractions is homogenized. On the other hand, large aspherical turbulence enhances the angular dispersion of mass fractions. Then, the dispersion becomes smaller as the turbulence weakens.

We explain the variation of O and Ne mass fraction distributions in phase III as an example. In 127 s ≤ t ≲ 150 s, strong turbulence is activated in the Si/O layer with (4–6) × 10⁸ cm (see also the right panel of Figure 3). The steep rise of the O mass fraction with radius in this region becomes less steep, and at the same time, the O mass fraction decreases with time. On the other hand, an increase in the angular dispersion is seen. The turbulence weakens after t ∼ 150 s, and this dispersion also becomes smaller (see the paragraph on Ne mass fraction). For the Ne mass fraction, we see the radial gradient in the O/Si/Ne layer at the beginning of phase III. The Ne mass fraction becomes gradually homogenized outward from r ∼ 7 × 10⁸ cm, and the angular dispersion becomes larger. The angular dispersion in the radially homogenized region is about 40% at t = 150 s. This dispersion gradually decreases after this time.

In what follows, we attempt to make an order-of-magnitude estimate to explore whether the 3D mixing (leading to the flattening of the initial composition distribution) could or could not be treated as a diffusion process, which an evolution code commonly does. The initial 1D data have typical values of the convective velocity v_cv ∼ 2 × 10⁷ cm s⁻¹, the mixing length l_cv = αH_p ∼ 5.5 × 10⁸ cm, and the diffusion coefficient D_cv = v_cvl_cv/3 ∼ 3.6 × 10¹⁵ cm² s⁻¹ in the O/Ne/Si layer of model 27L_A. Thus, the diffusion timescale to homogenize the composition gradient is estimated as τ_cv ∼ L²/D_cv ∼ 6.3 × 10² s, where the size of the composition gradient of Ne of L ∼ 1.5 × 10⁹ cm is used. This is much longer than the evolutionary timescale in this late phase, and it is why the original mass fraction distribution has a gradient within the convective region. However, our 3D simulation results in the radially homogeneous distributions of the angle-averaged Ne mass fraction in a timescale of ∼100 s by the turbulent mixing. We have noticed that this short timescale is rather compatible to the advection timescale of the turbulent flow, τ_turb ∼ L/v_turb ∼ 40 s. Much more systematic 3D stellar evolution modeling is preferable to seek for the method to mimic the multidimensional effects in mixing models in 1D stellar evolution models (e.g., Couch et al. 2015; Müller 2016; Müller et al. 2018; Yadav et al. 2020).

In what follows, we propose a useful method in order to extract the information about the typical eddy size of turbulence and its relation to the turbulent energy in the burning shells of our two 3D models. In our previous study (Yoshida et al. 2019), we performed a spectrum analysis of the turbulent velocity at the radius where the maximum Mach number was obtained (e.g., Figure 6 in that paper). If we follow the conventional treatment in model 27L_A, the location of the maximum Mach number is identified at ∼3 × 10⁹ cm (see the right panel of Figure 3). On the other hand, the O/Si/Ne layer with the similar amplitude of the Mach number extends in a wider region up to r ∼ 6 × 10⁹ cm near the upper boundary of the layer (the upper horizontal white dashed lines). In such a widely extended shell, it should be better to study the radial dependence of the spectrum c_ℓ(r) of the turbulent velocity, instead of estimating the power spectrum at one radial location with the maximum Mach number.

We calculate the power spectrum of the radial turbulent velocity as

$$\begin{eqnarray}&&{c}_{{\ell }}^{2}(r)=\sum _{m=-{\ell }}^{{\ell }}{\left|\int ({v}_{r}-\langle {v}_{r}\rangle ){Y}_{{\ell }m}^{* }(\theta ,\phi )d{\rm{\Omega }}\right|}^{2},\end{eqnarray} \tag{ 5 }$$

where ${Y}_{{\ell }m}^{* }(\theta ,\phi )$ is the (complex conjugate) spherical harmonics of degree ℓ and order m and ℓ is from 1 to 50. We investigate the peak mode ℓ_peak(r) as a function of radius. The (specific) turbulent energy, the energy of the radial turbulent velocity, is evaluated as

$$\begin{eqnarray}&&{\varepsilon }_{r,\mathrm{turb}}(r)=\displaystyle \frac{1}{8\pi }\sum _{{\ell }}{c}_{{\ell }}^{2}(r).\end{eqnarray} \tag{ 6 }$$

Figure 17 shows the radial profiles of ℓ_peak(r) (solid lines) and the turbulent energy ε_r,turb (dashed lines) at two representative times near (blue line) and at the final simulation time (red line) for models 22L (top panel) and 27L_A (bottom panel), respectively. From the top panel, one can see in model 22L that the value of ℓ_peak(r) (solid lines) is in the range of ∼3–5 in most of the Si/O layer (shaded region). The smaller turbulent mode ℓ_peak(r) ∼ 3 is associated with the higher turbulent velocity (the dashed lines), whereas near the outer boundary of the layer (e.g., the right edge of the shaded region) the turbulent mode becomes bigger (ℓ_peak(r) = 7) with the smaller turbulent energy (e.g., seen as a drop in the dashed lines). The high turbulent energy near the base of the Si/O layer (see also the right panel of Figure 3) is triggered by O shell burning. We consider that the low ℓ_peak(r) is a natural outcome of the growing large-scale turbulence near the base of the burning layer. One also sees the value of ℓ_peak ∼ 8 at the bottom of the Si/O layer (see the blue solid line). This corresponds to the start of the turbulence development by the ignition of the O shell burning. As a side remark, the high ℓ_peak(r) (≲12) below the lower convective boundary may be not surprising because of the weak turbulent activity and also the thinner width of the region.

Zoom In Zoom Out Reset image size

The above trend is similar also to model 27L_A (bottom panel of Figure 17). The minimum ℓ_peak(r) = 2 is obtained near the base of the O/Si/Ne layer, which is smaller than model 22L, whereas ℓ_peak(r) is found to be slightly higher (∼8) near the outer boundary than for model 22L (ℓ_peak(r) ∼ 5).

In order to characterize a typical scale of turbulent flows in the whole Si/O or O/Si/Ne layer (for models 22L and 27L_A, respectively), we define the average mode of ℓ_peak, which we obtain from ℓ_peak(r) by integrating with a weight of the turbulent energy as follows:

$$\begin{eqnarray} &&\langle {\ell }_{{\rm{peak}}} \rangle =\displaystyle \frac{\int {\ell }_{{\rm{peak}}}(r)\rho (r){\varepsilon }_{r,{\rm{turb}}}(r){r}^{2}{dr}}{\int \rho (r){\varepsilon }_{r,{\rm{turb}}}(r){r}^{2}{dr}},\end{eqnarray} \tag{ 7 }$$

where ρ(r) is the angle-averaged density as a function of radius. The obtained 〈ℓ_peak〉 values at the last simulation time are 3.75 and 3.56 for models 22L and 27L_A, respectively. Using these numbers, we touch on in the next section how much the difference in the turbulent modes (ℓ_peak), as well as the associated Mach number, could possible affect the explodability of these progenitors.

Bearing in mind an implementation of the precollapse inhomogeneities (obtained in this work) in future CC SN simulations, we present an analysis of the SSH mode of turbulent Mach number. Following Chatzopoulos et al. (2014), we evaluate the SSH decomposition of the radial Mach number (Ma_r) distribution in the burning shells of models 22L and 27L_A, as well as model 25M in our previous study (Yoshida et al. 2019).

Following the procedure in the Appendix of Chatzopoulos et al. (2014), one first needs to determine the inner (r_in) and outer (r_out) boundaries of the region of interest. To determine the inner and outer boundaries, we first determine the radius of the maximum turbulent Mach number ${r}_{\mathrm{Ma},\max }$ in the Si/O-rich layer. Then, the inner boundary is determined as the location where the turbulent Mach number becomes a minimum when we go inward from ${r}_{\mathrm{Ma},\max }$. The outer boundary is determined as the location where the turbulent Mach number becomes less than $1/3\langle {\mathrm{Ma}}^{2}{\rangle }_{\max }^{1/2}$ when we go outward from ${r}_{\mathrm{Ma},\max }$. There is no local minimum around the outer boundary of the Si/O-rich layer. The Mach number is $\lesssim 0.1\langle {\mathrm{Ma}}^{2}{\rangle }_{\max }^{1/2}$ even outside the Si/O-rich layer. Hence, we set the location of the outer boundary using the ratio to the maximum Mach number. The radii of the boundaries are listed in Table 1. We also set the shortest physical length scale λ_r for these models in Table 1. We set this scale to obtain the reconstruction factor, described later, below 0.1. With the given boundaries and shortest physical length, the maximum numbers of the n-modes, n_end, and ℓ-modes, ℓ_end, in the SSH power spectra are determined as

$$\begin{eqnarray}&&{n}_{\mathrm{end}}=\frac{2({r}_{\mathrm{out}}-{r}_{\mathrm{in}})}{{\lambda }_{r}}\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{{\ell }}_{\mathrm{end}}=\displaystyle \frac{\pi ({r}_{\mathrm{out}}+{r}_{\mathrm{in}})}{2{\lambda }_{r}}.\end{eqnarray} \tag{ 9 }$$

The values of n_end and ℓ_end are also listed in Table 1. Then, we evaluate the set of the components of the SSH decomposition a_nℓm with a mode n, ℓ, m from the radial Mach number distribution Ma_r(r, θ, ϕ) in the range between r_in and r_out. Details for the SSH decomposition are written in Appendix A. The SSH decomposition components a_nℓm of models 25M, 22L, and 27L_A are listed in Tables 2–4, respectively, in Appendix A. The machine-readable tables are also available.

Table 1.  Key Parameters for the SSH Decomposition of the Turbulent Mach Number

Model	rin	rout	λr	nend	ℓend	$\langle {\mathrm{Mach}}^{2}{\rangle }_{\max }^{1/2}$	$r(\langle {\mathrm{Ma}}^{2}{\rangle }_{\max }^{1/2})$	ℓave	ΔLcrit/Lcrit	$r{{\prime} }_{\mathrm{in}}$	$r{{\prime} }_{\mathrm{out}}$
	(108 cm)	(108 cm)	(108 cm)				(108 cm)		(%)	(108 cm)	(108 cm)
25M	2.30	12.99	0.5	42	48	0.172	5.48	4.99	8.1	2.49	11.40
22L	2.13	12.99	0.5	43	47	0.158	3.93	5.08	7.3	2.14	11.64
27LA	5.45	63.25	2.5	46	43	0.103	37.6	4.22	5.7	5.29	55.24

Download table as:  ASCIITypeset image

When one applies the spatial distribution of the radial Mach number to a 1D stellar evolution model having an Si/O-rich layer in the range between ${r}_{\mathrm{in}}^{{\prime} }$ and ${r}_{\mathrm{out}}^{{\prime} }$, one can calculate it using the modified SSH decomposition components $a{{\prime} }_{n{\ell }m}$ (see Equation (A10) in Appendix A). In Figure 18, we compare the original radial Mach number distribution (left panel) with the one that is reconstructed by the SSH method (right panel) for models 22L and 27L_A. Note that the radial range of the Si/O layer of the reconstructed one is fitted to the 1D model, which is listed in Table 1. As seen, we obtain a nice match of a global feature. To quantify the difference, we estimate the reconstruction factor, which we estimate in the same range as

$$\begin{eqnarray}&&{f}_{\mathrm{rec}}=\sqrt{\displaystyle \frac{\int {\left({\mathrm{Ma}}_{r}^{\prime} ({\boldsymbol{x}})-{\mathrm{Ma}}_{r}({\boldsymbol{x}})\right)}^{2}{dV}}{\int {\mathrm{Ma}}_{r}{\left({\boldsymbol{x}}\right)}^{2}{dV}}},\end{eqnarray} \tag{ 10 }$$

where Ma_r(x) and ${\mathrm{Ma}}_{r}^{\prime} ({\boldsymbol{x}})$ are the original and reconstructed radial Mach number, respectively (e.g., Chatzopoulos et al. 2014). The values of the reconstruction factor are 0.049, 0.056, and 0.065 for models 25M, 22L, and 27L_A, respectively. Note that in the volume integration in Equation (10) we do not include the regions of five meshes from the inner and outer boundaries. This is because the original Mach number has nonzero values at the boundary, whereas the reconstructed Mach number is set to be zero there by the imposed boundary conditions. Some differences between the reconstructed Mach number and the original one with ≳0.1 appear in a small region close to the boundary along the z-axis. Given this, we consider that the 5%–7% level of the mismatch is not a big concern.

Zoom In Zoom Out Reset image size

Finally, we give an exploratory discussion of how the inhomogeneities and turbulence obtained in this work may impact the perturbation-aided explosion. For this purpose, we estimate the reduction rate of the critical neutrino luminosity required for the onset of neutrino-driven explosion (Equation of (9) in Collins et al. 2018), namely, ΔL_crit/L_crit ∼ 2.34 Ma/ℓ, with Ma and ℓ denoting the Mach number and the mode in the burning shells, respectively. Here we take ℓ from ℓ_ave of Equation (7) and Ma from the maximum Mach number obtained in the computational domain (see Table 1). From this optimistic and crude estimate, we can expect ∼6%–8% levels of reduction in the critical neutrino luminosity. As already discussed in Müller et al. (2016) and Collins et al. (2018), the expected reduction of ∼5% could still make the turbulent perturbation in the burning shells one of several key ingredients for robust explosions. To put the final word, one needs to perform 3D SN simulations using 3D progenitors (Couch & Ott 2015; Müller et al. 2018; Bollig et al. 2020), which we leave for future study.