PROGENITOR-EXPLOSION CONNECTION AND REMNANT BIRTH MASSES FOR NEUTRINO-DRIVEN SUPERNOVAE OF IRON-CORE PROGENITORS

One of the main goals of the long-standing quest for a better understanding of the supernova (SN) explosion mechanism is the establishment of a theoretical connection between the properties of progenitor stars and those of SNe and their remnants. Understanding the progenitor-remnant systematics not only would mean a better definition of SNe as the end point of massive-star evolution but could also serve as input for continuative theoretical studies like binary population synthesis calculations. Moreover, it would provide indispensable information and the theoretical basis for deducing constraints from observations on the physical processes that trigger and power the onset of the explosion at the center of the SN. In particular, the explosion mechanism plays a crucial role for determining the energetic and dynamical properties of the SN blast wave. It thus controls the conditions for explosive nucleosynthesis and is of central importance for answering long-standing astrophysical questions: Which progenitors give birth to neutron stars (NSs), which ones to black holes (BHs)? How do SN properties depend on the progenitor metallicity? Are the observed masses of compact remnants compatible with modeling results? Do observed SNe provide evidence for more than one explosion mechanism? In particular, considering neutrino heating as the most widely favored process to trigger the explosions (Bethe & Wilson 1985), the specific question arises, which range of SN energies and nickel masses can be explained by the neutrino-driven mechanism?

A popular way to model the blast-wave initiation for studying SN nucleosynthesis and determining compact remnant masses are piston-driven explosions (e.g., Woosley & Weaver 1995; Zhang et al. 2008, and references therein), in which a piston is placed at a chosen position in the pre- or post-collapse stellar core and moved with a prescribed, time-dependent velocity to mimic the infall of the selected mass shell and the explosive expansion of the newly formed SN shock. Besides the initial mass cut—fallback of nonescaping material may shift the final mass cut—also the energy deposited by the piston is prescribed as a free parameter. Alternatively, internal energy "bombs" have been used, in which a defined amount of thermal energy is released by an artificial increase of the temperature in a certain volume (Aufderheide et al. 1991).

Investigating spherically symmetric stellar core collapse by hydrodynamic simulations with simplified neutrino treatment, O'Connor & Ott (2011) identified a single parameter, the "progenitor bounce compactness" ξ_2.5 of the innermost 2.5 M_☉ of a star, whose value above a certain threshold they considered to be indicative for an unlikely explosion and thus BH formation. Progenitors with ξ_2.5 ≳ 0.45 were found to reach the BH formation limit typically within ≲0.8 s after bounce because of their tremendous mass-accretion rates. Such stars also turned out to require a particularly high time-averaged neutrino-heating efficiency to develop an explosion against the huge ram pressure of the infalling core material. Nevertheless, the exact value of ξ_2.5, where the bifurcation of the behavior occurs, appears somewhat arbitrarily chosen and not supported by strong theoretical arguments. In particular, the BH formation limit could already take place at a lower value of ξ_2.5.

Belczynski et al. (2011) and Fryer et al. (2012) made theoretical predictions of the NS/BH mass distribution by applying a simple analytic concept for estimating the explosion time, explosion energy, and fallback mass on the basis of structural properties of the progenitor stars (Fryer 2006). Fryer & Kalogera (2001) derived theoretical BH mass distributions by employing the explosion energy as a function of the progenitor mass from a fit to results of a small set of four two-dimensional (2D) models of neutrino-driven explosions (Fryer 1999). They considered different fallback masses in dependence on the (parameterized) fraction of the energy released by the SN mechanism that is used to unbind the outer layers of the progenitor stars.

In this work we adopt a different approach. We perform hydrodynamical simulations in spherical symmetry for a large set of 101 progenitor models, following core collapse and bounce, postbounce accretion, the possible onset of an explosion and subsequent neutrino-driven wind phase of the cooling proto-neutron star (PNS), and the development of the SN blast wave over hours to days until the fallback mass is determined. In order to power the explosions, we refer to the neutrino-driven mechanism, whose viability seems to be supported by recent 2D simulations with detailed energy-dependent neutrino transport (Marek & Janka 2009) and general relativity (Müller et al. 2012a, 2012b), although an ultimate robustness of the model predictions may require three-dimensional simulations and may hinge on remaining uncertainties of the high-density neutrino opacities and equation of state (EOS). Since in our spherical simulations explosions by neutrino heating must be triggered artificially, we excise the high-density core of the nascent NS and replace it by a time-dependent, analytic cooling model. This core model yields the time evolution of the neutrino emission from the high-density interior of the PNS. Its free parameters are adjusted such that observational characteristics of SN 1987A, i.e., the explosion energy, E_87A ∼ 1.3 × 10⁵¹ erg, and ejected ⁵⁶Ni mass, M_{Ni, 87A} ∼ 0.07 M_☉, are reproduced for progenitor stars around 20 M_☉ (e.g., Tanaka et al. 1999; Utrobin & Chugai 2011). This implies the assumption that SN 1987A was an ordinary case of a neutrino-driven explosion in which, e.g., exceptionally rapid core rotation did not play a crucial role for the collapse dynamics. The NS core model with the same parameter values is then also applied to all other collapsing stars.

