One-, Two-, and Three-dimensional Simulations of Oxygen-shell Burning Just before the Core Collapse of Massive Stars

We calculate the 1D evolution of massive solar-metallicity stars with ZAMS masses between 9 and 40 M_⊙ up to the onset of collapse of the iron core. The calculations are performed using an up-to-date version of the 1D stellar evolution code, HOngo Stellar Hydrodynamics Investigator⁷ (HOSHI) code (e.g., Takahashi et al. 2016, 2018). In the code, a 300 species nuclear reaction network⁸ is included, the rates of which are taken from JINA REACLIB v1 (Cyburt et al. 2010) except for ¹²C(α, γ)¹⁶O (see Takahashi et al. 2016 for details).

The mass of the helium (He), carbon–oxygen (CO), and iron (Fe) cores as well as the advanced stage evolution depend on the treatment of convection. We use the Ledoux criterion for convective instability. Inside the convective region, we treat the chemical mixing by means of the MLT using the diffusion coefficients as described in Takahashi et al. (2018).

In order to take into account the chemical mixing by convective overshoot, an exponentially decaying coefficient,

$$\begin{eqnarray}&&{D}_{\mathrm{cv}}^{\mathrm{ov}}={D}_{\mathrm{cv},0}\exp \left(-2\displaystyle \frac{{\rm{\Delta }}r}{{f}_{\mathrm{ov}}{H}_{P0}}\right),\end{eqnarray} \tag{ 1 }$$

is included, where D_cv,0, Δr, f_ov, and H_P0 are the diffusion coefficients at the convective boundary, the distance from the convective boundary, the overshoot parameter, and the pressure scale height at the convective boundary, respectively (e.g., Takahashi et al. 2016). We consider the following four models to see the impacts of different overshoot parameters (f_ov) on the 1D stellar evolution. First, we consider two cases during the H- and He-core burning phases (f_ov = 0.03 or 0.01; see Table 1). The former and the latter values are determined based on the calibrations to early B-type stars in the Large Magellanic Cloud (LMC; Brott et al. 2011) and the main-sequence width observed for AB stars in open clusters of the Milky Way Galaxy (Maeder & Meynet 1989; Ekström et al. 2012), respectively. We name the former as model "L" after the LMC and the latter as model "M" after the Milky Way Galaxy. Similarly, in order to investigate the impact in more advanced stages, we test two different overshooting parameters of f_ov,A = 0 or 0.002 (see Table 1) for the advanced stages including the core carbon-burning phase. The convective overshoot during advanced stages is considered for both core and shell convective regions. The subscript "A" is added to the models with f_ov,A = 0.002. When a star model has a ZAMS mass of x M_⊙ and belongs to model L(_A) or M(_A), we set the model name to be xL(_A) or xM(_A). We also name the set of models contained within models xL(_A) or xM(_A) as Set L(_A) or Set M(_A).

We note that the stellar evolution and the final structure also depend on the metallicity and rotation. However, the main purpose of this study is to investigate precollapse inhomogeneities for canonical CCSNe with solar metallicity. We leave the investigation of the metallicity and rotation dependence to a future study.

The mass-loss rate at different evolution stages is important for determining the final mass and the He- and CO-core masses for high-mass stars. We adopt the Vink et al. (2001) mass-loss rate of a main-sequence star when the effective temperature is higher than log T_eff = 4.05, and the surface hydrogen mass fraction X_H is higher than or equal to 0.3. The mass-loss rate of Nugis & Lamers (2000) is adopted for Wolf–Rayet (WR) stars where the effective temperature is higher than log T_eff = 4.05 and the surface hydrogen mass fraction X_H is lower than 0.3. When the effective temperature is lower than log T_eff = 3.90, we adopt de Jager et al. (1988) mass-loss rate.

We compute 2D and 3D models with our hydrodynamics code, the 3 Dimensional nuclear hydrodynamic simulation code for Stellar EVolution (3DnSEV), which is a branch of the 3DnSNe code (see Nakamura et al. 2016; Takiwaki et al. 2016; Sasaki et al. 2017; Kotake et al. 2018 for recent code development). Similar to the base code, 3DnSEV solves Newtonian hydrodynamics equations using spherical polar coordinates as follows:

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{v}}\right)=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\partial }_{t}\left(\rho {\boldsymbol{v}}\right)+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{v}}{\boldsymbol{v}}+p{\boldsymbol{\delta }}\right)=-\rho {\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}{e}_{\mathrm{tot}}+{\rm{\nabla }}\cdot \left[({e}_{\mathrm{tot}}+p){\boldsymbol{v}}\right]=-\rho {\boldsymbol{v}}\cdot {\rm{\nabla }}{\rm{\Phi }}\\ \quad \ +\ \rho {\epsilon }_{\mathrm{burn}}+C,\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\partial }_{t}\rho {X}_{i}+{\rm{\nabla }}\cdot \left(\rho {X}_{i}{\boldsymbol{v}}\right)={\gamma }_{\mathrm{burn}},\end{eqnarray} \tag{ 5 }$$

