The Explosion of Helium Stars Evolved with Mass Loss

The code P-HOTB is a 1D, 2D, or 3D neutrino-hydrodynamics code for simulating the explosion of supernovae. The original version was derived from the Prometheus hydrodynamics code (Fryxell et al. 1989, 1991; Müller et al. 1991). It was first developed and employed in 1D and 2D simulations by Janka & Müller (1996) but later upgraded to include more microphysics and a more generalized numerical grid (e.g., Kifonidis et al. 2003; Scheck et al. 2006; Arcones et al. 2007; Ugliano et al. 2012; Ertl et al. 2016a). Since 2010, P-HOTB has also been applied to 3D supernova simulations by Hammer et al. (2010) and—supplemented for an axis-free Yin-Yang grid—by Wongwathanarat et al. (2010a, 2010b, 2013, 2015, 2017) and Gessner & Janka (2018).

The code P-HOTB is optimized to efficiently simulate neutrino-driven supernovae continuously from the onset of stellar core collapse, through core bounce and the onset of explosion, and onto the late stages of the ejecta expansion. In order to cover this large range of spatial and temporal scales for large sets of models, the treatment of the problem and the relevant microphysics is approximated in many aspects, but efforts are made to preserve the essential physics of the neutrino-driven mechanism.

The chief characteristics of P-HOTB are its use of an approximate, gray treatment of neutrino transport in the outer layers (optical depths below about ${ \mathcal O }(1000)$) of the newborn neutron star and the replacement of the central, high-density core of the neutron star by an inner grid boundary that shrinks with time. The movement of the grid boundary mimics the contraction of the cooling and deleptonizing compact remnant. Excising the core of the neutron star allows larger time steps, facilitating the efficient long-time simulation of the explosion including the central neutrino source. It also permits a choice of neutrino luminosities and mean spectral energies as inner boundary conditions that enable neutrino-driven explosions to be calculated with a chosen energy in spherical symmetry or higher dimensions. The replacement of the inner core of the nascent neutron star is justified by the fact that the high-density nuclear equation of state (EOS) remains uncertain despite considerable theoretical and experimental progress (see, e.g., Schneider et al. 2019).

At the onset of a calculation, the rest-mass density ρ, radial fluid velocity v, electron fraction Y_e, nuclear composition (in regions where nuclear statistical equilibrium (NSE) does not apply), and pressure P as functions of radius from the progenitor data are mapped onto the computational grid of P-HOTB. Using the pressure guarantees that hydrostatic equilibrium in the progenitor star and the developing core contraction at the onset of core collapse are well represented. The period from core collapse until core bounce is treated using the deleptonization scheme of Liebendörfer (2005), the Y_e(ρ) trajectory of Figure 1 in Ertl et al. (2016a), and the nuclear EOS of Lattimer & Swesty (1991) with an incompressibility modulus of 220 MeV. After core bounce, neutrinos and antineutrinos of all flavors are followed with the gray transport approximation.

About 10 ms after bounce (specifically, at the time the shock encloses 1.25 M_⊙), the central high-density core of the proto–neutron star (PNS) with a mass of 1.1 M_⊙ is excluded from the computational domain and replaced by a gravitating point mass. The gravitational force of this point mass and the self-gravity of the matter on the grid are corrected for general relativistic effects (see Arcones et al. 2007). For the hydrodynamics, the inner grid boundary corresponds to a Lagrangian mass coordinate at which the boundary conditions of the hydrodynamic variables ensure hydrostatic equilibrium. As long as the evolution of the PNS is tracked, the simulations of the present paper follow Arcones et al. (2007), Ugliano et al. (2012), Ertl et al. (2016a), and Sukhbold et al. (2016) in the contraction of the inner grid boundary and the boundary conditions for the neutrino properties at this location. The inner grid boundary, and with it the whole computational grid, is moved according to an exponential interpolation between initial radius R_ib,i and defined final radius R_ib,f: R_ib = R_ib,f + (R_ib,i − R_ib,f) · exp(−t/t₀). Here R_ib,i is the radius corresponding to an enclosed mass of 1.1 M_⊙ at about 10 ms after bounce. It is typically around 60 km. For the contraction timescale, we use t₀ = 0.4 s, and for the final radius, R_ib,f = 20 km (see also Ugliano et al. 2012).

The neutrino emission of the excised core, which defines the boundary condition for the neutrino transport at the inner grid boundary, is computed from an analytic one-zone model (see Ugliano et al. 2012; Sukhbold et al. 2016). The model equations are based on the requirements of energy conservation and the validity of the virial theorem. The relative weights of neutrinos and antineutrinos are constrained by the lepton number loss associated with the neutronization of the PNS core (Scheck et al. 2006). The values of the parameters in this analytic model are thus constrained by (global) physics arguments, but three of them are varied within small ranges to obtain explosions in 1D that are in agreement with observationally constrained energies and ⁵⁶Ni masses of the well-studied SN 1987A and SN 1054-Crab. These two observational cases are used to calibrate our "neutrino engine" that drives explosions in spherical symmetry. This choice is motivated, on the one hand, by the fact that the Crab supernova had a low explosion energy (only around 10⁵⁰ erg) with little nickel production (some 10⁻³ M_⊙), and the gaseous remnant has a high helium abundance and relatively oxygen-poor filaments. These characteristics are compatible with the explosion of a low-mass star with an extremely low core compactness value. On the other hand, SN 1987A had an explosion energy over 10⁵¹ erg and expelled about 0.07 M_⊙ of ⁵⁶Ni and roughly 1 M_⊙ of oxygen, classifying it as a more "normal" supernova of a massive star greater than 15 M_⊙. For a detailed discussion, see Sukhbold et al. (2016).

Following Sukhbold et al. (2016), we use a 9.6 M_⊙ star (model Z9.6) for a Crab progenitor and a set of 15, 18, and 20 M_⊙ blue supergiant stars (models W15, W18, W20, and N20), as well as a 19.8 M_⊙ red supergiant star (model S19.8), as SN 1987A progenitors. The values of the PNS core parameters for these neutrino engines, calibrated by comparison to Crab and SN 1987A, are listed in Table 3 of Sukhbold et al. (2016). In addition, we also tested several of the binary progenitor models proposed for SN 1987A by Menon & Heger (2017) and found that all of them could be exploded with at least one of our sets of engine-parameter values from single-star progenitors to yield reasonable agreement with the explosion energy and nickel production of SN 1987A (see Appendix A). Therefore, the binary progenitors of SN 1987A do not provide any new neutrino engines but can be considered as variants of the already available single-star engines.

In the present paper, we will present results of all engines used by Sukhbold et al. (2016), but for a detailed analysis of population-integrated supernova outcomes, we will mostly focus on Z9.6 in combination with W18. The latter provides a neutrino engine not much below the strongest ones (N20 and S19.8), and its choice allows for detailed comparisons with single-star results discussed by Sukhbold et al. (2016). In practice, we interpolate the PNS core model parameters between the Z9.6 and SN 1987A engines in terms of a parameter M₃₀₀₀, which measures the mass enclosed by a radius of 3000 km in the core of the precollapse star. This interpolation provides a smooth transition from "Crab-like" to "SN 1987A–like" cases. For details, see again Sukhbold et al. (2016).

In all simulations, the EOS of Lattimer & Swesty (1991) with K = 220 MeV is used for ρ ≥ 10¹¹ g cm⁻³. At lower densities, a multicomponent EOS is employed that includes radiation, electrons, positrons, and ions. The composition is followed using two approximations. Above a temperature of 7 × 10⁹ K, NSE is assumed, and the composition is tracked using a table that includes, as a function of ρ, T, Y_e, the 13 α-nuclei, neutrons, protons, and a "tracer nucleus" (Tr) that represents neutron-rich nuclei and replaces ⁵⁶Ni in regions where Y_e < 0.49 and neutron-rich species are expected to become abundant. Neutrino capture and electron and positron capture on neutrons and protons are followed so that Y_e evolves with time (for the neutrino treatment, see Scheck et al. 2006). For temperatures between 7 × 10⁹ and 10⁸ K, the composition is followed using a reaction network that links the 13 α-nuclei by the triple-α reaction for helium, subsequent (α, γ) reactions for nuclei with atomic mass numbers, A < 28, the faster (α, p)(p, γ) reaction channel compared to the latter for A ≥ 28, and heavy-ion reactions, ¹²C(¹²C, α)²⁰Ne, ¹²C(¹⁶O, α)²⁴Mg, ¹⁶O(¹⁶O, α)²⁸Si (original version by Müller 1986, upgraded by Kifonidis et al. 2003). In regions where lepton capture on nucleons has reduced Y_e below 0.49, the tracer nucleus Tr replaces ⁵⁶Ni in the network and is produced instead of it by the reaction ⁵²Fe(α, p)(p, γ)⁵⁶Ni. The small network approximately tracks the explosive burning of carbon, oxygen, and silicon as the shock moves out. Material that has been above 7 × 10⁹ K and is ejected has its composition followed in the small network as it expands and cools down, so that, e.g., the helium has an opportunity to recombine. In regions where silicon has been depleted, i.e., regions that have attained NSE, the final composition consists of a mixture of helium from photodisintegration, ⁵⁶Ni from regions where Y_e has remained greater 0.49 at all times, and Tr, a partially neutronized species whose actual identity reflects the history of weak interactions in the matter but is not well determined in P-HOTB.

As in Ertl et al. (2016a) and Sukhbold et al. (2016), the network and EOS are fully self-consistently coupled, thus taking into account the energy release by nuclear reactions for the hydrodynamical evolution of the explosion.

A number of improvements have been implemented in the P-HOTB code since the work of Sukhbold et al. (2016).

(1) In solving the neutrino-transport equation by an analytic method (see Scheck et al. 2006), the energy source term associated with neutrino–nucleon scattering (whose average energy transfer per scattering reaction is taken into account) is not described in a closed form as in Scheck et al. (2006) but is split into an absorption part and emission part of energy. This considerably improves the numerical stability of the transport integrator.

(2) For the long-time simulations of supernova explosions over several weeks, the radioactive decays of freshly produced ⁵⁶Ni to ⁵⁶Co (half-life 6.077 days) and ⁵⁶Co to stable ⁵⁶Fe (half-life 77.27 days) are included. This has some influence on the expansion velocities of the innermost ejecta shells (see Jerkstrand et al. 2018) and a minor effect on the late-time fallback.

(3) An adaptive grid-refinement algorithm is added to ensure that the density differences between neighboring grid cells in the PNS surface layers obey the constraint Δρ/ρ ≲ 6%. This improves the possibility of maintaining hydrostatic equilibrium for the increasingly steep density gradient at the PNS surface.

(4) For the long-time simulations during phases of homologous or nearly homologous expansion, a steerage algorithm moves the grid with the expanding matter. This reduces the numerical diffusion of, e.g., composition components, and it also reduces computational costs by allowing for larger time steps.

The first two points have some smaller consequences for the explosion results published by Ertl et al. (2016a) and Sukhbold et al. (2016), and in particular, point (1) leads to slightly more optimistic conditions for explosions. Therefore, with the improved P-HOTB code, we obtain a slightly larger number of exploding cases, in particular at the edges of ZAMS mass intervals where failures or fallback supernovae were obtained before. In Appendix A we provide a detailed comparison of explosion results obtained by Sukhbold et al. (2016) with the previous version of the P-HOTB code and results obtained for all neutrino engines with the present, upgraded version of the code.

The stellar collapse and explosion simulations are computed with a radial grid of up to 2000 zones. After core bounce, the grid is dynamically adjusted, and its number of zones varies between about 1300 and 2000. After usually 10–15 s, the neutrino cooling of the compact remnant is no longer followed, and the inner grid boundary is moved outward. Also, the outer grid boundary is relocated to larger radii as the supernova shock propagates through the progenitor star. We have confirmed by tests that the radial position of the inner grid boundary has no relevant effect on the computation of the fallback, which is determined by mass escaping through the inner grid boundary with negative velocities (see Ertl et al. 2016b).

The supernova explosions are typically followed for about 5 weeks. In order to perform these simulations beyond the breakout of the supernova shock from the stellar surface, the progenitor star is embedded in a dilute circumstellar medium. For this purpose, we assume a constant mass-loss rate ${\dot{M}}_{{\rm{w}}}$ and prescribe a wind velocity of ${v}_{{\rm{w}}}(r)={v}_{\infty }\cdot (1-{R}_{* }/r)$, where R_* is the radius of the progenitor star and ${v}_{\infty }$ is the asymptotic velocity at large distances. This yields a wind-density profile of ρ_w(r) = ${\dot{M}}_{{\rm{w}}}$ · [4πr²v_w(r)]⁻¹. The temperature in the wind is computed by the assumption of hydrostatic equilibrium. For the supernova simulations discussed in this work, we used ${\dot{M}}_{{\rm{w}}}$ = 10⁻⁴ M_⊙ yr⁻¹ and ${v}_{\infty }={10}^{8}$ cm s⁻¹. With these values, the density in the first computational cell exterior to the star is typically 1–2 orders of magnitude lower than the density in the outermost cell of the stellar model. The corresponding total mass of the circumstellar medium within the computational grid with an outer boundary radius of 10¹⁸ cm is lower than 0.04 M_⊙. This mass, in addition to the precollapse mass of the progenitor, is so low that it has no dynamical influence on the expanding supernova ejecta.

KEPLER is used to postprocess the explosions calculated by P-HOTB in order to capture the details of the nucleosynthesis and test for consistent results. This postprocessing is facilitated by the fact that the presupernova models used by P-HOTB were calculated using KEPLER, so the starting points are identical. It is complicated by the fact that P-HOTB uses an Eulerian grid, carries the PNS, and includes neutrino transport, while KEPLER is Lagrangian and does neither. Since regions outside what will eventually be the final mass separation have their dynamics modified considerably by neutrino-energy deposition, postprocessing is not so simple as just taking the history of a single mass shell in P-HOTB and using that as a piston in KEPLER. Particularly problematic is matter that passes through the accretion shock in P-HOTB and dwells for a considerable time between that shock and the neutrinosphere while experiencing an irregular radial history. It may take up to a second for the explosion to develop in P-HOTB. Attempting to model this delay with a piston in KEPLER results in multiple bounces and the artificial accumulation of high-density matter just outside the piston.

