A TWO-PARAMETER CRITERION FOR CLASSIFYING THE EXPLODABILITY OF MASSIVE STARS BY THE NEUTRINO-DRIVEN MECHANISM

Our basic modeling approach follows Ugliano et al. (2012) with a number of improvements. To trigger neutrino-driven explosions in spherically symmetric (1D) hydrodynamic simulations, we use a schematic model of the high-density core of the PNS as the neutrino source (for details, see Ugliano et al. 2012). This analytic description is applied to the innermost 1.1 M_⊙, which are excised from the computational domain, and it yields time-dependent neutrino luminosities that are imposed as boundary values at the contracting, Lagrangian inner grid boundary. On the numerical grid, where neutrino optical depths increase from initially ∼10 to finally several thousand, neutrino transport is approximated by the gray treatment described in Scheck et al. (2006) and Arcones et al. (2007). This allows us to account for the progenitor-dependent variations of the accretion luminosity.

Our approach replaces still uncertain physics connected to the equation of state (EOS) and neutrino opacities at high densities by a simple, computationally efficient PNS core model. The associated free parameters are calibrated by reproducing observational properties of SN 1987A. We emphasize that the neutrino emission is sensitive to the time- and progenitor-dependent mass accretion rate. Not only the accretion luminosity increases for progenitors with higher mass accretion rate of the PNS, but also the neutrino loss of the inner core rises with the accreted mass because of compressional work of the accretion layer on the core. Such dependences are accounted for in our modeling of neutron star (NS) core and accretion.

Our numerical realization improves the treatment by Ugliano et al. (2012) in several aspects. We use the high-density EOS of Lattimer & Swesty (1991) with a compressibility of K = 220 MeV and below ρ = 10¹¹ g cm⁻³ apply an e^±, photon, and baryon EOS (Timmes & Swesty 2000) for nuclear statistical equilibrium (NSE; K. Kifonidis 2004, private communication) with 16 nuclei for T > 7 × 10⁹ K and a 14-species alpha network (including an additional neutron-rich tracer nucleus of iron-group material) at lower temperatures (Müller 1986). The tracer nucleus is assumed to be formed in ejecta with Y_e < 0.49 and thus tracks the ejection of matter with neutron excess, when detailed nucleosynthesis calculations predict little production of ⁵⁶Ni (Thielemann et al. 1996).

The network is consistently coupled to the hydrodynamic modeling and allows us to include the contribution from explosive nuclear burning to the energetics of the SN explosions. The collapse phase until core bounce is modeled with the deleptonization scheme proposed by Liebendörfer (2005), using the Y_e(ρ) trajectory of Figure 1 for the evolution of the electron fraction Y_e as a function of density ρ (B. Müller 2013, private communication). This yields good overall agreement with full neutrino transport results and allows for a very efficient computation of large sets of post-bounce models.

Zoom In Zoom Out Reset image size

We perform collapse and explosion simulations for large progenitor sets of different metallicities, namely: the zero-metallicity z2002 set (30 models with ZAMS masses of 11.0–40.0 M_⊙), the low-metallicity (10⁻⁴ solar) u2002 series (247 models, 11.0–75.0 M_⊙), and the solar-metallicity s2002 series (101 models, 10.8–75.0 M_⊙) of Woosley et al. (2002) plus a 10.0 M_⊙ progenitor (S. E. Woosley 2007, private communication) and a 10.2 M_⊙ progenitor (A. Heger 2003, private communication); the solar-metallicity s2014 (151 models, 15.0–30.0 M_⊙) and sh2014 series (15 models, 30.0–60.0 M_⊙, no mass loss) of Sukhbold & Woosley (2014), supplemented by an additional 36 models with 9.0–14.9 M_⊙; the solar-metallicity s2007 series (32 models, 12.0–120.0 M_⊙) of S. E. Woosley et al. (2007, private communication); and the n2006 series (8 models, 13.0–50.0 M_⊙; Nomoto et al. 2006).

For the core-model parameter calibration we choose five different progenitors, namely, the (red supergiant) model s19.8 of the s2002 series as in Ugliano et al. (2012), and four blue supergiant pre-SN models of SN 1987A: w15.0 (ZAMS mass of 15 M_⊙; Woosley et al. 1988), w18.0 (18 M_⊙, evolved with rotation; S. E. Woosley et al. 2007, private communication), w20.0 (20 M_⊙; Woosley et al. 1997), and n20.0 (20 M_⊙; Shigeyama & Nomoto 1990). Compactness values and explosion and remnant parameters of these models are listed in Table 1.

Table 1.  Calibration Models with Explosion and Remnant Properties

Calibration	${\xi }_{1.5}$ a	${\xi }_{1.75}$ a	${\xi }_{2.0}$ a	${\xi }_{2.5}$ a	texpb	Eexpc	Mejd	Eexp/Mej	${M}_{{}^{56}\mathrm{Ni}}$ e	Mtracerf	Mnsg	Mwindh	Mfbi	tν,90j	Eν, totk
Model					(ms)	(B)	(${M}_{\odot }$)	(B ${M}_{\odot }$−1)	(${M}_{\odot }$)	(${M}_{\odot }$)	(${M}_{\odot }$)	(${M}_{\odot }$)	(${M}_{\odot }$)	(s)	(100 B)
s19.8 (2002)	1.03	0.35	0.22	0.14	750	1.30	12.98	0.100	0.072	0.034	1.55	0.096	0.00298	4.27	3.68
w15.0	0.34	0.09	0.03	0.01	580	1.41	13.70	0.103	0.045	0.046	1.32	0.088	0.00018	5.18	2.81
w18.0	0.76	0.26	0.16	0.10	730	1.25	15.42	0.081	0.056	0.036	1.48	0.081	0.00310	4.16	3.32
w20.0	0.98	0.35	0.18	0.06	620	1.24	17.81	0.070	0.063	0.027	1.56	0.089	0.00168	4.73	3.61
n20.0	0.87	0.36	0.19	0.12	560	1.49	14.84	0.100	0.036	0.052	1.55	0.117	0.00243	3.97	3.48

Notes.

^aCompactness evaluated for a central density of 5 × 10¹⁰ g cm⁻³. For w15.0 the value at core bounce is given because earlier data are not available. ^bPost-bounce time of onset of explosion, when shock expands beyond 500 km (1 B = 1 bethe = 10⁵¹ erg). ^cFinal explosion energy, including binding energy of preshock progenitor. ^dMass ejected in the explosion. ^eEjected ⁵⁶Ni mass produced by explosive burning with late-time fallback taken into account. ^fMass of neutron-rich tracer nucleus ejected in neutrino-driven wind material with neutron excess (fallback is taken into account). ^gFinal baryonic neutron star mass including late-time fallback. ^hNeutrino-driven wind mass measured by mass between gain radius at t_exp and preliminary mass cut before fallback. ⁱFallback mass. ^jEmission time for 90% of the radiated neutrino energy. ^kTotal radiated neutrino energy.

Download table as:  ASCIITypeset image

The calibration aims at producing the explosion energy and ejected ⁵⁶Ni mass of SN 1987A compatible with observations, for which the best values are E_exp = (1.50 ± 0.12) × 10⁵¹ erg (Utrobin 2005), E_exp ∼ 1.3 × 10⁵¹ erg (Utrobin & Chugai 2011), and M_Ni = 0.0723–0.0772 M_⊙ (Utrobin et al. 2014), but numbers reported by other authors cover a considerable range (see Handy et al. 2014 for a compilation). The explosion energy that we accept for an SN 1987A model in the calibration process is guided by the ejected ⁵⁶Ni mass (which fully accounts for short-time and long-time fallback) and a ratio of E_exp to ejecta mass in the ballpark of estimates based on light-curve analyses (see Table 1 for our values).

Because of the "gentle" acceleration of the SN shock by the neutrino-driven mechanism (also in 3D simulations; see Utrobin et al. 2014), it is difficult to produce this amount of ejected ⁵⁶Ni just by shock-induced explosive burning. ${M}_{{}^{56}\mathrm{Ni}}$ in Table 1 mainly measures this component but also contains ⁵⁶Ni from proton-rich neutrino-processed ejecta. However, also neutrino-processed ejecta and the neutrino-driven wind with a slight neutron excess could contribute significantly to the ⁵⁶Ni production. The electron fraction Y_e of these ejecta is set by ν_e and ${\bar{\nu }}_{e}$ interactions and depends extremely sensitively on the properties (luminosities and spectra) of the emitted neutrinos, which our transport approximation cannot reliably predict and which also depend on subtle effects connected to multidimensional physics and neutrino opacities. For these reasons we consider the ⁵⁶Ni as uncertain within the limits set by the true ⁵⁶Ni yield from our network on the low side and, in the maximal case, all tracer material added to that. We therefore provide as possible ⁵⁶Ni production of our models the range ${M}_{{}^{56}\mathrm{Ni}}\leqslant {M}_{{}^{56}\mathrm{Ni}}^{\mathrm{total}}\leqslant {M}_{{}^{56}\mathrm{Ni}}+{M}_{\mathrm{tracer}}$.

