REVISITING IMPACTS OF NUCLEAR BURNING FOR REVIVING WEAK SHOCKS IN NEUTRINO-DRIVEN SUPERNOVAE

We solve the hydrodynamic equations corresponding to the conservation of mass, momentum, and energy,

$$\begin{equation} \frac{d \rho }{d t} + \rho \nabla \cdot \mathbf {v} = 0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \rho \frac{d \mathbf {v}}{d t} = - \nabla p - \rho \nabla \Phi, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial e}{\partial t} + \nabla \cdot \left[ (e+p)\mathbf {v} \right] = - p \mathbf {v} \cdot \nabla \Phi + \rho (H - C + Q), \end{equation} \tag{ 3 }$$

where ρ is the mass density, v the fluid velocity, p the pressure, Φ the gravitational potential, and e the total energy density, respectively. The Lagrangian derivative is denoted by $d/dt \equiv \partial / \partial t + {\boldsymbol {v}} \cdot \nabla$. To treat Newtonian self-gravity, a monopole approximation is employed. The tabulated realistic EOS, based on the relativistic mean field theory (Shen et al. 1998), is implemented according to the prescription in Kotake et al. (2003). The term Q in Equation (3) denotes the net energy deposition rate by nuclear burning. The goal of this paper is to explore the effect of this term on shock revival. We compare two cases. (1) For the burning case, we estimate Q by calculating a simple nuclear reaction network, and (2) for the non-burning case, we do not solve the nuclear network (Q = 0), but adopt Shen EOS throughout our simulations.

For the burning case, we are keeping track of 13 species of α network (from ⁴He to ⁵⁶Ni) by solving a separate advection equation for each species. The nuclear reaction network is mainly based on the REACLIB database (Rauscher & Thielemann 2000). Experimentally determined masses (Audi & Wapstra 1995) and reactions (Angulo et al. 1999) are adopted if available. It should be noted that our network does not include the photodissociation of iron elements because Shen EOS adopted in this study takes these endothermic effects into account. It should also be noted that we solve the reaction network only for the grids where T < 5 × 10⁹ K, assuming that the local chemical composition is in nuclear statistical equilibrium above this temperature.

We also employ the LB scheme (Janka & Müller 1996), in which neutrino heating and cooling is adjusted parametrically to trigger explosions. Following Janka (2001) and Nordhaus et al. (2010), the neutrino heating (H) and cooling rates (C) in Equation (3) are given by,

$$\begin{eqnarray} H = 1.544 &\times & 10^{20} \left(\frac{L_{\nu _{\rm e}}}{10^{52} \, {\rm erg \, s^{-1}}}\right) \left(\frac{T_{\nu _{\rm e}}}{4 \, {\rm MeV}} \right)^2 \nonumber \\ & \times & \left(\frac{r}{100 \, {\rm km}} \right)^{-2} (Y_{\rm n} + Y_{\rm p}) \, e^{-\tau _{\nu _{\rm e}}} ({\rm \, erg}\, {\rm g}^{-1}\, {\rm s}^{-1}), \nonumber \\ \end{eqnarray} \tag{ 4 }$$

$$\begin{equation} C = 1.399 \times 10^{20} \left(\frac{T}{2 \, {\rm MeV}} \right)^6 (Y_{\rm n} + Y_{\rm p}) \, e^{-\tau _{\nu _{\rm e}}} ({\rm \, erg} {\rm \, g}^{-1}\, {\rm s}^{-1}), \end{equation} \tag{ 5 }$$

where $L_{\nu _{\rm e}}$ is the electron-neutrino luminosity that is assumed to be equal to the anti-electron neutrino luminosity ($L_{\bar{\nu }_{\rm e}} = L_{\nu _{\rm e}}$), $T_{\nu _{\rm e}}$ is the electron neutrino temperature assumed constant at 4 MeV, r is the distance from the center, T is the local fluid temperature, Y_n and Y_p are the neutron and proton fractions, respectively, and $\tau _{\nu _e}$ is the electron neutrino optical depth, which we estimate according to Equation (7) in Hanke et al. (2012).