In Sukhbold et al. (2016), the approach taken was to define a "special trajectory" in P-HOTB for use as a piston in KEPLER (see their Figures 10 and 11). This trajectory tracked the first mass shell to move outward vigorously once the stalled shock was revived. The radius of that shell contracted from a few thousand km to a minimum of 80–140 km (more typically 120–140 km) before rapidly moving out and being ejected. The mass inside this shell, its minimum radius, and the time at which the matter moved out defined the initial motion of a piston in the KEPLER calculation. The inward motion of that piston, from the time the core first collapsed, i.e., reached a collapse speed of 1000 km s⁻¹ in any zone, to the time it reached its minimum radius was approximated by fitting a parabolic curve to the initial and final radius; see Equation (7) of Sukhbold et al. (2016). Its subsequent outward motion was given by a ballistic trajectory in a modified gravitational field, i.e., a constant times the local Newtonian acceleration due to gravity. Given the assumed final radius of the piston, 10⁹ cm, the strength of this field, and hence the transit time to 10⁹ cm, was adjusted so that the kinetic energy of all ejecta at late time agreed with the P-HOTB simulation. By design then, the explosion energy in KEPLER and the initial mass separation were the same as in P-HOTB. Subsequent fallback was followed independently in both codes and usually agreed well. All remnant masses (e.g., those used in Section 6) were taken from the P-HOTB simulation and not KEPLER.

For the piston trajectory just described, the nucleosynthesis of the iron group and ⁵⁶Ni usually agreed quite well between the two codes. Most of the ⁵⁶Ni made in P-HOTB was produced outside of the special trajectory, so its synthesis was accurately captured by KEPLER. Between the special trajectory and the deeper final mass separation, however, there was matter in P-HOTB with an electron mole number, Y_e, that had been affected by neutrino capture during both the delayed explosion and the neutrino-powered wind that followed. As discussed in Section 3.1, this matter consisted, in P-HOTB, of helium from photodisintegration and iron-group elements, Tr, with an uncertain isotopic composition. For all models in Sukhbold et al. (2016), the mass in the form of Tr was much less than that of ⁵⁶Ni, so fitting the total iron (⁵⁶Ni plus other "neutronized" isotopes) well also resulted in fitting ⁵⁶Ni well. In the end, the total iron-group synthesis in KEPLER always fell between the values of Ni and Ni+Tr in all the matter ejected by P-HOTB. In some cases where this condition was not initially satisfied, the piston location was moved inward in KEPLER until it was. The adjustment was small, of order 0.01 M_⊙. Ultimately, the total ejected iron-group synthesis in KEPLER and P-HOTB agreed to within a few thousandths of a solar mass for models below 12 M_⊙ and about 0.02 M_⊙ for heavier models. The average difference in total iron ejected for the two codes was 24% or 0.014 M_⊙ (see Figure 12 and Tables 4 and 5 of Sukhbold et al. 2016).

While no additional modifications were made to tune the production of ⁵⁶Ni itself, it turned out that the agreement between the two codes was even better for that isotope. The ⁵⁶Ni production in KEPLER, on   average, was within 9.8% (0.005 M_⊙) of the P-HOTB values for just ⁵⁶Ni and 28% (0.016 M_⊙) for ⁵⁶Ni+Tr. In the worst case, the ⁵⁶Ni-to-⁵⁶Ni comparison was off by 0.021 M_⊙.

In the present paper, iron and ⁵⁶Ni synthesis in KEPLER were computed this "old way," but an additional set of models was calculated that used a deeper mass cut for the piston location. Except for a few cases with large amounts of late-time fallback, this location was equal to the final mass separation in the P-HOTB simulation. For those cases with massive fallback, the nucleosynthesis of iron did not matter, and the final remnant mass was taken from P-HOTB, so the special trajectory continued to be used. Moving the mass cut in like this increases the production of iron and ⁵⁶Ni in KEPLER. Because interactions with neutrinos may have reduced Y_e in these deepest layers, the ⁵⁶Ni production is a maximum. The timing of the bounce and the minimum radius in this other set were still that of the "special trajectory," and the same agreement in final kinetic energy was enforced, but the initial location of the piston was deeper and not varied.

Because of the central role of ⁵⁶Ni in producing the light curves of SNe I, we discuss here how our best estimates were calculated and their possible error. Ultimately, six possibilities were considered, each with its own special meaning.

Part of the variation is because KEPLER and P-HOTB both calculate iron-group nucleosynthesis. KEPLER uses a large network that is necessary for determining isotopic and trace element nucleosynthesis, while P-HOTB uses a 13-species network (α-nuclei from helium through ⁵⁶Ni) for temperatures below 7 × 10⁹ K, supplemented by a 15-species NSE solver for higher temperatures (Section 3.1; Kifonidis et al. 2003; Scheck et al. 2006; Ugliano et al. 2012; Ertl et al. 2016a). The dynamical history of zones near the mass cut is tracked more accurately by P-HOTB, which also calculates a 1D approximation to the neutrino-powered wind. To do so, it follows the evolution of ejecta near the mass cut for at least 10 s after core collapse. See Figure 10 of Sukhbold et al. (2016) for an example.

As discussed in Section 3.1, P-HOTB delineates its iron-group synthesis into three components, Ni, Tr, and α, and there is some ambiguity in the interpretation of each. The "Ni" component is the ⁵⁶Ni synthesized in the 13-isotope network and the 15-species NSE solver when Y_e is greater than 0.49. The "Tr" component is ⁵⁶Ni plus other iron-group species produced in both the network and the NSE solver when Y_e < 0.49, especially in the neutrino-powered wind. In the present studies, Y_e was typically ∼0.46–0.49 in the wind region, and thus Tr would contain essentially no ⁵⁶Ni. Other more realistic studies have shown, however, that Y_e in the wind should be close to or even greater than 0.50 (Martínez-Pinedo et al. 2012; Roberts 2012; Mirizzi et al. 2016). Studies by, e.g., Pruet et al. (2006) show the fraction of heavy species that are ⁵⁶Ni is 30%–80% for Y_e near 0.54. Consequently, we take the fraction of Tr that is ⁵⁶Ni here to be some uncertain fraction that might be near 0.5.

The α-particles are produced by photodisintegration in the shock but also in the wind at late times, where their mass fraction can approach 90%. The total mass of material that has achieved NSE in P-HOTB is bounded above by Ni+Tr+α. It is an upper bound, since the sum also includes a small amount of ⁵⁶Ni produced by incomplete explosive silicon burning for shock temperatures between 4 and 5 × 10⁹ K. Since we are interested not just in best estimates but also in upper bounds (Section 8), it is possible that some greater fraction of the α-particles might reassemble in a multidimensional model. Regardless of dimensionality, a comparable amount of matter needs to absorb energy from neutrinos in order to obtain the observed energy of the supernova. In 1D, this mass expands like a wind; in 3D, it will come from overturn, accretion, and re-ejection of matter after heating. An open question is whether the average entropy of the bottommost layers ejected in a 3D explosion is less than in 1D. This would promote the assembly of more α-particles into ⁵⁶Ni. Coupled with the desire to explore the sensitivity of the light curves to the ⁵⁶Ni mass, this motivates treating the α-particles calculated here by P-HOTB as a potential mixture of ⁵⁶Ni and α. Studies with KEPLER with the mass separation at the P-HOTB value (Section 3.2) show that for matter ejected solely by shock heating, the average ⁵⁶Ni mass ejected is 75% of the matter that achieved NSE. This suggests that 0.75 times the sum of Ni+Tr+α from the P-HOTB calculation is an upper bound to the ⁵⁶Ni synthesis. It is an upper bound because KEPLER does not include neutrino capture.

Figure 9, shown later in the paper in Section 4.3, illustrates the possibilities. For the neutrino-powered models calculated here (see Section 8.5 for other possibilities), the least amount of ⁵⁶Ni is mainly the iron-group material that did not experience any appreciable change in Y_e during the explosion as modeled in P-HOTB. This is the thin black line at the bottom of the figure. It also contains the ⁵⁶Ni made by explosive burning in shock-heated ejecta plus ⁵⁶Ni produced by nuclear recombination in neutrino-heated wind matter for Y_e ≥ 0.49. This last contribution is usually very small because the Y_e in neutrino-processed matter that ends up being ejected   is Y_e ∼ 0.46–0.49 in the P-HOTB models. Close to that, the blue line gives the KEPLER results for calculations using the special trajectory as the mass cut, the same approach used by Sukhbold et al. (2016; see Section 3.2). The two dashed lines are for Ni+Tr/2 and Ni+Tr as calculated by P-HOTB. As discussed in Section 3.1 it is probable that not all of the Tr is neutronized iron, but the fraction of ⁵⁶Ni is undetermined when the Y_e in the wind is lower than 0.49. The band between the two is our best estimate for ⁵⁶Ni synthesis in the current model set. The top solid black curve shows the KEPLER result when the mass cut is moved into the deepest possible value for the present models. It is very similar to the red curve, which is all of the ⁵⁶Ni, Tr, and helium from photodisintegration ejected by P-HOTB (i.e., essentially all the matter that has reached NSE) multiplied by 0.75. The factor 0.75 was chosen so that the two curves would match. The best value is thus in the gray band, and the two values explored as maxima and minima for the light curves (Section 8) are the red and blue lines. The results in the figure are for the W18 central engine, but those for S19.8 are similar.

