SYSTEMATIC STUDIES OF SHOCK REVIVAL AND THE SUBSEQUENT EVOLUTIONS IN CORE-COLLAPSE SUPERNOVAE WITH PARAMETRIC PROGENITOR MODELS

The methodology employed in this paper is essentially the same as that in our previous paper (Yamamoto et al. 2013). One single big difference is, as mentioned in Introduction, that multiple progenitor models are constructed systematically, as described in the next subsection in detail, and study the dependence of the energy and Ni mass in ejecta on the progenitor structure. In order to evaluate the asymptotic values of these quantities, we take the following four steps.

• 1.  
  Construct the toy models of pre-supernova stars with different masses of iron cores and Si+S layers.
• 2.  
  Run 1D hydrodynamical simulations of gravitational collapse of the pre-supernova models obtained in the first step, excising the central part and thus reducing the pressure support there by hand. This step is needed to obtain the mass accretion histories for these models, which are used in the following steps.
• 3.  
  Construct spherically symmetric, steady accretion flows through stalled shock waves for mass accretion rates obtained in the first step and neutrino luminosities, which are chosen to be close to the critical values for the assumed mass accretion rates.
• 4.  
  Run 1D and 2D hydrodynamical simulations for the accretion flows just constructed. Small perturbations are added in the 2D computations. The matter distributions obtained in the first step are employed outside the shock wave. The neutrino luminosities change in time in a prescribed way. Since the initial configurations are close enough to the critical states, shock revival takes place soon after the simulations are started. The ensuing evolutions are followed until the explosion energy is well settled.

It is noted that the determination of the critical curve (see Burrows & Goshy 1993; Yamasaki & Yamada 2005; Pejcha & Thompson 2012; Suwa et al. 2014, for details), which corresponds to the third step, is slightly changed from the previous work. We follow Murphy & Burrows (2008) to take account of the temporal change in the mass accreation rate which is commonly used by others these days.

It should be also mentioned that the simulations in the final step are terminated when the energy in the ejecta becomes almost unchanged. In 1D models this also means the settlement of the Ni mass in the ejecta by that time. This is not the case, however, in 2D computations unfortunately, since larger fractions of ejecta tend to fallback at later times and it takes much longer to get the Ni mass settled: in some cases the explosion energy reaches the final value in ∼1 s, whereas the Ni mass has not reached the asymptotic value yet even after ∼10–100 s. The simulations are terminated before the Ni mass is finally settled in those cases, since a further computation would require us to take into account the outer envelopes such as helium and hydrogen layers, which would not be a trivial thing. It is also stressed that we are concerned with the influences of the structures of the stellar central region involving the iron core and Si+S layers on the shck revival and the dynamics that follows. We hence decide to avoid unnecessary complications in this paper. The earlier termination of simulations would underestimate the fallback and hence overestimate the Ni mass. Then we should regard the numbers given in that case as upper limits, which would not spoil our arguments anyway.

As we have mentioned already, we do not use realistic progenitor models that are put in the public domain by researchers of stellar evolutions (Woosley et al. 2002; Limongi & Chieffi 2006; Paxton et al. 2011; Yoshida & Umeda 2011; Ekström et al. 2012), but instead construct toy models on our own. This is mainly because those realistic models have structures that are sensitive to the stellar mass as well as some physical processes, e.g., convections of various kinds, and their numerical implementations (Langer 2012; Sukhbold & Woosley 2014), which would make it difficult to decipher what is the consequence of which effect. Since our studies are of experimental nature, "realism" is not very important and toy models that capture the essential features in massive stars' structures but still have some degrees of freedom to change key parameters, e.g., the masses of Fe core and Si+S layer, rather arbitrarily, which is more useful. In this subsection, we present in detail the construction of such toy models.

It is well known that the masses of the Fe core and Si+S layer are critically important for shock revival (Buras et al. 2006; Suwa et al. 2014) because both of them one way or another affect the neutrino luminosity and mass accretion rate, the main controlling parameters of the stagnant shock wave. They are notoriously stochastic as a function of the stellar mass, though. Figure 1 illustrates this, showing the masses of the Fe core and Si+S layer for various pre-supernova models with 11–40 M_⊙ taken from Woosley et al. (2002, hereafter WHW2002). It is evident that neither the iron core mass, M_c, nor the Si+S mass, M_SiS, is a monotonic function of the progenitor mass.

Zoom In Zoom Out Reset image size

If we pay attention, however, to the entropy, S, and electron fraction, Y_e as a function of density inside the Fe core, which is defined hereafter to be the central region in the NSE, those apparently diverse progenitors can be nicely categorized into three groups, as demonstrated for the same pre-supernova models in Figure 2. The first group is characterized by the relatively low entropies, i.e., S < 1.75k_B at the edge of the Fe core and is plotted in the left panel with the label "L." Here k_B is the Boltzmann's constant. On the other hand, the right panel shows the models that have rather high entropies, 2.3k_B < S, at the edge and are referred to with the label "H." In between fall most of the pre-supernova models, as demonstrated in the middle panel labeled with "M." They are conveniently characterized by the moderate entropies of 1.75k_B ≤ S ≤ 2.3k_B at the edge of Fe core. The important thing is that the members of each group seem to obey approximately the same ρ–S and ρ–Y_e relations irrespective of their Fe core masses. One may complain that the boundary between groups "M" and "H" is rather blurred as long as the value of the entropy at the edge of Fe core is concerned. Note, however, the different locations of the jump in the entropy distributions between the two groups: it occurs at a lower density for group "H" than for group "M." It is also mentioned that the members of group "H" have larger values of compactness ξ_2.5 (O'Connor & Ott 2011; Sukhbold & Woosley 2014) and may form a black hole instead of a neutron star.

Zoom In Zoom Out Reset image size

It is found from Figure 2 that S and Y_e as functions of ρ have some features in common among these groups: plotted as in the figure, they are almost linearly increasing from the central point of star corresponding to the right ends of lines until a certain point, where S and Y_e jump nearly discontinuously; the positions of the latter points are different among the three groups, as mentioned already, and mark the locations at which the shell Si burning once took place. It induced violent convection outside the shell and produced almost homogeneous entropy and Y_e distributions there, which were later bent downwards at smaller radii by electron captures on nuclei during the subsequent phase. The resultant distributions of S and Y_e are approximated by parabolic functions quite well. Inside the Si burning shell, on the other hand, matter stratification is stable against convection mainly due to efficient neutrino cooling in the central region. As a result, the linear distributions just mentioned are obtained there for S and Y_e. These functional relations are later used to build the toy model.

We also point out the correlation between the central density and the mass of the Fe core as shown in Figure 3, although there are some scatters: massive cores tend to have smaller central densities just prior to collapse. This is mainly because of the generic trend that the more massive the core is, the higher the entropy becomes due to the shorter lifetime and, as a consequence, the less efficient neutrino cooling. We will use this correlation in the following when we construct the toy model.

Zoom In Zoom Out Reset image size

It is stressed again that the purpose of the toy model is to reproduce key features commonly observed in the massive star structures just prior to collapse, simultaneously leaving some degree of freedom to modify them rather arbitrarily for systematic studies. We thus might be able to suggest a feature or features that may be essential to robust shock revival but may have somehow eluded the current stellar evolution theory. In this series of papers, we will study the effects on the shock revival of various features in the massive star structures reproduced by the toy model one after another. In this first paper, we pay particular attention to the masses of the iron core, M_c, and of the Si+S layer, M_SiS, which certainly play major roles in the neutrino-heating mechanism, dictating the neutrino luminosity and accretion rate, the two key parameters in the mechanism. Other features may be equally or more important and will be investigated in the sequels (Y. Yamamoto & S. Yamada 2016, in preparation).

Among three groups we focus on group M in this paper, since it is the largest group, to which most of the progenitor models we show in Figure 2 belong, and it is hence supposed to represent the canonical progenitor of CCSNe. We then use the following parametrization for the entropy S and electron fraction Y_e as functions of ρ:

$$\begin{eqnarray}&&S(\rho )\;=\\ \left\{\begin{array}{lc}\displaystyle \frac{({S}_{1}-{S}_{0})}{(\mathrm{log}{\rho }_{1}-\mathrm{log}{\rho }_{0})}\cdot (\mathrm{log}\rho -\mathrm{log}{\rho }_{0})+{S}_{0} & (\rho \geqslant {\rho }_{1}),\\ \displaystyle \frac{({S}_{2}-{S}_{3})}{{(\mathrm{log}{\rho }_{2}-\mathrm{log}{\rho }_{3})}^{2}}\cdot {(\mathrm{log}\rho -\mathrm{log}{\rho }_{3})}^{2}+{S}_{3} & (\rho \leqslant {\rho }_{2}),\end{array}\right.\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{Y}_{e}(\rho )\;=\\ \left\{\begin{array}{lc}\displaystyle \frac{({Y}_{e1}-{Y}_{e0})}{(\mathrm{log}{\rho }_{1}-\mathrm{log}{\rho }_{0})}\cdot (\mathrm{log}\rho -\mathrm{log}{\rho }_{0})+{Y}_{e0} & (\rho \geqslant {\rho }_{1}),\\ \displaystyle \frac{({Y}_{e2}-{Y}_{e3})}{{(\mathrm{log}{\rho }_{2}-\mathrm{log}{\rho }_{3})}^{2}}\cdot {(\mathrm{log}\rho -\mathrm{log}{\rho }_{3})}^{2}+{Y}_{e3} & (\rho \leqslant {\rho }_{2}),\end{array}\right.\end{eqnarray} \tag{ 2 }$$

in which ρ₁ = 10^8.5 g cm⁻³ and ρ₂ = 10^8.4 g cm⁻³ specify the region where the jump occurs, and the other parameters determine the shapes of the functions and their values are set to ρ₀ = 10^9.8 g cm⁻³, ρ₃ = 10^7.1 g cm⁻³, S₀ = 0.70k_B, S₁ = 1.50k_B, S₂ = 2.00k_B, and S₃ = 2.30k_B and Y_e0 = 0.435, Y_e1 = 0.467, Y_e2 = 0.465, and Y_e3 = 0.482. In the transition region between ρ₁ and ρ₂, S and Y_e are interpolated linearly in $\mathrm{log}\rho$. Figure 4 displays the resultant relations (black dotted lines) between S and ρ in the left panel and between Y_e and ρ in the right panel together with the 30 pre-supernova models from WHW2002 (colored solid lines), which consist of group "M." It is clear that the essential features of the group are well reproduced.

Zoom In Zoom Out Reset image size

Next we shift our attention to outer layers. We begin with the Si+S layer, which neighbors the Fe core. The important feature of this layer is that nuclei are no longer in the NSE but the temperature is still high enough to maintain the so-called QSE (Hix & Thielemann 1999), in which the Fe group, Si group, and another group of lighter elements separately establish chemical equilibrium among the elements belonging to each group (see Appendix A).

The treatment of other outer envelopes, in which no chemical equilibrium is established even approximately, are simplified further: homogeneous distributions of entropy, electron fraction, and chemical abundance are assumed in each of these layers. This is not so bad an approximation, though, since convection once prevailed in most of these layers, making the entropy and chemical abundance uniform there. We take the typical values from WHW2002 (see Table 1). As explained shortly, since we solve the hydrostatic equations to build the toy models, we assume the continuity of pressure at their boundaries. The parameter values we choose are further constrained if one were to calculate the evolutions of chemical abundances later, since the timescale of nuclear reactions, ${\tau }_{{\rm{nuc}}}={\mathrm{min}}_{k}{Y}_{k}/\dot{{Y}_{k}}$, should be longer than the dynamical timescale, ${\tau }_{{\rm{ff}}}=\sqrt{{r}^{3}/{{GM}}_{r}}$, in every layer. Note that τ_nuc may be shorter than the typical time for the consumption of nuclear fuel owing to the rises of density by core contraction, especially in silicon and oxygen–neon–magnesium layers.

Table 1.  Parameter Sets used in Toy Pre-supernova Stage Construction

	Core	QSE	Si+S+O	O+Mg+Si	O+Ne+Mg	C+O	C+O+He	He
M (M⊙)	Mc	MSiS	0.36-MSiS	0.09	2.21-Mc	0.15	≳0.10	≲1.10
S$({k}_{{\rm{B}}}/\mathrm{nuc})$	Equation (1)	3.3	4.0	5.0	5.9	6.5	7.6	12.0
Chemical	⋯	⋯	XSi = 0.45,	XO = 0.80,	XO = 0.70,	XC = 0.20,	XC = 0.30,	XHe = 1.00
Abundance	NSE	QSE	XS = 0.35,	XMg = 0.15,	XNe = 0.25,	XO = 0.80	XO = 0.60,	⋯
Xk	⋯	⋯	XO = 0.20	XSi = 0.05	XMg = 0.05	⋯	XHe = 0.10	⋯

Note. In this paper, we pay attention mainly to the core mass, M_c = 1.3, 1.4, 1.5 M_⊙ and Si+S layer mass, M_SiS = 0.09, 0.18 M_⊙.

Download table as:  ASCIITypeset image

As already mentioned repeatedly, the toy models are obtained by solving the hydrostatic equations for non-rotating, spherically symmetric stars. Instead of solving the equations for energy generation and transport, we employ the functional relations described above for the Fe core and the constant S and Y_e distributions in the outer layers. As for the chemical abundances, NSE and QSE are assumed in the Fe core and Si+S layer, respectively. The uniform distributions of various elements are assumed in other outer layers. Then the basic equations for the NSE core are given as

$$\begin{eqnarray}&&\displaystyle \frac{{dr}}{{{dM}}_{r}}=\displaystyle \frac{1}{4\pi {r}^{2}\rho },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\rho }{{{dM}}_{r}}=\displaystyle \frac{-(1-{f}_{P})\cdot \displaystyle \frac{{{GM}}_{r}}{4\pi {r}^{4}}}{\left({\left(\displaystyle \frac{\partial P}{\partial \rho }\right)}_{S,{Y}_{e}}+{\left(\displaystyle \frac{\partial P}{\partial S}\right)}_{\rho ,{Y}_{e}}\displaystyle \frac{{dS}}{d\rho }+{\left(\displaystyle \frac{\partial P}{\partial {Y}_{e}}\right)}_{\rho ,S}\displaystyle \frac{{{dY}}_{e}}{d\rho }\right)},\end{eqnarray} \tag{ 4 }$$

where r, M_r, and P are the radius, mass coordinate, and pressure, respectively. f_P is a parameter to allow for some deviation from the complete hydrostatic equilibrium, which is commonly the case for realistic progenitor models provided by stellar evolution calculations, in which gravitational contraction has already sets in. It is assumed to be given by the following form:

$$\begin{eqnarray}&&{f}_{P}=0.05\mathrm{exp}(-0.5{(\mathrm{log}(\rho )-7.2)}^{2}),\end{eqnarray} \tag{ 5 }$$

which is found to fit the 30 realistic progenitor models with 11 ≤ M_ZAMS ≤ 40 that have entered the collapse phase. As for the outer envelopes, the basic equations are essentially the same except the piecewise constant S, Y_e, and element distributions, as well as f_p = 0. Since the layer boundaries are contact surfaces, pressure is assumed to be continuous there.

In this paper, we focus on the masses of the Fe core and Si+S layer. Six different combinations of three Fe core masses, M_c = 1.3 M_⊙, 1.4 M_⊙ and 1.5 M_⊙, and two Si+S layer masses, M_SiS = 0.09 M_⊙ and 0.18 M_⊙, are investigated. Given the Fe core mass, the central density is given by the following equation:

$$\begin{eqnarray}&&{\rho }_{c}({M}_{c})=6.5\times {10}^{9}\mathrm{exp}[-5\mathrm{ln}2({M}_{c}-1.40)],\end{eqnarray} \tag{ 6 }$$

which is a fit to the correlation suggested in Figure 3 (see the black dotted line in the figure). Then the hydrostatic equations are integrated from the center to the specified masses of the Fe core and Si+S layer to give the corresponding radii. The integration is actually continued to larger radii, with other outer layers being properly taken into account.

The six solutions are shown in Figure 5, which depicts the enclosed masses and entropies per baryon as functions of radius for the three Fe core masses with the Si+S layer mass of M_SiS = 0.09 M_⊙ in the left panel and with M_SiS = 0.18 M_⊙ in the right panel, respectively. Note that we adopt different masses of C+O+He and He layers for the models with M_c = 1.30 M_⊙ and with M_c = 1.40 M_⊙ in the case of M_SiS = 0.09 M_⊙ in order to avoid unrealistic He burnings, although this has essentially no influence on the explosion energy and Ni mass.

Zoom In Zoom Out Reset image size

It is evident that the rather small differences in the core mass, ΔM_c = 0.10 M_⊙, and the Si+S layer mass, ΔM_SiS = 0.09 M_⊙, cause noticeable changes in the R–M_r relations and hence the density profiles in outer layers (M_r ≳ 2 M_⊙): higher M_SiS tends to make the greater differences at large radii for the same ΔM_c. The compactness introduced by O'Connor & Ott (2011) also differs considerably among these models. The purpose of this paper is to study the consequences these differences may have on the shock revival, explosion energy, and Ni mass later in the dynamics of CCSNe.

It is noteworthy that the hierarchy in compactness depends on where it is measured: heavier Fe core models have higher compactnesses if they are measured in the silicon layer, whereas the hierarchy is inverted outside the ONeMg layer. This is because the entropy is higher at M_r ∼ 1.5–2 M_⊙ for lighter Fe cores, making this layer more bloated. In the O+Ne+Mg envelope at M_r ∼ 2–2.7 M_⊙, on the other hand, all models are assumed to have the same entropy again. Since the O+Ne+Mg envelope appears at lower densities in the models with heavier Fe cores, as can be seen in the lower panels of the figure, temperatures are higher in this layer for these models, which results in the more extended envelopes.

This feature was also reported in Baron & Cooperstein (1990; see their models "105" and "107" in Figure 9), who attempted to construct their own toy models, having similar parametric research in mind. They gave S and Y_e as functions of mass coordinate unlike in this paper, in which they are considered as functions of ρ. It seems that they paid little attention to the correlation between S and Y_e, both of which are influenced by the electron capture on heavy nuclei in a similar way. Note also that the two masses, M_c and M_SiS, cannot be chosen completely arbitrarily with other parameters being fixed, since we perform network calculations to follow the temporal evolution of chemical abundances, the initial chemical compositions should be stable, i.e., the timescale of nuclear reactions should be longer than the dynamical timescale. It is also found that very high M_c and/or M_SiS produce too large low-entropy regions and make it difficult to extend outer envelopes. Since we perform network calculations to follow the temporal evolution of chemical abundances, their models are not appropriate for our purpose.

The numerical method employed in this paper for hydrodynamical computations is essentially the same as that in our previous work (Yamamoto et al. 2013). The hydrodynamical evolutions are solved on the spherical coordinates with the axisymmetric Eulerian code ZEUS-2D (Stone & Norman 1992; Hayes et al. 2006), which is extended by the authors so that chemical abundance evolutions could be treated consistently with a non-NSE EOS. For the 1D collapse computation (step 2 in Section 2.1), 500 radial grid points are deployed to resolve 100 km ≤ r ≤ 20,000 km. For the simulations of shock revival and the subsequent evolution (steps 3 and 4), we employ 650 radial grid points in 1D and $({n}_{r},{n}_{\theta })=(600,76)$ mesh points in 2D, which cover the domain with 45 km ≤ r ≤ 20,000 km and 0 ≤ θ ≤ 2 π. The central portion of the core is excised and replaced by a point mass located at the center, which exerts gravitational forces to the matter in the computational region. The temporal change in its mass is calculated from the matter flow going into the inner boundary. The inner boundary condition is essentially the same as in the previous paper (Yamamoto et al. 2013): all variables but radial velocity are copied from the innermost active mesh point to the ghost mesh points. The radial velocity is fixed until outflow occurs, in which case it is treated just in the same way as other variables except for the upper limit of 10⁸ cm s⁻¹. The densities on the ghost mesh points are also modified so that the entropy per baryon should not exceed 100k_B. We use the same code both for 1D and 2D simulations; in the former of which we suppress motions in the θ direction.

We adopt the light bulb approximation instead of solving neutrino transfer and take into consideration only interactions of electron-type neutrinos and anti-neutrinos with nucleons when we conduct simulations in steps 3 and 4 for the shock revival and subsequent evolution. To approximately describe the neutrino-spheric radius, R_ν, and temperature, T_ν, that evolve in time due to PNS contractions, we employ the following analytical prescriptions for R_ν:

$$\begin{eqnarray}&&{R}_{{\nu }_{e}}({t}_{300})=31.00+16.69\cdot \mathrm{exp}\left(-\displaystyle \frac{{t}_{300}}{193}\right)\;{\rm{km}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{R}_{{\bar{\nu }}_{e}}({t}_{300})=27.00+16.45\cdot \mathrm{exp}\left(-\displaystyle \frac{{t}_{300}}{211}\right)\;{\rm{km}},\end{eqnarray} \tag{ 8 }$$

and for T_ν:

$$\begin{eqnarray}&&{T}_{{\nu }_{e}}({t}_{300})=6.00-0.67\cdot \mathrm{exp}\left(-\displaystyle \frac{{t}_{300}}{200}\right)\;{\rm{MeV}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{T}_{{\bar{\nu }}_{e}}({t}_{300})=7.00-0.67\cdot \mathrm{exp}\left(-\displaystyle \frac{{t}_{300}}{200}\right)\;{\rm{MeV}},\end{eqnarray} \tag{ 10 }$$

where t₃₀₀ denotes the time in millisecond elapsed from 300 ms post-bounce. The prescriptions are determined so that they should fit more realistic simulations (Nakamura et al. 2015). The luminosities of electron-type neutrinos and anti-neutrinos are assumed to be same, ${L}_{\nu }={L}_{{\nu }_{e}}={L}_{{\bar{\nu }}_{e}}$, gradually decay after the onset of explosion:

$$\begin{eqnarray}{L}_{\nu }(t)=\left\{\begin{array}{ll}{L}_{\nu } & \mathrm{before}\ \mathrm{explosion},\\ {L}_{\nu }\cdot \mathrm{exp}\left(-\displaystyle \frac{{t}_{300}}{800\;{\rm{ms}}}\right) & \mathrm{after}\ \mathrm{explosion}.\end{array}\right.\end{eqnarray} \tag{ 11 }$$

Another major improvement from the previous work is the nuclear physics employed. In order to calculate chemical compositions more precisely, we deploy 297 nuclei ($Z\quad \leqq \quad 32$) instead of 28 to describe NSE at T > 7.0 × 10⁹ K (see the list in Table 2). At T < 7 × 10⁹ K, where NSE no longer holds, on the other hand, the number is reduced back to 28: n, p, D, T, ³He, ⁴He, and 12 α-nuclei, i.e., ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti, ⁴⁸Cr, ⁵²Fe, ⁵⁶Ni, and 10 of their neutron-rich neighbors, that is, ${}^{27}\mathrm{Al}$, ${}^{31}{\rm{P}}$, ${}^{35}\mathrm{Cl}$, ${}^{39}{\rm{K}}$, ${}^{43}\mathrm{Sc}$, ${}^{47}{\rm{V}}$, ${}^{51}\mathrm{Mn}$, ${}^{53}\mathrm{Fe}$, ${}^{54}\mathrm{Fe}$, and ${}^{55}\mathrm{Co}$. The other 269 nuclei are represented by a single (imaginary) inert nucleus in the network calculations.

We present first the results of 1D computations of the collapse phase, i.e., Step 2 mentioned in Section 2.1. The time evolutions of the mass accretion rates, $\dot{M}$, measured at r = 100 km are depicted in Figure 6. After the onset of collapse, $\dot{M}$ first rises rapidly, reaching ∼100 M_⊙ s⁻¹ at maximum, and then drops rather quickly to ∼1 M_⊙ s⁻¹. The entire core shrinks to r ≤ 100 km at around the peak time, which roughly corresponds to the core-bounce time in realistic simulations and hence which we refer to as the bounce time in the following. The mass accretion rates decay more gradually thereafter with a few discernible steps, which correspond to the passages of layer boundaries. These behaviors are qualitatively consistent with Buras et al. (2006), Yamamoto et al. (2013), and Suwa et al. (2014). More detailed analyses are given in Appendix B.

Zoom In Zoom Out Reset image size

It is interesting to point out that if we horizontally shift the curves of $\dot{M}$ for the models with different core masses but with the same Si+S mass so that the bounce times should coincide with each other, then the entire curves also agree rather well with each other as demonstrated in the right panels of Figure 6, except that the steps occur at different post-bounce times. Comparing models with different M_SiS but with the same M_c, we find that the outer boundary of Si+S+O layer reaches r = 100 km at almost the same post-bounce time irrespective of M_SiS, although the outer edges of Si+S layer arrive at different times. This is mainly due to the fact that the enclosed mass at the outer boundary of the Si+S+O layer is common. It is interesting that the difference in M_SiS is reflected in the accretion of the Si+S layer alone, not affecting the infall of outer layers.

Note in particular the difference in the positions of the steps in the mass accretion rate as a function of time, which are produced by the passage of boundaries between different layers. This is a direct consequence of the different masses of Fe core and Si+S layer. As discussed later, it is reflected in the shape of the critical curve (see Section 3.2.1 and also Appendix C).

To explore the shock revival by neutrino heating and the subsequent evolution with hydrodynamical simulations in 1D and 2D, we first construct spherically symmetric, steady, shocked accretion flows for initial conditions (Step 3 in Section 2). They are characterized by three parameters, i.e., the mass accretion rate, neutrino luminosity, and mass of the PNS. As for the first and third parameters, we adopt the post-bounce time of the accretion histories in Step 2 at t_pb = 300 ms. We have three options for the second parameter, i.e., the neutrino luminosity, on the other hand: ${L}_{\nu }={L}_{{\nu }_{e}}={L}_{{\bar{\nu }}_{e}}$, L_ν = 2.0, 2,3, 2.5 × 10⁵² erg s⁻¹, or L_ν,52 = 2.0, 2.3, 2.5, which are applied to each of the six progenitor models constructed in Step 1. The properties of these 18 shocked accretion flows are summarized in Table 3.

Table 3.  The List of the Steady Solution Properties at t_pb = 300 ms in Different Neutrino Luminosity

Modela	$\dot{M}$ (M⊙ s−1)	Mgb (M⊙)	Rgc (cm)	Rshd (cm)	MPNSe (M⊙)
1.3S09Lnu2.0	5.20E−01	6.67E−03	8.61E+06	1.38E+07	1.49E+00
1.4S09Lnu2.0	4.94E−01	5.86E−03	8.56E+06	1.36E+07	1.56E+00
1.5S09Lnu2.0	5.13E−01	4.71E−03	8.43E+06	1.28E+07	1.65E+00
1.3S18Lnu2.0	4.47E−01	7.58E−03	8.79E+06	1.50E+07	1.50E+00
1.4S18Lnu2.0	4.50E−01	6.21E−03	8.64E+06	1.42E+07	1.58E+00
1.5S18Lnu2.0	4.71E−01	4.92E−03	8.48E+06	1.33E+07	1.66E+00
1.3S09Lnu2.3	5.20E−01	8.71E−03	8.84E+06	1.51E+07	1.48E+00
1.4S09Lnu2.3	4.94E−01	7.58E−03	8.75E+06	1.48E+07	1.56E+00
1.5S09Lnu2.3	5.13E−01	5.99E−03	8.55E+06	1.38E+07	1.64E+00
1.3S18Lnu2.3	4.47E−01	1.02E−02	9.09E+06	1.68E+07	1.49E+00
1.4S18Lnu2.3	4.50E−01	8.15E−03	8.86E+06	1.56E+07	1.58E+00
1.5S18Lnu2.3	4.71E−01	6.30E−03	8.61E+06	1.44E+07	1.66E+00
1.3S09Lnu2.5	5.20E−01	1.05E−02	9.04E+06	1.62E+07	1.48E+00
1.4S09Lnu2.5	4.94E−01	9.11E−03	8.91E+06	1.59E+07	1.56E+00
1.5S09Lnu2.5	5.13E−01	7.04E−03	8.66E+06	1.46E+07	1.64E+00
1.3S18Lnu2.5	4.47E−01	1.28E−02	9.34E+06	1.84E+07	1.49E+00
1.4S18Lnu2.5	4.50E−01	9.88E−03	9.05E+06	1.68E+07	1.57E+00
1.5S18Lnu2.5	4.71E−01	7.45E−03	8.75E+06	1.52E+07	1.66E+00

Notes.

^aThe model names is consisted of core mass, type of Si+S mass, and luminosity in order. ^bMass inside the gain region. ^cRadius at gain region. ^dRadius at shock front. ^eThe enclosed mass inside the inner boundary.

Download table as:  ASCIITypeset image

Using all of these 18 models as initial conditions, we run simulations assuming spherical symmetry. We also conduct 2D computations for eight models, which generate robust explosions in 1D. In these simulations, we employ in the region outside the stalled shock wave the profiles of the accreting envelopes at t_pb = 300 ms, which were obtained in the calculations of the accretion histories in Step 2. As a result of this, the mass accretion rate varies in time. We keep the neutrino luminosity constant in time until shock revival occurs. Once it occurs, we allow the neutrino luminosity to evolve in time according to Equation (11). The results are summarized in Table 4, which will be used frequently later for discussion. In this table, the mass accretion rates at the shock revival, ${\dot{M}}_{{\rm{rev}}}$, are sampled at r = 400 km and are essentially the same as those on the shock front. The diagnostic explosion energy, which is denoted by E_exp, is defined as

$$\begin{eqnarray}&&{E}_{\mathrm{exp}}=\displaystyle \sum _{{v}_{r}\gt 0,{e}_{\mathrm{tot}}\gt 0}^{}{e}_{\mathrm{tot}}\ {\rm{\Delta }}V,\end{eqnarray} \tag{ 12 }$$

in which the total energy density, e_tot, is given by the sum of kinetic (e_kin), internal (e_int), and gravitational (e_grav) energy densities as

$$\begin{eqnarray}&&{e}_{\mathrm{tot}}={e}_{\mathrm{kin}}+{e}_{\mathrm{int}}+{e}_{\mathrm{grav}};\end{eqnarray} \tag{ 13 }$$

the summation in Equation (12) extends over the grid points with positive radial velocities, v_r, and total energy densities; ΔV stands for the cell volumes for these grid points. Note that the rest mass energy is not included in the internal energy in the above definition. The diagnostic energy is also used in judging shock revival: we consider that the stalled shock wave revives if the diagnostic energy reaches 10⁴⁸ erg.

Table 4.  Summary of the Explosion Results in 1D and 2D

Model	Eexp (foe)	MNi (M⊙)	treva (s)	MPNS (M⊙)	${\dot{M}}_{{\rm{rev}}}$ b (M⊙ s−1)
1D
1.3S09Lnu2.0	8.40E−01	8.80E−02	7.40E−01	1.68E+00	3.09E−01
1.4S09Lnu2.0	6.70E−01	7.00E−02	8.30E−01	1.78E+00	2.54E−01
1.5S09Lnu2.0	4.40E−01	4.95E−02	1.00E+00	1.91E+00	2.16E−01
1.3S18Lnu2.0	7.50E−01	7.50E−02	6.00E−01	1.62E+00	2.77E−01
1.4S18Lnu2.0	5.30E−01	5.70E−02	7.00E−01	1.73E+00	2.22E−01
1.5S18Lnu2.0	4.20E−01	4.67E−02	7.90E−01	1.85E+00	1.92E−01
1.3S09Lnu2.3	1.22E+00	1.08E−01	5.70E−01	1.60E+00	3.63E−01
1.4S09Lnu2.3	8.20E−01	9.05E−02	6.40E−01	1.72E+00	3.18E−01
1.5S09Lnu2.3	6.10E−01	7.20E−02	7.80E−01	1.85E+00	2.48E−01
1.3S18Lnu2.3	1.37E+00	1.09E−01	4.30E−01	1.54E+00	4.06E−01
1.4S18Lnu2.3	1.13E+00	9.10E−02	4.70E−01	1.64E+00	3.60E−01
1.5S18Lnu2.3	8.14E−01	7.20E−02	7.00E−01	1.78E+00	3.01E−01
1.3S09Lnu2.5	1.64E+00	1.49E−01	4.40E−01	1.51E+00	4.64E−01
1.4S09Lnu2.5	1.43E+00	1.22E−01	4.80E−01	1.63E+00	4.27E−01
1.5S09Lnu2.5	7.30E−01	8.40E−02	6.70E−01	1.82E+00	3.10E−01
1.3S18Lnu2.5	1.40E+00	1.24E−01	3.70E−01	1.52E+00	4.32E−01
1.4S18Lnu2.5	1.45E+00	9.60E−02	4.30E−01	1.61E+00	3.83E−01
1.5S18Lnu2.5	1.07E+00	9.40E−02	4.70E−01	1.73E+00	3.66E−01
2D
1.3S09Lnu2.0	8.96E−01	3.95E−02	6.16E−01	1.63E+00	3.30E−01
1.4S09Lnu2.0	7.60E−01	6.08E−02	6.40E−01	1.72E+00	3.10E−01
1.3S18Lnu2.0	9.70E−01	3.45E−02	5.34E−01	1.60E+00	3.51E−01
1.4S18Lnu2.0	7.46E−01	4.03E−02	5.25E−01	1.72E+00	3.10E−01
1.3S09Lnu2.3	1.55E+00	1.01E−01	4.07E−01	1.53E+00	4.72E−01
1.4S09Lnu2.3	1.46E+00	1.08E−01	3.95E−01	1.61E+00	4.69E−01
1.3S18Lnu2.3	1.48E+00	1.04E−02	3.92E−01	1.50E+00	4.36E−01
1.4S18Lnu2.3	1.26E+00	6.27E−02	3.96E−01	1.63E+00	4.03E−01

Notes.

^aThe shock revival time measured from post-bounce. ^bThe mass accretion rate at t_rev.

Download table as:  ASCIITypeset image

3.2.1. Critical Luminosity

The luminosities, mass accretion rates, as well as PNS mass at shock revival are presented in Figure 7 for the 18 spherical models. The upper panels display the results for the models with M_SiS = 0.09 M_⊙, whereas the lower panels are the counterparts for the models with M_SiS = 0.18 M_⊙. The left panels exhibit the so-called critical curves, in which the neutrino luminosity at shock revival is regarded as a function of the mass accretion rate at the same time, ${\dot{M}}_{{\rm{rev}}}$. It is evident that critical luminosity is a monotonically increasing function of the mass accretion rate, which is consistent with other studies (Murphy & Burrows 2008; Nordhaus et al. 2010; Hanke et al. 2012; Couch 2013; Iwakami et al. 2014). It is interesting to point out that the curve is upward convex for the models with M_SiS = 0.09 M_⊙ while it is downward convex for the models with M_SiS = 0.18 M_⊙. This is related with the locations of layer boundaries, at which the mass accretion rate decreases rapidly, and will be discussed in more detail in Appendix C. It is also clear that the larger the core mass is, the higher the critical luminosity is, since the gravitational binding becomes stronger.

Zoom In Zoom Out Reset image size

In the right panels of Figure 7, we give the trajectories of the mass accretion rate and the PNS mass for different masses of the Fe core and Si+S layer. As time passes, the mass accretion rate gets lowered while the PNS mass becomes larger. Note that the mass accretion rates measured at r = 400 km are independent of the neutrino luminosities. The points, at which shock revival occurs for the 18 models, are marked with squares in the figure. As the luminosity decreases, shock revival is delayed and, as a result, the mass accretion rate becomes lower. Note that the positions of the squares are slightly off the trajectories. In fact, we employ the PNS masses at the end of simulations instead of those at shock revival in plotting these squares. Hence, their vertical deviations from the trajectories indicate the masses gained or lost by PNS owing to the inflow or outflow of matter through the inner boundary after shock revival. In the 1D case, the mass is always lost, since the outflow is dominant via neutrino winds. This is not always the case in 2D as shown later, in which the neutrino-driven wind occurs more strongly.

In Figure 8 we plot the time evolutions of the mass accretion rates and the shock revival points for the three core masses: M_c = 1.30 M_⊙ (left panel), 1.40 M_⊙ (middle panel), and 1.50 M_⊙ (right panel). It is evident from these plots that shock revival is delayed for more massive iron cores. This is simply because of the stronger gravitational attraction. The higher mass accretion rates for the model with M_SiS = 0.09 M_⊙ leads to the fact that it obtains shock revival later than the model with M_SiS = 0.18 M_⊙ irrespective of the core mass (see Appendix C). It may be interesting to see for the highest luminosity (L_ν,52 = 2.5) case in the model with M_SiS = 0.09 M_⊙ that shock revival occurs before the passage of the boundary between the Si+S+O layer and the O+Mg+Si layer for the lighter cores with M_c = 1.30 M_⊙ and 1.40 M_⊙, whereas it happens during the passage for the heaviest core with M_c = 1.50 M_⊙. As a result, the mass accretion rate at shock revival becomes smaller in this case than the counterpart in the model with M_SiS = 0.18 M_⊙. In fact, the layer boundary is located at a larger radius in the latter progenitor model and shock revival occurs before its passage even for the heaviest core mass, which is also the reason why the critical curve is almost straight in this case.

Zoom In Zoom Out Reset image size

3.2.2. Explosion Energy and Nickel Mass

We turn attention to the explosion energies and Ni yields in the 1D models, which are summarized in the second and third columns in Table 4. For fixed masses of Fe core and the Si+S layer, both the explosion energy and the mass of synthesized Ni monotonically decrease as the neutrino luminosity is reduced. This is consistent with the findings in our previous paper (Yamamoto et al. 2013). This is a consequence of the delay of shock revival for lower neutrino luminosities (see the fourth column in the table), which in turn results in smaller mass accretion rates (the last column of the table).

In Figure 9 the time evolutions of these quantities are presented. We show in addition the energies gained by neutrino heating, E_ν, and produced by recombinations and nuclear burnings, E_nuc, as well as the yields of some major elements as functions of time. The former two are evaluated, respectively, as

$$\begin{eqnarray}&&{E}_{\mathrm{nuc}}=\displaystyle \sum _{{v}_{r}\gt 0,{e}_{\mathrm{tot}}\gt 0}^{}{Q}_{\mathrm{nuc}}^{(n)}\ {\rm{\Delta }}V\ {\rm{\Delta }}{t}^{(n)},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{E}_{\nu }=\displaystyle \sum _{{v}_{r}\gt 0,{e}_{\mathrm{tot}}\gt 0}^{}{Q}_{\nu }^{(n)}\ {\rm{\Delta }}V\ {\rm{\Delta }}{t}^{(n)},\end{eqnarray} \tag{ 15 }$$

where ${\rm{\Delta }}{t}^{(n)}$ is the time interval of the nth time step in the simulations and ${Q}_{\mathrm{nuc}}^{(n)}$ and ${Q}_{\nu }^{(n)}$ are the local energy production rates via the nuclear reactions and neutrino-heating rate, respectively; the summation runs over the same grid points in Equation (12).

Zoom In Zoom Out Reset image size

We first compare the results for different neutrino luminosities, which are shown in Figure 9. In the top panel the temporal changes of the energies given above are displayed for the models with the Fe core mass of M_c = 1.40 M_⊙ and the Si+S layer mass of M_SiS = 0.18 M_⊙. In the bottom panels, on the other hand, the mass yields of ⁵⁶Ni, ⁵⁴Fe, ²⁸Si, ⁴He, and free neutron are presented. It is found from the top panels that lower luminosities indeed delay shock revival and result in smaller explosion energies. The same is true of the nickel mass yield, as is also evident from the bottom panels. The diagnostic explosion energy, E_exp, increases mainly by the nuclear energy releases from the recombinations to ⁴He and heavier elements at first and via the nuclear burnings later. The later slow decline in the diagnostic explosion energy is caused by the engulfing of gravitationally bound envelopes and by some fall backs. Rather small contributions from neutrino heating are deceptive. It should be noted in fact that these energies are calculated only after shock revival. Most of the energy given by neutrino heating is consumed to lift matter up from the gravitational well and make it ready for expansion, which is not included in the evaluation of E_exp, however.

The productions of ⁴He followed by ⁵⁴Fe and then ⁵⁶Ni are due to the recombinations, while ²⁸Si is mainly synthesized by oxygen burnings. The difference in the final yield of ⁵⁴Fe contributes considerably to the difference in E_nuc. It is also seen that not all ⁴He is consumed to assemble heavier elements and the so-called α-rich freeze out takes place. This sequence of events is essentially the same as those we found in our previous work (Yamamoto et al. 2013). It is mentioned that the wind contribution to E_exp is sometimes non-negligible and is the main reason why E_nuc is asymptotically larger than E_exp in some cases and smaller in others. We will come back to this issue in Sections 3.2.4 and 3.3.

It is interesting to point out that the times of shock revival, t_rev, are not very far apart between the models with L_ν = 2.5 × 10⁵² erg and 2.3 × 10⁵² erg, however, the explosion energies are substantially different. This corresponds to what we saw in the critical curve for the models with M_SiS = 0.18 M_⊙ (the lower left panel of Figure 7): the critical luminosities are different for not-so-different mass accretion rates.

The increase in the core mass affects the evolutions of the explosion energy and ejecta mass in a way similar to the decrease in the neutrino luminosity. For instance, for the fixed parameter sets (M_Si+S/ M_⊙, L_ν,52) = (0.18, 2.3), robust explosions with E_exp ≳ 1.0 × 10⁵¹ erg are obtained for smaller cores (M_c = 1.30 and 1.40 M_⊙). The yield of ⁵⁶Ni is also a decreasing function of the core mass. These are common trends shared by other models (except for the model with M_c = 1.30 M_⊙, M_SiS = 0.18 M, and L_ν = 2.5 × 10⁵² erg s⁻¹, in which shock revival occurs very early on (t_exp = 0.37 s), but the explosion energy is a little bit smaller than the model with the heavier core of M_c = 1.40 M_⊙). This may seem to be simply explained by different gravitational attractions by the central PNS but there is something more here, which we will discuss in Section 3.2.3.

The M_SiS dependence is rather complicated. In the case of L_ν,52 = 2.3, the higher M_SiS induces an earlier and more energetic explosion since the growth rate of radial perturbation is larger than that of the mass accretion rate. The behavior in E_nuc reflects the yield of ⁵⁴Fe as in the previous cases. The synthesized ⁵⁶Ni are nearly the same between the models with the different M_SiS. The PNS is more massive for M_SiS = 0.09 M_⊙, reflecting the delay of shock revival in that case. These trends are common to other core mass models (see Table 4). For the lower luminosity case (L_ν = 2.0 ×10⁵² erg s⁻¹), on the other hand, shock revival is still delayed for the lighter Si+S layer but the accretion rate at the shock revival is larger for M_SiS = 0.09 M_⊙ and, as a consequence, the explosion energy and the nickel yield are both greater for these models than those with M_SiS = 0.18 M_⊙. Note that the PNS is still more massive for the delayed shock revival. This happens again because high ${\dot{M}}_{{\rm{rev}}}$ leads to the canonical explosion energy (see more detailed discussions in the next section and Section 3.3). In the higher luminosity models (L_ν = 2.5 × 10⁵² erg s⁻¹), the shock revival time behaves similarly to the first case (L_ν = 2.3 × 10⁵² erg s⁻¹): it is delayed for M_SiS = 0.09 M_⊙. The accretion rates and hence the explosion energies and nickel masses are non-monotonic with the Si+S layer mass. In fact, they are larger for M_SiS = 0.09 M_⊙ than for M_SiS = 0.18 M_⊙ in the models with M_c = 1.30 and 1.40 M_⊙, whereas the opposite is true for M_c = 1.50 M_⊙, which has already explained in Section 3.2.1.

Figure 10 summarizes the explosion energies and nickel masses thus obtained in the E_exp–M_Ni plane with some observational data taken from Pejcha & Prieto (2015a). The upper left panel shows all 18 models, illuminating the clear correlation between them, which is well fitted by

$$\begin{eqnarray}&&{M}_{\mathrm{Ni}}\propto {E}_{\mathrm{exp}}^{0.70}.\end{eqnarray} \tag{ 16 }$$

This is consistent with our previous findings (Yamamoto et al. 2013). It is interesting, however, that the correlation still holds among different M_c and/or M_SiS models. The other three panels display the correlation for different neutrino luminosities separately. It is apparent from these results that as the neutrino luminosity becomes lower, the location of data points in the plane shifts roughly from the upper right to the lower left part on the line that expresses the suggested correlation, which implies that the explosions get weaker and the nickel yield becomes smaller. Looking into the individual cases more closely, we find that the models with L_ν = 2.0 × 10⁵² erg s⁻¹ are more strongly correlated, whereas the models with higher luminosities have larger dispersions. This is mainly due to larger contributions from the neutrino wind in these models. It is also mentioned that the shock revival for L_ν = 2.5 × 10⁵² erg s⁻¹ occurs around the Si/O interface, where $\dot{M}$ changes rapidly.

Zoom In Zoom Out Reset image size

It has been commonly supposed that the typical explosion of CCSNe produces E_exp ≳ 10⁵¹ erg and M_Ni ≲ 0.1 M_⊙ (Hamuy 2003; Nomoto et al. 2006; Kochanek et al. 2008; Utrobin & Chugai 2008; Smartt et al. 2009). This corresponds to the lower right quadrant in each panel. It is evident that not many models, actually only three of them, fall in this region. Recent observations (for example, the papers listed in Table 1 in Pejcha & Prieto 2015b) and their analysis (Pejcha & Prieto 2015a) suggest that the explosion energy and the nickel mass are correlated and occupy the region indicated in the figure. Note that the fitted line expressed by Equation (16) traverses only the upper part of this region. This means that it is difficult in 1D models to obtain the right combination of explosion energy and ⁵⁶Ni mass. As demonstrated later, the disagreement is resolved in the 2D models, which may be yet another support for the claim that multi-dimensionality is crucial for explosion.

3.2.3. Dependences of Explosion Energy on Model Parameters

Now we look into the dependences of the explosion energy on the parameters that characterize the progenitor structures. Although no clear relation is observed between the explosion energy and M_c or M_SiS, a tight correlation is found if we choose the mass accretion rate at shock revival is adopted, which is demonstrated in the left panels of Figure 11. Each panel plots for one of the three neutrino luminosities the results for six progenitor models. The solid line is the least-square fit to all 18 data points. It is rather surprising that all the models fall so closely to the line. The correlation is found to be fitted as

$$\begin{eqnarray}&&{E}_{\mathrm{exp}}\propto {\dot{M}}_{{\rm{rev}}}^{1.61}.\end{eqnarray} \tag{ 17 }$$

Larger $\dot{M}$ are favorable for strong explosions as long as shock revival occurs. In our 1D models, ${\dot{M}}_{{\rm{rev}}}\gtrsim 0.35\;{M}_{\odot }\;{{\rm{s}}}^{-1}$ is required to produce an explosion with E_exp ≳ 10⁵¹ erg.

Zoom In Zoom Out Reset image size

On the right panels of Figure 11, we examine possible correlations between the explosion energy and the final PNS mass. The solid line is again the least-square fit to all data. Although there may be a rough trend that the explosion energy gets smaller as the PNS mass becomes larger, the scatter around the fitting line is not so small. In fact, the following three combinations, $({M}_{c}/{M}_{\odot },{L}_{\nu ,52})$ $\quad =\quad (1.50,2.5)$, $(1.40,2.3),(1.30,2.0)$ for M_SiS = 0.09 M_⊙, give almost the same E_exp although the PNS masses are substantially different among them. They have almost same mass accretion rates at shock revival, on the other hand. It should be noted that although larger PNS masses generate stronger gravitational attractions and are not favorable for shock revival in this respect, they are often a consequence of greater mass accretion rates, which are favorable for robust explosions. These competing effects result in the non-monotonic behavior of the explosion energy versus PNS mass. It is hence better to consider the influence of the PNS mass on E_exp only through $\dot{M}$.

3.2.4. Axisymmetric 2D Explosions

We turn our attention to 2D models, in which we assume axisymmetry. It has been demonstrated by different authors (Dolence et al. 2013; Hanke et al. 2013; Murphy et al. 2013; Couch & O'Connor 2014; Fernández et al. 2014; Takiwaki et al. 2014; Fernández 2015; Müller 2015) that the nature of turbulence is qualitatively different between 2D and 3D and that the shock dynamics is affected accordingly. It has, however, also been shown that the difference between 2D and 3D is much smaller than that between 1D and 2D (Nordhaus et al. 2010; Hanke et al. 2012; Couch 2013). We hence expect that our 2D models will still be able to illuminate how the multi-dimensionality of dynamics modifies our findings so far obtained in the 1D models.

For 2D simulations, we employ only the models with M_c = 1.30 and 1.40 M_⊙, since the heaviest core (M_c = 1.50 M_⊙) models are likely to produce mostly under-energetic explosions except for the highest neutrino luminosity (L_ν = 2.5 × 10⁵² erg s⁻¹) case, which is also not considered in the following.

The critical neutrino luminosities are plotted in the left panels of Figure 12 together with the 1D results presented already in Figure 7. In most of the 2D models, shock revival sets in earlier at higher mass accretion rates than in the corresponding 1D cases. This is due to the hydrodynamical instabilities such as convection and SASI, which expedite shock revival by expanding the gain region and making the dwell time in the region longer, and is consistent with the previous findings by others (Murphy & Burrows 2008; Scheck et al. 2008). In the case of L_ν = 2.0 × 10⁵² erg s⁻¹, for example, the shock revivals occur immediately after the Si/O interface hits the shock front in 2D, since the multi-dimensional effects just mentioned push the post-shock flows almost to the critical points prior to the encounter with the layer boundary and these models do not need to wait further when the interface actually hits the stalled shock wave.

Zoom In Zoom Out Reset image size

The right panels present the mass accretion rates at shock revival together with the final values of the PNS mass just as the right panels of Figure 7. The solid lines are the trajectories for the designated models. The reason why the squares are not sitting on these lines is because the final values instead of the values at shock revival are employed for plotting these data (see also Figure 7). The vertical distance between the trajectory and data point represents the mass gained or lost by the PNS after shock revival. Since the fallback of matter is generally stronger in 2D, M_PNS tends to be heavier in 2D than in 1D.

The time evolutions of the explosion energy as well as the masses of some major elements in the ejecta are shown in Figure 13 for the 2D model with L_ν = 2.3 × 10⁵² erg s⁻¹, M_c = 1.40 M_⊙, and M_SiS = 0.09 M_⊙. For comparison, the corresponding 1D model is exhibited in the right panels. Due to the non-radial fluid instabilities, shock revival occurs much earlier in 2D. The higher mass accretion rate then results in the larger E_exp. It is noteworthy in this context that the curve of E_exp is always running below that of E_nuc in the top right panel for 1D, whereas the opposite is true in the 2D model. The contribution from the neutrino heating rises more slowly and continues longer in 2D. The mass of ⁵⁶Ni in the ejecta is gradually decreasing at late times in 2D. These facts are all consequences of stronger neutrino winds, i.e., outward matter flows from the inner boundary in our simulations. The neutrino wind mainly consists of hot nucleons earlier and α particles later, which will contribute to the explosion energy through recombinations. The addition of these fresh nucleons is the main reason for the more long-lasting neutrino heating in 2D. The later α particle ejection underestimates the neutrino heating and nuclear the recombination energy which should be provided outside the active calculation region (see also Section 3.3).

Zoom In Zoom Out Reset image size

The neutrino wind adds slightly neutron-rich matter to the ejecta, which later recombines to produce ⁵⁴Fe instead of ⁵⁶Ni. This is one of the reasons why the mass of synthesized ⁵⁶Ni is generally smaller in 2D than in 1D. It should be noted, however, that the neutrino wind is one of the most uncertain components in our simulations. In fact, it is actually dictated by the inner boundary condition imposed in our models (see also Yamamoto et al. 2013), although it should have been obtained from the self-consistent evolution of PNS. In fact, previous realistic simulations by other groups, e.g., Scheck et al. (2006), Marek & Janka (2009), and Bruenn et al. (2014), demonstrated that the neutrino-driven wind and its heating by neutrinos may last for more than ∼500 ms, and hence is one of the main contributors to the explosion energy.

We demonstrate in the left panels of Figure 14 that the explosion energy is tightly correlated with the mass accretion rate at shock revival. It is intriguing that the 2D models fall on the same line (see Equation (17)) with the 1D models. On the other hand, the correlation between E_exp and M_PNS seems to become stronger in the 2D cases, as shown in the right panels of Figure 14. It should be mentioned here that the PNS mass is more difficult to estimate in 2D, since the fallback is more violent and highly anisotropic and it takes place simultaneously with the neutrino wind. It will hence be necessary to investigate a larger number of models to justify the above statement.

Zoom In Zoom Out Reset image size

Finally, the explosion energies and masses of synthesized ⁵⁶Ni are plotted in Figure 15 for all models, including the 2D computations. One can easily recognize that the 2D results are qualitatively different from those for 1D unlike the correlation between E_exp and ${\dot{M}}_{{\rm{rev}}}$ demonstrated in Figure 14. In fact, the stronger fallback and neutrino wind in 2D both tend to reduce the mass of ⁵⁶Ni in the ejecta compared with the 1D counterparts and the data points for 2D seem to be fitted by a steeper line (red solid line), although there is a scatter for the low-energy explosions. We showed earlier that the 1D models tend to produce too much nickel if they explode in the first place (see also Figure 10) and that the correlation between the explosion energy and the nickel mass is qualitatively different from that suggested by the observations, which we claimed is yet another reason that 1D explosions are disfavored in the neutrino-heating mechanism. In the 2D case, the reduced nickel production puts both the explosion energy and the Ni mass in the right region (see references in Bruenn et al. 2014; Pejcha & Prieto 2015a). This is consistent with the findings in our previous paper (Yamamoto et al. 2013).

Zoom In Zoom Out Reset image size

From the figure it may be said that the mass of the Si+S layer affects the Ni synthesis more than the Fe core mass. It should be mentioned, however, that the precise estimation of the nickel mass in the ejecta is made difficult not only by the uncertainty in the neutrino wind but also by later interactions of the shock wave with the helium and/or hydrogen layers, which may cause further fallback. This issue will be discussed more in the last section.

We mention the accuracy of our numerical simulations. To assess it, we conducted four more 2D computations with a finer radial or angular mesh for the parameter sets of (M_c/ M_⊙, M_SiS/ M_⊙, L_ν,52) = (1.30, 0.18, 2.3) and (1.40, 0.18, 2.0). We find in computations with a larger number of radial grid points (150 up from 600) that both the explosion energy and the Ni mass are changed by less than 0.5% at most. In the case of the finer angular mesh with 124 grid points in the θ direction, the explosion energy is changed by ∼3%, whereas the Ni mass is affected more with a difference of ∼10%, which is mainly caused by a greater fallback. These numerical errors, though not completely negligible, are not large enough to qualitatively change the conclusions, either.

We have so far observed that the explosion energy, E_exp, is strongly correlated with the mass accretion rate at shock revival, ${\dot{M}}_{\mathrm{rev}}$, irrespective of the dimension of dynamics. In this section, we look into the origin of this correlation more in detail. For that purpose, we introduce a new quantity, the total mass of the matter, ${M}_{{T}_{{\rm{P}}}}$, with the peak temperature higher than T_P > 5 × 10⁹ K. Here the peak temperature is defined for each mass element as the highest temperature it attains during the expansion following shock revival. The threshold temperature, i.e., T_P = 5 × 10⁹ K, is a temperature around which α particles are recombined to iron group elements. ${M}_{{T}_{{\rm{P}}}}$ is hence the mass of the matter that experiences neutrino heating and recombinations, the two main contributors to the explosion energy. It is numerically evaluated as follows:

$$\begin{eqnarray}&&{M}_{{T}_{{\rm{P}}}}=\displaystyle \sum _{{T}_{P,9}^{({\rm{n}})}\gt 5.0,{v}_{r}^{({\rm{n}})}\gt 0}^{}{\rm{\Delta }}{m}^{({\rm{n}})},\end{eqnarray} \tag{ 18 }$$

in which the superscript n labels mass elements on the numerical mesh and T_P,9 = T_P/10⁹ K is a normalized temperature; the condition ${v}_{r}^{(n)}\gt 0$ implies that only expanding matter should be summed.

In Figure 16 we plot all the 1D (18 black crosses) and 2D (8 red crosses) results in the ${\dot{M}}_{\mathrm{rev}}\mbox{--}{M}_{{T}_{{\rm{P}}}}$ (left) and ${M}_{{T}_{{\rm{P}}}}\mbox{--}{E}_{\mathrm{exp}}$ (right) panels. It is evident that they are closely related with one another regardless of the dimension again. The least-square fits to them give the following power-law relations:

$$\begin{eqnarray}&&{M}_{{T}_{{\rm{P}}}}\propto {\dot{M}}_{\mathrm{rev}}^{1.47},\ {E}_{\mathrm{exp}}\propto {M}_{{T}_{{\rm{P}}}}^{1.07}.\end{eqnarray} \tag{ 19 }$$

It is important that the explosion energy is almost linearly correlated with ${M}_{{T}_{{\rm{P}}}}$, the fact suggesting that ${M}_{{T}_{{\rm{P}}}}$ is more directly related with the explosion energy. This will be endorsed by more detailed analyses in the following.

Zoom In Zoom Out Reset image size

We define the explosion energy per nucleon, ε_exp, as well as the contribution to it from the nuclear reactions, ε_nuc, and from the neutrino heating, ε_ν, as follows:

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{exp}}}={E}_{\mathrm{exp}}/{N}_{{T}_{{\rm{P}}}},\ {\varepsilon }_{{\rm{nuc}}}={E}_{{\rm{nuc}}}/{N}_{{T}_{{\rm{P}}}},\ {\varepsilon }_{\nu }={E}_{\nu }/{N}_{{T}_{{\rm{P}}}},\ \end{eqnarray} \tag{ 20 }$$

in which ${N}_{{T}_{{\rm{P}}}}$ is the number of the nucleons that consist of ${M}_{{T}_{{\rm{P}}}}$.

We also refer to the sum of the latter two as ε_src = ε_nuc + ε_ν. We also pay attention below to the total energy per nucleon of the outer envelope, which is defined as

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{bind}}}={E}_{{\rm{bind}}}/{N}_{{T}_{{\rm{P}}}},\end{eqnarray} \tag{ 21 }$$

with E_bind being given by

$$\begin{eqnarray}&&{E}_{{\rm{bind}}}=\displaystyle \sum _{500\;\mathrm{km}\lt {R}^{({\rm{n}})}}^{}{\rm{\Delta }}{V}^{({\rm{n}})}\cdot ({e}_{\mathrm{kin}}^{({\rm{n}})}+{e}_{\mathrm{int}}^{({\rm{n}})}+{e}_{g}^{({\rm{n}})}).\end{eqnarray} \tag{ 22 }$$

Note that E_bind is negative owing to the gravitational binding, although it also contains the kinetic and internal energies, which are positive.

We plot these quantities for all 1D and 2D models in Figure 17, in which the horizontal axis is chosen to be ${M}_{{T}_{{\rm{P}}}}$. First of all, it is clear that they are almost constant, indicating not only the total explosion energy but each contribution except for ε_bind is also nearly proportional to ${M}_{{T}_{{\rm{P}}}}$. If we look more closely, however, we can see some trends. In the left panel, in which the specific explosion energies are plotted, we can recognize that they increase slightly with ${M}_{{T}_{{\rm{P}}}}$, having values of 4.5 MeV ≲ ε_exp ≲ 6.0 MeV, although there are some scatterings. The reason for this increase is understood from the middle panel, in which ε_src and ε_bind are plotted for ${M}_{{T}_{{\rm{P}}}}$. It is found that they both decrease with ${M}_{{T}_{{\rm{P}}}}$. Since the binding energy of the envelope is almost unrelated to ${M}_{{T}_{{\rm{P}}}}$, the decrease of ε_bind is just a consequence of the increase of ${N}_{{T}_{{\rm{P}}}}$, the denominator in Equation (30), with ${M}_{{T}_{{\rm{P}}}}$. The decline in ε_src is much more moderate and, as a result, the total explosion energy per nucleon increases as shown in the left panel.

Zoom In Zoom Out Reset image size

The right panel shows the contents of ε_src, i.e., ε_nuc and ε_ν. It is discernible that both of them decrease slightly, as does ε_src in the middle panel, although 2D results appear to have an opposite trend. This is due to the following reason: as ${M}_{{T}_{{\rm{P}}}}$ increases, the mass of the post-shock matter that has larger α fractions becomes greater, since the photo-dissociations to α occur at larger radii. As a result, the energy gain by later recombinations gets smaller; since the nucleons are mainly responsible for the neutrino heating, ε_ν is also reduced. Note that the values of ε_nuc are ∼5–6 MeV, which are similar to those evaluated at r = 3000 km by Scheck et al. (2006), and are much larger than those of ε_ν, which is at most ∼1–2 MeV. This is partly because that those values are evaluated after shock revival whereas the neutrino heating has started much earlier and is consumed mainly to lift up matter from the gravitational well. It should also be mentioned that the neutrino heating is also used to compensate for the loss by photo-dissociations. In this sense, the recombination energy is actually originated from the neutrino heating. It is interesting to point out that ε_bind is roughly comparable to ε_ν.

The values of ${\varepsilon }_{{\rm{src}}}$ as well as ε_nuc and ε_ν are a bit smaller for 2D than for 1D, a fact which also seems to be the reason for the opposite trend we found in the middle and right panels between the 1D and 2D models. This happens because these quantities are evaluated only in the computational domain, that is, outside the artificial inner boundary. In fact, in 2D models some of the matter flowing into the computational domain through the inner boundary, to which we refer as the wind below, has already recombined to α particles, which then results in the reduction of ε_nuc and ε_ν. This is all the more true for the 2D models with smaller ${M}_{{T}_{{\rm{P}}}}$, producing the opposite trends mentioned above.

The nice thing with the correlation between E_exp and ${M}_{{T}_{{\rm{P}}}}$ is that the latter is fixed earlier than the former. In fact, ${M}_{{T}_{{\rm{P}}}}$ is determined by the time when the shock wave reaches ∼2000–3000 km in our models. It is true that the fallback occurs later after the shock reaches ∼5000 km, but this fallback matter does not have much energy as has been argued in Section 3.2.4. The correlation may hence be useful to estimate the explosion energy. There is one concern here, however. As mentioned above, the wind also contributes to the explosion energy, which is then expressed as

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{exp}}\sim {\varepsilon }_{{\rm{src}}}+{\varepsilon }_{{\rm{bind}}}+{\varepsilon }_{{\rm{wind}}}={\varepsilon }_{{\rm{nuc}}}+{\varepsilon }_{\nu }+{\varepsilon }_{{\rm{bind}}}+{\varepsilon }_{{\rm{wind}}}.\end{eqnarray} \tag{ 23 }$$