In this study, neutrino luminosity is assumed to evolve exponentially with time (Kifonidis et al. 2003) as

$$\begin{equation} L_{\nu _{\rm e}} = L_{\bar{\nu }_{\rm e}} = L_{\nu 0} \, {\rm exp} (- t_{\rm pb} /t_{d}), \end{equation} \tag{ 6 }$$

where L_ν0 denotes the initial luminosity, t_pb is the time measured after core bounce, and t_d is the decay time, and L_ν0 and t_d are treated as free parameters. Neutrino heating and cooling are switched on only after core bounce, according to the prescriptions (Equations (4) and (5)) assuming $L_{\nu _{\rm e}} = L_{\bar{\nu }_{\rm e}}$ and $T_{\nu _{\rm e}} = T_{\bar{\nu }_{\rm e}}$. Before bounce, we employ the Y_e prescription proposed by Liebendörfer (2005), in which Y_e is given simply as a function of density, and after that, we refrain from solving the change of Y_e, following Murphy & Burrows (2008), Nordhaus et al. (2010) and Hanke et al. (2012; see, however, Ohnishi et al. 2006). As for the hydro-solver, we employ the ZEUS-MP code (Hayes et al. 2006), which has been modified for core-collapse simulations (e.g., Iwakami et al. 2008, 2009). The computational grid is comprised of 300 logarithmically spaced, radial zones from the center up to 5000 km. For 2D models we adopt coarse mesh points (n_θ = 32 uniform grids) in the polar direction to make it possible to perform 174 models in 2D covering a wide range of the parameter region. For some selected models a finer resolution (n_θ = 128) is taken.

In order to induce non-spherical instability after the stall of the prompt bounce shock, we added a radial velocity perturbation, δv_r(r, θ, ϕ), to the steady spherically symmetric flow, according to the following equation,

$$\begin{equation} v_r(r, \theta, \phi) = v_r^{0}(r) +\delta v_r(\theta, \phi), \end{equation} \tag{ 7 }$$

with

$$\begin{equation} \delta v_r = 0.01 \times {\rm rnum} \times \, v_r^0 (r,\theta) \,, \end{equation} \tag{ 8 }$$

where $v_r^0 (r,\theta)$ is the unperturbed radial velocity and δv_r is the random multi-mode perturbation with a random number −1 < rnum < 1.

In this study, we employ four progenitor models: three for 15.0 M_☉ stars of Limongi & Chieffi (2006, hereafter LC15), Woosley & Weaver (1995, hereafter WW15), and Woosley et al. (2002, hereafter WHW15) and one for an 11.2 M_☉ star of Woosley et al. (2002, hereafter WHW11). For all the models, Figure 1 shows the pre-collapse composition profiles from near the outer edge of the iron core to the outside. As we will explain in the next section, burning the oxygen shell behind the (weakly propagating) shock plays an important role in assisting the shock expansion. Therefore, the earlier the oxygen layer touches the (stalling) shock after bounce, the better it could work. Among the three variants of the 15 M_☉ progenitors, the inner edge of the oxygen layer (seen as a sharp decline in solid red lines of Figure 1) is positioned much closer to the center for models LC15 (closest, top left panel) and WHW15 (next closest, top right panel) compared to model WW15 (bottom left panel). Table 1 shows a summary of the precollapse abundance distributions, in which each quantity, from the left to right column, corresponds to the different progenitor models, the progenitor mass, the mass of the iron core, the outer edge of the iron core, the mass of the silicon layer, the outer edge of the silicon layer, and the mass of the oxygen layer, respectively. The edge between each layer is defined as the radius where the most abundant element shifts from one to another (see Figure 1). The mass of oxygen layer for the 15 M_☉ models of LC15 and WHW15 (M_O in the table) is larger than the other progenitors (i.e., WW15 and WHW11) and their oxygen layers (denoted by R_Si/O) are positioned much closer to the center, so that they can touch the SN shock in a shorter timescale after bounce (before the neutrino luminosity gets smaller with time). As one would anticipate, the impacts of nuclear burning are most remarkable for the LC15 progenitor as we will show in the later sections.