Provided that SN 1987A was such a normal case and neutrino-energy deposition was responsible for its explosion, our approach thus allows us to explore the response of other progenitor stars to the input of neutrino energy released from their forming compact remnants, and in the case of successful explosions to determine the blast wave and remnant properties of the stars. In the following we will give more information on our modeling strategy and tools, present our results, and discuss the implications.

We employ the explicit, finite-volume, Eulerian multi-fluid code Prometheus (Fryxell et al. 1989) to solve the one-dimensional (1D) hydrodynamics equations on a spherical grid. We use typically 1000 geometrically spaced radial zones with the outer boundary at a radius R_ob. The latter is placed at a large distance (initially at 150,000 km) that is not directly affected by the simulated stellar collapse and SN explosion. In particular, it is sufficiently large to prevent the SN shock from leaving the computational domain during the simulated time. The advection of nuclear species is treated by the consistent multi-fluid advection scheme of Plewa & Müller (1999), and self-gravity of the stellar plasma includes corrections of the Newtonian gravitational potential for general relativistic effects as described in detail in Scheck et al. (2006) and Arcones et al. (2007).

The inner core of the PNS with a mass of M_c = 1.1 M_☉ and densities well above those of the neutrinospheric layer is excised and replaced by a point mass at the coordinate origin. The shrinking of the PNS is mimicked by a retraction of the closed inner boundary at R_ib together with the radial grid, applying the function defined in Arcones et al. (2007). For the simulations discussed in this paper we use the values Rⁱ_ib = R(1.1 M_☉  at  t ∼ 10 ms p.b.) ∼60 km, R^f_ib = 20 km, and t_ib = 0.4 s for the initial and final boundary radii and the exponential contraction timescale, respectively. The chosen values yield a boundary contraction similar to the "standard case" of Scheck et al. (2006). Hydrostatic equilibrium is assumed at the boundary radius R_ib between excised core and accretion envelope of the nascent NS.

The cooling of the PNS interior and corresponding boundary luminosities are described by the neutrino core-cooling model presented in Section 2.1. The temperatures of the spectra of inflowing neutrinos (assumed to have Fermi–Dirac shape) are set equal to the gas temperature T_i in the innermost grid zone, T_ν(R_ib) = T_i, consistent with the fact that neutrinos in highly optically thick regions are close to local thermal equilibrium. The spectral degeneracy parameters are chosen to be $\eta _{\nu _e} = \eta _{\bar{\nu }_e} = \eta _{\bar{\nu }_x} = 0$ for electron neutrinos ν_e, electron antineutrinos $\bar{\nu }_e$, and heavy-lepton neutrinos ν_x (= ν_μ, $\bar{\nu }_\mu$, ν_τ, $\bar{\nu }_\tau$), respectively. (Using the values for local chemical equilibrium for $\eta _{\nu _e}$ and $\eta _{\bar{\nu }_e}$ does not cause any significant differences.) Neutrino transport and neutrino-matter interactions in the computational domain are approximated as in Scheck et al. (2006) by radial integration of the 1D (spherical), gray transport equations for neutrino number and energy, assuming the neutrino spectra to have Fermi–Dirac shape with a calculated local temperature T_ν.

The high-density plasma is described by the tabulated microphysical EOS of Janka & Müller (1996) including arbitrarily degenerate and arbitrarily relativistic electrons and positrons, photons, and four predefined nuclear species (neutrons, protons, α-particles, and a representative Fe-group nucleus) in nuclear statistical equilibrium (NSE). At temperatures and densities below NSE the thermodynamics are treated by the Helmholtz-EOS of Timmes & Swesty (2000).

The NS cooling and its neutrino emission are followed for 8–25 s until the shock has crossed 10¹⁰ cm. At that time the power of the neutrino-driven wind, which pushes the dense ejecta shell enclosed by the expanding SN shock and the wind-termination shock, has decayed to a dynamically insignificant level. The models are then mapped to a grid extending to larger radii, and the inner boundary is moved from below the neutrinosphere to 10⁹ cm in order to ease the time step constraint. The gravity potential of the central object is thereby corrected for the gas mass in the larger excised region, and free-outflow conditions are applied at the new inner boundary. When the SN blast is followed beyond shock breakout for several days, the star is embedded in a wind-like medium with a density ρ∝r⁻² and a constant temperature of T = 1100 K, and the grid is shifted outward once again to cover radii from R_ib = 10¹⁰ cm to R_ob up to 10¹⁵ cm.

During the long-time simulations of the SN dynamics, we follow approximately the explosive nucleosynthesis by solving a small alpha-reactions network, similar to the network described in Kifonidis et al. (2003). It consists of the 13 α-nuclei, ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti, ⁴⁸Cr, ⁵²Fe, ⁵⁶Ni, and an additional tracer nucleus, which is produced via the reaction rate for ⁵²Fe(α,γ)⁵⁶Ni within grid cells whose electron fraction Y_e is below 0.49. This allows us to keep track of element formation in regions with neutron excess. The network is solved in grid cells whose temperature is between 10⁸ K and 7 × 10⁹ K. We assume that at T > 7 × 10⁹ K all nuclei are photo-disintegrated to α-particles. Such a composition is consistent with the NSE yields that are produced by our network solver in the high-temperature limit because the burning network contains α-particles as the only representatives of light nuclear species. A feedback from the network composition to the EOS and thus to the hydrodynamic evolution is neglected.

2.1. Proto-neutron Star Core Model

