A High-resolution Study of Presupernova Core Structure

As a massive star below the pair-instability threshold (∼80 M_⊙; Woosley 2017) evolves through its final stages of nuclear burning, its central regions cool by neutrino emission, become degenerate, and tend to decouple from the overlying layers and evolve as separate stars. Although never becoming completely detached except for the lowest-mass stars, the presupernova core takes on a structure similar to that of a white dwarf, mostly composed of iron, surrounded by a dense mantle of oxygen and intermediate-mass elements. As numerous studies have shown, the structure of this configuration, and especially the rate at which the density declines outside the iron core, is strongly correlated with the difficulty of blowing the star up (e.g., Burrows & Lattimer 1987; Fryer 1999). Recent studies have sought to capture this complex structure in just one or two parameters that might predict, albeit approximately, whether a star with a given mass blows up or collapses to a black hole simply from looking at one-dimensional models for stellar evolution (O'Connor & Ott 2011, 2013; Ugliano et al. 2012; Pejcha & Thompson 2015; Ertl et al. 2016; Müller et al. 2016; Sukhbold et al. 2016).

Sukhbold & Woosley (2014) discussed how the advanced burning stages, especially convective carbon and oxygen burning, sculpt this structure and provided a library of presupernova stars consisting of 503 models in the mass range 12–65 M_⊙ to demonstrate the systematics and its dependence on input physics. Müller et al. (2016) prepared a larger grid of over 2000 models between 10 and 32.5 M_⊙ and corrected some errors in Sukhbold & Woosley (2014). Their results showed finer structure but similar global systematics (see especially their Figure 6).

More recently, Farmer et al. (2016), using the MESA code, found that the choice of mass resolution, i.e., zoning, "dominates the variations in the structure of the intermediate convection zone and secondary convection zone during core and shell hydrogen burning, respectively" and greatly affects the structure of presupernova stars. They found that a minimum mass resolution of ∼0.01 M_⊙ was necessary to achieve convergence in the final helium core mass at the ∼5% level. They also found ∼30% variations in the central electron fraction and mass locations of the main nuclear-burning shells, and that a minimum of ∼127 isotopes was needed to attain convergence of these values at the ∼10% level.

Renzo et al. (2017) also used the MESA code to explore the sensitivity of massive star evolution due to variations in the mass loss rate for stars with initial masses between 15 and 35 M_⊙. They found variations in the presupernova core-compactness parameter (Section 4.3) of ∼30%, depending upon the choice of the algorithm. In a limited study of resolution in one model, they found roughly 9% variation in final core compactness, smaller than that resulting from uncertainties in the mass loss prescription.

In this paper, we present a new survey similar to Sukhbold & Woosley (2014) and Müller et al. (2016) but using much finer zoning, a greater number of models, and several different mass loss rates. The key differences of the new survey are listed in Section 2, including the correction to an erroneous pair-neutrino loss rate used by Sukhbold & Woosley (2014). Section 2 also offers a detailed discussion of the effects due to zoning, network size, and boundary pressure. In general, we find that the pattern of final core compactness seen in these previous studies is unaltered by finer zoning or a larger reaction network, though the mass limits for different behaviors are shifted by about 10% when the pair-neutrino loss rate is corrected. In Sections 3 and 4, we discuss the general characteristics of the new survey, including its observable properties. With a much greater number of models, clear evidence emerges for multivalued solutions to the presupernova structure of stars in the 14–19 and 22–24 M_⊙ ranges (see also Müller et al. 2016). In Section 5, we offer an interpretation for this feature through the physics of the advanced-stage evolution in the cores of massive stars. Large variations in core structure are expected for stars of these masses no matter what resolution, mass loss rate, or reaction network is employed in the calculation (for the effects of such variations in individual models, see Farmer et al. 2016; Renzo et al. 2017). Future surveys of supernova, as well as presupernova, models should thus take care to examine, when feasible, a large number of masses in order to sample an outcome that is essentially statistical in nature. More specific results and effects on explodability and remnant masses are discussed in Section 6. Finally, in Section 7, we offer our conclusions.

With several important exceptions, the code, physics, and input parameters used here are the same as in Sukhbold & Woosley (2014) and Müller et al. (2016). The (solar) initial composition is from Asplund et al. (2009), as was used in Müller et al. (2016), but with an appropriate correction for the ratio ¹⁵N/¹⁴N (Meibom et al. 2007). Convective and semiconvective settings, nuclear reaction rates, opacities, and mass loss rates are the same as in both of the earlier works. The mass range studied and lack of rotation are the same as well.

Nuclear burning was handled as usual, using a 19 isotope approximation network until the central oxygen mass fraction declines below 0.03 and a silicon quasi-equilibrium network with 121 isotopes (Weaver et al. 1978) thereafter. Numerous studies that compare this treatment with the energy generation and bulk nucleosynthesis obtained by "coprocessing" with a network of several hundred to over 1000 isotopes have shown excellent agreement (e.g., Woosley et al. 2002; Woosley & Heger 2007), at least up to central oxygen depletion, where the switch to the quasi-equilibrium network is made. By coprocessing, we mean carrying a large network in each zone using the same time step, temperature, and density but not coupling the energy generation from that large network directly to the iterative loop within a time step. After oxygen depletion, the quasi-equilibrium network is more stable, contains the same weak interaction physics, and is roughly an order of magnitude faster. An exception is the treatment of oxygen and silicon burning in stars lighter than about 11 M_⊙ (Woosley & Heger 2015). For such light stars, it is important to follow neutronization in multiple off-center shells where the quasi-equilibrium approximation has questionable validity and can be unstable. Such light stars are not part of the present survey, which starts at 12 M_⊙. Section 2.2 discusses the issue of network sensitivity in greater detail and confirms the validity of the approximation network plus quasi-equilibrium for representative cases.

Key differences in the new study are as follows.

• 1.  
  The models in this paper employ much finer mass resolution (Figures 1, 2, and 4). For the lightest stars considered (12 M_⊙), the increase is roughly a factor of four. Since an effort was made to keep fine resolution even in the larger stars, the factor for the highest masses considered, around 60 M_⊙, was closer to 15. In the most massive models, the total number of mass shells was approximately 16,000. Additional studies of individual cases of zoning sensitivity in 15 and 25 M_⊙ models (Section 2.1) showed little systematic difference. In all cases, zones contained less than about 0.01 M_⊙ everywhere but were about 10 times smaller than that in the heavy element core (Figure 4). It is this core that is most critical in determining the final presupernova core properties, such as compactness. Surface zoning was a few times 10⁻⁵ M_⊙, and all temperature and density gradients were well resolved (Figure 2).
• 2.  
  A major, long-standing coding error in the KEPLER code was repaired. The error affected the axial vector component of the pair-neutrino losses. In particular, the error resulted in the accidental zeroing of the second term involving ${C}_{{\rm{A}}}^{{\prime} 2}$ and ${Q}_{\mathrm{pair}}^{-}$ in the expression (Itoh et al. 1996, Equation (2.1))

  $$\begin{eqnarray}&&{Q}_{\mathrm{pair}}=\displaystyle \frac{1}{2}[({C}_{V}^{2}+{C}_{A}^{2})+n({C}_{V}^{{\prime} 2}+{C}_{A}^{{\prime} 2})]{Q}_{\mathrm{pair}}^{+}\end{eqnarray} \tag{ 1 }$$

  $$\begin{eqnarray}&&+\displaystyle \frac{1}{2}[({C}_{V}^{2}-{C}_{A}^{2})+n({C}_{V}^{{\prime} 2}-{C}_{A}^{{\prime} 2})]{Q}_{\mathrm{pair}}^{-}.\end{eqnarray} \tag{ 2 }$$

  See Itoh et al. (1996) for the definitions of the quantities. The consequence of this error was that the pair-neutrino loss rate was underestimated by close to a factor of two during carbon burning and somewhat less during later, higher-temperature stages. The error, introduced through inadequate checking of a routine provided by email, was included in 2001 and affected all KEPLER calculations published through 2014. In particular, it affected the often-cited calculations of Woosley et al. (2002), Woosley & Heger (2007), and Sukhbold & Woosley (2014). Because some of the models from Sukhbold & Woosley (2014) were used by Sukhbold et al. (2016), it also affected the outcome of that work for masses larger than 14 M_⊙. The bug was repaired, however, in the works of Woosley & Heger (2015) on 9–11 M_⊙ stars and Woosley (2017) on pulsational pair instability supernovae. Most relevant to this present work, the bug was also repaired for the work of Müller et al. (2016). Because of the strong sensitivity of neutrino and nuclear reaction rates to temperature, the effect of the bug was a slight shift upward in the burning temperature for the late stages of stellar evolution. For a given model, the change in presupernova structure was not great and, as we shall see, within the "noise" of other uncertainties. For the lightest models (<14 M_⊙), the difference is hardly noticeable, but it did systematically shift the outcome for higher-mass stars significantly downward. For example, the peak in compactness that occurred for Sukhbold & Woosley (2014) at about 24 M_⊙ is shifted downward in the work of Müller et al. (2016) and in the present work (Section 4) by about 3 M_⊙.
• 3.  
  The surface boundary pressure is much less than in Müller et al. (2016), which significantly affects the final red supergiant (RSG) properties but does not significantly alter the core structure. Whereas the mass loss rates are varied in the present study, the new calculations do not include Wolf–Rayet models where the envelope is completely lost. All stars studied here retained a substantial hydrogen envelope when they died.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.1. The Effects of Zoning

2.2. The Effects of Networks

Uncertainty in nuclear energy generation and neutronization during the various episodes of burning can be another source of variability in presupernova models. In this section, the effect of using either the standard 19 isotope approximation network plus silicon quasi-equilibrium (Section 2) or a much larger nuclear reaction network is explored. It should be noted that while the standard network carries the abundances of only 19 species, it has the power of a network roughly twice that size, since reactions through trace species like ²³Na, ²⁷Al, etc. are carried in a steady-state approximation (Weaver et al. 1978). Reactions coupling ⁵⁴Fe to ⁵⁶Ni are also included, as are approximations to energy generation by the CNO and pp cycles and to ¹⁴N(α, γ)¹⁸F(${e}^{+}\nu {)}^{18}{\rm{O}}$(α, γ)²²Ne.⁷

The large network used for comparison here is drawn from an extensive nuclear database maintained and used for nucleosynthesis studies for decades (e.g., Woosley et al. 2002). Key reaction rates, such as ¹⁴N(p, γ)¹⁵O, 3α, ¹²C(α, γ)¹⁶O, the (α, γ) and (α, p) reactions on α-nuclei, and the weak interaction rates are the same in both networks. The network is "adaptive" (Rauscher et al. 2002) in the sense that isotopes are added or subtracted at each step to accommodate the reaction flow (e.g., there is no need to include a detailed network for the iron group during hydrogen burning). Usually, the settings on the network size are quite liberal, resulting in from 700 to 1200 isotopes being carried in the early and late stages of evolution, respectively (Woosley & Heger 2007). Here, because we are only interested in stellar structure and not, e.g., the s-process, a smaller number was carried, typically 300 isotopes extending up through the element selenium. All of the iron-group isotopes in the quasi-equilibrium network were included, and many more, so that core neutronization during and after silicon burning was equivalently calculated.

Network sensitivity was examined for four different masses of stars, 15.00, 15.01, 25.00, and 25.01 M_⊙. Two studies were carried out to test network sensitivity. Table 2 edits the evolution of the 15.00 star, similar to Model 15.00D in Table 1, at different times in its evolution using three approximations to the energy generation and neutronization. The Approx case corresponds to using the 19 isotope approximation network and quasi-equilibrium network, as is standard in the rest of this paper. The Coproc case carries the large network of about 300 isotopes along with the approximation network in "coprocessing mode." That is, the approximation and quasi-equilibrium networks are used for energy generation in the stellar structure calculation, but the network is also carried along in passive mode. The same time step is used at the same temperature and density to evolve, zone by zone, a much larger number of isotopes in parallel with the approximation network. Subcycling is done within a given time step to follow the abundances of trace isotopes accurately. Coprocessing is only performed in the code up to the point where the transition to quasi-equilibrium occurs in a zone, typically at oxygen depletion X(¹⁶O = 0.03). The output from the large coprocessing network is used to continually update the electron mole number, Y_e, which feeds back into the structure calculation. Because Y_e evolves slowly, no iteration is required. The energy generation from the large network can also be used to check the validity of the approximation network, though no correction is applied. In the third case, Full, the big network is coupled directly to the structure calculation throughout the entire life of the star, including silicon burning and core collapse. This might be deemed the most accurate approach, but it is slow, requiring an order of magnitude more computer resources, and prone to instability during silicon burning because of strongly coupled flows.

Table 2.  Central Properties for Different Network Treatments for the 15.00 M_⊙ Model

Network	$\mathrm{log}{T}_{{\rm{c}}}$	$\mathrm{log}{\rho }_{{\rm{c}}}$	${Y}_{{\rm{e}},{\rm{c}}}$	$\mathrm{log}{S}_{{\rm{n}}}$	$\mathrm{log}| {S}_{\nu }|$
	(K)	(g cm−3)		(erg g−1 s−1)
H burning
Approx	7.564	0.822	0.70042	5.104	3.910
Coproc	7.564	0.822	0.70196	5.104	3.910
			0.70196	5.107
Full	7.565	0.820	0.69991	5.109	3.914
H depletion
Approx	7.676	1.138	0.50581	5.376	4.114
Coproc	7.676	1.138	0.50616	5.376	4.114
			0.50616	5.469
Full	7.677	1.136	0.50489	5.387	4.129
He burning
Approx	8.249	3.127	0.49998	5.600	1.375
Coproc	8.250	3.123	0.49932	5.621	1.383
			0.49932	5.631
Full	8.249	3.122	0.49932	5.621	1.371
He depletion
Approx	8.406	3.426	0.49998	5.624	2.560
Coproc	8.406	3.426	0.49932	5.627	2.562
			0.49932	5.667
Full	8.407	3.422	0.49932	5.667	2.577
C ignition
Approx	8.703	4.544	0.49998	0.405	4.750
Coproc	8.704	4.546	0.49933	0.434	4.756
			0.49933	0.368
Full	8.700	4.513	0.49933	0.159	4.729
C depletion
Approx	9.081	6.602	0.49998	6.193	8.198
Coproc	9.080	6.598	0.49894	6.181	8.202
			0.49894	6.108
Full	9.079	6.580	0.49897	6.063	8.212
O ignition
Approx	9.176	6.895	0.49998	10.63	9.123
Coproc	9.176	6.899	0.49896	10.60	9.116
			0.49896	10.58
Full	9.176	6.944	0.49899	10.57	9.188
O depletion
Approx	9.343	7.056	0.49998	11.53	11.05
Coproc	9.346	7.037	0.49242	11.60	11.12
			0.49242	11.64
Full	9.352	6.990	0.49309	11.81	11.27
Si ignition
Approx	9.477	8.068	0.48679	12.14	11.21
Coproc	9.477	8.087	0.48609	12.18	11.19
Full	9.477	8.084	0.47945	11.51	11.29
Si depletion
Approx	9.583	7.699	0.47060	12.71	13.01
Coproc	9.571	7.706	0.46987	12.67	12.97
Full	9.477	7.739	0.46809	12.40	12.92
Presupernova
Approx	9.885	9.747	0.43848	(16.11)	16.24
Coproc	9.880	7.762	0.43747	(16.09)	16.20
Full	9.869	10.10	0.43266	(16.83)	16.94

Note. Here S_n is the energy generation rate excluding neutrino losses and is negative during the presupernova stage due to photodisintegration. Values in parentheses denote $\mathrm{log}| {S}_{{\rm{n}}}|$. The two values of Y_e,c and S_n given for coprocessing calculations until core oxygen depletion (until quasi-equilibrium) are from the approximation network (top) and coprocessed big network (bottom).

Download table as:  ASCIITypeset image

In Table 2, up until oxygen depletion, two different numbers are given for the coprocessing run (Coproc). The first is the energy generation from the approximation network, which is used to calculate the stellar model, and the second is the output from the passively carried big network. Hydrogen burning and depletion correspond to central hydrogen mass fractions of 0.4 and 0.01. Helium burning and depletion are when the central helium mass fraction is 0.5 and 0.01. Carbon "ignition" is actually evaluated during the Kelvin–Helmholtz contraction between helium depletion and real carbon burning when the central temperature is 5 × 10⁸ K, and carbon depletion is when the central carbon mass fraction is 1%. Similarly, "oxygen ignition" is when the central temperature during the Kelvin–Helmholtz contraction is 1.5 × 10⁹ K and oxygen has yet to burn, and oxygen depletion is when the central oxygen mass fraction is 0.03. Silicon ignition is at 3.0 × 10⁹ K, and silicon depletion is when the silicon mass fraction is 1%. Presupernova is when the core-collapse speed exceeds a maximum of 900 km s⁻¹.

The table shows near-exact agreement in energy generation and Y_e calculated using all three approaches up until at least oxygen ignition. The slight differences in energy generation at carbon ignition do not matter because neutrino losses (far right column) dominate at the time examined. By oxygen depletion, though, things are starting to mildly diverge. This divergence has three causes. One is the electron capture that goes on in the late stages of oxygen burning and starts to appreciably affect Y_e. This change is not followed by the Approx network. Carrying Y_e in the coprocessing run addresses most of this divergence; the values of Y_e in the coprocessing and full network runs are nearly the same. Another effect, more difficult to disentangle, is the fact that, at late times, the stellar structure in the different calculations starts to diverge. It diverges, in part, because of the different nuclear physics, but even more because of the chaotic nature of the evolution. As was seen in the zoning study (Section 2.1), in certain initial mass ranges, two models with even slightly different physics will differ appreciably in final outcome. The third effect is the instantaneous difference due solely to the differing treatments of nuclear physics. This is what we are actually trying to study, but it is difficult to separate from cumulative effects due to structural changes. Even so, at silicon depletion, the three values of Y_e agree very well, and the energy generation differs by a factor of two, which is mostly due to the higher temperature in the Approx case and the very sensitive temperature dependence of the silicon-burning reactions.

These small differences probably cause very little change in the presupernova model because the star has a certain amount of fuel to burn and the total energy release is set by known nuclear binding energies. The different burning rates only modestly affect the temperature at which the silicon burns in steady state (Woosley et al. 2002). We conclude that the approximation–plus–quasi-equilibrium approach is adequate for our survey and introduces errors that are small compared with other variations that occur when other uncertain quantities, or even the zoning, change.

Table 3 examines the effect of using either the large network or the approximation network plus quasi-equilibrium to study stars of 15.00, 15.01, 25.00, and 25.01 M_⊙, the same masses considered in the zoning sensitivity study (Section 2.1). Only one zoning was considered here, but the resolution was about three times greater in all four cases than in Sukhbold & Woosley (2014). Fewer time steps were taken during hydrogen burning than in the zoning sensitivity study. An overall limit of 5 × 10¹⁰ s was placed on the time step, as opposed to 1 × 10¹⁰ s in the survey and zoning sensitivity study. Here the emphasis is on comparing runs with different networks, not on the highest spatial and temporal resolution. Consequently, the 15.00 model in Table 3 is not directly comparable with the one in Table 1.

Table 3.  Presupernova Properties with Varying Network Size

MZAMS	Network	Nzones	Nsteps	MHe	MCO	12Cign.	MFe	ξ2.5	M4	μ4	Ye,c
(M⊙)				(M⊙)	(M⊙)		(M⊙)		(M⊙)
15.00	Approx	4233	32013	4.319	2.949	0.215	1.527	0.187	1.760	0.087	0.438
15.00	Coproc	4227	33294	4.319	2.946	0.215	1.537	0.179	1.742	0.084	0.437
15.00	Full	4185	39361	4.383	3.014	0.202	1.407	0.139	1.469	0.068	0.433
15.01	Approx	4235	33091	4.308	2.937	0.211	1.373	0.143	1.449	0.088	0.434
15.01	Coproc	4248	35736	4.298	2.925	0.216	1.387	0.135	1.421	0.097	0.432
15.01	Full	4194	36083	4.390	2.998	0.214	1.511	0.195	1.761	0.089	0.437
25.00	Approx	4717	30080	8.663	6.711	0.179	1.604	0.311	1.913	0.107	0.444
25.00	Coproc	4721	30143	8.670	6.722	0.180	1.580	0.306	1.891	0.106	0.443
25.00	Full	4589	42642	8.179	6.322	0.179	1.806	0.240	1.807	0.084	0.438
25.01	Approx	4696	30317	8.633	6.720	0.176	1.617	0.312	1.912	0.108	0.444
25.01	Coproc	4681	30120	8.640	6.670	0.176	1.564	0.289	1.854	0.110	0.443
25.01	Full	4531	43698	8.411	6.513	0.180	1.835	0.252	1.837	0.088	0.438

Note. All quantities are measured at the presupernova stage, except ¹²C_ign (Section 2.2). Here N_zones is slightly larger at the beginning of calculation, and Y_e,c is evaluated at the center.

Download table as:  ASCIITypeset image

Table 3 gives, in addition to the presupernova zoning and time steps, the masses of the iron, carbon–oxygen, and helium cores; the central mass fraction of carbon evaluated just prior to carbon ignition; the central value of the electron mole number, Y_e, in the presupernova star; and the descriptors of presupernova core structure, ξ_2.5, μ₄, and M₄. Though not given, the presupernova masses, radii, luminosities, and effective temperatures are in the same range as the same-mass stars in Table 1.

Based solely on the changes in the core structure and iron core mass for the 15.00 and 15.01 M_⊙ models, one might conclude that using a large network made a substantial difference, and all future runs should take much greater care with the nuclear physics. Though the models do indeed differ, this would not be a valid inference. As shown in Section 2.1 and discussed elsewhere in the paper, changing anything for models in some mass ranges can give qualitatively different answers for core structure. The differences between runs with large and small networks is within the range of differences resulting from increased or decreased zoning (Table 1) or a small change in the star's mass. More telling is the fact that the central carbon abundances, the helium and carbon–oxygen core masses, and, most critically, the central value of Y_e all agree very well. The approximation network is doing a fine job representing the nuclear physics. It is just that other, uncontrollable factors introduce large variations in the core structure of a 15 M_⊙, and to a lesser extent, 25 M_⊙ star. We conclude that the 19 isotope approximation network and quasi-equilibrium hypothesis are quite adequate for surveys like this. Instead, the focus should be placed on studying the impact of stellar physics on the statistical properties over the entire mass range.

2.3. Boundary Pressure

Especially during its post-main-sequence evolution, the radius of a model star is sensitive to surface boundary conditions. Different codes treat the surface in different ways. Many use a boundary condition on pressure or density that sometimes includes a reduction of gravity by the Eddington factor and the inertia term from wind acceleration (e.g., Appenzeller 1970; Heger 1998; Paxton et al. 2015). Others solve the wind equation for conditions at the sonic radius (Grassitelli, private communication, 2018). Still others fit the surface to a stellar atmosphere of varying complexity (e.g., plane-parallel gray or tables; e.g., Chieffi & Straniero 1989; Paxton et al. 2011). KEPLER uses a constant boundary pressure, P_bound, that does not vary during the evolution. The advantage of such an approach is its simplicity, but care must be taken that P_bound is not so large as to greatly alter the solution.

Ideally, P_bound only influences the structure of a tiny mass near the surface. No reasonable value of P_bound affects main-sequence evolution, for example, because the pressure gradient there is very steep. The gradient becomes shallower during the RSG phase, though, and even slightly different radii there can significantly alter the mass loss rate. The timing and development of the surface convection zone is also affected. Both can have a significant effect on the star's location in the HR diagram and, to a lesser extent, the helium core mass and presupernova structure. Traditionally, studies with KEPLER have used P_bound = 50–100 dyne cm⁻² (e.g., Woosley et al. 2002; Woosley & Heger 2007; Sukhbold & Woosley 2014). An exception is the work of Müller et al. (2016), where a much larger value, ∼3500 dyne cm⁻², was used.

In the present study, a value of 50 dyne cm⁻² is employed. Is this low enough? To explore the sensitivity of the results to boundary pressure, a standard 15.00 M_⊙, Model S15.00D in Table 1, was calculated using a variety of boundary pressures from 10 to 6400 dyne cm⁻² maintained throughout the evolution. The resulting HR diagrams and radial histories are shown in Figure 3. For P_bound ≲ 500 dyne cm⁻², the trajectory in the HR diagram is essentially identical for six different choices of ${P}_{\mathrm{bound}}$. The presupernova radius and luminosity do not vary. For larger values of P_bound, though, especially P_bound ≳ 1000 dyne cm⁻², significant variation is found.

Zoom In Zoom Out Reset image size

Even the models with small surface pressures show some variation in final properties. For example, for P_bound between 10 and 400 dyne cm⁻², the final star mass varied between 12.59 and 12.61 M_⊙, while the helium core mass varied between 3.61 and 3.67 M_⊙. Some of this variation is a consequence of the irreducible noise in running any model twice with even small variations in the physics or mass, but part could be a residual sensitivity to ${P}_{\mathrm{bound}}$.

Closely related to the treatment of surface boundary conditions is surface zoning. In our recent studies (Sukhbold & Woosley 2014; Müller et al. 2016), rezoning was allowed to continue all the way to the surface of the star throughout its evolution. This often resulted in coarse zoning near the photosphere during the RSG stage. Fine zoning on the main sequence was lost to dezoning when the temperature and density gradients were shallow. See Figure 1 for a sample comparison with the older model. In the new models, rezoning was not allowed in the outer 250 zones, typically 0.3 M_⊙. Thus, the fine surface zoning shown in Figure 1 was maintained throughout the evolution for all models. Typical zoning near the surface is less than 5 × 10⁻⁴ M_⊙.

The use of finer surface zoning and a reduced boundary pressure generally gave larger radii and increased mass loss during the RSG stage, especially as compared with Müller et al. (2016). Because the helium core mass was not appreciably affected, though, the final core compactness, nucleosynthesis, and remnant masses are not altered by this change.

Three grids of models were calculated for initial main-sequence masses in the range 12–60 M_⊙ using the zoning described in Section 2 and three choices of mass loss rate. Since all stars retained at least a small envelope mass as presupernova stars, the relevant mass loss rate was that of Nieuwenhuijzen & de Jager (1990). The first set of models uses that rate unmodified (${\dot{M}}_{{\rm{N}}}$). This results in so much of the envelope being lost by stars above about 27 M_⊙ that the residual envelope, typically a few M_⊙, becomes unstable and develops density inversions of more than an order of magnitude. We suspect that massive stars with extremely large radii and low-mass envelopes become unstable in nature at this point and lose their remaining hydrogen rapidly (e.g., Sanyal et al. 2015, 2017). Mass loss is not the focus here, though. In order to have a greater range of helium core masses, a second set of models used half of the standard mass loss (${\dot{M}}_{{\rm{N}}}/2$) appropriate for solar metallicity and encountered no instability up to 40 M_⊙. Finally, a third set was calculated with one-tenth the standard mass loss (${\dot{M}}_{{\rm{N}}}/10$) to provide a sparse mass grid between 12 and 60 M_⊙. Considering that mass loss prescriptions are probably uncertain (Renzo et al. 2017; Beasor & Davies 2018), the factors of 2 and 10 could reflect a sensitivity study, or they might be appropriate for stars with reduced metallicity.

Altogether, 1,499 models were calculated in the mass range 12–27 M_⊙ using the standard mass loss rate and 2799 models from 12 to 40 M_⊙ using half that value. The grid spacing was 0.01 M_⊙ for both sets. The third set, with one-tenth the standard mass loss rate, consisted of only 49 models between 12 and 60 M_⊙ with a mesh size of 1 M_⊙. Lower-metallicity stars were not considered, but, except for zero, metallicity would have affected the answer chiefly by changing the mass loss and hence hydrogen envelope mass of the presupernova star. As we shall see, the helium core mass is not entirely independent of the remaining envelope mass, but it is insensitive.

The code was compelled to spend a much greater number of time steps during hydrogen and helium burning by imposing a limit on the maximum time step of 10¹⁰ s. A minimum of roughly 10,000 steps was thus spent burning hydrogen on the main sequence. This helped to weaken the impact of random semiconvective mixing events that affected the final envelope mass and, to a smaller extent, the helium core mass. It also improved the accuracy of the treatment of convection and burning as "split" operations (convection and burning are not implicitly coupled in the code). Very little difference was noted in a few test cases when this maximum step was increased to 5 × 10¹⁰ s.

Minor changes to improve the stability of the code when transitioning to silicon quasi-equilibrium and the convergence criteria were also incorporated. These had insignificant effects on the outcome but improved code performance. Models for the main survey were all calculated on identical processors using the same version of the code and compiler so as to reduce any possible noise introduced by machine architecture (Ohio Supercomputer Center 1987).

One major goal of this work was to explore the effect of resolution on calculations of presupernova structure. The number of zones used in the main survey was typically 4000 × (M_ZAMS/12 M_⊙). Stars more massive than 30 M_⊙ thus had over 10,000 initial mass shells. In all cases, the maximum mass zone anywhere in the star was, at all times, ≲0.01 M_⊙, but zoning was by no means uniform. The initial zoning (Figure 4), which continued to characterize most of the star until the presupernova stage, had much finer zoning in the helium core, down to about 0.001 M_⊙ at the star's center.

Zoom In Zoom Out Reset image size

Except as noted above for hydrogen burning, the time-step criteria were not varied from the previous studies. Nevertheless, the finer resolution of the new runs did result in taking more time steps as smaller zones experienced mass loss or were cycled through the large density contrast at the base of the hydrogen convective shell. Abundances within a given zone convectively coupled to many other zones also required more time steps to change. The typical number of steps for a new calculation is 45,000, about twice the earlier studies' ∼25,000.

The presupernova structure and composition tables for all models from the main sets with ${\dot{M}}_{{\rm{N}}}$, ${\dot{M}}_{{\rm{N}}}/2$, and ${\dot{M}}_{{\rm{N}}}/10$ are available through the Harvard Dataverse archive.⁸

4.1. Observable Properties

The observable properties of the new models are summarized in Figures 5 and 6 and Table 4. Final (presupernova) masses are sensitive to both the uncertain mass loss rate prescription (Nieuwenhuijzen & de Jager 1990) and the uncertain history of the stellar radius. The results vary considerably. For lower-mass stars, the spread in final mass is smaller, since the star only loses a small amount of mass. For bigger stars, though, a large fraction of the hydrogen envelope may be lost. Most of the mass loss occurs during helium burning, and most of that occurs during the late stages when the central helium mass fraction is less than 0.5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Properties of Three Main Model Sets

MZAMS	MpreSN	MHe	MCO	logLpreSN	Teff,preSN	Teff,He dep	RpreSN
(M⊙)	(M⊙)	(M⊙)	(M⊙)	(erg s−1)	(K)	(K)	(1013 cm)
${\dot{M}}_{{\rm{N}}}$ (binned by 1 M⊙)
12	11.19	3.35	2.20	38.39	3440	3670	4.97
13	11.72	3.79	2.53	38.47	3420	3610	5.49
14	12.22	4.20	2.85	38.54	3410	3570	5.97
15	12.73	4.61	3.19	38.60	3405	3545	6.41
16	13.23	5.01	3.53	38.65	3405	3535	6.80
17	13.74	5.42	3.87	38.69	3410	3535	7.15
18	14.12	5.84	4.24	38.73	3420	3530	7.45
19	14.60	6.24	4.60	38.77	3425	3530	7.72
20	14.87	6.66	4.98	38.81	3430	3535	8.07
21	14.97	7.09	5.35	38.85	3430	3535	8.51
22	15.09	7.49	5.71	38.90	3430	3540	8.95
23	15.11	7.90	6.05	38.94	3435	3550	9.36
24	14.72	8.34	6.44	38.98	3450	3565	9.73
25	14.49	8.75	6.80	39.01	3470	3590	10.01
26	14.27	9.15	7.15	39.04	3500	3620	10.20
${\dot{M}}_{{\rm{N}}}/2$ (binned by 1 M⊙)
12	11.82	3.39	2.22	38.40	3460	3685	4.95
13	12.59	3.82	2.55	38.47	3450	3630	5.44
14	13.35	4.24	2.88	38.54	3440	3595	5.90
15	14.12	4.65	3.22	38.60	3440	3580	6.31
16	14.87	5.07	3.58	38.65	3440	3570	6.71
17	15.62	5.49	3.93	38.70	3450	3565	7.04
18	16.37	5.91	4.30	38.74	3460	3565	7.32
19	17.01	6.36	4.69	38.78	3465	3560	7.65
20	17.66	6.79	5.08	38.82	3465	3560	8.04
21	18.27	7.22	5.45	38.87	3460	3560	8.49
22	18.69	7.67	5.85	38.92	3455	3560	9.00
23	19.25	8.07	6.20	38.96	3455	3565	9.42
24	19.58	8.52	6.59	39.00	3455	3565	9.86
25	19.95	8.94	6.94	39.03	3455	3570	10.23
26	20.14	9.38	7.31	39.06	3460	3575	10.60
27	20.36	9.80	7.67	39.09	3465	3585	10.92
28	20.66	10.19	7.99	39.11	3470	3590	11.19
29	20.61	10.65	8.38	39.14	3480	3600	11.50
30	20.72	11.08	8.73	39.17	3495	3610	11.75
31	20.83	11.49	9.10	39.19	3510	3625	11.97
32	20.85	11.95	9.48	39.21	3525	3640	12.18
33	20.98	12.37	9.86	39.23	3540	3655	12.37
34	21.03	12.85	10.27	39.26	3560	3675	12.55
35	21.08	13.33	10.70	39.28	3580	3690	12.74
36	21.36	13.69	11.01	39.29	3590	3705	12.86
37	21.40	14.18	11.46	39.31	3610	3725	12.99
38	21.63	14.56	11.80	39.33	3625	3735	13.15
39	21.80	14.98	12.18	39.34	3635	3750	13.29
${\dot{M}}_{{\rm{N}}}/10$ (not binned)
12	11.88	3.21	2.09	38.36	3480	3725	4.69
15	14.76	4.46	3.07	38.58	3465	3615	6.05
20	19.48	6.66	4.96	38.81	3500	3590	7.76
25	24.00	8.86	6.89	39.02	3495	3605	9.88
30	27.75	11.39	8.95	39.17	3510	3615	11.73
35	31.55	13.89	11.05	39.29	3535	3625	13.17
40	35.25	16.23	13.19	39.37	3575	3660	14.23
45	36.65	21.41	18.22	39.53	3610	3650	16.73
50	41.40	21.21	17.87	39.52	3635	3710	16.28
55	44.27	23.73	20.20	39.58	3655	3715	17.26
60	47.31	25.64	22.00	39.62	3680	3740	17.79

Note. All quantities were measured at the presupernova stage except T_{eff,He dep}, which was measured when the helium mass fraction in the core was depleted to 1%. All T_eff values were rounded to multiples of 5.

Download table as:  ASCIITypeset image

Since the luminosity does not vary greatly during helium burning, the amount of mass lost and its uncertainty depends mostly upon the radial history of the star: R^0.8 in the Nieuwenhuijzen & de Jager (1990) formalism, i.e., the amount of time the helium-burning star spends as a blue supergiant (BSG) or RSG. It is well known that semiconvective and overshoot mixing play a key role in determining the ratio of time spent as each (Lauterborn et al. 1971). In the present calculations, small changes in individual semiconvective mixing episodes affect the timing of the development of a deep surface convection zone, and hence the movement of the star to the red. With more semiconvection, the star spends a greater fraction of its helium-burning lifetime as a BSG and hence loses less mass, ending its life with a larger hydrogen envelope still intact. However, RSGs have the converse behavior. Less semiconvection, e.g., strictly Ledoux convection, favors a longer time as an RSG and thus reduces the threshold for a massive star to lose its envelope. Though the spread in Figure 5 looks large, the uncertainty is usually a small fraction of the total mass lost. For example, a 25 M_⊙ star with the standard mass loss rate might end up as a 12 or 15 M_⊙ star. That is, it might lose 10–13 M_⊙, a range of about 25%. Lower-mass models have a smaller variation.

The final masses for stars calculated using ${\dot{M}}_{{\rm{N}}}$/2, of course, are larger. Bigger stars on the main sequence then retain more envelope by the time they die, and we were thus able to determine the core structure of more massive stars. For the most massive models considered, the final hydrogen envelope mass was roughly 4 M_⊙ for the ${\dot{M}}_{{\rm{N}}}$ models and 5 M_⊙ for the ${\dot{M}}_{{\rm{N}}}/2$ models. For smaller envelope masses, the radius expanded beyond 1.5 × 10¹⁴ cm and became unstable in the KEPLER code due to recombination. In these cases, the envelope could only be retained on the star by the addition of a large, unrealistic boundary pressure that resulted in gross density inversions in the outer envelope. We suspect that this is a real instability in nature, that once the envelope mass decreases below a critical value in a very massive star and the radius extends beyond 10 au, the remaining envelope is lost through nonsteady processes (Sanyal et al. 2015, 2017). The luminosity in the hydrogen envelope is a substantial fraction of the Eddington value. This interesting possibility is deferred for a later study.

Since a larger pressure can give a smaller radius and reduced mass loss, different choices of boundary pressure can lead to significant differences in the final mass. The boundary pressure used in Müller et al. (2016), for example, several thousand dyne cm⁻², was larger than what we now believe reasonable and much larger than that used in the present study, 50 dyne cm⁻². This probably accounts, at least partly, for the larger final masses found in the Müller et al. (2016) study and their ability to study higher-mass stars using a single mass loss prescription.

This boundary pressure does not appreciably affect the compactness of the cores for stars of given main-sequence masses, however. This is because the large spread in final star masses (Figure 5) is not reflected in the helium core mass (Section 4.2) or the luminosity of the star (Figure 6). Since the core structure, which is the main focus of this paper, is chiefly sensitive to the helium core mass and not the envelope mass, the spread in final masses seen in Figure 5 does not affect our major conclusions.

Figure 6 and Table 4 give the final luminosities and effective temperatures of our presupernova models. The luminosity of the presupernova star depends chiefly on the helium core mass and is approximately given by L_preSN ≈ 5.77 ×10³⁸ (M_He/6 M_⊙)^1.5 erg s⁻¹. The effective temperature, essentially bounded by the Hayashi limit, remains pegged at close to ∼3500 K for all presupernova stars. Thus, the radius of the star increases as ${M}_{\mathrm{He}}^{0.75}$. Except for a few presupernova stars, most RSGs will be observed during their helium-burning stage, where their effective temperatures will be hotter. Our calculations show a systematic ∼200 K offset between the effective temperatures at helium depletion and presupernova stars. Earlier in helium burning, the temperature will be even hotter; thus, T_eff ≈ 3550 K should be a lower bound to what is seen. Larger values are also expected for subsolar metallicity and lower values for presupernova stars. These numbers are in good accord with measurements of the hottest RSGs (Levesque et al. 2005; Davies et al. 2013).

4.2. Core Masses

Figure 7 shows the helium, carbon–oxygen (CO), and iron core masses for all of our presupernova models with ${\dot{M}}_{{\rm{N}}}$ and ${\dot{M}}_{{\rm{N}}}/2$. The dispersion in helium and CO core masses for a given main-sequence mass is small, much less than the hydrogen envelope masses inferred in Figure 5. This implies that the helium and CO core masses are only slightly affected by the assumed mass loss rate and probably not very sensitive to metallicity. The helium core mass as a function of main-sequence mass is approximately ${M}_{\mathrm{He}}\approx 6.46\ {({M}_{\mathrm{ZAMS}}/20)}^{1.27}$, which can be combined with the previous relation in Section 4.1 to give L_preSN ≈ 6.5 × 10³⁸ (M_ZAMS/20 M_⊙)^1.92 erg s⁻¹. Since presupernova core compactness is chiefly a function of helium core mass (Sukhbold & Woosley 2014), these results suggest a near-universal dependence of presupernova core structure on initial mass, provided mass loss does not remove the entire hydrogen envelope.

Zoom In Zoom Out Reset image size

Iron core masses increase, on average, with increasing stellar mass, reaching a maximum of about 2.0 M_⊙ for the most massive stars studied (>40 M_⊙). Still larger iron cores, up to about 2.5 M_⊙, characterize more massive stars in the pulsational pair-instability range (70–140 M_⊙; Woosley 2017). For stars below 23 M_⊙, the iron core mass is markedly multivalued for stars with nearly the same initial mass. This reflects the operation of multiple shells of carbon and oxygen burning, as will be discussed further in Section 5. The two major branches of iron core masses below 20 M_⊙, which are most of the stars that leave neutron star remnants, might result in bimodality in the neutron star mass function.

4.3. Core Structure

4.3.1. Measures of "Explodability"

Early theoretical studies of supernovae noted a strong correlation of the rapidly declining density outside the iron core with the degree of difficulty encountered in trying to explode the star using the neutrino energy transport (e.g., Burrows & Lattimer 1987; Fryer 1999). O'Connor & Ott (2011) introduced a simple, single-parameter measure of this density decline called the "compactness parameter,"

$$\begin{eqnarray}&&{\left.{\xi }_{{\rm{M}}}=\displaystyle \frac{M/{{\rm{M}}}_{\odot }}{R({M}_{\mathrm{bary}}=M)/1000\mathrm{km}}\right|}_{{t}_{\mathrm{bounce}}},\end{eqnarray} \tag{ 3 }$$

where R(M_bary = M) is the radius of the Lagrangian mass shell enclosing mass M in the presupernova star. The fiducial mass is often chosen as the innermost 2.5 M_⊙, so that for a wide range of initial masses, it not only encloses the iron core but samples enough of the overlying shell material around it. Though the parameter is defined to be evaluated at the time of bounce, it is more often measured at the time of presupernova (when the collapse begins), since the systematics are insensitive to slightly different fiducial choices (Sukhbold & Woosley 2014). Subsequent studies by Ugliano et al. (2012), O'Connor & Ott (2013), and Sukhbold et al. (2016) showed strong correlation between the ease with which model stars exploded and the ξ parameter, in the sense that stars with small ξ, i.e., steep density gradients outside the iron core, were easier to explode using a standard, albeit approximate, set of supernova physics. Although this is a useful starting point, a single parameter conveys limited information about the structure of the core, and more physics-based representations followed.

In particular, Ertl et al. (2016) suggested an alternative two-parameter characterization based upon M₄, the mass coordinate in solar masses, where the entropy per baryon reaches 4.0 k_B, and the radial gradient, μ₄, of the density at that point. In practice, μ₄ was obtained by evaluating the change in mass over the change in radius between several mass shells separated by 0.3 M_⊙ in the vicinity of M₄, i.e.,

$$\begin{eqnarray}&&{\mu }_{4}=\displaystyle \frac{{dm}/{{\rm{M}}}_{\odot }}{{dr}/1000\,\mathrm{km}},\end{eqnarray} \tag{ 4 }$$

where dm = 0.3 M_⊙. Smaller values of μ₄ thus imply steeper density gradients (less change in enclosed mass when the radius changes). The quantity M₄ has long been used to locate the steep density decline often associated with a strong oxygen-burning shell in the presupernova star (e.g., Woosley & Heger 2007). Where the entropy per baryon abruptly rises at nearly constant temperature, the density declines. When the core collapses, the arrival at the neutrinosphere of lower-density matter reduces the "ram pressure" and facilitates the launch of an outgoing shock. Besides its location, it is helpful to characterize the rate of the density change, which is the role played by μ₄.

It is reasonable that higher neutrino luminosities and smaller accretion rates on the proto-neutron star should favor explosion. Ertl et al. (2016) made the case that μ₄, multiplied by its radius (i.e., the radius at which M₄ is found) and divided by a collapse timescale, is a surrogate for the accretion rate at the time of explosion. The product μ₄M₄ is similarly a surrogate for the accretion luminosity. This assumes, as the models suggest, that the radius where the neutrino luminosity is generated does not vary greatly with mass. Explosion is thus favored by small μ₄ (accretion rate) and large μ₄M₄ (luminosity). Ertl et al. (2016) determined a simple condition expressed as a straight line,

$$\begin{eqnarray}&&{\mu }_{4}={k}_{1}({\mu }_{4}{M}_{4})\,+\,{k}_{2},\end{eqnarray} \tag{ 5 }$$

such that models that lay below the line (lower μ₄) were more likely to explode than those above. The values k₁ and k₂ were determined from several large-scale simulations and leading models for SN 1987A. For one representative set, k₁ = 0.194 (M_⊙)⁻¹ and k₂ = 0.058 ("N20" model in their Table 2).

More recently, Müller et al. (2016) adopted a different approach using a semi-analytic model for the explosion based upon presupernova properties. By using a fuller description of the presupernova star than afforded by just one or two parameters, they were able to approximately estimate not only the success or failure of the explosion but also the explosion energy, ejected nickel mass, and compact remnant mass. This approach does not include, in its present form, several important pieces of explosion physics, such as proto-neutron star cooling, self-consistent fallback, and explosive nuclear burning, and thus it is not as powerful as the parameterized simulations such as Sukhbold et al. (2016). Its main power lies is its ability to rapidly calculate general trends in explosion properties that depend on presupernova structure. Their approach uses five physically motivated parameters to determine the outcome. Here we analyze our new presupernova models using their standard values: β_expl = 4, ζ = 0.7, α_out = 0.5, α_turb = 1.18, and τ_1.5 = 1.2 (see their Table 1). For the definitions of these parameters and further technical details, see Müller et al. (2016).

4.3.2. Results for the Compactness Parameter

Figure 8 shows the compactness parameter, ξ_2.5, i.e., ξ for our three model sets, measured at a fiducial mass of 2.5 M_⊙ at the time the presupernova collapse speed reaches 900 km s⁻¹. One notable feature seen in Sukhbold & Woosley (2014) is the relatively large "scatter" of the compactness parameter between 18 and 22 M_⊙, which is now shifted to roughly 14 and 19 M_⊙ (panel (a)) with updated neutrino losses (Section 2). As has long been noted (e.g., Timmes et al. 1996), ∼20 M_⊙ marks the transition region between exoergic convective central carbon burning at lower masses and endoergic radiative burning at higher ones. That is, at high mass, the energy produced by carbon fusion does not exceed the neutrino losses in the star's center, and it proceeds from helium depletion to oxygen ignition on a short Kelvin–Helmholtz timescale. As a result, subsequent burning—especially in shells, where the net energy generation is positive—proceeds differently. The carbon convective shells move inward (with increasing initial mass) until this transition point and are weaker and greater in number (Barkat 1994).

Zoom In Zoom Out Reset image size

Quite noticeable in the new results and the prior work of Müller et al. (2016) are concentrations of points in the compactness parameter plot, which were not clearly seen in our earlier studies due to much coarser increments in mass space. Points are not randomly scattered between some local maximum and minimum value. Note, for example, the existence of two, and possibly three, discrete solutions for helium core masses near 4.8 M_⊙ (panel (b)). Multiple branches also exist at other masses but with less clarity. For a helium core mass of 5.6 M_⊙, the compactness parameter for stars with nearly identical masses varies by a factor of three, from ξ_2.5 = 0.1 to 0.3. Such a large variation spans the range of stars that might explode or collapse into black holes and suggests that both outcomes are possible for relatively low-mass stars of almost the same initial mass.

When all three new sets are plotted as a function of initial mass (panel (c)), the models with lower mass loss rates are slightly shifted with respect to models with higher mass loss rates. The shifts are hardly noticeable at a lower initial mass (<20 M_⊙), but they grow with increasing mass, and for our most massive models, the ${\dot{M}}_{{\rm{N}}}/10$ set is shifted by a few M_⊙ with respect to the ${\dot{M}}_{{\rm{N}}}/2$ set. The primary reason for this shift is that for a given initial mass, models with lower mass loss rates attain slightly larger He cores and slightly different compositions. There is also a general underlying scatter, similar to that of the total star mass shown in Figure 5, some of which results from the slight change of He core mass due to the semiconvective effects discussed in Section 4.1.

Panel (d) of Figure 8 shows the compactness as a function of He core mass for all of our models. Since the initial mass scales cleanly with He core mass (Figure 7), the structure of the curve stays nearly identical; however, the abovementioned shifts and scatters are mostly gone. The He core mass and, within it, the CO core mass are the chief determinants of core structure, not the final supernova mass, which includes the hydrogen envelope. Uncertainties in mass loss rates are thus not particularly relevant to the compactness, except for determining the main-sequence mass above which the entire envelope is lost.

As was noted by O'Connor & Ott (2011) and studied many times since, the compactness plotted against mass is nonmonotonic and highly structured, with two distinct peaks in the vicinity of 21 and 35 M_⊙ (Sukhbold & Woosley 2014). A shallow dip occurs near 50 M_⊙, followed by a slow rise at still higher masses (Sukhbold & Woosley 2014), until finally, the pulsational pair domain is reached near 70 M_⊙. As discussed in Sukhbold & Woosley (2014) and Sukhbold (2016), this structure is crafted out of an overall gradually rising curve by the operation of an advanced-stage shell-burning episode of one fuel modulating the core-burning episode of the next fuel. In particular, the structure near the first peak is primarily driven by the effect of shell carbon burning on core oxygen burning, while the structure near the second peak is driven by the effect of shell oxygen burning on core silicon burning. At very high mass, oxygen burning ultimately ignites a strong shell outside of 2.5 M_⊙, and therefore the compactness stays high.

Figure 9 compares the new results with those of Sukhbold & Woosley (2014) and Müller et al. (2016) and illustrates these points. Most notable is a shift of the new models above ∼14 M_⊙ downward in mass by about 2–3 M_⊙ compared with that of Sukhbold & Woosley (2014). This is a consequence of fixing the bug in the neutrino losses as discussed in Section 2 and not due to differing resolution. Models below 14 M_⊙ are affected very little, but the more massive stars change appreciably. The region of variability that was between 18 and 22 M_⊙ in Sukhbold & Woosley (2014) is now shifted down to 14–19 M_⊙. If anything, that variability is now greater with the finer resolution, perhaps due to the larger number of models surveyed.

Zoom In Zoom Out Reset image size

Indeed, the comparison with Müller et al. (2016), who used the corrected neutrino loss rate and sampled many more models, is excellent, both in the location of structures and the range of ξ_2.5 within those structures. This agreement persists despite the use of 4–10 times more mass shells in each of the new models and a radical decrease in the surface boundary pressure. The latter affected the mass lost by the star but did not appreciably affect the helium core masses. There is no reason to believe that still finer zoning, smaller time steps, or a different reaction network will greatly alter these results, unless the code physics itself (reaction rates, semiconvection, rotation, etc.) is changed. The studies of Müller et al. (2016) and the present work are mutually confirming.

4.3.3. New Results and the Ertl Parameterization

An important subsidiary question is whether ξ_2.5 is really the best measure for presupernova core structure. Might some of the variability seen in Figure 8 be simply because of the choice of a single arbitrary point in the star to sample its structure? Perhaps other parameterizations might give less variability and more reliable predictions? In particular, ξ_2.5 is sensitive to recent shell activity in the vicinity of 2.5 M_⊙ that might not always describe well what went on deeper inside.

The Ertl parameterization of our results is given in Figure 10. As previously discussed (Section 4.3), the points beneath the dashed line represent models that are more likely to explode in a simple neutrino-transport scheme. Multivalued solutions are clearly seen, especially for stars with μ₄M₄ less than 0.25. In some ranges, the difference between solutions is enough to significantly affect the probable outcome of the explosion. Stars with μ₄M₄ less than 0.25 typically have masses less than 20 M_⊙. For more massive stars—and, in particular, for μ₄M₄ greater than 0.25—explosion seems quite unlikely. The reasons for the multivalued solution are discussed in Section 5.

Zoom In Zoom Out Reset image size

Also striking is the much tighter clustering of points in the Ertl representation of our new models (Figure 10), in contrast with an equivalent plot of "compactness" (Figure 8). The results are less noisy and show less overall variability. This is due to the obvious correlation between two Ertl parameters on one hand, but it also presumably reflects the better representation of core structure by a two-parameter, physics-based representation than a single parameter anchored to a single point in the core. It also suggests that much of the "noise" in Figure 8, especially for initial masses of 14–19 M_⊙, may not be so much a consequence of "nonconvergence" of the models but rather a poor representation of the results.

Why is the core structure of presupernova stars nonmonotonic and multivalued for some ranges of mass? The short answer is that the advanced-stage evolution of massive stars involves two to three carbon-burning shells (plus central core burning for models below ∼19 M_⊙) and one or two oxygen shells (plus core burning). Combinations of these shells lead to variable outcomes but not a continuum of all possibilities, because the shells have finite sizes and their number is an integer. The transitions are abrupt, and small changes upstream in the strength or extent of one or more shells can send a star down one path or another. This is especially true for stars below 19 M_⊙, where there are more carbon-burning shells.

Globally, ξ_2.5, μ₄, and M₄ all increase gradually with mass. Higher-mass stars have greater entropy in their middles and are less degenerate in their final stages. Greater degeneracy leads to an increased central concentration of the mass and "core convergence." For lighter stars, the presupernova structure more resembles a white dwarf embedded in a low-density envelope, where the density declines rapidly at the edge of the central "white dwarf." This effect is clearly at work in the lightest stars surveyed. Below 13 M_⊙, compactness defined at a mass of 2.5 M_⊙ has little meaning as a measure of explodability, since the fiducial point is not even inside of the helium-burning shell. For stars in this mass range, the helium-burning shell always lies outside 2 × 10⁹ cm, thus guaranteeing a small compactness parameter. For a 10 M_⊙ presupernova star, 2.5 M_⊙ is even outside the helium core and in the hydrogen envelope. For all parameterizations of core structure considered in this paper, these light stars will always explode and need no further discussion here.

Degeneracy remains an important consideration, though, even for the larger stars. A 12 M_⊙ presupernova star has a core that is degenerate (degeneracy parameter η ≡ μ/kT > 0) out to 1.46 M_⊙. For 30 M_⊙, degeneracy extends out to 1.87 M_⊙. Even for 100 M_⊙, the iron core is degenerate out to 2.1 M_⊙. Most of the increase in degeneracy occurs between carbon depletion and oxygen depletion in the core when the core is efficiently cooled by the neutrinos (Sukhbold & Woosley 2014). Degeneracy affects the core structure, making the (thermally adjusted) Chandrasekhar mass relevant for all stars that are likely to explode as common supernovae. On the other hand, the shells outside what will be the iron core in the presupernova star—and, in particular, the carbon- and oxygen-burning shells—are always nondegenerate. Strong burning thus leads to expansion that affects both ξ_2.5 and μ₄.

The upper panel of Figure 11 shows a generally strong correlation between the compactness parameter, ξ_2.5, and M₄, except at very high and low masses. A correlation is expected, since a sharp density decline close to the iron core requires a large radius to enclose 2.5 M_⊙, thus small ξ_2.5. The sharp density decline implies a deeper location where the entropy per baryon exceeds 4.0 k_B and thus a smaller M₄ and a steeper mass gradient at that point as well. The reversal of the plot for the most massive models is a consequence of both the increasing central entropy and the outward migration of the first oxygen shell. This migration causes a nonmonotonic dependence on mass for both M₄ and ξ_2.5, as shown in Figure 8, but M₄ starts to decline at a slightly lower initial mass due to increasing entropy. Therefore, when the location of M₄ recedes from ∼2.5 to 2 M_⊙, the ξ stays roughly constant and results in the reversal near ∼35 M_⊙.

Zoom In Zoom Out Reset image size

The lower panel of Figure 11 shows that the behavior of compactness and M₄ for stars lighter than about 20 M_⊙ is clearly multivalued. Three distinct branches are apparent for the lighter stars. There is no reason to believe that a more dense grid of calculations would randomly fill out the spaces between the branches.

Figure 12 delves deeper into this correlation between compactness and M₄ and helps to understand why both are multivalued relations of mass for moderate-mass stars (13 M_⊙ < M_ZAMS < 19 M_⊙). The first panel shows M₄ in the presupernova star as a function of helium core mass. A cleaner pattern results from using the helium core mass instead of the main-sequence mass, since it eliminates the variations due to envelope mass loss (Figure 5). As the figure shows, M₄ is usually pegged to the "oxygen shell," by which we mean the location in the presupernova star where the energy generation from oxygen fusion (excluding neutrino losses) is at maximum. This is frequently, though not always, interior to the boundary of the "silicon core," inside of which the silicon mass fraction is greater than the oxygen fraction. Though one might naively assume that the oxygen-burning shell is at the edge of the silicon core, this is not generally true. There are often gradients on the silicon-to-oxygen ratio left behind by receding convection or radiative burning. In these cases, the oxygen shell lies inside the silicon core. On the other extreme, in some stars, the oxygen shell can be so active as to merge with the carbon- and neon-burning shells, producing one large convective shell where all three are burning vigorously and the silicon core and oxygen shell are coincident. When this occurs, the compactness is usually small (Figure 16 of Sukhbold & Woosley 2014). Many of the pink, brown, and green points in the top panel of Figure 12 are of this sort.

Zoom In Zoom Out Reset image size

It is not surprising that the oxygen-burning shell, the strongest burning shell during the late stages of stellar evolution, is usually the location of a jump in entropy. The oxygen fusion rate is very temperature-sensitive, so oxygen burns at a nearly constant temperature. The overlying star expands, decreasing the local density and thus increasing the entropy. The quantity μ₄ in the Ertl parameterization is a measure of the strength of this burning.

The middle panel of Figure 12 shows that the behavior of ξ_2.5 in this same mass range (equivalent to panel (b) of Figure 8) correlates reasonably well with that of M₄, including, approximately, the location and extent of clusters of solutions. Within a cluster, there is usually a quasi-linear relation with a well-defined slope. A larger value of M₄ implies a larger value of ξ_2.5. The farther out in the star the strong burning shell, the less centrally concentrated its density. This correlation is not perfect, though. Sometimes long lines in the M₄ plot, e.g., the long string of points from M_He = 4.1 to 5.2 M_⊙, break into several segments with different ξ_2.5. This is partly due to the arbitrary pinning of ξ_2.5 to a single mass shell, but also because specifying M₄ alone does not measure the strength of the burning there.

Instead, as might be expected, the bottom panel of Figure 12 shows that ξ_2.5 correlates better with μ₄, the mass gradient at M₄. This correlation is very strong, since both quantities are sensitive to the strength of the most active burning shell in the last days of the star's life. Why, though, are M₄ and μ₄ not randomly scattered between their extrema? Why the patterns? The top panel of Figure 12 suggests a reason. Here M₄ traces the location of the strongest oxygen-burning shell. Oxygen burns there because it is at the deepest location that has not already depleted oxygen by prior shell burning. That location is, in turn, set by the extent of oxygen-burning shell(s) during the previous evolution, which in turn depends upon the entropy structure set up during carbon shell burning.

Sukhbold & Woosley (2014) pointed out this correlation between the core compactness and the extent of the first oxygen-burning convective shell (their Figure 14). Figure 13 shows a related quantity for the new model set. The figure shows the mass of the "silicon core," the mass interior to which oxygen has burned out, at the time silicon burning ignites in the star's center. Here a specific central temperature, 3.0 × 10⁹ K, was chosen in order to make the plot, but the conditions reflect what the oxygen-burning shell (or shells) have accomplished prior to silicon ignition. Visual inspection of the convective history of these 700 models shows that the silicon core mass at silicon ignition is very nearly equal to the maximum extent of the second convective oxygen-burning shell for stars below ∼14.6 M_⊙ (helium core mass 4.2 M_⊙) and the first oxygen convective shell for more massive stars up to at least 20 M_⊙. In this plot, we see the clearest evidence yet for regular but nonmonotonic and occasionally multivalued behavior. Two helium cores of very nearly 4.2 M_⊙ can give rise to presupernova cores with structures in one of two well-defined states. Slight shifts in mass, composition, or even numerical approach (including zoning) can send the star one way or another.

Zoom In Zoom Out Reset image size

Further analysis of the convective histories reveals the systematics behind this behavior—at least in the 1D code, if not in nature.

• 1.  
  The little jump at M_He = 3.7 M_⊙ (M_ZAMS = 13.2 M_⊙) from M_Si = 1.7 to 1.6 M_⊙ reflects the operation of the third convective carbon-burning shell. Below this mass, the third shell ignites inside the former full extent of the second shell and outside the effective Chandrasekhar mass. Thus, the core oxygen burning starts to burn in a smaller extent while this outermost carbon shell operates. For masses above this value, the third shell ignites at the outer boundary of the previous shell and within the Chandrasekhar mass. Core oxygen burning now has to wait until this third shell is complete; therefore, the base of this third shell sets (approximately) the extent of the first oxygen convective shell and thus the base of the second one (top left panel of Figure 14).
• 2.  
  The much larger decrease in M_Si from 1.7 to 1.4 M_⊙ at M_He = 4.2 M_⊙ (M_ZAMS = 14.6 M_⊙) is due to the diminished significance of the second oxygen-burning shell. The silicon core shrinks to the extent of the first oxygen-burning shell, though the second shell continues to be sporadically important for a time (bottom left panel of Figure 14).
• 3.  
  The jump and wild variations starting at M_He = 5.2 M_⊙ and extending up to 5.7 M_⊙ (M_ZAMS = 17.2–19.0 M_⊙) mark the transition from convective carbon burning to radiative burning. After the transition near ∼19 M_⊙, carbon no longer burns exoergically at the center of the star, and there are only two convective shells before oxygen burning rather than three or more. During this transition, carbon core and shell burning vary greatly in location and extent. The large spread in M_Si and, ultimately, in compactness reflects the irregularity of this transition (right panels of Figure 14).

Zoom In Zoom Out Reset image size

Other shell interactions result in the weaker multiple branches seen in higher-mass stars in Figure 8. Despite the specific masses given in the list above, none of these transitions are abrupt, and, given slight nudges, the star may oscillate from one to the other solution when it is close to boundaries. This leads to the "multivalued" behavior. Although complicated and probably sensitive to the one-dimensional treatment of the problem, the conclusion is that the presupernova structure in stars from 13 to 19 M_⊙ results from an interplay of convective carbon- and oxygen-burning shells after carbon ignition.

Putting these various factors together, a deterministic picture emerges. Multiple carbon-burning stages—core burning below 19 M_⊙ and two or more episodes of shell burning—act to sculpt an entropy distribution in the core so that oxygen shell burning, when it occurs, ignites at and extends to various mass shells. Whether the number of carbon shells is two or three or four and the number of oxygen shells is one or two strongly affects the core structure of a presupernova star in the mass range where a significant number of events are observed. These changes in the extent of convective shells, which amplify small differences in the earlier evolution, may end up determining whether the star explodes or makes a black hole.

Since the new models differ from those of Sukhbold & Woosley (2014) and Sukhbold et al. (2016) and, to a lesser extent, Müller et al. (2016), it is worth revisiting some of the conclusions of those papers regarding compact remnants using the new models.

Figure 15 shows the models now expected to explode based on the criteria of Ertl et al. (2016) and Müller et al. (2016). This figure can be compared with Figure 4 of Ertl et al. and Figure 6 of Müller et al., which it very closely resembles. The "N20" engine parameterization was used to make the comparison with Ertl et al. (2016), and the "standard" choices of five parameters were used for the Müller et al. (2016) method.

Zoom In Zoom Out Reset image size

The most significant difference with the Ertl criterion is the shift of the pattern for models above about 14 M_⊙ to slightly (∼10%) lower masses. Changes in average quantities, like the remnant masses, are small, however. Figure 16 shows a comparison of IMF-weighted remnant mass distributions of imploding models assuming a Salpeter (1955) initial mass function with a power of −2.35. The usual two cases are considered: (a) collapse of the helium core but ejection of the hydrogen envelope and (b) collapse of the entire presupernova star. Case (a) is more appropriate for stars where the envelope is lost to a binary companion or, for very massive stars, to winds. The envelope might also be lost to a very weak explosion that did not unbind the helium core (Quataert & Shiode 2012; Lovegrove & Woosley 2013; Fuller 2017). Case (b) is for more robust explosions. In each case, there is an element of uncertainty because the final mass depends on the mass loss rate, e.g., for RSGs in case (b) and Wolf–Rayet stars in case (a). In making Figure 16, but not in computing the averages, only the mass range covered by our ${\dot{M}}_{{\rm{N}}}$ mass loss survey, 12–27 M_⊙, was included.

Zoom In Zoom Out Reset image size

Astronomers, of course, observe black holes coming from all masses of stars, not just 12–27 M_⊙. To compute the averages for solar-metallicity stars (and the answer will be sensitive to metallicity), the new results were supplemented with the prior models of Sukhbold et al. (2016) for stars more massive than 27 M_⊙. For just the limited mass range 12–27 M_⊙, the average black hole masses from Sukhbold et al. (2016) were previously 6.95 and 15.9 M_⊙, respectively, for the helium core and full star assumptions. The corresponding new numbers are 6.50 and 14.3 M_⊙. Considering the entire mass range of stars that experience iron core collapse, and using the N20 parameters in Sukhbold et al. (2016), the previous averages were 9.35 and 13.7 M_⊙. The equivalent new numbers are 8.25 and 13.7 M_⊙. A less than 1% adjustment has been applied to Table 6 of Sukhbold et al. (2016) based upon a slightly different way of interpolating the grid.

Similar corrections can be estimated for the neutron star gravitational masses. Lacking an explosion model, it is assumed, based on prior experience (i.e., Ertl et al. 2016), that the baryonic mass of the resulting neutron star is usually equal to M₄. After appropriately correcting for neutrino mass loss (same as in Müller et al. 2016), the new average neutron star gravitational mass in the range 12–27 M_⊙ is 1.45 M_⊙. Explosions below 12 M_⊙ were calculated by Sukhbold et al. (2016) using presupernova models in which the neutrino bug had been fixed. They can thus be combined with the current set. Assuming neutron star production above 27 M_⊙ is negligibly small; the new global average neutron star mass is 1.38 M_⊙. In Table 6 of Sukhbold et al. (2016), the corresponding global average was 1.41 M_⊙.

The agreement of Figure 15 with the earlier work (Figure 6 of Müller et al. 2016) is excellent, and it is thus expected that the explosion outcomes will also be very similar. Indeed, plots of neutron star mass (Figure 17), black hole mass, and explosion energy as a function of initial mass (not shown) are virtually indistinguishable from the panels in their Figure 2. This suggests that the difference in envelope structure and zoning between the series of Müller et al. (2016) and the present paper did not have much influence on the core structure and statistical explosion properties.

Zoom In Zoom Out Reset image size

Figure 15 shows that a significantly larger fraction of stars in the interesting mass range 15–20 M_⊙ explode using the Müller et al. (2016) formalism, even though "N20" is one of the more energetic formulations of the Ertl et al. (2016) model. Conservatively, this could be regarded as an uncertainty in outcome until more realistic simulations of the actual explosion can be done. The various energies in the Ertl et al. (2016) model came, however, from calibrating a central engine with 1D neutrino transport, a shrinking proto-neutron star, and fallback to SN 1987A using various presupernova models and might, for now, be considered the more realistic of the two. This calibration to SN 1987A has its own uncertainties, though, especially since the structure of the presupernova star could be different in more realistic binary merger models (e.g., Menon & Heger 2017 and references therein).

Figure 17 gives the expected neutron star masses that result when our new presupernova models are analyzed using the formalism and standard parameters of Müller et al. (2016). Fallback is neglected in their analysis, so the baryonic mass of the proto-neutron star is equal to the mass coordinate where the neutrino-driven engine shuts off a few hundred ms after shock revival. As mentioned before, this correlates strongly with M₄ (top panel). There are branches of neutron star masses above the M₄ line, which primarily originate from less massive progenitors (M_ZAMS < 15 M_⊙), reflecting the operation of a prominent second oxygen-burning shell and the resulting shallow entropy profile. The most massive neutron stars have a baryonic mass that is bounded by the base of the convective carbon-burning shell, which is located much further out than M₄ in the presupernova star. The bottom panel of Figure 17 shows the expected neutron star (gravitational) mass as a function of main-sequence mass. These results agree quite well with those of Müller et al. (2016; see their Figure 2(c)), including the existence of very massive neutron stars produced for main-sequence stars of 14–15 M_⊙ and the highly variable nature of the solution between 14 and 19 M_⊙. For the majority of models, the resulting neutron star mass is tightly correlated with M₄ (Figure 12) and the oxygen-depleted core at the time of silicon ignition (Figure 13), and in certain initial mass ranges, the explosion of models with nearly identical mass (or slightly different input physics) can result in very different neutron stars.

It is not expected that the nucleosynthesis and light curves calculated by Sukhbold et al. (2016) will be significantly altered by using the new models, though further exploration is certainly encouraged. In particular, the deficiency of light s-process elements seen in Sukhbold et al. (2016) will persist in the new model set, since most of the production of these elements is due to the most massive stars (Brown & Woosley 2013), which still fail to explode. The fraction of solar-metallicity stars above 9 M_⊙ exploding as supernovae was 74% based on the N20 parameterization (Table 6 of Sukhbold et al. 2016) and is now 65%. Above 18 M_⊙, the fraction of supernovae reduces to only 8%, which actually is in slightly better accord with the results from Smartt (2009, 2015), as compared to the previous study. Since the envelope masses have not changed appreciably and the explosion energies are also expected to be unchanged, the light curves will be unaltered. The fractions given in their Table 6 for supernovae above 12, 20, and 30 M_⊙ will also probably not change within the variations already seen for the different central engine characteristics.

The full evolution of over 4000 massive stars of solar metallicity in the mass range 12–60 M_⊙ has been studied using unprecedented resolution, both in zones per star and in number of stars within a given mass range. The mass loss rate was varied, and an important bug in the neutrino loss routine was repaired that caused a significant variation in the presupernova properties from those calculated by Sukhbold & Woosley (2014) and prior works. Our chief conclusions are as follows.

• 1.  
  The pattern of core compactness seen in previous studies (e.g., Sukhbold & Woosley 2014; Müller et al. 2016) is robust. The "noise" in these studies was not a consequence of inadequate zoning but reflected real variability. The range of variation and the location of peaks (Figure 8) in the new study are virtually identical to those seen by Müller et al. (2016), even though the new models use 4–10 times the zoning and a much smaller surface boundary pressure. The results are also qualitatively similar to those of Sukhbold & Woosley (2014), but important peaks in the compactness plot are shifted downward by about 10% in initial mass (Figure 9) due to the corrected neutrino loss rate. For the stellar physics used and spherically symmetric nature of the calculation, the variations and peaks seen in Figure 8 are now well determined.
• 2.  
  The large variation in core structure seen for stars between 14 and 19 M_⊙ by Sukhbold & Woosley (2014) is not completely random. For a larger set of models with finer spacing in initial mass, several branches of solutions emerge. This behavior was also seen by Müller et al. (2016). The branches apparently result from variations in the locations of the oxygen-burning shells in the presupernova star, which in turn "remember" the location of several carbon-burning shells of variable extent. Some noise is introduced by the fact that the carbon abundance and size of the carbon–oxygen core are not precisely monotonic with mass (Figure 7). This, in turn, reflects the operation of semiconvection at the boundary of the hydrogen and helium convective core. In the mass ranges 14–19 and 22–24 M_⊙, the presupernova core structure is more sensitive to events in the last year of a star's life (and sometimes the last hours) than to the star's initial mass. Presupernova structure in these mass ranges does not necessarily follow the Vogt–Russell theorem (Russell 1917; Vogt 1926). It may be more comparable to weather on Earth.
• 3.  
  Even accounting for these systematics, the branches of the solutions are still noisy. The level of noise is reduced if one characterizes the presupernova star by its helium core mass, not its starting mass or final total mass. The scatter is also reduced in the Ertl two-parameter characterization of core structure rather than the O'Connor-Ott compactness parameter.
• 4.  
  Given the large number of models, it is possible to give statistically meaningful results for the radius, luminosity, and effective temperature of supernova progenitors (Section 4.1).
• 5.  
  The mass distributions of neutron stars and black holes resulting from supernovae in the mass range studied are not greatly altered from our earlier surveys (Section 6). The average gravitational mass of neutron stars, including all masses of supernovae, is now 1.38 M_⊙. The average black hole mass, if only the helium core implodes, is 8.61 M_⊙. If the entire presupernova star collapses, the average black hole mass is 13.5 M_⊙. The fraction of stars above 9 M_⊙ that explode rather than collapsing is estimated to be 65%. The fraction above 18 M_⊙ is 8%. Nucleosynthesis and light curves are basically unchanged.

Binary mass exchange and rotation undoubtedly play key roles but have been neglected in this study. To first order, stars that end up with the same final helium core mass will have similar presupernova compactness and fates. Similar systematic variations and multivalued solutions are expected to persist. The statistical averages of compact remnant masses may vary, however, and certainly the light curves will be different. Our models are publicly available to those wanting to estimate outcomes based upon their own distribution of helium core masses, carbon–oxygen core masses, etc. In the future, we will consider rotating models, but this was principally a study of how resolution affects the solution to a well-defined, frequently studied problem.

Our results suggest a slightly different strategy to the study of presupernova evolution and supernova modeling than sometimes used in the past. Given the variation in outcome for stars of nearly the same mass, or the same mass with different codes, a statistically meaningful sample of models must be calculated before drawing strong conclusions about supernova mass ranges, remnant mass distributions, nucleosynthesis, etc. Historically, researchers have sometimes focused on the calculation of just a few masses, e.g., 15, 20, and 25 M_⊙, and sought to test the sensitivity to changes in physics in just those cases. Here we see that calculating a statistically meaningful sample may be just as important as getting the physics precisely right. The size of such a sample depends upon the need to resolve regions of rapid variability found with a given code and physics, but a minimum initial mass grid of 0.1 M_⊙ is reasonable within such regions.

Some of the models calculated here had merged carbon and oxygen convective shells at the end. Many did not. Others were only separated by a single thin zone from being coupled. Similar to the numerical artifact seen when helium burning is calculated with too little semiconvection, these mathematical bifurcations may not be real and might be overcome, or at least smoothed out (Alexakis et al. 2004), by increasing the convective overshoot or doing a multidimensional calculation (e.g., Meakin & Arnett 2007; Viallet et al. 2013; Arnett & Meakin 2016; Cristini et al. 2017; Jones et al. 2017). In general, linked convective shells give a structure with a smaller compactness parameter, ξ, and favor explosions (see, e.g., Collins et al. 2018). Further study with other representations of convective overshoot mixing and multidimensional codes are thus encouraged.

Semiconvection and overshoot mixing remain major uncertainties in studies of this sort and play important roles in determining the helium and carbon–oxygen core masses. How KEPLER treats semiconvection is described in Weaver & Woosley (1993). It is expected that other studies using different treatments will find results that differ in important details from those presented here. The overall pattern, multipeaked structure, and range of variability of the compactness parameter and other measures of core structure should persist, however. Clusters of solutions due to variable numbers and extents of the carbon and oxygen shells should also be a common feature. Further exploration is again encouraged.

All of the new presupernova models presented here are available online. Also online are plots of the convective histories of 800 models between 12 and 20 M_⊙ of the $\dot{M}$ set (see footnote 9). Other auxiliary presupernova models (e.g., those presented in Sections 2.1 and 2.2) are available upon request to the authors.

We thank Raphael Hirschi for a careful review of this work and for providing many useful suggestions. We also thank Bernhard Müller for helping us to process our new progenitor results through his semi-analytical explosion modeling and feedback on this manuscript. Thomas Janka and Thomas Ertl provided valuable insight into the supernova explosion models and their dependence on progenitor properties. We thank Todd Thompson and John Beacom for useful comments on the manuscript. We also appreciate the efforts of Gang Gao in discovering the neutrino bug in the KEPLER code. All numerical KEPLER calculations presented in this work were performed on the RUBY cluster at the Ohio Supercomputer Center (Ohio Supercomputer Center 1987). TS is partly supported by NSF PHY-1404311 to John Beacom. SW is supported by NASA NNX14AH34G. AH is supported by an ARC Future Fellowship, FT120100363, and, in part, by the National Science Foundation under grant No. PHY-1430152 (JINA Center for the Evolution of the Elements).