PUSHING CORE-COLLAPSE SUPERNOVAE TO EXPLOSIONS IN SPHERICAL SYMMETRY. I. THE MODEL AND THE CASE OF SN 1987A

We make use of the general relativistic hydrodynamics code AGILE in spherical symmetry (Liebendörfer et al. 2001). For the stellar collapse, we apply the deleptonization scheme of Liebendörfer (2005). For the neutrino transport, we employ the isotropic diffusion source approximation (IDSA) for the electron neutrinos ${\nu }_{e}$ and electron anti-neutrinos ${\bar{\nu }}_{e}$ (Liebendörfer et al. 2009), and an advanced spectral leakage scheme (ASL) for the heavy-lepton flavor neutrinos ${\nu }_{x}={\nu }_{\mu },{\bar{\nu }}_{\mu },{\nu }_{\tau },{\bar{\nu }}_{\tau }$ (Perego et al. 2014). We discretize the neutrino energy using 20 geometrically increasing energy bins, in the range $3\;\mathrm{MeV}\leqslant {E}_{\nu }\leqslant 300\;\mathrm{MeV}$. The neutrino reactions included in the IDSA and ASL scheme are summarized in Table 1. They represent the minimal set of the most relevant weak processes in the post-bounce phase, particularly up to the onset of an explosion. Note that electron captures on heavy nuclei and neutrino scattering on electrons, which are relevant in the collapse phase (see, for example, Mezzacappa & Bruenn 1993c, 1993b), are not included explicitly in the form of reaction rates, but as part of the parameterized deleptonization scheme. In the ASL scheme, we omit nucleon–nucleon bremsstrahlung, $N+N\leftrightarrow N+N+{\nu }_{x}+{\bar{\nu }}_{x}$, where N denotes any nucleon (see, for example, Hannestad & Raffelt 1998; Bartl et al. 2014). We have observed that its inclusion would overestimate μ and τ neutrino luminosities during the PNS cooling phase due to the missing neutrino thermalization provided by inelastic scattering on electrons and positrons at the PNS surface. However, we have also tested that the omission of this process does not significantly change the μ and τ neutrino luminosities predicted by the ASL scheme before the explosion sets in. Thus, neglecting N–N bremsstrahlung is only relevant for the cooling phase, where it improves the overall behavior when compared to simulations obtained with detailed Boltzmann neutrino transport (e.g., Fischer et al. 2010).

Table 1.  Relevant Neutrino Reactions

Reactions	Treatment	Reference
${e}^{-}+p\leftrightarrow n+{\nu }_{e}$	IDSA	a
${e}^{+}+n\leftrightarrow p+{\bar{\nu }}_{e}$	IDSA	a
$N+\nu \leftrightarrow N+\nu$	IDSA & ASL	a
$(A,Z)+\nu \leftrightarrow (A,Z)+\nu$	IDSA & ASL	a
${e}^{-}+{e}^{+}\leftrightarrow {\nu }_{\mu ,\tau }+{\bar{\nu }}_{\mu ,\tau }$	ASL	a, b

Note. Nucleons are denoted by N. The nucleon charged current rates are based on Bruenn (1985), but effects of mean-field interactions (Reddy et al. 1998; Martínez-Pinedo et al. 2012; Roberts et al. 2012; Hempel 2015) are taken into account.

References. (a) Bruenn (1985); (b) Mezzacappa & Bruenn (1993a).

Download table as:  ASCIITypeset image

The equation of state (EOS) of Hempel & Schaffner-Bielich (2010; HS) that we are using includes various light nuclei, such as alphas, deuterons, or tritons (details below). However, the inclusion of all neutrino reactions for this detailed nuclear composition would be beyond the standard approach implemented in current supernova simulations, where only scattering on alpha particles is typically included. To not completely neglect the contributions of the other light nuclei, we have added their mass fractions to the unbound nucleons. This is motivated by their very weak binding energies and, therefore, by the idea that they behave similarly as the unbound nucleon component.

In all our models we use 180 radial zones which include the progenitor star up to the helium shell. This corresponds to a radius of R ≈ (1.3–1.5) × 10¹⁰ cm. With this setup, we model the collapse, bounce, and onset of the explosion. The grid of AGILE is adaptive, with more resolution where the gradients of the thermodynamic variables are steeper. Thus, in the post-bounce and explosion phases, the surface of the PNS and the shock are the better resolved regions. The simulations are run for a total time of 5 s, corresponding to $\gtrsim 4.6$ s after core bounce. At this time, the shock has not yet reached the external edge of our computational domain.

For the high-density plasma in nuclear statistical equilibrium (NSE) the tabulated microphysical EOS HS(DD2) is used. This supernova EOS is based on the model of Hempel & Schaffner-Bielich (2010). It uses the DD2 parametrization for the nucleon interactions (Typel et al. 2010), the nuclear masses from Audi et al. (2003), and the Finite Range Droplet Model (Möller et al. 1995a). 8140 nuclei are included in total, up to Z = 136 and to the neutron drip line. The HS(DD2) was first introduced in Fischer et al. (2014), where its characteristic properties were discussed and general EOS effects in core-collapse supernova simulations were investigated. Fischer et al. (2014) showed that the HS(DD2) EOS gives a better agreement with constraints from nuclear experiments and astrophysical observations than the commonly used EOSs of Lattimer & Swesty (1991) and Shen et al. (1998). Furthermore, additional degrees of freedom, such as various light nuclei and a statistical ensemble of heavy nuclei, are taken into account. The nucleon mean-field potentials, which are used in the charged-current rates, have been calculated consistently (Hempel 2015). The maximum mass of a cold neutron star for the HS(DD2) EOS is 2.42 ${{{\rm M}}}_{\odot }$ (Fischer et al. 2014), which is well above the limits from Demorest et al. (2010) and Antoniadis et al. (2013).

The EOS employed in our simulations includes an extension to non-NSE conditions. In the non-NSE regime the nuclear composition is described by 25 representative nuclei from neutrons and protons to iron-group nuclei. The chosen nuclei are the alpha-nuclei ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti, ⁴⁸Cr, ⁵²Fe, ⁵⁶Ni, complemented by ¹⁴N and the following asymmetric isotopes: ³He, ³⁶S, ⁵⁰Ti, ⁵⁴Fe, ⁵⁶Fe, ⁵⁸Fe, ⁶⁰Fe, ⁶²Fe, and ⁶²Ni. With these nuclei it is possible to achieve a mapping of the abundances from the progenitor calculations onto our simulations which is consistent with the provided electron fraction, i.e., mantaining charge neutrality. All the nuclear masses M_i are taken from Audi et al. (2003). To advect the nuclear composition inside the adaptive grid, we implement the Consistent Multi-fluid Advection method by Plewa and Müller (1999). For given abundances, the non-NSE EOS is calculated based on the same underlying description used in the NSE regime (Hempel & Schaffner-Bielich 2010), but with the following modifications: excited states of nuclei are neglected, excluded volume effects are not taken into account, and the nucleons are treated as non-interacting Maxwell–Boltzmann gases. Such a consistent description of the non-NSE and NSE phases prevents spurious effects at the transitions between the two regimes.

Outside of NSE, an approximate α-network is used to follow the changes in composition. Explosive Helium, Carbon, Neon, and Oxygen burning are currently implemented in the simulation. Note that the thermal energy generation by the nuclear reactions is fully incorporated via the detailed non-NSE treatment. We do not have to calculate explicitly any energy liberation, but just the changes in the abundances. Within our relativistic treatment of the EOS (applied both in the non-NSE and NSE regime) energy conservation means that the specific internal energy ${e}_{\mathrm{int}}$ is not changed by nuclear reactions. This is due to fact that ${e}_{\mathrm{int}}$ includes the specific rest mass energy ${e}_{\mathrm{mass}}$, where ${e}_{\mathrm{mass}}$ is given by the sum over the masses of all nuclei weighted with their yield ${Y}_{i}={X}_{i}/{A}_{i}$,

$$\begin{eqnarray}&&{e}_{\mathrm{mass}}=\displaystyle \sum _{i}{Y}_{i}{M}_{i}.\end{eqnarray} \tag{ 1 }$$

However, if we define the specific thermal energy ${e}_{\mathrm{th}}$ as

$$\begin{eqnarray}&&{e}_{\mathrm{th}}={e}_{\mathrm{int}}-{e}_{\mathrm{mass}},\end{eqnarray} \tag{ 2 }$$

the nuclear reactions will decrease the rest mass energy (i.e., increase the binding) and consequently increase the thermal energy. This treatment of the non-NSE EOS is consistent with the convention used in all high-density NSE EOSs.

Due to limitations of our approximate α-network, and because we do not include any quasi-statistical equilibrium description, we apply some parameterized burning for temperatures between 0.3 and 0.4 MeV. A temperature dependent burning timescale is introduced, which gradually transforms the initial non-NSE composition toward NSE. For temperatures of 0.4 MeV and above, the non-NSE phase always reaches a composition close to NSE and its thermodynamic properties become very similar to the NSE phase of the HS(DD2) EOS. Even though the two phases are based on the same input physics, small, but unavoidable differences can remain, due to the limited set of nuclei considered in non-NSE. To assure a smooth transition for all conditions, we have introduced a transition region as an additional means of thermodynamic stability. We have chosen a parameterization in terms of temperature and implement a linear transition in the temperature interval from 0.40 to 0.44 MeV. We check that the basic thermodynamic stability condition ${ds}/{dT}\gt 0$ is always fulfilled.

For this study, we use solar-metallicity, non-rotating stellar models from the stellar evolution code KEPLER (Woosley et al. 2002). Our set includes 16 pre-explosion models with zero-age main sequence (ZAMS) mass between 18.0 ${{{\rm M}}}_{\odot }$ and 21.0 ${{{\rm M}}}_{\odot }$ in increments of 0.2 ${{{\rm M}}}_{\odot }$. These models have been selected to have ZAMS mass around 20 ${{{\rm M}}}_{\odot }$, similar to the progenitor of SN 1987A (e.g., Podsiadlowski et al. 2007). We label the models by their ZAMS mass. In Figure 1, the density profiles of the progenitor models are shown. For each of them the compactness parameter ${\xi }_{M}$ is defined following O'Connor & Ott (2011) by the ratio of a given mass M and the radius R(M) which encloses this mass:

$$\begin{eqnarray}&&{\xi }_{M}\equiv \displaystyle \frac{M/{{{\rm M}}}_{\odot }}{R(M)/1000\;\mathrm{km}}.\end{eqnarray} \tag{ 3 }$$

Typically, either ${\xi }_{1.75}$ or ${\xi }_{2.5}$ are used. The compactness can be computed at the onset of collapse or at bounce, as suggested by O'Connor & Ott (2011). For our progenitors, the difference in the compactness parameter between these two moments is not significant for our discussions. Thus, for simplicity, in the following we will use ${\xi }_{1.75}$ computed at the onset of the collapse. The progenitor models considered here fall into two distinct families of compactness: low compactness (${\xi }_{1.75}\lt 0.45;$ LC models) and high compactness (${\xi }_{1.75}\gt 0.45$; HC models), see Table 2. Figure 2 shows the compactness as function of ZAMS mass for the progenitors of this study. The non-monotonous behavior is a result of the evolution before collapse. The mass range between 19 and 21 ${{{\rm M}}}_{\odot }$ is particularly prone to variations of the compactness. For a detailed discussion of the behavior of the compactness as function of ZAMS mass see Sukhbold & Woosley (2014).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Progenitor Properties

${M}_{\mathrm{ZAMS}}$	${\xi }_{1.75}$		${\xi }_{2.5}$		${M}_{\mathrm{prog}}$	${M}_{\mathrm{Fe}}$	${M}_{\mathrm{CO}}$	${M}_{\mathrm{He}}$	${M}_{\mathrm{env}}$
(${{{\rm M}}}_{\odot }$)	At collapse	At bounce	At collapse	At bounce	(${{{\rm M}}}_{\odot }$)	(${{{\rm M}}}_{\odot }$)	(${{{\rm M}}}_{\odot }$)	(${{{\rm M}}}_{\odot }$)	(${{{\rm M}}}_{\odot }$)
18.2	0.37	0.380	0.173	0.173	14.58	1.399	4.174	5.395	9.186
18.6	0.365	0.375	0.170	0.170	14.85	1.407	4.317	5.540	9.313
18.8	0.357	0.366	0.166	0.166	15.05	1.399	4.390	5.613	9.435
19.6	0.282	0.288	0.118	0.117	13.37	1.461	4.959	6.243	7.125
19.8	0.334	0.341	0.135	0.135	14.54	1.438	4.867	6.112	8.428
20.0	0.283	0.287	0.125	0.125	14.73	1.456	4.960	6.215	8.517
20.2	0.238	0.241	0.104	0.104	14.47	1.458	5.069	6.342	8.125
18.0	0.463	0.485	0.199	0.199	14.50	1.384	4.104	5.314	9.187
18.4	0.634	0.741	0.185	0.185	14.82	1.490	4.238	5.459	9.366
19.0	0.607	0.715	0.191	0.191	15.03	1.580	4.461	5.693	9.341
19.2	0.633	0.737	0.191	0.192	15.08	1.481	4.545	5.760	9.325
19.4	0.501	0.535	0.185	0.185	15.22	1.367	4.626	5.860	9.365
20.4	0.532	0.594	0.192	0.192	14.81	1.500	5.106	6.376	8.433
20.6	0.742	0.95	0.278	0.279	14.03	1.540	5.260	6.579	7.450
20.8	0.726	0.904	0.271	0.272	14.34	1.528	5.296	6.609	7.735
21.0	0.654	0.764	0.211	0.212	13.00	1.454	5.571	6.969	6.026

