PUSHing Core-collapse Supernovae to Explosions in Spherical Symmetry. II. Explodability and Remnant Properties

After the hydrostatic burning stages stars more massive than 8 M_⊙ undergo core collapse. Up to about 10 M_⊙ (the limit is uncertain, see, e.g., Poelarends et al. 2008; Woosley & Heger 2015; Doherty et al. 2017), the collapse occurs due to electron capture (EC) reactions on C-burning products in the O–Ne–Mg core, resulting in fast contraction and the transition of the core matter to NSE followed by the explosion of the star (dubbed EC supernovae). Beyond 10 M_⊙ the end of the life of massive stars is initiated by the collapse of the central Fe-core, which occurs after central Si-burning when the core exceeds its temperature-dependent Chandrasekhar mass. This collapse of the stellar core marks the onset of a core-collapse supernova (CCSN), a violent explosion that disrupts the star, allows for nuclei beyond iron to be formed and ejected into the interstellar medium, and ultimately leaves behind a neutron star (NS) or black hole (BH) as a remnant.

The major open questions are related to the explosion mechanism itself and how the transition from a "regular" CCSN with the formation of a NS to the formation of a central BH occurs in nature (if there exists such a clear region of transition, see, e.g., Nomoto et al. 2013). This depends strongly on the properties of the stellar progenitor, among others the compactness, i.e., the central mass concentration. Detailed reviews and recent findings on the present understanding of the end stages of massive stars are given in, e.g., Heger & Woosley (2010), Chieffi & Limongi (2013, 2017), Nomoto et al. (2006, 2013), and Nomoto (2017). So-called hypernovae/long-duration gamma-ray bursts (GRBs) possibly occur after BH formation in the case of fast rotation and strong magnetic fields (see, e.g., MacFadyen & Woosley 1999; MacFadyen et al. 2001). If fast rotation and strong magnetic fields are not present, BH formation leads to a failed or faint supernova (see, e.g., Lovegrove & Woosley 2013). For smaller initial stellar masses a NS with magnetic fields as high as 10¹⁵ Gauss, known as a magnetar, can be formed. In the present paper, we will only address models without rotation and magnetic fields, i.e., neutrino-driven CCSNe.

The full solution to the CCSN problem in a self-consistent way is still not converged yet. Multi-D simulations are a well-suited and necessary tool to investigate the underlying mechanism of CCSNe. There exists a growing set of 2D and 3D CCSN explosion models (see, e.g., Janka 2012; Burrows 2013; Nakamura et al. 2015; Bruenn et al. 2016; Janka et al. 2016; Burrows et al. 2018). Note that while 2D simulations tend to explode more easily, in 3D the resulting explosion energies may be higher, as is shown, e.g., in Takiwaki et al. (2014), Lentz et al. (2015), Melson et al. (2015a, 2015b), Müller (2015), Janka et al. (2016), and Hix et al. (2016). Multidimensional simulations are also computationally intensive. They are prohibitively expensive if large spatial domains are required and/or if simulations are required to follow the evolution for long timescales beyond the onset of the explosion. Thus, given the present status, it is still too early to provide complete predictions from multi-D simulations for an extended sample of progenitor stars.

Computationally more affordable models would be a viable tool for studies of large samples of stars but self-consistent state-of-the-art spherically symmetric simulations of CCSNe do not lead to explosions via the delayed neutrino-driven mechanism (see, e.g., Liebendörfer et al. 2004). The reason is simply the lack of dimensionality, which is the basis for many of the proposed explosion mechanisms, e.g., the convective or SASI-aided neutrino-driven (see references given above) and the magnetorotational mechanism (Winteler et al. 2012; Mösta et al. 2014, 2015; Nishimura et al. 2015, 2017; Halevi & Mösta 2018). An exception is the phase-transition mechanism (Sagert et al. 2009; Fischer et al. 2011), which leads to successful SN explosions in one-dimensional models. Given the complexities and computational costs of multi-D approaches, a number of approximate 1D approaches have been proposed, which try to mimic the net effects of multi-D simulations. Such parameterized explosions induced in spherically symmetric models are still a pragmatic and valuable approach to study large numbers of stellar progenitors, from the onset of the explosion up to several seconds after core bounce.

In the past, some simplified approaches have been used to predict supernova nucleosynthesis, which artificially induce explosions with estimated typical explosion energies (e.g., Woosley & Weaver 1995; Thielemann et al. 1996, and many more). Interim approaches beyond piston or thermal bomb models attempt to mimic multi-D neutrino heating in a spherical approach, in order to obtain more appropriate predictions of the explosion energy, the mass cut between neutron star and ejecta, as well as nucleosynthesis (including the effects of neutrinos on Y_e). This includes the "neutrino light-bulb method"—where the proto-neutron star (PNS) is excised and replaced with an inner boundary condition that contains an analytical prescription for the neutrino luminosities (e.g., Yamasaki & Yamada 2005; Iwakami et al. 2008; Yamamoto et al. 2013)—and the "absorption methods"—where the neutrino energy deposition is increased by modifying the neutrino opacities in spherically symmetric models with detailed Boltzmann neutrino transport (Fröhlich et al. 2006a, 2006b; Fischer et al. 2010).

More recent spherically symmetric approaches try to mimic the effect of multi-D neutrino transport in a way that is adapted more consistently to core collapse and PNS accretion. These approaches still need calibrations that can be provided by comparison with a variety of observations of explosion energies, deduced ejected ⁵⁶Ni masses, progenitor properties, as well as the outcome of multi-D studies. Ugliano et al. (2012) presented a more sophisticated light-bulb method to explode spherically symmetric models, using neutrino energy deposition in post-shock layers. They used an approximate, gray neutrino transport and replaced the innermost 1.1 M${}_{\odot }$ of the PNS by an inner boundary. The evolution of the neutrino luminosity at the boundary was based on an analytic cooling model of the PNS, which depends on a set of free parameters. These parameters (within the so-called Prometheus-Hot Bubble (P-HOTB) approach) are set by fitting observational properties of SN 1987A for progenitor masses around 20 M${}_{\odot }$ (see also Ertl et al. 2016; Sukhbold et al. 2016).

Perego et al. (2015, hereafter our Paper I) utilized the energy in muon and tau neutrinos as an additional energy source that approximately captures the essential effects of multi-D neutrino transport (the PUSH method). PUSH relies on the so-called delayed neutrino-driven mechanism as the central engine of CCSNe. In particular, it provides an artificially enhanced neutrino energy deposition inside the gain region in spherically symmetric models by depositing a fraction of the energy carried by heavy flavor neutrinos behind the shock (see the detailed discussion in Section 2). This allows us to trigger explosions in 1D simulations without modifying ν_e and ${\bar{\nu }}_{e}$ luminosities nor changing charged current reactions. This increases the accuracy of treating the electron fraction for the innermost ejecta, which is a crucial ingredient for nucleosynthesis calculations. The PUSH method is also still parameterized, in the sense that the additional energy deposition is calibrated by comparing the obtained explosion energies and nucleosynthesis yields with observations of nearby supernovae (in particular, SN 1987A).