In this expression, the wind energy is denoted by ε_wind (see Marek & Janka 2009; Bruenn et al. 2014). Although this is not a dominant contributor in our models, the correlation may be degraded if it is, since we cannot expect that it is strongly correlated with ${M}_{{T}_{{\rm{P}}}}$. Our models have admittedly no predictive power on the issue, since the wind is rather sensitive to the inner boundary condition we impose by hand. As a matter of fact, the wind contribution seems to be substantial in more realistic simulations (Burrows et al. 2007; Marek & Janka 2009; Fischer et al. 2010; Müller et al. 2012), and we are afraid that this may be the reason why the explosion energy is fixed rather early on in our models compared with other realistic simulations (Bruenn et al. 2014; Nakamura et al. 2015; Takiwaki et al. 2014; Müller 2015). Systematic studies will be conducted on the wind in our forthcoming paper.

In the models considered in this paper, there is a clear tendency regardless of spatial dimensions that the explosion energy gets larger as the core mass becomes smaller. This may seem at odds with some recent realistic simulations (Bruenn et al. 2014; Nakamura et al. 2015; Takiwaki et al. 2014; Müller 2015) for light progenitors. In fact, they generically produced weak explosions of E_exp ≲ 0.1 × 10⁵¹ erg, e.g., for model s11.2 of 11.2 M_⊙ by WHW2002. It is true that the time of shock revival is self-consistently determined in these realistic simulations and that the post-shock flow is not steady in reality but we believe that they are not the main cause for the discrepancy mentioned above but the progenitor models are.