Different from Ugliano et al. (2012), we reduce the compression parameter of the NS core model by a linear relation $\zeta ^{\prime} \propto {\xi }_{1.75,{\rm{b}}}$ (the value of ${\xi }_{1.75,{\rm{b}}}$ is measured at core bounce) for progenitors with ≤13.5 M_⊙, i.e., we use the function

$$\begin{eqnarray}&&\zeta ^{\prime} =\zeta \left(\displaystyle \frac{{\xi }_{1.75,{\rm{b}}}}{0.5}\right)\ \ \mathrm{for}\ \ M\leqslant 13.5\;{M}_{\odot },\end{eqnarray} \tag{ 4 }$$

with ζ being the value determined from the SN 1987A calibration for a considered progenitor model of that SN. We note that the values of ${\xi }_{1.75,{\rm{b}}}$ are less than 0.5 for all progenitors below 13.5 M_⊙ and close to 0.5 for M ∼ 13.5 M_⊙, for which reason Equation (4) connects smoothly to the ζ value applied for stars above 13.5 M_⊙ according to the SN 1987A calibration.

The modification of Equation (4) accounts for the reduced burden of the small mass of the accretion layer of these stars with their extremely low compactnesses. Such a modification allows us to reproduce the trend to weak explosions obtained in sophisticated 2D and 3D simulations for low-mass iron-core progenitors (Janka et al. 2012; Melson et al. 2015b; Müller 2015). We point out that the Crab supernova SN 1054 is considered to be connected to the explosion of a ∼10 M_⊙ star (e.g., Nomoto et al. 1982; Smith 2013), and its explosion energy is estimated to be up to only ∼10⁵⁰ erg (e.g., Yang & Chevalier 2015). This fact lends support to the results of recent, self-consistent 1D and multidimensional SN models of ≲10 M_⊙ stars (e.g., Kitaura et al. 2006; Fischer et al. 2010; Melson et al. 2015b), whose low explosion energies and low nickel production agree with the Crab observations.

Although as a consequence of our ζ reduction the explosion times, t_exp, tend to be late for stars in the 10.5–12.5 M_⊙ range (Figure 3), this behavior also seems to be compatible with self-consistent, multidimensional simulations of stellar explosions in the 11–12 M_⊙ range by Müller et al. (2013), Müller (2015), and Janka et al. (2012), where these stars were found to have a long-lasting phase of accretion and simultaneous mass outflow after a quite inert onset of the explosion. A more detailed discussion and justification of our modified treatment of low-ZAMS mass cases will be provided in Section 2.3.4 and can also be found in a follow-up paper by Sukhbold et al. (2015), where the values of all PNS core-model parameters are tabulated for all calibrations.

It must be emphasized that in the context of the present work the detailed treatment of the low-mass stars is not overly important. These stars usually explode fairly easily, independent of the treatment of the PNS core model with the original or with our revised calibration. Therefore, these stars lie far away from the boundary curve that separates exploding from nonexploding models and whose determination will be our main goal in Section 3. For this reason exactly the same separation line is obtained when the core-model parameter values from the SN 1987A calibration are applied to all stars.

In Table 1 and the rest of our paper, time-dependent structural parameters of the stars (like compactness values, μ₄ of Equation (3), the iron-core mass M_Fe) are measured when the stars possess a central density of 5 × 10¹⁰ g cm⁻³, unless otherwise stated. This choice of reference density defines a clear standard for the comparison of stellar profiles of different progenitors (see Appendix A of Buras et al. 2006). Different from the moment of core bounce, which was used in other works, our reference density has the advantage of being still close to the initial state of the pre-collapse models provided by stellar evolution modeling and therefore to yield values of the structural parameters that are more similar to those of the pre-collapse progenitor data. Our calibration model w15.0, however, must be treated as an exception. Because pre-collapse profiles of this model are not available any more, all structural quantities for this case are given (roughly) at core bounce.

2.3.1. Status of "Ab Initio" SN Modeling

2.3.2. Modeling Recipes in Recent Literature

2.3.3. Motivation of Modeling Assumptions of This Work

The analytic NS-core model introduced by Ugliano et al. (2012) in combination with their approximate transport solver for treating the accretion component of the neutrino luminosity, as well as neutrino cooling and heating between the PNS and shock, is an attempt to realize the tight coupling of accretion behavior and explodability in close similarity to what is found in current 2D simulations (e.g., those of Marek & Janka 2009; Müller et al. 2012a, 2012b, 2013; Müller & Janka 2014). Since 1D models with elaborate neutrino physics and a fully self-consistent calculation of PNS cooling miss the critical condition for explosions by far, it is not the goal of Ugliano et al. (2012) to closely reproduce the neutrino emission properties of such more sophisticated calculations. Rather, it is the goal to approximate the combined effects of neutrino heating and multidimensional postshock hydrodynamics by a simple and computationally efficient neutrino-source model, which allows for the fast processing of large progenitor sets including the long-time evolution of the SN explosion to determine also the shock breakout and fallback evolution.

Free parameters in the NS core model and the prescribed contraction behavior of the inner grid boundary are calibrated by matching basic observational features (explosion energy, ⁵⁶Ni yield, total release of neutrino energy) of SN 1987A. This is intended to ensure that the overall properties of the neutrino-source model are anchored on empirical ground. Of course, the setting of the parameter values cannot be unambiguous when only a few elements of a single observed SN are used for deriving constraints. However, the approximate nature of the neutrino source treatment as a whole does not require the perfectly accurate description of each individual model component in order to still contain the essence of the physics of the system like important feedback effects between accretion and outflows and neutrinos, which govern the progenitor-dependent variations of explodability and SN properties. A reasonable interplay of the different components is more relevant than a most sophisticated representation of any single aspect of the neutrino-source model.

In detail, the basic features of our 1D realization of the neutrino-driven mechanism along the lines of Ugliano et al. (2012) are the following.

• 1.  
  The possibility of an explosion is coupled closely to the progenitor-dependent strength and evolution of the post-bounce accretion. This is achieved not only by taking into account the accretion luminosity through the approximate neutrino transport scheme but also through the response of the PNS core to the presence of a hot accretion mantle. The evolution of the latter is explicitly followed in our hydrodynamic simulations, which track the accumulation of the collapsing stellar matter around the inner PNS core. The existence of the mantle layer enters the analytic core model in terms of the parameter m_acc for the mass of this layer and the corresponding accretion rate ${\dot{m}}_{\mathrm{acc}}$.
• 2.  
  The inner 1.1 M_⊙ core of the PNS is cut out and replaced by a contracting inner grid boundary and a corresponding boundary condition in our model. This inner core is considered to be the supranuclear high-density region of the nascent NS, whose detailed physics is still subject to considerable uncertainties. This region is replaced by an analytic description, whose parameters Γ, R_c(t), and n (see Ugliano et al. 2012) are set to the same values for all stars. This makes sense because the supranuclear phase is highly incompressible, for which reason it can be expected that the volume of the core is not largely different during the explosion phase for different PNS masses. Moreover, the neutrino diffusion timescale out of this core is seconds, which implies that its neutrino emission is of secondary importance during the shorter post-bounce phase when the explosion develops. Despite its simplicity, our core treatment still includes progenitor- and accretion-dependent variations through the mass m_acc of the hot accretion mantle of the PNS and the mass accretion rate ${\dot{m}}_{\mathrm{acc}}$, whose influence on the inner core is accounted for in Equations (1)–(4) of Ugliano et al. (2012) for describing the energy evolution of the core model.
• 3.  
  The onset of the explosion is considerably delayed (typically between several hundred milliseconds and about a second) with a slow (instead of abrupt) rise of the explosion energy during the subsequent shock acceleration phase, when an intense neutrino-driven wind ejects matter and delivers power to the explosion. Neutrino heating cannot deposit the explosion energy impulsively, because the ejected matter needs to absorb enough energy from neutrinos to be accelerated outward. The rate of energy input to the explosion is therefore limited by the rate at which matter can be channeled trough the heating region. A long-lasting period (hundreds of milliseconds to more than a second) of increasing energy is characteristic of neutrino-driven explosions (Figure 2) and is observed as gradual growth of the explosion energy also in 2D explosion models, e.g., by Scheck et al. (2006) and Bruenn et al. (2014). To achieve this behavior in our 1D models, the core-neutrino source needs to keep up high neutrino luminosities for a more extended period of time than found in fully self-consistent SN simulations, where the rapid decline of the mass accretion rate at the surface of the iron core and at the interface of silicon and silicon-enriched oxygen layers leads to a strong decrease of the accretion luminosity. The longer period of high neutrino emission is compensated by a somewhat underestimated early post-bounce neutrino luminosity (Figure 2) in order to satisfy the energy constraints set by the total gravitational binding energy of the forming NS.