Zoom In Zoom Out Reset image size

Relying on the LB scheme in this study, the destiny of the stalling bounce shock (i.e., whether or not it will revive) depends on two parameters: the input neutrino luminosity L_ν0 and the decay time t_d (see Equation (6)).⁹ We characterize models as (L_{ν0, 52}, t_d) = (x, y) for convenience, in which the luminosity and the decay time is x × 10⁵² (erg s⁻¹) and y (s), respectively.

Figure 2 shows comparisons of the postbounce shock evolution in 1D LC15 models depending on the two parameters (L_{ν0, 52}, t_d) and nuclear energy feedback from α network calculation. Three sets of parameters were chosen, (L_{ν0, 52}, t_d) = (2.2, 5.0), (2.0, 5.0), and (2.2, 2.0), and are shown as a solid line, a dotted line, and a dash-dotted line, respectively. All the models with nuclear burning, marked with red lines in Figure 2, present a shock expansion leading to explosions, whereas among the three models without nuclear burning (blue lines) only the model with relatively higher luminosity and longer decay time ((L_{ν0, 52}, t_d) = (2.2, 5.0), solid blue line) exhibits a shock revival. Note that in all six cases in Figure 2, the bounce shock stalls and then transits to a passive shock, which presents negative radial velocity behind the shock. Only after the transition, the additional energy gain due to nuclear burning acts to bifurcate the path of the passive shock. Two of three models without burning (blue lines in Figure 2) experiences shock recession afterward, whereas the shock of the burning models (for red lines) expands with different revival timescales depending on the input neutrino parameters.

Zoom In Zoom Out Reset image size

As seen in Figure 2, larger input neutrino luminosity and shorter decay timescale unsurprisingly leads to easier explosions. More importantly, by comparing the dotted red line with the dotted blue line ((L_{ν0, 52}, t_d) = (2.0, 5.0)), the shock is shown to shift from recession to expansion by the inclusion of nuclear burning. In the case of higher neutrino luminosity (L_{ν0, 52} = 2.2), the trajectories of the shock are observed to be rather similar when the effect of nuclear burning is taken into account (compare the solid red with dashed red line).