Three different outcomes of stellar core collapse can be distinguished: supernova explosions with neutron star formation, supernova explosions with black hole formation by massive fallback ("fallback supernovae"), and failed explosions with "direct" black hole formation, i.e., continuous accretion of the transiently existing neutron star until collapse to a black hole takes place ("Implosion and BH" in Figures 3 and 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 3 and 4 show our results from P-HOTB simulations for all of the central engines employed for the set of helium stars with a standard mass-loss rate and the set with an enhanced mass-loss rate, respectively. The engines are named by the combination of Crab and SN 1987A progenitors, whose sets of parameter values are interpolated with a dependence on the core compactness (more specifically, M₃₀₀₀; see Section 3.1). The strongest neutrino engine is shown   at the top, the weakest at the bottom. It is evident that the number of cases with supernova explosions and neutron star formation (red) decreases from top to bottom, whereas the number of cases with black hole formation, either by failed explosions (black) or fallback supernovae (blue), increases. While below a final helium star mass of about 6 M_⊙, essentially all progenitors blow up and give birth to neutron stars (with very few exceptions for the weakest engine), there is a mix of outcomes with neutron star or black hole formation for final helium star masses between ∼6 and nearly 12 M_⊙, where the strongest engine produces mostly successful explosions, and the weakest engine produces mostly black holes. Above a final helium star mass of ∼12 M_⊙ for the progenitors with standard mass loss and ∼11.3 M_⊙ for enhanced mass loss, only black hole formation by continuous or fallback accretion takes place.

A subset of the massive helium stars is especially interesting because they produce very massive PNSs before an explosion sets in (the baryonic neutron star masses at this time are between ∼2.3 and ∼2.9 M_⊙). Later, fallback triggers black hole formation (see Figure 2). This is the case for final helium star masses between ∼12 and 14 M_⊙ (ZAMS mass between ∼58 and 68 M_⊙) for standard mass loss, in particular with the three strongest neutrino engines, and for final helium star masses between ∼11.3 and 12.5 M_⊙ (ZAMS mass between ∼77 and 85 M_⊙) for enhanced mass loss and all neutrino engines. Below these mass windows, there is only a single exception of a lower-mass helium star with a similar behavior, namely a progenitor with a final helium star mass of nearly 8 M_⊙ (ZAMS mass ∼39 M_⊙) for the standard mass loss and one with a final helium star mass close to 7.7 M_⊙ (ZAMS mass ∼54 M_⊙) for the higher mass loss. In all of these cases, M₄, i.e., the mass where the dimensionless entropy per baryon reaches a value of s = 4, is large (>1.85 M_⊙). But the explosions set in even later, typically later than 1.5 s after bounce, at which time the neutron star masses are considerably bigger than M₄. Consequently, the final explosion energies are low, less than ∼0.5 × 10⁵¹ erg.

In these cases, M₄ underestimates the neutron star mass at the time the explosion sets in, because the star blows up only when matter whose entropy per baryon equals about 6 in the presupernova star collapses into the stalled shock. This roughly corresponds to the outer boundary of the convective oxygen-burning shell in the presupernova star and is often close to strong carbon- and neon-burning shells in cases where these shells have merged. At this location, the mass accretion rate makes another drop by about a factor of 2. Successful explosions in these cases occur because M₆ (which is the mass enclosed by the radius where s = 6) is still smaller than in other nonexploding progenitors, namely between about 2.3 and 2.9 M_⊙. Using M₄ instead of M₆ as a proxy of the neutron star mass and μ₄ as a parameter scaling with the mass accretion rate places these cases above the two-parameter curve separating neutron star formation from black hole formation in the M₄μ₄–μ₄ plane (see Appendix B and Ertl et al. 2016a for the theoretical background), i.e., on the side of black hole formation. This is appropriate because it reflects the fact that these stars do not blow up when the infalling M₄ shell arrives at the shock, but they explode marginally only later with low energies and several solar masses of stellar material falling back   onto the compact remnant.

In the following discussion of explosion results and their astrophysical consequences, we will focus on the two engines based on the W18 and S19.8 progenitors of SN 1987A. The reason is that the former allows for the most detailed comparison with the single-star results discussed by Sukhbold et al. (2016). The latter is our strongest neutrino engine, marginally beating N20, and thus is suitable to demonstrate the range of possibilities on the optimistic side for explosions. When this aspect is of relevance, we will also refer to our results with the S19.8 engine.

In our P-HOTB simulations, neutrino-energy loss from the newly formed compact remnant is taken into account by the binding energy of the analytic inner-core model plus the neutrino emission from the matter that gravitationally settles in the outer PNS layers followed directly by our hydrodynamics and neutrino-transport simulations. The description is very approximate. Nevertheless, detailed comparisons   (D. Kresse et al. 2019, in preparation) reveal that the total release of neutrino energy, E_ν,tot, overestimates the binding energy computed from Equation (36) of Lattimer & Prakash (2001; see also Equations (9) and (10) of Sukhbold et al. 2016) only by a modest 10%–20% when neutron star radii of 11–12 km are used in these formulae. This difference not only stems from the rough approximations in our neutron star model but also results from the fact that the PNS is transiently more massive than the neutron star after its mass loss by the neutrino-driven wind and before fallback, i.e., M_b,NS(t_exp) > ${M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}$ (see Figure 2). It therefore radiates transiently more neutrinos than expected by the final neutron star mass. Such time-dependent effects of the dynamical supernova evolution in spherical symmetry are not captured by the formulae of Lattimer & Prakash (2001).

The difference implies a small underestimation of the gravitational neutron star mass by typically 0.02–0.04 M_⊙ in the P-HOTB simulations compared to the estimates derived on the grounds of Lattimer & Prakash (2001). For example, for neutron stars with baryonic masses of 1.4, 1.6, 1.8, 2.0, and 2.2 M_⊙, our P-HOTB simulations yield about 1.24, 1.40, 1.56, 1.71, and 1.90 M_⊙ for the gravitational masses, whereas Equation (36) of Lattimer & Prakash (2001) gives 1.27, 1.43, 1.60, 1.75, and 1.90 M_⊙ with R_NS = 12 km and 1.26, 1.42, 1.58, 1.73, and 1.88 M_⊙ with R_NS = 11 km. In the worst case and for a very massive neutron star, the gravitational mass of the cold remnant estimated with our neutrino loss in the P-HOTB simulation is lower by ∼0.1 (0.07) M_⊙ for a neutron star radius of 11 (10) km, corresponding to an underestimation of the gravitational neutron star mass by 4%–5% compared to the values derived from Equation (36) of Lattimer & Prakash (2001).

However, the gravitational wave and kilonova observations of GW170817 have set new constraints on neutron star radii, which reduce the uncertainties of the nuclear EOS in cold neutron stars. The numerical factors in Equations (35) and (36) of Lattimer & Prakash (2001) are averages over a wide range of possibilities, some of which are not compatible with the new radius constraints. Abbott et al. (2018) concluded that the radii of the two merger components (whose masses are most likely between ∼1.1 and ∼1.7 M_⊙) are 11.9_−1.4^+1.4 km for EOSs that permit neutron stars with masses larger than 1.97 M_⊙ as required by observations (Antoniadis et al. 2013). Bauswein et al. (2017) argued that the radius of a nonrotating cold neutron star with a gravitational mass of 1.6 M_⊙ is larger than 10.68_−0.04^+0.15 km, and the radius of the maximum mass configuration must be larger than 9.60_−0.03^+0.14 km. Such small radii tend to favor higher binding energies than obtained with the numerical factors used in Equations (35) and (36) of Lattimer & Prakash (2001). The EOSs that are still compatible with the new constraints—for example, WFF2, AP3, AP4, ENG, and, on the very compact side, WFF1—are above the average line in Figure 8 of Lattimer & Prakash (2001). Therefore, they yield binding energies that are roughly 5%–25% higher than the mean value when neutron stars have gravitational masses exceeding ∼1.4 M_⊙.

For all of these reasons, we will use our P-HOTB results for the neutrino-energy loss of the newborn neutron stars in the whole paper. They are displayed in Figures 5 and 7, and they are employed to convert baryonic masses to gravitational masses of the compact remnants as given in Figures 6 and 8. For reference, we will also provide the numbers derived with Equation (36) of Lattimer & Prakash (2001) when we discuss initial mass function (IMF)–averaged masses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Different from our approach in the previous works on single stars (Ertl et al. 2016a; Sukhbold et al. 2016; see also Appendix A), the greater number of present cases with massive fallback motivates us now to also take into account the release of the gravitational binding energy associated with the accretion of fallback matter. Because the neutrino-energy release calculated with P-HOTB and the values computed from Equation (36) of Lattimer & Prakash (2001) exhibit close agreement, we use the latter formula for estimating the energy loss from accretion. Several cases need to be distinguished.

(1) In the case of direct black hole formation by continuous accretion without supernova explosion, the neutron star is assumed to collapse to a black hole when it reaches a baryonic mass of M_b,NS^max = 2.75 M_⊙.

(1a) If this limiting mass is reached during the P-HOTB simulation, we use the neutrino-energy output E_ν,tot from the simulation to estimate the gravitational mass of the black hole by

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

Here M_b,BH is the total mass of the stellar matter that collapses into the black hole, and we assume that the loss of energy from matter accreted after black hole formation is negligible. This assumption holds well for radial accretion, in which case neutrinos (and, of course, photons as well) do not have enough time to efficiently escape from the fast inward flow. Moreover, we make the approximation that thermal stabilization of the hot neutron star has a minor effect; i.e., we take our assumed mass limit for cold neutron stars to also set the limiting mass of the accreting remnant. This is valid only when the accretion proceeds slowly and the neutron star survives for a long period of time to cool efficiently. However, if the neutron star collapses on a short timescale compared to the cooling timescale, our approximation is acceptable as well, because the neutrino-energy loss in this case is fairly small, and therefore the gravitational-to-baryonic mass difference accounts only for an insignificant correction to the black hole mass.

(1b) If the accreting neutron star does not reach the limiting mass of ${M}_{{\rm{b}},\mathrm{NS}}^{\max }=2.75$ M_⊙ until our P-HOTB simulations are terminated at t_end = 10 s after bounce, we extrapolate its further energy loss by using the Lattimer & Prakash (2001) estimate. The gravitational black hole mass is therefore calculated as

$$\begin{eqnarray}&&{M}_{{\rm{g}},\mathrm{BH}}={M}_{{\rm{b}},\mathrm{BH}}-\displaystyle \frac{1}{{c}^{2}}\,\max \left\{{E}_{\nu ,\mathrm{tot}}({t}_{\mathrm{end}}),{E}_{{\rm{b}},10\,\mathrm{km}}^{\mathrm{LP}01}({M}_{{\rm{b}},\mathrm{NS}}^{\max })\right\},\end{eqnarray} \tag{ 2 }$$

where ${E}_{{\rm{b}},10\,\mathrm{km}}^{\mathrm{LP}01}({M}_{{\rm{b}},\mathrm{NS}}^{\max })$ is the binding energy of a cold, maximum mass neutron star of 10 km radius according to Equation (36) of Lattimer & Prakash (2001). In practice, its value is higher than the neutrino-mass decrement obtained in the P-HOTB simulation at t_end.

(2) In the case of black hole formation by fallback accretion, we compute the gravitational black hole mass according to

$$\begin{eqnarray}&&{M}_{{\rm{g}},\mathrm{BH}}={M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}+{M}_{\mathrm{fb}}-\displaystyle \frac{1}{{c}^{2}}\,{E}_{{\rm{b}},\mathrm{NS}}^{\max },\end{eqnarray} \tag{ 3 }$$

where M_fb is the fallback mass and ${M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}$ the baryonic mass of the transiently stable neutron star before fallback, whose neutrino emission has triggered the supernova explosion. The maximum binding energy of a neutron star is assumed to be

$$\begin{eqnarray}&&{E}_{{\rm{b}},\mathrm{NS}}^{\max }={E}_{\nu ,\mathrm{tot}}({t}_{\mathrm{end}})+{E}_{{\rm{b}},\mathrm{fb}}^{\mathrm{LP}01}.\end{eqnarray} \tag{ 4 }$$

Here we estimate the neutrino-energy release associated with the fallback, ${E}_{{\rm{b}},\mathrm{fb}}^{\mathrm{LP}01}$, by employing Equation (36) from Lattimer & Prakash (2001) for computing the difference between the binding energy of the final neutron star when it collapses to a black hole and the initial neutron star mass before fallback accretion:

$$\begin{eqnarray}&&{E}_{{\rm{b}},\mathrm{fb}}^{\mathrm{LP}01}={E}_{{\rm{b}},10\,\mathrm{km}}^{\mathrm{LP}01}({M}_{{\rm{b}},\mathrm{NS}}^{\max })-{E}_{{\rm{b}},12\,\mathrm{km}}^{\mathrm{LP}01}({M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}).\end{eqnarray} \tag{ 5 }$$

This assumes that fallback accretion proceeds slowly enough, and therefore the accreting neutron star remains stable for a sufficiently long period of time, to allow all of the binding energy to be radiated away.

Neutrino- or radiation-energy loss by further radial (in 1D) accretion after black hole formation is again assumed to be negligible. In the few cases of the standard mass-loss set of progenitors where ${M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}\gt {M}_{{\rm{b}},\mathrm{NS}}^{\max }=2.75$ M_⊙ (see Figure 2, left panel), we expect that the neutron star collapses to a black hole due to neutrino emission after 10 s of simulated post-bounce evolution. Therefore, in these cases, we also consider no additional neutrino-energy loss during the fallback accretion by the black hole.

(3) For consistency, we also include mass decrement corrections due to neutrino and radiation loss from fallback matter accreted by the stable neutron stars in our simulations, although these corrections are minuscule in the vast majority of cases because of the usually small fallback masses. In order to do that, we follow a procedure analog to the one applied in point (2). The corrected gravitational mass of the neutron star is

$$\begin{eqnarray}&&{M}_{{\rm{g}},\mathrm{NS}}={M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}+{M}_{\mathrm{fb}}-\displaystyle \frac{1}{{c}^{2}}\,{E}_{{\rm{b}},\mathrm{NS}}({M}_{{\rm{b}},\mathrm{NS}}),\end{eqnarray} \tag{ 6 }$$

where ${M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}$ is the baryonic mass of the neutron star before fallback, and M_b,NS is the final baryonic neutron star mass including fallback. The binding energy of the final neutron star is taken to be

$$\begin{eqnarray}&&{E}_{{\rm{b}},\mathrm{NS}}({M}_{{\rm{b}},\mathrm{NS}})={E}_{\nu ,\mathrm{tot}}({t}_{\mathrm{end}})+{E}_{{\rm{b}},\mathrm{fb}}^{\mathrm{LP}01},\end{eqnarray} \tag{ 7 }$$

with

$$\begin{eqnarray}&&{E}_{{\rm{b}},\mathrm{fb}}^{\mathrm{LP}01}={E}_{{\rm{b}},12\,\mathrm{km}}^{\mathrm{LP}01}({M}_{{\rm{b}},\mathrm{NS}})-{E}_{{\rm{b}},12\,\mathrm{km}}^{\mathrm{LP}01}({M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}).\end{eqnarray} \tag{ 8 }$$

Again, we estimate the neutrino-energy release associated with the fallback, ${E}_{{\rm{b}},\mathrm{fb}}^{\mathrm{LP}01}$, by employing Equation (36) from Lattimer & Prakash (2001), but now we compute the difference between the binding energy of the final neutron star and the binding energy of the neutron star before fallback by adopting a radius of 12 km both times.

We finally repeat that only in a tiny subset of our present simulations is the assumed baryonic neutron star mass limit of 2.75 M_⊙ exceeded at the time the explosion sets in, i.e., M_b,NS(t_exp) > 2.75 M_⊙ (see Figure 2, left panel). However, the overshoot is so small that we can safely assume that thermal pressure stabilizes the neutron stars long enough for them to power the associated weak fallback supernova explosions by their neutrino emission. We also treat these neutron stars such that they survive for the full simulation time of 10 s to radiate neutrinos. This is a crude assumption. But since these cases have high fallback masses and give birth to black holes of ∼10 M_⊙ or more, this implies only minor errors in our estimates of the black hole masses.

Figures 5–8 summarize the P-HOTB simulations of helium stars with standard and enhanced mass-loss rates that use the Z9.6 and W18 engine. Cases with neutron star formation are shown by red bars, fallback supernovae with black hole formation by blue bars, and black hole formation without explosion by black bars. The direct black hole formation cases are also marked by short vertical black dashes in the upper halves of all panels, and the fallback black hole formation cases are marked by short vertical blue dashes in some panels.

Overall, the explosion energies and timescales, Ni+Tr masses, and remnant masses for the helium stars are quite similar to those for single stars; see Sukhbold et al. (2016) and Appendix A. The minimum values of the explosion energy, E_exp, with the Z9.6 and W18 engine are near 10⁵⁰ erg, and the maximum values are near 1.8 × 10⁵¹ erg. The corresponding Ni+Tr masses range from ∼0.007 M_⊙ for the weakest explosions to ∼0.10 M_⊙ for the most energetic ones. A single outlier in the model set with enhanced mass loss produces slightly more than 0.15 M_⊙ (Figure 7). This star, with an initial helium star mass of 17.5 M_⊙, exploded unusually early with a high energy of 1.75 × 10⁵¹ erg and a correspondingly high yield of Ni+Tr. These special explosion conditions resulted because the base of the oxygen shell was characterized by an exceptionally high entropy jump (roughly from s ≈ 2.5 to a value near 6), leading to a more dramatic drop of the mass accretion rate than for neighboring stars with similar helium star masses and values of M₄.

Zoom In Zoom Out Reset image size

The lowest fallback masses of less than 10⁻⁴ M_⊙ were obtained for the helium stars with the smallest precollapse masses and the lowest explosion energies, while the more typical cases had fallback masses between ∼10⁻³ and roughly 10⁻² M_⊙. Only for progenitors with initial helium star masses above about 15 M_⊙ did fallback significantly greater than 10⁻² M_⊙ become more frequent.

The estimated gravitational masses of the neutron stars ranged from about 1.24 M_⊙ up to ∼2.23 M_⊙ for a progenitor of roughly 57 M_⊙ ZAMS mass (22.50 M_⊙ initial helium star mass) in the standard mass-loss set. For this heavy case, fallback of nearly 1.0 M_⊙ lifted the baryonic mass of the compact remnant to ∼2.71 M_⊙, which was just below our assumed mass limit for cold neutron stars. A similarly high value of M_g,NS ≈ 2 M_⊙ (M_b,NS ≈ 2.4 M_⊙) came from another progenitor of ∼75 M_⊙ ZAMS mass (31.50 M_⊙ initial helium star mass) for the model set with enhanced mass loss.

The black hole masses range from values just above the threshold mass for black hole formation, to ∼19.5 M_⊙ for the most massive helium stars with standard mass loss, to ∼10.7 M_⊙ for the set with a higher mass-loss rate. At the lower end, the black hole mass was set by fallback; at the upper end, the entire presupernova star collapsed.

The most obvious difference of the helium star explosions compared to the single-star models of Sukhbold et al. (2016; see also Appendix A) is the greater number of cases with significant fallback. These fill in a "gap" in the remnant mass distribution that might have existed between gravitational masses of about 2–6 M_⊙. Most cases where black holes are formed by fallback are associated with progenitors that explode very late (t_exp ∼ 2 s) at the time the infalling point with entropy per baryon 6.0 reaches the shock. Such cases form massive PNSs (see Section 4.1) that explode with relatively low energies (about (1–4) × 10⁵⁰ erg) and eject no iron-group material (in the absence of mixing). The fallback includes much of the progenitor mass. Typically, 2–3 M_⊙ is still ejected in models with standard mass loss but only 1–1.5 M_⊙ for the most massive helium stars with enhanced mass loss. We estimate the neutrino-energy release in such cases (including the energy loss from fallback accretion) to be up to ∼10⁵⁴ erg (or about 0.55 M_⊙), very similar to black hole formation by continuous slow accretion.

Also, the black holes with the lowest masses originate from fallback accretion. In the standard mass-loss set, these are just above the gravitational mass threshold for black hole formation, namely 2.31, 2.29, and 2.56 M_⊙ for final helium star masses of 10.91, 11.23, and 11.34 M_⊙, respectively. In the progenitor set with enhanced mass loss, there is a similar case with 2.89 M_⊙ for the fallback black hole in a star with a final helium star mass of 11.38 M_⊙.

In the context of this work, the production of radioactive nickel during the explosion is of particular interest. A detailed discussion of the associated uncertainties, which are connected to our simplifications of the explosion modeling, was provided in Section 3.3. For a subset of stars with initial helium star masses between about 4.6 and 10 M_⊙ and standard mass loss exploded with the Z9.6 and W18 neutrino engine, Figure 9 displays lower and upper estimates of the nickel yields compared to our best estimates. Five different quantities are shown: ⁵⁶Ni calculated using the same approach as Sukhbold et al. (2016), Ni from P-HOTB, Ni+Tr/2 from P-HOTB, Ni+Tr, and 0.75 times Ni+Tr+α from P-HOTB. The thin solid black line for Ni from P-HOTB represents the results shown by the red parts of the histogram bars in the third panel of Figure 5, and the thick dashed black line for Ni+Tr corresponds to the whole histogram bars in the same panel. The analysis of Sukhbold et al. (2016) gives values that we now believe are close to the minimum value, "Ni only" in the figure, but not far below the best estimate, Ni+Tr/2. The upper bound, 0.75 times Ni+Tr+α, is much larger than any of these and is probably a gross upper bound, though more extreme variations, probably not appropriate for a purely neutrino-powered explosion, are discussed in Section 8.5.

In neutrino-driven explosions, independent of the dimensionality of the modeling, the explosion energy, E_exp, is provided by neutrino-heated matter that is expelled from the close vicinity or surface of the neutron star after having absorbed energy from neutrinos. Neutrino heating basically transfers the energy to lift this matter to a gravitationally marginally bound state, i.e., to specific energy _tot = _int + _grv + _kin ≲ 0 (with _int, _grv, and _kin being the specific internal, gravitational, and kinetic energies, respectively). The excess energy that fuels the supernova explosion mainly stems from nuclear recombination in the neutrino-heated matter as it expands and cools down from its initial conditions in NSE (see discussion in Appendix C of Scheck et al. 2006; see also Marek & Janka 2009; Müller 2015; Bruenn et al. 2016; Janka 2017). This leads to the simple relation (Scheck et al. 2006)

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{\exp }}{{10}^{51}\,\mathrm{erg}}\approx \displaystyle \frac{{M}_{\mathrm{rec}}}{0.1\,{M}_{\odot }}\cdot \displaystyle \frac{{\epsilon }_{\mathrm{rec}}}{5\,\mathrm{MeV}},\end{eqnarray} \tag{ 9 }$$