where $\rho ,{\boldsymbol{v}},p,{e}_{\mathrm{tot}},{\rm{\Phi }}$ are the density, velocity, pressure, total energy density (sum of the internal energy and kinetic energy), and gravitational potential, respectively. X_i denotes the mass fraction of the ith isotopes and _burn is the energy generated by the change in composition, γ_burn, due to nuclear burning. C is the energy loss by neutrino emission. The subgrid-scale physics is handled by implicit numerical diffusion instead of solving filtered hydrodynamic equations and by creating a subgrid model for the dissipation of kinetic energy as a large-eddy simulation. A piecewise linear method with the geometric correction of the spherical coordinates is used to reconstruct variables at the cell edge, where a modified van Leer limiter is employed to satisfy the condition of total variation diminishing (TVD; Mignone 2014). The numerical flux is basically calculated by an HLLC solver (Toro et al. 1994). For the numerical flux of isotopes, the consistent multifluid advection method of Plewa & Müller (1999) is used. The models are computed on a spherical polar coordinate grid with a resolution of n_r × n_θ × n_ϕ = 512 × 64 × 128 (3D) and n_r × n_θ = 512 × 128 (2D) zones. The radial grid is logarithmically spaced and covers the center up to the outer boundary of 10¹⁰ cm. For the polar and azimuthal angles, the grid covers all 4π sr. To focus on the convective activity mainly in the oxygen shell, the inner 10⁸ cm is solved in spherical symmetry. We include self-gravity assuming a spherically symmetric (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.

We use the "Helmholtz" equation of state (EOS; Timmes & Swesty 2000). Neutrino cooling is taken into account (Itoh et al. 1996) as a sink term in the energy equation. A nuclear reaction network of 21 isotopes⁹ (aprox21; Paxton et al. 2011) is implemented, where the inclusion of ⁵⁴Fe, ⁵⁶Fe, and ⁵⁶Cr is crucial for treating low-electron fractions Y_e ≳ 0.43 in the presupernova stage. The network is as large as that of Couch et al. (2015) and a little larger than the 19 isotopes of Müller et al. (2016). When the temperature is higher than 5 × 10⁹ K, the chemical composition is assumed to be in nuclear statistical equilibrium (NSE). To avoid the temperature variations caused by numerical instability, we set an artificial upper bound in our multi-D runs, in such a way that the (absolute) sum of the local energy generation rates by thermonuclear reactions and weak interactions does not exceed 100 times the local neutrino cooling rate. To correctly treat 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. Because the NSE region appears mainly in the Fe core may not significantly affect the convection in the O layer on which we focus.

About 100 s before the onset of collapse, the 1D evolution models are mapped to our multi-D hydrodynamics code, and we follow the 2D and 3D evolution for ∼100 s until the onset of collapse. When we start the multi-D 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 the 2D/3D runs when the central temperature exceeds 9 × 10⁹ K, because the core is dynamically collapsing at this time.

In total, 100 stellar evolution models are calculated in 1D. Four sets of models are constructed, to which different overshoot parameters are applied (Table 1). Each set consists of 25 models to cover the initial mass range of 9–40 M_⊙. These 1D models are evolved from the ZAMS stage up to the onset of collapse, which is determined using the threshold central temperature of T_C ∼ 10^9.9 K. Details of our 1D evolution models (e.g., comparison with reference stellar evolution codes) are given in the Appendix.

Figure 1 shows the total mass (red) and masses of the He (blue), CO (green), and Fe (magenta and cyan) cores at the onset of collapse as a function of the ZAMS mass (M_ZAMS). The top panel shows results for Sets L and L_A, while results for Sets M and M_A are shown in the bottom panel. We note that Sets L and L_A result in very similar total, He-core, and CO-core masses, and differ only in the Fe-core mass. This is because the He- and CO-core masses are mostly determined by the size of the core convection during the H- and He-burning phases, respectively, and are largely independent of the overshoot during more advanced stages. The same is true for Sets M and M_A. The He-core mass is defined as the largest enclosed mass with hydrogen mass fraction less than 10⁻³. Similarly, the CO-core mass is defined as the largest enclosed mass with He mass fraction less than 0.1, and the Fe-core mass is defined as the largest enclosed mass with the sum of the mass fractions of Z ≥ 21 elements larger than 0.5.

Zoom In Zoom Out Reset image size

The total mass at the collapse is determined by the mass-loss history. Because the mass loss is relatively weak, the total mass monotonically increases with the ZAMS mass for models below M_ZAMS ≲ 20–25 M_⊙ in both Sets L and M. The mass-loss rate increases with increasing luminosity, and thus with increasing ZAMS mass. The increasing mass-loss rate explains the flat and even decreasing trends seen in the 20–30 M_⊙ models in Set L. At the same time, models in the same mass range show a stochastic trend for Set M. This is caused by the bistability jump of the mass-loss rate, which results from the discontinuous rate increase along the decreasing effective temperature (e.g., Vink et al. 2000). For more massive models above M_ZAMS ≳ 30 M_⊙, the mass-loss rate becomes so efficient to remove most of the H envelope during the He-burning phase. Therefore, the total masses of these models coincide with their He-core masses. This is why the total mass again shows a monotonic increase in this massive end of the ZAMS mass range. The most massive models (the 32, 35, and 38 M_⊙ models of Set L and the 35, 38, and 40 M_⊙ models in Set M) finally retain only a small amount of hydrogen of 0.26–0.29 M_⊙ in their envelopes, which will correspond to these being observed as late-type WN stars (Crowther 2007).

As an exception, model 40L (M_ZAMS = 40 M_⊙ model in Set L) has lost not only the entire H envelope but also most of the He layer. This is due to the even stronger WR wind mass loss during the helium- and carbon-burning phases. The He mass remaining on the surface is 0.24 M_⊙. We apply the mass-loss rate of Nugis & Lamers (2000) for the H-deficient stars. However, there is a large uncertainty in the estimation of the WR wind mass-loss rate. Among H-deficient stars especially, the mass-loss rate of He-deficient WC stars can be larger than the rate in Nugis & Lamers (2000) by a factor of ∼10 (Yoon 2017). The remaining He mass can be even less if we consider more efficient WR wind mass loss. Therefore, we expect that the star will most probably be observed as a He-deficient WC star, and moreover, it will be observed as a Type Ic supernova when this star explodes.

The features of the distributions of the final stellar mass and the He- and CO-core masses as a function of ZAMS mass are also seen in the results obtained in previous works using KEPLER code (Woosley & Heger 2007; e.g., Figure 4 of Ebinger et al. 2019 for a concise summary).

The He- and CO-core masses monotonically increase with ZAMS mass except for model 40L_A, which is affected by the strong WR wind during the He-burning phase. The mass of the helium layer, which is shown as the difference between the He-core and the CO-core masses, also increases with the ZAMS mass. As mentioned earlier, the He- and CO-core masses are insensitive to the overshoot parameter after the He burning, f_ov,A. Thus, the difference in the CO-core mass between Sets L_A and L (and similarly between Sets M_A and M) is less than 0.7%. Furthermore, the difference in the He-core mass is less than 0.1%. Note that model 18M_A exceptionally forms a CO-core mass about 3% larger than model 18M. This results from the emergence of a narrow convection in the outer layer of the CO core, in which a small amount of He is contained. Because this narrow convection is activated only after the oxygen-core depletion, the structure outside the He layer of model 18M_A is mostly the same as that of model 18M, like other models.

We will compare Set L with Set M. Below M_ZAMS ≲ 25 M_⊙, the final mass is not that sensitive to the overshoot parameter for the H- and He-core convection. On the other hand, models above 25 M_⊙ show a scatter within a factor of ∼0.3. Set L tends to show larger He- and CO-core masses than Set M. This is simply due to the larger overshoot parameter applied during the H- and He-burning phases. For most of the models, the ratio of the He-core masses is about a factor of 1.2–1.3, and the CO-core mass ratio is somewhat larger than the He-core mass ratio. As an exception, the He-core mass ratio reaches 2.16 for the 9 M_⊙ models. This is due to the merging of He layer to the H envelope by the second dredge-up. A larger overshoot for Set L brings about more effective convective mixing to make more massive He and CO cores. Small core-mass ratios for 40 M_⊙ models are due to the strong mass loss occurring in model 40L.

To select models that will show strong convective activities in the SiO-rich layer, we utilize two measures, the SiO-coexistence parameters of f_M,SiO and f_V,SiO. Both of them are defined based on mass fraction distributions of ¹⁶O and/or ²⁸Si. The f_M,SiO is a product of the mass fractions of ¹⁶O and ²⁸Si weighted by the enclosed mass between 10⁸ cm and 10⁹ cm:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{M,\mathrm{SiO}} & = & {c}_{M}{\displaystyle \int }_{1}^{10}X{(}^{16}{\rm{O}})X{(}^{28}\mathrm{Si})\\ & & \times \ {\rm{\Theta }}(X{(}^{16}{\rm{O}})-0.1){\rm{\Theta }}(X{(}^{28}\mathrm{Si})-0.1)\rho {r}_{8}^{2}d({r}_{8}),\end{array}\end{eqnarray} \tag{ 6 }$$

where c_M is a scaling coefficient, X(^AZ) is the mass fraction of isotope ^AZ, and Θ(x) is the step function, which satisfies Θ(x) = 1 for x ≥ 0 and 0 for x < 0, ρ is the density, and r₈ is the radius in units of 10⁸ cm. Therefore, the value becomes large in a model that has a layer mainly composed of both oxygen and silicon. Such a layer would be the most preferable site to host strong turbulence powered by oxygen-shell burning. This definition of f_M,SiO has an uncertainty on how the local density strength of the turbulence depends. Hence, we also test another indicator, f_V,SiO, in which the product of the mass fraction is weighted not by the enclosed mass but by the enclosed volume instead:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{V,\mathrm{SiO}} & = & {c}_{V}{\displaystyle \int }_{1}^{10}X{(}^{16}{\rm{O}})X{(}^{28}\mathrm{Si})\\ & & \times \ {\rm{\Theta }}(X{(}^{16}{\rm{O}})-0.1){\rm{\Theta }}(X{(}^{28}\mathrm{Si})-0.1){r}_{8}^{2}d({r}_{8}),\end{array}\end{eqnarray} \tag{ 7 }$$

where c_V is a scaling coefficient. The scaling coefficients are arbitrarily chosen. We calculate these two measures at every time step from 120 to 80 s before the last step of the calculations to see the characteristics at times close to the onset of the multidimensional simulations.

The result is shown in Figure 2, wherein c_M = 3.2 × 10⁻¹⁰ and c_V = 0.025 are applied. We do not see clear dependencies among different treatments of the overshoot. Some models in the ZAMS mass range ≳22 M_⊙ show large (≳0.6) f_M,SiO values. In the volume-weighted case, the ZAMS mass range showing the models having large (≳0.9) f_V,SiO values is 13–28 M_⊙. From this result, we selected 11 models, in which either f_M,SiO or f_V,SiO, or possibly both of them, shows a large value. Models showing the seven highest f_M,SiO values are models 28M, 23L_A, 25M, 28L_A, 27L_A, 27M, and 22L. Models showing the six highest f_V,SiO values are models 13L_A, 28M, 21M_A, 25M, 16M_A, and 18M_A. Among them, models 25M and 28M show large values for both f_M,SiO and f_V,SiO. The actual values of the parameters are shown in the second and third columns of Table 2.

Zoom In Zoom Out Reset image size

Table 2.  2D Model Properties and SiO-coexistence Parameters

Model	fM,SiO	fV,SiO	$\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$	r($\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$)	Layer	ℓmax	dc/HP
				(108 cm)
Low Ma
13LA	0.27–0.28	0.95–0.96	0.018	11.6	O/Si	12	6.22
16MA	0.24–0.24	0.90–0.91	0.015	3.9	O/Si	4	3.20
18MA	0.57–0.58	0.91–0.91	0.131	3.1	Si/O	14	1.06
21MA	0.47–0.47	0.91–0.95	0.134	3.0	Si/O	8	4.42
23LA	0.75–0.80	0.78–0.80	0.069	11.5	O/Si	4	5.20
High Ma
22L	0.57–0.61	0.77–0.82	0.108	9.4	Si/O	2	2.50
25M	0.75–0.79	0.91–0.94	0.160	5.8	Si/O	3	3.65
27LA	0.59–0.66	0.76–0.76	0.179	45.0	O/Si	2	4.56
27M	0.58–0.65	0.37–0.40	0.134	4.7	Si/O	10	2.44
28LA	0.60–0.68	0.37–0.42	0.117	5.3	Si/O	8	1.81
28M	0.83–0.90	0.90–0.95	0.369	14.6	O/Si	2	4.08

Note. SiO-coexistence parameters f_M,SiO and f_V,SiO are obtained from the result of the 1D evolution simulations. $\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$ represents the maximum turbulent Mach number obtained at a radius of r($\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$) at the end of the 2D simulations. "Layer" represents the composition of the convective region. ℓ_max represents the ℓ value where ${c}_{{\ell }}^{2}$ has a peak (see Equation (9)). d_c/H_P represents the width of the convective region normalized by the local scale height. These quantities are all estimated at the last step of the simulations. See the text for a more detailed definition.

Download table as:  ASCIITypeset image

For later convenience, we separate the SiO-rich layer into the "Si/O" layer and the "O/Si" layer. The "Si/O" layer has a larger Si mass fraction than O mass fraction in the layer, i.e., X(²⁸Si) ≥ X(¹⁶O) and X(¹⁶O) ≥ 0.1, whereas the "O/Si" layer has the relation 0.1 ≤ X(²⁸Si) < X(¹⁶O). Then, we can classify these 11 models into two groups having different structures of the SiO-rich layer. One group has an extended O/Si layer instead of the O/Ne layer above the Si/Fe layer. The other group has a Si/O layer between the inner Si/Fe layer and the outer O/Ne layer. The former consists of models 13L_A, 16M_A, 18M_A, 21M_A, 23L_A, 27L_A, and 28M. The top panel of Figure 3 shows the mass fraction distribution of model 13L_A as a function of radius. The radius of the outer boundary of the O/Si layer is ∼3 × 10⁹ cm. This layer was originally formed as an O/Ne layer. Neon burning has started after the core silicon-burning phase, transforming neon into oxygen and silicon. In the case of models 18M_A, 21M_A and 23L_A, a thin Si/O layer exists between the Si/Fe layer and the O/Si layer with a width of less than 3 × 10⁸ cm.

Zoom In Zoom Out Reset image size

Models in the latter group have a layered structure, in which the innermost Fe core is surrounded by the Si/Fe, Si/O, and O/Ne layers. Models 22L, 25M, 27M, and 28L_A comprise this group. The bottom panel of Figure 3 shows the mass fraction distribution of model 25M, an example of the latter group. The model has a Si/O layer between 3 × 10⁸ cm and 1.1 × 10⁹ cm. For other models in this group, the width of the Si/O layer is typically several times 10⁸ cm. We show the mass fraction distributions as a function of radius for models other than 13L_A and 25M in Figure 17.

In order to investigate the convective activities in a multidimensional space, we perform 2D hydrodynamics simulations of oxygen-shell burning. In the previous subsection, we picked up 11 models that show large SiO-coexistence parameters, f_M,SiO and/or f_V,SiO. Profiles of these models at ∼100 s before the end of the 1D calculations are taken as the initial conditions. The 2D calculations are continued until the central temperature reaches 9 × 10⁹ K, by which point the stars have started runaway collapse due to the gravitational instability.

Following Müller et al. (2016), we evaluate the angle-averaged turbulent Mach number as an indicator of the turbulence strength,

$$\begin{eqnarray}&&\langle {{Ma}}^{2}{\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{ 8 }$$

where ρ is the density; v_r, v_θ, and v_ϕ are the radial, tangential, and azimuthal velocities; $\langle {v}_{r}\rangle$ is the angle-averaged radial velocity; c_s is the sound velocity; and Ω is the solid angle. The maxima of $\langle {{Ma}}^{2}{\rangle }^{1/2}(r)$ evaluated at the end of the simulations $\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$ are shown in the fourth column of Table 2, and $r(\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2})$ represents the radii where $\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$ are obtained.

Based on the Mach number, we divide our 2D models into two groups, either showing "low Ma" or "high Ma." The criterion of high Ma is set as $\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}\geqslant 0.1$, because the turbulence with such a high Mach number potentially fosters the perturbation-aided explosion (Müller & Janka 2015; Müller et al. 2016). It is noted that models 18M_A and 21M_A are exceptionally classified into low Ma despite their "large" Mach numbers. We discuss this later in this section. The column "Layer" in Table 2 represents the dominant chemical composition in the convective layer, i.e., the Si/O layer or the O/Si layer. We also show the mass fraction distribution of models 13L_A and 25M in Figure 3. For other models, the mass fraction distributions are shown in Figure 17 in the Appendix. See the last part of the previous subsection for the definition of the layer.