Note. ZAMS mass, compactness ${\xi }_{1.75}$ and ${\xi }_{2.5}$ at the onset of collapse and at bounce, total progenitor mass at collapse (${M}_{\mathrm{prog}}$), mass of the iron core (${M}_{\mathrm{Fe}}$), carbon–oxygen core (${M}_{\mathrm{CO}}$), helium core (${M}_{\mathrm{He}}$), and mass of the hydrogen-rich envelope (${M}_{\mathrm{env}}$) at collapse for all the progenitor models included in this study. The top part of the table includes the low-compactness progenitors (LC; ${\xi }_{1.75}\lt 0.4$ at collapse), the bottom part includes the high-compactness progenitors (HC; ${\xi }_{1.75}\gt 0.45$ at collapse).

Download table as:  ASCIITypeset image

2.4.1. Rationale

2.4.2. Implementation

The additional energy deposition, that represents the main feature of PUSH, is achieved by introducing a local heating term, ${Q}_{\mathrm{push}}^{+}(t,R)$ (energy per unit mass and time), given by

$$\begin{eqnarray}&&{Q}_{\mathrm{push}}^{+}(t,r)=4\;{\mathcal{G}}(t){\displaystyle \int }_{0}^{\infty }{q}_{\mathrm{push}}^{+}(r,E)\;{dE},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{q}_{\mathrm{push}}^{+}(r,E)\equiv {\sigma }_{0}\ \displaystyle \frac{1}{4\;{m}_{b}}{\left(\displaystyle \frac{E}{{m}_{e}{c}^{2}}\right)}^{2}\displaystyle \frac{1}{4\pi {r}^{2}}\left(\displaystyle \frac{{{dL}}_{{\nu }_{x}}}{{dE}}\right){\mathcal{F}}(r,E),\end{eqnarray} \tag{ 5 }$$

with

$$\begin{eqnarray}&&{\sigma }_{0}=\displaystyle \frac{4{G}_{F}^{2}{\left({m}_{e}{c}^{2}\right)}^{2}}{\pi {({\hslash }c)}^{4}}\approx 1.759\times {10}^{-44}\;{\mathrm{cm}}^{2}\end{eqnarray} \tag{ 6 }$$

being the typical neutrino cross-section, ${m}_{b}\approx 1.674\times {10}^{-24}{{\rm g}}$ an average baryon mass, and $({{dL}}_{{\nu }_{x}}/{dE})/(4\pi {r}^{2})$ the spectral energy flux for any single ${\nu }_{x}$ neutrino species with energy E. Note that all four heavy neutrino flavors are treated identically by the ASL scheme, and contribute equally to ${Q}_{\mathrm{push}}^{+}$ (see the factor 4 appearing in Equation (4)).

The term ${\mathcal{F}}(r,E)$ in Equation (5) defines the spatial location where ${Q}_{\mathrm{push}}^{+}(t,r)$ is active:

$$\begin{eqnarray}{\mathcal{F}}(r,E)=\;\left\{\begin{array}{ll}0 & {{\rm if}}\quad {ds}/{dr}\gt 0\quad {{\rm or}}\quad {\dot{e}}_{{\nu }_{e},{\bar{\nu }}_{e}}\lt 0\\ \mathrm{exp}(-{\tau }_{{\nu }_{e}}(r,E)) & {{\rm otherwise}}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

where ${\tau }_{{\nu }_{e}}$ denotes the (radial) optical depth of the electron neutrinos, s is the matter entropy, and ${\dot{e}}_{{\nu }_{e},{\bar{\nu }}_{e}}$ the net specific energy rate due to electron neutrinos and anti-neutrinos. The two criteria above are a crucial ingredient in our description of triggering CCSN explosions: PUSH is only active where electron neutrinos are heating (${\dot{e}}_{{\nu }_{e},{\bar{\nu }}_{e}}\gt 0$) and where neutrino-driven convection can occur (${ds}/{dr}\lt 0$).

The term ${\mathcal{G}}(t)$ in Equation (4) determines the temporal behavior of ${Q}_{\mathrm{push}}^{+}(t,r)$. Its expression reads

$$\begin{eqnarray}{\mathcal{G}}(t)={k}_{\mathrm{push}}\times \left\{\begin{array}{ll}0 & t\leqslant {t}_{\mathrm{on}}\\ \left(t-{t}_{\mathrm{on}}\right)/{t}_{\mathrm{rise}} & {t}_{\mathrm{on}}\lt t\leqslant {t}_{\mathrm{on}}+{t}_{\mathrm{rise}}\\ 1 & {t}_{\mathrm{on}}+{t}_{\mathrm{rise}}\lt t\leqslant {t}_{\mathrm{off}}\\ \left({t}_{\mathrm{off}}+{t}_{\mathrm{rise}}-t\right)/{t}_{\mathrm{rise}} & {t}_{\mathrm{off}}\lt t\leqslant {t}_{\mathrm{off}}+{t}_{\mathrm{rise}}\\ 0 & t\gt {t}_{\mathrm{off}}+{t}_{\mathrm{rise}}\end{array}\right.,\end{eqnarray} \tag{ 8 }$$

and it is sketched in Figure 3. Note that throughout the article we always measure the time relative to bounce, if not noted otherwise.

Zoom In Zoom Out Reset image size

The cumulative energy deposited by PUSH, ${E}_{\mathrm{push}}$, can be calculated from the energy deposition rate ${{dE}}_{\mathrm{push}}/{dt}$ as

$$\begin{eqnarray}&&{E}_{\mathrm{push}}(t)={\displaystyle \int }_{{t}_{\mathrm{on}}}^{t}\left(\displaystyle \frac{{{dE}}_{\mathrm{push}}}{{dt}}\right)\;{dt}\prime ={\displaystyle \int }_{{t}_{\mathrm{on}}}^{t}\left({\displaystyle \int }_{{V}_{\mathrm{gain}}}{Q}_{\mathrm{push}}^{+}\;\rho \;{dV}\right){dt}\prime .\end{eqnarray} \tag{ 9 }$$

where ${V}_{\mathrm{gain}}$ is the volume of the gain region. Both these quantities have to be distinguished from the corresponding energy and energy rate obtained by IDSA:

$$\begin{eqnarray}&&{E}_{\mathrm{idsa}}(t)={\displaystyle \int }_{0}^{t}\left(\displaystyle \frac{{{dE}}_{\mathrm{idsa}}}{{dt}}\right){dt}\prime ={\displaystyle \int }_{0}^{t}\left({\displaystyle \int }_{{V}_{\mathrm{gain}}}{\dot{e}}_{{\nu }_{e},{\bar{\nu }}_{e}}\;\rho \;{dV}\right){dt}.\end{eqnarray} \tag{ 10 }$$

The definition of ${\mathcal{G}}(t)$ introduces a set of (potentially) free parameters:

• 1.  
  ${k}_{\mathrm{push}}$ is a global multiplication factor that controls directly the amount of extra heating provided by PUSH. The choices of ${\sigma }_{0}$ as reference cross-section and of the μ and τ neutrino luminosity as energy reservoir suggest ${k}_{\mathrm{push}}\gtrsim 1$.
• 2.  
  ${t}_{\mathrm{on}}$ sets the time at which PUSH starts to act. We relate ${t}_{\mathrm{on}}$ to the time when deviations from spherically symmetric behavior appear in multi-dimensional models. Matter convection in the gain region sets in once the advection time scale ${\tau }_{\mathrm{adv}}$ and the convective growth time scale ${\tau }_{\mathrm{conv}}$ satisfy ${\tau }_{\mathrm{adv}}/{\tau }_{\mathrm{conv}}\gtrsim 3$ (Foglizzo et al. 2006). For all the models we have explored, this happens around t = 0.06–0.08 s. In the above estimates, ${\tau }_{\mathrm{adv}}={\dot{M}}_{\mathrm{shock}}/{M}_{\mathrm{gain}}$, where ${\dot{M}}_{\mathrm{shock}}$ is the accretion rate at the shock and ${M}_{\mathrm{gain}}$ the mass in the gain region, and ${\tau }_{\mathrm{conv}}={f}_{{{\rm B}}-{{\rm V}}}^{-1}$, where ${f}_{{{\rm B}}-{{\rm V}}}$ is the Brunt–Väisäla frequency. Considering that ${\tau }_{\mathrm{conv}}\sim 4{{\rm -}}5\;\mathrm{ms}$, we expect ${t}_{\mathrm{on}}\sim 0.08{{\rm -}}0.10\;{{\rm s}}$, in agreement with recent multi-dimensional simulations (see, for example, the 1D-2D comparison of the shock position in Bruenn et al. 2013).
• 3.  
  ${t}_{\mathrm{rise}}$ defines the time scale over which ${\mathcal{G}}(t)$ increases from zero to k_push. We connect ${t}_{\mathrm{rise}}$ with the time scale that characterizes the growth of the largest multi-dimensional perturbations between the shock radius (${R}_{\mathrm{shock}}$) and the gain radius (${R}_{\mathrm{gain}}$) (e.g., Janka & Mueller 1996). Foglizzo et al. 2006 showed that convection in the gain region can be significantly stabilized by advection, especially if ${\tau }_{\mathrm{adv}}/{\tau }_{\mathrm{conv}}$ only marginally exceeds 3, and that the growth rate of the fastest growing mode is diminished. Thus, ${t}_{\mathrm{rise}}\gg {\tau }_{\mathrm{conv}}$. On the other hand, a lower limit to ${t}_{\mathrm{rise}}$ is represented by the overturn time scale, ${\tau }_{\mathrm{overturn}}$, defined as

  $$\begin{eqnarray}&&{\tau }_{\mathrm{overturn}}\sim \displaystyle \frac{\pi ({R}_{\mathrm{shock}}-{R}_{\mathrm{gain}})}{{\langle v\rangle }_{\mathrm{gain}}},\end{eqnarray} \tag{ 11 }$$

  where $\langle v{\rangle }_{\mathrm{gain}}$ is the average fluid velocity inside the gain region. In our simulations, we have found ${\tau }_{\mathrm{overturn}}\approx 0.05\;{{\rm s}}$ around and after ${t}_{\mathrm{on}}$. In case of a contracting shock, SASIs are also expected to develop around $0.2{{\rm -}}0.3\;{{\rm s}}$ after bounce (Hanke et al. 2013). Hence, we assume $0.05\;{{\rm s}}\lesssim {t}_{\mathrm{rise}}\lesssim \left(0.30\;{{\rm s}}-{t}_{\mathrm{on}}\right)$.
• 4.  
  ${t}_{\mathrm{off}}$ sets the time after which PUSH starts to be switched off. We expect neutrino driven explosions to develop for $t\lesssim 1\;{{\rm s}}$ due to the fast decrease of the luminosities during the first seconds after core bounce. Hence, we fix ${t}_{\mathrm{off}}=1\;{{\rm s}}.$ PUSH is not switched off suddenly at the onset of the explosion, but rather starts decreasing naturally even before 1 s after core bounce due to the decreasing neutrino luminosities and due to the rarefaction of the gain region above the PNS. The subsequent injection of energy by neutrinos in the accelerating shock is qualitatively consistent with multi-dimensional simulations, where accretion and explosion can coexist during the early stages of the shock expansion. The decrease of ${{dE}}_{\mathrm{push}}/{dt}$ on a time scale of a few hundreds of milliseconds after the explosion has been launched makes our results largely independent of the choice of ${t}_{\mathrm{off}}$ for explosions happening not too close to it.

While ${t}_{\mathrm{on}}$ is relatively well constrained and ${t}_{\mathrm{off}}$ is robustly set, ${t}_{\mathrm{rise}}$ and especially ${k}_{\mathrm{push}}$ are still undefined. We will discuss their impact on the model and on the explosion properties in detail in Sections 3.1.1–3.1.3. Ultimately, we will fix them using a calibration procure detailed in Section 3.4.

For the analysis of our results we determine several key quantities for each simulation. These quantities are obtained from a post-processing approach. We distinguish between the explosion properties, such as the explosion time, the mass cut,or the explosion energy, and the nucleosynthesis yields. The former are calculated from the hydrodynamics profiles. The latter are obtained from detailed nuclear network calculations for extrapolated trajectories.

2.5.1. Accretion Rates and Explosion Properties

For the accretion process, we distinguish between the accretion rate at the shock front, ${\dot{M}}_{\mathrm{shock}}={dM}({R}_{\mathrm{shock}})/{dt}$, and the accretion rate on the PNS, ${\dot{M}}_{\mathrm{PNS}}={dM}({R}_{\mathrm{PNS}})/{dt}$. In these expressions, M(R) is the baryonic mass enclosed in a radius R, ${R}_{\mathrm{shock}}$ is the shock radius, and ${R}_{\mathrm{PNS}}$ is the PNS radius that satisfies the condition $\rho ({R}_{\mathrm{PNS}})={10}^{11}\;{{\rm g}}\;{\mathrm{cm}}^{-3}$.

We consider the explosion time ${t}_{\mathrm{expl}}$ as the time when the shock reaches $500\;\mathrm{km}$, measured with respect to core bounce (see Ugliano et al. 2012). In all our models, the velocity of matter at the shock front has turned positive at that radius and the explosion has been irreversibly launched. There is no unique definition of ${t}_{\mathrm{expl}}$ in the literature and some other studies (see Janka & Mueller 1996; Handy et al. 2014) use different definitions, e.g., the time when the explosion energy increases above 10⁴⁸ erg. However, we do not expect that the different definitions give qualitatively different explosion times.

For the following discussion, we will use the total energy of the matter between a given mass shell m₀ up to the stellar surface:

$$\begin{eqnarray}&&{E}_{\mathrm{total}}({m}_{0},t)=-{\displaystyle \int }_{M}^{{m}_{0}}{e}_{\mathrm{total}}(m,t)\;{dm}.\end{eqnarray} \tag{ 12 }$$

M is the enclosed baryonic mass at the surface of the star and m₀ is a baryonic mass coordinate ($0\leqslant {m}_{0}\leqslant M$). ${e}_{\mathrm{total}}$ is the specific total energy, given by

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

i.e., the sum of the (relativistic) internal, kinetic, and gravitational specific energies. For all these quantities we make use of the general-relativistic expressions in the laboratory frame (Fischer et al. 2010). The integral in Equation (12) includes both the portion of the star evolved in the hydrodynamical simulation and the outer layers, which are considered as stationary profiles from the progenitor structure.

The explosion energy emerges from different physical contributions (see, for example, the appendix of Scheck et al. 2006 and the discussion in Ugliano et al. 2012). In our model, we are taking into account: (i) the total energy of the neutrino-heated matter that causes the shock revival; (ii) the nuclear energy released by the recombination of nucleons and alpha particles into heavy nuclei at the transition to non-NSE conditions; (iii) the total energy associated with the neutrino-driven wind developing after the explosion up to the end of the simulation; (iv) the energy released by the explosive nuclear burning in the shock-heated ejecta; and (v) the total (negative) energy of the outer stellar layers (also called the "overburden"). We are presently not taking into account the variation of the ejecta energy due to the appearance of late-time fallback. This is justified as long as the fallback represents only a small fraction of the total ejected mass.

To compute the explosion energy, we assume that the total energy of the ejecta with rest-masses subtracted eventually converts into kinetic energy of the expanding supernova remnant at $t\gg {t}_{\mathrm{expl}}$. The quantity ${e}_{\mathrm{total}}$ includes the rest mass contribution via ${e}_{\mathrm{int}}$, see Equation (2). Instead, if we want to calculate the explosion energy, we have to consider the thermal energy ${e}_{\mathrm{th}}$. Therefore, we define the specific explosion energy as

$$\begin{eqnarray}&&{e}_{\mathrm{expl}}={e}_{\mathrm{th}}+{e}_{\mathrm{kin}}+{e}_{\mathrm{grav}},\end{eqnarray} \tag{ 14 }$$

and the time- and mass-dependent explosion energy for the fixed mass domain between m₀ and M as

$$\begin{eqnarray}&&{H}_{\mathrm{expl}}({m}_{0},t)=-{\displaystyle \int }_{M}^{{m}_{0}}{e}_{\mathrm{expl}}(m,t)\;{dm}.\end{eqnarray} \tag{ 15 }$$

This can be interpreted as the total energy of this region in a non-relativistic EOS approach, where rest masses are not included.

The actual explosion energy (still time-dependent) is given by

$$\begin{eqnarray}&&{E}_{\mathrm{expl}}(t)={H}_{\mathrm{expl}}({m}_{\mathrm{cut}}(t),t),\end{eqnarray} \tag{ 16 }$$

i.e., for the matter above the mass cut ${m}_{\mathrm{cut}}(t)$.

To identify the mass cut, we consider the expression suggested by Bruenn in Fischer et al. (2010):

$$\begin{eqnarray}&&{m}_{\mathrm{cut}}(t)=m\left(\mathrm{max}({H}_{\mathrm{expl}}(m,t))\right),\end{eqnarray} \tag{ 17 }$$

where the maximum is evaluated outside the homologous core ($m\gtrsim 0.6$ ${{{\rm M}}}_{\odot }$), which has large positive values of the specific explosion energy ${e}_{\mathrm{expl}}$ once the PNS has formed due to the high compression. In the outer stellar envelope, before the passage of the shock wave, ${e}_{\mathrm{expl}}$ is dominated by the negative gravitational contribution. However, it is positive in the neutrino-heated region and in the shocked region above it. Hence, the above definition of ${m}_{\mathrm{cut}}$ locates essentially the transition from gravitationally unbound to bound layers. The final mass cut is obtained for $t={t}_{\mathrm{final}}$.

Our final simulation time ${t}_{\mathrm{final}}\gtrsim 4.6$ s is always much larger than the explosion time and, as we will show later, it allows ${E}_{\mathrm{expl}}(t)$ to saturate. Thus, we consider ${E}_{\mathrm{expl}}(t={t}_{\mathrm{final}})$ as the ultimate explosion energy of our models. In the following, if we use ${E}_{\mathrm{expl}}$ without the time as argument, we mean this final explosion energy.

2.5.2. Nucleosynthesis Yields

To predict the composition of the ejecta, we perform nucleosynthesis calculations using the full nuclear network Winnet (Winteler et al. 2012). We include isotopes up to ²¹¹Eu covering the neutron-deficient as well as the neutron-rich side of the valley of β-stability. The reaction rates are the same as in Winteler et al. (2012). They are based on experimentally known rates where available and predictions otherwise. The n, p, and alpha-captures are taken from Rauscher & Thielemann (2000), who used known nuclear masses where available and the Finite Range Droplet Model (Möller et al. 1995b) for unstable nuclei far from stability. The β-decay rates are from the nuclear database NuDat2.⁴

We divide the ejecta into different mass elements of 10⁻³ ${{{\rm M}}}_{\odot }$ each and follow the trajectory of each individual mass element. As we are mainly interested in the amounts of ⁵⁶Ni, ⁵⁷Ni, ⁵⁸Ni, and ⁴⁴Ti, we only consider the 340 innermost mass elements above the mass cut, corresponding to a total mass of 0.34 ${{{\rm M}}}_{\odot }$. The contribution of the outer mass elements to the production of those nuclei is negligible.

For $t\lt {t}_{\mathrm{final}}$, we use the temperature and density evolution from the hydrodynamical simulations as inputs for our network. For each mass element we start the nucleosynthesis post-processing when the temperature drops below 10 GK, using the NSE abundances (determined by the current electron fraction Y_e) as the initial composition. For mass elements that never reach 10 GK we start at the moment of bounce and use the abundances from the approximate α-network at this point as the initial composition. Note that for all tracers the further evolution of Y_e in the nucleosynthesis post-processing is determined inside the Winnet network.

At the end of the simulations, i.e., $t={t}_{\mathrm{final}}$, the temperature and density of the inner zones are still sufficiently high for nuclear reactions to occur ($T\approx 1$ GK and $\rho \approx 2.5\times {10}^{3}$ g cm⁻³). Therefore, we extrapolate the radius, density, and temperature up to ${t}_{\mathrm{end}}=100$ s using:

$$\begin{eqnarray}&&r(t)={r}_{\mathrm{final}}+{{tv}}_{\mathrm{final}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\rho (t)={\rho }_{\mathrm{final}}{\left(\displaystyle \frac{t}{{t}_{\mathrm{final}}}\right)}^{-3},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&T(t)=T[{s}_{\mathrm{final}},\rho (t),{Y}_{e}(t)],\end{eqnarray} \tag{ 20 }$$

where r is the radial position, v the radial velocity, ρ the density, T the temperature, s the entropy per baryon, and Y_e the electron fraction of the mass zone. The temperature is calculated at each timestep using the EOS of Timmes & Swesty (2000). The prescription in Equations (18)–(20) corresponds to a free expansion for the density and an adiabatic expansion for the temperature (see, for example, Korobkin et al. 2012).

3.1.1.  ${k}_{\mathrm{push}}$

The parameter with the most intuitive and strongest impact on the explosion is k_push. Its value directly affects the amount of extra heating which is provided by PUSH. As expected, larger values of k_push (assuming all other parameters to be fixed) result in the explosion being more energetic and occurring earlier. In addition, a faster explosion implies a lower remnant mass, as there is less time for the accretion to add mass to the forming PNS.

Beyond these general trends with k_push, the detailed behavior depends also on the compactness of the progenitor. For all 16 progenitor models in the 18–21 ${{{\rm M}}}_{\odot }$ ZAMS mass range, we have explored several PUSH models, varying ${k}_{\mathrm{push}}$ between 0.0 and 4.0 in increments of 0.5 but fixing ${t}_{\mathrm{on}}=80$ ms and ${t}_{\mathrm{rise}}=150$ ms. For ${k}_{\mathrm{push}}\leqslant 1$, none of the models explode and for ${k}_{\mathrm{push}}=1.5$ only the lowest compactness models explode. Figure 4 shows the explosion energy, the explosion time, and the (baryonic) remnant mass as function of the progenitor compactness for ${k}_{\mathrm{push}}=1.5,2.0,3.0,4.0$. A distinct behavior between LC and HC models is seen. The LC models (${\xi }_{1.75}\lt 0.4$) result in slightly weaker and faster explosions, with less variability in the explosion energy and in the explosion time for different values of k_push. Even for relatively large values of k_push, the explosion energies remain below 1 Bethe (1 Bethe, abbreviated as 1 B, is equivalent to 10⁵¹ erg). On the other hand, the HC models (${\xi }_{1.75}\gt 0.45$) explode stronger and later, with a larger variation in the explosion properties. In this case, for high enough values of k_push ($\gtrsim 3.0$), explosion energies of $\gtrsim 1$ Bethe can be obtained.

Zoom In Zoom Out Reset image size

The HC models also lead to a larger variability of the remnant masses, even though this effect is less pronounced than for the explosion time or energy. For the values of k_push used here, we obtain (baryonic) remnant masses of approximately 1.4–1.9 ${{{\rm M}}}_{\odot }$. The differences of LC and HC models will be investigated further in Section 3.3.

There are three models with $0.37\lesssim {\xi }_{1.75}\lesssim 0.50$ (corresponding to ZAMS masses of 18.0 (HC), 18.2 (LC), and 19.4 ${{{\rm M}}}_{\odot }$ (HC)) which do not follow the general trend. In particular, we find the threshold value of ${k}_{\mathrm{push}}$ for successful explosions to be higher for these models. A common feature of these three models is that they have the lowest Fe-core mass of all the models in our sample and the highest central densities at the onset of collapse.

The choice of t_rise does not affect the observed trends with k_push: similar behaviors are also seen for $50\;\mathrm{ms}\lesssim {t}_{\mathrm{rise}}\lesssim 250\;\mathrm{ms}$.

3.1.2.  ${t}_{\mathrm{on}}$

3.1.3.  ${k}_{\mathrm{push}}$ and ${t}_{\mathrm{rise}}$

In Sections 3.1.1 and 3.1.2, we have investigated the dependency of the model on the single parameters k_push and ${t}_{\mathrm{on}}$. Now, we explore the role of t_rise in combination with k_push. For this, we approximately fix the explosion energy to the canonical value of ∼1 B for the HC models (corresponding, for example, to the previously examined models with ${k}_{\mathrm{push}}=3.0$ and ${t}_{\mathrm{rise}}=150\;\mathrm{ms}$), and investigate which other combinations of k_push and t_rise result in the desired explosion energy. We restrict our explorations to a sub-set of progenitor models (18.0 ${{{\rm M}}}_{\odot }$, 18.6 ${{{\rm M}}}_{\odot }$, 19.2 ${{{\rm M}}}_{\odot }$, 19.4 ${{{\rm M}}}_{\odot }$, 19.8 ${{{\rm M}}}_{\odot }$, 20.0 ${{{\rm M}}}_{\odot }$, 20.2 ${{{\rm M}}}_{\odot }$, and 20.6 ${{{\rm M}}}_{\odot }$) that spans the ${\xi }_{1.75}$-range of all 16 progenitors. Figure 5 summarizes the explosion energies, explosion times, and remnant masses for various combinations of k_push and t_rise for progenitors of different compactness. The required constraint can be obtained by several combinations of parameters, which lie on a curve in the k_push–t_rise plane. As a general result, a longer t_rise requires a larger k_push to obtain the same explosion energy. This can be understood from the different roles of the two parameters: while k_push sets the maximum efficiency at which PUSH deposits energy from the reservoir represented by the ${\nu }_{\mu ,\tau }$ luminosity, t_rise sets the time scale over which the mechanism reaches this maximum. Together, they control the slope of ${\mathcal{G}}(t)$ in the rising phase (see Figure 3). A model with a longer rise time reaches its maximum efficiency later, at which time the luminosities have already decreased and a part of the absorbed energy has been advected on the PNS or re-emitted in the form of neutrinos. To compensate for these effects, a larger k_push is required for a longer t_rise. This is seen in Figure 6, where we plot the cumulative neutrino contribution $({E}_{\mathrm{push}}+{E}_{\mathrm{idsa}})$ and its time derivative for four runs of the 18.0 ${{{\rm M}}}_{\odot }$ progenitor model, but with different combinations of t_rise and k_push. Runs with larger parameter values require PUSH to deposit more energy (see $({E}_{\mathrm{push}}+{E}_{\mathrm{idsa}})$ at $t\approx {t}_{\mathrm{expl}}$), and the corresponding deposition rates are shifted toward later times. Moreover, for increasing values of ${t}_{\mathrm{rise}}$, the explosion time ${t}_{\mathrm{expl}}$ becomes larger, but the interval between $({t}_{\mathrm{on}}+{t}_{\mathrm{rise}})$ and ${t}_{\mathrm{expl}}$ decreases. Despite the significant variation of ${k}_{\mathrm{push}}$ between different runs, the peak values of $d({E}_{\mathrm{push}}+{E}_{\mathrm{idsa}})/{dt}$ at the onset of the shock revival that preceeds the explosion are very similar in all cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.4.  ${t}_{\mathrm{off}}$

In the following, we discuss the contributions to and the sources of the explosion energy, i.e., we investigate how the explosion energy is generated. This is done in several steps: first, we have a closer look at the neutrino energy deposition. Then we show how it relates to the increase of the total energy of the ejected layers, and finally how this increase of the total energy transforms into the explosion energy. For this analysis, we have chosen the 19.2 and 20.0 ${{{\rm M}}}_{\odot }$ ZAMS mass progenitor models as representatives of the HC and LC samples, respectively. We consider their exploding models obtained with ${t}_{\mathrm{on}}=80$ ms, ${t}_{\mathrm{rise}}=150$ ms, and ${k}_{\mathrm{push}}=3.0$. A summary of the explosion properties can be found in Table 3.

Table 3.  Explosion Properties for Two Reference Runs

Quantity		HC	LC
ZAMS	(${{{\rm M}}}_{\odot }$)	19.2	20.0
${\xi }_{1.75}$	(−)	0.637	0.283
${t}_{\mathrm{on}}$	(ms)	80
${t}_{\mathrm{rise}}$	(ms)	150
${k}_{\mathrm{push}}$	(−)	3.0
${t}_{\mathrm{expl}}$	(ms)	307	206
${M}_{\mathrm{remn}}$	(${{{\rm M}}}_{\odot }$)	1.713	1.469
${E}_{\mathrm{expl}}$ (${t}_{\mathrm{final}}$)	(B)	1.36	0.57
${E}_{\mathrm{push}}$ (${t}_{\mathrm{off}}+{t}_{\mathrm{rise}}$)	(B)	3.51	1.08
${E}_{\mathrm{idsa}}$ (${t}_{\mathrm{off}}+{t}_{\mathrm{rise}}$)	(B)	2.76	1.01
${E}_{\mathrm{idsa}}$ (${t}_{\mathrm{final}}$)	(B)	4.10	2.11

Note. These two runs are used to compare the HC and LC samples.

Download table as:  ASCIITypeset image

The table shows that for both models neutrinos are required to deposit a net cumulative energy $({E}_{\mathrm{push}}+{E}_{\mathrm{idsa}})$ much larger than ${E}_{\mathrm{expl}}$ to revive the shock and to lead to an explosion that matches the expected energetics. For the two reference runs, when the PUSH contribution is switched off ($t={t}_{\mathrm{off}}+{t}_{\mathrm{rise}}$), the cumulative deposited energy is ∼4 times larger than ${E}_{\mathrm{expl}}$. This can also be inferred from Figure 6 for other runs. That ratio increases further up to ∼5.5 at $t={t}_{\mathrm{final}}$, due to the neutrino energy deposition happening at the surface of the PNS which generates the ν-driven wind. According to Equations (9) and (10), ${E}_{\mathrm{push}}$ and ${E}_{\mathrm{idsa}}$ are the total energies which are deposited in the (time-dependent) gain region. This neutrino energy deposition increases the internal energy of the matter flowing in that region. However, since the advection timescale is much shorter than the explosion timescale, a large fraction of this energy is advected onto the PNS surface by the accreting mass before the explosion sets in, and hence does not contribute to the explosion energy. Only the energy deposited by neutrinos in the region above the final mass cut will eventually contribute to the explosion energy.

To identify this relevant neutrino contribution, in Figure 7 we show the time evolution of the integrated net neutrino energy deposition ${E}_{\nu }({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$ within the domain above the fixed mass ${m}_{\mathrm{cut}}^{\mathrm{fin}}={m}_{\mathrm{cut}}({t}_{\mathrm{final}})$. We choose ${m}_{\mathrm{cut}}^{\mathrm{fin}}$ to include all the relevant energy contributions to the explosion energy, up to the end of our simulations. Despite the significant differences in magnitudes, the two models show overall similar evolutions. If we compare ${E}_{\nu }({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$ at late times with $({E}_{\mathrm{push}}({t}_{\mathrm{off}}+{t}_{\mathrm{rise}})+{E}_{\mathrm{idsa}}({t}_{\mathrm{final}}))$ from Table 3, we see that it is significantly smaller. About two thirds of the energy originally deposited in the gain region are advected onto the PNS and hence do not contribute to the explosion energy.

Zoom In Zoom Out Reset image size

In addition to the neutrino energy deposition, in Figure 7 we also show the variation of the total energy for the domain above ${m}_{\mathrm{cut}}^{\mathrm{fin}}$, i.e., ${{\rm \Delta }}{E}_{\mathrm{total}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$ = ${E}_{\mathrm{total}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$ − ${E}_{\mathrm{total}}({m}_{\mathrm{cut}}^{\mathrm{fin}},{t}_{\mathrm{initial}})$, where ${t}_{\mathrm{initial}}$ is the time when we start our simulation from the stage of the progenitor star. The variation of the total energy can be separated into the net neutrino contribution and the mechanical work at the inner boundary, ${{\rm \Delta }}{E}_{\mathrm{total}}={E}_{\nu }+{E}_{\mathrm{mech}}$. We note that in our general relativistic approach the variation of the gravitational mass due to the intense neutrino emission from the PNS is consistently taken into accout. It is visible in Figure 7, that the net deposition by neutrinos makes up the largest part of the change of the total energy. The transfer of mechanical energy ${E}_{\mathrm{mech}}$ is negative because of the expansion work performed by the inner boundary during the collapse and the PNS shrinking. However, it is significantly smaller in magnitude than ${E}_{\nu }$.

Next, we investigate the connection between the variation of the total energy and the explosion energy. In Figure 7, we show the variation of the explosion energy above the fixed mass ${m}_{\mathrm{cut}}^{\mathrm{fin}}$, i.e., ${{\rm \Delta }}{H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)={H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)-{H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},{t}_{\mathrm{initial}})$, together with the relative variation of the time-dependent explosion energy, ${{\rm \Delta }}{E}_{\mathrm{expl}}(t)={E}_{\mathrm{expl}}(t)-{H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},{t}_{\mathrm{initial}})$. It is obvious from Equations (2), (12), and (15) that the difference between ${{\rm \Delta }}{H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$ and ${{\rm \Delta }}{E}_{\mathrm{total}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$ is given by the variation of the integrated rest mass energy, ${{\rm \Delta }}{H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)={{\rm \Delta }}{E}_{\mathrm{total}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)-{{\rm \Delta }}{E}_{\mathrm{mass}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$. In Figure 7, $-{{\rm \Delta }}{E}_{\mathrm{mass}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$ can thus be identified as the difference between the long-thin and the short-thick dashed lines. We find that the overall rest mass contribution to the final explosion energy is positive, but much smaller than the neutrino contribution. Figure 7 also makes evident the conceptual difference between ${H}_{\mathrm{expl}}$ and ${E}_{\mathrm{expl}}$, and, at the same time, shows that ${H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)\to {E}_{\mathrm{expl}}(t)$ for $t\to {t}_{\mathrm{final}}$, since we have chosen ${m}_{\mathrm{cut}}^{\mathrm{fin}}={m}_{\mathrm{cut}}({t}_{\mathrm{final}})$. It also reveals that the explosion energy ${E}_{\mathrm{expl}}$ has practically saturated for $t\gtrsim 1\;{{\rm s}}$, while ${E}_{\nu }$ (and, consequently, ${{\rm \Delta }}{E}_{\mathrm{tot}}$ and ${{\rm \Delta }}{H}_{\mathrm{expl}}$) increases up to ${t}_{\mathrm{final}}$, when ${m}_{\mathrm{cut}}^{\mathrm{fin}}$ is finally ejected. However, this energy provided by neutrinos is mostly spent to unbind matter from the PNS surface. Thus, the late ν-driven wind, which occurs for several seconds after 1 s, still increases ${E}_{\mathrm{expl}}$, but at a relative small, decreasing rate.

To summarize, the variation of the explosion energy above ${m}_{\mathrm{cut}}^{\mathrm{fin}}$ can be expressed as

$$\begin{eqnarray}{{\rm \Delta }}{H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t) & = & {{\rm \Delta }}{E}_{\mathrm{total}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)-{{\rm \Delta }}{E}_{\mathrm{mass}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)\\ & = & {E}_{\nu }({m}_{\mathrm{cut}}^{\mathrm{fin}},t)+{E}_{\mathrm{mech}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)\\ & & -{{\rm \Delta }}{E}_{\mathrm{mass}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t).\end{eqnarray} \tag{ 21 }$$

The quantity $-{{\rm \Delta }}{E}_{\mathrm{mass}}$ is positive, but significantly smaller than ${E}_{\nu }({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$. ${E}_{\mathrm{mech}}$ is negative and also smaller than ${E}_{\nu }({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$. Therefore, we conclude that in our models the explosion energy is mostly generated by the energy deposition of neutrinos in the eventually ejected layers, especially within the first second after bounce.

To give further insight, in Figure 8 we show the time evolution of all energies which contribute to the explosion energy together with the explosion energy itself, for both the HC (left panel) and the LC model (right panel). We present ${E}_{\mathrm{int}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$, $-{E}_{\mathrm{mass}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$, ${E}_{\mathrm{grav}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$, and ${E}_{\mathrm{kin}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)$, which together give a complete decomposition of the explosion energy, i.e.,

$$\begin{eqnarray}{H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t) & = & {E}_{\mathrm{kin}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)+{E}_{\mathrm{grav}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)\\ & & +{E}_{\mathrm{int}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)-{E}_{\mathrm{mass}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t).\end{eqnarray} \tag{ 22 }$$

Compared to Figure 7 we are now not dealing with variations anymore, but with absolute values. Gravitational energy initially dominates (${H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t)\lt 0$), meaning that the portion of the star above ${m}_{\mathrm{cut}}^{\mathrm{fin}}$ is still gravitationally bound. The HC model is initially more bound than the LC model (for example, ${H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t=0.1\;{{\rm s}})\approx -0.54\;{{\rm B}}$, versus ${H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},t=0.1\;{{\rm s}})\approx -0.40\;{{\rm B}}$, respectively). Before providing positive explosion energy, neutrinos have to compensate for this initial negative binding energy as well as for the negative ${E}_{\mathrm{mech}}$. This can be seen explicitly by expressing Equation (21) as:

$$\begin{eqnarray}{H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},{t}_{\mathrm{final}}) & \sim & {H}_{\mathrm{expl}}({m}_{\mathrm{cut}}^{\mathrm{fin}},{t}_{\mathrm{initial}})+{E}_{\nu }({m}_{\mathrm{cut}}^{\mathrm{fin}},{t}_{\mathrm{final}})\\ & & +{E}_{\mathrm{mech}}({t}_{\mathrm{final}}),\end{eqnarray} \tag{ 23 }$$

where we have neglected ${{\rm \Delta }}{E}_{\mathrm{mass}}$.

Zoom In Zoom Out Reset image size

In the following, we discuss the evolution of the relevant energies and, in particular, of the rest mass energy (see Section 2.2 for the description of the (non-) NSE EOS and of the related definitions of the internal, thermal and rest mass energies). The innermost part of the ejecta (corresponding to ∼0.15 ${{{\rm M}}}_{\odot }$ and ∼0.07 ${{{\rm M}}}_{\odot }$ above ${m}_{\mathrm{cut}}^{\mathrm{fin}}$ for the 19.2 ${{{\rm M}}}_{\odot }$ and 20.0 ${{{\rm M}}}_{\odot }$ model, respectively) is initially composed of intermediate mass nuclei (mainly silicon and magnesium). In the first part of the evolution, during the gravitational collapse, no significant changes of ${E}_{\mathrm{int}}$ and ${E}_{\mathrm{mass}}$ are observed in Figure 8. However, when this matter enters the shock, it is quickly photodissociated into neutrons, protons, and alpha particles. This process increases the rest mass energy, as is visible in Figure 8 between roughly 200 and 300 ms for the HC model and between 100 and 200 ms for the LC model. At the same time, the release of gravitational energy of the still infalling matter and the dissipation of kinetic energy happening at the shock, together with the large and intense neutrino absorption on free nucleons, increase ${E}_{\mathrm{int}}$. Later, once neutrino heating has halted the collapse and started the explosion, the expanding shock decreases its temperature and free neutrons and protons inside it recombine first into alpha particles and then into iron group nuclei. At the same time, fresh infalling layers are heated by the shock to temperatures above 0.44 MeV, and silicon and magnesium are converted into heavier nuclei and alpha particles under NSE conditions, leading to an alpha-rich freeze-out from NSE. The production of alpha particles, which are less bound than the heavy nuclei initially present in the same layers, limits the amount of rest mass energy finally released. Thus, these recombination and burning processes liberate in a few hundred milliseconds after ${t}_{\mathrm{expl}}$ an amount of rest mass energy larger but comparable to the energy spent by the shock to photodissociate the infalling matter during shock revival and early expansion. We have checked in post-processing that the full nucleosynthesis network Winnet confirms these results.

The distributions of the explosion energy and explosion time obtained with PUSH, as well as their variations in response to changes of the model parameters, suggest a possible distinction between HC and LC progenitors. In the following, we investigate how basic properties of the models (e.g., the accretion history or the neutrino luminosities), ultimately connected with the compactness, relate to differences in the explosion process and properties. For a similar discussion in self-consistent 1D and 2D supernova simulations, see Suwa et al. (2014). Again, we choose the 19.2 and 20.0 ZAMS mass progenitor runs with ${t}_{\mathrm{rise}}=150\;\mathrm{ms}$ and ${k}_{\mathrm{push}}=3.0$, as representatives of the HC and LC samples, respectively.

In Figure 9, we show the temporal evolution of several quantities of interest for both the 19.2 ${{{\rm M}}}_{\odot }$ and 20.0 ${{{\rm M}}}_{\odot }$ models, with and without PUSH. The evolution before ${t}_{\mathrm{on}}$ follows the well known early shock dynamics in CCSNe (see, for example, Burrows & Goshy 1993). In both models, a few tens of milliseconds after core bounce, the expanding shock turns into an accretion front, and the mantle between the PNS surface and the shock reaches a quasi-stationary state. In this accretion phase, ${\dot{M}}_{\mathrm{shock}}$ and ${\dot{M}}_{\mathrm{PNS}}$ are firmly related. However, the two different density profiles already affect the evolution of the shock. Since ${\rho }_{19.2}/{\rho }_{20.0}\gtrsim 1.2$ outside the shock and up to a radius of $2\times {10}^{8}$ cm (while the infalling velocities of the unshocked matter are initially almost identical), ${\dot{M}}_{\mathrm{shock}}$ (and in turn also ${\dot{M}}_{\mathrm{PNS}}$) starts to differ between the two models around ${t}_{\mathrm{pb}}\approx 30$ ms.

Zoom In Zoom Out Reset image size

The difference in the accretion rates has a series of immediate consequences. For the HC case, (i) neutrino luminosities are larger (Figure 9(c)); (ii) the shock is subject to a larger ram pressure (i.e., a larger momentum transfer provided by the collectively infalling mass flowing through the shock), and, as visible in the case without PUSH, shock stalling happens earlier and at a smaller radius (Figure 9(b)); (iii) the PNS mass grows faster. Since the mass of the PNS at bounce is almost identical for the two models (${M}_{\mathrm{PNS}}\approx 0.63$ ${{{\rm M}}}_{\odot }$), the stronger gravitational potential implied by (iii) increases the differences in the accretion rates even further by augmenting the ratio of the radial velocities inside the gain region (larger by 12–15% at $t\approx {t}_{\mathrm{on}}$ for the 19.2 ${{{\rm M}}}_{\odot }$ case).

For $t\gt {t}_{\mathrm{on}}$, the differences between the two runs amplify as a result of the PUSH action. In the LC case, due to the lower accretion rate, a relatively small energy deposition by PUSH in the gain region (smaller than or comparable to the energy deposition by ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ from IDSA, as visible in Figure 10) is able to revive the shock expansion a few milliseconds after ${t}_{\mathrm{on}}$. Later, the increasing ${{dE}}_{\mathrm{push}}/{dt}$ triggers an explosion in a few tens of milliseconds, even before ${\mathcal{G}}(t)$ reaches its maximum (Figure 9(b)). In the HC case, the energy deposition by neutrinos is more intense from the beginning due to the larger neutrino luminosities and harder neutrino spectra (Figures 9(c) and (d)) and due to the higher density inside the gain region. However, because of the larger accretion rate, the extra contribution provided by PUSH is initially only able to prevent the fast shock contraction observed in the model without PUSH. During this shock stalling phase, the accretion rate and the luminosity decrease, but only marginally and very similarly to the non-exploding case. When PUSH reaches its maximum energy deposition rate ($t\approx {t}_{\mathrm{on}}+{t}_{\mathrm{rise}}$), the shock revives and the explosion sets in (Figure 9(b)).

Zoom In Zoom Out Reset image size

In Figure 11, we plot the ratio of the ram pressure just above the shock front (${P}_{\mathrm{ram}}({R}_{\mathrm{shock}}^{+})=\rho {v}^{2}$ calculated at ${R}_{\mathrm{shock}}^{+}={R}_{\mathrm{shock}}+1\;\mathrm{km}$) to the thermal pressure just inside it (${P}_{\mathrm{th}}({R}_{\mathrm{shock}}^{-})$ where ${R}_{\mathrm{shock}}^{-}={R}_{\mathrm{shock}}-1\;\mathrm{km}$). In the non-exploding runs (i.e., without PUSH), both these pressures decrease with time, but their ratio stays always well above unity. On the other hand, in runs with PUSH, the more efficient energy deposition by neutrinos reduces the decrease of the thermal pressure inside the shock. The corresponding drop in the pressure ratio below unity determines the onset of the explosion.

Zoom In Zoom Out Reset image size

In both runs, once the explosion has been launched, the density in the gain region decreases and the PUSH energy deposition rate reduces accordingly. The conversion from an accreting to an expanding shock front decouples ${\dot{M}}_{\mathrm{shock}}$ from ${\dot{M}}_{\mathrm{PNS}}$. The latter drops steeply, together with the accretion neutrino luminosities (Figures 9(a) and (c)), while ${\dot{M}}_{\mathrm{shock}}$ decreases first but then stabilizes around an almost constant (slightly decreasing) value. In the case where the shock expansion velocity is much larger than the infalling matter velocity at ${R}_{\mathrm{shock}}$, ${\dot{M}}_{\mathrm{shock}}$ can be re-expressed as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{shock}}\approx 4\pi {R}_{\mathrm{shock}}^{2}\;\rho ({R}_{\mathrm{shock}})\;{v}_{\mathrm{shock}},\end{eqnarray} \tag{ 24 }$$

where ${v}_{\mathrm{shock}}={{dR}}_{\mathrm{shock}}/dt\propto {R}_{\mathrm{shock}}^{\delta }$. For $R\gt {R}_{\mathrm{shock}}$ we have in good approximation $\rho (R)\propto {R}^{-2}$, and thus

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{shock}}\propto {R}_{\mathrm{shock}}^{\delta }.\end{eqnarray} \tag{ 25 }$$

The stationary value of ${\dot{M}}_{\mathrm{shock}}$ implies that $\delta \approx 0$. Thus, after an initial exponential expansion, the shock velocity is almost constant during the first second after the explosion.

Despite the larger difficulties to trigger an explosion, the HC model explodes more energetically than the LC model. According to the analysis performed in Section 3.2, the difference in the explosion energy between the HC and the LC model depends ultimately on the different amount of energy deposited by neutrinos. Since the HC model requires a larger energy deposition to overcome the ram pressure and the gravitational potential, the total energy of the corresponding ejecta (and in turn the explosion energy) will be more substantially increased. In addition, after the explosion has been triggered, the larger neutrino luminosities and densities that characterize the HC model inject more energy in the expanding shock compared with the LC model.

The ultimate goal of core-collapse supernova simulations is to reproduce the properties observed in real supernovae. So far we have only focused on the dependence of dynamical features of the explosion (e.g., the explosion energy) on the parameter choices in the PUSH method. However, the ejected mass of radioactive nuclides (such as ⁵⁶Ni) is an equally important property of the supernova explosion. Here, we describe how we calibrate the PUSH method by reproducing the explosion energy and mass of Ni ejecta of SN 1987A for a progenitor within the expected mass range for this supernova.

3.4.1. Observational Constraints from SN 1987A

The analysis and the modeling of the observational properties of SN 1987A just after the luminosity peak have been the topics of a long series of works (e.g., Woosley 1988; Arnett et al. 1989; Shigeyama & Nomoto 1990; Kozma & Fransson 1998a, 1998b; Blinnikov et al. 2000; Fransson & Kozma 2002; Utrobin & Chugai 2005; Seitenzahl et al. 2014, and references therein). They provide observational estimates for the explosion energy, the progenitor mass, and the ejected masses of ⁵⁶Ni, ⁵⁷Ni, ⁵⁸Ni, and ⁴⁴Ti, all of which carry rather large uncertainties. In Table 4, the values used for the calibration of the PUSH method are summarized.

Table 4.  Observational Properties of SN 1987A

${E}_{\mathrm{expl}}$	$(1.1\pm 0.3)\times {10}^{51}$ erg
${m}_{\mathrm{prog}}$	18–21 ${{{\rm M}}}_{\odot }$
$m{(}^{56}\mathrm{Ni})$	$(0.071\pm 0.003)$ ${{{\rm M}}}_{\odot }$
$m{(}^{57}\mathrm{Ni})$	$(0.0041\pm 0.0018)$ ${{{\rm M}}}_{\odot }$
$m{(}^{58}\mathrm{Ni})$	0.006 ${{{\rm M}}}_{\odot }$
$m{(}^{44}\mathrm{Ti})$	$(0.55\pm 0.17)\times {10}^{-4}$ ${{{\rm M}}}_{\odot }$

Note. The nucleosynthesis yields are taken from Seitenzahl et al. (2014) except for ⁵⁸Ni which is taken from Fransson & Kozma (2002). No error estimates were given for ⁵⁸Ni. The explosion energy is adapted from Blinnikov et al. (2000). For the progenitor range we chose typical values found in the literature, see e.g., Shigeyama & Nomoto (1990) and Woosley (1988).

Download table as:  ASCIITypeset image

The ZAMS progenitor mass is assumed to be between 18 ${{{\rm M}}}_{\odot }$ and 21 ${{{\rm M}}}_{\odot }$, corresponding to typical values reported in the literature for the SN 1987A progenitor, see, e.g., Woosley (1988) and Shigeyama & Nomoto (1990). For the explosion energy we consider the estimate reported by Blinnikov et al. (2000), ${E}_{\mathrm{expl}}=(1.1\pm 0.3)\times {10}^{51}$ erg (for a detailed list of explosion energy estimates for SN 1987 A, see for example Table 1 in Handy et al. 2014). This value was obtained assuming ∼14.7 ${{{\rm M}}}_{\odot }$ of ejecta and an hydrogen-rich envelope of ∼10.3 ${{{\rm M}}}_{\odot }$. The uncertainties in the progenitor properties and in the SN distance were taken into account in the error bar. The employed values of the total ejecta and of the hydrogen-rich envelope are compatible (within a 15% tolerance) with a significant fraction of our progenitor candidates, especially for ${M}_{\mathrm{ZAMS}}\lt 19.6$ ${{{\rm M}}}_{\odot }$ (see Table 2, where the total ejecta can be estimated subtracting 1.6 ${{{\rm M}}}_{\odot }$ from the mass of the star at the onset of the collapse). Explosion models with larger ejected mass (i.e., less compatible with our candidate sample) tend to have larger explosion energies (see, for example, Utrobin & Chugai 2005). Finally, we consider the element abundances for ${}^{\mathrm{56,57}}\mathrm{Ni}$ and ⁴⁴Ti provided by Seitenzahl et al. (2014), which were obtained from a least squares fit of the decay chains to the bolometric lightcurve. For ⁵⁸Ni we use the value provided by Fransson & Kozma (2002).

3.4.2. Fitting Procedure

We calibrate the PUSH method by finding a combination of progenitor mass, k_push, and t_rise which provides the best fit to the all observational quantities of SN 1987A mentioned above. The weight given to each quantity is related to the uncertainty. For example, due to the large uncertainty in the ⁴⁴Ti mass, this does not provide a strong constraint on selecting the best fit.

Figure 12 shows the explosion energy and ejected mass of ⁵⁶Ni, ⁵⁷Ni, ⁵⁸Ni, and ⁴⁴Ti for different cases of k_push and t_rise and for four select HC progenitors used to calibrate the PUSH method. We do not consider the LC progenitors because of their generally lower explosion energies, see Figure 4. The different cases of k_push and t_rise span a wide range of explosion energies around 1 Bethe. For all parameter combinations shown, at least one progenitor in the 18–21 ${{{\rm M}}}_{\odot }$ range fulfills the requirement of an explosion energy between 0.8 Bethe and 1.4 Bethe.

Zoom In Zoom Out Reset image size

There is a roughly linear correlation between the explosion energy and the synthesized ⁵⁶Ni mass. However, this correlation is not directly compatible with the observations, as the ejected ⁵⁶Ni is systematically larger than expected (up to a factor of ∼2 for models with an explosion energy around 1 Bethe). There is a weak trend that models with higher ${t}_{\mathrm{rise}}$ tend to give lower nickel masses for given explosion energy. Among the parameter combinations that produce robustly high explosion energies (i.e., ${k}_{\mathrm{push}}\geqslant 3$), ${k}_{\mathrm{push}}=3.5$ with the high value of ${t}_{\mathrm{rise}}$ of 200 ms gives the lowest ⁵⁶Ni mass for similar explosion energies, but still much too high.

Our simulations can be reconciled with the observations by taking into account fallback from the initially unbound matter. Since we do not model the explosion long enough to see the development of the reverse shock and the appearance of the related fallback when the shock reaches the hydrogen-rich envelope, we have to impose it, removing some matter from the innermost ejecta.⁵ With a value of ∼0.1 ${{{\rm M}}}_{\odot }$ we can match both the expected explosion energy and ⁵⁶Ni ejecta mass, see Figure 13. In this way we have fixed the final mass cut by observations. However, we point out that we are able to identify the amount of late-time fallback only because we also have the dynamical mass cut from our hydrodynamical simulations. This is not possible in other methods such as pistons or thermal bombs. Our value of ∼0.1 ${{{\rm M}}}_{\odot }$ of fallback in SN 1987A will be further discussed and compared with other works in Section 4.3.

Zoom In Zoom Out Reset image size

The observed yield of ⁵⁶Ni provides a strong constraint on which parameter combination would fit the data. From the observed yields of ⁵⁷Ni and ⁵⁸Ni, only the 18.0 and 19.4 progenitors remain viable candidates. Without fallback our predicted ⁴⁴Ti yields are compatible with the observed yields (see Figure 12). However, if we include fallback (which is needed to explain the observed Ni yields), ⁴⁴Ti becomes underproduced compared to the oberved value. Since this behavior is true for all out models, we exclude the constraint given by ⁴⁴Ti from our calibration procedure. From the considered parameter combinations, we obtained the best fit to SN 1987A for the 18.0 ${{{\rm M}}}_{\odot }$ progenitor model with ${k}_{\mathrm{push}}=3.5$, ${t}_{\mathrm{rise}}=200$ ms, and a fallback of 0.1 ${{{\rm M}}}_{\odot }$. These parameters are summarized in Table 5. In Figure 14, we show the temporal evolution of the accretion rates, of the relevant radii, and of the neutrino luminosities and mean energies for our best fit model. For comparison purposes, we present also the results obtained for the same model without PUSH. Note that in this non-exploding case the ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ luminosities stay almost constant over several ∼100 ms after core bounce, despite the decreasing accretion rate. This is due to the relatively slow variation of ${\dot{M}}_{\mathrm{PNS}}$ (for example, compared with the variation obtained in the 19.2 ${{{\rm M}}}_{\odot }$ model, Figure 9) and due to the simultaneous increase of the PNS gravitational potential, proportional to ${M}_{\mathrm{PNS}}/{R}_{\mathrm{PNS}}$ (see, for example, Fischer et al. 2009). A summary of the most important results of the simulations using this parameter set for the different progenitors in the 18–21 ${{{\rm M}}}_{\odot }$ window is given in Table 6. For the remnant mass and for the ⁵⁶Ni yields of our best-fit model, we provide both the values obtained with and without assuming a fallback of 0.1 ${{{\rm M}}}_{\odot }$.

Zoom In Zoom Out Reset image size

Table 5.  Parameter Values for Best Fit to SN 1987A

kpush	trise	${t}_{\mathrm{on}}$	${t}_{\mathrm{off}}$
(−)	(ms)	(ms)	(s)
3.5	200	80	1

Note. We identified the 18.0 ${{{\rm M}}}_{\odot }$ model as the progenitor which fits best, whereas we had to impose a late-time fallback of 0.1 ${{{\rm M}}}_{\odot }$.

Download table as:  ASCIITypeset image

Table 6.  Summary of Simulations for ${k}_{\mathrm{push}}=3.5$ and ${t}_{\mathrm{rise}}=200$ ms

ZAMS	${E}_{\mathrm{expl}}$	${t}_{\mathrm{expl}}$	${M}_{\mathrm{remnant}}^{B}$	${M}_{\mathrm{remnant}}^{G}$	M(56Ni )
(${{{\rm M}}}_{\odot }$)	(Bethe)	(s)	(${{{\rm M}}}_{\odot }$)	(${{{\rm M}}}_{\odot }$)	(${{{\rm M}}}_{\odot }$)
18.0	1.092	0.304	1.563	1.416	0.158
18.2	0.808	0.249	1.509	1.371	0.110
18.4	1.358	0.318	1.728	1.549	0.144
18.6	0.702	0.239	1.529	1.388	0.090
18.8	0.721	0.236	1.522	1.382	0.093
19.0	1.366	0.317	1.716	1.54	0.161
19.2	1.356	0.318	1.724	1.546	0.152
19.4	1.15	0.326	1.608	1.452	0.158
19.6	0.371	0.230	1.584	1.433	0.04
19.8	0.661	0.225	1.523	1.383	0.088
20.0	0.613	0.222	1.474	1.342	0.085
20.2	0.379	0.224	1.554	1.408	0.039
20.4	0.743	0.263	1.674	1.506	0.094
20.6	1.005	0.277	1.781	1.592	0.141
20.8	0.959	0.277	1.764	1.578	0.135
21.0	1.457	0.316	1.733	1.554	0.198
18.0 (fb)	1.092	0.304	1.663	1.497	0.073

Note. For the model 18.0 (fb), which is our best fit to SN 1987A, we have included 0.1 ${{{\rm M}}}_{\odot }$ of fallback, determined from obervational constraints. See the text for more details.

Download table as:  ASCIITypeset image

Figures 12 and 13 show that the composition of the ejecta is highly dependent on the progenitor model, especially for the amount of ⁵⁷Ni and ⁵⁸Ni ejected. From the four HC progenitors shown, two (18.0 ${{{\rm M}}}_{\odot }$ and 19.4 ${{{\rm M}}}_{\odot }$) produce a fairly high amount of those isotopes, while the other two (19.2 ${{{\rm M}}}_{\odot }$ and 20.6 ${{{\rm M}}}_{\odot }$) do not reach the amount observed in SN 1987A. A thorough investigation of the composition profile of the ejecta reveals that ⁵⁷Ni and ⁵⁸Ni are mainly produced in the slightly neutron-rich layers (${Y}_{e}\lt 0.5$), where the alpha-rich freeze-out leads to nuclei only one or two neutron units away from the N = Z line. A comparison of the Y_e and composition profiles for the 18.0 ${{{\rm M}}}_{\odot }$ and the 20.6 ${{{\rm M}}}_{\odot }$ progenitors is shown in Figure 15. For the 18.0 ${{{\rm M}}}_{\odot }$ model, the cut-off mass is 1.56 ${{{\rm M}}}_{\odot }$ and a large part of the silicon shell is ejected. In this shell, the initial matter composition is slightly neutron-rich (due to a small contribution from ⁵⁶Fe) with ${Y}_{e}\simeq 0.498$ (dotted line in top left graph) and the conditions for the production of ⁵⁷Ni and ⁵⁸Ni are favorable. The increase in Y_e around 1.9 ${{{\rm M}}}_{\odot }$ marks the transition to the oxygen shell. The same transition for the 20.6 ${{{\rm M}}}_{\odot }$ model happens around 1.74 ${{{\rm M}}}_{\odot }$, i.e., inside the mass cut. Therefore, this model ejects less ⁵⁷Ni and ⁵⁸Ni (see also Thielemann et al. 1990). In all our models, ⁴⁴Ti is produced within the innermost 0.15 ${{{\rm M}}}_{\odot }$ of the ejecta (see Figure 15). Since we assume 0.1 ${{{\rm M}}}_{\odot }$ fallback onto the PNS, most of the synthesized ⁴⁴Ti is not ejected in our simulations.

Zoom In Zoom Out Reset image size

While post-processing the ejecta trajectories for nucleosynthesis, Y_e is evolved by the nuclear network independently of the hydrodynamical evolution. This leads to a discrepancy at later times between the electron fraction in the initial trajectory (${Y}_{e}^{\mathrm{hydro}}$) and in the network (${Y}_{e}^{\mathrm{nuc}}$). In order to estimate the possible error in our nucleosynthesis calculations arising from this discrepancy, we have performed reference calculations using ${Y}_{e}^{\mathrm{hydro}}(t={t}_{\mathrm{final}})$ instead of ${Y}_{e}^{\mathrm{hydro}}(T=10\ \mathrm{GK})$ as a starting value for the network (see Section 2.5.2). The results are shown in Figure 15 for two progenitors: 18.0 ${{{\rm M}}}_{\odot }$ and 20.6 ${{{\rm M}}}_{\odot }$. The label "standard" refers to the regular case which uses ${Y}_{e}^{\mathrm{hydro}}(T=10\ \mathrm{GK})$ as input. The calculation using ${Y}_{e}^{\mathrm{hydro}}(t={t}_{\mathrm{final}})$ as input is labeled "alternative" and is represented by the dashed lines. The point in time at which the Y_e profile is shown is indicated by the supplements "input" (before the first timestep) and "final" (at t = 100 s). The corresponding nuclear compositions of the ejecta, each at the final calculation time of 100 s, are shown in the bottom panels. For the alternative Y_e profile of the 18.0 ${{{\rm M}}}_{\odot }$ progenitor (top left) the minimum around 1.59 ${{{\rm M}}}_{\odot }$ disappears, leading to an increase in ⁵⁶Ni in this region at the expense of ⁵⁷Ni and ⁵⁸Ni (bottom left). For the 20.6 ${{{\rm M}}}_{\odot }$ progenitor the situation is similar, with only a very small region just above 1.8 ${{{\rm M}}}_{\odot }$ showing significant differences. In general, we observe that the uncertainties in Y_e in our calculations are only present up to 0.05 ${{{\rm M}}}_{\odot }$ above the mass cut. The resulting uncertainties in the composition of the ejecta are very small or even inexistent in the scenarios where we consider fallback.

The radioactive isotope ⁴⁴Ti can be detected in supernovae and supernova remnants. Several groups have used different techniques to estimate the ⁴⁴Ti yield (Chugai et al. 1997; Fransson & Kozma 2002; Jerkstrand et al. 2011; Larsson et al. 2011; Grebenev et al. 2012; Grefenstette et al. 2014; Seitenzahl et al. 2014). The inferred values span a broad range, $(0.5{{\rm -}}4)\times {10}^{-4}$ ${{{\rm M}}}_{\odot }$. Traditional supernova nucleosynthesis calculations (e.g., Woosley & Weaver 1995; Thielemann et al. 1996) typically predict too low ⁴⁴Ti yields. Only very few models predict high ⁴⁴Ti yields: Thielemann et al. (1990) report ⁴⁴Ti yields around 10⁻⁴ and above in the best fits of their artificial SN explosions to SN 1987A. Rauscher et al. (2002) argue that the yields of ⁵⁶Ni and ⁴⁴Ti are very sensitive to the "final mass cut" (as we have shown, too), which is often determined by fallback. Ejecta in a supernova may be subject to convective overturn. To account for this, we can assume homogeneous mixing in the inner layers up to the outer boundary of the silicon shell before cutting off the fallback material (see, for example, Umeda & Nomoto 2002 and references therein). For our best-fit model, the ejected ⁴⁴Ti mass increases to $2.70\;\times \;{10}^{-5}$ ${{{\rm M}}}_{\odot }$, if this prescription is applied. Comparing to the previous yield of $1.04\;\times \;{10}^{-5}$ ${{{\rm M}}}_{\odot }$, we observe that the effect of homogeneous mixing is considerable, but not sufficient to match the observational values. The ejected ${}^{56{{\rm -}}58}\mathrm{Ni}$ masses also show a slight increase. However, there are also uncertainties in the nuclear physics connected to the production and destruction of ⁴⁴Ti. The final amount of produced ⁴⁴Ti depends mainly on two reactions: ⁴⁰Ca($\alpha ,\gamma$)⁴⁴Ti and ⁴⁴Ti($\alpha ,p$)⁴⁷V. Recent measurements of the ⁴⁴Ti($\alpha ,p$)⁴⁷V reaction rate within the Gamow window concluded that it may be considerably smaller than previous theoretical predictions (Margerin et al. 2014). In this study, an upper limit cross-section is reported that is a factor of 2.2 smaller than the cross-section we have used in our calculations (at a confidence level of 68%). Using this smaller cross-section for the ⁴⁴Ti($\alpha ,p$)⁴⁷V reaction, our yield of ejected ⁴⁴Ti for our best-fit model (18.0 ${{{\rm M}}}_{\odot }$ progenitor, ${k}_{\mathrm{push}}=3.5$, ${t}_{\mathrm{rise}}=200$ ms) rises to $1.49\;\times \;{10}^{-5}$ ${{{\rm M}}}_{\odot }$ with fallback and $5.65\;\times \;{10}^{-5}$ ${{{\rm M}}}_{\odot }$ without fallback. This corresponds to a relative increase of 43% with fallback and 48% without fallback. If we include both the new cross-section and homogeneous mixing, the amount of ⁴⁴Ti in the ejecta is $3.99\;\times \;{10}^{-5}$ ${{{\rm M}}}_{\odot }$ including fallback. This value, however, is still below the expected value derived from observations, but within the error box.

In the analysis of the nucleosynthesis yields above, we have used a mass resolution of 0.001 ${{{\rm M}}}_{\odot }$ for the tracers. This is too coarse to resolve the ejecta of the late neutrino-driven wind. Note that in our best-fit approach, where no mixing is assumed, none of the neutrino-driven wind is ejected because it is part of the fallback. Nevertheless, in the following we report briefly on the properties of the wind obtained by our detailed neutrino-transport scheme. For our best-fit model, the 18.0 ${{{\rm M}}}_{\odot }$ progenitor, at ${t}_{\mathrm{final}}$ we find an electron fraction around 0.32, entropies up to 80 k_B per baryon, and fast expansion velocities ($\sim {10}^{9}$ cm s⁻¹). Similar conditions are also found for the other progenitors. They are not sufficient for a full r-process (see, for example, Farouqi et al. 2010). On the other hand, we have found that the entropy is still increasing and the electron fraction still decreasing in the further evolution. The high asymmetries are only obtained if we include the nucleon mean-field interaction potentials in the neutrino charged-current rates (Martínez-Pinedo et al. 2012). However, they are much higher than found in other long-term simulations which also include these potentials (Martínez-Pinedo et al. 2012, 2014; Roberts et al. 2012). This could be related to the missing neutrino-electron scattering in our neutrino transport, which is an important source of thermalization and down-scattering, especially for the high energy electron anti-neutrinos at late times, see Fischer et al. (2012). More detailed comparisons are required to identify the origin of the found differences which will be addressed in a future study.

To reconcile our models with the nucleosynthesis observables of SN 1987A we need to invoke 0.1 ${{{\rm M}}}_{\odot }$ of fallback (see Section 3.4.2). The variation in the amount of synthesized Ni isotopes between runs obtained with different PUSH parameters (Figure 12) suggests that a smaller ${t}_{\mathrm{rise}}$ (and consequently smaller ${k}_{\mathrm{push}}$) could also be compatible with SN 1987A observables, if a larger fallback is assumed. On the one hand, assuming that ${t}_{\mathrm{rise}}$ ranges between 50 and 250 ms, fallback for the 18.0 ${{{\rm M}}}_{\odot }$ model compatible with observations is between 0.14 ${{{\rm M}}}_{\odot }$ (for ${t}_{\mathrm{rise}}=50$ ms) and 0.09 ${{{\rm M}}}_{\odot }$ (for ${t}_{\mathrm{rise}}=250$ ms). On the other hand, if the amount of fallback has been fixed, the observed yields (especially of ⁵⁶Ni) reduce the uncertainty in ${t}_{\mathrm{rise}}$ to $\lesssim 50$ ms.

Our choice of 0.1 ${{{\rm M}}}_{\odot }$ is compatible with the fallback obtained by Ugliano et al. (2012) in exploding spherically symmetric models for progenitor stars in the same ZAMS mass window. Moreover, Chevalier (1989) estimated a total fallback around 0.1 ${{{\rm M}}}_{\odot }$ for SN 1987A, which is supposed to be an unusually high value compared to "normal" SNe II. Recent multi-dimensional numerical simulations by Bernal et al. (2013) and Fraija et al. (2014) confirmed this scenario and furthermore showed that such a hypercritical accretion can lead to a submergence of the magnetic field, giving a natural explanation why the neutron star (possibly) born in SN 1987A has not been found yet.

From the observational side, the compact remnant in SN 1987A is still obscure. From the neutrino signal (see, e.g., Arnett et al. 1989; Koshiba 1992; Vissani 2015 for a recent detailed analysis) one can conclude that a PNS star was formed and that it lasted at least for about 12 s. The mass cut in our calibration run is located at an enclosed baryon mass of 1.56 ${{{\rm M}}}_{\odot }$ without fallback. If we include the 0.1 ${{{\rm M}}}_{\odot }$ of late-time fallback required to fit the observed nickel yields and the explosion energy, we have a final baryonic mass of 1.66 ${{{\rm M}}}_{\odot }$. For the employed HS(DD2) EOS this corresponds to a gravitational mass of a cold neutron star of 1.42 ${{{\rm M}}}_{\odot }$ (without fallback) or 1.50 ${{{\rm M}}}_{\odot }$ (with fallback). The CCSN simulations with artificial explosions of Thielemann et al. (1990), where a final kinetic energy of 1 Bethe was obtained by hand and where the mass cut was deduced from a ⁵⁶Ni yield of $(0.07\pm 0.01)$ ${{{\rm M}}}_{\odot }$, lead to a similar baryonic mass of $(1.6\pm 0.045)$ ${{{\rm M}}}_{\odot }$. These authors also wrote that uncertainties in the stellar models could increase this value to 1.7 ${{{\rm M}}}_{\odot }$ which would also be fully compatible with our result.

The prediction of the neutron star mass has important consequences. From the observations of Demorest et al. (2010) and Antoniadis et al. (2013) it follows that the maximum gravitational mass of neutron stars has to be above two solar masses. The maximum mass of the HS(DD2) EOS is 2.42 ${{{\rm M}}}_{\odot }$, corresponding to a baryonic mass of 2.92 ${{{\rm M}}}_{\odot }$. If the compact remnant in SN 1987A was a black hole, and not a neutron star, it means that at least ∼0.5 ${{{\rm M}}}_{\odot }$ of additional accreted mass were required, if we just take the two solar mass limit. If we use the maximum baryonic mass of HS(DD2) we even have to accrete ∼1.3 ${{{\rm M}}}_{\odot }$ of additional material. Obviously, if such a huge amount of material would be accreted onto the neutron star, our predictions for the explosion energy and the nucleosynthesis would not apply any more.

Nevertheless, we have the impression that it would be difficult to fit the SN 1987A observables and obtain a black hole as the compact remnant at the same time. For spherical fallback, it is certainly excluded. The only possibility could be a highly anisotropic explosion and aspherical accretion, which we cannot address with our study. To show if such a scenario can be realized remains a task for future multi-dimensional studies. In the 2D simulations of Yamamoto et al. (2013) the remnant mass is decreasing with the explosion energy and an explosion energy above 1 Bethe would result in neutron stars below ∼2 ${{{\rm M}}}_{\odot }$ baryonic mass. Note that Kifonidis et al. (2006) already came to the same conclusion that the formation of a black hole in SN 1987A "is quite unlikely," based on 2D simulations with a 15 ${{{\rm M}}}_{\odot }$ progenitor.

Another possibility was proposed by Chan et al. (2009). These authors argued that the time delay of ∼5 s observed for the neutrino signal by the IMB detector could be related to a collapse to a quark star. Due to the proposed faster neutrino cooling of quark stars, this would give a natural explanation why it has not been observed until today. The end of our simulations is also around 5 s, thus we can make statements about the conditions at which the phase transition to quark matter took place in SN 1987A, if the scenario of Chan et al. (2009) was true. We have a central mass density of $4.56\times {10}^{14}$ g cm⁻³ corresponding to ${n}_{B}^{c}=0.272$ fm⁻³ or ${n}_{B}^{c}=1.83\;{n}_{B}^{0}$, a temperature of 23.2 MeV, and an electron fraction of 0.24. Some simplified models for quark matter predict that the phase transition in symmetric matter is shifted to higher densities compared with supernova conditions (Fischer et al. 2011). Under that hypothesis, a phase transition around 2 ${\rho }_{0}$ and 20 MeV cannot be excluded.

A simpler explanation is given by the possibility that a pulsar in the SN 1987A remnant is simply not (yet) observable. Ögelman and Alpar (2004) and Graves et al. (2005) showed that the non-detection of any compact remnant puts important limits on the magnetic field the NS can have (either unusually low or very high, in the realm of magnetars). Furthermore, for both cases (NS and BH) Graves et al. (2005) put severe constraints on currently ongoing accretion scenarios, e.g., spherical accretion is almost ruled out. Graves et al. (2005) conclude that "it seems unlikely that the remnant of SN 1987A currently harbors a pulsar." Our simulations would be in line with the option of a neutron star with a very low magnetic field or with a "normal" magnetic field which is still (partly) buried in the crust due to the late time fallback, similar to what is observed for neutron stars in binary systems. In this respect, recent high-resolution radio observations of the remnant indicate the presence of a compact source or a pulsar wind nebula (Zanardo et al. 2013, 2014). Future observations will be able to clarify the nature of this emission.

As a byproduct of exploring the 18–21 ${{{\rm M}}}_{\odot }$ window and the fitting procedure to SN 1987A we have found interesting correlations between different quantities, which we will discuss here. In Figure 16, we plot the explosion energy, the explosion time, and the (baryonic) remnant mass as function of the progenitor compactness. The results obtained with the calibrated runs indicate a general trend with progenitor compactness for ${E}_{\mathrm{expl}}$. The explosion time, ${t}_{\mathrm{expl}}$, is almost constant within each the LC and the HC group, while the difference between the two groups is related to the difference between how LC and HC models explode (discussed in Section 3.3). The remnant mass increases with compactness, as expected. Nevertheless, we notice significant deviations from the described trends: for ${E}_{\mathrm{expl}}$ and ${t}_{\mathrm{expl}}$ in the HC sample, for ${M}_{\mathrm{rem}}$ mainly in the LC sample.

Zoom In Zoom Out Reset image size

Figure 17 shows explosion times and explosion energies for all the exploding runs in our sample. We can identify a correlation between ${t}_{\mathrm{expl}}$ and ${E}_{\mathrm{expl}}$ for a given progenitor: the larger ${t}_{\mathrm{expl}}$ the lower is ${E}_{\mathrm{expl}}$. This correlation is more pronounced for the HC models than for the LC models. It means that the explosion in PUSH cannot set in too late, if the observed explosion energy should be achieved.

Zoom In Zoom Out Reset image size

In the context of CCSNe, the heating efficiency η is often defined as the ratio between the volume-integrated, net energy deposition inside the gain region and the sum of the ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ luminosities at infinity:

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\displaystyle \int }_{{V}_{\mathrm{gain}}}\rho \;{\dot{e}}_{{\nu }_{e},{\bar{\nu }}_{e}}{dV}}{{L}_{{\nu }_{e}}+{L}_{{\bar{\nu }}_{e}}},\end{eqnarray} \tag{ 26 }$$

see, e.g., Murphy & Burrows (2008), Marek & Janka (2009), Müller et al. (2012a), Suwa et al. (2013), and Couch & O'Connor (2014). In non-exploding, spherically symmetric simulations, η usually rises within a few tens of milliseconds after core bounce and reaches its maximum around $\eta \sim 0.1$ at $t\approx 100\;\mathrm{ms}$, when the shock approaches its maximum radial extension. As soon as the shock starts to recede and the volume of the gain region decreases, η diminishes quickly to a few percents (see, for example, the long-dashed lines in Figure 18).

Zoom In Zoom Out Reset image size

In multi-dimensional simulations, where the shock contraction is delayed or even not happening, energy deposition is expected to be slightly more efficient (η ∼ 0.10–0.15 at maximum) and to decrease more slowly, within a few hundreds of milliseconds after bounce or at the onset of an explosion (see, for example, Murphy & Burrows 2008; Müller et al. 2012a; Couch 2013; Couch & O'Connor 2014). These differences arise not only because the gain region does not contract, but also because neutrino-driven convection efficiently mixes low and high entropy matter between the neutrino cooling and the heating regions below the shock front. Furthermore, convective motion and SASIs are expected to increase significantly the residence time of fluid particles inside the gain region during which they are subject to intense neutrino heating (see, e.g., Murphy & Burrows 2008; Handy et al. 2014). Since the increase of the particle internal energy is given by the time integral of the energy absorption rate over the residence time, this translates to a larger energy variation (Handy et al. 2014).

In spherically symmetric models, the imposed radial motion does not allow the increase of the residence time. This constraint limits the energy gain of a mass element traveling through the gain region. In models exploded using the light-bulb approximation, a large enough internal energy variation is provided by increasing the neutrino luminosity above a critical value, which depends on the mass accretion rate and on the dimensionality of the model (e.g., Burrows & Goshy 1993; Yamasaki & Yamada 2005; Iwakami et al. 2008; Murphy & Burrows 2008; Iwakami et al. 2009; Nordhaus et al. 2010; Hanke et al. 2012; Couch 2013; Dolence et al. 2013; Handy et al. 2014; Suwa et al. 2014). Since in our model the neutrino luminosities are univocally defined by the cooling of the PNS and by the accretion rate history, we increase the energy gain by acting on the neutrino heating efficiency. This effect can be made visible by defining a heating efficiency that takes the PUSH contribution into account, ${\eta }_{\mathrm{tot}}$:

$$\begin{eqnarray}&&{\eta }_{\mathrm{tot}}=\eta +{\eta }_{\mathrm{push}}=\displaystyle \frac{{\displaystyle \int }_{{V}_{\mathrm{gain}}}\rho \left({\dot{e}}_{{\nu }_{e},{\bar{\nu }}_{e}}+{\dot{Q}}_{\mathrm{push}}^{+}\right){dV}}{{L}_{{\nu }_{e}}+{L}_{{\bar{\nu }}_{e}}}.\end{eqnarray} \tag{ 27 }$$

In Figure 18, we plot ${\eta }_{\mathrm{tot}}$ as a function of time for our SN 1987A calibration model, with PUSH (${k}_{\mathrm{push}}=3.5$) and without it (${k}_{\mathrm{push}}=0$). We first notice that the heating efficiency provided by ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ can differ between exploding (short-thick dashed lines) and non-exploding models (long-thin dashed lines). In the case of the exploding model, PUSH provides an increasing contribution to ${\eta }_{\mathrm{tot}}$. It continues to increase steeply up to $t\approx {t}_{\mathrm{on}}+{t}_{\mathrm{rise}}$, but also later, up to $t\approx {t}_{\mathrm{expl}}$, due to the shock expansion preceding the explosion. Thus, the increasing heating efficiency in our spherically symmetric models can be interpreted as an effective way to include average residence times longer than the advection timescale.

In Figure 19, we collect the average and the maximum heating efficiencies for all the models obtained with the set of parameters resulting from the fit procedure (Table 5). Both the average and the maximum values are computed within the interval ${t}_{\mathrm{on}}\leqslant t\leqslant {t}_{\mathrm{expl}}$. We plot them as a function of the compactness and we distinguish between η and ${\eta }_{\mathrm{tot}}$. The maximum of η is usually realized at $t\approx {t}_{\mathrm{on}}$, while the maximum of ${\eta }_{\mathrm{tot}}$ is reached around $t\approx {t}_{\mathrm{expl}}$ (see also Figure 18). Since the explosion sets in later for HC models, when ${t}_{\mathrm{expl}}\gtrsim {t}_{\mathrm{on}}+{t}_{\mathrm{rise}}$, the PUSH factor ${\mathcal{G}}$ has time to rise to k_push for these models. This increases not only the maximum but also the average ${\eta }_{\mathrm{tot}}$ compared with the LC cases. We notice that all four quantities show a correlation with ${\xi }_{1.75}$, but much weaker in the case of η than in the case of ${\eta }_{\mathrm{tot}}$. Moreover, in the HC region, we recognise deviations from monotonic behaviors which reproduce the irregularities already observed in the explosion properties.

Zoom In Zoom Out Reset image size

In the following, we discuss alternative measures of the explosion energy used in the literature for reasons of comparison. We investigate their behaviors at early simulation times and their general rate of convergence. The diagnostic energy ${E}^{+}(t)$, see e.g., Bruenn et al. (2013), is given by the integral of the specific explosion energy ${e}_{\mathrm{expl}}$ over regions where it is positive (again, excluding the PNS core, see Section 2.5.1). The quantity ${E}^{+}(t)$ is often used in multi-dimensional simulations as an estimate of the explosion energy at early simulation times, see e.g., Buras et al. (2006), Suwa et al. (2010), Müller et al. (2012b), Couch & O'Connor (2014), and Takiwaki et al. (2014).

The overburden ${E}_{\mathrm{ov}}(t)$, see Bruenn et al. (2013), is given by the integral of the specific explosion energy of the still gravitationally bound regions between the expanding shock front and the surface of the progenitor star. If we define ${E}_{\mathrm{ov}}^{+}(t)$ as the sum of the overburden and of the diagnostic energy, we recover a measure of the explosion energy equivalent to the one defined in Equation (16):

$$\begin{eqnarray}&&{E}_{\mathrm{expl}}(t)\equiv {E}_{\mathrm{ov}}^{+}(t)={E}^{+}(t)+{E}_{\mathrm{ov}}(t).\end{eqnarray} \tag{ 28 }$$

For long enough simulation times, all matter above the mass-cut should get positive specific explosion energies, and thus the overburden should approach zero and the diagnostic energy should become equal to the explosion energy ${E}_{\mathrm{expl}}(t)$.

Finally, an upper limit for the explosion energy is obtained by also taking into account the "residual recombination energy" ${E}_{\mathrm{rec}}(t)$ (Bruenn et al. 2013):

$$\begin{eqnarray}&&{E}_{\mathrm{ov},r}^{+}(t)={E}_{\mathrm{ov}}^{+}(t)+{E}_{\mathrm{rec}}(t),\end{eqnarray} \tag{ 29 }$$

where ${E}_{\mathrm{rec}}(t)$ is the energy that would be released if all neutron–proton pairs and all ⁴He recombined to ⁵⁶Ni in the regions of positive specific explosion energy. We call it residual recombination energy to make clear that this is energy which is not liberated in our simulations, in contrast to the energy of the recombination processes which we identified in Section 3.2.

In Figure 20, we investigate the behavior of the diagnostic energy ${E}^{+}(t)$, and we compare it with our estimate of the explosion energy ${E}_{\mathrm{expl}}(t)\equiv {E}_{\mathrm{ov}}^{+}(t)$ and with its upper limit represented by ${E}_{\mathrm{ov},r}^{+}(t)$. We want to emphasize that these quantites are obtained from mass integrals above the time-dependent mass cut, in contrast to most of the energies investigated in Section 3.2, where a fixed mass domain was considered.

Zoom In Zoom Out Reset image size

While ${E}_{\mathrm{ov}}^{+}(t)$ and ${E}_{\mathrm{ov},r}^{+}(t)$ have already saturated to a constant value at $t\approx 1.5$ s, even at $t\approx 4.6$ s the diagnostic energy has not yet converged. ${E}_{\mathrm{ov}}^{+}(t)$ and ${E}_{\mathrm{ov},r}^{+}(t)$ approach their asymptotic values from below, and any late time increase ($t\gtrsim 1.5$ s) is due to the energy carried by the neutrino-driven wind ejected from the PNS surface. On the other hand, ${E}^{+}(t)$ reaches its maximum around $\lesssim 1$ s after ${t}_{\mathrm{expl}}$, when the neutrino absorption and the nuclear recombination have released most of their energy in the expanding shock wave, and then it decreases toward ${E}_{\mathrm{ov}}^{+}(t)$, since matter with negative total specific energy is accreted at the shock. The difference between ${E}_{\mathrm{ov}}^{+}(t)$ and ${E}^{+}(t)$ is mainly represented by the gravitational binding energy of the stellar layers above the shock front. Thus, the rate of convergence of the diagnostic energy depends on the amount of gravitational binding energy contained in the outer envelope of the star and on the relative speed at which the shock propagates inside it. Since the gravitational binding energy of the outer layers is similar between the two explored models, the different rate of convergence depends mostly on the different expansion velocity of the shock wave, which is larger for more energetic HC models.

Yamamoto et al. (2013) found for a 15 ${{{\rm M}}}_{\odot }$ progenitor that the diagnostic energy saturates and thus reaches the asymptotic explosion energy already between 1 and 2 s post-bounce. This difference to what we find is related to the different progenitors used and, in particular, to the different binding energy of the outer envelopes, which is expected to be much smaller for a 15 ${{{\rm M}}}_{\odot }$ progenitor than for a ∼20 ${{{\rm M}}}_{\odot }$ progenitor (see, for example, Figure 5 of Burrows 2013). Nevertheless, we conclude that the diagnostic energy is in general (i.e., without further considerations) not suited to give an accurate estimate of the explosion energy at early times.

A similar fitting to SN 1987A energetics has been done for multi-dimensional simulations (2D and 3D) using a light-bulb scheme for the neutrinos by Handy et al. (2014). As initial conditions they used a post-collapse model based on the 15 ${{{\rm M}}}_{\odot }$ blue supergiant progenitor model of Woosley (1988). Even if they did not provide the corresponding compactness, the values of the accretion rate ($\sim 0.2{{\rm -}}0.3$ ${{{\rm M}}}_{\odot }$ s⁻¹) and of the electron neutrino luminosity ($\sim 1.8{{\rm -}}3.5\times {10}^{51}$ erg s⁻¹) at the onset of the explosion are more compatible with our LC models. In their fitting, only the diagnostic explosion energy ${E}^{+}$ was used at a time of ${t}_{\mathrm{pb}}=1.5$ s when it is expected to have saturated to ${E}_{\mathrm{expl}}$ (see Yamamoto et al. 2013). But no estimates for the nucleosynthesis yields were given. The time when the shock reaches 500 km (which corresponds for us to ${t}_{\mathrm{expl}}$) is significantly lower in their models (90–140 ms after bounce), mainly due to the different extension and evolution of the shock during the first 100 ms after core bounce. A more detailed quantitative comparison (albeit limited by the different dimensionality of the two models) requires to use a more similar progenitor. However, the advection timescale and the mass in the gain region are larger than the corresponding values we have obtained in all our models, as expected from the larger average residence time resulting from multi-dimensional hydrodynamical effects.

Ugliano et al. (2012) also calibrated their spherically symmetric exploding models with the observational constraints from SN 1987A, and used progenitor models identical to the ones we have adopted (Woosley et al. 2002). They also found that the remnant mass and the properties of the explosion exhibit a large variability inside the narrow 18–21 ${{{\rm M}}}_{\odot }$ ZAMS mass window (they even found some non-exploding models). However, they did not find any clear trend with progenitor compactness (for example, their calibration model is represented by the 19.8 ${{{\rm M}}}_{\odot }$ ZAMS mass progenitor which belongs to the LC sample). The explosion timescales for models in the 18–21 ${{{\rm M}}}_{\odot }$ ZAMS mass interval are much longer in their case (${t}_{\mathrm{expl}}\sim 0.3{{\rm -}}1$ s), while their range for the explosion energy (0.6–1.6 Bethe) is relatively compatible with ours (0.4–1.6 Bethe). Clearly, all these differences are related to the numerous diversities between the two models.

A possible relation between explosion properites and progenitor compactness has been first pointed out by O'Connor & Ott (2010), who searched for a minimum enhanced neutrino energy deposition in spherically symmetric models. Similarly to us, they found that more compact progenitors require larger heating efficiency to explode. However, they do not investigate the explosion energy of their models. Moreover, they consider it to be unlikely that a model which requires $\eta \gtrsim 0.23$ (${\xi }_{2.5}\gtrsim 0.45$) will explode in nature. In our analysis, we have interpreted a large neutrino heating efficiency in spherically symmetric models as an effective way to take into account longer residence time inside the gain region. We have pointed out that HC models, characterized by larger ${\eta }_{\mathrm{tot}}$, are required to obtain the observed properties of SN 1987A. However, these models still have ${\xi }_{2.5}\lt 0.45$ and our average heating efficiency are below the critical value of O'Connor & Ott (2010).

A clear correlation between explosion properties and progenitor compactness has been recently discussed by Nakamura et al. (2014). They performed systematic 2D calculations of exploding CCSNe for a large variety of progenitors, using the IDSA to model ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ transport. Due to computational limitations and due to the usage of only a NSE EOS, their simulations were limited to ∼1 s after core bounce. Thus, they could not ensure the convergence of the diagnostic energy and could not directly compare their results with CCSN observables. However, they found trends with compactness similar to the ones we have found in our reduce sample.

Other authors have also compared the predicted explosion energy and Ni yield from their models to the observational constraints. For example, Yamamoto et al. (2013), using the neutrino light-bulb method to trigger explosions in spherical symmetry, found a similar trend between explosion energies and nickel masses as we found (see Table 6). They also compared to a thermal bomb model with similar explosion energies and mass cut, and found that the neutrino heating mechanism leads to systematically larger ⁵⁶Ni yields. They related it to higher peak temperatures, which appear because a greater thermal energy is required to unbind the accreting envelope. They also concluded that the neutrino-driven mechanism is more similar to piston-driven models by comparing with Young & Fryer (2007). The problem of overproducing ⁵⁶Ni is lessened in the 2D simulations of Yamamoto et al. (2013) because of slightly lower peak temperatures and the occurrence of fallback.

The conclusions drawn in Section 3.2 about the contributions of nuclear reactions to the explosion energy are somewhat opposite to what can be found in other works in the literature. For example, Yamamoto et al. (2013) state that the contribution of the nuclear reactions to the explosion energy is comparable to or greater than that of neutrino heating. Furthermore, they identify the recombinations of nucleons into nuclei in NSE as the most important nuclear reactions. However, they also point out that this "recombination energy eventually originates from neutrino heating." We think that this aspect is crucial for understanding the global energetics. Indeed, if we had started the analysis presented in Figure 8 not at bounce but at ${t}_{\mathrm{expl}}$ we would also have identified a strong contribution from the nuclear reactions, given roughly by the difference between $-({E}_{\mathrm{mass}}-{E}_{\mathrm{mass},0})$ at ${t}_{\mathrm{expl}}$ (which is close to the minimum) and the final value. However, as is clear from the figure, roughly the same amount of energy was actually taken from the thermal energy before ${t}_{\mathrm{expl}}$. The dominant net contribution to the explosion energy originates from neutrino heating, as is evident from Figure 7 and as we haved discussed in detail in Section 3.2.