where M_rec is the neutrino-driven ejecta mass releasing the recombination energy, and _rec is the mean recombination energy per nucleon. Typical values of _rec are around 5–6 MeV. These are below the maximum value of 8.8 MeV set free when nucleons recombine to iron-group nuclei because recombination can be incomplete (producing α-particles). Moreover, some fraction of the released energy is usually still needed to overcome the gravitational binding energy of the ejected gas.

In our 1D simulations, neutrino-heated matter has Y_e values between ∼0.45 and 0.49. Therefore, neutrons and protons present at NSE recombine to tracer (Tr) material and α-particles, depending on the entropy and expansion timescale of the ejecta, i.e., M_rec = M(Tr) + M(α) (see Section 3.1). At conditions of Y_e ≳ 0.5, a significant mass fraction of ⁵⁶Ni would be assembled.

Because of Equation (9), the explosion energy correlates with M_rec, and since the mass of explosively produced ⁵⁶Ni in the shock-heated progenitor layers with Y_e = 0.5 should also correlate with the explosion energy, we expect a strong correlation between E_exp and the summed masses of Ni+Tr+α. This correlation is visible in the top and bottom panels of the two plots in Figure 10 and even more tightly in the middle panels. In our simulations, the masses of ⁵⁶Ni and Tr are roughly equal (see Figure 9 and third panels of Figures 5 and 7). Since Ni is formed in nuclear reactions of oxygen and silicon, releasing only ∼1 MeV per nucleon, its fusion contributes only a little to the explosion energy.

Note that in Figure 10, M(Ni+Tr+α) is measured before fallback, which means that it constitutes the mass of these nuclear species that contributes to the energy provided by the explosion mechanism through nuclear recombination and fusion. The total energy that the mechanism would have to make available to unbind the complete star and explode it with energy E_exp is denoted by E_mech in Figure 10. It is defined as

$$\begin{eqnarray}&&{E}_{\mathrm{mech}}={E}_{\exp }-{E}_{\mathrm{bind}}(m\gt {M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}};{t}_{\exp }),\end{eqnarray} \tag{ 10 }$$

where ${E}_{\mathrm{bind}}(m\gt {M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}};{t}_{\exp })$ is the (negative) total energy (internal plus gravitational plus kinetic), measured at the time the explosion sets in, of all matter exterior to the baryonic mass of the neutron star before fallback. Consistent with our hydrodynamic modeling, we take general relativistic corrections into account in calculating the gravitational energy of all matter, $m\gt {M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}$. The bounding mass shell, $m={M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}$, is located below the surface of the PNS at time t_exp, and the overlying shells are blown out by neutrino heating in the neutrino-driven wind to deliver the energy to the explosion as described above. Therefore, E_mech is much larger than the final explosion energy E_exp. On the one hand, this is so because E_mech includes the energy that the explosion mechanism has to spend on lifting the matter that contains M(Ni+Tr+α) out of the gravitational potential trough of the neutron star and unbinding all of the overlying stellar material ahead of the supernova shock. On the other hand, fallback removes some of the energy initially transferred to the transiently expanding stellar matter, which can lead to a reduction of E_exp. Measuring ${E}_{\mathrm{bind}}(m\gt {M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}};{t}_{\exp })$ at t_exp is an approximation because the neutron star contracts with time and its surface layers become gravitationally more strongly bound as time goes on. However, our approximation yields a useful proxy because the neutrino-driven outflow is strongest in the phase when the explosion sets in and shortly afterward.

The correlation between E_mech and M(Ni+Tr+α) in the middle panels of Figure 10 is considerably tighter than between E_exp and M(Ni+Tr+α) in the top two panels. In the latter panels, the uppermost data points describe the relation of Equation (9) very closely. These are the normal explosions with very little fallback (typically much less than ∼10⁻² M_⊙; see Figures 6 and 8). The cases that lie well below this imagined correlation line in the top panels, while their correspondents in the middle panels are well integrated in the main data band, are supernova explosions with more significant fallback. In these cases, the binding energy of the collapsed stellar core is considerably larger, and a bigger fraction of E_mech is needed to lift the ejecta out of the gravitational potential trough. Therefore, the resulting explosion energies are noticeably lower than the typical values.

Zoom In Zoom Out Reset image size

The cluster of data points below the main data band at M(Ni+Tr+α) between ∼0.05 and ∼0.10 M_⊙ in the top panels of Figure 10 is even more extreme. These are our cases of fallback supernovae with black hole formation caused by very large fallback masses as discussed above. They are characterized by collapsed stellar cores that possess particularly high binding energies because the explosion sets in very late and the neutron star is extremely massive and compact. These conditions permit only marginal explosions because (a) the binding energy in E_mech is much bigger than usual and (b) M(Ni+Tr+α) is rather low when the explosion starts late. Since the mechanism cannot provide the energy to unbind the whole star, massive fallback is the consequence, and E_exp becomes small. Therefore, the data points are located close to the bottom of the top panels in Figure 10.

Figure 10 demonstrates that the mass of Ni+Tr+α produced by neutrino-driven explosions is limited due to its correlation with the explosion energy. The latter is not found to exceed ∼2 × 10⁵¹ erg even for our strongest neutrino engine, Z9.6 and S19.8 (see right plot in Figure 10). For this reason, we do not obtain an M(Ni+Tr+α) higher than ∼0.2 M_⊙. For the Z9.6 and W18 engine, the limit is even only around ∼0.15 M_⊙ (corresponding to an extreme upper bound on ⁵⁶Ni of ∼0.12 M_⊙, indicated by the red line in Figure 9). In most cases, fallback reduces the ejected mass of Ni+Tr+α, M_ej(Ni+Tr+α), only slightly compared to the nucleosynthesized mass represented by M(Ni+Tr+α). However, the data points for the cases with more massive fallback are shifted considerably in the two bottom panels of Figure 10, where M_ej(Ni+Tr+α) is plotted on the abscissa instead of M(Ni+Tr+α). This concerns most of the points with small explosion energies in the lower parts of the two top panels. They are moved closer to the correlation line describing Equation (9). This leads to a narrower band of data points in the bottom panels and a correlation similarly as tight as the one seen in the two middle panels. Outliers are only the (few) cases in which all or most of the nucleosynthesized Ni+Tr+α falls back. The corresponding data points cluster in the left part of the bottom panels, where the ejected mass M_ej(Ni+Tr+α) is very small or zero. Since mixing is not included in our 1D simulations, it is unclear whether multidimensional effects during the long-time evolution of the supernova can change this result.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the new helium star models, a considerably larger number of cases show appreciable fallback than in our previous studies of single stars (Ertl et al. 2016a; Sukhbold et al. 2016). This fallback results in neutron stars with high gravitational masses up to more than 2 M_⊙ and black holes with smaller gravitational masses between 2.3 and 6 M_⊙. Such cases were absent in our previous works on single-star explosions. Why?

Figure 11 compares the fallback conditions of the helium stars with the explosions of the single-star progenitors considered by Sukhbold et al. (2016). All simulations, including those of single stars, were calculated using the upgraded version of P-HOTB. As discussed in Section 4.3, large fallback often happens when a high energy, E_mech (see Equation (10)), is required in order to gravitationally unbind the whole star and provide an energy of E_exp. This quantity can be used to measure the probability of massive fallback in the explosion. In the left panels of Figure 11, the fallback mass is displayed as a function of the ratio of E_mech to the absolute value of the total energy E_bind as used and defined in Equation (10). In the right panels of Figure 11, ${E}_{\mathrm{mech}}/| {E}_{\mathrm{bind}}|$ is replaced by the corresponding ratio of energies when the binding energy of the star is not computed from the structure at t_exp but from the precollapse profile of the star. This means that E_bind^prog is the (negative) total energy of all progenitor shells above the initial mass cut, i.e., of all mass shells $m\gt {M}_{{\rm{b}},\mathrm{NS}}^{\mathrm{bfb}}$. With that, ${E}_{\mathrm{mech}}^{\mathrm{prog}}={E}_{\exp }-{E}_{\mathrm{bind}}^{\mathrm{prog}}$. The upper panel of each plot shows the results for helium stars, and the lower panel shows single stars for comparison.

Zoom In Zoom Out Reset image size

The figure shows, as expected, that massive fallback anticorrelates with these energy ratios. In all panels of Figure 11, there is a narrow main band that contains the majority of all data points. Outliers from this main band, where most of the models cluster, are only low-mass helium star and single-star progenitors that explode with very little fallback and are therefore located below the main band. For helium stars, fallback masses of more than about 0.1 M_⊙ become frequent when ${E}_{\mathrm{mech}}/| {E}_{\mathrm{bind}}| \lesssim 1.1$ or ${E}_{\mathrm{mech}}^{\mathrm{prog}}/| {E}_{\mathrm{bind}}^{\mathrm{prog}}| \lesssim 2.4$. This implies that in explosions with high fallback masses, the final explosion energy is less than ∼10% of the energy needed to gravitationally unbind the initial ejecta, or E_exp is less than 1.4 times the gravitational binding energy of the progenitor shells exterior to the initial mass cut. These values are lower for single-star progenitors, where M_fallback ≳ 0.1 M_⊙ for ${E}_{\mathrm{mech}}/| {E}_{\mathrm{bind}}| \lesssim 1.07$ or ${E}_{\mathrm{mech}}^{\mathrm{prog}}/| {E}_{\mathrm{bind}}^{\mathrm{prog}}| \lesssim 1.7$, which implies ${E}_{\exp }\lesssim 0.07\,| {E}_{\mathrm{bind}}|$ and ${E}_{\exp }\lesssim 0.7\,| {E}_{\mathrm{bind}}^{\mathrm{prog}}|$. Since E_bind is similar for single and helium stars, this means that the former explode less energetically when the fallback masses are ≳0.1 M_⊙. This seems to be a consequence of the fact that there is not only a reverse shock from the C+O/He interface but also a second reverse shock from the interface between the helium core and the hydrogen envelope. Both together damp the supernova blast more strongly than just the reverse shock from the C+O/He interface.