In the following, we elaborate on how and why the shock expansion is affected by nuclear burning as observed in Figure 2. Figures 3 and 4 show the mass-shell trajectory of models LC15 with the different parameter set of (L_{ν0, 52}, t_d) = (2.2, 2.0), and (L_{ν0, 52}, t_d) = (2.0, 5.0), respectively. The former and latter case corresponds to the dashed and dotted line in Figure 2. Without nuclear burning (left panels in Figures 3 and 4), the stalled shock oscillates but does not turn into expansion (see also Figure 2). With nuclear burning (right panels of Figures 3 and 4), the shock expansion can be seen to take place when the shock front passes through the Si-rich layer (see the behavior of the thick red line in the green region in the right panel of Figure 3) or when it touches the O-rich layer (e.g., the shock in the red region in the right panel of Figure 4). For the latter case, the bounce shock firstly stalls, similar to the non-burning model (compare the left with the right panel in Figure 4), but the shock front deviates from the non-burning case when it encounters the O-rich layer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For more detail about how the nuclear burning contributes to the shock acceleration, Figure 5 shows the radial velocity profiles and the composition distributions for model LC15 with the parameter set (L_{ν0, 52}, t_d) = (2.2, 2.0) (which is the same parameter set as in Figure 3). At t_pb = 150 ms (top left panel), the shock front is in the progenitor Si-rich layer. Behind the shock front, heavier elements are synthesized as shown. The nuclear energy released by silicon burning heats the material behind the shock, giving it a small positive velocity (compare the velocity profiles with and without nuclear burning in the top left panel). The difference between the velocity profiles with versus without nuclear burning becomes outstanding when the O-rich layer starts to touch the shock front (t_pb > 250 ms). Figure 6 shows the evolution of the diagnostic (explosion) energy for burning (red line) and non-burning (blue line) cases as well as the net energy released via nuclear reactions (green line). As in Suwa et al. (2010), we define the diagnostic energy that refers to the integral of the energy over all zones that have a positive sum of the specific internal, kinetic, and gravitational energy. It is impossible to calculate the final energy of an explosion that is still occurring at this early post-bounce stage. After silicon burning starts to feed the energy behind the shock, in addition to neutrino heating (in the gain region, e.g., t_pb = 150 ms, see the top left panel in Figure 5), the diagnostic energy deviates from that without burning (compare the red and blue lines of Figure 6), which is also clearly visible in the shock evolution (Figure 3). From Figure 6, the total amount of 3.1 × 10⁵⁰ erg is shown to be released through nuclear burning in this case, raising the diagnostic energy to 5.0 × 10⁵⁰ erg.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As we already mentioned, oxygen burning predominantly triggers the shock expansion for the parameter set taken for Figure 4. In addition, the silicon layer is shown to be burned as a heating source (top left panel of Figure 7), which is the reason that the shock position becomes larger compared to the non-burning model (Figure 4). When the shock front begins to swallow the oxygen layer at ∼400 ms postbounce (the right panel of Figure 4), the fresh fuel supplies energy to assist the shock expansion (see from top right, bottom left, to bottom right panels of Figure 7). If not for the energy gain, the stalled shock would not revive earlier than t = 750 ms, as seen from the left panel of Figure 4. Even with the aid of nuclear burning, the explosion for this model (Figure 8) is weaker (≲ 10⁵⁰ erg) compared to the more luminous models (Figures 6 and 9). This suggests that nuclear burning has a secondary impact on the explosion mechanism—it can assist explosions only when the neutrino heating is working enough to push the weak shock to the fuel layers.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10 shows a parameter map on the (L_{ν0, 52}, t_d) plane for each progenitor. As can be seen, higher neutrino luminosity and/or longer decay time leads to easier explosions (colored by red and denoted as "Explosion"), whereas it is opposite for smaller neutrino luminosity and/or shorter decay time (colored by light-blue and denoted as "No explosion" in the figure). For the LC15 progenitor model (top left panel), a parameter region colored by yellow can be seen between the exploding and the non-exploding regime, in which an explosion is obtained only when nuclear burning is included in the hydrodynamics simulations. The emergence of the yellow region means that the minimum neutrino luminosity necessary to drive an explosion is reduced by taking into account energy feedback from nuclear burning. The burning-mediated regime is clearly visible only for the LC15 progenitor. As already mentioned in Section 2.2, this is because this model possesses a massive oxygen layer and the oxygen shell is positioned closest to the center among the progenitors taken in this study.

Zoom In Zoom Out Reset image size

The area of the yellow region in Figure 10 is not so large even for the LC15 progenitor, which suggests that nuclear burning has a secondary importance. In the case of energetic explosions, for example (L_{ν0, 52}, t_d) = (3.0, 0.8), a difference of diagnostic energy between models with and without nuclear burning is ∼5 × 10⁴⁹ erg (Table 2). For marginal weak explosions with E_dia ≲ 10⁵⁰ erg (which are often the case in recent first-principle CCSN simulations), however, it should be emphasized that the inclusion of nuclear burning could increase the diagnostic energy up to about ∼0.6 × 10⁵¹ erg.

The critical luminosity for explosions can be read from the y-axis in Figure 10 in the limit of long t_d (namely, approaching to a constant neutrino luminosity), which corresponds to 2.7 (WHW15), 2.0 (LC15), 1.9 (WW15), and 1.2 (WHW11) in unit of 10⁵² erg s⁻¹, respectively. The WHW15 model shows the highest critical neutrino luminosity among our models. This is because of the higher mass accretion rate of the WHW15 model (blue solid line in Figure 11), which makes the impact of nuclear burning relatively small. The critical luminosity becomes smallest for model WHW11, mainly owing to the compactness of the pre-collapse core and small mass accretion rate coming from its tenuous envelope, as shown in a blue dotted line in Figure 11. The mass accretion rate, averaged between 200 ms and 600 ms after bounce for each model, is 0.33 (WHW15), 0.23 (LC15), 0.21 (WW15), and 0.08 M_☉ s⁻¹ (WHW11), respectively, which is roughly proportional to the critical luminosity except for the WHW11 progenitor model. When the input luminosity is taken below the critical curves, nuclear burning cannot drive explosions alone because the shock needs to expand far away from the central protoneutron star by neutrino heating (i.e., the shock revival due to neutrino heating is preconditioned to enjoy the assistance from nuclear burning).