• 4.  
  The timescale and duration of the growth of the explosion energy in multidimensional models of neutrino-driven SNe are connected to an extended period of continued accretion and simultaneous shock expansion that follows after the revival of the stalled shock (see Marek & Janka 2009). Persistent accretion thereby ensures the maintenance of a significant accretion luminosity, while partial re-ejection of accreted and neutrino-heated matter boosts the explosion energy. In our 1D simulations the physics of such a two-component flow cannot be accurately accounted for. In order to approximate the consequences of this truly multidimensional phase, our 1D models are constructed with two important properties: On the one hand, they refer to a high level of the PNS-core luminosity for about 1 s. On the other hand, they are set up to possess a more extended accretion phase that precedes the onset of the delayed explosion before the intense neutrino-driven wind pumps energy into the explosion. The power and mass loss in this wind are overestimated compared to sophisticated neutrino-cooling simulations of PNSs. However, this overestimation of the wind strength has its justification: The early wind is supposed to mimic the mass ejection that is fed in the multidimensional case by the inflow and partial re-ejection of matter falling toward the gain radius during the episode of simultaneous accretion and shock expansion. The enhanced wind mass counterbalances the extra mass accretion by the PNS during the long phase before the shock acceleration is launched, and this enhanced wind mass is of crucial importance to carry the energy of the neutrino-powered blast.

Zoom In Zoom Out Reset image size

Figure 2 shows the post-bounce evolution of the stalled SN shock, the onset of the explosion, energy evolution, and the time evolution of the neutrino emission properties for three representative progenitors, namely, s14, s21, and s27 of the s2014 series, which explode successfully with the w18.0 calibration. The shock stagnation at a radius of approximately 200 km lasts between ∼700 and 900 ms and can exhibit the well-known oscillatory expansion and contraction phases, which signal proximity to the explosion (see, e.g., Buras et al. 2006; Murphy & Burrows 2008; Fernández 2012). The explosion sets in shortly after M₄ (Equation (2)) has fallen through the shock and well before the mass shell corresponding to ${M}_{4}+0.3\;$ has collapsed. High neutrino luminosities are maintained by high mass accretion rates and, after the onset of the explosion, by a contribution from the core emission (dashed lines in the luminosity panel of Figure 2) that grows until roughly 1 s. The current models underestimate the surface luminosity of heavy-lepton neutrinos compared to more sophisticated simulations because of the chosen modest contraction of the inner boundary of the computational grid (which leads to underestimated temperatures in the accretion layer of the PNS) and because neutrino-pair production by nucleon-nucleon bremsstrahlung is not taken into account. We did not upgrade our treatment in this respect because ν_μ and ν_τ are not of immediate relevance for our study since the explosion hinges exclusively on the heating by ν_e and ${\bar{\nu }}_{e}$.

Replacing the inner core of the PNS by a contracting inner boundary of the computational mesh introduces a number of free parameters, whose settings allow one to achieve the desired accretion and neutrino-emission behavior as detailed above. On the one hand, our model contains parameters for the prescription of the contraction of the grid boundary; on the other hand, there are parameters for the simple high-density core model (see Ugliano et al. 2012 and references therein). While the core-model parameters (Γ, ζ, ${R}_{{\rm{c}}}(t)$, n) regulate the neutrino-emission evolution of the excised, high-density core of the PNS, the prescribed grid-boundary radius, ${R}_{\mathrm{ib}}(t)$, governs the settling of the hot accretion mantle of the PNS. Because of partially compensating influences and dependences, not all of these parameters have a sensitive impact on the outcome of our study. Again, more relevant than a highly accurate description of individual components of the modeling is a reasonable reproduction of the overall properties of the accretion and neutrino emission history of the stalled SN shock and mass-accumulating PNS. For example, the moderate increase of the mean neutrino energies with time and their regular hierarchy ($\langle {\epsilon }_{{\nu }_{e}}\rangle \lt \langle {\epsilon }_{{\bar{\nu }}_{e}}\rangle \lt \langle {\epsilon }_{{\nu }_{x}}\rangle ;$ Figure 2) are not compatible with the most sophisticated current models (see, e.g., Marek & Janka 2009; Müller & Janka 2014). They reflect our choice of a less extreme contraction of the 1.1 M_⊙ shell than found in simulations with soft nuclear EOSs for the core matter, where the PNS contracts more strongly and its accretion mantle heats up to higher temperatures at later post-bounce times (compare Scheck et al. 2006 and see the discussion by Pejcha & Thompson 2015). Our choice is motivated solely by numerical reasons (because of less stringent time-step constraints), but it has no immediate drawbacks for our systematic exploration of explosion conditions in large progenitor sets. Since neutrino-energy deposition depends on ${L}_{\nu }\langle {\epsilon }_{\nu }^{2}\rangle$, the underestimated mean neutrino energies at late times can be compensated by higher neutrino luminosities L_ν of the PNS core.

The neutrino emission from the PNS-core region is parameterized in accordance with basic physics constraints. This means that the total loss of electron-lepton number is compatible with the typical neutronization of the inner 1.1 M_⊙ core, whose release of gravitational binding energy satisfies energy conservation and the virial theorem (see Ugliano et al. 2012). Correspondingly chosen boundary luminosities therefore ensure a basically realistic deleptonization and cooling evolution of the PNS as a whole and of the accretion mantle in particular, where much of the inner-boundary fluxes are absorbed and reprocessed. Again, a proper representation of progenitor-dependent variations requires a reasonable description of the overall system behavior but does not need a very high sophistication of all individual components of the system.

2.3.4. Calibration for Low-ZAMS Mass Range

Agreement with the constraints from SN 1987A employed in our work (i.e., the observed explosion energy, ⁵⁶Ni mass, total neutrino energy loss, and the duration of the neutrino signal) can be achieved with different sets of values of the PNS core-model parameters. Using only one observed SN case, the parameter set is underconstrained and the choice of suitable values is ambiguous. It is therefore not guaranteed that the calibration works equally well in the whole mass range of investigated stellar models.

In particular, stars in the low-ZAMS mass regime (M_ZAMS ≲ 12–13 M_⊙) possess properties that are distinctly different from those of the adopted SN 1987A progenitors and in the mass neighborhood of these progenitors. Stars with M_ZAMS ≲ 12–13 M_⊙ are characterized by very small values of compactness (Equation (1)), binding energy outside of the iron core and outside of M₄ (Equation (2)), and mass derivative μ₄ (Equation (3)). The progenitor of SN 1054 giving birth to the Crab remnant is considered to belong to this mass range, more specifically to have been a star with mass around 10 M_⊙ (Nomoto et al. 1982; Smith 2013; Tominaga et al. 2013). Because of their structural similarities and distinctive differences compared to more massive stars, Sukhbold et al. (2015) call progenitors below roughly 12–13 M_⊙ "Crab-like," in contrast to stars above this mass limit, which they term "SN 1987-like."

Stars below ∼10 M_⊙ were found to explode easily in self-consistent, sophisticated 1D, 2D, and 3D simulations (Kitaura et al. 2006; Janka et al. 2008; Fischer et al. 2010; Wanajo et al. 2011; Janka et al. 2012; Melson et al. 2015b) with low energies (less than or around 10⁵⁰ erg = 0.1 B) and little nickel production (<0.01 M_⊙), in agreement with observational properties concluded from detailed analyses of the Crab remnant (e.g., Yang & Chevalier 2015). We therefore consider the results of these state-of-the-art SN models together with the empirical constraints for the Crab SN as important benchmarks that should be reproduced by our approximate 1D modeling of neutrino-powered explosions.

The results of Ugliano et al. (2012) revealed a problem in this respect, because they showed far more energetic explosions of stars in the low-mass domain than expected on grounds of the sophisticated simulations and from observations of Crab. Obviously, the neutrino-source calibration used by Ugliano et al. (2012) is not appropriate to reproduce "realistic" explosion conditions in stars with very dilute shells around the iron core. Instead, it leads to an overestimated power of the neutrino-driven wind and therefore overestimated explosion energies. In particular, the strong and energetic wind is in conflict with the short period of simultaneous postshock accretion and mass ejection after the onset of the explosion in ≲10 M_⊙ stars. Since the mass accretion rate is low and the duration of the accretion phase is limited by the fast shock expansion, the energetic importance of this phase is diminished by the small mass that is channeled through the neutrino-heating layer in convective flows (Kitaura et al. 2006; Janka et al. 2008; Wanajo et al. 2011; Janka et al. 2012; Melson et al. 2015b). In order to account for these features found in the most refined simulations of low-mass stellar explosions, the neutrino-driven wind power of our parametric models has to be reduced.

We realize such a reduction of the wind power by decreasing the parameter ζ, which scales the compression work exerted on the inner (excised) core of the PNS by the overlying accretion mantle (see Equations (1)–(4) in Ugliano et al. 2012), in proportionality to the compactness parameter ξ_1.75,b, which drops strongly for low-mass progenitors (see Equation (4)). This procedure can be justified by the much lighter accretion layers of such stars, which implies less compression of the PNS core by the outer weight. Such a modification reduces the neutrino emission of the high-density core and therefore the mass outflow in the early neutrino-driven wind. As a consequence, the explosion energy falls off toward the low-mass end of the investigated progenitor sets. This can be seen in the top panel of Figure 3, which should be compared to the top left panel of Figure 5 in Ugliano et al. (2012).