As discussed in Section 2.2, the progenitor models can be classified into three groups labeled with "L," "M," and "H." In this paper, we focus on group M, since the majority of progenitor models belongs to this group. We then employ the entropy and Y_e distributions, Equations (1) and (2), suitable for this group. Note that the lowest Fe core mass models are also based on these relations. On the other hand, some of the light progenitor models obtained by realistic evolution calculations are known to be affiliated with group L. Model s11.2 of 11.2 M_⊙ by WHW2002 is indeed one of such models.

To see the difference between groups L and M, we show in Figure 18 the density and entropy distributions for their members. In the left panel, the horizontal axis is the enclosed mass, whereas the radius is adopted in the right panel. Note that only the entropy distributions, which are depicted with dashed lines, in the core as a function of density are different between the two models and the other parameters such as the masses of the Fe core and Si+S layer are identical. It is interesting that the density distributions in the Fe core are more or less the same between the two models. Instead, the difference manifests itself in the outer layers as is particularly evident in the left panel.

Zoom In Zoom Out Reset image size

This may be understood as follows. In both models, the cores are in NSE and are mostly supported by degenerate electrons. As a result, the difference in entropy has a limited influence on the core structure. This is not the case, however, near the core surface as well as in the Si+S layer. Since the entropy is smaller at the core surface for group L than for group M but it is identical in the adjacent Si+S layer for both groups, the entropy jump at the interface is greater in group L. As the densities are nearly the same at the core surface between the two groups, this implies that the pressure is smaller for group L. Since the gravitational attraction is almost identical at the core surface as the mass and radius of the core are nearly the same, the pressure scale height is shorter for group L accordingly. This is reflected in the density scale height in the Si+S layer directly because the entropy is assumed to be identical for both groups there. This is the reason why the density is gradient steeper for group L than for group M as observed in the right panel of Figure 18.

If such an entropy difference in the core is indeed responsible for the discrepancy of the explosion energies mentioned above between the recent realistic simulations and our experimental computations, it may suggest that we would obtain canonical explosions if we could obtain via evolution calculations progenitor models with relatively high entropies in small cores. Detailed analyses of groups L and H will be reported in the forthcoming paper.