Time evolution of convective motion—In Figure 4, we show the time evolution of the turbulent Mach number and the Si mass fraction for representative models from low Ma (13L_A, top) and high Ma (25M, middle, and 27L_A, bottom). The color visualizes the angle-averaged turbulent Mach number $\langle {{Ma}}^{2}{\rangle }^{1/2}$ (left) and the ²⁸Si mass fraction X(Si) (right). Note that the outer radial frame of the panels for model 27L_A is set to 8 × 10⁹ cm in order to show how the outer edge of the convective region keeps moving outward and reaches this radius at the end.

Zoom In Zoom Out Reset image size

The model 13L_A has no Si/O layer. This star has an Fe core at the central region of R ≲ 2 × 10⁸ cm, which is surrounded by the convective Si/Fe layer (R ∼ 2–4 × 10⁸ cm) and the convective O/Si layer (R ≳ 4 × 10⁸ cm). The Si mass fraction at the Si/Fe layer is ∼0.5 (see the top panel of Figure 3). As shown in the top-right panel of Figure 4, the Si mass fraction is small compared to that of the 25M model that is shown in the middle right panel of Figure 4. Reflecting this structure, the turbulent Mach number is lower than 0.1 in the inner Si/Fe layer and in the outer O/Si layer throughout the simulation (see the top-left panel of Figure 4). However, oxygen burning slightly enhances the ²⁸Si mass fraction in the base region of the O/Si layer of ∼4–8 × 10⁸ cm. Note that $\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$ in Table 2 is estimated at the end of the simulations, and it does not refer to the peak seen at ∼4 × 10⁸ cm at ∼30 s.