Zoom In Zoom Out Reset image size

The ζ scaling of Equation (4) is introduced as a quick fix in the course of this work and is a fairly ad hoc measure to cure the problem of overestimated explosion energies for low-mass SN progenitors. In Sukhbold et al. (2015) a different approach is taken, in which the final value of the core-radius parameter of the one-zone model describing the supranuclear PNS interior as the neutrino source (see Ugliano et al. 2012) is modified. This procedure can directly be motivated by the contraction of the PNS found in self-consistent cooling simulations with microphysical high-density EOSs and detailed neutrino transport. Mathematically, the modification of the core radius has a similar effect on the core-neutrino emission to the ζ scaling employed here. While we refer the reader to Sukhbold et al. (2015) for details, we emphasize that the consequences for the overall explosion behavior of the low-mass progenitors are very similar for both the ζ reduction and the core-radius adjustment applied by Sukhbold et al. (2015). They lead to considerably lower explosion energies for M_ZAMS ≲ 12 M_⊙ stars and a further drop of the explosion energies below ∼9.5 M_⊙.

As a drawback of this modification, the explosions of some of the low-mass progenitors between ∼10.5 and ∼12.5 M_⊙ set in rather late (>1 s post-bounce; see Figure 3).⁵ This, however, is basically compatible with the tendency of relatively slow shock expansion and late explosions that are also found in sophisticated multidimensional simulations of such stars, which, in addition, reveal long-lasting phases of simultaneous accretion and mass ejection after the onset of the explosion (Müller et al. 2013; Müller 2015). This extended accretion phase has only moderate consequences for the estimated remnant masses because the mass accretion rate of these progenitors reaches a low level of ≲0.05–0.1 M_⊙ s⁻¹ after a few hundred milliseconds post-bounce, and some or even most of the accreted mass is re-ejected in the neutrino-driven wind.

Besides providing information on the explosion energies and the onset times of the explosion (defined by the time the outgoing shock reaches 500 km), Figure 3 also displays for exploding models the ejected masses of ⁵⁶Ni and the iron-group tracer element, the baryonic and gravitational remnant masses, the fallback masses, and the total energies radiated by neutrinos. Overall, these results exhibit features very similar to those discussed in detail by Ugliano et al. (2012) for a different progenitor series and a different calibration model. We point out that the fallback masses in the low-mass range of progenitors were overestimated by Ugliano et al. (2012) owing to an error in the analysis (more discussion will follow in Section 3.5; an erratum on this aspect is in preparation).

Figure 4 shows ξ_2.5 versus ZAMS mass with BH-formation cases indicated by gray and explosions by red bars for the s2002 and s2014 series and all calibrations. The irregular pattern found by Ugliano et al. (2012) for the s2002 progenitors is reproduced and appears similarly in the s2014 set. High compactness ξ_2.5 exhibits a tendency to correlate with BHs. But also other parameters reflect this trend, for example, ξ_1.5, the iron-core mass M_Fe (defined as the core where ${\sum }_{\{i| {A}_{i}\gt 46\}}{X}_{i}\gt 0.5$ for nuclei with mass numbers A_i and mass fractions X_i), and the enclosed mass at the bottom of the O-burning shell. All three of them are tightly correlated; see Figure 5 as well as Figure 4 of Ugliano et al. (2012). Also, high values of the binding energy ${E}_{{\rm{b}}}(m\gt {M}_{\mathrm{Fe}})$ outside of M_Fe signal a tendency for BH formation, because this energy correlates with ξ_2.5 (see Figure 4 in Ugliano et al. 2012). However, for none of these single parameters does a sharp boundary value exist that discriminates between explosions and nonexplosions. For all such choices of a parameter, the BH formation limit tends to vary (nonmonotonically) with M_ZAMS, and in a broad interval of values either explosion or BH formation can happen. Pejcha & Thompson (2015) tried to optimize the choice of M for ξ_M, but even their best case achieved only 88% correct predictions. Since in the cases of ξ_2.5 and ${E}_{{\rm{b}}}(m\gt {M}_{\mathrm{Fe}})$, for example, the threshold value for BH formation tends to grow with higher ZAMS mass, one may hypothesize that a second parameter could improve the predictions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Placing the progenitors in a two-parameter space spanned by ξ_1.5 and ξ_2.5 or, equally good, M_Fe and ξ_2.5, begins to show a cleaner separation of successful and failed explosions: SNe are obtained for small values of ξ_2.5, whereas BHs are formed for high values of ξ_2.5, but the value of this threshold increases with ξ_1.5 and M_Fe. For given ξ_1.5 (or M_Fe) there is a value of ξ_2.5 above which only BHs are formed. However, there is still a broad overlap region of mixed cases.

This beginning separation can be understood in view of the theoretical background of the neutrino-driven mechanism, where the expansion of the SN shock is obstructed by the ram pressure of infalling stellar-core matter and shock expansion is pushed by neutrino-energy deposition behind the shock. For neutrino luminosities above a critical threshold ${L}_{\nu ,\mathrm{crit}}(\dot{M})$, which depends on the mass accretion rate $\dot{M}$ of the shock, shock runaway and explosion are triggered by neutrino heating (see Figure 6, left panel, and, e.g., Burrows & Goshy 1993; Janka 2001, 2012; Murphy & Burrows 2008; Nordhaus et al. 2010; Fernández 2012; Hanke et al. 2012; Pejcha & Thompson 2012; Müller & Janka 2015). M_Fe (or ξ_1.5) can be considered as a measure of the mass M_ns of the PNS as an accretor, which determines the strength of the gravitational potential and the size of the neutrino luminosities. Such a dependence can be concluded from the proportionality ${L}_{\nu }\propto {R}_{\nu }^{2}{T}_{\nu }^{4}$, where R_ν is the largely progenitor-independent neutrinosphere radius and the neutrinospheric temperature T_ν increases roughly linearly with M_ns (see Müller & Janka 2014). On the other hand, the long-time mass accretion rate of the PNS grows with ξ_2.5, which is higher for denser stellar cores. For each PNS mass explosions become impossible above a certain value of $\dot{M}$ or ξ_2.5.

Zoom In Zoom Out Reset image size

If the initial mass cut at the onset of the explosion develops at an enclosed mass M = m(r) of the progenitor, we can choose M = m(r) as a suitable proxy of the initial PNS mass, M_ns. A rough measure of the mass accretion rate $\dot{M}$ by the stalled shock around the onset of the explosion is then given by the mass gradient $m^{\prime} (r)\equiv {dm}(r)/{dr}\quad =\quad 4\pi {r}^{2}\rho (r)$ at the corresponding radius r. This is the case because $m^{\prime} (r)$ can be directly linked to the free-fall accretion rate of matter collapsing into the shock from initial radius r according to

$$\begin{eqnarray}\dot{M} & = & \displaystyle \frac{{dm}}{{{dt}}_{\mathrm{ff}}}=\displaystyle \frac{{dm}(r)}{{dr}}{\left(\displaystyle \frac{{{dt}}_{\mathrm{ff}}(r)}{{dr}}\right)}^{-1}\\ & = & \displaystyle \frac{2m^{\prime} (r)}{{t}_{\mathrm{ff}}[(3/r)-m^{\prime} (r)/m(r)]}\approx \displaystyle \frac{2}{3}\;\displaystyle \frac{r}{{t}_{\mathrm{ff}}}\;m^{\prime} (r),\end{eqnarray} \tag{ 5 }$$

where ${t}_{\mathrm{ff}}=\sqrt{{r}^{3}/[{Gm}(r)]}$ is the free-fall timescale (Suwa et al. 2014) and the last, approximate equality is justified by the fact that ${(m^{\prime} /m)}^{-1}={dr}/d\mathrm{ln}(m)\gg r$ outside of the dense stellar core.