The cooling of the PNS core is described by an analytic model that couples the excised core region to the surrounding accretion layer, whose evolution is followed on the computational grid. We assume the dense PNS core with mass M_c and radius R_c to be approximately homogeneous and its EOS to be an ideal Γ law, P = (Γ − 1)e (with P being the pressure, e the internal energy density). Combining the total core energy, E_c = E_g + E_i, and the virial theorem, E_g + 3(Γ − 1)E_i + S = 0, we can replace the integrated internal energy, E_i, and express E_c in terms of the Newtonian gravitational energy, E_g = −2/5GM²_c/R_c, and the surface term, S = −4πR³_cP_s, for pressure P_s at R_c:

$$\begin{equation} E_\mathrm{c} = \frac{3\Gamma -4}{3(\Gamma -1)}\,E_\mathrm{g} - \frac{S}{3(\Gamma -1)} \,. \end{equation} \tag{ 1 }$$

The energetic evolution of the core is given by its loss of neutrino energy and compression work done on its surface as

$$\begin{equation} \dot{E}_\mathrm{c} \equiv \frac{{d}E_\mathrm{c}}{{d}t} = - L_{\nu ,\mathrm{c}} - 4\pi P_\mathrm{s} R_\mathrm{c}^2 \dot{R}_\mathrm{c} \,, \end{equation} \tag{ 2 }$$

where $\dot{E}_\mathrm{c}$ can be computed as time derivative of Equation (1), L_{ν, c} is the total neutrino luminosity, and the second term results from the time derivative of the core volume. Instead of setting P_s equal to the boundary pressure on the hydro grid, we prefer to link it to overall properties of the accretion layer. This prescription is intended to capture the nature of the core–mantle coupling but not to constrain the freedom to tune the parameters of the simple PNS-core model. We therefore consider hydrostatic equilibrium in terms of the mass coordinate m(r), dP/dm = −GM/(4πr⁴), and linearize both sides to obtain

$$\begin{equation} P_\mathrm{s} = \zeta \, \frac{G M_\mathrm{c} m_\mathrm{acc}}{4\pi R_\mathrm{c}^4} \,, \end{equation} \tag{ 3 }$$

where m_acc is the mass of the accretion layer that surrounds the PNS core, ζ > 0 is a numerical factor of order unity, and we assumed P₀ ≪ P_s for the pressure P₀ outside of the accretion layer. Moreover, in performing the time derivative of Equation (1) we assumed M_c and Γ to be constant. Combining Equations (1)–(3), we arrive at

$$\begin{equation} L_{\nu ,\mathrm{c}} = \frac{3\Gamma -4}{3(\Gamma -1)}\, (E_\mathrm{g} + S)\,\frac{\dot{R}_\mathrm{c}}{R_\mathrm{c}} - \frac{\zeta }{3(\Gamma -1)}\,\frac{\delta E_\mathrm{acc}}{\delta t}, \end{equation} \tag{ 4 }$$

with S = −ζ GM_cm_acc/R_c and $\delta E_\mathrm{acc}/\delta t \equiv GM_\mathrm{c}\dot{m}_\mathrm{acc}/R_\mathrm{c}$. While the first term on the right-hand side of Equation (4) describes the luminosity increase due to the deepening of the gravity potential and surface work in the case of PNS compression, the second term accounts for the higher core pressure (and internal energy) needed when the accretion layer grows in mass.

Equation (4) is used to prescribe the boundary luminosities, $L_{\nu _i,\mathrm{c}} \equiv \ell _{\nu _i} L_{\nu ,\mathrm{c}}$, of neutrinos of all kinds with $\ell _{\nu _e} = 0.20$, $\ell _{\bar{\nu }_e} = 0.15$, $\ell _{\nu _x} = 0.1625$. (This choice corresponds to a certain loss of lepton number from the PNS core (see Scheck et al. 2006) and ensures a reasonable evolution of Y_e in the PNS mantle and surface layers but does not have much relevance for the dynamical evolution.)

In our simulations m_acc is taken to be the mass between the inner grid boundary and a density of ρ₀ = 10¹⁰ g cm⁻³ at radius r₀, where we define $\dot{m}_\mathrm{acc} = -4\pi r_0^2 v_0 \rho _0$ (accretion means a velocity v₀ < 0 and $\dot{m}_\mathrm{acc} > 0$). The core radius is assumed to contract according to R_c(t) = R_{c, f} + (R_{c, i} − R_{c, f})/(1 + t)ⁿ with R_{c, i} = Rⁱ_ib and R_{c, f} being the initial and final radius, respectively, and t is measured in seconds. With Γ = 3, n = 3, and ζ = 0.6 a choice of R_{c, f} = 6 km allows us to reproduce E_87A and M_{Ni, 87A} of SN 1987A for progenitors in the 20 M_☉ range. For the simulations discussed below the 19.8 M_☉ progenitor serves for the calibration of the PNS-core model, but the overall results are similar when neighboring stars or an SN 1987A blue supergiant progenitor (Woosley et al. 1988) are used for the calibration.