Readers may be confused by the models 18M_A and 21M_A because they have the Si/O layer in Table 2 but they are classified into the low-Ma group. Actually, models 18M_A and 21M_A have a thin Si/O layer and $\langle {{Ma}}^{2}{\rangle }^{1/2}$ ∼ 0.13 in the Si/O layer, but only after the last ∼10 s of the simulations. This is because the turbulence is triggered by gravitational contraction, which amplifies the temperature at the bottom of the Si/O layer, enhancing the oxygen-burning rate. The turbulence powered by gravitational contraction has too short a time to form an extended convective region, which is in contrast to the shell convection powered by hydrostatic burning. This is why we selected these models as members of low Ma.

Models 22L, 25M, 27L_A, 27M, 28L_A, and 28M are categorized into models with high Ma. Convective motion with such a strong turbulence develops by oxygen burning in the Si/O layer in these models. We pick out two models in which it is easy to explain the typical dynamics of the convection.

Model 25M consists of the central Fe core (R ≲ 2 × 10⁸ cm), the Si/Fe layer (R ∼ 2–3 × 10⁸ cm), the Si/O layer (R ∼ 3–10 × 10⁸ cm), and the O/Ne layer (R ≳ 10 × 10⁸ cm). It is noteworthy that, despite the Si/Fe layer seeming to have a homogeneous chemical composition (see the lower panel of Figure 3), the outer part of R ∼ 2.5–3 × 10⁸ cm is actually composed of a small amount of oxygen with X(¹⁶O) < 0.01. The oxygen-free region of R ∼ 2–2.5 × 10⁸ cm becomes convective within a short timescale of ∼10 s from the start of the simulation, though the shell convection is not extended further. At the bottom of the outer Si/O layer, hydrostatic oxygen-shell burning takes place. In this case, nuclear burning drives high turbulent velocity with $\langle {{Ma}}^{2}{\rangle }^{1/2}\gt 0.1$, which is sustained for 20–110 s (see the middle-left panel of Figure 4). Accordingly, turbulent mixing homogenizes the ²⁸Si mass fraction in the region of R = 3–10 × 10⁸ cm. Furthermore, oxygen burning also takes place in the oxygen-containing outer region of the Si/Fe layer. Despite the large mean molecular weight, the heating due to oxygen burning is strong enough to lift the silicon-rich material up into the surrounding Si/O layer. As a result, the silicon mass fraction in the Si/O layer is significantly enhanced. This silicon enhancement repetitively takes place at ∼50 and 80 s, which accompanies the enhancement of the convective Mach number as well. It seems that repetitive mixing follows the oscillation of the outer edge of the Si/Fe convection. This will be because, with the small oxygen mass fraction, the temperature fluctuation originally caused by the oscillation is enhanced by the O burning, resulting in a large density fluctuation that triggers convection in the Si/O layer. Indeed, the temperature rise at the bottom of the Si/O layer where the O mass fraction is ∼0.08–0.1 reaches 8% at maximum, which is much higher than the temperature change solely due to the oscillation, which is less than ∼1%, measured in the outer region of the Si/F and Si/O layers.