Following the critical-luminosity concept now points the way to further improvements toward a classification scheme of explosion conditions: the ${L}_{\nu }\mbox{--}\dot{M}$ dependence of the neutrino-driven mechanism suggests that the explodability of the progenitors may be classified by the parameters M = m(r) and $\dot{M}\propto m^{\prime} (r)$ because the accretion luminosity ${L}_{\nu }^{\mathrm{acc}}\propto {{GM}}_{\mathrm{ns}}\dot{M}/{R}_{\mathrm{ns}}\propto {Mm}^{\prime} (r)$ accounts for a major fraction of the neutrino luminosity of the PNS at the time of shock revival (Müller & Janka 2014, 2015), and, in particular, it is the part of the neutrino emission that reflects the main progenitor dependence. It is important to note that the time-evolving NS and neutrinospheric radii, ${R}_{\mathrm{ns}}\sim {R}_{\nu }$, are nearly the same for different progenitors and only weakly time dependent when the explosions take place rather late after bounce. This is true for our simulations (where ${t}_{\mathrm{exp}}\gtrsim 0.5\;{\rm{s}}$ with few exceptions; Figure 3), as well as for self-consistent, sophisticated models (see, e.g., Figure 3 in Müller & Janka 2014). In both cases the spread of the NS radii and their evolution between ∼0.4 and 1 s after bounce accounts for less than 25% variation around an average value of all investigated models. In Section 3.6 we will come back to this argument and give reasons why the NS radius has little influence on the results discussed in this work. Moreover, the neutrino loss from the low-entropy, degenerate PNS core, whose properties are determined by the incompressibility of supranuclear matter, should exhibit a progenitor dependence mostly through the different weight of the surrounding accretion mantle, whose growth depends on $\dot{M}$. Such a connection is expressed by the terms depending on m_acc and ${\dot{m}}_{\mathrm{acc}}$ in the neutrino luminosity of the high-density PNS core in Equation (4) of Ugliano et al. (2012). We therefore hypothesize, and demonstrate below, a correspondence of the ${L}_{\nu }\mbox{--}\dot{M}$ space and the Mm'–m' parameter plane and expect that the critical luminosity curve ${L}_{\nu ,\mathrm{crit}}(\dot{M})$ maps to a curve separating BH formation and successful explosions in the ${Mm}^{\prime} \mbox{--}m^{\prime}$ plane.⁶

In our simulations neutrino-driven explosions set in around the time or shortly after the moment when infalling matter arriving at the shock possesses an entropy s ∼ 4. We therefore choose M₄ (Equation (2)) as our proxy of the PNS mass, ${M}_{4}\propto {M}_{\mathrm{ns}}$, and ${\mu }_{4}\equiv m^{\prime} (r){[{M}_{\odot }/(1000\mathrm{km})]}^{-1}{| }_{s=4}$ (Equation (3)) as a corresponding measure of the mass accretion rate at this time, ${\mu }_{4}\propto \dot{M}$. The product ${Mm}^{\prime}$ is therefore repesented by ${M}_{4}{\mu }_{4}$. Tests showed that replacing M₄ by the iron-core mass, M_Fe, is similarly good and yields results of nearly the same quality in the analysis following below (which points to an underlying correlation between M₄ and M_Fe). In practice, we evaluate Equation (3) for μ₄ by the average mass gradient of the progenitor just outside of s = 4 according to

$$\begin{eqnarray}{\mu }_{4} & \equiv & {\left.\displaystyle \frac{{\rm{\Delta }}m/{M}_{\odot }}{{\rm{\Delta }}r/1000\mathrm{km}}\right|}_{s=4}\\ & = & \displaystyle \frac{({M}_{4}+{\rm{\Delta }}m/{M}_{\odot })-{M}_{4}}{[r({M}_{4}+{\rm{\Delta }}m/{M}_{\odot })-r(s=4)]/1000\;\mathrm{km}},\end{eqnarray} \tag{ 6 }$$

with ${\rm{\Delta }}m=0.3\;{M}_{\odot }$ yielding optimal results according to tests with varied mass intervals ${\rm{\Delta }}m$. With the parameters M₄ and μ₄ picked, our imagined mapping between critical conditions in the ${L}_{\nu }\mbox{--}\dot{M}$ and ${Mm}^{\prime} \mbox{--}m^{\prime}$ spaces transforms into such a mapping relation between the ${L}_{\nu }\mbox{--}\dot{M}$ and ${M}_{4}{\mu }_{4}\mbox{--}{\mu }_{4}$ planes as illustrated by Figure 6.

Figure 7 demonstrates the strong correlation of the mass accretion rate $\dot{M}={dm}/{dt}$ with the parameter μ₄ as given by Equation (6) (panel (c)), as well as the tight correlations between M₄μ₄ and the sum of ν_e and ${\bar{\nu }}_{e}$ luminosities (panel (a)) and the summed product of the luminosities and mean squared energies of ν_e and ${\bar{\nu }}_{e}$ (panel (b)). It is important to note that the nonstationarity of the conditions requires us to average the quantities plotted on the abscissas over time from the moment when the s = 4 interface passes through the shock until either the models explode (i.e., the shock radius expands beyond 500 km; open circles) or the mass shell $({M}_{4}+0.3){M}_{\odot }$ has fallen through the shock, which sets an end point to the time interval within which explosions are obtained (nonexploding cases marked by filled circles). The time-averaging is needed not only because of evolutionary changes of the preshock mass accretion rate (as determined by the progenitor structure) and corresponding evolutionary trends of the emitted neutrino luminosities and mean energies. The averaging is necessary, in particular, because the majority of our models develop large-amplitude shock oscillations after the accretion of the s = 4 interface, which leads to quasi-periodic variations of the neutrino emission properties with more or less pronounced, growing amplitudes (see the examples in Figure 2). Panel (d) demonstrates that exploding models (open circles) exceed a value of unity for the ratio of advection timescale, t_adv, to heating timescale, t_heat, in the gain layer, which was considered as a useful critical threshold for diagnosing explosions in many previous works (e.g., Janka & Keil 1998; Janka et al. 2001; Thompson 2001; Thompson et al. 2005; Buras et al. 2006; Marek & Janka 2009; Fernández 2012; Müller et al. 2012b; Müller & Janka 2015). The exploding models also populate the region toward low mass accretion rates (as visible in panels (a)–(c), too), which confirms our observation reported in Section 3.1. In contrast, nonexploding models cluster, clearly separated, in the left, upper area of panel (d), where ${t}_{\mathrm{adv}}/{t}_{\mathrm{heat}}\lesssim 1$ and the mass accretion rate tends to be higher. For the calculation of the timescales we follow the definitions previously used by, e.g., Buras et al. (2006), Marek & Janka (2009), Müller et al. (2012b), and Müller & Janka (2015):

$$\begin{eqnarray}&&{t}_{\mathrm{heat}}=\left({\displaystyle \int }_{{R}_{{\rm{g}}}}^{{R}_{{\rm{s}}}}(e+{\rm{\Phi }})\rho \;{dV}\right){\left({\displaystyle \int }_{{R}_{{\rm{g}}}}^{{R}_{{\rm{s}}}}{\dot{q}}_{\nu }\rho {dV}\right)}^{-1},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{adv}}={\displaystyle \int }_{{R}_{{\rm{g}}}}^{{R}_{{\rm{s}}}}\displaystyle \frac{1}{| {v}_{r}| }\;{dr}.\end{eqnarray} \tag{ 8 }$$

Here the volume and radius integrals are performed over the gain layer between gain radius R_g and shock radius R_s. e is the sum of the specific kinetic and internal energies, Φ the (Newtonian) gravitational potential, ρ the density, ${\dot{q}}_{\nu }$ the net heating rate per unit of mass, and v_r the velocity of the flow. Again, because of the variations of the diagnostic quantities associated with the time evolution of the collapsing star and the oscillations of the gain layer, the mass accretion rate and timescale ratio are time-averaged from the moment when the s = 4 interface passes the shock until either 300 ms later or until the model explodes (shock radius exceeding 500 km).⁷ Panels (e) and (f) lend support to the concept of a critical threshold luminosity in the ${L}_{\nu }\mbox{--}\dot{M}$ (and ${L}_{\nu }\langle {\epsilon }_{\nu }^{2}\rangle \mbox{--}\dot{M}$) space mentioned above. The two plots show a separation of exploding (open circles) and nonexploding (filled circles) models in a plane spanned by the time-averaged values of the preshock mass accretion rate on the one hand and the sum of ν_e and ${\bar{\nu }}_{e}$ luminosities (panel (e)) or the summed product of ν_e and ${\bar{\nu }}_{e}$ luminosities times their mean squared energies on the other hand (panel (f)). Again the same time-averaging procedure as for panel (d) is applied. After the s = 4 progenitor shell passes the shock, the time-averaged conditions of the exploding models reach the lower halves of these panels, whereas the time-averaged properties of the nonexploding models define the positions of these unsuccessful explosions in the upper halves. A separation appears that can be imagined to resemble the critical luminosity curve ${L}_{\nu ,\mathrm{crit}}(\dot{M})$ sketched in the left panel of Figure 6.

Zoom In Zoom Out Reset image size

Because of the strong time dependence of the postshock conditions and of the neutrino emission during the phase of dynamical shock expansion and contraction after the accretion of the s = 4 interface, it is very difficult to exactly determine the critical luminosity curve that captures the physics of our exploding 1D simulations. Although we do not consider such an effort as hopeless if the governing parameters are carefully taken into account (see, e.g., the discussions in Müller & Janka 2015; Pejcha & Thompson 2015) and their time variations are suitably averaged, we think that the clear separation of successful and failed explosions visible in panels (e) and (f) of Figure 7 provides proper and sufficient support for the notion of such a critical curve (or, more general: condition) in the ${L}_{\nu ,\mathrm{crit}}\mbox{--}\dot{M}$ space. Figure 6 illustrates our imagined relation between the evolution tracks of collapsing stellar cores that explode or fail to explode in this ${L}_{\nu }\mbox{--}\dot{M}$ space on the left side and the locations of SN-producing and BH-forming progenitors in the M₄μ₄–μ₄ plane on the right side. The sketched evolution paths are guided by our results in panels (e) and (f) of Figure 7. In the following section we will demonstrate that exploding and nonexploding simulations indeed separate in the M₄μ₄–μ₄ parameter space.