Massive fallback occurs for helium stars with final helium star masses between ∼10 and ∼14 M_⊙, whereas in the set of single stars, only a few cases with final helium core masses around 10 M_⊙ show this behavior. There are two reasons for more fallback cases in our current P-HOTB simulations of helium stars. First, in the set of single-star progenitors used by Sukhbold et al. (2016), there were many fewer models with final helium core masses from 10.5 to 13 M_⊙, where both the helium stars and the single stars make fallback supernovae. Unfortunately, our model set of 2016 was sparse in exactly this region (see Table 2 in Sukhbold et al. 2016). Second, the helium stars have lower compactness values for precollapse helium star masses between ∼10 and ∼12 M_⊙, namely ξ_2.5 ≈ 0.15–0.2 (Figure 12, upper two panels), whereas the single stars possess ξ_2.5 ≈ 0.2–0.3. This favors direct black hole formation for single stars and more fallback supernovae for helium stars.

Zoom In Zoom Out Reset image size

In the panels that show single-star results in Figure 11 for both neutrino engines (Z9.6 and W18 and Z9.6 and S19.8), one can see a corresponding gap of explosion cases around ${E}_{\mathrm{mech}}/| {E}_{\mathrm{bind}}| \sim 1.05$ and ${E}_{\mathrm{mech}}^{\mathrm{prog}}/| {E}_{\mathrm{bind}}^{\mathrm{prog}}| \sim 1.5$, which separates the few cases with the largest fallback from the continuum of data with fallback masses M_fallback ≲ 0.1 M_⊙. For helium stars, such a gap is absent (for the Z9.6 and W18 engine) or much less pronounced (for the Z9.6 and S19.8 engine).

Figure 12 displays the core-collapse outcomes, i.e., supernova explosion and neutron star formation (red), explosion and black hole formation (blue), or collapse with black hole formation (black), for the Z9.6 and W18 neutrino engine, with a dependence on the compactness parameter ξ_2.5 as a function of the final helium star mass for both model sets with standard and enhanced mass loss (top panels; a similar pattern for single stars was found by Pejcha & Thompson 2015; Müller et al. 2016). In comparison, the bottom panel shows the separation of neutron star formation cases (red) from black hole formation cases (black and blue) by the Ertl et al. (2016a) two-parameter criterion for the same neutrino engine. Compactness values ξ_2.5 ≲ 0.2 favor explosions, compactness values above this threshold characterize mostly core collapse with black hole formation, but for values between 0.1 and 0.2, both outcomes are possible.

In contrast, the Ertl et al. (2016a) criterion separates neutron star formation (red) from black hole formation (blue and black) more reliably. This can be seen in Figure 13, where the patterns for direct simulation results with P-HOTB and core-collapse outcomes predicted by the Ertl et al. (2016a) criterion agree very well for helium stars with both standard and enhanced mass loss. Essentially, only fallback supernovae, which are weak explosions with high fallback masses and black hole formation, are hard to predict with the Ertl et al. (2016a) criterion.

Zoom In Zoom Out Reset image size

Figure 14 compares the results for our current survey of mass-losing helium stars to our previous results (Sukhbold et al. 2016). The x-axis in each case is the presupernova mass for the helium stars or the mass of the helium and heavy-element core of the presupernova models for single stars. The mass ranges for the two are slightly different. Here it is assumed that helium stars below 2.50 M_⊙ (presupernova masses below 2.1 M_⊙) result in neon–oxygen white dwarfs. There could be a narrow range of electron-capture supernovae below this mass, but their existence and mass range are uncertain, and they are not included here. The single stars are bounded on the lower end by a ZAMS mass of 9.0 M_⊙ that has a helium core of 1.57 M_⊙ when the iron core collapses. This limit is also set by the neglect of electron-capture supernovae. The limits are the same in the sense that they are both the lightest models to experience iron core collapse. They are different because radius expansion in the helium star models is thought to lead to a second phase of binary mass transfer.

Zoom In Zoom Out Reset image size

The models also differ at higher masses. Above about 35 M_⊙, the single stars lose their hydrogen envelope and experience WR mass loss, similar to the present models. But a different, larger mass-loss rate was used in 2016, so the final masses were smaller. For example, the presupernova mass for a 100 M_⊙ star in Sukhbold et al. (2016) was 6.04 M_⊙. Here it would be 22 M_⊙. Both stars collapsed to black holes with gravitational masses approximately equal to their presupernova masses.

Figure 14 shows a remarkable overall similarity between the earlier single-star results and the present study. To first order, it is the presupernova mass of the helium and heavy-element core (aka "the helium core") that determines its explosion energy, remnant mass, and ⁵⁶Ni nucleosynthesis. This has important implications. The minimum and maximum iron yields derived here probably characterize all massive stars, whether single or in binaries, even though the core masses at the time of collapse may have resulted from very different evolutionary paths and main-sequence masses.

In detail, however, they are different. The compactness of the core, however measured (Figure 12), is different in the two sets (Figure 14). This reflects the different entropy and composition of the two stars, one being derived from a helium core that lost mass—including, at the high end, most of its helium—and the other from a helium core that grew until the end. As noted by Woosley (2019), there is an offset in the mass giving rise to the first peak in compactness. For single stars, the peak is at about 8 M_⊙; for the present models, it is more like 9 M_⊙. This reflects in part the larger carbon abundance following helium depletion in the helium stars (Woosley 2019) and ultimately is a consequence of central carbon burning transitioning to the radiative regime at higher progenitor mass in the binary models (see Sukhbold & Adams 2020).

As a result, there is a gap in neutron star production around 8 M_⊙ in the old models that is shifted to 9 M_⊙ in the new ones.

Of the many values in Table 3, we suggest the range 1.35–1.38 M_⊙ as the most appropriate for the median gravitational mass of neutron stars in close binary systems. This is the range of results for the median value using the W18 central engine and either the P-HOTB or Lattimer and Prakash corrections to the baryonic mass.   It is the average for systems of all initial masses. The table also gives the partition of this overall average into intervals of initial helium star mass: 3–5, 5–8, and >8 M_⊙. As expected, the average increases with progenitor mass as the core structure becomes less centrally condensed (Figure 15). The weighted birth function for neutron stars is given as a function of their gravitational mass in Figure 16. The mean mass (as opposed to median mass) overall is 1.37–1.40 M_⊙ with a standard deviation of 0.11 M_⊙.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The lightest neutron star in the standard (W18) set has a baryonic mass of 1.332 M_⊙, corresponding to a gravitational mass of 1.239 M_⊙ using the P-HOTB correction for neutrino losses (Section 4.2). If one uses the Lattimer & Prakash (2001) approximation instead, the gravitational mass is 1.214 M_⊙ for an assumed radius of 12 km and 1.205 M_⊙ for a radius of 11 km. Looking at all possible central engines, the overall minimum is for the 2.50 M_⊙ helium star using the W20 central engine. This combination gives a gravitational mass of 1.238 M_⊙ using the P-HOTB correction and 1.214 M_⊙ (1.204 M_⊙) using the Lattimer & Prakash (2001) approximation and an assumed radius of 12 km (11 km). Slightly smaller values are, in principle, possible using an artificial mass separation. For the 2.70 M_⊙ helium star and W20 central engine, the iron core is 1.288 M_⊙, and the equivalent gravitational mass is 1.181 M_⊙ (1.167–1.178 M_⊙ for the Lattimer–Prakash mass with a radius of 12 or 11 km). This is unrealistically small, since mass separations this deep are not found in modern models. We conclude that a gravitational mass of 1.24 M_⊙ is a reasonable lower bound for our models, and 1.20 M_⊙ is the boundary of what might be reasonable with the physics employed in this study. Less massive neutron stars would require special conditions for their creation.

The most massive neutron star comes from the theoretical upper limit, taken here to be about 2.25 M_⊙ (Section 4). Thus, model He22.50 produces a neutron star with gravitational mass 2.23 M_⊙, while model He22.00 makes a black hole with mass 2.29 M_⊙. This continuum of masses is due to the action of fallback after the explosion, and low-mass black holes, while rare, should exist, as should a smattering of neutron stars up to the maximum mass. All neutron stars with gravitational mass over 1.65 M_⊙ come from helium stars over 20.0 M_⊙. All neutron stars with masses below 1.40 M_⊙ come from helium stars lighter than 10 M_⊙, as do a lot of more massive neutron stars. All neutron stars less massive than 1.30 M_⊙ come from helium stars less than 3.5 M_⊙ (Figure 15). Neutron stars between 1.24 and 1.7 M_⊙ are formed promptly in the explosion, but heavier neutron stars are produced by fallback (e.g., the six most massive points in Figure 15).

Figure 16 also suggests that the neutron star birth function has structure with a hint of peaks around 1.35 and 1.5 M_⊙. A larger set of models that explores variable metallicity, mass loss, explosion characterization, and binary parameters will be needed before this structure is statistically significant, but an exciting prospect is that the neutron star mass function reflects the nonmonotonic nature of the compactness parameter (Figure 12), which, in turn, is sensitive to the location of critical shell-burning episodes in the presupernova star.

The values derived here for median masses (1.35–1.38 M_⊙ favored) are in good agreement with those determined from observations. Lattimer (2012) gave an (error-weighted) observational mean of 1.368 M_⊙ for X-ray and optical binaries, 1.402 M_⊙ for neutron star binaries, and 1.369 M_⊙ for neutron star–white dwarf binaries. These are not far from the average neutron star mass found by Schwab et al. (2010), 1.325 ± 0.056 M_⊙.

Özel & Freire (2016) gave the range of accurately determined neutron star masses in 2016 as 1.17–2.0 M_⊙. Our calculated range, 1.24–2.25 M_⊙, is in reasonable agreement. As already remarked, a 1.17 M_⊙ neutron star may be difficult to produce. They also gave 1.33 M_⊙ with a dispersion of 0.09 M_⊙ as the average for double neutron stars; 1.54 M_⊙ with a dispersion of 0.23 M_⊙ for the recycled neutron stars; and 1.49 M_⊙ with a dispersion of 0.19 M_⊙ for the slow pulsars. A recent study by Antoniadis et al. (2016) raised the possibility of two peaks (see also Valentim et al. 2011) within the recycled millisecond pulsar population, with the first peak at 1.388 M_⊙ with a dispersion of 0.058 M_⊙ and the second peak at 1.814 with a dispersion of 0.152 M_⊙ (see single-star comparison in Raithel et al. 2018).

Also given in Table 3 are the neutron star masses resulting from taking the mass separation to arbitrarily occur at the edge of the iron core. This gives a lower limit neutron star mass and a maximum production of ⁵⁶Ni that is of interest for the light curve. These masses are consistently small compared with the observations, suggesting that such a deep mass cut is, on average, probably not realistic.

The neutron star masses based on the same approach as in Sukhbold et al. (2016) are also given for comparison. Since the "special trajectory" (Section 3.2) is more shallow, the neutron star masses are slightly higher.

Stellar collapses that fail to create a strong outward moving shock before 10 s in P-HOTB are assumed to form black holes. While Woosley (2019) gave results for helium cores up to 120 M_⊙ (presupernova masses up to 60 M_⊙), only those that had a reasonable likelihood of leaving a neutron star were followed with P-HOTB. It is expected that for helium stars of up to 60 M_⊙ (presupernova up to 30 M_⊙), the presupernova mass will collapse to a black hole with little mass ejection. For helium cores of 60–70 M_⊙ (presupernova 30–35 M_⊙), a mild pulsational instability is encountered that does not greatly affect the mass of the star when its iron core collapses. To good approximation, these stars also make black holes with masses equal to their presupernova mass. For helium stars initially above 70 M_⊙, the pair instability becomes an important consideration. As far as remnant masses go, the chief effect of the pair instability is to reduce the mass of the remnant. Woosley (2019) derived an upper limit to black hole masses coming from pulsational pair instability supernovae of 46 M_⊙. To good approximation, the black hole mass for masses bigger than 20 M_⊙ for the assumed mass-loss rate is given by

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}\ =\ 0.463{M}_{\mathrm{He},{\rm{i}}}+1.49\,{{\rm{M}}}_{\odot },\end{eqnarray} \tag{ 11 }$$

where M_He,i is the initial mass of the helium star, and

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}\ =\ 0.232{M}_{\mathrm{ZAMS}}-1.23\,{{\rm{M}}}_{\odot },\end{eqnarray} \tag{ 12 }$$

where M_ZAMS is the main-sequence mass of the star.

Figure 17 shows the resulting IMF-weighted birth function for black holes using the P-HOTB results and an extrapolation from 19 to 46 M_⊙ using the above equations (Woosley 2019). A more thorough study that explicitly includes the pulsational pair-instability supernovae results is planned.

Zoom In Zoom Out Reset image size

The smaller-mass black holes are made by fallback after the initial launch of a successful shock. The lowest-mass black holes made in any models for the W18 central engine were 2.29, 2.31, and 2.56 M_⊙, which came from stars with initial helium star masses of 22.0, 21.25, and 22.25 M_⊙, respectively. The presupernova masses of these systems were 11.23, 10.91, and 11.34 M_⊙, respectively, so most of the mass was ejected in the explosion. There is no low-mass "gap" in the black hole mass distribution in the present study. Fallback produces a continuous spectrum of masses from the heaviest neutron stars to the lightest black holes. These low-mass black holes are rare, however, compared with their higher-mass counterparts.

The smallest black hole to be made directly in a failed explosion with no outgoing shock has a baryonic mass of 6.95 M_⊙ and a gravitational mass of 6.42 M_⊙. The low-mass peak in the birth function (Figure 17) below ∼6 M_⊙ is due to fallback; the remainder is mostly black holes produced by the direct implosion of the progenitor star. From 5 to 15 M_⊙, the distribution shows structure correlated with the compactness of the presupernova stellar core (Figure 12). A pronounced minimum at ∼10 M_⊙ results from the island of explodability near that mass. This suppression of black hole formation is robust for different choices of central engine, though its location may vary depending on the details of the presupernova model, especially the rate for ¹²C(α, γ)¹⁶O and convection physics. It should be a target of future gravitational wave surveys.