In model 27L_A, the turbulent activity in the O/Si layer starts to increase at ∼45 s (see the bottom-left panel of Figure 4). During the simulation time of ≳200 s, the high Mach number region extends outward, finally reaching ∼6 × 10⁸ cm, which roughly corresponds to the composition jump between the CO-rich and the He-rich layers. At the same time, the turbulent Mach number grows with time in almost the entire region in the convective layer. As a consequence, the turbulent Mach number exceeds ∼0.15 in the wide outer region of R ∼ 30–50 × 10⁸ cm at the end of the simulation. Initially, the ²⁸Si mass fraction decreases with radius in the O/Si layer (see the mass fraction distribution of model 27L_A in Figure 17). The convection powered by the oxygen burning mixes material in the slightly silicon-enriched region, which is initially located below ∼8 × 10⁸ cm, into the outer, slightly silicon-poor region after ∼70 s. However, the convective mixing in the 2D simulation is still not efficient enough to achieve the homogeneous chemical distribution. This is due to the limitation of the calculation time, because the total time of this simulation covers only about one convection-turnover time. Model 28M shows similar convective properties to model 27L_A.

High/low Mach number and the chemical distribution—Models 13L_A, 16M_A, 18M_A, 21M_A, and 23L_A belong to low Ma, i.e., they do not show strong turbulence in their convective regions during the simulations. These low-Ma models have characteristic chemical composition profiles. A main characteristic is no or thin Si/O layer. Models 13L_A and 16M_A do not have a Si/O layer and have an extended O/Si layer on the Si layer (see the top panel of Figure 3 and top-left panel of Figure 17, respectively). Model 23L_A does not have a Si/O layer at the beginning of the 2D simulation. The turbulent Mach number of these models in the O/Si layer is low. Models 18M_A and 21M_A have a Si/O layer at the beginning of the 2D calculations but the width is less than ∼1 × 10⁸ cm. Although the turbulent Mach number exceeds 0.1 at the bottom of the Si/O layer for a few seconds before the termination of the simulations, the turbulence in this layer does not develop before this time.