The existence of a separation line between BH-forming and SN-producing progenitors in the M₄μ₄–μ₄ plane is demonstrated by Figure 8, which shows the positions of the progenitors for all investigated model series in this two-dimensional parameter space. For all five calibrations successful explosions are marked by colored symbols, whereas BH formation is indicated by black symbols. The regions of failed explosions are underlaid in gray. They are bounded by straight lines with fit functions as indicated in the panels of Figure 8,

$$\begin{eqnarray}&&{y}_{\mathrm{sep}}(x)={k}_{1}\cdot x+{k}_{2},\end{eqnarray} \tag{ 9 }$$

where $x\equiv {M}_{4}{\mu }_{4}$ and $y\equiv {\mu }_{4}$ are dimensionless variables with M₄ in solar masses and μ₄ computed by Equation (6). The values of the dimensionless coefficients k₁ and k₂ as listed in Table 2 are determined by minimizing the numbers of outliers.

Zoom In Zoom Out Reset image size

Table 2.  BH–SN Separation Curves for All Calibration Models

Calibration Model	k1a	k2a	${M}_{4}$ b	${\mu }_{4}$ b	${M}_{4}{\mu }_{4}$ b
s19.8 (2002)	0.274	0.0470	1.529	0.0662	0.101
w15.0c	0.225	0.0495	1.318	0.0176	0.023
w18.0	0.283	0.0430	1.472	0.0530	0.078
w20.0	0.284	0.0393	1.616	0.0469	0.076
n20.0	0.194	0.0580	1.679	0.0441	0.074

Notes.

^aFit parameters of separation curve (Equation (9)) when x and y are measured for a central stellar density of 5 × 10¹⁰ g cm⁻³. ^bMeasured for a central stellar density of 5 × 10¹⁰ g cm⁻³. ^c ${M}_{4}$ and ${\mu }_{4}$ measured roughly at core bounce, because pre-collapse data are not available.

Download table as:  ASCIITypeset image

The stellar models of all progenitor sets populate a narrow strip in the x–y plane of Figure 8, left panels. BH-formation cases are located in the upper left part of the x–y plane. The inclination of the separation line implies that the explosion limit in terms of μ₄ depends on the value of the M₄μ₄ and therefore a single parameter would fail to predict the right behavior in a large number of cases. Denser cores outside of M₄ with high mass accretion rates (larger μ₄) prevent explosions above some limiting value. This limit grows for more massive cores and thus higher M₄ because not only do larger mass accretion rates hamper shock expansion by higher ram pressure but larger core masses and bigger accretion rates also correlate with an increase of the neutrino luminosity of the PNS as expressed by our parameter M₄μ₄. The evolution tracks of successful explosion cases in the left panel of Figure 6 indicate that for higher ${L}_{\nu }\propto {M}_{\mathrm{ns}}\dot{M}$ the explosion threshold, ${L}_{\nu ,\mathrm{crit}}(\dot{M})$, can be reached for larger values of $\dot{M}$.

Explosions are supported by the combination of a massive PNS, which is associated with a high neutrino luminosity from the cooling of the accretion mantle, on the one hand, and a rapid decline of the accretion rate, which leads to decreasing ram pressure, on the other hand. A high value of M₄ combined with a low value of μ₄ is therefore favorable for an explosion because a high accretion luminosity (due to a high accretor mass ${M}_{\mathrm{ns}}\approx {M}_{4}$) comes together with a low mass accretion rate (and thus low ram pressure and low binding energy) exterior to the s = 4 interface (see Figure 9). Such conditions are met, and explosions occur readily, when the entropy step at the s = 4 location is big, because a high entropy value outside of M₄ correlates with low densities and a low accretion rate. M₄ is usually the base of the oxygen shell and a place where the entropy changes discontinuously, causing (or resulting from) a sudden decrease in density due to burning there. This translates into an abrupt decrease in $\dot{M}$ when the mass M₄ accretes. Figure 14 of Sukhbold & Woosley (2014) shows a strong correlation between compactness ξ_2.5 and location of the oxygen shell. The decrease of the mass accretion rate is abrupt only if the entropy change is steep with mass, for which μ₄ at M₄ is a relevant measure.

Zoom In Zoom Out Reset image size

Progenitors with ${M}_{\mathrm{ZAMS}}\lesssim 22\;{M}_{\odot }$ that are harder to explode often have relatively small values of M₄ and an entropy ledge above s = 4 on a lower level than the entropy reached in more easily exploding stars. The lower neutrino luminosity associated with the smaller accretor mass in combination with the higher ram pressure can prohibit shock expansion in many of these cases. Corresponding to the relatively small values of M₄ and relatively higher densities outside of this mass, these cases stick out from their neighboring stars with respect to the binding energy of overlying material, namely, nonexploding models in almost all cases are characterized by local maxima of ${E}_{{\rm{b}},4}={E}_{{\rm{b}}}(m/{M}_{\odot }\gt {M}_{4})$ (see Figure 9).

In view of this insight, it is not astonishing that exploding and nonexploding progenitors can be seen to start separating from each other in the two-parameter space spanned by M₄ and the average entropy value $\langle s{\rangle }_{4}$ just outside of M₄ (Figure 10). Averaging s over the mass interval [M₄, M₄ + 0.5] turns out to yield the best results. Exploding models cluster toward the side of high $\langle s{\rangle }_{4}$ and low M₄, while failures are found preferentially for low values of $\langle s{\rangle }_{4}$. The threshold for success tends to grow with M₄. However, there is still a broad band where both types of outcomes overlap. The disentanglement of SNe and BH-formation events is clearly better achieved by the parameter set of M₄μ₄ and μ₄, which, in addition, applies correctly not only for stars with ${M}_{\mathrm{ZAMS}}\geqslant 15\;{M}_{\odot }$ but also for progenitors with lower masses.

Zoom In Zoom Out Reset image size

Out of 621 simulated stellar models for the s19.8, w15.0, w18.0, w20.0, and n20.0 calibrations, only 9, 14, 16, 11, and 9 models, respectively, do not follow the behavior predicted by their locations on the one or the other side of the separation line in the M₄μ₄–μ₄ plane (see the zoom-ins in the right column of Figure 8). But most of these cases lie very close to the boundary curve, and their explosion or nonexplosion can be affected by fine details and will certainly depend on multidimensional effects. A small sample of outliers is farther away from the boundary line. The w20.0 calibration is the weakest driver of neutrino-powered explosions in our set and tends to yield the largest number of such more extreme outliers.

These cases possess unusual structural features that influence their readiness to explode. On the nonexploding side of the separation line, model s20.8 of the s2014 series with (M₄μ₄, μ₄) ≈ (0.142, 0.0981) is one example of a progenitor that blows up with all calibrations except w20.0, although it is predicted to fail (see Figure 8). In contrast, its mass neighbor, s20.9 with (M₄μ₄, μ₄) ≈ (0.123,0.085), and its close neighbor in the M₄μ₄–μ₄ space, s15.8 of the s2002 series with (M₄μ₄, μ₄) ≈ (0.140, 0.096), both form BHs as expected. The structure of these pre-SN models in the s = 4 region is very similar to M₄= 1.45, 1.45, and 1.46 for s20.8, s20.9, and s15.8, respectively. Although s20.8 reaches a lower entropy level outside of s = 4 than the other two cases and therefore is also predicted to fail, its explosion becomes possible when the next entropy step at an enclosed mass of 1.77 M_⊙ reaches the shock. This step is slightly farther out (at 1.78 M_⊙) in the s20.9 case and comes much later (at ∼1.9 M_⊙) in the s15.8 model. Both the earlier entropy jump and the lower preceding entropy level enable the explosion of s20.8, because the associated higher density maintains a higher mass accretion rate and therefore higher neutrino luminosity until the entropy jump at 1.77 M_⊙ falls into the shock. The abnormal structure of the progenitor therefore prevents that the explosion behavior is correctly captured by our two-parameter criterion for the explodability.