The chosen parameter values lead to typical PNS neutrino-cooling times (t₉₀ for 90% of the total neutrino-energy release) of 3.5–5.5 s (Figure 1). This is shorter than the ∼10 s of emission inferred from the SN 1987A neutrino events of Kamiokande II. However, this detector reported a 7 s gap after eight events in the first two seconds, and the last three events were very close to the detection threshold (Hirata et al. 1987). It is interesting to note that the neutrino signal in all three experiments (Kamiokande II, Irvine–Michigan–Brookhaven, and Baksan) is compatible with a PNS cooling period (exponential cooling timescale) of only 4–5 s (Raffelt 1996; Loredo & Lamb 2002; Pagliaroli et al. 2009). We emphasize that our 1D simulations, which ignore postshock convection and yield standard values for the mean energies of radiated ν_e and $\bar{\nu }_e$ (〈〉 ≈ 10–17 MeV), naturally require overestimated postbounce neutrino fluxes to trigger explosions. Since the neutrino emission of our PNS model is roughly compatible with the SN 1987A data, we expect that our calibration by SN 1987A explosion properties reflects trends of neutrino-driven explosions whose validity for the investigated progenitors in relative comparison holds beyond the considered 1D setup. We also stress that our description of the core behavior accounts for the presence of an accretion layer surrounding the (excised) high-density core. The growth of this layer depends on the progenitor-specific stellar structure and corresponding mass-infall rate. Its accretion luminosity, which adds to the analytically modeled core-neutrino luminosity, is also included in our simulations by means of a simple, gray transport approximation.

Zoom In Zoom Out Reset image size

The progenitor structure therefore influences the postbounce evolution in different ways. The density profile of the outer Fe-core, Si-, and O-layers not only determines the mass-accretion rate and thus the accretion luminosity but also governs the shock stagnation radius because high/low accretion rates damp/favor shock expansion. Moreover, according to Equation (4) our PNS-core model connects the neutrino emission of the PNS core to the growth of the accretion mantle of the PNS. Core and accretion luminosities conspire in reviving and powering the shock wave, for which reason the neutrino emission properties of the PNS (core and accretion contributions) are crucial for the explosion. It is therefore clear that the calibration of the PNS-core parameters by the explosion properties of a chosen progenitor model cannot be independent of the structure of this progenitor star. The fact that basic aspects of our results exhibit some robustness against variations of the calibration model is assuring. However, this does not imply that the predicted explosion properties could not change if the core structure of the SN 1987A progenitor was considerably different from the ∼20 M_☉ solar-metallicity stars and the blue supergiant model used in our calibration tests. Metallicity effects, interaction in a binary scenario, or rotation and convective boundary effects could have affected the mass-loss history and thus could have altered the growth of the stellar core compared to the progenitor models considered by us for calibration.

2.2. Progenitor Stars

Density profiles versus radius and mass for a subset of the investigated 101 solar-metallicity progenitors with iron cores (Woosley et al. 2002) are displayed in Figure 2. The models are given in 0.2 M_☉ steps between 10.8 M_☉ and 28 M_☉ and further up to 40 M_☉ in 1.0 M_☉ steps. We also include two stars with 10.0 M_☉ and 10.2 M_☉. While the density profiles are labeled with the zero-age main-sequence (ZAMS) masses, the corresponding pre-collapse masses are visible in Figure 3. The images show that mass loss is moderate up to 19.6 M_☉ (except for the 10.0 M_☉ case) but grows considerably beyond that value.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The nonmonotonic star-to-star scatter is a consequence of the turn-on and -off of convective shells during the late evolution stages. The scatter can be seen in Figure 4 for different structural parameters as functions of the ZAMS mass: the iron-core mass, M_Fe, is identified with the deleptonized core, i.e., the volume where the electron faction Y_e < 0.497; the binding energy of the overlying shells, E_b, is evaluated at the onset of core collapse; and the enclosed mass at the inner edge of the oxygen-burning shell, $M_{\mathrm{O}_<}$, is defined by the point where the dimensionless entropy per nucleon reaches a value of 4. Finally, the compactness parameter ξ_2.5 is defined as in Equation (10) of O'Connor & Ott (2011) by the ratio of mass M = 2.5 M_☉ and radius R(M) that encloses this mass:

$$\begin{equation} \xi _{2.5} \equiv \frac{M/M_\odot }{R(M)/1000\,\mathrm{km}}\,, \end{equation} \tag{ 5 }$$