Table 4 gives the average properties of the distribution. Black holes are produced by 21% of all models. The median gravitational mass (after all neutrino losses) for the W18 central engine when all black holes up to 46 M_⊙ are included is 10.88 M_⊙. This is influenced by the extrapolation of the birth function from 20 to 46 M_⊙ in Figure 17 and is probably an overestimate. If only helium stars up to 40 M_⊙ (main-sequence masses up to 80 M_⊙) are considered, the median is 8.96 M_⊙. Both the IMF for stars above 80 M_⊙ on the main sequence and their mass-loss history, especially as luminous blue variables, are poorly known. A higher mass-loss rate would act to suppress the high-mass tail in Figure 17. Neither the median nor mean is a well-defined characteristic of the structure seen in Figure 17, especially if it is bimodal. There may not even be any black holes with masses near 10 M_⊙. Using a more energetic central engine, S19.8, actually increases the median value because more of the lower-mass models explode while the robustly collapsing massive models remain unchanged.

Given these uncertainties, our results compare favorably with black hole masses observed in X-ray binaries. Özel et al. (2010) gave a mean mass of 7.8 ± 1.2 M_⊙. The observed mean is consistent, but the narrow distribution is not. Farr et al. (2011) gave a lower limit at the 99% confidence level of about 4.5 M_⊙. This is inconsistent with the presence of a substantial number of less massive black holes in our results (Figure 17). These black holes are all made by fallback, though, and their masses are uncertain. On the other hand, Kreidberg et al. (2012) pointed out that errors in the inclination angle may have led to the overestimate of black hole masses in the gap and that black holes may exist all the way down to the maximum neutron star mass, as our analysis suggests. Their analysis also suggests a peak around 7–8 M_⊙ and a sharp drop-off above 10 M_⊙. Figure 17 shows such a peak and a falloff, but the production remains substantial at higher masses. As remarked earlier, the abundance of such very massive black holes is sensitive to the treatment of mass loss and the IMF. All models with presupernova masses over 12 M_⊙ collapse directly to black holes. The issue is only how many such presupernova stars there are.

Fryer et al. (2012) also explored the black hole birth function using a different approach to modeling the explosion and found results that were sensitive to the mass loss. For solar metallicity, they found a maximum black hole mass around 10 M_⊙ when using presupernova models from Woosley et al. (2002), which had a larger mass-loss rate than the present study (see their Figure 9). On the other hand, using presupernova models from Limongi & Chieffi (2006) that employed a smaller mass-loss rate gave larger black hole masses. Clearly, an accurate birth function will require careful consideration of the effects of metallicity and mass loss. For now, we note that the outcome—explosion or collapse—and the mass of the collapsing object are mainly determined by the presupernova mass. It is thus feasible to estimate the black hole birth function to good accuracy from any choice of IMF, metallicity, mass-loss rate, binary parameters, etc. simply from an estimate of the masses of the stars at carbon ignition and the results shown in Figures 5–8.

For the models presented here (solar metallicity, W18 central engine, standard mass loss), 21% made black holes. This includes all stars with final masses above about 12 M_⊙ (initial helium star mass 25 M_⊙). This is substantially less than the 33% that make black holes in single stars for the same central engine (Table 6 of Sukhbold et al. 2016). This reflects both the smaller presupernova masses for stars of a given main-sequence mass when evolved in a mass-exchanging binary and their smaller compactness parameters (Figure 10 of Woosley 2019).

Since all presupernova stars with masses above 12 M_⊙ collapse, there is a maximum luminosity of any SN Ib or Ic progenitor, 10^5.6 L_⊙. The equivalent value for single stars is 10^5.3 L_⊙ or a birth mass of about 23 M_⊙ (Sukhbold et al. 2018), though the current observational limit is 10^5.1 L_⊙ or a birth mass of about 18 M_⊙ (Smartt 2009, 2015).