On the exploding side, model s15.3 of the s2014 series with (M₄μ₄, μ₄) ≈ (0.146,0.0797) is expected to blow up according to the two-parameter criterion, but does not do so for all calibrations (Figure 8). Similarly, s15.0 of the s2007 series with (M₄μ₄, μ₄) ≈ (0.137, 0.0749) fails with the w15.0 and w20.0 calibrations, although success is predicted. We compare their structure with the nearby successful cases of s25.4 (s2014 series, (M₄μ₄, μ₄) ≈ (0.150, 0.0820)), s25.2, s25.5 (both from the s2014 series), and s25.8 (s2002 series), all of which group around (M₄μ₄, μ₄) ≈ (0.143, 0.0783). The successfully exploding models all have similar entropy and density structures, namely, fairly low entropies (s ≲ 3) and therefore high densities up to 1.81–1.82 M_⊙, where the entropy jumps to s ≳ 6. The high mass accretion rate leads to an early arrival of the s = 4 interface at the shock (∼300 ms after bounce) and high accretion luminosity. Together with the strong decline of the accretion rate afterward, this fosters the explosion. In contrast, the two models that blow up less easily have higher entropies and lower densities so that the s = 4 mass shells (at ∼1.8 M_⊙ in s15.0 and at ∼1.82 M_⊙ in the s15.3) arrive at the shock much later (at ∼680 and ∼830 ms post-bounce, respectively), at which time accretion contributes less neutrino luminosity. Moreover, both models have a pronounced entropy ledge with a width of ∼0.05 M_⊙ (s15.0) and ∼0.08 M_⊙ (s15.3) before the entropy rises above s ∼ 5. This ledge is much narrower than in the majority of nonexploding models, where it stretches across typically 0.3 M_⊙ or more. The continued, relatively high accretion rate prohibits shock expansion and explosion. This is obvious from the fact that model s15.0 with the less extended entropy ledge exhibits a stronger tendency to explode and for some calibrations indeed does, whereas s15.3 with the wider ledge fails for all calibrations. Our diagnostic parameter μ₄ to measure the mass derivative in an interval of ${\rm{\Delta }}m=0.3\;{M}_{\odot }$, however, is dominated by the high-entropy level (low-density region) above the ledge and therefore underestimates the mass accretion rate in the ledge domain, which is relevant for describing the explosion conditions. Again the abnormal structure of the s15.0 and s15.3 progenitors prevents our two-parameter classification from correctly describing the explosion behavior of these models.

In Figure 11 colored symbols show the positions of the progenitors of the s2014 series and those of the supplementary low-mass models with ${M}_{\mathrm{ZAMS}}\lt 15\;{M}_{\odot }$ in the x–y-plane relative to the separation lines ${y}_{\mathrm{sep}}(x)$ of exploding and nonexploding cases. In the top panel the color-coding corresponds to M_ZAMS, in the middle panel to the iron-core mass, M_Fe, and in the bottom panel to the binding energy of matter outside of the iron core. M_Fe is taken to be the value provided by the stellar progenitor model at the start of the collapse simulation in order to avoid misestimation associated with our simplified nuclear burning network and with inaccuracies from the initial mapping of the progenitor data. Since we use the pressure profile of the progenitor model instead of the temperature profile, slight differences of the derived temperatures can affect the temperature-sensitive shell burning and thus the growth of the iron-core mass.

Zoom In Zoom Out Reset image size

While low-mass progenitors with small iron cores and low binding energies populate the region toward the lower left corner with significant distance to the separation curve, stars above 20 ${M}_{\odot }$ with bigger iron cores and high binding energies can be mostly found well above the separation curve. However, there are quite a number of intermediate-mass progenitors above the line and higher-mass cases below. In particular, a lot of stars with masses between ∼25 and 30 M_⊙ cluster around y_sep(x) in the x ∼ 0.13–0.15 region. These stars are characterized by M_Fe ∼ 1.4–1.5 M_⊙ and high exterior binding energies. Some of them explode, but most fail (see Figure 9). The ones that group on the unsuccessful side are mostly cases with smaller iron cores, whose neutrino luminosity is insufficient to create enough power of neutrino heating to overcome the ram pressure of the massive infall.

Figure 12 displays the BH-formation cases of the s2014 series without associated SNe by black crosses in the x–y-plane. Successful SN explosions of this series plus additional ${M}_{\mathrm{ZAMS}}\lt 15\;{M}_{\odot }$ progenitors are shown by color-coded symbols, which represent, from top to bottom, the final explosion energy (E_exp, with the binding energy of the whole progenitor taken into account), the explosion time (t_exp, measured by the time the outgoing shock reaches 500 km), the gravitational mass of the remnant (with fallback taken into account), the ejected mass of ⁵⁶Ni plus the tracer element (see Section 2; fallback also taken into account), and the fallback mass. The left panel shows the results of our w18.0 calibration, while the right panel is for the n20.0 calibration. The gravitational mass of the NS remnant, ${M}_{\mathrm{ns},{\rm{g}}}$, is estimated from the baryonic mass, ${M}_{\mathrm{ns},{\rm{b}}}$, by subtracting the rest-mass equivalent of the total neutrino energy carried away in our simulations:

$$\begin{eqnarray}&&{M}_{\mathrm{ns},{\rm{g}}}={M}_{\mathrm{ns},{\rm{b}}}-\displaystyle \frac{1}{{c}^{2}}{E}_{\nu ,\mathrm{tot}}.\end{eqnarray} \tag{ 10 }$$

Our estimates of NS binding energies, ${E}_{\mathrm{ns},{\rm{b}}}={E}_{\nu ,\mathrm{tot}}$, are roughly compatible with the Lattimer & Yahil (1989) fit of ${E}_{\mathrm{ns},{\rm{b}}}^{\mathrm{LY}}=1.5\times {10}^{53}$ ${({M}_{\mathrm{ns},{\rm{g}}}/{M}_{\odot })}^{2}$ $\mathrm{erg}=0.084\;{M}_{\odot }{c}^{2}$ ${({M}_{\mathrm{ns},{\rm{g}}}/{M}_{\odot })}^{2}$. Blue symbols in the middle and bottom panels mark fallback BH-formation cases, for which gravitational and baryonic masses (the latter in parentheses) are listed in the panels. We ignore neutrino-energy losses during fallback accretion for both NS- and BH-forming remnants.

Zoom In Zoom Out Reset image size

Progenitors in Figure 12 that lie very close to the separation line tend to produce weaker explosions that set in later than those of progenitors with a somewhat greater distance from the line. Moreover, there is a tendency of more massive NSs to be produced higher up along the x–y band where the progenitors cluster, i.e., bigger NS masses are made at higher values of x = M₄μ₄. Also, the largest ejecta masses of ⁵⁶Ni and ⁵⁶Ni plus the tracer are found toward the right side of the displayed progenitor band just below the boundary of the BH-formation region.

Fallback masses tend to decrease toward the lower left corner of the x–y-plane in Figure 12, far away from the separation curve, where predominantly low-mass progenitors are located, besides five progenitors around 20 M_⊙, which lie in this region because they have extremely low values of ${E}_{{\rm{b}}}(m\gt {M}_{4})$ (see Figure 9) and exceptionally small values of μ₄ (Figure 11) and develop fast and strong explosions with small NS masses, large masses of ejected ⁵⁶Ni plus the tracer, and very little fallback. Closer to the separation line the fallback masses are higher, but for successfully exploding models they exceed ∼0.05 M_⊙ only in a few special cases of fallback SNe (see Figure 3), where the fallback mass can amount to up to several solar masses. We point out that the fallback masses in particular of stars below ∼20 M_⊙ were massively overestimated by Ugliano et al. (2012). The reason was an erroneous interpretation of the outward reflection of reverse-shock-accelerated matter as a numerical artifact connected to the use of the condition at the inner grid boundary. The reverse shock, which forms when the SN shock passes the He/H interface, travels backward through the ejecta and decelerates the outward-moving matter to initially negative velocities. This inward flow of stellar material, however, is slowed down and reflected back outward by the large negative pressure gradient that builds up in the reverse-shock-heated inner region. With this outward reflection, which is a true physical phenomenon and not a boundary artifact, the matter that ultimately can be accreted by the NS is diminished to typically between some 10⁻⁴ M_⊙ and some 10⁻² M_⊙ of early fallback (see Figure 3).

In a handful of high-mass s2014 progenitors (s27.2 and s27.3 for the w18.0 calibation and s27.4, s29.0, s29.1, s29.2, and s29.6 for the n20.0 set) the SN explosion is unable to unbind a large fraction of the star, so that fallback of more than a solar mass of stellar matter is likely to push the NS beyond the BH-formation limit. Such fallback SN cases cluster in the vicinity of x ∼ 0.13–0.14 and y ∼ 0.080–0.081 on the explosion side of the separation line. In the ZAMS mass sequence they lie at interfaces between mass intervals of successfully exploding and nonexploding models, or they appear isolated in BH-formation regions of the ZAMS mass space (see Figures 4 and 9). Their fallback masses and estimated BH masses are listed in the corresponding panels of Figure 12. Naturally, they stick out also by their extremely low ejecta masses of ⁵⁶Ni and tracer elements, late explosion times (around 1 s post-bounce or later), and relatively low explosion energies (∼0.3–0.5 B).