where we consider M as the baryonic mass. For our solar-metallicity progenitors it makes no noticeable difference whether we compute ξ_2.5 at the onset of collapse or at the time of bounce (as suggested by O'Connor & Ott 2011), in contrast to the situation for low-metallicity stars, which experience less mass loss during their evolution and which develop more compact cores.

Zoom In Zoom Out Reset image size

The gray histogram bars mark the cases where we found BH formation without an SN explosion. No clear correlation can be observed with extrema of any of the parameters that characterize the density profile of the progenitor core. (We note that the binding energy at bounce instead of the value at the onset of collapse, though different in detail, does not reveal any better correlation.) While many of the BH formation cases are located in a local peak of M_Fe, E_b, $M_{\mathrm{O}_<}$, and ξ_2.5 between ∼22 M_☉ and ∼26 M_☉, there are exceptions in this interval. Moreover, there are BH formation cases below and above this ZAMS mass window with clusterings around 15–16 M_☉, ∼18 M_☉, and 21–22 M_☉. Also the C+O-core and O-shell masses, which exhibit a nearly monotonic increase up to a progenitor ZAMS mass of 28 M_☉ and scattering around a high value beyond, do not show any correlations with BH formation cases.

The clearest indication of possible BH formation seems to be high values of ξ_2.5, at least local maxima (as visible between 15 and 16 M_☉ and around 18 M_☉). But there are exceptions like progenitors near 21 M_☉, which have relatively high values of ξ_2.5 and E_b. Despite this fact, they explode, correlated with possessing lower values of $M_{\mathrm{O}_<}$ than BH formation cases in the close neighborhood. We therefore conclude that the dynamics leading to a successful SN blast does not depend on a single property of the progenitor star but on different ones or even a combination of them.

Although of all quantities displayed in Figure 4 the compactness parameter ξ_2.5 in the form of local maxima exhibits the tightest correlation with BH formation cases, the corresponding threshold value considered by O'Connor & Ott (2011), i.e., ξ_2.5 > 0.45, significantly underestimates the BH formation probability compared to our results. The progenitor variability seen in Figure 4, however, defies the specification of a better and reliable number for the BH formation limit. Instead, there are a wide range of ξ_2.5 values between ∼0.15 on the one side and nearly 0.35 on the other side, for which NS and BH formation cases occur. Of course, any such conclusions can be affected by the 1D nature of our simulations. It cannot be excluded that multi-dimensional effects will shift the ξ_2.5 limit for BH formation to larger values.

The stellar core collapse until ∼10–15 ms after bounce is computed with the Prometheus-Vertex code including sophisticated neutrino transport (Rampp & Janka 2002; Buras et al. 2006a). At that time the bounce shock has converted to a slowly expanding accretion shock. We then continue the simulations of the subsequent accretion phase and the possible explosion with the Prometheus version described in Section 2. The mapping, excision of the NS core, and approximate neutrino treatment do not cause any worrisome transients.

3.1. Explosion Properties

Explosions can develop in the case of a favorable interplay of mass-accretion rate and neutrino luminosities (e.g., Burrows & Goshy 1993; Janka 2001; Fernández 2012). In all successful cases compared to failed explosions of neighboring progenitors, the mass-accretion rate either is lower during a long postbounce period or decreases rapidly when a composition–shell interface arrives at the shock. Shock revival occurs when the neutrino luminosity is still sufficiently high (and thus neutrino heating strong enough) at this time. In a large number of successful and unsuccessful models the decreasing mass-accretion rate triggers shock oscillations, which indicate the proximity to runaway conditions (Buras et al. 2006b; Murphy & Burrows 2008; Fernández 2012) and whose amplification also leads to large-amplitude pulses of the accretion component of the driving neutrino luminosity (see Buras et al. 2006b). In some stars the explosion is fostered by the Si/O interface reaching the shock relatively soon after bounce, either due to its location at a smaller mass coordinate or because of higher mass-accretion rates at earlier times, corresponding to a more compact Si-layer. In this case the high accretion luminosity seems to be supportive. (More information on the time evolution, dynamics, and the neutrino emission of our models will be provided in a separate paper.) In summary, the destiny of a collapsing star does not hinge on a single parameter but depends on the overall structure of the stellar core.

Figure 5 gives an overview of the results of our whole model set. All displayed quantities exhibit considerable scatter even in narrow mass windows, which is a consequence of the nonmonotonicities of the progenitor structure. Failed explosions with BH formation seem to be possible for progenitors below 20 M_☉, and successful SNe with NS formation are also found between 20 and 40 M_☉. While below 15 M_☉ all core collapses produce NSs, the investigated progenitor set yields several "islands" with preferred BH creation above 15 M_☉. A discussion how BH formation cases correlate or do not correlate with the density structure and characteristic quantities of the progenitor cores can be found in Section 2.2.

Zoom In Zoom Out Reset image size

The energies of the neutrino-driven explosions do not exceed 2 × 10⁵¹ erg, and ⁵⁶Ni production up to 0.1–0.15 M_☉ can be expected. Note that our determination of nickel yields is uncertain because the neutrino transport approximation does not provide very accurate information of the Y_e in the neutrino-heated ejecta and in the neutrino-driven wind, where a sizable contribution to the ejected nickel may come from (orange histogram bars in the bottom left panel of Figure 5 denote the corresponding production masses of the neutron-rich tracer nucleus in the network). The onset of the SN blast (roughly defined by the moment the shock passes 500 km) varies in time between ∼0.2 and 1.2 s after bounce, thus including "early" as well as "late" explosions. Later explosions tend to be less energetic because less mass is available in the more dilute gain layer for absorbing energy from neutrinos. This simple relation for each progenitor is partly masked by the large differences of the mass-infall rates in the different stars: a dense stellar core allows for a high mass-accretion rate by the shock and a relatively massive gain layer even at later postbounce times. Moreover, the neutrino-driven wind, which consists of mass shed off the PNS surface by neutrino heating, adds a contribution to the explosion power that does not depend on the density structure of the silicon and oxygen layers.

The blast wave energy of the SN emerges from a combination of different positive contributions: (a) the total (internal plus kinetic plus gravitational) energy of the neutrino-heated postshock matter when it begins expansion away from the gain radius after shock revival; (b) the energy released by the recombination of nucleons and α-particles to heavy nuclei in the expanding postshock material; (c) energy from neutrino absorption in matter that is accreted through the shock and channeled toward the gain radius to be partly ejected again during a phase of simultaneous accretion and outflow after the onset of the explosion; (d) energy carried outward by the neutrino-driven wind from the PNS surface after the accretion has ended; (e) energy released from nuclear burning in the shock-heated ejecta. Two negative effects reduce the energy that finally escapes with the SN ejecta: (f) the binding energy of the stellar layers swept up by the accelerating explosion shock; (g) energy extracted from ejection by mass fallback to the compact remnant. In our explosion models we fully account for points (a), (b), (d), (f), and (g), while point (c) was estimated to be potentially relevant in the multi-dimensional situation (Marek & Janka 2009), and effect (e) typically makes only a small contribution to the SN explosion energy (the burning of 0.1 M_☉ of Si and O releases roughly 10⁵⁰ erg). Also term (a) is usually unimportant because simulations show that the matter in the postshock layer is marginally unbound (i.e., its total energy is near zero) when shock revival sets in and the gas begins outward expansion behind the accelerating blast wave (Scheck et al. 2006; Fernández 2012). The dominant positive contributions to the explosion energy in our 1D models are therefore given by terms (b) and (d). The power of the neutrino-driven wind, point (d), can account for ∼30%–70% of the SN explosion energy (with higher relative importance in more energetic explosions; see the detailed analysis and especially Figure C.5 in Scheck et al. 2006).⁴

3.2. Remnant Properties

The NS baryonic masses are in the range of ∼1.2–2 M_☉. For the observable gravitational masses these numbers need to be reduced by the mass defect associated with the total neutrino-energy loss during NS formation, corresponding to ΔE_tot/c². Our model estimates of ΔE_tot are given in the middle right panel of Figure 5, implying corrections roughly between 0.11 M_☉ and 0.23 M_☉ for the least massive up to the most massive NSs. However, these numbers are approximative—though in a reasonable range—and exact binding energies depend on the NS EOS and have to be determined from detailed general relativistic NS structure models. The smallest BH contains 6.5 M_☉ and forms in the 37 M_☉ star by late fallback of 4.5 M_☉. Since the progenitor has a mass of ∼9.7 M_☉ at collapse, the SN explosion ejects roughly 3.2 M_☉. All other BHs originate from failed explosions and absorb the total mass of the progenitor at collapse (>8.5 M_☉). Although the remnant mass follows a trend of growth with the enclosed mass at the bottom of the oxygen-burning shell (Figure 5), this location is not a reliable indicator for the fate of the star because some models with relatively small silicon core do not explode. However, BH formation clusters in regions where local maxima of the NS mass occur.

Note that we do not terminate the neutrino emission from the excised core when the mass of the compact remnant exceeds the BH formation limit. Instead, the accreting central object is considered to remain gravitationally stable, for which reason its total energy loss can reach ΔE_tot > 5 × 10⁵³ erg (middle right panel of Figure 5) and its cooling time can become t₉₀ ≳ 6 s (Figure 1). However, none of our successful explosions were triggered by this unphysical, long-time release of neutrino energy because all NSs have a baryonic mass of ⩽2 M_☉. This also holds for the compact object that is temporarily present (i.e., until massive fallback occurs long after the launch of the explosion) in the 37 M_☉ progenitor. The value of the baryonic mass of ⩽2 M_☉ is safely below the gravitational mass limit of 1.97 ± 0.04 M_☉ set by the recent precision measurement of Demorest et al. (2010), according to which objects up to at least this mass do not collapse to BHs.

A gap that seems to be present in the observed remnant distribution between NS and BH masses (e.g., Özel et al. 2010; Valentim et al. 2011; Farr et al. 2011) is compatible with our results. This is visible in Figure 6, which displays the normalized (baryonic) remnant mass distribution originating from our neutrino-driven explosion models, weighted by Salpeter's (1955) stellar birthrate function (with power-law index −2.35). We obtain a relatively wide distribution of NS masses (see inset). Note that these are birth masses resulting from single, nonrotating (or at most slowly rotating) stars, for which effects like mass transfer and mass loss associated with possible binary interaction during the progenitor evolution have been ignored. Accretion of the NSs in binaries will, of course, also modify this NS birth-mass distribution.

Zoom In Zoom Out Reset image size

The fallback mass that adds to the initial NS mass is shown in blue in the middle left panel of Figure 5. Fallback is bigger (up to ∼0.2 M_☉) for the lower-mass progenitors where an extended hydrogen envelope leads to a significant deceleration of the explosion shock and a stronger reverse shock. The distribution of Figure 6 implies that 22.5% of the solar-metallicity progenitors produce BHs at the end of their lives. We point out that in the considered sample of 101 solar-metallicity progenitors only the 37 M_☉ star produces a fallback SN with delayed BH formation. In this case the BH swallows a baryon mass of ∼6.5 M_☉ while the explosion ejects the remaining ∼3.2 M_☉ of the pre-collapse mass of the star.

Compared to the neighboring, exploding stars of 36 M_☉ and 38 M_☉, the 37 M_☉ progenitor possesses a higher enclosed mass m(r) and higher gravitational binding energy E_b(r) in the central 15,000 km. This is reflected by its larger values of E_b(M_Fe), $M_{\mathrm{O}_<}$, and ξ_2.5 in Figure 4, although M_Fe is slightly smaller. For this reason the 37 M_☉ model explodes only very late (t_exp ∼ 1 s) and with low explosion energy (Figure 5). Because of the late onset of the blast wave acceleration, the neutrino energy injected to the explosion and the shock velocity become considerably lower than in the 36 and 38 M_☉ cases, in fact too low to unbind the major part of the star. It is interesting to note that the 39 M_☉ progenitor is even more extreme in all of its characteristic parameters displayed in Figure 4. Correspondingly, it does not explode and collapse to a BH as a whole. In contrast, the 34 M_☉ model is also locally extreme in some of its structural parameters, but its compactness ξ_2.5 is lower and it explodes much earlier than the 37 M_☉ star. The latter therefore defines an intermediate, seemingly rare situation in which the combination of a high binding energy and core compactness of the star and a late onset of the explosion and low explosion energy favors the fallback of a large fraction of the progenitor mass.

We performed simulations of neutrino-powered explosions in spherical symmetry for a large set of solar-metallicity, nonrotating stars with Fe-cores. We calibrated an analytic neutrino-cooling model of the PNS core such that the observed explosion properties of SN 1987A were reproduced for progenitors around 20 M_☉. This yields numbers for the total energy release and cooling timescale of the newly formed NS in rough agreement with the SN 1987A neutrino measurement, thus confirming the overall consistency of our approach.

The neutrino-driven explosions turned out to be fostered by "steps" in the stellar density and entropy profiles, which reduce the mass-infall rate and ram pressure of the accretion flow and thus allow the shock to expand. Explosions set in either by a continuous acceleration of the shock or after a phase of shock oscillations with growing amplitude. The onset times range from only 0.1 s to more than 1 s after core bounce, meaning that, depending on the progenitor, the explosion can occur relatively early or with a considerable delay. This conflicts with the claim that neutrino-driven explosions should generally be launched soon after bounce in order to account for the observed gap between highest NS masses and smallest BH masses (Belczynski et al. 2011).

Reflecting the "erratic" structure differences of the progenitor models, the explosion properties and remnant masses exhibit highly nonmonotonic variability even in narrow ranges of the ZAMS mass. This might explain the considerable differences of the explosion energies and nickel masses inferred for observed SNe from light curve and spectral analyses as a function of the ZAMS mass of the progenitors (e.g., Utrobin & Chugai 2011; Smartt et al. 2009; Smartt 2009) without invoking additional pre-collapse degrees of freedom like rapid core rotation. However, the large variations of explosion and remnant properties even for small changes of the progenitor ZAMS mass should be taken with caution. The mixing-length or mixing-length-plus-overshooting description of convection used in the employed progenitor models is extremely unstable to numerical effects in the late burning stages. While the trends we found in dependence of iron-core mass, core-binding energy, and compactness are likely to be valid, the large structural differences of progenitor models with similar ZAMS masses probably do not reflect physical reality. Instead, late-stage convection is subject to substantial dynamic instability, for which reason there may be stochastic variations in the core properties at any given mass. An exploration of this possibility, however, has to be deferred until multi-dimensional progenitor models, self-consistently evolved to the onset of core collapse, become available.

Neutrino heating seems hardly able to provide explosion energies significantly higher than 2 × 10⁵¹ erg and ⁵⁶Ni yields in excess of 0.1–0.2 M_☉. Observed SNe with higher explosion energies and nickel production (e.g., Tanaka et al. 1999) are therefore likely to be driven by a mechanism that is different from neutrino heating and could involve magnetohydrodynamic processes.

The distribution of NS birth masses (including fallback) spans from ∼1.2 M_☉ to ∼2.0 M_☉ in baryonic mass with a broad maximum between 1.4 M_☉ and 1.7 M_☉. BH formation occurs in certain mass intervals and seems to be possible for progenitors with ZAMS masses not only beyond 20 M_☉ but also between 15 M_☉ and 20 M_☉. With a single exception (the 37 M_☉ progenitor) all BHs swallow the progenitor completely. This suggests that fallback SNe could be rare cases in solar-metallicity environments although roughly 23% of the progenitor population end as BHs. Because of the large mass loss of stars with M_ZAMS > 15M_☉, the fallback is strongest (up to ∼0.2 M_☉) for stars below ∼20 M_☉. Correspondingly, we find a wide gap between the maximum NS birth mass (around 1.8 M_☉ gravitational mass) and the minimum BH mass (roughly 6 M_☉), compatible with conclusions from observed compact remnant masses (e.g., Özel et al. 2010; Valentim et al. 2011; Farr et al. 2011). It should, however, be noted that our predictions do not account for binary effects prior to and after the SN explosion, nor do we include NSs formed from ONeMg-core progenitors. A direct comparison with observations therefore requires caution.

The compactness parameter ξ_2.5 introduced by O'Connor & Ott (2011) provides the best indicative quantity based on the pre-collapse structure for the final fate of a star, i.e., whether it is likely to end as an NS or BH. However, our models fail to explode significantly below the discriminating value of ξ_2.5 = 0.45 considered by O'Connor & Ott (2011). If this result were valid more generally, in particular also for lower-metallicity progenitors, it would have implications for the discussion of gamma-ray burst progenitors (O'Connor & Ott 2011; Dessart et al. 2012). Instead of depending on an exact bifurcation value of ξ_2.5, we see BH formation occurring preferentially (but not exclusively) in regions with local maxima of ξ_2.5 and NS as well as BH formation to be possible for ξ_2.5 between ∼0.15 and ∼0.35. The success or failure of the explosion mechanism obviously depends on the progenitor structure (determining the time-dependent mass-accretion rate, shock radius, and accretion luminosity) in a complex way and cannot be predicted exactly on grounds of a single parameter.

So far we have investigated only solar-metallicity progenitors from one modeling group, but we plan to extend our studies also on sets of progenitors from other groups and for different metallicities. Naturally, the spherical symmetry of our simulations is a major constraint, since the explosion mechanism has been recognized to be intrinsically multi-dimensional (e.g., Janka et al. 2007; Burrows et al. 2007). Nonradial fluid flows in the gain layer have a supportive influence and thus reduce the critical luminosity for driving the explosion (e.g., Janka & Müller 1996; Murphy & Burrows 2008; Nordhaus et al. 2010; Hanke et al. 2012; Murphy & Meakin 2011). With a higher neutrino-heating efficiency lower values of the early postbounce luminosities from our analytic PNS core-cooling model would be needed to reproduce the explosion properties of SN 1987A. This, in turn, would stretch the PNS cooling timescale, which in our simulations tends to be on the lower side of the range compatible with the SN 1987A neutrino data.

Nevertheless, this would only mean a recalibration of the time decay of the core-neutrino source adopted for our modeling, but it is unlikely to fundamentally alter qualitative aspects of the behavior of the different progenitors in relative comparison to each other. Therefore, we are hopeful that our basic conclusions will retain their validity also in the multi-dimensional case even if quantitative changes (for example, lower explosion energies by less powerful neutrino-driven winds in low-mass progenitors) could be the consequence. It also cannot be excluded that some of the extreme sensitivity of 1D explosions to details of the progenitor structure might be moderated or disappear in multi-dimensional models. It is interesting to note, however, that the relatively early explosions of the 11.2 M_☉ and 27 M_☉ progenitors of our model set and the later explosion of progenitors around 15 M_☉ are fully compatible with recent results of 2D models with sophisticated neutrino transport (cf. Marek & Janka 2009; Müller et al. 2012a, 2012b).

Our findings challenge canonical thinking, e.g., of the progenitor masses leading to BH formation (e.g., Zhang et al. 2008; O'Connor & Ott 2011), and they raise doubts about claims in the literature that the possible gap between NS and BH masses requires explosions that set in relatively early after core bounce (Belczynski et al. 2011; Fryer et al. 2012). We have shown that a highly variable onset time of the explosion and little fallback for massive progenitors (i.e., a low occurrence rate of fallback SNe) naturally lead to the mass gap. On the other hand, a long delay of the beginning of the explosion is questioned by the recent Bayesian analysis of the NS mass distribution in double NS systems by Pejcha et al. (2012). These authors conclude that for solar-metallicity stars the explosion most likely develops near the edge of the iron core at a specific entropy of S/(N_Ak_B) ≈ 2.8. It is hard to understand how the explosion mechanism could dynamically couple to this specific pre-collapse value of the entropy, because the core density and entropy profiles are highly variable for the different progenitors, they are modified by the nuclear burning during the collapse phase, and the neutrino-driven mechanism in spherical symmetry as well as in the multi-dimensional case is sensitive to the mass-infall rate rather than the pre-collapse entropy structure. For these reasons the mass-accretion rate and the accretion luminosity exhibit large differences between different progenitors when compared for specific values of the entropy, for which reason a wide range of explosion times appears plausible. Interestingly, however, Pejcha et al. (2012) favor an NS formation scenario with no or very little mass fallback, which is consistent with the tendency of our models compared to piston-driven explosion models. Correspondingly, the probability distribution computed by Pejcha et al. (2012) with our predicted remnant masses has little power for NS masses above ∼1.6 M_☉ even when our fallback masses are included. Nevertheless, the theoretical distribution is shifted by ∼0.1 M_☉ to higher masses relative to the observational double NS data. If this mismatch cannot be attributed to an overestimated delay time of the explosion in our SN models as argued above, it might be connected to the low-number statistics of the empirical data or might point to smaller iron cores than present in the solar-metallicity progenitors of our study.

Certainly, our simulations can mean only a very first step of systematically exploring the consequences of neutrino-powered explosions for a wide variety of progenitors. Nevertheless, our results already shed new light on a number of important astrophysical questions connected to how the properties of SN explosions and their remnants connect to those of the progenitor stars. While the fundamental dependences of explosion and remnant properties on the progenitor structure are likely to possess more general validity, the detailed functional variation with the progenitor mass is subject to a number of caveats. These caveats are, in particular, the employed set of progenitor models with their large structural variations even for similar ZAMS masses, the calibration of the neutrino-cooling models of the PNS core based on a certain choice of the SN 1987A progenitor model, and the constraint of the explosion modeling to spherical symmetry. Future work will have to explore the consequences of these constraints and deficiencies.

We thank Lorenz Hüdepohl and Annop Wongwathanarat for support. At Garching, this work was supported by the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich/Transregio 27 "Neutrinos and Beyond," Sonderforschungsbereich/Transregio 7 "Gravitational-Wave Astronomy," and the Cluster of Excellence EXC 153 "Origin and Structure of the Universe." A.A. acknowledges funding by the Helmholtz-University Young Investigator Grant No. VH-NG-825. The computations were partially performed at the Rechenzentrum Garching.