Models 22L, 25M, 27L_A, 27M, 28L_A, and 28M belong to high Ma. The main characteristic of the chemical profiles is an extended Si/O layer. Models 22L, 25M, and 28M have a Si/O layer with the width of ∼8 × 10⁸ cm. Note that the Si/O layer of model 28M has merged to the O/Ne layer before the end of the simulation. Models 27M and 28L_A also have a Si/O layer, although their width is thinner than that of the three models.

Model 27L_A is an exception. This model does not have a Si/O layer but the chemical composition profile is similar to the last step of model 28M.

We briefly discuss the relation to the SiO-coexistence parameters. All models in high Ma are selected using a large f_M,SiO value. Although models 25M and 28M are also selected using a large f_V,SiO value, they also have a large f_M,SiO. The models having a large f_M,SiO rather than a large f_V,SiO are associated with high Ma.

2D distribution—Figure 5 shows the 2D distributions of the radial turbulent Mach number, Ma_r = ${v}_{r}/{c}_{s}-\langle {v}_{r}\rangle /\langle {c}_{s}\rangle$, and the ²⁸Si mass fraction taken at the last step of the simulations for models 13L_A, 25M, and 27L_A. The turbulent Mach number of model 13L_A (top panels in Figure 5) develops only within the level of Ma_r ∼ 0.01. The spherical boundary is clearly observed at r ∼ 3 × 10⁸ cm, where the turbulent Mach number becomes almost zero. Inside the boundary, convection is developed in the Si/Fe layer. A more extended but even weaker turbulent motion is also developed in the O/Si layer above the boundary. The outer boundary of the O/Si convection may be defined at ∼6 × 10⁸ cm, but there is only a diffuse Ma_r ∼ 0 region in this case. A thin and nearly spherical band with X(²⁸Si) ∼ 0.2 surrounding the inner boundary between the Si/Fe and the O/Si layers exists. As a result of the low-velocity turbulence in the O/Si layer, this silicon-rich material is slowly mixed into the inner region of the O/Si layer at R ≲ 6 × 10⁸ cm.

Zoom In Zoom Out Reset image size

Model 25M (middle panels in Figure 5) develops convective motion with a high turbulent Mach number in the Si/O layer ranging from R ∼ 3 × 10⁸ cm to R ∼ 10 × 10⁸ cm. At the end of the simulation, outflows stream in three directions: the northern pole direction, ∼45° from the polar axis, and ∼135° from the polar axis, and inflows are sandwiched by the outflows. These convective flows have turbulent Mach numbers larger than ∼0.1. The ²⁸Si mass fraction is roughly homogenized inside the convective region, having the value of 0.3–0.4, though some fluctuations are observed especially near the outer boundary.

Model 27L_A (bottom panels in Figure 5) has an extended O/Si layer distributed from R ∼ 5 × 10⁸ to 50 × 10⁸ cm. A large-scale convective motion is developed in this layer; a broad conical outflow with the opening angle of ∼45° is formed in both polar regions, and between them, a thick inflow is formed around the equatorial plane. The convective Mach number reaches ∼0.12. The large-scale outflow mixes the silicon-rich material into the O/Si layer. The silicon mass fraction in most parts of this layer is initially ∼0.1, while the outflow has a higher fraction of ∼0.16.

Width of the convective region—We briefly discuss the width of the convective region divided by the local scale height. We define a convective region as a region having $\langle {{Ma}}^{2}{\rangle }^{1/2}\gt 1/3\,\times \langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$ including $r(\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2})$. We determined the factor of 1/3 to avoid including the neighboring convective region because we obtain a small turbulent Mach number even at the convection boundary. The width of the convective region is compared with the pressure scale height H_P at $r(\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2})$. This results in 1.8–4.6 for the high-Ma models, which are listed in Table 2. We should note that the above definition does not specify the width of the convective region correctly for low-Ma models. In models 16M_A and 18M_A, the specified region contains the Si layer inside the O/Si and Si/O layer, respectively. In model 21M_A, the calculated region contains part of the Fe core, the Si and Si/O layer, and part of the O/Si layer. The turbulent Mach number at the boundary determined by the abundance distribution is not small enough compared with $\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2}$ to specify the boundary of the convective region for these models.

Typical scale of the convection—In addition to the Mach number, the dominant angular wave number in spherical harmonics also characterizes the convection. This is related to the typical size of the convective flow. This quantity is important because a large-scale convective flow can amplify the explodability of core-collapse supernovae (Müller et al. 2016). The power spectrum of the radial turbulent velocity at $r(\langle {{Ma}}^{2}{\rangle }_{\max }^{1/2})$ is calculated as

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

where Y_ℓm(θ) is the spherical harmonics function of degree ℓ and order m. ℓ_max in the table represents the ℓ value at which ${c}_{{\ell }}^{2}$ has a peak.