In Table 2, the time of explosion, t_exp, and diagnostic energy, E_dia, are listed for some chosen sets of the neutrino parameters. The time of explosion is defined as the moment when a shock reaches an average radius of 4500 km, while non-exploding models are denoted by a "—" symbol. The diagnostic energy is plotted as a function of neutrino luminosity for a various decay time scale in Figure 12. In fact, the diagnostic energy is shown to be remarkably enhanced in the case of marginal explosions (i.e., low neutrino luminosity L_{ν, 0} and/or short decay time t_d), and the difference gets smaller for large L_{ν, 0} and t_d, in which explosions are predominantly triggered by neutrino heating. As previously mentioned, these features, which are due to nuclear burning, are only remarkable in the LC15 progenitor. Nevertheless, it is worth mentioning that even for the WW95 progenitor the shock extent, for which the impact of nuclear burning is relatively small (see the bottom left panel of Figure 10), becomes bigger for models with nuclear burning compared to those without (Figure 13).

We move on to discuss axi-symmetric 2D models and examine the effects of nuclear burning in the same manner as in the previous section. To clearly see the impacts of nuclear burning in our 2D simulations, we choose to employ the LC15 in the following.

Figures 14 and 15 show entropy evolution (left-hand side in each panel) with the mass fraction of silicon (right-hand side) for two sets of neutrino parameters at selected postbounce epochs (t_pb = 100 and 200 ms postbounce). Small- and large-scale inhomogeneities in the entropy plots come from neutrino-driven convection and the SASI, both of which lead to more easier explosions in 2D than 1D (e.g., Marek & Janka 2009; Murphy & Burrows 2008; Ohnishi et al. 2006).

Reflecting the stochastic motions of the expanding shocks, the way the (anisotropic) shock surfaces touch the nuclear fuel (in the shape of spherical shells) changes from model to model in 2D. In the case of L_{ν0, 52} = 2.4 (Figure 14), the expanding shock first reaches to the silicon layer near in the vicinity of the north pole at t ∼ 150 ms. Simultaneously, heavy elements like nickel are synthesized, which helps to push the burning material preferentially along the direction. In a less luminous case, assuming smaller luminosity (L_{ν0, 52} = 2.2) and the same decay time (Figure 15), the shock encounters the silicon layer closer to the center, where a mass accretion rate is effectively higher, resulting in the longer explosion time.

As shown from Figure 16, nuclear burning does assist 2D explosions similar to 1D, but the energy difference (here ∼0.2 × 10⁵¹ erg) is generally smaller in 2D than in 1D (compare with Figures 6 and 9). The comparison of the energy gain due to nuclear burning between 1D and 2D models is more clearly shown in Table 2 and Figure 17. The difference of the diagnostic energy with and without nuclear burning is larger in 1D than in 2D. This may be because neutrino-driven convection and the SASI in 2D models (as indicated by entropy distributions in Figures 14 and 15) enhances the neutrino heating efficiency, which makes the contribution of nuclear burning relatively smaller compared to 1D models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally Figure 18 is the parameter map in 2D for the LC15 progenitor. As expected, 2D hydrodynamics leads to easier explosions compared to 1D (see the dashed lines that are the critical curves in 1D). More importantly, the yellow region still exists in 2D models for the LC15 progenitor. Nuclear energy released in 2D models reduces the critical luminosity by 0.1–0.5 × 10⁵² erg s⁻¹ depending on t_d, as well as 1D models, although its impact on the diagnostic energy is weaker than for 1D models. It would be interesting to perform multi-D (radiation-hydro) simulations with a nuclear network calculation for the previously untested progenitor model.