The opacity consists of two parts, electron scattering and other processes that depend on details of the composition, ionization, velocity shear, and level populations. This second part will be referred to loosely as the "line opacity" or the "additive opacity." Electron scattering is calculated very well in the KEPLER code. A Saha equation is solved in every zone at every time step using the current temperature, density, and composition (Ensman & Woosley 1988). All ionization stages of all of the even-Z elements up to iron are calculated assuming thermal equilibrium. The resulting electron density is used to determine the electron scattering opacity, which is corrected for Compton scattering (including the change of a photon's energy in a scattering event), degeneracy, and the presence of pairs.

The line opacity, on the other hand, is treated poorly. It is represented by a single constant added to the electron scattering opacity at all times and places. In calculations of SNe IIp, this number is given a small value, 10⁻⁵ cm² g⁻¹, that does not significantly affect the light curve. For SNe Ia, where it is generally recognized that, owing to the large metal content, lines play a dominant role, the constant frequently used is κ_a ∼ 0.1 cm² g⁻¹, where "a" stands for "additive." Figure 24 of Dessart et al. (2015) shows the actual opacity for a typical SN Ib model calculated using the CMFGEN code (Hillier & Dessart 2012), a state-of-the-art supernova radiation transport code. The non–electron scattering part is by no means constant but is larger near the center, where the nickel is concentrated, and at earlier times, when the density is higher. If one were to attempt to assign a single value to the difference in the dashed and solid lines in their figure (electron scattering and total opacity, respectively), it would be close to κ_a ∼ 0.2 at early times for velocities less than 6000 km s⁻¹ but much smaller near the peak of the light curve. There it ranges from about 0.07 near the center (v less than 3000 km s⁻¹) to at most a few hundredths cm² g⁻¹ farther out. Since we are more interested in the behavior at peak than on the rise, κ_a was chosen to be 0.03 cm² g⁻¹. This low value may slightly underestimate the rise time compared with that of Dessart et al. (2012) but approximately captures the behavior near peak.

To test this calibration, the light curve of model 4.41 in Dessart et al. (2012) was recalculated using the KEPLER code and their assumption for mixing, which, in their most extreme case, is very similar to ours (Section 7.2). This entailed a running boxcar average with a box size of 0.2 M_⊙ passed through the model of Yoon et al. (2010) four times, a procedure identical to that used by Dessart et al. (2012). Several choices for κ_a were explored (Figure 18). In their Figure 7, the peak for their model 4.41 × 4 is 10^42.47 erg s⁻¹, while ours is 10^42.48 erg s⁻¹ for κ_a = 0.03 cm² g⁻¹ and 10^42.42 erg s⁻¹ for κ_a = 0.05 cm² g⁻¹. Our two light curves peak at 27.5 and 30.6 days after explosion, while theirs peaks at 34.2 days. The discrepancy probably reflects the larger opacity early on in their models. The decline rate at late times is sensitive to the treatment of gamma-ray trapping, which is captured better in their full Monte Carlo simulation. All in all, though, if a single constant must be used, κ_a = 0.03 cm² g⁻¹ seems to be a good value for calculating L_peak.

Zoom In Zoom Out Reset image size

The principal effect of mixing is to disperse ⁵⁶Ni from the central tenth or so solar masses where it is made to masses and speeds farther out in the star. The mixing in SNe I has its origin in the nonspherical nature of the central explosion, not so much the reverse shock as in SNe II. Mixing leads to a light curve that rises earlier, is a little broader, and declines earlier. Except for filling in a gap in emission at earlier times, the modification of the light curve near peak is not great. The effect on the spectrum (not calculated here), though, is very important (Dessart et al. 2012, 2015).

Mixing is inherently a 3D phenomenon and does not imply the homogenization at the atomic scale of any ejecta. Rather, clumps of one composition, which may carry radioactivity, are mixed out in velocity. As noted above, the average effects of mixing are incorporated in the KEPLER code using a "running boxcar average." A specified interval of mass (the "boxcar") is thoroughly mixed, and the process is repeated as this boxcar is moved out, zone by zone, through the star. The degree of mixing depends on the size of the boxcar and the number of times it passes through the star. The resulting distribution is quasi-exponential in mass.

Dessart et al. (2012) found that appreciable mixing was essential to obtain the correct colors and spectra for SNe Ib and Ic. Indeed, the degree to which ⁵⁶Ni was mixed into helium was the determining factor in whether the supernova was of Type Ib or Ic. Lacking the tamping influence of a massive hydrogen envelope, mixing in an almost stripped core extends farther out in velocity space than the 3500 km s⁻¹ typical of models for SN 1987A (Wongwathanarat et al. 2015; Utrobin et al. 2019). There is a lack of 3D calculations of mixing in supernovae without hydrogen. The only exceptions are simulations of ultrastripped supernovae by Müller et al. (2018) and an SN IIb model of Wongwathanarat et al. (2017). The latter model, which had a hydrogen envelope of only 0.3 M_⊙ and therefore did not experience much deceleration by a reverse shock from the hydrogen-rich layer, was used to calibrate mixing in the present study.

The model they considered was the stripped-down core of a 15 M_⊙ red supergiant with a helium core of 4.5 M_⊙. For comparison, we used our 6.0 M_⊙ model, which had a final helium–plus–heavy-element core of 4.45 M_⊙. Figure 7 of Wongwathanarat et al. (2017) shows the angle-averaged ⁵⁶Ni mixing from their 3D run. The Wongwathanarat et al. (2017) model may have experienced more mixing than most others (A. Wongwathanarat 2019, private communication), so it is treated as an upper bound. The prescription finally employed in KEPLER used a boxcar width of 0.15 M_ej, where M_ej is the mass ejected in the supernova (presupernova mass minus remnant mass). This mixing region was moved through all of the ejecta three times, and then a final fourth mixing was applied using an interval half as great, i.e., 0.075 M_ej. The comparison with the model from Wongwathanarat et al. (2017) is shown in Figure 19. The result of the same mixing applied to the 2.5, 6.0, and 19.75 M_⊙ models is shown in Figure 20.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Models with initial helium core masses less than about 3.0 M_⊙ (3.3 M_⊙ if silicon flashes are strong) and standard mass loss experience significant radius expansion during their final stages of evolution (see Table 3 and Figure 5 of Woosley 2019). These stars have presupernova masses less than 2.45 M_⊙ (Table 1) and eject less than 1 M_⊙. Their explosion energies and ⁵⁶Ni production are smaller than the more massive models (Figure 5), so the radioactive portion of their light curves is fainter and faster. Some representative light curves are shown in Figure 21.

Zoom In Zoom Out Reset image size

These are not what we usually think of as common SNe Ib or Ic. Their early display is dominated by envelope recombination. There are often two peaks, though the first is very blue. Some of these might be associated with "fast blue optical transients" (Kleiser et al. 2018; Woosley 2019), for example, models He2.60 and He2.70 with presupernova masses of 2.15 and 2.22 M_⊙ and radii of 7.8 × 10¹² and 6.5 × 10¹² cm, respectively. The properties of these two supernovae are summarized in   Figure 21 and Table 6. These large presupernova radii, small ejected masses, and extended epochs of high effective temperature (Figure 21, lower panel) associated with shock breakout are also similar to what was inferred for SN 2013ge by Drout et al. (2016).

Particularly intriguing are the light curves of the low-mass models that experience a silicon deflagration and eject substantial matter well in advance of core collapse. As discussed in Woosley (2019), the amount of matter ejected and its timing depend on the poorly understood 3D aspects of the silicon flash and are uncertain. In cases where a lot of mass is lost months to years before collapse, the resulting supernova can be very bright.

Figure 22 and Table 7 show a representative range of possibilities with the brightness of the light curve being mostly determined by the timing and energy of the silicon flash. Typically, the flash occurs several weeks prior to iron core collapse. The mean velocity of the ejecta is v_ej,Si ∼ (2 KE_ej,Si/M_ej,Si)^1/2, and models where the shell expands to v_ej,Sit_ej,Si ∼ 10¹⁴–10¹⁵ cm will be especially bright. This criterion explains the great brightness of models He3.0 and He3.1, which have v_ej,Sit_ej,Si = 0.85 and 0.32 × 10¹⁵ cm, respectively. The amount of mass ejected and the radius reached by that ejecta are not so extreme in the other two models (only to a few times 10¹³ cm), so they are brighter than typical SNe Ib and Ic only during the first few days after the explosion. At their second bright peaks, models He3.0 and He3.1 had peak luminosities of 10^44.54 and 10^44.88 erg s⁻¹ at days 4.5 and 8.3. The effective temperatures there were not well determined, since they were due to shock interaction, but they were calculated to be 21,000 and 32,000 K. Models He2.5 and He3.2, on the other hand, had weak secondary maxima due to radioactivity of 10^41.90 and 10^42.34 erg s⁻¹ at days 22 and 15. Their effective temperatures then were only 5900 and 6600 K.

Zoom In Zoom Out Reset image size

The mass and energy of the ejected shell and the delay time are determined mostly by the amount of silicon that burns explosively in the flash and not directly by the mass of the star that explodes. Brighter explosions could be generated by 2.5 and 3.2 M_⊙ models and fainter ones by 3.0 and 3.1 M_⊙ models. The four choices here are a very limited sampling meant to illustrate a broad possibility of outcomes. In cases where the silicon flash is weak or absent, these supernovae would resemble the models of Section 8.1. See Woosley (2019) for further details and a comparison to SN 2014ft (De et al. 2018).

For the standard mass-loss rate, the most likely candidates for common SNe Ib and Ic have initial helium star masses between 3.0 and 8.0 M_⊙. These produce presupernova stars with masses of 2.45–5.63 M_⊙. Lower-mass stars experience radius expansion, produce little ⁵⁶Ni, and do not look like common events. The subinterval 3.0–3.2 M_⊙ might make normal SNe Ib and Ic but is complicated by the possibility of strong silicon flashes. More massive stars expand too slowly and make broad, faint events (Section 8.6). Both boundaries, 3.0 and 8.0 M_⊙, are imprecise and depend on the assumed mass-loss rate.

Figure 23 shows the expected light curves for all of the models in this mass range that do not experience a strong silicon flash. The interval of initial helium core masses between 3.0 and 8.0 M_⊙ has been further divided into two subintervals, 3.0–5.0 M_⊙ (2.5 < M_preSN < 3.8) and 5.0–8.0 M_⊙ (3.8 < M_preSN < 5.6), to highlight some systematic variations with mass. More realistic calculations of radiation transport than possible with KEPLER will be needed to tell if, spectroscopically, these two regions correspond to SNe Ib and Ic. From Table 1, one expects that all of the models in this mass range would be Type Ib, since Y_s remains close to 1, but it may not be necessary to remove all of the helium to make a Type Ic. Similar presupernova models result from higher-mass stars with greater mass loss; Table 1 shows that if the mass-loss rate is multiplied by 1.5, stars with presupernova masses over 5 M_⊙ will lose their helium.

Zoom In Zoom Out Reset image size

Figure 24 shows the mass ejected by supernovae in approximately this same mass range. These calculations used the W18 engine. Including every star that exploded, the median mass ejected was 2.21 M_⊙, and the average was 2.97 M_⊙, but for the limited range shown here and relevant to SNe Ib and Ic, the median and mean masses were 2.31 and 2.30 M_⊙. This does not include the low-mass explosions with radius expansion or the high-mass ones that produce unusual light curves (Section 8.6). These values are similar to those for observed SNe Ib and Ic (Drout et al. 2011; Lyman et al. 2016; Prentice et al. 2019). In particular, Prentice et al. (2019) reported observed median and mean ejected masses for Type Ib of 2.0 and 2.20 M_⊙. For Type Ic, the median is 2.2 M_⊙ and the mean is 3.2 M_⊙, though the error bars in all cases are large (of order 1 M_⊙). Figure 24 shows that this may reflect an actual large spread of masses in nature.

Zoom In Zoom Out Reset image size

For a typical explosion energy of 10⁵¹ erg, our masses imply average expansion speeds of order 7000 km s⁻¹, but substantial material moves much faster than this average and dominates the spectrum at early times. The lower panels of Figure 24 show the effective temperature and the speed of the zone at the photosphere at the time of peak light. Low-mass models below 3.5 M_⊙ have not been included because of the lingering effects of radius expansion. The 8.0 M_⊙ model was also dropped because of an anomalously low explosion energy (7 × 10⁵⁰ erg). That is not to say that these explosions do not exist, just that they are atypical. The remaining velocities at peak light lie between 6500 and 10,000 km s⁻¹. Speeds above the photosphere are larger and may affect the spectrum. Generally, though, these predictions agree with what is observed for Types Ib and Ic at peak (e.g., Table 6 of Lyman et al. 2016). Photospheric temperatures are between 5000 and 6000 K, and that also agrees reasonably well with observations (Prentice et al. 2019), though perhaps on the cool side.

Generally speaking, the models also have light-curve shapes—rise times and durations—similar to observations (Figure 25; Table 8) but are too faint at peak, on average, by a factor of about 1.5, even using our most optimistic ⁵⁶Ni yields (0.75 times Ni+Tr+α; Figure 9; Section 3.3). Prentice et al. (2016) gave a median luminosity at peak for common SNe Ib and Ic of 10^42.50 and 10^42.51 erg s⁻¹, respectively. Lyman et al. (2016) gave a similar value for Type Ib but 10^42.6 erg s⁻¹ for Type Ic (our Figure 23; their Figure 3). Generally, only our brightest, most energetic supernovae approach these medians. The discrepancy is worse if the ⁵⁶Ni yield is estimated more conservatively using the approach of Sukhbold et al. (2016). A more recent observational study by Prentice et al. (2019) gives substantially reduced values, 10^42.3±0.2 erg s⁻¹ for the mean of SNe Ib and 10^42.3±0.3 erg s⁻¹ for narrow-lined Type Ic. These values are within reach of our more optimistic (highest-energy, greatest ⁵⁶Ni mass) models but are only for the emission from 4000 to 10000 Å. A larger value is expected when the ultraviolet and infrared emission is included (S. Prentice 2019, private communication). All in all, though, it seems that neutrino-powered explosions as parameterized here have trouble making enough ⁵⁶Ni to explain the light curves of many of the brighter ordinary SNe Ib and Ic reported in the literature. A similar conclusion was reached by Anderson (2019). This flies in the face of current convention, which states that all such supernovae come from binaries, are neutrino-powered, and have radioactive-powered light curves, and warrants some discussion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simplest, traditional way to make a SN I brighter is to increase its ⁵⁶Ni yield. Though numerically convenient, such variations are actually severely limited by physics and other observations. Whatever is done for Type Ib and Ic models must also be done for Type IIp models, which have very similar core structure and explosion characteristics (Section 4.5 and Appendix A). There are limits from iron nucleosynthesis (Section 5), the ⁵⁶Ni produced in SN 1987A (Section 5), the tails of the light curves of SNe IIp (Anderson 2019), and the masses of neutron stars (Section 6.1) that prohibit moving the mass cut in too deep. Thus, other explanations of the bright light curves of SNe Ib and Ic are also explored here: an embedded magnetar and circumstellar interaction.

The required average ⁵⁶Ni production for common SNe Ib and Ic, based upon observations and assuming the luminosity at peak is derived solely from radioactivity and given by Arnett's rule (a questionable assumption; Dessart et al. 2015), is 0.14^+0.04_−0.04 M_⊙ for Type Ib and ${0.16}_{-0.10}^{+0.03}$ M_⊙ for Type Ic, according to Prentice et al. (2016). Lyman et al. (2016) gave slightly larger values, 0.17 ± 0.16 M_⊙ for Type Ib and 0.22 ± 0.16 M_⊙ for Type Ic. Anderson (2019) also gave similar values based on a literature survey, 0.163 M_⊙ for Type Ib and 0.155 M_⊙ for Type Ic. In mild tension with these numbers, recent studies by Prentice et al. (2019) gave mean ⁵⁶Ni production of 0.09 ± 0.06 M_⊙ for Type Ib and 0.11 ± 0.09 M_⊙ for Type Ic, but as noted previously, these values are based on observations in a limited wavelength range and are lower bounds. The larger values for the average ⁵⁶Ni masses (i.e., those other than Prentice et al. 2019) are appreciably greater than computed here for models in the initial mass range 3.0–8 M_⊙, the ones in our study deemed most likely to make these sorts of supernovae (Figure 26). As noted in Section 5 and Table 2, the average ⁵⁶Ni production for the W18 central engine for stars likely to become SNe Ib or Ic ranges from 0.034 M_⊙ for the lower-mass stars, weaker engine, and conservative approach to evaluating ⁵⁶Ni production to 0.13 M_⊙ for the highest-mass, most energetic models evaluated in the most optimistic way. Our best estimate would be between Ni and Ni+Tr/2 in Table 2 for the W18 and S19.8 central engines, which is between 0.034 and 0.097 M_⊙, with the smaller value appropriate for the lower-mass supernovae and weaker engine.

Since a large fraction of the observations are explained by models within our error bar, it is natural to ask if still larger values might be possible. Physically, the ⁵⁶Ni yield is constrained by the amount of matter that attains a peak temperature in excess of about 5 × 10⁹ K and is ejected with a small neutron excess. This is an upper bound because other iron-group isotopes will be produced, and the fraction of helium from the alpha-rich freeze out is not small. A much smaller contribution to ⁵⁶Ni also comes from matter experiencing explosive silicon burning between 4 and 5 × 10⁹ K. A popular ansatz for the explosion temperature, T_exp, is given by (Woosley et al. 2002)

$$\begin{eqnarray}&&\displaystyle \frac{4}{3}\pi {r}^{3}{{faT}}_{\exp }^{4}\approx {E}_{\exp },\end{eqnarray} \tag{ 13 }$$

where E_exp is the energy of the explosion, usually taken to be the final kinetic energy; r is the radius at the time the supernova shock begins to move outward from the zone that eventually achieves T_exp; a is the radiation constant; and f is a factor that accounts for the presence of electron–positron pairs. Often the radius of the given mass shell is evaluated using a presupernova model, assuming that the adjustment before the shock arrives is negligible and f is taken to be 1.

In the present situation, even an approximate result requires some improvements. First, the abundance of pairs is not negligible. For the temperatures of relevance, 4 × 10⁹ K and higher, f ≈ 11/4, similar to its value in the early Big Bang (Peebles 1993). This has been verified for our models using the KEPLER EOS. Second, the settling of mass shells located only a few thousand km out as the explosion develops cannot be neglected. The gravitational acceleration at 4000 km is about 10⁹ cm s⁻², and the time from the presupernova stage to the launch of the shock is typically 0.5–1 s. Even though the layers are not far from hydrostatic equilibrium, a contraction of 35% is typical. Finally, the full energy of the explosion may not have been developed when the shock temperature declines to 5 × 10⁹ K. Since the explosion energy only enters as E^1/3 and most of the explosion energy has usually been deposited, this last factor is not critical.

Still, appropriately corrected, Equation (13) works very well. For an explosion energy of 1.2 × 10⁵¹ erg and f = 2.75, the equation says that a temperature greater than 5 × 10⁹ K is achieved inside a radius of 2800 km. For f = 1, the boundary radius is increased to 3900 km. Examining the runs themselves, 90% of ⁵⁶Ni is typically synthesized inside a radius that, if the structure is evaluated when the shock is launched, has a value of 2970 ± 9 km. The same 90% criterion gives a radius of 4470 ± 5 km using radii in the presupernova star. Of course, it is not realistic to neglect pairs, but it turns out that the correction to the radius for pairs, (11/4)^1/3 ≈ 1.40, very nearly cancels the contraction factor, 4470/2970 ≈ 1.50. Thus, Equation (13) can be used approximately with either f = 1 and r (the presupernova radius) or f = 2.75 and r (the radius in the star when the shock is launched). This Lagrangian mass shell in the presupernova star where the shock temperature declines below 5 × 10⁹ K, the value required to achieve NSE on a hydrodynamic timescale, will be referred to as M_T9=5.

It is important to note that M_T9=5 is mostly a function of presupernova structure. It depends only weakly on explosion energy and not on the mass cut. The mass of the ejecta that is heated above 5 × 10⁹ K is then just M_T9=5 − M_rem, where M_rem is the baryonic mass of the gravitationally bound remnant. Figure 27 shows the good agreement of this estimate with the mass of ejected material that actually achieved NSE in those models likely to make SNe Ib or Ic. The agreement suggests that the ⁵⁶Ni yield is not sensitive to details of how the explosion is calculated but mostly to the final mass cut and, secondarily, the explosion time.

Zoom In Zoom Out Reset image size

This semianalytic model allows a qualitative exploration of how ⁵⁶Ni synthesis might change as the explosion energy and mass cut are varied for a given presupernova star. Figure 28 shows the results if the explosion energy is varied between 1 and 5 × 10⁵¹ erg or the mass cut is moved in to the edge of the iron core. These are extreme excursions. None of our models had explosion energies exceeding 2.0 × 10⁵¹ erg, even with the energetic S19.8 engine, and most had less (Figure 5). No modern calculation of a purely neutrino-powered model that we know of has greatly exceeded that (e.g., Fryer 1999; Fryer & Kalogera 2001; Scheck et al. 2006; Bruenn et al. 2016; Müller et al. 2017). Similarly, the mass cut in neutrino-powered models rarely, if ever, occurs as deep as the edge of the iron core. The reduction in remnant mass might also cause an unacceptable decrease in the mean neutron star mass (Section 6.1) and could make the contribution of SNe Ia to iron synthesis a minor component. Nevertheless, neglecting these constraints, the synthesis of iron group plus alphas might be pushed to ∼0.20 M_⊙. For reasonable explosion energies, the placement of the mass cut has more leverage on the yield than the explosion energy. The maximum ⁵⁶Ni mass, even in models that have unreasonable mass cuts, is 0.33 M_⊙. Some past studies, e.g., Dessart et al. (2012, 2015), arbitrarily placed the mass cut at the edge of the iron core and thus overestimated the ⁵⁶Ni yield for purely neutrino-powered explosions.

Zoom In Zoom Out Reset image size

In summary, larger ⁵⁶Ni yields capable of resolving the discrepancy between the models and the observed brightness of SNe Ib and Ic might conceivably be possible, and certainly can be prepared in parameterized explosions, but they are not achievable in today's neutrino-powered models. The yields may violate limits on iron nucleosynthesis, and the neutron stars may be too small compared to observations. The same approximations applied to SN 1987A would overproduce ⁵⁶Ni there. An additional source of energy seems to be required to explain supernovae brighter than about 10^42.5 erg s⁻¹.

8.5.1. Circumstellar Interaction

The brightness of the supernova could be amplified by colliding with a circumstellar gas, presumably the outcome of mass loss during the final year of the star's life (e.g., Shiode & Quataert 2014; Woosley 2017). On the positive side, ample energy exists in the outer edge of the ejecta, and the required circumstellar mass (CSM) is not large. In a typical Ib progenitor, e.g., model He6.0, the outer 0.001 M_⊙ moves at speeds in excess of 0.1c, implying a kinetic energy in this small mass of ∼10⁴⁹ erg, more than enough to explain the full light curve of a typical SN Ib. This would only need to encounter ∼0.001 M_⊙ of CSM matter inside a radius of a few times 10¹⁵ cm.

On the other hand, getting the time structure of the light curve is problematic and requires the introduction of many uncertain parameters (Chatzopoulos et al. 2013; Vreeswijk et al. 2017). The spectrum from such high-velocity material might not resemble common SNe Ib and Ic. In the simplest case of a steady wind in which ρ ∝ r⁻², the luminosity

$$\begin{eqnarray}&&{L}_{\mathrm{csm}}\approx 0.5\dot{M}{v}_{\mathrm{shock}}^{3}/{v}_{\mathrm{wind}},\end{eqnarray} \tag{ 14 }$$

where $\dot{M}$ is the mass-loss rate, v_wind is the wind speed, and v_shock is the speed of the fast-moving supernova ejecta. For $\dot{M}$ = 2 × 10⁻⁴ M_⊙ yr⁻¹, v_wind = 500 km s⁻¹, and v_shock = 30,000 km s⁻¹, the luminosity is 3 × 10⁴² erg s⁻¹. This luminosity would persist and slowly decline until the fast-moving ejecta encountered about half its mass, or ∼20 days. But the light curve would commence immediately and decline slowly, quite unlike what is observed.

Assuming a CSM shell with constant density only improves things marginally. Now the luminosity

$$\begin{eqnarray}&&{L}_{\mathrm{csm}}\approx 2\pi {r}^{2}\rho {v}_{\mathrm{shock}}^{3}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\approx \ 2\pi \rho {v}_{\mathrm{shock}}^{5}{t}^{2}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\approx \ 4\times {10}^{42}{\left(\displaystyle \frac{{v}_{\mathrm{shock}}}{3\times {10}^{9}}\right)}^{5}\left(\displaystyle \frac{\rho }{{10}^{-17}}\right)\ \mathrm{erg}\,{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 17 }$$

increases quadratically with time and declines once the edge of the shell is reached. A density structure, ρ(r), could be found to match any given light curve (Chatzopoulos et al. 2013), but its construction seems artificial. Given the lack of narrow CSM lines in the spectra of common SNe Ib, CSM interaction seems an unlikely estimate for their luminosity excess.

8.5.2. A Magnetar

A significant fraction, substantially more than 10%, of all neutron stars are born with magnetic fields (0.3–1) × 10¹⁵ G (Beniamini et al. 2019). This is roughly consistent with the fraction of core-collapse supernovae that are of Types Ib and Ic. It is thus possible that a substantial fraction of this class of supernovae makes magnetars. In one of the few studies of Type Ib progenitors that included rotation, Yoon et al. (2010) found periods of ∼10–20 ms, rapid enough for rotational energy to affect the light curve of a supernova, though not its explosion energy. Loss of the hydrogen envelope diminishes the rotational braking the helium core would have experienced in a red giant, though angular momentum is still lost to winds. The explosion would still be neutrino-powered as calculated here, but the light curve would be dominated by the magnetar.

As an example, consider model He6.0. The unmodified version of this supernova, using the S19.8 central engine, exploded with a kinetic energy of 1.2 × 10⁵¹ erg and produced 0.083 M_⊙ of ⁵⁶Ni. Figure 29 shows the effect on the light curve of embedding a magnetar with field strength 5, 7, and 10 × 10¹⁴ G and initial rotational energy 2.5, 3.0, and 3.5 × 10⁴⁹ erg. The prescription for energy deposition here was the same as that of Woosley (2010), which is very similar to that of Kasen & Bildsten (2010) but uses slightly different field strengths to achieve the same luminosity. Energy from radioactive decay was included in the three magnetar models.

Zoom In Zoom Out Reset image size

The corresponding initial rotational periods are 28, 26, and 24 ms. Other parameters could be chosen to give different peak luminosities and slopes on the tail, but these three peaked near 10^42.7 erg s⁻¹, similar to a brighter-than-average SN Ib or Ic (Prentice et al. 2016). This compares with 10^42.21 erg s⁻¹ for the same model with no magnetar. The rise time from 50% to peak luminosity (t_−1/2) for the three models varied from 10.4 to 12.4 days, with the models with greater rotational energy but weaker magnetic fields rising more slowly. The decline time from peak to 50% of peak luminosity (t_+1/2) was 19.0–28.5 days for the same models. For the model without the magnetar, the rise and decline timescales were 10.4 and 23.4 days. The ejected mass here was 2.816 M_⊙, and the neutron star remnant was 1.417 M_⊙ (gravitational). The radioactive contribution came from the decay of 0.084 M_⊙ of ⁵⁶Ni.

The idea that many "normal" SNe Ib and Ic might not be powered at peak by radioactive decay is novel and will have important implications for the role of magnetars in many diverse kinds of supernovae. Since they are much more common than superluminous supernovae and a more regular class of events, they might be a good laboratory for the study of magnetar birth. After 100 yr, the magnetars in Figure 29 would have periods of 1.26, 0.88, and 0.63 s and luminosities of 4 × 10³⁶, 8 × 10³⁶, and 1.6 × 10³⁷ erg s⁻¹, assuming the field strength remains constant.

The tails of the light curves of SNe Ib and Ic might be studied carefully to ascertain ⁵⁶Ni abundances, but this could be complicated by the fact that gamma rays are already escaping shortly after peak in the radioactive model. For some field strengths, the tail of the magnetar models also declines at a similar rate as the radioactive model (Figure 29). On the other hand, light curves declining slower than the 77.3 day half-life of ⁵⁶Co, i.e., full gamma-ray trapping, might be indicative of a magnetar energy source. If the iron abundance could be accurately inferred from late-time spectra, the magnetar model predicts a much smaller abundance than needed to explain the brightness of the supernova.

Helium stars with presupernova masses greater than 6 M_⊙ require extraordinary explosion energies or gamma-ray escape to avoid producing light curves that are too faint and broad to be common SNe Ib and Ic (Ensman & Woosley 1988), yet our present formalism predicts many such explosions (Table 5; Figure 3). This is particularly true when the more energetic S19.8 central engine is employed. Figure 30 shows a selection of light curves. Though some are as bright as their lower-mass counterparts in Figure 23, the curves all peak appreciably later and decline much more slowly. The decline time to 50% peak luminosity (t_+1/2) takes more than 40 days for the lightest models in this group and longer for the more massive ones (Table 5). The three broadest, faintest light curves in Figure 30 came from helium stars with initial masses of 19.25, 19.50, and 19.75 M_⊙ that had presupernova masses from 10.07 to 10.29 M_⊙. These supernovae took 2 months to decline to half their peak luminosity. Their effective temperatures at peak were near 3500 K, so they would be quite red. Since they lost most of their helium prior to dying, these supernovae would probably be Type Ic.

Zoom In Zoom Out Reset image size

We are not aware of any observed supernovae with these characteristics. Perhaps they have been missed because they are faint, slow, and red, or these massive progenitors, essentially all of the stripped stars with presupernova masses greater than about 6 M_⊙, collapsed to black holes. This is an issue that gravitational radiation studies might some day address. Are there indications in the black hole birth function of cores with masses greater than 6 M_⊙ exploding and not leaving black hole remnants?

Or, our models may be deficient. It is hard to see how anything other than a major alteration of the neutrino-powered, radioactive-illuminated paradigm could greatly increase the luminosities or decrease the light curves' widths for such massive progenitors. The alterations again include circumstellar interaction or extra energy from a rotationally powered central engine.

In favor of a magnetar explanation is the fact that, faint or not, many of these massive models do explode (Table 5). A successful neutrino-powered explosion during the first second after the iron core collapses facilitates the subsequent formation and operation of a magnetar by evacuating a region around the PNS and allowing it time to cool, contract, and spin up. Since the initial explosion makes only ≲0.15 M_⊙ of ⁵⁶Ni and the magnetar would make none, the luminosity both at peak and, usually, on the tail would not be due to radioactivity (see, e.g., Woosley 2010), but the iron lines might be strong in the remnant. While the calculation of rotating presupernova models is deferred, we note the tendency of more massive stars, especially those that lose their envelopes early on, to produce rapidly rotating neutron stars (Heger et al. 2005). Most of the angular momentum loss in the core occurs, in single stars, during helium burning, when the star's hydrogen envelope expands to become a giant. The essentially stationary envelope saps angular momentum from the rapidly rotating core, so stars stripped of their envelopes, as these are, end up rotating more rapidly at death.

There are two sorts of magnetar models: those where the magnetar's initial rotational energy is small compared with the kinetic energy of the neutrino-powered explosion and those where it is not. In the former case, the magnetar just illuminates a coasting configuration whose dynamics were determined by the neutrinos. These are the magnetars that were studied in Section 8.5.2. In the latter case, energy input by the magnetar appreciably modifies the velocity profile of the supernova. These two cases are illustrated with models He10.00 (Figure 31) and He19.75 (Figure 32). Both used the W18 explosion model (Table 5) for their initial explosions. Model He10.00, which had a presupernova mass of 6.74 M_⊙, is a bit heavier than usually deemed responsible for common SNe Ibc. Model He19.75, with a presupernova mass of 10.3 M_⊙, was in the island of high-mass explodability in Figure 13 and one of the heaviest models to explode in the current study.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The five magnetars introduced in model He10.00 all had magnetic fields of 5 × 10¹⁴ G and rotational energies of 0.1, 0.2, 0.5, 1.0, and 2 × 10⁵¹ erg. These rotational energies were chosen to span a range where the rotational energy was either a small or major adjustment to the energy of the neutrino-powered supernova, 7.6 × 10⁵⁰ erg. Given the single value of the magnetic field considered, the light curves are all very similar on their tails (Figure 31) but vary greatly in peak luminosity and rise time. All are substantially brighter than the purely radioactive case. For the smaller rotational energies, the light curves are similar to those for the lower-mass models (Figure 29) but broader due to the larger mass. They could potentially account for some of the brighter, broader SNe Ib and Ic in Figure 26. For larger magnetic fields, the light curve would decline more rapidly. The light curves in Figure 31 assume the full trapping of the magnetar-deposited energy. Depending on the spectrum of the magnetar and the pulsar wind interaction inside the supernova, a substantial fraction of the energy might escape at late times and not contribute to the optical light curve. The column depths at 200 days for the five magnetar energies shown are 28, 22, 14, 9, and 6 g cm⁻². If the opacity were like gamma rays from radioactivity, ∼0.05 cm² g⁻¹, then the light curve at 200 days for the most energetic case would be about three times fainter. The optical depth at earlier times would scale approximately as t⁻².

The light curves for the more energetic cases bear some similarity to SNe Ic-BL. A potential problem with such a scenario, though, is that, except for a small mass near the surface (Figure 31), the velocities are quite low, hence no "broad" lines. Typical SNe Ic-BL are estimated to have velocities characteristic of 4 M_⊙ exploding with 7 × 10⁵¹ erg with ⁵⁶Ni masses near 0.3 M_⊙ (Taddia et al. 2019), or over 10,000 km s⁻¹. Drout et al. (2011) adopted a velocity of 20,000 km s⁻¹ at peak luminosity. Both this extreme energy and mass of ⁵⁶Ni are well out of reach of our standard models. There are three possible solutions to this difficulty. (1) The masses of SNe Ic-BL may be lighter or their explosion energies greater than considered here. A different explosion mechanism would need to be invoked to raise the energy. (2) The light at peak could come from circumstellar interaction. The outer 0.015 M_⊙ of model He10.00 moves faster than 12,000 km s⁻¹ and has kinetic energy 3 × 10⁴⁹ erg, more than enough to provide the entire light curve of most SNe Ic-BL, if the matter interacted with a comparable mass of CSM at a few × 10¹⁵ cm. The problems with this scenario are explaining why the necessary CSM is there with the right density distribution to give the observed light-curve shape (Section 8.5.1), and why narrow lines from the CSM are not usually seen in SNe Ic-BL. This may be a more promising model for SNe Ibn (Pastorello et al. 2016; though typical "broad-line" components are slower for this class). (3) The light curve and most of the explosion energy could come from a magnetar or accretion into a black hole. The explosion would then probably be asymmetric.

Also shown in Figure 31 are the velocity profiles resulting from the various magnetars. If plotted versus radius instead of mass, the regions of near-constant coasting speed would show up as thin, dense shells. For large magnetar energy and 1D simulations, most of the mass of the explosion is piled up in this thin shell. That the shell is thin and dense is an artifact. In a multidimensional simulation, the hot pulsar wind for the more energetic magnetars would break through the slower-moving denser material outside (Chen et al. 2016) and surround the supernova with an envelope of faster-moving material. The explosion would be asymmetric. The problem here is that a definitive calculation of the spectrum is difficult. At late times, the typical velocities would be slower, more like the broad flat spots in Figure 31. Velocities of 4000–5000 km s⁻¹ are typical of some long duration Type Ic super luminous supernova (SLSN; Jerkstrand et al. 2017) and may be similar to what is seen in some SNe Ic-BL. Maurer et al. (2010) found that the velocities for SNe Ic-BL are highly variable and might depend on observing angle. They derived an average of about 6000 km s⁻¹.

As a second example, consider model 19.75. Without a magnetar, this model exploded, using the W18 central engine, with a terminal kinetic energy 8.4 × 10⁵⁰ erg and ejected 0.088 M_⊙ of ⁵⁶Ni. The presupernova mass was 10.28 M_⊙, of which 8.5 M_⊙ was ejected. Only 0.25 M_⊙ of that was helium. This would have been a SN Ic. The original ZAMS mass of the star before it lost its envelope in the binary would have been 51 M_⊙. This explosion was recomputed with an embedded magnetar with an initial rotational energy of 3.8 × 10⁵¹ erg (rotational period 2.3 ms) and magnetic field strengths of 1, 2, and 5 × 10¹⁴ G. Energy deposition was calculated according to Woosley (2010), and full absorption was assumed. The case with the smallest field strength is similar to magnetar parameters used by Vreeswijk et al. (2017) to fit the Type I SLSN iPTF12dam and is rather typical of the large spread observed in SLSN (Nicholl et al. 2017; Villar et al. 2018). Figure 32 shows the light curves both with and without the magnetar. With the magnetar, the model essentially replicates the results of Vreeswijk et al. (2017) and observations of iPTF12dam. As noted by those authors, the magnetar model, even with some approximation for escape, is too bright at very late times compared with the observations, but magnetar energy transport and breakout are poorly understood. The improvement here over Vreeswijk et al. (2017) is that the explosion is derived from an actual stellar evolution and explosion calculation.

The light curve shown in Figure 32 might also rise more rapidly than would actually be observed in the optical because the magnetar pushes the ejected matter into a thin shell that has an unphysically short diffusion time. The bolometric light curve also peaks earlier than the optical light curve because of its high effective temperature (Nicholl et al. 2017).