Figure 6 shows the power spectrum ${c}_{{\ell }}^{2}$ at three different times for models 13L_A (top), 25M (middle), and 27L_A (bottom). For model 13L_A, ${c}_{{\ell }}^{2}$ has a maximum at ℓ = 12, but the spectrum is rather flat and less energetic. The radius of the highest Mach number in the O/Si layer is 1.16 × 10⁹ cm. Although convective mixing occurs in the inner region of R ≲ 6 × 10⁸ cm, the turbulent velocity is lower than that in the outer region at the last step. Large-scale convection in the O/Si layer is not developed probably because of small turbulent motion. For model 16M_A, ℓ_max is equal to 4 and the trend of the power spectrum is similar to model 13L_A.

Zoom In Zoom Out Reset image size

For the other low-Ma models, models 18M_A and 21M_A show large ℓ_max (see Table 2), and they have a thin Si/O layer. This trend is roughly consistent with the analysis of the convective eddy scale relating to the scale of the convective layer and the typical radius of the layer (e.g., Müller et al. 2016). In these models, the SiO-rich layer is thin compared to the radius of the layer.

In models 25M and 27L_A, which are high-Ma models, the power spectrum peaks at ℓ_max = 3 and 2, respectively, and the spectrum decreases with increasing ℓ above that. Models 22L and 28M show a similar power spectrum to model 27L_A. For these models, large-scale convective eddies have developed. On the other hand, models 27M and 28L_A indicate larger ℓ values probably owing to the thin SiO-rich layer. Indeed, the previous three models show a larger width of the convective region normalized by the scale height compared to the latter two models (see Table 2). Note that these models develop shell convection in the Si/Fe layer as well. However, the convective region is always confined inside the layer. This will be because the timescale of the silicon burning is shorter than the convective turnover time, so that the mean molecular weight of the convective blob soon increases, suppressing the convective motion.

From the results shown above, it is discerned that f_M,SiO will be a more suitable measure than f_V,SiO to discriminate a model that develops convection with high turbulent velocity and a small ℓ_max. First, we have shown that high-Ma models are selected based on the high f_M,SiO values. Moreover, models showing small ℓ_max ≤ 3 are all selected based on f_M,SiO (models 22L, 25M, 27L_A, and 28M), and only one of the two models showing ℓ_max = 4 (model 16M_A) is selected based on f_V,SiO.

A 3D hydrodynamics simulation is conducted using the 1D model 25M as the initial condition. The size of the convective region of model 25M would be suitable for investigating the multidimensional effects of the structure of a presupernova star on the supernova explosion.¹⁰ The hydrodynamical evolution is followed for ∼100 s until the central temperature reaches 9 × 10⁹ K, at which point the Fe core is unstable enough to collapse.

Figure 7 shows the time evolution of the ²⁸Si mass fraction distribution of model 25M. The initial distribution of the ²⁸Si mass fraction is spherically symmetric (top-left panel). After the start of the simulation, the convection in the Si/O layer develops from the inner region. We see that Si-enriched plumes go up into the Si/O layer (top-right panel). The convective motion reaches a steady flow by ∼20 s. The inhomogeneous ²⁸Si mass fraction distribution introduced by the convection at 30 s is shown in the middle-left panel. After a while, the turbulent velocity becomes small. We will discuss the mechanism of this weakening later.

At ∼70 s, the Si/O layer gradually contracts, triggering the strong oxygen-shell burning at the bottom of the Si/O layer. This strong burning drives high-velocity turbulence and expands the Si/O layer. We see some ²⁸Si-enriched plumes flow from the inner region of the Si/O layer at 75 s (see red region in the middle right panel). As a result, the high-velocity turbulent flow mixes with the surroundings and increases the Si mass fraction in the whole Si/O layer. The ²⁸Si mass fraction in the Si/O layer slightly increases from about 0.36 at 30 s to 0.38 at 90 s (bottom-left panel). The convective motion in the Si/O layer continues until the last step of the simulation. We see inhomogeneous ²⁸Si mass fraction distribution at the last step (bottom right panel).

In Figure 8, the time evolution of $\langle {{Ma}}^{2}{\rangle }^{1/2}$ and the ²⁸Si mass fraction obtained from the 3D simulation is shown. As shown by the left yellow region in the top panel, the convective motion with $\langle {{Ma}}^{2}{\rangle }^{1/2}\sim 0.1$ is obtained in the Si/O layer at ∼20 s. The bottom panel shows that the ²⁸Si mass fraction is enhanced in this layer by that time. After a while, the Si/O layer slightly expands and the convective motion weakens. The averaged Si mass fraction in the Si/O layer does not change significantly from 20 to 70 s.

Zoom In Zoom Out Reset image size

By ∼70 s, the Fe core contracts and silicon-shell burning as well as oxygen-shell burning is enhanced. Because of this, the turbulent Mach number becomes large not only at R ∼ 3–12 × 10⁸ cm but also at R ∼ 2 × 10⁸ cm. In addition, the Si-rich material at the bottom of the Si/O layer is carried into the Si/O layer through this convective motion, and the ²⁸Si mass fraction in this region increases during ∼70–80 s.

From ∼90 s, the convective motion becomes weak again. The whole Si/O layer gradually contracts toward core-collapse. Meanwhile, the base temperature of the Si/O layer increases, boosting the oxygen-shell burning. This results in the enhancement of the convective motion at the bottom of the Si/O layer. The Si/O layer at the end of the simulation is distributed from 3.0 × 10⁸ cm to 10.5 × 10⁸ cm. The maximum radial convective Mach number is $\langle {{Ma}}^{2}{\rangle }^{1/2}=0.087$, which is obtained at R = 3.6 × 10⁸ cm.