As mentioned above, a major open question is whether core collapse eventually leads to a supernova explosion with a NS remnant or whether the final outcome is a central BH. Using P-HOTB, Sukhbold et al. (2016) show that both possible outcomes can occur within the same mass interval, mainly dependent on the pre-collapse stellar model and its compactness (O'Connor & Ott 2011). The best observed core-collapse supernova to date is SN 1987A, which is commonly used as the standard against which numerical supernova models are tested. For SN 1987A, observations provide its explosion energy and ejected masses of ${}^{\mathrm{56,57,58}}$Ni and ⁴⁴Ti as well as progenitor mass and metallicity. Besides the well-observed SN 1987A, many other CCSNe are known. In addition, observational evidence for failed CCSN explosions has recently become available. To date, the LIGO-VIRGO collaboration has detected gravitational wave signals from four different BH–BH merger events (Abbott et al. 2016a, 2016b, 2017; The LIGO Scientific Collaboration et al. 2017). The individual pre-merger BH masses in these events range from ${7}_{-2}^{+2}\,{M}_{\odot }$ (GW170608) to ${36}_{-4}^{+5}{M}_{\odot }$ (GW150914). In two of these events, both pre-merger BHs have masses above 25 M_⊙. One possible formation mechanism of these stellar mass BHs is through failed SNe of low-metallicity massive stars (Abbott et al. 2016b). A complementary observation has been recently made by Adams et al. (2017) where they used Hubble Space Telescope imaging to confirm the optical disappearance of a 25 M_⊙ red supergiant. This suggests that stars with initial masses around 25 M_⊙ may end their lives in a failed explosion, ultimately resulting in a BH. In a recent review, Smartt (2015) argues, based on a distance-limited sample, that supernovae in the local universe are, on the whole, not produced by high-mass stars. They interpret their result in such a way that the missing high-mass SN-progenitor stars are the high-compactness stars that die as failed SNe or directly collapse to a BH. However, Davies & Beasor (2018) argue that the cutoff mass for progenitors of type II-P and II L CCSNe may be 25 M_⊙, which is higher than the previously assumed 17 M_⊙ due to a source of systematic error in converting pre-explosion photometry into an initial ZAMS mass. Full-scale multi-D simulations of core-collapse obtaining resulting BH formation for 40 M_⊙ and 70 M_⊙ ZAMS masses (Kuroda et al. 2018; Pan et al. 2018).

Improved modeling approaches in spherical symmetry like P-HOTB and PUSH represent for the first time treatments that provide supernova explosion energies and resulting NS masses, as well as the transition to BH formation as a function of ZAMS mass. This modeling also permits the investigation of the transition from regular CCSNe to faint supernovae, if they exist, or failed supernovae as a function of stellar mass.

In Paper I, the PUSH method was calibrated to SN 1987A, using pre-explosion models representing red giants in the mass range of 18–21 M_⊙ such that the explosion energy and the yields of ^56–58Ni and ⁴⁴Ti are consistent with their corresponding values derived from observations as proof of principle. Here, in Paper II, we discuss how PUSH is extended to a broad range of progenitor masses. We apply the refined PUSH method to large samples of solar metallicity pre-explosion models from Woosley et al. (2002) and Woosley & Heger (2007), with masses between 10.8 and 120 M_⊙ and we present the explosion properties and the progenitor-remnant connection. A detailed discussion of the nucleosynthesis obtained with PUSH is presented in a paper by Curtis et al. (2018, hereafter Paper III).

The paper is organized as follows: In Section 2, our numerical setup and the input models are discussed. The focus is on the changes compared to Paper I. Next, based on general features of CCSN observations, we introduce requirements that the PUSH method should fulfill. These requirements are formulated in the form of constraints on the parameterized neutrino heating in Section 3. In Section 4, we present the resulting systematics of explosion properties across the ZAMS-mass range and trends in compactness. Furthermore, we show the predicted NS and BH mass distributions resulting from this approach in Section 5 and summarize our results in Section 6.

The PUSH method provides a computationally efficient and physically motivated framework to explode massive stars in spherical symmetry. It allows us to study multiple aspects related to CCSNe that require modeling of the explosion for several seconds after its onset and for extended sets of progenitors. The method, based on its calibration, is also well-suited to explore the effects of the shock passage through the star and to predict the NS and BH mass distribution as a function of stellar mass. PUSH relies on the so-called delayed neutrino-driven mechanism as the central engine of CCSNe. In particular, it provides an artificially enhanced neutrino energy deposition inside the gain region in spherically symmetric models, which do not explode in self-consistent simulations. This more efficient energy deposition is inspired by the increase of the net neutrino heating that a fluid element experiences due to the presence of multidimensional effects. Unlike other methods that employ external energy sources or that use modified electron flavor neutrino luminosities to trigger artificial explosions, PUSH deposits a fraction of the energy otherwise carried away by heavy flavor neutrinos (${\nu }_{x}={\nu }_{\mu },{\bar{\nu }}_{\mu },{\nu }_{\tau },{\bar{\nu }}_{\tau }$) behind the shock to ultimately provide successful explosion conditions. In self-consistent core-collapse models, ν_x represent a marginal dependence on the temporal evolution of the accretion rate (e.g., Liebendörfer et al. 2004), and their contribution to the energy deposition inside the gain region is negligible. However, their usage in PUSH presents a number of advantages. The properties of the ν_x emission, which includes dynamical feedback from accretion history as well as cooling properties of the forming compact object, correlate significantly with the main features of the ν_e and ${\bar{\nu }}_{e}$ emission (O'Connor & Ott 2013). Moreover, the accretion luminosity depends not only on the accretion rate but also on the evolution of the mass and radius of the PNS, which is treated accurately and self-consistently in this method. This allows us to trigger explosions in 1D simulations without modifying ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ luminosities or changing charged current reactions. This increases the accuracy of treating the electron fraction for the innermost ejecta, which is a crucial ingredient for nucleosynthesis calculations. In addition, unlike the electron (anti-)neutrino luminosities, which decrease suddenly once the shock has been revived in spherically symmetric models, ν_x luminosities are only marginally affected by the development of an explosion. This allows PUSH to continue injecting energy inside the expanding shock for a few hundreds of milliseconds after the explosion has set in. A first implementation of the PUSH method and its calibration strategy was extensively documented in Paper I (Perego et al. 2015).

2.1. Simulation Framework of PUSH

The simulations presented in this paper have been performed with the same numerical setup as described in Perego et al. (2015). Hence, we only provide the essential aspects here and refer the reader to Paper I for more details. We use the general relativistic hydrodynamics code Agile, which uses an adaptive mesh to achieve a good resolution at the shock front and at the PNS surface (Liebendörfer et al. 2001). For the neutrino transport, we use IDSA for the electron flavor neutrinos (Liebendörfer et al. 2009) and ASL for the heavy flavor neutrinos (Perego et al. 2016). Dense and hot nuclear matter in NSE conditions is described by the HS(DD2) nuclear, finite temperature EOS (Hempel & Schaffner-Bielich 2010; Typel et al. 2010). For the non-NSE regime, we use an extension to ideal gas coupled with an approximative alpha-network (Hempel & Schaffner-Bielich 2010; Perego et al. 2015). The PUSH method relies on additional (parameterized) heating from heavy-lepton flavor neutrinos. It was introduced to take increased heating efficiencies based on multidimensional effects into account and it enables exploding models in spherical symmetry. The additional energy deposition is given by the local heating term

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

where

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

The term $({{dL}}_{{\nu }_{x}}/{dE})/(4\pi {r}^{2})$ is the spectral energy flux for any single ν_x neutrino species with energy E, σ₀ the typical neutrino cross-section, and ${m}_{b}\approx 1.674\times {10}^{-24}{\rm{g}}$ an average baryon mass. The function ${ \mathcal G }(t)$ (which contains the two free parameters ${k}_{\mathrm{push}}$ and ${t}_{\mathrm{rise}}$) controls the temporal evolution of PUSH, as shown in Figure 1.

Zoom In Zoom Out Reset image size

The spatial dependence of PUSH is described by

$$\begin{eqnarray}{ \mathcal F }(r,E)=\left\{\begin{array}{ll}0 & \mathrm{if}\quad r\gt {R}_{s}\quad \mathrm{or}\quad {\dot{e}}_{{\nu }_{e},{\overline{\nu }}_{e}}\lt 0\\ \exp (-{\tau }_{{\nu }_{e}}(r,E)) & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 3 }$$

where R_s is the shock radius, ${\dot{e}}_{{\nu }_{e},{\bar{\nu }}_{e}}$ the net specific energy deposition rate due to electron neutrinos and anti-neutrinos, and ${\tau }_{{\nu }_{e}}$ is the optical depth for electron neutrinos. Note that the entropy-gradient criterion from Paper I (where PUSH only provides extra heating in regions where ${ds}/{dr}\lt 0$) is no longer included in the spatial heating term. The effect of lifting the entropy-gradient criterion on the simulations and on our calibration procedure will be discussed in detail in Section 2.3.1.

We stress that we do not modify the electron neutrino and antineutrino transport (hence, preserving a consistent Y_e-evolution) and that the mass cut emerges naturally from the simulations. Moreover, we want to emphasize that we can self-consistently follow the simulation through to the formation of a BH (${\rho }_{{\rm{c}}}\gt {10}^{15}$ g cm⁻³) since we include the full PNS in our computational domain.

2.2. Initial Models

We use two sets of nonrotating pre-explosion models with solar metallicity from the stellar evolution code KEPLER presented in Woosley et al. (2002) and Woosley & Heger (2007). From here on, we will refer to these sets of pre-explosion models as WHW02 and WH07, respectively. The two sets of models span zero-age main sequence (ZAMS) masses between 10.8 M${}_{\odot }$ and 75 M_⊙ (WHW02) and between 12 M_⊙ and 120 M_⊙ (WH07). Table 1 summarizes all the pre-explosion models used, including their ZAMS mass and a unique label for each model.

Table 1.  Pre-explosion Models Used in This Study

Series	Label	Min Mass	Max Mass	Δm	References
		(M⊙)	(M⊙)	(M⊙)
WHW02	s	10.8	28.2	0.2	1
		29.0	40.0	1.0	1
		75.0			1
WH07	w	12.0	33.0	1.0	2
		35.0	60.0	5.0	2
		70.0			2
		80.0	120.0	20.0	2

Note. All models have solar metallicity. Note that the pre-explosion models are irregularly spaced in mass. The Δm indicates the mass spacings of available pre-explosion models for each mass interval Min Mass to Max Mass.

References. (1) Woosley et al. (2002); (2) Woosley & Heger (2007).

Download table as:  ASCIITypeset image

It is well-known that the initial stellar mass is a poor predictor of the final fate (successful or failed explosion) of massive stars. Hence, we make use of the compactness parameter (O'Connor & Ott 2011)

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

Compactness values evaluated for M between 1.5 and 3.0 M_⊙ are commonly used, see, e.g., O'Connor & Ott (2011, 2013), Nakamura et al. (2015), and Burrows et al. (2018). In our Paper I, where we investigated exploding models that were aimed at reproducing SN 1987A, we used M = 1.75 M_⊙. In this work, we investigate a much broader ZAMS-mass range of progenitors, which involves the formation of BHs and more massive BHs than in the previous study. To account for the larger possible mass included in the layers of the star that are crucial for the explosion mechanism as well as the larger possible NS masses and BH formation, we chose compactness values (at bounce) evaluated at 2 M_⊙. The compactness parameter for all WHW02 and WH07 models is shown in Figure 2 where solid lines denote the compactness at bounce and dashed lines at the onset of collapse. The difference in compactness between the two sets illustrates the uncertainties involved in the pre-explosion models themselves. The compactness curves reflect the structure of the pre-explosion models, as illustrated in Figures 3 and 4, which show the Fe-core mass (defined as the layers with Y_e < 0.495), the carbon–oxygen (CO) core mass (enclosed mass with X_He ≤ 0.2, i.e., up to the beginning of the He-shell), He-core mass (mass regions with X_H ≤ 0.2, i.e., up to the beginning of the H-shell), and the total mass including the H-envelope at the end of the stellar evolution simulation for the WHW02 and WH07 models. For models with ZAMS masses up to ∼35 M_⊙ (WHW02) and ∼40 M_⊙ (WH07), respectively, the CO-core mass grows continuously with increasing initial mass, while the envelope exhibits a different behavior: The helium envelope remains approximately constant up to 20–25 M_⊙ ZAMS mass. For higher ZAMS masses, mass loss decreases the envelope mass (hydrogen and helium) until at approximately 30–40 M_⊙ ZAMS mass, the models are almost completely stripped of their envelope. The decrease of the H/He-envelope mass is mirrored in the compactness. Beyond 30 M_⊙ the rise in compactness reflects the increasing CO-core mass, as these models have lost virtually all of their H- and He-envelope. For initial masses above 40 M_⊙, the evolution is quite uncertain, in particular, because mass loss is poorly understood. We include these models for completeness; however, we caution the reader that they may be quite uncertain. Moreover, successful explosions of models beyond 30 M_⊙ result in Type Ib or Ic SNe due to the loss of the H- (and He-) envelope. More generally, the efficiency of semiconvection and overshooting, together with the metallicity, leads to a rather complex interplay of different burning shells (Rauscher et al. 2002), which contributes to the rapid, nonmonotonic variations in compactness with ZAMS mass (Sukhbold & Woosley 2014; Sukhbold et al. 2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.3. Hydrodynamic Simulations

Using the framework described in Section 2.1 and the initial models from Section 2.2, we perform hydrodynamic simulations of the collapse, bounce, and post-bounce phases. We include the pre-explosion model up to the He-layer in the hydrodynamic simulations with the exception of the low-mass model s11.2, where the computational domain ends just below the He-layer and a few models with ZAMS masses above 30 M_⊙ with large pre-explosion mass loss. The total run time of the simulation depends on its outcome, which we categorize based on the explosion energy (representing the kinetic energy of the ejected matter) and the central density (indicating black hole formation if in excess of ∼10¹⁵g cm⁻³ and ultimately diverging). We compute the explosion energy as in Paper I: we assume that the total energy of the ejecta with rest-masses subtracted eventually converts into kinetic energy of the expanding supernova remnant at late times and we include the gravitational binding energy of the progenitor (see also Equations (12)–(17) in Paper I). For exploding models, we define the explosion time t_expl as the time when the shock reaches 500 km, measured with respect to core bounce. We distinguish between successful explosions, BH formation, and failed explosions, defined as follows: (i) If a simulation completes the set total simulation time of t_final = 5 s (approximately 4.6 s–4.8 s post-bounce) and/or has a saturated positive explosion energy it corresponds to an explosion. (ii) If a simulation forms a BH in the 5 s run time it corresponds to a nonexploding model (BH formation). (iii) If a simulation has a negative explosion energy at times when PUSH is no longer active (or at t_final) it corresponds to a failed explosion and hence a nonexploding model that will eventually form a BH.

2.3.1. Entropy-gradient Criterion

As indicated in Section 2.1, the current version of the PUSH heating term does not include the entropy-gradient criterion anymore. Unlike simulations in spherical symmetry, multidimensional simulations exhibit simultaneous downflow and outflow of matter. These regions are at similar distances from the PNS, but have quite different entropies (see, for example, Pan et al. 2016). This is a fundamentally multidimensional behavior and is difficult to reproduce in spherically symmetric simulations. The spatial heating term ${ \mathcal F }(r,E)$ in Paper I included a criterion (the entropy-gradient criterion) that restricted heating to regions where the one-dimensional system fulfills the condition for convective instability according to the Schwarzschild criterion (i.e., to zones with ${ds}/{dr}\lt 0$). This proved to be too restrictive in spherically symmetric simulations and did not allow for the extended heating seen in multidimensional systems. It introduces a self-canceling effect of the increased heating efficiency on a relatively short timescale. This severely limits the heating and considerably confines the allowed parameter space. This self-canceling effect of the heating due to the entropy-gradient criterion required relatively small values of ${t}_{\mathrm{rise}}$ in Paper I and hence led to a rapid evolution of the entropy profiles. When combined with the observational constraints from SN 1987A (explosion energy and Ni ejecta yields), this resulted in relatively early explosions. The explosion time impacts the location of the mass cut and hence the ejecta mass. The mass cut determines the total amount of Ni ejecta and also how much material with ${Y}_{e}\lt 0.5$ is ejected (and hence the relative amounts of ⁵⁷Ni and ⁵⁸Ni). For early explosions, the mass cut resides deeper in the pre-explosion structure while for late explosions the mass cut is further out in the star. With the entropy criterion on (Paper I), it was only possible to reproduce SN 1987A if we imposed 0.1 M_⊙ of fallback by hand. Generally, fallback consists of two components: the early fallback, which can be only determined in multidimensional simulations, and the late fallback, which requires simulation times far longer than what is feasible with the presented PUSH setup. The late fallback, due to shock reflection at density jumps of outer shell boundaries is not included in this study as its effect is considered to be small (Sukhbold et al. 2016). In this work, we relaxed the entropy-gradient criterion. In this case, a broader range of values for ${t}_{\mathrm{rise}}$ allows for explosion with energies around 1 Bethe. This can be seen in Figure 5, which shows explosion energy as a function of ${t}_{\mathrm{rise}}$ and ${k}_{\mathrm{push}}$ with (triangles) and without (circles) the entropy-gradient criterion for the s18.8 model. Points within the shaded region represent simulations that are consistent with observational constraints of explosion energy. Also note that larger values for ${t}_{\mathrm{rise}}$ result in lower explosion energies and in later explosions. For the largest value of ${t}_{\mathrm{rise}}$ considered here, only simulations without the entropy criterion result in explosion energies of around 1 Bethe. This leaves the choice between early explosions with the necessity to impose fallback or later explosions without the need for additional fallback. This choice only plays a very minor role for the outcome of the simulations (successful explosion with a NS or failed explosion with a BH), as in both cases the final ejecta mass is very similar (the difference in mass cut is similar to the amount of imposed fallback that is required to reproduce the observed nucleosynthesis yields). However, it is not obvious how much fallback should be imposed for models with progenitor masses outside of the range of SN 1987A. We find that the new setup of this paper (without the entropy-gradient criterion) enables explosions that are in better agreement with the temporal evolution of the shock radius, neutrino heating rates, and entropy profiles of multidimensional simulations by Pan et al. (2016), Bruenn et al. (2016), Ebinger (2017), and that have nucleosynthesis yields and explosion energy consistent with SN 1987A. In the following, all the results presented are for simulations using this new setup.

Zoom In Zoom Out Reset image size

2.4. Nucleosynthesis Post-processing

We use a post-processing approach to predict the detailed composition of the ejecta. The ejecta are divided into tracer shells of 10⁻³ M_⊙ each. We follow the trajectory of each mass element that reaches a peak temperature ≥1.75 GK using a nuclear reaction network that includes isotopes from free nucleons to ²¹¹At. For details, we refer to Paper III.

Within each tracer, we record the evolution of fundamental hydrodynamical and neutrino quantities, including radius, density, temperature, electron fraction, neutrino luminosities, and mean energies. For each mass element of the ejecta, we start the nucleosynthesis post-processing when the temperature drops below 10 GK (∼0.86 MeV), using the NSE abundances (determined by the current electron fraction) as the initial composition. For the tracer particles that never reach 10 GK, we start the post-processing at the beginning of the hydrodynamic simulation and use the abundances from the approximate α-network at this point (WHW02) or from the progenitor model (WH07) as the initial composition. For the times t < t_final, we use the hydrodynamic evolution. If the conditions at t_final are still sufficiently hot for nucleosynthesis to occur, we extrapolate radius, temperature, and density as follows:

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

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

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

where r is the radial position, v is the radial velocity, ρ is the density, T is the temperature, s is the entropy per baryon, and Y_e is the electron fraction of the tracer. The subscript "final" indicates the end time of the hydrodynamical simulation. Equations (5)–(7) correspond to a free expansion for the density and an adiabatic expansion for the temperature (e.g., Korobkin et al. 2012). For the extrapolation, we calculate the temperature at each time-step using the equation of state of Timmes & Swesty (2000).

In this section, we develop PUSH into a method that can be applied to a broad set of pre-explosion models to predict remnant properties (this paper) and detailed nucleosynthesis (Paper III). The direct application of the best-fit parameters from Paper I to the whole progenitor range may sound straight forward but we found that a more variable approach is necessary to reproduce the available observational properties of CCSNe. This complication does not come as a big surprise, since, after all, the understanding of the nature of the explosion mechanism of CCSNe has been elusive for decades.

3.1. Calibration Against SN 1987A

In Paper I, the calibration of PUSH was performed using solar metallicity models between 18 and 21 M_⊙ from Woosley et al. (2002). This resulted in two possible candidates, s18.0 and s19.4, both of which were able to reproduce the observed Ni yields when we imposed 0.1 M_⊙ of fallback by hand. In this paper, we use an updated spatial heating term that does not include the entropy criterion anymore (see Section 2.3.1). Therefore, we repeat the calibration procedure as described in Paper I using the pre-explosion models from the WHW02 set with ZAMS masses between 18.0 and 21.0 M_⊙ in order to find a suitable candidate able to reproduce the observed properties of SN 1987A. We find that the s18.8 model is in good agreement with the observed explosion energy and yields of Ni and Ti for k_push = 4.3 and t_rise = 400 ms (Table 2) without the need for any additional fallback. The observational estimates for SN 1987A are summarized in Table 3. These values are the same as those used in Paper I except for ⁴⁴Ti where we now use the newer results from Boggs et al. (2015) instead of Seitenzahl et al. (2014). For comparison, Table 3 also gives the corresponding values from our calibrated model for SN 1987A.

Table 2.  Parameter Values from Calibration

kpush (−)	trise (ms)	ton (ms)	toff (s)
4.3	400	80	1

Download table as:  ASCIITypeset image

Table 3.  Observed and Calculated Properties of SN 1987A

Quantity	SN 1987A	PUSH
	(observed)	(s18.8)
Eexpl $({10}^{51}\,\mathrm{erg})$	1.1 ± 0.3	1.2
Mprog (M${}_{\odot }$)	18–21	18.8
${}^{56}\mathrm{Ni}$ (M${}_{\odot }$)	(0.071 ± 0.003)	0.069
57Ni (M${}_{\odot }$)	(0.0041 ± 0.0018)	0.0027
58Ni (M${}_{\odot }$)	0.006	0.0066
44Ti (M${}_{\odot }$)	$(1.5\pm 0.3)\times {10}^{-4}$	3.05 × 10−5

Note. The s18.8 model was identified as the progenitor that reproduces SN 1987A for ${k}_{\mathrm{push}}=4.3$ and t_rise = 400 ms. The nucleosynthesis yields for SN 1987A are taken from Seitenzahl et al. (2014) except for ⁵⁸Ni, which is taken from Fransson & Kozma (2002), and ⁴⁴Ti, which reflects the value from Boggs et al. (2015). The explosion energy is adapted from Blinnikov et al. (2000). The corresponding values from our calibration model are shown for comparison.

Download table as:  ASCIITypeset image

3.2. Black Hole Formation

An important question in the investigation of CCSNe across the whole progenitor mass range is whether a collapsing star ultimately leads to a successful supernova explosion or whether it fails to explode and forms a BH. In this section, we turn our attention to failed SNe and the formation of BHs. To do this, we investigate the behavior of spherically symmetric CCSN models without the application of the extra heating from PUSH (i.e., setting k_push = 0). These simulations are not expected to explode. Rather, they will collapse to BHs. The timescale on which the different models undergo collapse can depend on the progenitor structure and on the choice of the EOS. We use a subset of models that samples the mass range up to 40 M_⊙ of both sets of pre-explosion models (WHW02 and WH07) and perform simulations for two equations of state (SFHO and HS(DD2); Hempel & Schaffner-Bielich 2010; Steiner et al. 2013; Fischer et al. 2014) until the formation of a central black hole, which starts to occur for central densities around ∼10¹⁵ g cm⁻³.

In Figure 6, the BH formation times for several pre-explosion models from WHW02 and WH07 are given. The differences in BH formation time between the progenitors can be related to different accretion rates, which are correlated to compactness (see Equation (4)). Figure 7 shows the temporal evolution of the central density of the two 40 M_⊙ models from the two pre-explosion model sets (WHW02 in blue; WH07 in green) and for two different equations of state (HS(DD2) solid lines, SFHO dashed lines). It is evident that the BH formation time strongly depends on the EOS used (indicated by the colored areas; models using the SFHO EOS collapse faster than the ones using the HS(DD2) EOS) and even stronger on the progenitor model (difference between green and blue lines). Thus, both aspects can have an impact on the explodability of numerical models since they set the estimated upper limit for the timescale on which a delayed mechanism should revive the stalled shock.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the collapse timescales (see Figure 6), we can draw first conclusions on the explodability of the different progenitors. As expected, progenitors with higher compactness have a shorter collapse time. In addition, the WH07 models collapse to BHs on shorter timescales overall. This relation between collapse time, or BH formation timescale, and compactness is also discussed in other works (see, e.g., O'Connor & Ott 2011). In our simulations, the BH candidates have compactness ξ_2.0 > 0.5. The BH formation timescales from simulations with k_push = 0 can be interpreted as upper bound for the timescale on which the neutrino-driven mechanism can yield successful explosions. From this, we expect that (at least some) pre-explosion models with compactness ξ_2.0 > 0.5 do not explode (i.e., form black holes), which can later be summarized as a constraint on the PUSH method when applied to an entire set of pre-explosion models.

3.3. Standard Neutrino-driven CCSNe

The PUSH method represents a parameterization of the neutrino-driven mechanism intended for the investigation of CCSNe in spherically symmetric simulations. Hence, we focus on standard neutrino-driven CCSNe (between ∼10 and 21 M_⊙ where 21 M_⊙ represents the upper mass limit of SN 1987A which we assume to be a standard neutrino-driven explosion) and the putative faint/failed SN branch (above 21 M_⊙). Note that a definitive consensus on the existence of such a "faint supernova branch" is still debated based on the explosion observations and progenitor mass determination (see Appendix A).

Based on general features of CCSN observations, we require that the PUSH method should

• (i)  
  reproduce the observed properties of SN 1987A for a suitable pre-explosion model,
• (ii)  
  allow for the possibility of BH formation, and
• (iii)  
  result in lower explosion energies for stars with masses of ≲13 M_⊙ ("Crab-like SNe")

when applied across the entire mass range.

In the following, we quantify the general requirements and discuss how they are included in setting the parameters of the PUSH method such that it ultimately can be applied to an entire mass range of pre-explosion models and is in agreement with general observational properties of CCSNe across the ZAMS-mass range.

As a first step, we applied the parameters from the calibration (see Table 2) to all WHW02 pre-explosion models. This approach leads to robust explosions for all considered stars if reasonable explosion energies of the order of 1 Bethe are supposed to be achieved for suitable candidates of SN 1987A. The resulting explosion energies cannot explain the observations and they even partially behave opposite to the expectations (see Section 3.2 and also Appendix A). With a constant ${k}_{\mathrm{push}}$ for the entire mass range ruled out, a compactness-dependent ${k}_{\mathrm{push}}$ factor is the next obvious choice for the calibration of our effective model such that it can fulfill all three requirements we imposed.

Any calibration of PUSH has to include the best fit for SN 1987A, i.e., have k_push = 4.3 for ξ_2.0 = 0.245 (model s18.8). The compactness values of progenitors for ZAMS masses below 21 M_⊙ are ξ_2.0 ≤ 0.4–0.5. To emulate a transition from the standard convective neutrino-driven mechanism to a regime of less efficient convective neutrino-driven mechanism and eventually to BH formation, we tune down the ${k}_{\mathrm{push}}$ parameter above ξ_2.0 = 0.4–0.5. This leads to the constraint for ${k}_{\mathrm{push}}$ to approach zero for compactness values above ξ_2.0 ∼ 0.5 and to be set to zero for compactness values ξ_2.0 ≥ 0.7. The observed explosion energies of progenitors at the lower end of the mass range indicate weaker explosions (and hence lower values of ${k}_{\mathrm{push}}$). These lower ZAMS-mass progenitors coincide with lower compactness values. Thus, we obtain the constraint for ${k}_{\mathrm{push}}$ to be smaller for lower compactness values. We choose k_push = 2.5 for ξ_2.0 = 0.0 as the third constraint that the compactness-dependent PUSH parameter function ${k}_{\mathrm{push}}$ (ξ) ultimately has to fulfill. In summary, the three points defining the dependence of ${k}_{\mathrm{push}}$ on the pre-explosion model compactness are:

• (i)  
  k_push = 2.5 at ${\xi }_{2.0}=0.0$ ("Crab-like SNe"),
• (ii)  
  ${k}_{\mathrm{push}}=4.3$ at ξ_2.0 = 0.245 (SN 1987A calibration model s18.8), and
• (iii)  
  k_push = 0.0 at ξ_2.0 ≥ 0.7 (BH formation).

We assume a polynomial dependence of k_push on the compactness and, given that we have three fixed points, we consider a parabolic dependence: ${k}_{\mathrm{push}}(\xi )=a{\xi }^{2}\,+\,b\xi +c$ (see also Figure 8). The resulting parabola through the three points given above is described by $a=-23.99$, b = 13.22, and c = 2.5.

Zoom In Zoom Out Reset image size

In addition to the parabola, we also tested piece-wise linear dependences of ${k}_{\mathrm{push}}$ on the compactness (using the same three fixed points as for the parabola) and comparable values for the calibration to SN 1987A. The results are very similar, confirming that our results are not sensitive to the functional dependence of ${k}_{\mathrm{push}}$ on the compactness.

Now, we apply this compactness-dependent ${k}_{\mathrm{push}}$ to both series of pre-explosion models. The left panel of Figure 9 shows the resulting landscapes of explodability for all the progenitors in the WHW02 and WH07 sets in comparison with observed explosion energies (see Appendix A and Bruenn et al. 2016 for the selection of observed SNe and a discussion of the associated uncertainties). For both sets, the explosion energy increases from "Crab-like SNe" at the lowest-mass end to robust explosions (with E_expl ≈ 0.8—1.6 Bethe) between 15 M_⊙ and 21 M_⊙. Above about 25 M_⊙, no observational data are available. We find that models between 22 M_⊙ and 26 M_⊙ form BHs (see the vertical dashes at the bottom of the figure in this mass range). Above 26 M_⊙, we find a mixture of successful explosions (models that have lost most of their hydrogen and helium envelopes) and BHs (models with higher CO-core and Fe-core masses). Due to the relative sparseness of pre-explosion models at higher ZAMS masses, it is not clear whether we obtain many faint explosions, or if we instead go back and forth between regular explosions and BH formation (cf. Figures 9 and 10). The models with ZAMS mass above 50 M_⊙ explode for both sets. When looking at the Ni ejecta masses for the same models together with the explosion energy (Figure 9, right panel), our models lie within the same range of explosion energies and Ni ejecta masses as the observations.

Zoom In Zoom Out Reset image size

We also investigated a somewhat different interpretation of the observational constraints by using the compactness ξ_1.75 instead of ξ_2.0.⁷ The compactness ξ_1.75 is more strongly dependent on the iron core mass and shows a slightly different behavior than the compactness ξ_2.0. The left panel of Figure 10 shows the parabolic dependence of ${k}_{\mathrm{push}}$ for the second calibration and the right panel of Figure 10 shows the resulting explosion energies compared to observational data for this second calibration. In Figure 11, we summarize the outcomes of both calibrations for both progenitor samples. It is evident that the second calibration results in a lower overall explodability and leads to a larger fraction of faint explosions and black holes. The first calibration, which represents our standard calibration, is in good agreement with observations (SN 1987A, putative faint SN branch, and BH formation) and is well suited to investigate the full range of CCSN progenitors. The relation between compactness and explodability (i.e., explosion energies and BH formation times; first investigated in O'Connor & Ott 2011), is expected to be one of the relevant aspects of the neutrino-driven mechanism. We showed that with a relatively simple calibration in compactness we can fulfill the constraints we formulated and obtain an effective model to investigate the properties of CCSNe. From here on, we discuss and use the standard calibration of PUSH introduced in this section unless otherwise indicated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we present and discuss the explosion properties of our simulations for the entire WHW02 and WH07 samples (combined including ZAMS masses of 10.8 – 120 M_⊙) based on the standard calibration presented in the previous section. These same models are used in the detailed discussion of the corresponding nucleosynthesis yields and their possible impact an the galactic chemical evolution in Paper III.

4.1. Explosion Properties

Figure 12 gives an overview of the resulting explosion and the predicted remnant properties for all pre-explosion models considered in this study. From top to bottom, we show the explosion energy, the explosion time, the ejected Ni mass, the total ejecta mass, and the remnant mass (baryonic mass) as a function of ZAMS mass for the WHW02 (left column) and WH07 (right column) pre-explosion models. Note that the WH07 series consists of a smaller number of models than the WHW02 series. The general features are similar for both series. With the standard calibration, we obtain explosion energies from less than 0.5 to 1.7 Bethe, gravitational NS masses from 1.2 to 1.8 M_⊙ (see also Section 5), and ejected nickel masses from 0.02 to 0.16 M_⊙. The lowest explosion energies and Ni masses are obtained for the lowest-mass progenitors (around 11–12 M_⊙), while the highest explosion energies and Ni masses result from pre-explosion models around 15 M_⊙, 18 M_⊙, and 21 M_⊙ ZAMS mass.

Zoom In Zoom Out Reset image size

For the most massive pre-explosion models (≳30 M_⊙), there is more variation both between the two series used in this work as well as in comparison with other works. The pre-SN evolution of such massive stars is much less certain (e.g., Hirschi 2015) and the outcome of core collapse is much more sensitive to the details of the prescription (see, for example, Figures 8 and 9 in Sukhbold et al. 2016). If these stars explode, they are of Type Ib or Ic. Overall, the WH07 series has lower explosion energies and is more prone to BH formation (as discussed in Section 3.2). Coupled to that, this series results in somewhat more massive NSs. The WH07 series also has ejecta masses that are slightly larger for the same ZAMS mass when compared to the WHW02 series. Note that for lower metallicity, massive stars experience less mass loss during their evolution and hence their mass at collapse is closer to their initial ZAMS mass. Above ∼30 M_⊙ they also have a higher compactness than their counterparts with solar metallicity. From this, we expect pre-explosion models at low metallicity to not have any successful explosion above ∼30–35 M_⊙, constituting a potential origin of the BHs seen by LIGO-VIRGO.

Our results for the solar metallicity pre-explosion models show a moderate trend of the explosion properties with ZAMS mass up to about 15 M_⊙ as expected from our calibration. This is the same mass range where the compactness exhibits an increasing trend with ZAMS mass (and where the CO-core mass grows with increasing ZAMS mass, see Figures 2–4). Beyond about 15M_⊙, there is no obvious trend with ZAMS mass. In this mass range, the relation between ZAMS mass and compactness is also more complex, as already discussed in Sukhbold & Woosley (2014). Not surprisingly, we see a correlation between the explosion energy, the ejected Ni mass, and to a lesser degree the remnant mass. The total ejected mass is dominated by the mass at collapse, so we see a similar trend of the total ejecta mass with ZAMS mass as seen in Figure 3.

Figure 13 shows explosion energy, remnant mass, and explosion time against compactness ξ_2.0. This reveals several trends. The strongest correlation is seen between the remnant mass (neutron stars only) and the compactness. More compact models have a higher mass accretion rate and hence more matter is accreted onto the PNS before the successful explosion (Perego et al. 2015). This trend to higher remnant mass would have continued to higher compactness if we had applied a constant value of ${k}_{\mathrm{push}}$ for all models, without forming any BHs. Instead, with our standard calibration for PUSH these models do not explode and eventually collapse to BHs, which are not included in this Figure. The explosion time shows a similar trend with compactness; however, the distribution is broader and the lowest-compactness models do not follow the general trend. We note that for the models of intermediate compactness (ξ_2.0 ≈ 0.3), where the PUSH heating reaches maximum values, the explosion time becomes comparable to the value of t_rise. The behavior of the explosion energy with compactness has two distinct features. The highest explosion energies for each model set show a parabolic dependence on the compactness, as expected from our ${k}_{\mathrm{push}}$ dependence. In addition, we observe a large scatter (up to 0.5 B) in the explosion energies for models with compactness between 0.2 and 0.45. In this intermediate compactness region, the explosion energies and the explosion times fall into two groups. In Figure 2, one can see a peak in compactness value around 24.2 M_⊙ (WHW02) and 23.0 M_⊙ (WH07). Models to the left of the peak (i.e., with lower ZAMS mass) and to the right of the peak (i.e., with higher ZAMS mass), have similar compactness values of 0.2 to 0.45. Nevertheless, they are different in their behavior. In Figure 13, the models to the left of the peak in compactness are indicated by open markers and the models to the right of the peak in compactness by filled markers. Lower ZAMS-mass stars before the peak in compactness ultimately take longer to form explosions and have a higher explosion energy than higher ZAMS-mass stars of the same compactness but different ZAMs and CO-core masses. For the remnant mass, this split in two groups is less pronounced.

Zoom In Zoom Out Reset image size

Similar studies of CCSN properties have been done in Ugliano et al. (2012), Ertl et al. (2016), and Sukhbold et al. (2016) based on P-HOTB. By calibrating their spherically symmetric models to SN 1987A and, in the case of Ertl et al. (2016) and Sukhbold et al. (2016), to lower explosion energies for less massive progenitors, the explosion properties and yields of CCSNe were investigated. Another recent study on the explodability of CCSNe without the use of hydrodynamic simulations has been presented in Müller et al. (2016). This work uses physically motivated scaling laws and differential equations to describe crucial quantities of CCSNe like the shock propagation, the neutron star contraction and the heating conditions. Note that with the exception of Ugliano et al. (2012) these other works use different progenitor samples and range from 9 to 120 M_⊙, i.e., the smallest mass range is covered in Müller et al. (2016), which investigates progenitors between 10 and 32.5 M_⊙ ZAMS mass. Our explosion energies are compatible with values found in observations and similar to other studies (Ugliano et al. 2012; Pejcha & Thompson 2015; Ertl et al. 2016; Müller et al. 2016; Sukhbold et al. 2016). Exceptions are the relatively high explosion energies that Ugliano et al. (2012) and Pejcha & Thompson (2015) obtain for low-mass progenitors. In Sukhbold et al. (2016), slightly lower explosion energies and nickel ejecta masses are obtained with a calibration for "Crab-like" SNe for ZAMS masses below 12 M_⊙. The results presented in Pejcha & Thompson (2015) and Müller et al. (2016) yield the largest range of explosion energies, spanning from  0.2 to  6 B, and from a few 0.01 B to above 2 B, respectively. A prominent feature is the region of nonexplodability in the vicinity of 25 M_⊙ ZAMS mass. In our simulations, we find BH forming models in similar mass regions to other works at 20 M_⊙, around 25 M_⊙ and between 30 and 100 M_⊙. Similar to Müller et al. (2016), we do not find many BH forming models between 15 and 20 M_⊙. We want to stress that in other studies (Ertl et al. 2016; Müller et al. 2016; Sukhbold et al. 2016) different progenitors have been used and thus a perfect agreement is not expected.

4.2. Key Isotopes

We now turn our attention to the behavior of four key isotopes (⁵⁶Ni, ⁵⁷Ni, ⁵⁸Ni, and ⁴⁴Ti) with compactness. Figure 14 shows the ejected mass of ⁵⁶Ni (top left), ⁵⁷Ni (top right), ⁵⁸Ni (bottom left), and ⁴⁴Ti (bottom right) for all models in both series (WHW02 with blue circles; WH07 with green stars; s18.8 with a red triangle) as a function of compactness. The amount of ejected ⁵⁶Ni shows a strong linear correlation with the compactness. A similar correlation with compactness can be seen for ⁴⁴Ti. Both of these isotopes are symmetric, N = Z isotopes. For these isotopes, ⁵⁶Ni in particular, the amount ejected depends on the explosion energy and how much of the Si-shell gets processed and ejected. This trend is not shared by ⁵⁷Ni and ⁵⁸Ni. Instead, the correlation broadens toward high compactness (at ξ_2.0 there is almost a factor of three difference between the lowest and highest value of ejected ⁵⁷Ni). A similar, but stronger, trend is seen for ⁵⁸Ni. While the amount of ⁵⁶Ni and ⁴⁴Ti mainly depends on the explosion energy, the yields of ⁵⁷Ni and ⁵⁸Ni depend strongly on the local electron fraction and whether regions with slightly lower Y_e are ejected or not. In general, for models of similar compactness, the models with lower ZAMS mass (open symbols in Figures 13 and 14) eject more ^57,58Ni and ⁴⁴Ti. An extended discussion of the nucleosynthesis yields and the trend with the local Y_e of the exploding models presented here can be found in Paper III.

Zoom In Zoom Out Reset image size

Besides explosion energies and ejected nickel masses, there are additional observables against which we can compare our results, for example, the mass of the compact remnant. Observationally, the distribution of NS masses in slow pulsars and nonrecycled high-mass eclipsing binaries has a mean of 1.28 M_⊙ (dispersion of 0.24 M_⊙), which is consistent with the formation in CCSNe (Özel et al. 2012). With our simulation setup (no excised mass and realistic nuclear equation of state) we can follow the full evolution of the PNS and obtain the baryonic mass of the resulting newly born hot neutron star. The gravitational mass of the corresponding cooled NS (that may be observed) is smaller due to neutrino losses in the cooling phase. We compute the zero-temperature gravitational mass of the hot neutron star formed in our simulations using the HS(DD2) nuclear equation of state. We arrive at a distribution of gravitational birth masses of NSs by weighting the predicted NS masses as a function of ZAMS mass with the initial mass function (IMF) of stars according to Salpeter (1955), which is sufficiently accurate for stars with masses above 10 M_⊙. In Figure 15, we show the gravitational birth mass distribution of cool NSs for the two sets of pre-explosion models from our standard calibration. The different colors indicate different ranges of ZAMS masses of the pre-explosion models. The two samples combined include stars with masses between 10.8 and 120 M_⊙ (see also Table 1). The resulting NS masses are between 1.2 and 1.8 M_⊙. The lowest NS masses between 1.2 and 1.4 M_⊙ are only noteworthily populated by the WHW02 sample. In Figure 15, we see that lower ZAMS-mass stars contribute most of the lower NS masses around 1.4 M_⊙ and the higher ZAMS-mass stars are the main contribution to the NS masses in the vicinity of 1.6 M_⊙. The resulting distribution of NS birth masses for the second (less energetic) calibration does not differ too much from our standard calibration, as can be seen in Figure 15. In comparison to the reported observed NS distribution, the mass range of our predicted distribution is somewhat shifted to higher masses. Note that due to a lack of progenitor models below 10.8 M${}_{\odot }$, remnants from even lighter progenitors are not included, which due to the applied Salpeter-IMF also would considerably contribute to the lightest NSs and would arguably reduce the lower limit of the predicted NS mass distribution. Also note that for the WHW02 sample we did not include the outlier with 75 M_⊙ ZAMS mass. The results are mass distributions of single-star systems and do not consider the possible effects present in binary systems, such as mass accretion or loss due to a companion, that can lead to a very different evolution of the stars. Therefore, a statistically significant comparison of our models with the observed NS masses of double systems is beyond the scope of this paper. Overall, our NS masses are 0.1 to 0.2 M_⊙ higher than those in Ugliano et al. (2012) and in Pejcha & Thompson (2015) who used the same WHW02 pre-explosion models. Our results are in general agreement with Ertl et al. (2016), who use different pre-explosion models below 30 M_⊙. Nakamura et al. (2015) performed 2D axisymmetric simulations up to 1 s post-bounce of the same WHW02 pre-explosion models. Similar to us, they find NS masses in the range of 1.2 to 2.1 M_⊙. It should be noted that none of their simulations failed to explode on the timescale of their simulations. Since we do not take into account fallback, which has been found to not have a major effect in other studies (Sukhbold et al. 2016), all our exploding models result in NSs. This is a consequence of the 1D treatment, where it is not possible to have infalling and outgoing matter at the same time. Thus, we do not obtain any fallback for exploding models in our simulations, unlike in multi-D simulations where it is possible to have simultaneous in- and outflow, which allows for substantial fallback in 2D and 3D simulations (e.g., Marek & Janka 2009; Chan et al. 2018; Ott et al. 2018). In our framework, simulations of stellar collapse that run beyond the time in which PUSH is active and ultimately fail to explode and simulations that directly form a black hole contribute to the resulting birth mass distribution of BHs. Our failed explosions for the (nonrotating) WHW02 and WH07 pre-explosion models correspond to BH masses from failed neutrino-driven CCSNe of nonrotating (or weakly rotating) stars. The final mass that collapses to a BH depends on the amount of mass stripping and hence for solar metallicity progenitors can be as low as the CO-core mass. For the two samples considered here, the final stellar mass at the onset of collapse does not exceed ∼17 M_⊙ due to wind mass loss (see also Figures 3 and 4). Again, we weight the resulting BH masses with an IMF. Figure 16 shows the resulting BH mass distributions for our standard calibration. The different shaded regions correspond to different stellar cores that collapsed to a BH, to illustrate how mass loss may effect the final BH masses. For the WHW02 and WH07 progenitor samples, we find BH formation for stars between 20 and 30 M_⊙, centered around 25 M_⊙ ZAMS mass, resulting in BH masses centered around ∼14 M_⊙. Stars above 30 M_⊙ that form BHs are mainly found in the WH07 sample, which overall consists of models that are more likely to collapse to BHs. For our second calibration of PUSH, considerably more BHs are found, which shifts the resulting BH mass distribution to slightly lower mass (see Figure 16). The gap between possible low and high mass BHs is not as strongly present as before. For both calibrations we can compute the fraction of stars that ultimately form BHs. We do so by considering a mass range from 8 to 150 M_⊙ for the estimate and assuming that stars between 8 M_⊙ and the lowest ZAMS mass in each sample successfully explode and leave behind an NS as a remnant. Furthermore, the fate of the star with the highest ZAMS mass within each sample is continued up to 150 M_⊙. Again, we use the Salpeter IMF for our estimate. For the WHW02 progenitor sample, ∼5% of progenitor stars leave behind a BH for the standard calibration and ∼16% for the second calibration. The WH07 progenitor sample leads to more BHs: ∼8% of stars have a BH as a remnant for the standard calibration and ∼21% for the second calibration. Our results are broadly consistent with the observationally determined BH mass distribution ((7.8 ± 1.2) M_⊙, Özel et al. 2010), when we assume that the helium core mass sets the BH mass (Kochanek 2014). The BH masses found from the solar metallicity pre-explosion models (WHW02 and WH07) are not massive enough to explain the BH from the recent LIGO/VIRGO observations. These BHs can originate from low-metallicity stars, which experience less (or no) mass loss during their evolution and hence collapse with almost their entire ZAMS mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The investigation of the CCSN mechanism remains an intriguing and unsolved problem. The solution most likely requires multidimensional hydrodynamical simulations, including general relativity, a nuclear equation of state, sophisticated neutrino transport, magnetic fields and rotation, and asphericity of the progenitor structure. The high computational cost of such models still limits the number of models that can be investigated and still motivates the usage of effective spherically symmetric models for extended progenitor and outcome studies.

In this work, we have improved the PUSH method from Perego et al. (2015), which is used to artificially trigger parameterized core-collapse supernova explosions in spherical symmetry. The PUSH method provides a computationally affordable and robust framework to study aspects of CCSNe that require modeling many different progenitors for several seconds after the onset of the explosion. Here, we focused on applying it to two sets of nonrotating, solar metallicity pre-explosion models between 10.8 M_⊙ and 120 M_⊙ to obtain explosions and predict compact remnant properties.

The PUSH method has two free parameters (${k}_{\mathrm{push}}$ and ${t}_{\mathrm{rise}}$) that need to be determined from external constraints. We have used three general constraints from observations of core-collapse supernovae to set the parameters for any pre-explosion model a priori. We required that the PUSH method (i) reproduces the observed properties of SN 1987A for a suitable pre-explosion model, (ii) allows for the formation of black holes, and (iii) results in lower explosion energies for the lowest-mass progenitors ("Crab-like SNe"). These requirements led to a compactness-dependent value for ${k}_{\mathrm{push}}$ (and a fixed value for ${t}_{\mathrm{rise}}$).

Using this setup, we have simulated the death of 133 SN-progenitor models as either successful neutrino-driven explosion or as failed explosion or direct collapse to a black hole. This study has led to several interesting predictions and conclusions:

• 1.  
  As a whole, the resulting explosion energies from the PUSH method are in good agreement with the explosion energies of observed CCSNe. In addition, the same models also match observations simultaneously for ⁵⁶Ni ejecta and explosion energy.
• 2.  
  We have shown and discussed that it is possible to infer several interesting trends of explosion properties with compactness (and CO-core mass). The compactness is a better indicator for the expected outcome than the ZAMS mass; however, it does not tell the entire story. For example, we found that for models of similar compactness a degeneracy exists that can be partially broken with the CO-core mass. This is consistent with the finding in recent multidimensional simulations, where the outcome of the CCSN simulation is related to the binding energy of the outer envelope exterior of a given mass, which has a strong correlation with compactness (Burrows et al. 2018).
• 3.  
  We found linear trends of ejected ⁵⁶Ni and ⁴⁴Ti yields with compactness (and explosion energy). The yields of ^57,58Ni do not follow the same simple correlation. Instead, the local electron fraction has a strong impact on the final yields of these isotopes.
• 4.  
  The predicted outcome (neutron star or black hole) is in agreement with predictions from other works that employ a comparable parameterized approach based on the neutrino-driven CCSN mechanism. In particular, we also find a region of nonexplodability around ∼25 M_⊙ initial stellar mass. Note that some of these studies use different progenitor sets below 30 M_⊙.
• 5.  
  We predict NS mass and BH mass distributions that are broadly consistent with observations. However, we do not find any BHs with m_BH ≳ 18  M_⊙ since none of the pre-explosion models have a final mass at collapse above this due to mass loss. We expect that lower metallicity pre-explosion models (which experience less mass loss) will result in more massive BHs.
• 6.  
  This set of CCSN models is used for a detailed nucleosynthesis study in Paper III (Curtis et al. 2018).

The work at NC State (K.E., S.C., and C.F.) was supported by the Department of Energy through an Early CAREER Award (DOE grant No.SC0010263). C.F. acknowledges support from the Research Corporation for Science Advancement through a Cottrell Scholar Award. The effort at the Universität Basel was supported by the Schweizerischer Nationalfond and by the ERC Advanced Grant "FISH." K.E. acknowledges the support from GSI. A.P. acknowledges support from the INFN project "High Performance data Network" funded by the Italian CIPE.

Software: Agile Liebendörfer et al. (2001), CFNET Fröhlich et al. (2006a), Matplotlib (Hunter 2007).

Besides the well-observed SN 1987A, many other CCSNe are known with a range of progenitor masses and explosion energies. One sees a rising trend in explosion energy and ejected nickel mass from ∼10 to ∼20 M_⊙, which corresponds to the transition from weaker "Crab-like" SNe to standard neutrino-driven SNe. Beyond 20 M_⊙ ZAMS mass, the properties of observed CCSNe are uncertain (see also Section 1). Nomoto et al. (2003) argue that they can be separated into two branches ("hypernova (HN) branch" and "faint SN branch") with different magnitudes of explosion energy and ejected nickel mass, and arguably different mechanisms powering the explosions. We refer the reader, for example, to Figure 2 in Nomoto et al. (2003; see also Janka 2012; Nomoto et al. 2013) for the bifurcation into two branches. Note that a definitive consensus on the existence of such a "faint supernova branch" has not yet been reached by observers.

The HN branch consists of very energetic γ-ray burst SNe (GRB-SNe) and of hypernovae (HNe; Woosley & Bloom 2006; Janka 2012; Nomoto et al. 2013). For these explosions, rapid stellar rotation and strong magnetic fields are thought to be crucial. The term HNe originates from the exceptional brightness caused by a large production of nickel in these hyperenergetic explosions with explosion energies of ≳10⁵² erg (Paczyński 1998; Nomoto et al. 2004, 2006). A possible central engine of these HNe and GRB-SNe are rapidly rotating BHs (collapsars), where accreting matter around the central compact object radiates energy in neutrinos, electromagnetic Poynting flux, and mass outflow. An alternative scenario for the central engine of HNe and GRB-SNe is a rapidly spinning neutron star with a strong magnetic field (B ≳  10¹⁵ Gauss) that is formed during stellar collapse. In this case, the HN or GRB-SN is powered by rotational energy, which is converted into explosion energy by the magnetic field. The resulting explosion strengths and amount of Ni produced are higher than what is typically expected for the neutrino-driven mechanism (Nomoto et al. 2006). This is another indication that the HN branch is likely caused by a different mechanism, distinct from the standard neutrino-driven mechanism.

The existence of the other branch, the faint SN branch, is observationally less certain. It represents CCSNe with low explosion energies. These SNe mark a transitional region to failed SNe, which indicates a mechanism that no longer can efficiently power explosions for stars above a certain ZAMS mass (the value is quite uncertain, with observations indicating 18 M_⊙ (Smartt 2015) to 25 M_⊙ (Davies & Beasor 2018)). These failed explosions lead to the formation of a BH either directly without an associated explosion or through the continued fallback of material onto the central object in very weak (and eventually failing) explosions (see, e.g., Lovegrove & Woosley 2013). The change in typical explosion energies within the faint SN branch could indicate a transition from efficient, strongly convective neutrino-driven SNe to an inefficient neutrino-driven mechanism. Thus, the observation of two branches (faint SN branch and HN branch), together with the constraint that HNe and GRB-SNe are likely due to a different and somewhat rarer mechanism, suggests that the faint SN branch represents neutrino-powered SNe.

In Table 4, we list the CCSNe included in our Figures 9 and 10 together with their relevant observational properties. This selection follows Bruenn et al. (2016) and Nomoto et al. (2013; only SNe with explosion energies ≤10⁵² erg). These CCSNe are in the mass range between 9 M_⊙ and 30 M_⊙ and can be considered to be caused by the neutrino-driven mechanism. However, it is important to emphasize that the observationally derived ZAMS masses of the CCSN progenitors have relatively large uncertainties depending on the method used. The first line for each SN indicates the values used in this paper (taken from Bruenn et al. 2016 or Nomoto et al. 2013 as indicated in the References column). Additional lines indicate that significantly different values are available in the literature. We have added these for the reader's convenience. However, a detailed discussion of the observational uncertainties is beyond the scope of this paper.

Table 4.  Observational Properties of CCSNe

Event	MZAMS	Eexpl	m(56Ni)		References
	(M⊙)	(1051 erg)	(M⊙)
			BL; S; L; O/Va
SN 1987A	18–21	1.1 ± 0.3	0.071 ± 0.003	L	(1), (2)
SN 1993J	12–17	1–2	0.06–0.09	L	(3), (4), (5), (6), (7)
			0.07–0.11b	L	(8)
SN 2004A	12.0 ± 2.1	0.76–1.3	${0.046}_{-0.017}^{+0.031}$	O/V	(3), (9), (10), (11), (12)
SN 2004dj	12–15	0.7–0.9	0.020 ± 0.002	O/V	(3), (9), (10), (13), (14), (15), (16), (17)
SN 2004et	12–15	1.1–1.8	0.062 ± 0.02	S	(3), (9), (10), (18), (19)
	25–29b	2.0–2.6b			(20)
			0.062 ± 0.02b	BL	(18)
	${8}_{-1}^{+5}$ b				(21)
SN 2005cs	${9}_{-2}^{+3}$	0.27–0.39	0.009 ± 0.003	BL	(3), (9), (10), (22), (23)
	17.2–19.2b	0.41 ± 0.03b			(24)
SN 2009kr	11–20	1.6–3			(3), (9), (10), (25)
SN 2012aw	14–18	1.0–1.7	0.06 ± 0.01	O/V	(3), (9), (10), (18), (26), (27)
		1–2b	0.06 ± 0.01b		(27)
		1.5b	0.05 ± 0.06b		(28)
SN 2012ec	14–22	0.6–1.9	⋯		(3), (9), (10), (29)
SN 1994I	∼13–15	∼1.0	∼0.07	O/V	(30), (31)
SN 2005bf	∼25–30	∼1.0–1–5	∼0.32	O/V	(30), (32), (33)
	8.3b	∼2b			(34)

Notes. With the exception of the values for SN 1987A, this table mainly adapts values presented in Bruenn et al. (2016) and Nomoto et al. (2013), and references therein.

^aThe abbreviations given in the column of the reported nickel values denote values obtained with the BL-method (BL), values obtained with the S-method (S), values obtained with modeling of the light curve (L), and values obtained by other methods or with a combination of various methods (O/V). For the event SN 1994I the values given here are taken from Nomoto et al. (2013, 1994). See Bruenn et al. (2016) for a discussion on the selection procedure. ^bValues in rows following a given SN name represent alternative ZAMS mass, explosion energy, and the amount of ejected nickel determinations for the readers convenience.

References. (1) Blinnikov et al. (2000); (2) Seitenzahl et al. (2014); (3) Bruenn et al. (2016); (4) Shigeyama et al. (1994); (5) Woosley et al. (1994); (6) Young et al. (1995); (7) Bartunov et al. (1994); (8) Freedman et al. (1994); (9) Poznanski (2013); (10) Dessart et al. (2010); (11) Hendry et al. (2006); (12) Maund et al. (2014); (13) Chugai et al. (2005); (14) Zhang et al. (2006); (15) Maíz-Apellániz et al. (2004); (16) Wang et al. (2005); (17) Vinkó et al. (2009); (18) Jerkstrand et al. (2012); (19) Sahu et al. (2006); (20) Utrobin & Chugai (2009); (21) Crockett et al. (2011); (22) Maund et al. (2005); (23) Takáts & Vinkó (2006); (24) Utrobin & Chugai (2008); (25) Fraser et al. (2010); (26) Jerkstrand et al. (2014); (27) Bose et al. (2013); (28) Dall'Ora et al. (2014); (29) Maund et al. (2013); (30) Nomoto et al. (2013); (31) Nomoto et al. (1994); (32) Janka (2012); (33) Tominaga et al. (2005); (34) Folatelli et al. (2006).

Download table as:  ASCIITypeset image

The explosion properties obtained from the standard PUSH calibration for the solar metallicity progenitor samples WHW02 and WH07 (see Figure 12) are provided as machine-readable tables. Only the exploding models are included in the table. A portion of the table is shown here for guidance regarding its form and content (see Table 5). The table is published in its entirety in machine-readable format.

Table 5.  WHW02 and WH07 Models (Standard Calibration)

MZAMS	Mcollapse	Eexpl	texpl	M${}_{{}^{56}\mathrm{Ni}}$	M${}_{\mathrm{remn}}$
(M${}_{\odot }$)	(M${}_{\odot }$)	(B)	(s)	(M⊙)	(M⊙)
18.0	14.50	1.45	0.44	1.12E–01	1.62
18.2	14.59	1.24	0.39	7.98E–02	1.55
18.4	14.83	1.58	0.42	9.72E–02	1.78
18.6	14.86	1.22	0.36	6.88E–02	1.56
18.8	15.05	1.20	0.36	6.90E–02	1.55
19.0	15.04	1.58	0.43	1.12E–01	1.78

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image