Basically, the accretion luminosity, which is given by ${L}_{\nu }^{\mathrm{acc}}\propto {{GM}}_{\mathrm{ns}}\dot{M}/{R}_{\mathrm{ns}}$, depends not only on the PNS mass, M_ns, and the mass accretion rate, $\dot{M}$, but also on the PNS radius, R_ns. One may wonder whether our two-parameter criterion is able to capture the essential physics, although we disregard the radius dependence when using x = M₄μ₄ as a proxy of the accretion luminosity.

For this reason, we tested a redefinition of x = M₄μ₄ by including a factor $\langle {R}_{\mathrm{ns}}{\rangle }_{300\;\mathrm{ms}}^{-1}$, i.e., using $\tilde{x}\equiv {M}_{4}{\mu }_{4}\langle {R}_{\mathrm{ns}}{\rangle }_{300\;\mathrm{ms}}^{-1}$ instead.⁸ Doing so, we found essentially no relevant effects on the location of the boundary curve. In fact, the separation of exploding and nonexploding models in the $\tilde{x}$–y-plane is even slightly improved compared to the x–y plane, because $\langle {R}_{\mathrm{ns}}{\rangle }_{300\;\mathrm{ms}}^{-1}$ for most nonexploding models is larger than for the far majority of exploding ones. As a consequence, the nonexploding cases tend to be shifted to the left away from the boundary line, whereas most of the exploding cases are shifted to the right, also increasing their distances to the boundary line. This trend leads to a marginally clearer disentanglement of both model groups near the border between explosion and nonexplosion regions. A subset of the (anyway few and marginal) outliers can thus move to the correct side, while very few cases can become new, marginal outliers. It might therefore even be possible to improve the success rate for the classification of explodability by a corresponding (minor) relocation of the boundary curve. The improvement, however, is not significant enough to justify the introduction of an additional parameter into our two-parameter criterion in the form of $\langle {R}_{\mathrm{ns}}{\rangle }_{300\mathrm{ms}}$, which has the disadvantage of not being based on progenitor properties and whose exact, case-dependent value cannot be predicted by simple arguments.

Because the line separating exploding and nonexploding models did not change in our test, the criterion advocated in this paper captures the basic physics, within the limitations of the modeling. While we report here this marginal sensitivity of our two-parameter criterion to the NS radius as a result of the present study, a finally conclusive assessment of this question would require repeating our set of model calculations for different prescriptions of the time-dependent contraction of the inner boundary of our computational grid. The chosen functional behavior of this boundary movement with time determines how the PNS contraction proceeds during the crucial phase of shock revival. In order to avoid overly severe numerical time-step constraints, which can become a serious handicap for our long-time simulations with explicit neutrino transport over typically 20 s, we follow Ugliano et al. (2012) in using a relatively slow contraction of the inner grid boundary. It would be highly desirable to perform model calculations also for faster boundary contractions, which is our plan for future work. In view of this caveat, the arguments and test results discussed in this section should still be taken with a grain of salt.

On grounds of the discussion of our results in the x–y plane, one can actually easily understand why the definition of the separation curve of exploding and nonexploding models in the present paper did not require us to take into account a possible dependence of the accretion luminosity on the NS radius. Instead, we could safely ignore such a dependence when we coined our ansatz that ${L}_{\nu }\propto {Mm}^{\prime} (r)\propto {M}_{4}{\mu }_{4}$. There are two reasons for this. On the one hand, the separation line in the ${M}_{4}{\mu }_{4}$–μ₄ plane is fairly flat. A variation of the NS radius corresponds to a horizontal shift of the location of data points in the x–y plane. However, with the contraction behavior of the NSs obtained in our simulations, only for (relatively few) points in the very close vicinity of the separation curve is such a horizontal shift sufficiently big to potentially have an influence on whether the models are classified as nonexploding or exploding. On the other hand, the models near the separation line typically blow up fairly late (t_exp ≳ 0.6 s), and the NS masses lie in a rather narrow range between roughly 1.4 and 1.6 M_⊙ for the gravitational mass (see Figure 12). For such conditions the variation of the NS radius is of secondary importance (see Figure 3 in Müller & Janka 2014). Low-mass progenitors with less massive NSs, whose radius (at the same time) can be somewhat larger (Figure 3 in Müller & Janka 2014), however, are located toward the lower left corner of the x–y plane and therefore far away from the separation curve (compare Figures 11 and 12). An incorrect horizontal placement of these cases (due to the omission of the dependence of the neutrino luminosity on the NS radius) does not have any relevance for the classification of these models.

Our result of a complex pattern of NS- and BH-forming cases as a function of progenitor mass was previously found by Ugliano et al. (2012), too, and was confirmed by Pejcha & Thompson (2015). Aside from differences in details depending on the use of different progenitor sets and different SN 1987A calibration models, the main differences of the results presented here compared to those of Ugliano et al. (2012) are lower explosion energies for progenitors with M ≲ 13 M_⊙ (see discussion in Section 2.3.4) and lower fallback masses as mentioned in Section 3.5. Based on simple arguments (which, however, cannot account for the complex dynamics of fallback), Pejcha & Thompson (2015) already expected (in particular for their parameterization (a)) that cases with significant fallback—in the sense that the remnant masses are significantly affected—should be rare for solar-metallicity progenitors. Our results confirm this expectation, although the ZAMS masses with significant fallback are different and less numerous than in the work by Pejcha & Thompson (2015). As pointed out by these authors, fallback has potentially important consequences for the remnant mass distribution, and the observed NS and BH masses seem to favor little fallback for the majority of SNe.

As discussed in detail by Ugliano et al. (2012), our explosion models (as well as parameterization (a) of Pejcha & Thompson 2015) predict many more BH-formation cases and more mass intervals of nonexploding stars than O'Connor & Ott (2011), who made the assumption that stars with compactness ξ_2.5 > 0.45 do not explode. Roughly consistent with our results and those of Ugliano et al. (2012), Horiuchi et al. (2014) concluded on the basis of observational arguments (comparisons of the SN rate with the star formation rate; the red supergiant problem as a lack of SNe IIP from progenitors above a mass of ∼16 M_⊙) that stars with an "average" compactness of ξ_2.5 ∼ 0.2 should fail to produce canonical SNe (see Figure 15 in Sukhbold et al. 2015).

A correlation of explosion energy and ⁵⁶Ni mass as found by Pejcha & Thompson (2015) and Nakamura et al. (2015) (both, however, without rigorous determination of observable explosion energies at infinity) and as suggested by observations (see Pejcha & Thompson 2015 for details) is also predicted by our present results (but not by Ugliano et al. 2012 with their erroneous estimate of the fallback masses and the overestimated explosion energies of low-mass progenitors; see Section 2.3.4). Our models yield low explosion energies and low nickel production toward the low-mass side of the SN progenitors (see also Sukhbold et al. 2015). In contrast, the modeling approach by Pejcha & Thompson (2015) predicts a tendency of lower explosion energies and lower ⁵⁶Ni masses toward the high-mass side of the progenitor distribution (because of the larger binding energies of more massive stars), although there is a large mass-to-mass scatter in all results. This difference could be interesting for observational diagnostics. The modeling approach by Pejcha & Thompson (2015) seems to yield neutrino-driven winds that are considerably stronger, especially in cases of low-mass SN progenitors, than those obtained in the simulations of Ugliano et al. (2012) and, in particular, than those found in current, fully self-consistent SN and PNS cooling models, whose behavior we attempt to reproduce better by the Crab-motivated recalibration of the low-mass explosions used in the present work. Moreover, for the ⁵⁶Ni–E_exp correlation reported by Pejcha & Thompson (2015), a mass interval between ∼14 and ∼15 M_⊙ with very low explosion energies and very low Ni production, which does not exist in our models, also plays an important role. Nakamura et al. (2015) found a positive correlation of the explosion energy and the ⁵⁶Ni mass with compactness ξ_2.5. We can confirm this result for the nickel production but not for the explosion energy (see Sukhbold et al. 2015, Figure 15 there). A possible reason for this discrepancy could be the consideration of "diagnostic" energies by Nakamura et al. (2015) at model-dependent final times of their simulations instead of asymptotic explosion energies at infinity, whose determination requires seconds of calculation and the inclusion of the binding energy of the outer progenitor shells (see Figure 3 and Figure 6 in Sukhbold et al. 2015, for the evolution of the energies in some exemplary simulations).

Since a comprehensive discussion of explosion energies, nickel production, and remnant masses is not in the focus of the present work, we refer the reader interested in these aspects to the follow-up paper by Sukhbold et al. (2015). For the same reason we also refrain here from more extended comparisons to the progenitor-dependent explosion and remnant predictions of other studies, in particular those of O'Connor & Ott (2011), Nakamura et al. (2015), and Pejcha & Thompson (2015). A detailed assessment of the different modeling methodologies and their underlying assumptions would be needed to understand and judge the reasons for differences of the results and to value their meaning in the context of SN predictions based on the neutrino-driven mechanism. Such a goal reaches far beyond the scope of our paper.