Figure 9 shows the power spectrum of the radial convective velocity ${c}_{{\ell }}^{2}$ at R = 5.8 × 10⁸ cm taken at three different times. The radius, which is located in the middle of the convective Si/O layer, is determined from the radius where the maximum $\langle {{Ma}}^{2}{\rangle }^{1/2}$ is obtained in the 2D simulation. Similar to the 2D case, the power spectrum in the 3D case is calculated as

$$\begin{eqnarray}&&{c}_{{\ell }}^{2}=\displaystyle \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{ 10 }$$

The power spectrum ${c}_{{\ell }}^{2}$ in the 3D simulation peaks at ℓ = 2. This is consistent with the finding in Müller et al. (2016) for the 18 M_⊙ star. The small maximum mode means that the convective motion is dominated by a large-scale flow. The main difference between 2D and 3D power spectra is the weaker turbulence in the 3D simulation. Note that the stronger turbulence in 2D compared to 3D is not surprising because the turbulent energy cascade could occur artificially, as previously identified, from small to large scales along the coordinate symmetry axis. This is mainly due to the reduced degree of freedom in which the energy can dissipate.

Zoom In Zoom Out Reset image size

Finally, we briefly present the result of the neutrino emission at the precollapse stage of the 3D simulation. Using the same method in Yoshida et al. (2016), we calculate the time evolution of the luminosity and the spectrum of neutrinos emitted via the pair-neutrino process for model 25M. In Figure 10, the emission rate and the average energy of neutrinos obtained from the 1D and the 3D simulations are compared.

Zoom In Zoom Out Reset image size

The overall features of the neutrino spectra for the 1D and 3D simulations are in common. The neutrino emission rate and the average temperature increase with time toward the core-collapse after ∼30 s in both simulations. The emission rates of ν_e and ${\bar{\nu }}_{e}$ are larger than those of ν_μ, τ, and ${\bar{\nu }}_{\mu ,\tau }$. The average energy of ${\bar{\nu }}_{e}$ is slightly smaller than that of ${\bar{\nu }}_{\mu ,\tau }$. In the 3D simulation, however, the decrease in the emission rate and the average energy of neutrinos is observed in 70–90 s.

Furthermore, we evaluate the time evolution of the event rate of ${\bar{\nu }}_{e}$ using KamLAND (e.g., Gando et al. 2013), assuming that the initial mass of 25 M_⊙ at the distance of 200 pc is at the ongoing core-collapse phase, emitting neutrinos via the pair-neutrino process.

KamLAND is a one-kiloton size liquid-scintillation-type neutrino detector (see, e.g., Gando et al. 2013). We take the neutrino oscillation into account in a simple manner: the survival probability of ${\bar{\nu }}_{e}$ is set to be 0.675 and 0.024 in the normal and inverted mass ordering, respectively (Yoshida et al. 2016). The live-time-to-runtime ratio and the total detection efficiency are set to be 0.903 and 0.64, respectively (Yoshida et al. 2016). Figure 11 shows the evolution of the neutrino event rate using KamLAND. As an overall feature, the rate increases with time independent of mass ordering. In the 3D simulation, however, the decrease in the event rate due to the oxygen-shell burning at the bottom of the Si/O layer is seen at 70–90 s. The time variation in the precollapse neutrino detection thus can be used as the indicator of the multidimensional convective activity deep inside the star, although it is practically impossible to detect such a time variation with current one-kiloton-size neutrino detectors.

Zoom In Zoom Out Reset image size

Next, we examine the detectability of precollapse neutrinos by Hyper-Kamiokande. Hyper-Kamiokande is a planned water-Čerenkov-type neutrino detector, which has a huge fiducial mass of 190 kton (Abe et al. 2018). If a supernova explodes in the vicinity of Earth, high-quality data of supernova neutrinos with well-resolved time and energy bins are expected to be obtained (e.g., Takahashi et al. 2001; Mirizzi et al. 2016). However, because the threshold energy of Hyper-Kamiokande is expected to be at the same level as Super-Kamiokande (∼4.79 MeV corresponding to 3.5 MeV for recoil electrons; Sekiya 2013) and is higher than that of a liquid-scintillation-type detector, Hyper-Kamiokande is less preferable for observations of neutrinos produced through the pair-neutrino process during the precollapse stage. Yoshida et al. (2016) discussed the possibility of neutrino observations using delayed γ-rays from Gd in Gd-loaded Hyper-Kamiokande (Beacom & Vagins 2004), because the energy threshold will be reduced to 1.8 MeV, an energy similar to the case of KamLAND, in this case. Considering a moderate detection efficiency of 50%, the detection rate is about 178 times as large as the rate of KamLAND. When the enhancement factor of 178 for the detection rate is applied, the event rate is expected to be 3–14 s⁻¹ in the normal mass ordering. So, the time variation due to the convective motion with a timescale of seconds could be observed by Hyper-Kamiokande.

We should note that we consider only the pair-neutrino process for presupernova neutrinos in this study. For a few minutes before core-collapse, the main source of ${\bar{\nu }}_{e}$ will be β⁻ decays rather than the pair-neutrino process (Kato et al. 2017). β⁻ decays mainly take place in the innermost Fe core, and the ${\bar{\nu }}_{e}$ energy is similar to pair neutrinos. We expect that the neutrino event rate by β⁻ decays also decreases for 70–90 s because the central temperature and density decrease during this period. So, even when we take into account the neutrino event rate by β⁻ decays, the time variation of presupernova neutrino events would give us information about the dynamics in the SiO-rich layer or a collapsing Fe core.