High Resolution Study of Presupernova Compactness
–– Oct, 2017
Preprint typeset using L A TEX style emulateapj v. 12/16/11
HIGH RESOLUTION STUDY OF PRESUPERNOVA CORE STRUCTURE
Tuguldur Sukhbold , S. E. Woosley , and Alexander Heger (Accepted – —, 2018) – Oct, 2017 ABSTRACTThe density structure surrounding the iron core of a massive star when it dies is known to have amajor effect on whether or not the star explodes. Here we repeat previous surveys of presupernovaevolution with some important corrections to code physics and four to ten times better mass resolutionin each star. The number of presupernova masses considered is also much larger. Over 4,000 modelsare calculated in the range from 12 to 60 M (cid:12) with varying mass loss rates. The core structure is notgreatly affected by the increased spatial resolution. The qualitative patterns of compactness measuresand their extrema are the same, but with the increased number of models, the scatter seen in previousstudies is replaced by several localized branches. More physics-based analyses by Ertl et al. (2016) andM¨uller et al. (2016) show these branches with less scatter than the single parameter characterizationof O’Connor & Ott (2011). These branches are particularly apparent for stars in the mass ranges14 - 19 M (cid:12) and 22 - 24 M (cid:12) . The multi-valued solutions are a consequence of interference betweenseveral carbon and oxygen burning shells during the late stages of evolution. For a relevant rangeof masses, whether a star explodes or not may reflect more the small, almost random differences inits late evolution than its initial mass. The large number of models allow statistically meaningfulstatements about the radius, luminosity, and effective temperatures of presupernova stars, their corestructures, and their remnant mass distributions.
Subject headings: stars: supernovae, evolution, black holes; nucleosynthesis; hydrodynamics INTRODUCTION
As a massive star below the pair-instability threshold( ∼
80 M (cid:12) ; Woosley 2017) evolves through its final stagesof nuclear burning, its central regions cool by neutrinoemission, become degenerate, and tend to decouple fromthe overlying layers and evolve as separate stars. Al-though never becoming completely detached except forthe lowest mass stars, the presupernova core takes on astructure similar to that of a white dwarf, mostly com-posed of iron, surrounded by a dense mantle of oxygenand intermediate mass elements. As numerous studieshave shown, the structure of this configuration, and es-pecially the rate at which the density declines outsidethe iron core is strongly correlated with the difficultyof blowing the star up (e.g., Burrows & Lattimer 1987;Fryer 1999). Recent studies have sought to capture thiscomplex structure in just one or two parameters thatmight predict, albeit approximately, whether a star witha given mass blows up or collapses to a black hole simplyfrom looking at one-dimensional models for stellar evo-lution (O’Connor & Ott 2011, 2013; Ugliano et al. 2012;Pejcha & Thompson 2015; Ertl et al. 2016; Sukhbold etal. 2016; M¨uller et al. 2016). Department of Astronomy, The Ohio State University,Columbus, OH 43210, USA, [email protected] Center for Cosmology and AstroParticle Physics, The OhioState University, Columbus, OH 43210, USA Department of Astronomy and Astrophysics, University ofCalifornia, Santa Cruz, CA 95064, [email protected] Monash Center for Astrophysics, Monash University, Vic3800, Australia, [email protected] School of Physics & Astronomy, University of Minnesota,Minneapolis, MN 55455, U.S.A. Department of Astronomy, Shanghai Jiao-Tong University,Shanghai 200240, P. R. China.
Sukhbold & Woosley (2014) discussed how the ad-vanced burning stages, especially convective carbon andoxygen burning, sculpt this structure, and provided a li-brary of presupernova stars consisting of 503 models inthe mass range 12 to 65 M (cid:12) to demonstrate the sys-tematics and its dependence on input physics. M¨ulleret al. (2016) prepared a larger grid of over 2,000 modelsbetween 10 and 32.5 M (cid:12) and corrected some errors inSukhbold & Woosley (2014). Their results showed finerstructure, but similar global systematics (see especiallytheir Fig. 6).More recently, Farmer et al. (2016) using the
MESA code, found that the choice of mass resolution, i.e., zon-ing, “ dominates the variations in the structure of theintermediate convection zone and secondary convectionzone during core and shell hydrogen burning, respec-tively ” and greatly affects the structure of presupernovastars. They found that a minimum mass resolution of ∼ .
01 M (cid:12) was necessary to achieve convergence in thefinal helium core mass at the ∼ ∼
30 % variations in the central electron fractionand mass locations of the main nuclear burning shells,and that a minimum of ∼
127 isotopes was needed to at-tain convergence of these values at the ∼
10 % level.Renzo et al. (2017) also used the
MESA code to explorethe sensitivity of massive star evolution due to varia-tions in the mass loss rate for stars with initial massesbetween 15 and 35 M (cid:12) . They found variations in thepresupernova core compactness parameter ( § ∼
30% depending upon the choice of the algorithm. In a lim-ited study of resolution in one model, they found roughly9% variation in final core compactness, smaller than thatresulting from uncertainties in mass loss prescription.In this paper, we present a new survey similar toSukhbold & Woosley (2014) and M¨uller et al. (2016), but a r X i v : . [ a s t r o - ph . H E ] M a r using much finer zoning, a greater number of models andseveral different mass loss rates. The key differences ofthe new survey are listed in §
2, including the correctionto an erroneous pair-neutrino loss rate used by Sukhbold& Woosley (2014). § § § (cid:12) and 22 - 24 M (cid:12) ranges (seealso M¨uller et al. 2016). In § §
6. Finally, in § CODE PHYSICS AND ASSUMPTIONS
With several important exceptions the code, physics,and input parameters used here are the same as inSukhbold & Woosley (2014) and M¨uller et al. (2016).The (solar) initial composition is from Asplund et al.(2009), as was used in M¨uller et al. (2016), but with anappropriate correction for the ratio N / N (Meibom etal. 2007). Convective and semi-convective settings, nu-clear reaction rates, opacities, and mass loss rates arethe same as in both of the earlier works. The mass rangestudied and lack of rotation are the same as well.Nuclear burning was handled as usual, using a 19 iso-tope approximation network until the central oxygenmass fraction declines below 0 .
03 and a silicon quasi-equilibrium network with 121 isotopes (Weaver, Zimmer-man, & Woosley 1978) thereafter. Numerous studies thatcompare this treatment with the energy generation andbulk nucleosynthesis obtained by “co-processing” with anetwork of several hundred to over 1,000 isotopes haveshown excellent agreement (e.g., Woosley et al. 2002;Woosley & Heger 2007), at least up to central oxygendepletion where the switch to the quasi-equilibrium net-work is made. By co-processing, we mean carrying alarge network in each zone using the same time step,temperature, and density, but not coupling the energygeneration from that large network directly to the itera-tive loop within a time step. After oxygen depletion, thequasi-equilibrium network is more stable, contains thesame weak interaction physics, and is roughly an orderof magnitude faster. An exception is the treatment ofoxygen and silicon burning in stars lighter than about11 M (cid:12) (Woosley & Heger 2015). For such light stars,it is important to follow neutronization in multiple off- center shells where the quasi-equilibrium approximationhas questionable validity and can be unstable. Such lightstars are not part of the present survey which starts at12 M (cid:12) . § • The models in this paper employ much finer massresolution (Fig. 1, Fig. 2, and Fig. 4). For thelightest stars considered (12 M (cid:12) ), the increase isroughly a factor of four. Since an effort was madeto keep fine resolution even in the larger stars, thefactor for the highest masses considered, around60 M (cid:12) , was closer to fifteen. In the most mas-sive models the total number of mass shells wasapproximately 16,000. Additional studies of in-dividual cases of zoning sensitivity in 15 M (cid:12) and25 M (cid:12) models ( § .
01 M (cid:12) everywhere, but were about tentimes smaller than that in the heavy element core(Fig. 4). It is this core that is most critical in de-termining the final presupernova core properties,such as compactness. Surface zoning was a fewtimes 10 − M (cid:12) , and all temperature and densitygradients were well resolved (Fig. 2). • A major, long standing coding error in the
KEPLER code was repaired. The error affected the axial vec-tor component of the pair neutrino losses. In par-ticular, the error resulted in the accidental zeroingof the second term involving C (cid:48) and Q − pair in theexpression (Itoh et al. 1996, eq. 2.1)Q pair = 12 [(C + C ) + n(C (cid:48) + C (cid:48) )]Q +pair (1)+ 12 [(C − C ) + n(C (cid:48) − C (cid:48) )]Q − pair (2)See Itoh et al. (1996) for the definitions of quan-tities. The consequence of this error was thatthe pair neutrino loss rate was underestimated byclose to a factor of two during carbon burningand somewhat less during later, higher tempera-ture stages. The error, introduced through an in-adequate checking of a routine provided by email,was included in 2001, and affected all KEPLER calcu-lations published through 2014. In particular it af-fected the often cited calculations of Woosley et al.(2002), Woosley & Heger (2007), and Sukhbold &Woosley (2014). Because some of the models fromSukhbold & Woosley (2014) were used by Sukhboldet al. (2016), it also affected the outcome of thatwork for masses larger than 14 M (cid:12) . The bug wasrepaired, however, in the works of Woosley & Heger(2015) on 9 M (cid:12) – 11 M (cid:12) stars and Woosley (2017)on pulsational pair instability supernovae. Mostrelevant to this present work, the bug was also re-paired for the work of M¨uller et al. (2016). Becauseof the strong sensitivity of neutrino and nuclear re-action rates to temperature, the effect of the bugwas a slight shift upwards in the burning temper-ature for the late stages of stellar evolution. Fora given model the change in presupernova struc-ture was not great and, as we shall see, withinthe “noise” of other uncertainties. For the lightestmodels ( <
14 M (cid:12) ) the difference is hardly notice-able, but it did systematically shift the outcome fora higher mass stars significantly downwards. Forexample, the peak in compactness that occurredfor Sukhbold & Woosley (2014) at about 24 M (cid:12) is shifted downwards in the work of M¨uller et al.(2016) and in the present work ( §
4) by about 3 M (cid:12) . • The surface boundary pressure is much less than inM¨uller et al. (2016), which significantly affects thefinal red supergiant (RSG) properties, but does notsignificantly alter the core structure. Whereas themass loss rates are varied in the present study, thenew calculations do not include Wolf-Rayet mod-els where the envelope is completely lost. All starsstudied here retained a substantial hydrogen enve-lope when they died.
The Effects of Zoning
One criterion for adequate zoning is that key variableslike temperature and density be well resolved, that is,that they do not vary greatly in going from zone “i” to“i+1”. Stellar evolution codes solve linear approxima-tions to non-linear differential equations. Fig. 2 showsthe resolution for a standard M ZAMS = 15 M (cid:12) presu-pernova star with 4225 zones. The figure gives the scaleheights in Lagrangian units (solar masses here) for thetemperature and density. Pressure and other derivedquantities, though not plotted, show similar variation.The scale height, e.g., for density is defined as M ρ =d m/ d ln ρ . In places, this quantity can abruptly be-come artificially small because of discontinuous changesin mean atomic weight at the edges of convective shells.Some of the prominent downward spikes in the iron coreare where partial photodisintegration has changed ¯ A ap-preciably. These abrupt changes in ¯ A are responsiblefor most of the discontinuous spikes in the figure. Largeupward spikes indicate regions of near constant tempera-ture or density. An especially large spike near the centerof the presupernova star reflects a temperature inversion(note that the scale height can be either positive or nega-tive; the absolute value is plotted). Where the derivativechanges sign and passes through zero, a spike results.The figure shows that there are roughly 100 to 1000zones per scale height everywhere in the star, exceptin the steep gradient at the edge of the helium core,where the actual density varies by six orders of mag-nitude over an interval of ∼ (cid:12) . Here the zoning isworst. In one location, the density changes by a max-imum of 30% from one zone to the next. Zone masseshere were ∼ × − M (cid:12) . Because this location moves inmass as the star evolves, still finer zoning would have sig-nificantly lengthened the calculation. We conclude thatthe relevant physical quantities are well resolved in thenew study. In fact, except at the edge of the helium core,they are over-resolved. This was done in order to explorethe sensitivity of outcome, e.g., the extent of convectiveshells, and to see if fine zoning alone would lead to a“converged” answer. (cid:12) ] − . − . − . − . − . l og m z o n e [ M (cid:12) ]
15 M (cid:12)
ZAMS
SW14This Work0 2 4 6 8 10 12Interior Mass [M (cid:12) ] − . − . − . − . − . l og m z o n e [ M (cid:12) ]
15 M (cid:12) presupernova
SW14This Work
Fig. 1.—
The masses of individual zones as a function of in-terior mass for a 15 M (cid:12) zero age main sequence star (top) andthe corresponding presupernova star (bottom). The gray curvesare the zoning for the earlier model of Sukhbold & Woosley (2014,SW14) which was very similar to that of Woosley & Heger (2007),Sukhbold et al. (2016) and M¨uller et al. (2016); the blue curvesare for the new survey reported here. The previous ZAMS modelhad 1,068 zones, the new one, 4,257. Stars are continually rezonedas they run to accommodate changing gradients in key quantities,but the zoning shown here did not vary greatly, with the exceptionof fine zones added at the base of the hydrogen envelope around4.4 M (cid:12) , and a few regions of finer zoning at what were once theboundaries of convective shells inside 2 M (cid:12) . The surface zoningwas ∼ g in all studies ( § Convergence assumes the existence of a well-defined so-lution to the stellar structure equations. As the numer-ical resolution is increased, if the subgrid physics (e.g.,convection) is coded in a zoning independent way, theanswer from a given calculation should approach this so-lution and give a constant answer. Certainly it is possi-ble for a stellar model to be inadequately zoned. Envi-sion Fig. 2 with zoning of 0.1 M (cid:12) . No one today wouldthink of trying to use just a few hundred zones in a pre-supernova model for a 25 M (cid:12) star, though substantialprogress was once made that way (Weaver, Zimmerman,& Woosley 1978). But given the power of modern ma-chines, how many zones are enough?Typical
KEPLER models by Sukhbold & Woosley (2014)and M¨uller et al. (2016) used about 1,200 zones, roughlyindependent of the star’s total mass, though the zon- − − − − l og M a ss [ M (cid:12) ]
15 M (cid:12) presupernova M ρ M T m zone (cid:12) ] − . − . − . − . − . − . − . − . l og | d ρ / ρ |
15 M (cid:12) presupernova ¯ A Fig. 2.—
Density and temperature resolution in the same 15 M (cid:12) presupernova model shown in the lower panel of Fig. 1. The toppanel shows the scale heights ( M x ≡ d m/ d ln x ) expressed in solarmasses for the density ( M ρ ; grey curve) and temperature ( M T ;red curve). The blue dashed curve shows the zoning. The bot-tom panel shows the actual variation in density between zones andmean atomic weight, ¯ A , plotted also as a function of interior mass.Discontinuities in ¯ A occur at the boundries of active and fossilconvective shells causing abrupt changes in density that the codeattempts to resolve. The large spike in M T at 0.15 M (cid:12) reflects atemperature inversion at the star’s center (where d T/ d m is nega-tive) and a small region of nearly constant temperature boundingit. Both panels show that temperature and density gradients arevery well resolved throughout the star, except perhaps in the verysteep density gradient at the edge of the helium core (4.33 M (cid:12) ),where changes of up to 30% occur for the density in a few zones. ing was by no means uniform. Zones were concentratedat the center, where temperature-dependent burning re-quired fine resolution, in the steep density gradient atthe base of the red giant convective envelope, and atthe surface. Continuous rezoning kept gradients in den-sity, temperature, and zone radius well resolved. For alarge range of masses, characteristics of the presupernovastar, such as its helium, carbon, and iron core masses,its luminosity and radius, and its “compactness” variedsmoothly with mass and were consistent with numerouspast studies. In some mass ranges, however, especially inthe range 14 M (cid:12) to 19 M (cid:12) , the compactness parameteroften varied wildly. These variations were attributed bySukhbold & Woosley (2014) to the interaction of convec-tive carbon and oxygen burning shells. The character- istics of those burning shells were irregular in location,extent, duration, and intensity, and so, for some masses,was the final presupernova core structure. But mightsome of that variation also be attributed to inadequateresolution? Would a set of finer-zoned models show lessvariation, or reveal structure in the noise?Most of this paper is about the results of repeating thesurvey of Sukhbold & Woosley (2014) with finer zoning(Fig. 1). In this section, however, we briefly examine, forjust four masses, the question: “What is enough?” Inpart, the answer must have a pragmatic aspect. If run-ning with 3,000 and 30,000 zones gives answers that differin some important quantity by 5 %, but if changing thestellar mass by 0.01 M (cid:12) , or some bit of uncertain stel-lar physics by a small fraction of its error bar alters theanswer by 50 %, perhaps 3,000 zones is “good enough.”Resources would be better spent studying the variationsthat depend on these other variables.There is a deeper issue though when the answer is de-terministic, but unpredictable. It is an inherent aspectof chaos that tiny changes upstream produce large differ-ences in outcome. As the noted chaos theorist EdwardLorentz once said “ Chaos exists when the present de-termines the future, but the approximate present doesnot approximately determine the future ” . As noted bySukhbold & Woosley (2014) and M¨uller et al. (2016) and,as will be explored in greater depth here, for some rangesof mass, presupernova evolution is like that. Seeminglyminuscule changes in initial conditions may determinewhether a star in an important mass range explodes asa supernova or collapses to a black hole. How differentcould two outcomes be?Table 1 gives the final presupernova properties of foursets of models with varying resolution for initial massesof 15.00, 15.01, 25.00, and 25.01 M (cid:12) . The spatial res-olution ranges from about 1,200 to 16,000 mass shells.The most finely resolved models used roughly twice asmany timesteps to reach the presupernova stage as theleast resolved models. The “compactness” parameter, ξ . (O’Connor & Ott 2011) is inversely proportional tothe radius enclosing innermost 2.5 M (cid:12) of the presuper-nova core, while M and µ (Ertl et al. 2016, furtherdiscussed in § s = 4 k B and the derivative of mass at that location evaluated overa mass interval of 0.3 M (cid:12) . The motivation for using theseparameters is discussed in § (cid:12) case in Fig. 1 corre-sponds to the moderately zoned case ( ) in Table1. Initial zoning in other 15.00 M (cid:12) and 15.01 M (cid:12) caseswas scaled at all locations in the star by roughly fac-tors of two. For the 25.00 M (cid:12) and 25.01 M (cid:12) modelsthe number of zones was proportionately larger. Thatis Model had roughly 25 /
15 as many zones asModel , etc. There was no Model , in-stead there was a Model which had about halfthe zoning of . The continuous rezoning parame-ters were set so as to preserve, within about 10 %, theinitial zoning. The number of zones given in Table 1 isfor the presupernova star.Two closely adjacent masses were studied at 15 M (cid:12) mpe2013.org/2013/03/17/chaos-in-an-atmosphere-hanging-on-a-wall TABLE 1Presupernova properties with varying resolution and mass M ZAMS N zones N steps M presSN L preSN R preSN T eff , preSN M He M Fe ξ . M µ [M (cid:12) ] [M (cid:12) ] [10 ergs s − ] [10 cm] [K] [M (cid:12) ] [M (cid:12) ] [M (cid:12) ]
874 38558 16.047 7.88 8.56 3505 8.342 1.609 0.303 1.915 0.103
876 38416 16.138 7.88 8.56 3505 8.269 1.606 0.301 1.911 0.102
Note . — All quantities are measured at the presupernova stage. N zones is slightly larger at the beginning of calculation. and 25 M (cid:12) in order to compare the effects of fine zon-ing with slight variations in any other parameter of theproblem. One might have chosen instead to vary timestep, overshoot, the rate for C( α, γ ) O, semiconvec-tion, mass loss, or rotation. Although no substitutefor a full survey of such dependencies, varying the massslightly in a situation where the solution is chaotic mightbe expected to send the calculation down one path oranother. Stated another way, if varying the mass by lessthan 0.1 % results in an answer significantly differentthan one obtained by increasing the zoning by a factor oftwo, perhaps there is no reason to use much finer zoningfor that mass. On the other hand, if all models give thesame answer, one might have confidence in the unique-ness of that model, at least for an assumed set of stellarphysics.The two choices, 15 M (cid:12) and 25 M (cid:12) , illustrate thesetwo possibilities nicely. In all cases the observable prop-erties of the presupernova star, its luminosity and radius,and hence its effective temperature, are well determined.Changing opacities, composition, etc., would certainlycause variation, but the bulk observables of our mod-els are well determined. There is a weak, but noticeabletrend to produce smaller presupernova stellar masses andlarger helium core masses when the zoning is finer. Aswill be discussed later, this is a consequence of variationin the time spent as a red or blue supergiant, which, inturn, is sensitive to an uncertain treatment of semicon-vection. A small number of episodic mixing events atthe edge of convective hydrogen burning core can tem-porarily “rejuvinate” the core and by affecting the grav-itational potential of hydrogen burning shell, eventuallycausing the transition from blue to red to occur at dif-ferent times. This trend does not emerge because themodels are converging with higher resolution, but ratheras a result of enhanced mixing due to the uncertain treat-ment of convection physics.The measures of presupernova core structure on the other hand vary by factors of two for the 15 M (cid:12) models,but scarcely at all for the 25 M (cid:12) models. The reasonsfor the difference will be discussed in §
5, but we notethat any reasonable zoning suffices to get the structureof the 25 M (cid:12) correct, but no reasonable zoning showsconvergence for the 15 M (cid:12) model. Not knowing before-hand what the outcome will be for an arbitrary mass,we used here the maximum zoning that time and thedesire to survey thousands of models allowed. This ismost like the “D” models in Table 1. There is little mo-tivation, however, for carrying out still finer resolutionstudies when a trivial change in mass, 15.00 M (cid:12) to 15.01M (cid:12) , causes more variation than an order of magnitudeincrease in zoning.
The Effects of Networks
Uncertainty in nuclear energy generation and neutron-ization during the various episodes of burning can beanother source of variability in presupernova models. Inthis section the effect of using either the standard 19-isotope approximation network plus silicon quasiequilib-rium ( §
2) or a much larger nuclear reaction network isexplored. It should be noted that while the standardnetwork carries the abundances of only 19 species, it hasthe power of a network roughly twice that size since re-actions through trace species like Na, Al, etc. arecarried in a steady-state approximation (Weaver, Zim-merman, & Woosley 1978). Reactions coupling Feto Ni are also included, as are approximations toenergy generation by the CNO- and pp-cycles and to N( α, γ ) F( e + ν ) O( α, γ ) Ne .The large network used for comparison here is drawnfrom an extensive nuclear data base maintained and usedfor nucleosynthesis studies for decades (e.g., Woosley etal. 2002). Key reaction rates such as N( p, γ ) O, 3 α , In the code, N is converted to Ne at an appropriate rateconserving mass. C( α, γ ) O, the ( α, γ ) and ( α, p ) reactions on α − nucleiand the weak interaction rates are the same in both net-works. The network is “adaptive” (Rauscher et al. 2002)in the sense that isotopes are added or subtracted at eachstep to accommodate the reaction flow (e.g., there is noneed to include a detailed network for the iron groupduring hydrogen burning). Usually the settings on thenetwork size are quite liberal, resulting in from 700 to1,200 isotopes being carried in the early and late stagesof evolution respectively (Woosley & Heger 2007). Here,because we are only interested in stellar structure andnot, e.g., the s -process, a smaller number was carried,typically 300 isotopes extending up through the elementSelenium. All of the iron-group isotopes in the quasi-equilibrium network were included, and many more, sothat core neutronization during and after silicon burningwas equivalently calculated.Network sensitivity was examined for four differentmasses of star, 15.00, 15.01, 25.00, and 25.01 M (cid:12) . Twostudies were carried out to test network sensitivity. Ta-ble 2 edits the evolution of the 15.00 star, similar toModel in Table 1, at different times in its evolu-tion using 3 approximations to the energy generation andneutronization. The Approx case corresponds to usingthe 19 isotope approximation network and quasiequilib-rium network, as is standard in the rest of this paper.The
Coproc case carries the large network of about 300isotopes along with the approximation network in “co-processing mode”. That is, the approximation and quasi-equilibrium networks are used for energy generation inthe stellar structure calculation, but the network is alsocarried along in passive mode. The same time step isused at the same temperature and density to evolve, zoneby zone, a much larger number of isotopes in parallel withthe approximation network. Sub-cycling is done withina given time step to follow the abundances of trace iso-topes accurately. Coprocessing is only performed in thecode up to the point where the transition to quasiequi-librium occurs in a zone, typically at oxygen depletion X ( O = 0.03). The output from the large co-processingnetwork is used to continually update the electron molenumber, Y e , which feeds back into the structure calcula-tion. Because Y e evolves slowly, no iteration is required.The energy generation from the large network can also beused 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 cal-culation throughout the entire life of the star, includingsilicon burning and core collapse. This might be deemedthe most accurate approach, but it is slow, requiring anorder of magnitude more computer resources, and proneto instability during silicon burning because of stronglycoupled flows.In Table 2, up until oxygen depletion, two differentnumbers are given for the co-processing run (
Coproc ).The first is the energy generation from the approximationnetwork, which is used to calculate the stellar model, andthe second is the output from the passively carried bignetwork. Hydrogen burning and depletion correspond tocentral hydrogen mass fractions of 0.4 and 0.01. Heliumburning and helium depletion are when the central he-lium mass fraction is 0.5 and 0.01. Carbon “ignition” isactually evaluated during the Kelvin-Helmholtz contrac-tion between helium depletion and real carbon burning when the central temperature is 5 × K, and carbondepletion is when the central carbon bass fraction is 1 %.Similarly, “oxygen ignition” is when the central temper-ature during Kelvin-Helmholtz contraction is 1 . × K and oxygen has yet to burn, and oxygen depletion iswhen the central oxygen mass fraction is 0.03. Silicon ig-nition is at 3 . × K, and silicon depletion is when thesilicon mass fraction is 1 %. Presupernova is when thecore collapse speed exceeds a maximum of 900 km s − .The table shows near exact agreement in energy gen-eration and Y e calculated using all three approaches upuntil at least oxygen ignition. The slight differences inenergy generation at carbon ignition do not matter be-cause neutrino losses (far right column) dominate at thetime examined. By oxygen depletion though, things arestarting to mildly diverge. This divergence has threecauses. One is the electron capture that goes on in thelate stages of oxygen burning and starts to apprecia-bly affect Y e . This change is not followed by the Ap-prox network. Carrying Y e in the co-processing run ad-dresses most of this divergence; the values of Y e in theco-processing and full network runs are nearly the same.Another effect, more difficult to disentangle, is the factthat, at late times, the stellar structure in the differentcalculations starts to diverge. It diverges, in part, ofthe different nuclear physics, but even more because ofthe chaotic nature of the evolution. As was seen in thezoning study ( § Y e agree very well and the energy gen-eration differs by a factor of two, which is mostly due tothe higher temperature in the Approx case and the verysensitive temperature dependence of the silicon burningreactions.These small differences probably cause very littlechange in the presupernova model because the star hasa certain amount of fuel to burn and the total energyrelease is set by known nuclear binding energies. Thedifferent burning rates only affects modestly the temper-ature at which the silicon burns in steady state (Woosleyet al. 2002). We conclude that the approximation plusquasiequilibrium approach is adequate for our survey andintroduces errors that are small compared with othervariations that occur when other uncertain quantities,or even the zoning change.Table 3 examines the effect of using either the largenetwork or the approximation network plus quasiequilib-riumin to study stars of 15.00, 15.01, 25.00, and 25.01M (cid:12) , the same masses considered in the zoning sensitiv-ity study ( § × swas placed on the time step as opposed to 1 × s in thesurvey and zoning sensitivity study. Here, the emphasisis on comparing runs with different networks, not on thehighest spatial and temporal resolution. Consequentlythe 15.00 model in Table 3 is not directly comparable TABLE 2Central properties for different networktreatments for 15.00 M (cid:12) model network log T c log ρ c Y e , c log S n log | S ν | [K] [g cm − ] [ergs g − s − ]H burning Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Note . — S n is the energy generation rate excludingneutrino-losses, and is negative during the presupernovastage due to photodisintegration. Values in parentheses de-note log | S n | . The two values of Y e , c and S n given forco-processing calculations until core oxygen depletion (un-til quasi-equilibrium) are form approximation network (top)and co-processed big network (bottom). with the one in Table 1.Table 3 gives, in addition to the presupernova zoningand time steps, the masses of the iron, carbon-oxygen,and helium cores, the central mass fraction of carbonevaluated just prior to carbon ignition, the central valueof the electron mole number, Y e in the presupernova star,and the descriptors of presupernova core structure, ξ . , µ , and M . Though not given, the presupernova masses, radii, luminosities, and effective temperatures are in thesame range as for the same mass stars in Table 1.Based solely on the changes in the core structure andiron core mass for the 15.00 and 15.01 M (cid:12) models, onemight conclude that using a large network made a sub-stantial difference and all future runs should take muchgreater care with the nuclear physics. Though the mod-els do indeed differ, this would not be a valid inference.As shown in § anything for models in some mass ranges cangive qualitatively different answers for core structure.The differences between runs with large and small net-works is within the range of differences resulting from in-creased or decreased zoning (Table 1) or a small changein the star’s mass. More telling is the fact that the cen-tral carbon abundances, the helium and carbon-oxygencore masses, and most critically the central value of Y e all agree very well. The approximation network is doinga fine job representing the nuclear physics. It is just thatother, uncontrollable factors introduce large variations inthe core structure of a 15 M (cid:12) , and to a lesser extent 25M (cid:12) star. We conclude that the 19-isotope approxima-tion network and quasiequilibrium hypothesis are quiteadequate for surveys like this. Instead, focus should beplaced on studying the impact of stellar physics on thestatistical properties over the entire mass range. Boundary Pressure
Especially during its post-main sequence evolution, theradius of a model star is sensitive to surface boundaryconditions. Different codes treat the surface in differ-ent ways. Many use a boundary condition on pressure ordensity that sometimes includes a reduction of gravity bythe Eddington factor and by the inertia term from windacceleration (e.g., Appenzeller 1970; Heger 1998; Paxtonet al. 2015). Others solve the wind equation for condi-tions at the sonic radius (Grassitelli, private communi-cation, 2018). Still others fit the surface to a stellar at-mosphere of varying complexity (e.g. plane-parallel grayor tables) (e.g. Paxton et al. 2011; Chieffi & Straniero1989).
KEPLER uses a constant boundary pressure, P bound that does not vary during the evolution. The advantageof such an approach is its simplicity, but care must betaken that P bound is not so large as to greatly alter thesolution.Ideally, P bound only influences the structure of a tinymass near the surface. No reasonable value of P bound af-fects main sequence evolution, for example, because thethe pressure gradient there is very steep. The gradientbecomes shallower during the RSG phase though, andeven slightly different radii there can significantly alterthe mass loss rate. The timing and development of thesurface convection zone is also affected. Both can have asignificant affect on the star’s location in the HR-diagramand, to a lesser extent, the helium core mass and presu-pernova structure. Traditionally, studies with KEPLER have used P bound = 50 to 100 dyne cm − (e.g., Woosleyet al. 2002; Woosley & Heger 2007; Sukhbold & Woosley2014). An exception is the work of M¨uller et al. (2016),where a much larger value, ∼ − was used.In the present study a value of 50 dyne cm − is em-ployed. Is this low enough? To explore the sensitivityof results to boundary pressure, a standard 15.00 M (cid:12) ,Model S15.00D in Table 1, was calculated using a vari- TABLE 3Presupernova properties with varying network size M ZAMS network N zones N steps M He M CO 12 C ign . M Fe ξ . M µ Y e , c [M (cid:12) ] [M (cid:12) ] [M (cid:12) ] [M (cid:12) ] [M (cid:12) ]15.00 Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Approx
Coproc
Full
Note . — All quantities are measured at the presupernova stage, except C ign ( § N zones is slightlylarger at the beginning of calculation. Y e , c is evaluated at the center. ety of boundary pressures from 10 to 6400 dyne cm − maintained throughout the evolution. The resulting HRdiagrams and radial histories are shown in Fig. 3. For P bound < ∼
500 dyne cm − the trajectory in the HR-diagramis essentially identical for six different choices of P bound .The presupernova radius and luminosity do not vary. Forlarger values of P bound though, especially P bound > ∼ − significant variation is found.Even the models with small surface pressures showsome variation in final properties. For example, for P bound between 10 and 400 dyne cm − the final star massvaried between 12.59 M (cid:12) and 12.61 M (cid:12) , while the he-lium core mass varied between 3.61 M (cid:12) and 3.67 M (cid:12) .Some of this variation is a consequence of the irreduciblenoise in running any model twice with even small varia-tions in the physics or mass, but part could be a residualsensitivity to P bound .Closely related to the treatment of surface bound-ary conditions is surface zoning. In our recent studies(Sukhbold & Woosley 2014; M¨uller et al. 2016), rezon-ing was allowed to continue all the way to the surfaceof the star throughout its evolution. This often resultedin coarse zoning near the photosphere during the RSGstage. Fine zoning on the main sequence was lost to de-zoning when the temperature and density gradients wereshallow. See Fig. 1 for a sample comparison with theolder model. In the new models rezoning was not al-lowed in the outer 250 zones, typically 0.3 M (cid:12) . Thusthe fine surface zoning shown in Fig. 1 was maintainedthroughout the evolution for all models. Typical zoningnear the surface is less than 5 × − M (cid:12) .The use of finer surface zoning and a reduced boundarypressure generally gave larger radii and increased massloss during the RSG stage, especially as compared withM¨uller et al. (2016). Because the helium core mass wasnot appreciably affected though, the final core compact-ness, nucleosynthesis, and remnant masses are not al-tered by this change. THE NEW SURVEY
Three grids of models were calculated for initial mainsequence masses in the range 12 through 60 M (cid:12) usingthe zoning described in § M N ). This re-sults in so much of the envelope being lost by stars aboveabout 27 M (cid:12) , that the residual envelope, typically a fewsolar masses, becomes unstable and develops density in-versions of more than an order of magnitude. We sus-pect that massive stars with extremely large radii andlow mass envelopes become unstable in nature at thispoint and lose their remaining hydrogen rapidly (e.g.,Sanyal et al. 2015, 2017). Mass loss is not the focus herethough. In order to have a greater range of helium coremasses, a second set of models used one-half of the stan-dard mass loss ( ˙ M N /
2) appropriate for solar metallicityand encountered no instability up to 40 M (cid:12) . Finally,a third set was calculated with one-tenth the standardmass loss ( ˙ M N /
10) to provide a sparse mass grid be-tween 12 M (cid:12) and 60 M (cid:12) . Considering that mass lossprescriptions are probably uncertain (Renzo et al. 2017;Beasor & Davies 2018), the factors of two and ten couldreflect a sensitivity study, or they might be appropriatefor stars with reduced metallicity.Altogether 1,499 models were calculated in the massrange 12 to 27 M (cid:12) using the standard mass loss rate and2,799 models from 12 to 40 M (cid:12) using half that value.The grid spacing was 0.01 M (cid:12) for both sets. The thirdset with one-tenth the standard mass loss rate, consistedof only 49 models between 12 and 60 M (cid:12) with a meshsize of 1 M (cid:12) . Lower metallicity stars were not consid-ered but, except for zero, metallicity would have affectedthe answer chiefly by changing the mass loss and hencehydrogen envelope mass of the presupernova star. As weshall see, the helium core mass is not entirely indepen-dent of the remaining envelope mass, but it is insensitive.The code was compelled to spend a much greater num-ber of time steps during hydrogen and helium burning byimposing a limit on the maximum time step of 10 s. Aminimum of roughly 10,000 steps was thus spent burninghydrogen on the main sequence. This helped to weakenthe impact of random semiconvective mixing events thataffected the final envelope mass and, to a smaller extent,the helium core mass. It also improved the accuracy ofthe treatment of convection and burning as “split” oper-ations (convection and burning are not implicitly coupled . . . . . . . log T eff [ K ] . . . . . . . . . l og L [ e r g ss − ] P bound . : [dyne cm − ] . . . . . .
00 12 .
25 12 .
50 12 .
75 13 .
00 13 .
25 13 . Time [Myr] R e ff . [ c m ] Fig. 3.— (Top panel:) The HR-diagram for variations of Model15.00D (Table 1) in which the surface boundary pressure is variedfrom 10 dyne cm − to 6400 dyne cm − in roughly factor of twosteps. The inset shows that the evolution is insensitive to thisboundary condition until the star becomes a red supergiant. Themain figure shows nearly identical results for pressures < ∼
500 dynecm − , but significant variations for larger values. (Bottom panel:)Radius as a function of time during the final evolution of the samemodels. Only the last 1.4 My when the stars become RSGs isshown. Models with boundary pressures less than 500 dyne cm − have nearly identical radii as supergiants until very close to the endwhile those with bigger pressures vary. The time spent as RSG’sis similar. in the code). Very little difference was noted in a few testcases when this maximum step was increased to 5 × s. Minor changes to improve the stability of the codewhen transitioning to silicon quasiequilibrium and theconvergence criteria were also incorporated. These hadinsignificant effects on the outcome, but improved codeperformance. Models for the main survey were all calcu-lated on identical processors using the same version of thecode and compiler so as to reduce any possible noise in-troduced by machine architecture (Ohio SupercomputerCenter 1987).One major goal of this work was to explore the effectof resolution on calculations of presupernova structure.The number of zones used in the main survey was typi-cally 4 , × ( M ZAMS /
12 M (cid:12) ). Stars more massive than
10 20 30 40 50 60 M ZAMS [M (cid:12) ]200040006000800010000120001400016000 N u m b e r o f Z o n e s p r e s u p e r n o v a H e − c o r e C O − c o r e F e − c o r e ˙ M N ˙ M N / M N / Fig. 4.—
Total number of zones in the presupernova star as afunction of main sequence mass for the model sets with standardmass loss rate ˙ M N (blue), with reduced rate of ˙ M N / M N /
10 (black). In order to maintain fine reso-lution in stars of increasing mass, the number of shells becomesgreater. The lower curves show the zoning inside the helium core,the carbon-oxygen core and iron-core. These new models retainmore zones only in the iron-core alone as the earlier studies did inthe whole star.
30 M (cid:12) thus had over 10,000 initial mass shells. In allcases, the maximum mass zone anywhere in the star was,at all times, (cid:46) .
01 M (cid:12) , but zoning was by no meansuniform. The initial zoning (Fig. 4), which continuedto characterise most of the star until the presupernovastage, had much finer zoning in the helium core, down toabout 0.001 M (cid:12) at the stars center.Except as noted above for hydrogen burning, the timestep criteria were not varied from the previous studies.Nevertheless the finer resolution of the new runs did re-sult in taking more time steps as smaller zones experi-enced mass loss or were cycled through the large densitycontrast at the base of the hydrogen convective shell.Abundances within a given zone convectively coupledto many other zones also required more time steps tochange. The typical number of steps for a new calcula-tion is 45,000, about twice the earlier studies ∼ M N , ˙ M N / M N /
10 are available through the Harvard-Dataversearchive . RESULTS
Observable Properties
The observable properties of the new models are sum-marized in Fig. 5, Fig. 6, and Table 4. Final (presuper-nova) masses are sensitive both to the uncertain massloss rate prescription (Nieuwenhuijzen & de Jager 1990),and to the uncertain history of the stellar radius. The re-sults vary considerably. For lower mass stars, the spreadin final mass is smaller since the star only loses a smallamount of mass. For bigger stars, though, a large frac-tion of the hydrogen envelope may be lost. Most ofthe mass loss occurs during helium burning, and mostof that, during the late stages when the central heliummass fraction is less than 0.5. http://doi.org/10.7910/DVN/VOEXDE R . in the Nieuwenhuijzen & de Jager (1990) for-malism, i.e., the amount of time the helium burning starspends as a blue supergiant (BSG) or RSG. It is wellknown that the semi-convective mixing and overshootmixing play a key role in determining the ratio of timespent as each (Lauterborn et al. 1971). In the presentcalculations, small changes in individual semiconvectivemixing episodes affect the timing of the development of adeep surface convection zone, and hence the movement ofthe star to the red. With more semiconvection, the starspends a greater fraction of its helium burning lifetimeas a BSG and hence loses less mass, ending its life witha larger hydrogen envelope still intact. RSGs have theconverse behaviour. Less semiconvection, e.g., strictlyLedoux convection, favors a longer time as a RSG andthus reduces the threshold for a massive star to lose itsenvelope. Though the spread in Fig. 5 looks large, theuncertainty is usually a small fraction of the total masslost. For example a 25 M (cid:12) star with the standard massloss rate might end up as a 12 or 15 M (cid:12) star. That is itmight lose 10 to 13 M (cid:12) , a range of about 25 %. Lowermass models have a smaller variation.The final masses for stars calculated using ˙ M N /2, ofcourse, are larger. Bigger stars on the main sequencethen retain more envelope by the time they die, and wewere thus able to determine the core structure of moremassive stars. For the most massive models consideredthe final hydrogen envelope mass was roughly 4 M (cid:12) forthe ˙ M N models and 5 M (cid:12) for the ˙ M N / . × cm and became unstable, in the KEPLER code,due to recombination. In these cases, the envelope couldonly be retained on the star by the addition of a large,unrealistic boundary pressure that resulted in gross den-sity inversions in the outer envelope. We suspect thatthis is a real instability in nature, that once the envelopemass decreases below a critical value in a very massivestar and the radius extends beyond 10 AU, the remain-ing envelope is lost through non-steady processes (Sanyalet al. 2015, 2017). The luminosity in the hydrogen en-velope 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 andreduced mass loss, different choices of boundary pressurecan lead to significant differences in the final mass. Theboundary pressure used in M¨uller et al. (2016), for exam-ple, several thousand dyne cm − , was larger than whatwe now believe reasonable and much larger than usedin the present study, 50 dyne cm − . This probably ac-counts, at least partly, for the larger final masses foundin the M¨uller et al. (2016) study and accounts for theirability to study higher mass stars using a single mass lossprescription.This boundary pressure does not appreciably affect thecompactness of the cores for stars of given main sequencemasses, however. This is because the large spread in finalstar masses (Fig. 5) is not reflected in the helium coremass ( §
12 14 16 18 20 22 24 26 M ZAMS [M (cid:12) ]1214161820222426 S t a r M a ss [ M (cid:12) ] He − ign . He − burn . He − depl . preSN .
15 20 25 30 35 40 M ZAMS [M (cid:12) ]10121416182022 F i n a l M a ss [ M (cid:12) ] ˙ M N ˙ M N / Fig. 5.—
Stellar masses during helium burning and for presu-pernova stars. The upper panel shows results for standard massloss ( ˙ M N ) at helium ignition (pink) when 1 % of the helium hasburned to carbon; helium burning (orange) when central heliummass fraction is 50 %; helium depletion (red) when the central he-lium mass fraction is 1 %; and the presupernova star (blue). Thepresupernova masses and helium depletion masses are almost indis-tinguishable. Most of the mass loss, and most of the dispersion inthe presupernova mass, arises due to radius expansion after sub-stantial helium has already burned. The lower panel shows thelarger presupernova masses expected when the mass loss rate isreduced by a factor of two. Models with reduced mass loss gener-ally have lower dispersion, but both show the dispersion increasewith increasing initial mass (thus increasing amount of mass lost).The presupernova helium core masses show much less variation (seeFig. 7). envelope mass, the spread in final masses seen in Fig. 5does not affect our major conclusions.Fig. 6 and Table 4 give the final luminosities and ef-fective temperatures of our presupernova models. Theluminosity of the presupernova star depends chiefly onthe helium core mass and is approximately given by L preSN ≈ . × ( M He / (cid:12) ) . ergs s − . The ef-fective temperature, essentially bounded by the Hayashilimit remains pegged at close to ∼ M . . Except for a few presupernova stars, most RSGswill be observed during their helium burning stage wheretheir effective temperatures will be hotter. Our calcula-tions show a systematic ∼
200 K offset between the ef-fective temperatures at helium depletion and presuper-1 . . . . . .
56 log T eff , preSN [K]38 . . . . . . . . l og L p r e S N [ e r g ss − ] M Z A M S [ M (cid:12) ]
12 14 16 18 20 22 24 26 M ZAMS [M (cid:12) ]335034003450350035503600365037003750 T e ff [ K ] p r e s up e r n o v a H e − d e p l e t i o n Fig. 6.— (Top panel:) Location in the HR-diagram for the presu-pernova stars that used the standard mass loss rate ( ˙ M N ). (Bottompanel:) T eff at helium depletion (red) and for the presupernovastars (blue). The presupernova radii vary almost linearly from5 × cm to 10 × cm as the mass increases from 13 to 26M (cid:12) . The stars all die with a nearly constant effective temperaturebetween 3400 and 3600 K. Observations would select RSGs prior tohelium depletion, so the helium-depletion curve (red) in the lowerpanel is a lower bound to what is likely to be observed for T eff . nova stars. Earlier in helium burning the temperaturewill be even hotter and thus T eff ≈ ,
550 K should bea lower bound to what is seen. Larger values are alsoexpected for sub-solar metallicity and lower values forpresupernova stars. These numbers are in good accordwith measurements of the hottest RSGs (Davies et al.2013; Levesque et al. 2005).
Core Masses
Fig. 7 shows the helium, carbon-oxygen (CO), andiron-core masses for all our presupernova models with˙ M N and ˙ M N /
2. The dispersion in helium and CO coremasses for a given main sequence mass is small, much lessthan the hydrogen envelope masses inferred in Fig. 5.This implies that the helium and CO core masses areonly slightly affected by the assumed mass loss rate andare probably not very sensitive to metallicity. The he-lium core mass as a function of main sequence mass isapproximately M He ≈ .
46 ( M ZAMS / . which can
15 20 25 30 35 40246810121416 C o r e M a ss [ M (cid:12) ] C O C o r e s H e C o r e s ˙ M N ˙ M N /
215 20 25 30 35 40 M ZAMS [M (cid:12) ]1 . . . . . . . I r o n C o r e M a ss [ M (cid:12) ] ˙ M N ˙ M N / Fig. 7.—
Helium core, CO core, and Fe core masses for all pre-supernova stars from sets with ˙ M N and ˙ M N /
2. Despite significantvariations in mass loss (Fig. 5), the final helium and carbon-oxygencores are well determined by the star’s initial mass and a standardchoice of stellar physics. Note multiple branches for the iron coremass below 19 M (cid:12) . be combined with the previous relation in § L preSN ≈ . × ( M ZAMS /
20 M (cid:12) ) . ergs s − . Sincepresupernova core compactness is chiefly a function ofhelium core mass (Sukhbold & Woosley 2014), these re-sults suggest a near universal dependence of presuper-nova core structure on initial mass, provided mass lossdoes not remove the entire hydrogen envelope.Iron core masses increase, on the average with increas-ing stellar mass reaching maximum of about 2.0 M (cid:12) forthe most massive stars studied ( >
40 M (cid:12) ). Still largeriron cores, up to about 2.5 M (cid:12) characterize more mas-sive stars in the pulsational-pair instability range (70 -140 M (cid:12) ; Woosley 2017). For stars below 23 M (cid:12) , theiron core mass is markedly multi-valued for stars withnearly the same initial mass. This reflects the operationof multiple shells of carbon and oxygen burning as willbe discussed further in §
5. The two major branches ofiron core masses below 20 M (cid:12) , which are most of thestars that leave neutron star remnants might result inbimodality in the neutron star mass function.
Core Structure TABLE 4Properties of three main model sets M ZAMS M preSN M He M CO log L preSN T eff , preSN T eff , He dep R preSN [M (cid:12) ] [M (cid:12) ] [M (cid:12) ] [M (cid:12) ] [ergs s − ] [K] [K] [10 cm]˙ M N (binned by 1 M (cid:12) )12 11.19 3.35 2.20 38.39 3440 3670 4.9713 11.72 3.79 2.53 38.47 3420 3610 5.4914 12.22 4.20 2.85 38.54 3410 3570 5.9715 12.73 4.61 3.19 38.60 3405 3545 6.4116 13.23 5.01 3.53 38.65 3405 3535 6.8017 13.74 5.42 3.87 38.69 3410 3535 7.1518 14.12 5.84 4.24 38.73 3420 3530 7.4519 14.60 6.24 4.60 38.77 3425 3530 7.7220 14.87 6.66 4.98 38.81 3430 3535 8.0721 14.97 7.09 5.35 38.85 3430 3535 8.5122 15.09 7.49 5.71 38.90 3430 3540 8.9523 15.11 7.90 6.05 38.94 3435 3550 9.3624 14.72 8.34 6.44 38.98 3450 3565 9.7325 14.49 8.75 6.80 39.01 3470 3590 10.0126 14.27 9.15 7.15 39.04 3500 3620 10.20˙ M N / (cid:12) )12 11.82 3.39 2.22 38.40 3460 3685 4.9513 12.59 3.82 2.55 38.47 3450 3630 5.4414 13.35 4.24 2.88 38.54 3440 3595 5.9015 14.12 4.65 3.22 38.60 3440 3580 6.3116 14.87 5.07 3.58 38.65 3440 3570 6.7117 15.62 5.49 3.93 38.70 3450 3565 7.0418 16.37 5.91 4.30 38.74 3460 3565 7.3219 17.01 6.36 4.69 38.78 3465 3560 7.6520 17.66 6.79 5.08 38.82 3465 3560 8.0421 18.27 7.22 5.45 38.87 3460 3560 8.4922 18.69 7.67 5.85 38.92 3455 3560 9.0023 19.25 8.07 6.20 38.96 3455 3565 9.4224 19.58 8.52 6.59 39.00 3455 3565 9.8625 19.95 8.94 6.94 39.03 3455 3570 10.2326 20.14 9.38 7.31 39.06 3460 3575 10.6027 20.36 9.80 7.67 39.09 3465 3585 10.9228 20.66 10.19 7.99 39.11 3470 3590 11.1929 20.61 10.65 8.38 39.14 3480 3600 11.5030 20.72 11.08 8.73 39.17 3495 3610 11.7531 20.83 11.49 9.10 39.19 3510 3625 11.9732 20.85 11.95 9.48 39.21 3525 3640 12.1833 20.98 12.37 9.86 39.23 3540 3655 12.3734 21.03 12.85 10.27 39.26 3560 3675 12.5535 21.08 13.33 10.70 39.28 3580 3690 12.7436 21.36 13.69 11.01 39.29 3590 3705 12.8637 21.40 14.18 11.46 39.31 3610 3725 12.9938 21.63 14.56 11.80 39.33 3625 3735 13.1539 21.80 14.98 12.18 39.34 3635 3750 13.29˙ M N /
10 (not binned)12 11.88 3.21 2.09 38.36 3480 3725 4.6915 14.76 4.46 3.07 38.58 3465 3615 6.0520 19.48 6.66 4.96 38.81 3500 3590 7.7625 24.00 8.86 6.89 39.02 3495 3605 9.8830 27.75 11.39 8.95 39.17 3510 3615 11.7335 31.55 13.89 11.05 39.29 3535 3625 13.1740 35.25 16.23 13.19 39.37 3575 3660 14.2345 36.65 21.41 18.22 39.53 3610 3650 16.7350 41.40 21.21 17.87 39.52 3635 3710 16.2855 44.27 23.73 20.20 39.58 3655 3715 17.2660 47.31 25.64 22.00 39.62 3680 3740 17.79
Note . — All quantities are 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. Measures of “Explodability”
Early theoretical studies of supernovae noted a strongcorrelation of a rapidly declining density outside the ironcore with the degree of difficulty encountered in tryingto explode the star using the neutrino energy transport(e.g., Burrows & Lattimer 1987; Fryer 1999). O’Connor& Ott (2011) introduced a simple, single parameter mea-sure of this density decline called the “compactness pa-rameter”: ξ M = M/ M (cid:12) R ( M bary = M ) / (cid:12)(cid:12)(cid:12) t bounce , (3)where R ( M bary = M ) is the radius of the Lagrangianmass shell enclosing mass M in the presupernova star.The fiducial mass is often chosen as the innermost 2.5M (cid:12) so that for a wide range of initial masses it notonly encloses the iron-core but samples enough of theoverlying shell material around it. Though the param-eter 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 in-sensitive 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 withwhich model stars exploded and the ξ parameter, in thesense that stars with small ξ , i.e., steep density gradientsoutside the iron core were easier to explode using a stan-dard, albeit approximate, set of supernova physics. Al-though this is a useful starting point, a single parameterconveys limited information about the structure of thecore and more physics-based representations followed.In particular, Ertl et al. (2016) suggested an alternativetwo-parameter characterisation based upon M , the masscoordinate, in solar masses, where the entropy per baryonreaches 4 . k B , and the radial gradient, µ , of the densityat that point. In practice, µ , was obtained by evaluatingthe change in mass over the change in radius betweenseveral mass shells separated by 0.3 M (cid:12) in the vicinityof M , i.e., µ = dm/ M (cid:12) dr/ dm = 0 . (cid:12) . Smaller values of µ thus implysteeper density gradients (less change in enclosed masswhen the radius changes). The quantity M has longbeen used to locate the the steep density decline oftenassociated with a strong oxygen-burning shell in the pre-supernova star (e.g., Woosley & Heger 2007). Where theentropy per baryon abruptly rises at nearly constant tem-perature, the density declines. When the core collapses,the arrival at the neutrinosphere of lower density matterreduces the “ram pressure” and facilitates the launch ofan outgoing shock. Besides its location, it is helpful tocharacterise the rate of the density change, which is therole played by µ .It is reasonable that higher neutrino luminosities andsmaller accretion rates on the proto-neutron star shouldfavor explosion. Ertl et al. (2016) made the case that µ ,multiplied by its radius (i.e., the radius at which M isfound) and divided by a collapse time scale, is a surro-gate for the accretion rate at the time of explosion. Theproduct µ M is similarly a surrogate for the accretion luminosity. This assumes, as the models suggest, that theradius where the neutrino luminosity is generated doesnot vary greatly with mass. Explosion is thus favoredby small µ (accretion rate) and large µ M (luminos-ity). Ertl et al. (2016) determined a simple conditionexpressed as a straight line µ = k ( µ M ) + k (5)such that models that lay below the line (lower µ ) weremore likely to explode than those above. The values k and k were determined from several large scale simula-tions and leading models for SN 1987A. For one repre-sentative set, k = 0.194 (M (cid:12) ) − and k = 0.058 (“N20”model in their Table 2).More recently, M¨uller et al. (2016) adopted a differentapproach using a semi-analytic model for the explosionbased upon presupernova properties. By using a fullerdesription of the presupernova star than afforded by justone or two parameters, they were able to approximatelyestimate not only the success or failure of the explosion,but also the explosion energy, ejected nickel mass, andcompact remnant mass. This approach does not include,in its present form, several important pieces of explo-sions physics, such as proto-neutron star cooling, self-consistent fallback and explosive nuclear burning, andthus it is not as powerful as the parameterised simula-tions such as Sukhbold et al. (2016). Its main power liesis its ability to rapidly calculate general trends in explo-sion properties that depend on presupernova structure.Their approach uses five physically motivated parametersto determine the outcome. Here we analyze our new pre-supernova models using their standard values - β expl =4 , ζ = 0 . , α out = 0 . , α turb = 1 . , and τ . = 1 . Results for the Compactness Parameter
Fig. 8 shows the compactness parameter, ξ . , i.e., ξ forour three model sets measured at a fiducial mass of 2.5M (cid:12) at the time the presupernova collapse speed reaches900 km s − . One notable feature seen in Sukhbold &Woosley (2014) is the relatively large “scatter” of com-pactness parameter between 18 and 22 M (cid:12) , which arenow shifted to roughly 14 and 19 M (cid:12) (Panel a) withupdated neutrino-losses ( § ∼
20 M (cid:12) marks the transitionregion between exoergic convective central carbon burn-ing, at lower masses, to endoergic radiative burning athigher ones. That is, at high mass the energy producedby carbon fusion does not exceed neutrino losses in thestar’s center and it proceeds from helium depletion tooxygen ignition on a short Kelvin-Helmholtz time scale.As a result subsequent burning especially in shells, wherethe net energy generation is positive, proceeds differently.The carbon convective shells move inwards (with increas-ing initial mass) up until this transition point and areweaker and more in number (Barkat 1994).Quite noticeable in the new results and also in theprior work of M¨uller et al. (2016) are concentrations ofpoints in the compactness parameter plot, which werenot clearly seen in our earlier studies due to much coarserincrements in mass-space. Points are not randomly scat-tered between some local maximum and minimum value.Note for example, the existence of two, and possibly three4
13 14 15 16 17 18 19 M ZAMS [M (cid:12) ]0 . . . . . . . ξ . (a) ˙ M N
10 20 30 40 50 60 M ZAMS [M (cid:12) ]0 . . . . . . . ξ . (c) ˙ M N ˙ M N / M N / . . . . . . M He [M (cid:12) ]0 . . . . . . . ξ . (b) ˙ M N M He [M (cid:12) ]0 . . . . . . . ξ . (d) ˙ M N ˙ M N / M N / Fig. 8.— (a): The compactness parameter for a fiducial mass of 2.5 M (cid:12) evaluated at the presupernova stage shown for models below < (cid:12) as a function of initial mass for the ˙ M N set. Even where rapid variations occur, the scatter is not wholly random. Points tend toconcentrate in one of several solutions. (b): Same points as in (a) now shown as a function of helium core mass. The branching solutionsare more clearly seen. The range of ξ . seen in this mass range includes stars that might explode or implode, and thus both outcomes arepossible for stars with nearly identical initial mass. (c): Compactness shown for all three values of mass loss. Slight shifts between the setsare due to the growth of the helium core by hydrogen shell burning, resulting in a slightly larger compactness for the models with reducedmass loss. (d): same points as in (c) are now shown as a function of helium core mass. Much of the shift seen in (c) is gone since thehelium core mass is a chief determinant of the core structure. discrete solutions for helium core masses near 4.8M (cid:12) (Panel b). Multiple branches also exist at other massesbut with less clarity. For a helium core mass of 5.6 M (cid:12) ,the compactness parameter for stars with nearly identi-cal masses varies by a factor of three, from ξ . = 0 . <
20 M (cid:12) ), but they grow with increasingmass and for our most massive models the ˙ M N /
10 set isshifted by a few M (cid:12) with respect to ˙ M N / § (cid:12) (Sukhbold5& Woosley 2014). A shallow dip, occurs near 50 M (cid:12) followed by a slow rise at still higher masses (Sukhbold& Woosley 2014) until finally the pulsational-pair do-main is reached near 70 M (cid:12) . As discussed in Sukhbold& Woosley (2014) and Sukhbold (2016), this structureis crafted out of an overall gradually rising curve by theoperation of an advanced stage shell burning episode ofone fuel modulating the core burning episode of the nextfuel. In particular, the structure near the first peak isprimarily driven by the effect of shell carbon burning oncore oxygen burning, while the structure near the secondpeak is driven by the effect of shell oxygen burning oncore silicon burning. At very high mass, oxygen burningultimately ignites a strong shell outside of 2.5 M (cid:12) , andtherefore the compactness stays high.Fig. 9 compares the new results with those of Sukhbold& Woosley (2014) and M¨uller et al. (2016), and illustratesthese points. Most notable is a shift of the new modelsabove ∼
14 M (cid:12) downwards in mass by about 2 to 3M (cid:12) compared with Sukhbold & Woosley (2014). Thisis a consequence of fixing the bug in the neutrino lossesas discussed in §
2, and not due to differing resolution.Models below 14 M (cid:12) are affected very little, but themore massive stars change appreciably. The region ofvariability which was between 18 and 22 M (cid:12) in Sukhbold& Woosley (2014) is now shifted down to 14 to 19 M (cid:12) .If anything that variability, is now greater with the finerresolution, perhaps due to the larger number of modelssurveyed.Indeed, the comparison with M¨uller et al. (2016), whoused the corrected neutrino loss rate and sampled manymore models is excellent, both in the location of struc-tures and the range of ξ . within those structures. Thisagreement persists despite the use of 4 to 10 times moremass shells in each of the new models and a radical de-crease in the surface boundary pressure. The latter af-fected the mass lost by the star, but not appreciably thehelium core masses. There is no reason to believe thatstill finer zoning, smaller time steps, or a different reac-tion network will greatly alter these results, unless thecode physics itself (reaction rates, semiconvection, rota-tion, etc.) is changed. The studies of M¨uller et al. (2016)and the present work are mutually confirming. New Results and the Ertl Parametrization
An important subsidiary question is whether ξ . isreally the best measure for presupernova core structure.Might some of the variability seen in Fig. 8 be simplybecause of the choice of a single arbitrary point in the starto sample its structure? Perhaps other parametrizationsmight give less variability and more reliable predictions?In particular, ξ . is sensitive to recent shell activity inthe vicinity of 2.5 M (cid:12) that might not always describewell what went on deeper inside.The Ertl parametrization of our results is given inFig. 10. As previously discussed ( § µ M less than 0.25. In some ranges, the difference be-tween solutions is enough to significantly affect the prob-able outcome of the explosion. Stars with µ M lessthan 0.25 typically have masses less than 20 M (cid:12) . Formore massive stars, and, in particular, for µ M greater
12 14 16 18 20 22 24 26 M ZAMS [M (cid:12) ]0 . . . . . ξ . this workSW1415 20 25 30 35 40 M ZAMS [M (cid:12) ]0 . . . . . . . ξ . ˙ M N , high res . ˙ M N / , high res . M¨ueller + 2016 , low res . Fig. 9.—
Presupernova compactness compared with previous re-sults for solar metallicity stars obtained by Sukhbold & Woosley(2014) (top frame) and M¨uller et al. (2016) (bottom frame). Thenew models differ little from Sukhbold & Woosley (2014) below 14M (cid:12) , but the peak in compactness formerly at 23 M (cid:12) is shifteddownwards to about 20.5 M (cid:12) . The “noise” noted by Sukhbold& Woosley (2014) from 17 to 22 M (cid:12) is also shifted downwardsto 14 to 19 M (cid:12) and now with finely incremented models it showsseveral concentrations of points suggestive of a multi-valued solu-tion. The lower panel shows comparison with the work of M¨ulleret al. (2016). Since the ˙ M N and ˙ M N / M N set for lowermass and ˙ M N / than 0.25, explosion seems quite unlikely. Reasons forthe multi-valued solution are discussed in § (cid:12) may not6 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 .
30 0 . µ M . . . . . . . µ e x p l o s i o n i m p l o s i o n M Z A M S [ M (cid:12) ] Fig. 10.— “Explodability” according to sample “engine” modelfrom Ertl et al. (2016). The Ertl parameters, µ and µ M , areshown for the new models with standard mass loss rates. See textfor definitions. Points above the dashed line are more likely toexplode than those below. Multiple solutions are clearly visible for µ M below 0.25 (M (cid:12) ) − The scatter in this figure is a lot lessthan Fig. 9. be so much a consequence of “non-convergence” of themodels, but as a poor representation of the results. INTERPRETATION
Why is the core structure of presupernova stars non-monotonic and multi-valued for some ranges of mass?The short answer is that the advanced stage evolutionof massive stars involves two to three carbon burningshells (plus central core burning for models below ∼ (cid:12) ) 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 shellshave finite sizes and their number is an integer. Thetransitions are abrupt and small changes upstream inthe strength or extent of one or more shells can send astar down one path or another. This is especially true forstars below 19 M (cid:12) where there are more carbon burningshells.Globally, ξ . , µ , and M all increase gradually withmass. Higher mass stars have greater entropy in theirmiddles and are less degenerate in their final stages.Greater degeneracy leads to an increased central concen-tration of the mass and “core convergence.” For lighterstars, the presupernova structure resembles more a whitedwarf embedded in a low density envelope, where thedensity declines rapidly at the edge of the central “whitedwarf’. This effect is clearly at work in the lightest starssurveyed. Below 13 M (cid:12) , compactness defined at a massof 2.5 M (cid:12) has little meaning as a measure of explod-ability, since the fiducial point is not even inside of thehelium burning shell. For stars in this mass range, thehelium burning shell always lies outside 2 × cm, thusguaranteeing a small compactness parameter. For a 10M (cid:12) presupernova star, 2.5 M (cid:12) is even outside the heliumcore and in the hydrogen envelope. For all parametriza-tions of core structure considered in this paper, theselight stars will always explode and need no further dis-cussion here.Degeneracy remains an important considerationthough, even for the larger stars. A 12 M (cid:12) presupernova . . . . . . M [M (cid:12) ]0 . . . . . . . ξ . M Z A M S [ M (cid:12) ] . . . . . . M [M (cid:12) ]00 . . . ξ . . . . . . . . M H e [ M (cid:12) ] Fig. 11.—
Compactness enclosing innermost 2.5 M (cid:12) at the timeof presupernova for ˙ M N / M N /
10 set of models are plottedas a function of M , the Lagrangian coordinate of the first massshell where the entropy per baryon exceeds 4 . k B . For lower massmodels three major branches are clearly apparent below M ≈ (cid:12) . The lower panel shows a zoomed in version of the toppanel, characterizing models with M ZAMS <
20 M (cid:12) , and withcolor contours denoting their corresponding helium core masses. star has a core that is degenerate (degeneracy parameter η ≡ µ/kT >
0) out to 1.46 M (cid:12) . For 30 M (cid:12) , degener-acy extends out to 1.87 M (cid:12) . Even for 100 M (cid:12) , the ironcore is degenerate out to 2.1 M (cid:12) . Most of the increase indegeneracy occurs between carbon depletion and oxygendepletion in the core when the core is efficiently cooledby the neutrinos (Sukhbold & Woosley 2014). Degen-eracy affects the core structure making the (thermallyadjusted) Chandrasekhar mass relevant for all stars thatare likely to explode as common supernovae. On theother hand, the shells outside what will be the iron corein the presupernova star, and in particular the carbonand oxygen burning shells, are always non-degenerate.Strong burning thus leads to expansion that affects both ξ . and µ .The upper panel of Fig. 11 shows a generally strongcorrelation among the compactness parameter, ξ . , and M , except at very high and low masses. A correlationis expected since a sharp density decline close to the ironcore requires a large radius to enclose 2.5 M (cid:12) , thus small ξ . . The sharp density decline implies a deeper location7where the entropy per baryon exceeds 4 . k B and thus asmaller M and a steeper mass gradient at that point aswell. The reversal of the plot for the most massive modelsis a consequence of both the increasing central entropyand the outward migration of the first oxygen shell. Thismigration causes a non-monotonic dependence on massfor both M and ξ . as shown in Fig. 8, but M starts todecline at a slightly lower initial mass due to increasingentropy. Therefore when the location of M recedes from ∼ (cid:12) , the ξ stays roughly constant, and resultsin the reversal near ∼
35 M (cid:12) .The lower panel of Fig. 11 shows the behavior of com-pactness and M for stars lighter than about 20 M (cid:12) isclearly multi-valued. Three distinct branches are appar-ent for the lighter stars. There is no reason to believethat more dense grid of calculations would randomly fillout the spaces between the branches.Fig. 12 delves deeper into this correlation between com-pactness and M and helps to understand why both aremulti-valued relations of mass for moderate mass stars(13 < M ZAMS <
19 M (cid:12) ). The first panel shows M inthe presupernova star as a function of helium core mass.A cleaner pattern results from using the helium core massinstead of the main sequence mass since it eliminates thevariations due to envelope mass loss (Fig. 5). As the fig-ure 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 (exclud-ing neutrino losses) is a maximum. This is frequently,though not always interior to the boundary of the “siliconcore”, inside of which the silicon mass fraction is greaterthan the oxygen fraction. Though one might naively as-sume that the oxygen burning shell is at the edge of thesilicon core, this is not generally true. There are oftengradients on the silicon to oxygen ratio left behind byreceding convection or radiative burning. In these cases,the oxygen shell lies inside the silicon core. On the otherextreme, in some stars the oxygen shell can be so ac-tive as to merge with the carbon and neon burning shellsproducing one large convective shell where all three areburning vigorously and the silicon core and oxygen shellare coincident. When this occurs, the compactness isusually small (Fig. 16 of Sukhbold & Woosley 2014).Many of the red points in the top panel of Fig. 12 are ofthis sort.It is not surprising that the oxygen burning shell, thestrongest burning shell during the late stages of stellarevolution, is usually the location of a jump in entropy.The oxygen fusion rate is very temperature sensitive, sooxygen burns at a nearly constant temperature. Theoverlying star expands, decreasing the local density andthus increasing the entropy. The quantity µ in the Ertlparametrization is a measure of the strength of this burn-ing.The middle panel of Fig. 12 shows that the behaviorof ξ . in this same mass range (equivalent to Panel b ofFig. 8) correlates reasonably well with that of M , includ-ing, approximately, the location and extent of clusters ofsolutions. Within a cluster, there is usually a quasi-linearrelation with a well-defined slope. A larger value of M implies a larger value of ξ . . The farther out in the starthe strong burning shell, the less centrally concentratedis its density. This correlation is not perfect though.Sometimes long lines in the M plot, e.g. the long string . . . . . . M He [M (cid:12) ]1 . . . . . . . I n t e r i o r M a ss [ M (cid:12) ] M M Si M O − shell . . . . . . M He [M (cid:12) ]0 . . . . . . . ξ . . . . . . . M He [M (cid:12) ]0 . . . . . . . µ Fig. 12.— (Top:) M vs Helium core mass for presupernova starswith main sequence masses between 13 and 19 M (cid:12) . Light pinkpoints which usually underly green stars or red pentagon pointsindicate M . Green stars mark the location of the most activeoxygen burning zone and red pentagons are where the silicon massfraction first exceeds the oxygen mass fraction, going inward. Forclarity, red and green points not closely affiliated to a M within0.01 M (cid:12) have not been plotted. Usually the M point is foundat or near the maximum oxygen burning zone. Several clustersof points sometimes lying in nearly straight lines are apparent.(Middle:) The compactness parameter, ξ . plotted for the samehelium cores, equivalent to Panel b of Fig. 8. Note also a clusteringof points here and a high degree of correlation with M in the toppanel. (Bottom:) µ as a function of helium core mass. Note thealmost completely congruent pattern with ξ . , since both measurethe most active burning shell at the time of presupernova. M He = 4.1 M (cid:12) to 5.2 M (cid:12) , break into sev-eral segments with different ξ . . This is partly due tothe arbitrary pinning of ξ . to a single mass shell, butalso because specifying M alone does not measure thestrength of the burning there.Instead, as might be expected, the bottom panel ofFig. 12 shows that ξ . correlates better with µ , themass gradient at M . This correlation is very strongsince both quantities are sensitive to the strength of themost active burning shell in the last days of the star’s life.Why though are M and µ not randomly scattered be-tween their extrema? Why the patterns? The first panelof Fig. 12 suggests a reason. M traces the location ofthe strongest oxygen burning shell. Oxygen burns therebecause it is at the deepest location that has not alreadydepleted oxygen by prior shell burning. That location is,in turn, set by the extent of oxygen burning shell(s) dur-ing the previous evolution, which in turn depends uponthe entropy structure set up during carbon shell burning.Sukhbold & Woosley (2014) pointed out this corre-lation between the core compactness and the extent ofthe first oxygen burning convective shell (their Fig. 14).Fig. 13 shows a related quantity for the new model set.The figure shows the mass of the “silicon core”, the massinterior to which oxygen has burned out, at the time sil-icon burning ignites in the star’s center. Here a specificcentral temperature, 3 . × K, was chosen in orderto make the plot, but the conditions here reflect whatthe oxygen burning shell (or shells) have accomplishedprior to silicon ignition. Visual inspection of the con-vective history of these 700 models shows that the sil-icon core mass at silicon ignition is very nearly equalto the maximum extent of the second convective oxygenburning shell for stars below ∼ . (cid:12) (helium coremass 4.2 M (cid:12) ) and of the first oxygen convective shellfor more massive stars up to at least 20 M (cid:12) . In thisplot, we see the clearest evidence yet for regular, butnon-monotonic and occasionally multi-valued behavior.Two helium cores of very nearly 4.2 M (cid:12) can give riseto presupernova cores with structures in one of two well-defined states. Slight shifts in mass, composition, or evennumerical approach (including zoning) can send the starone way or another.Further analysis of the convective histories reveals thesystematics behind this behavior - at least in the 1Dcode, if not in nature: • The little jump at M He = 3.7 M (cid:12) ( M ZAMS = 13.2M (cid:12) ) from M Si = 1.7 M (cid:12) to 1.6 M (cid:12) , reflects theoperation of the third convective carbon burningshell. Below this mass, the third shell ignites in-side the former full extent of the second shell andoutside the effective Chandrasekhar mass. Thusthe core oxygen burning start to burn in a smallerextent while this outermost carbon shell operates.For masses above this value, the third shell ig-nites at the outer boundary of the previous shelland within the Chandrasekhar mass. Core oxygenburning now has to wait until this third shell iscomplete and therefore the base of this third shellsets (approximately) the extent of the first oxygenconvective shell, and thus the base of the secondone (Panel 13.45 M (cid:12) of Fig. 14). • The much larger decrease in M Si from 1.7 M (cid:12) to 1.4M (cid:12) at M He = 4.2 M (cid:12) ( M ZAMS = 14.6 M (cid:12) ) is dueto the diminished significance of the second oxygenburning shell. The silicon core shrinks to the extentof the first oxygen burning shell though the secondshell continues to be sporadically important for atime (Panel 15.01 M (cid:12) of Fig. 14). • The jump and wild variations starting at M He = 5.2M (cid:12) and extending up to 5.7 M (cid:12) ( M ZAMS = 17.2M (cid:12) to 19.0 M (cid:12) ) mark the transition from convec-tive carbon burning to radiative burning. After thetransition near ∼
19 M (cid:12) , carbon no longer burnsexoergically at the center of the star, and there areonly two convective shells before oxygen burningrather than three or more. During this transition,carbon core and shell burning vary greatly in lo-cation and extent. The large spread in M Si , andultimately in compactness, reflects the irregularityof this transition (Panels 17.90 and 18.81 M (cid:12) ofFig. 14).Other shell interactions result in the weaker multiplebranches seen in higher mass stars in Fig. 8. Despite thespecific masses given in the list above, none of these tran-sitions are abrupt and, given slight nudges, the star mayoscillate from one to the other solution when it is closeto boundaries. This leads to the “multivalued” behavior.Although complicated and probably sensitive to the one-dimensional treatment of the problem, the conclusion isthat the presupernova structure in stars from 13 to 19M (cid:12) results from an interplay of convective carbon andoxygen burning shells after carbon ignition.Putting these various factors together, a determinis-tic picture emerges. Multiple carbon burning stages -core burning below 19 M (cid:12) , and two or more episodes ofshell burning - act to sculpt an entropy distribution inthe core so that oxygen shell burning, when it occurs, ig-nites at and extends to various mass shells. Whether thenumber of carbon shells is two or three or four and thenumber of oxygen shells one or two strongly affects thecore structure of a presupernova star in the mass rangewhere a significant number of events are observed. Thesechanges in the extent of convective shells, which amplifysmall differences in the earlier evolution, may end up de-termining whether the star explodes or makes a blackhole. COMPACT REMNANTS AND EXPLOSIONS
Since the new models differ from those of Sukhbold& Woosley (2014) and Sukhbold et al. (2016), and to alesser extent, from M¨uller et al. (2016), it is worth revis-iting some of the conclusions of those papers regardingcompact remnants using the new models.Fig. 15 shows the models now expected to explodebased on the criteria of Ertl et al. (2016) and M¨ulleret al. (2016). These figures can be compared with Fig. 4of Ertl et al. and Fig. 6 of M¨uller et al., which theyvery closely resemble. The “N20” engine parameteri-zation was used to make the comparison with Ertl etal. (2016) and the “standard” choices of five parameterswere used for M¨uller et al. (2016) method.The most significant difference with the Ertl criterionis the shift of the pattern for models above about 149 . . . . . . M He [M (cid:12) ]1 . . . . . . . . M S i ( a t S ii g n . ) [ M (cid:12) ] Fig. 13.—
The size of the oxygen depleted core (oxygen less thansilicon or iron) as a function of helium core mass at silicon ignition(when the central temperature first reaches 3 . × K). This isapproximately the maximum extent of the first oxygen convectiveshell. The function has multiple trajectories that correlate withthe major features in the plot of M vs helium core mass in thetop panel of Fig. 12. M (cid:12) to slightly ( ∼
10 %) lower masses. Changes in aver-age quantities like the remnant masses are small, how-ever. Fig. 16 shows a comparison of IMF-weighted rem-nant mass distributions of imploding models assuminga Salpeter (1955) initial mass function with a power of-2.35. The usual two cases are considered: a) collapse ofthe 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 toa binary companion or, for very massive stars, to winds.The envelope might also be lost to a very weak explosionthat did not unbind the helium core (Quataert & Shiode2012; Fuller 2017; Lovegrove & Woosley 2013). Case b)is for more robust explosions. In each case there is anelement of uncertainty because the final mass dependson the mass loss rate, e.g., for red supergiants in Case b)and Wolf-Rayet stars in Case a). In making this Fig. 16,but not in computing the averages, only the mass rangecovered by our ˙ M N mass loss survey, 12 - 27 M (cid:12) , wasincluded.Astronomers, of course, observe black holes comingfrom all masses of stars, not just 12 - 27 M (cid:12) . To com-pute the averages for solar metallicity stars (and the an-swer will be sensitive to metallicity), the new results weresupplemented with the prior models of Sukhbold et al.(2016) for stars more massive than 27 M (cid:12) . For just thelimited mass range 12 - 27 M (cid:12) , the average black holemasses from Sukhbold et al. (2016) were previously 6.52and 15.3 M (cid:12) respectively for the helium core and fullstar assumptions. The corresponding new numbers are6.60 and 14.2 M (cid:12) . Considering the entire mass range ofstars that experience iron core collapse, and using theN20 parameters in Sukhbold et al. (2016), the previousaverages were 9.25 M (cid:12) and 13.7 M (cid:12) . The equivaluentnew numbers are 8.61 M (cid:12) and 13.5 M (cid:12) . A less than 1%adjustment has been applied to Table 4 of Sukhbold etal. (2016) based upon a slightly different way of interpo-lating 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 neutronstar is usually equal to M . After appropriately cor-recting for neutrino mass loss (same as in M¨uller et al.2016), the new average neutron star gravitational massin the range 12 - 27 M (cid:12) is 1.45 M (cid:12) . Explosions below12 M (cid:12) were calculated by Sukhbold et al. (2016) usingpresupernova models in which the neutrino bug had beenfixed. They can thus be combined with the current set.Assuming neutron star production above 27 M (cid:12) is neg-ligibly small, the new global average neutron star massis 1.38 M (cid:12) . In Table 4 of Sukhbold et al. (2016) thecorresponding global average was 1.41 M (cid:12) .The agreement of Fig. 15 with the earlier work (Fig.6 of M¨uller et al. 2016) is excellent, and it is thus ex-pected that the explosion outcomes will also be very sim-ilar. Indeed, plots of neutron star mass (Fig. 17), blackhole mass, and explosion energy as a function of initialmass (not shown) are virtually indistinguishable from thepanels in their Figure 2. This suggests that the differ-ence in envelope structure and zoning between the seriesof M¨uller et al. (2016) and the present paper did nothave much influence on the core structure and statisticalexplosion properties.Fig. 15 shows that a significantly larger fraction of starsin the interesting mass range 15 to 20 M (cid:12) explode usingthe M¨uller et al. (2016) formalism, even though ”N20” isone of the more energetic formulations of the Ertl et al.(2016) model. Conservatively, this could be regarded asan uncertainty in outcome until more realisic simulationsof the actual explosion can be done. The various ener-gies in the Ertl et al. (2016) model came, however, fromcalibrating a central engine with 1D neutrino transport,a shrinking protoneutron star, and fallback to SN 1987Ausing various presupernova models and might, for now,be considered the more realistic of the two. This cali-bration to SN 1987A has its own uncertainties though,especially since the structure of the presupernova starcould be different in more realistic binary merger models(e.g., Menon & Heger 2017, and references therein).Fig. 17 gives the expected neutron star masses thatresult when our new presupernova models are analyzedusing the formalism and standard parameters of M¨ulleret al. (2016). Fallback is neglected in their analysis,so the baryonic mass of the proto-neutron star is equalto the mass coordinate where the neutrino-driven en-gine shuts off a few hundred ms after shock revival. Asmentioned before, this correlates strongly with M (toppanel). There are branches of neutron star masses abovethe M line, which primarily originate from less mas-sive progenitors ( M ZAMS <
15 M (cid:12) ), reflecting the op-eration of a prominent second oxygen burning shell andthe resulting shallow entropy profile. The most massiveneutron stars have a baryonic mass that is bounded bythe base of the convective carbon burning shell, which islocated much further out than M in the presupernovastar. The second panel of Fig. 17 shows the expectedneutron star (gravitational) mass as a function of mainsequence mass. These results agree quite well with thoseof M¨uller et al. (2016, see their Fig. 2c) including the ex-istence of very massive neutron stars produced for mainsequence stars 14 - 15 M (cid:12) and the highly variable natureof the solution between 14 and 19 M (cid:12) . For the majority0 Fig. 14.—
Convection histories of four sample models from the range of initial masses from 14 to 19 M (cid:12) , representing key structuralchanges that are responsible for the significant variations in the final presupernova core properties. Each panel shows the evolution of theinnermost 4.5 M (cid:12) material roughly during the last thousand years of its life, i.e., from core carbon burning until presupernova. Coloredshades denote energy generation (red) and energy loss (blue) gradients. Hatched black regions are convective episodes. The x -axis is shownas the log of time until core collapse, t cc . The initial mass of each model is denoted on the top left corner of panels. Two inserted mini-plotsillustrate the compactness parameter (purple) and the lagrangian location of entropy per baryon equal 4 . k B point (green) respectively,both as a function of helium core mass corresponding to the above mentioned initial mass range. The crosses inside each mini-plots denotethe values corresponding to the model. (13.45 M (cid:12) :) The lowest mass stars have 4 or more convective carbon burning episodes followedby 3 oxygen burning episodes (including central core burning episode). The shell helium burning lies within the 2.5 M (cid:12) location, thepoint on which the compactness parameter, ξ . , is measured, and thus these models have small compactness and small M . (15.01 M (cid:12) :)Central carbon burning weakens with increasing mass, and as a result the shell carbon burning episodes gradually “migrate” inwards (withincreasing initial mass). Once the third shell carbon burning ignites within the effective Chandrasekhar mass, the core oxygen burning hasto wait until the overlying carbon burning episode is finished, which leads to significantly weakened or, in some cases absent, immediatesecond oxygen shell burning. In these models oxygen eventually burns later when silicon is already ignited in the core. This marks thetransition of M to the “second” branch. Due to its sensitivity to the last major shell burning episode, the compactness parameter hasmultiple solutions depending on the final configuration of oxygen, neon and carbon burning (i.e. together burn vigorously in one shell,or separately). (17.90 M (cid:12) :) Around ∼
18 M (cid:12) it becomes increasingly harder for carbon to burn convectively near the center. The tinyconvective central burning episode is followed by a long lasting radiative flame, and as a result, with increasing initial mass, the base of theconvective carbon burning shells start to migrate outwards. (18.81 M (cid:12) :) Eventually at high enough initial mass all central carbon burningis radiative. Models between 17 and 19 M (cid:12) have rapidly changing core structures, and thus rapidly varying ξ . and M . Same formatplots for all 800 models between 12.00 and 19.99 M (cid:12) are available online (see Footnote 9). of models, the resulting neutron star mass is tightly cor-related with M (Fig. 12) and oxygen depleted core atthe time of silicon ignition (Fig. 13), and in certain initialmass ranges the explosion of models with nearly identi-cal mass (or slightly different input physics) can resultin very different neutron stars.It is not expected that the nucleosynthesis and lightcurves calculated by Sukhbold et al. (2016) will be signif-icantly altered by using the new models, though further exploration is certainly encouraged. In particular, thedeficiency of light s -process elements seen in Sukhbold etal. (2016) will persist in the new model set, since mostof the production of these elements is due to the mostmassive stars (Brown & Woosley 2013) which still fail toexplode. The fraction of solar metallicity stars above 9M (cid:12) exploding as supernovae was 74% based on the N20parameterization (Table 4 of Sukhbold et al. 2016) andis now 65 %. Above 18 M (cid:12) , the fraction of SN reduces1
12 14 16 18 20 22 24 260 . . . . . ξ . Ertl + 2016 criterion explosionimplosion12 14 16 18 20 22 24 26 M ZAMS [M (cid:12) ]0 . . . . . ξ . M¨ueller + 2016 criterion explosionimplosion
Fig. 15.—
The compactness parameter shown in Fig. 8 is plottedagain twice, color coded as to the success or failure of the explosionbased upon the parameterization of Ertl et al. (2016) (top) andthe semi-analytical method of M¨uller et al. (2016) (bottom). Redsymbols denote successfull explosions by these criteria, and blacksymbols, failures. The very good agreement between these twodifferent approches suggests that both are good representations ofcore structure, though the criterion of M¨uller et al. (2016) favorsslightly more explosions for the parameters chosen. to only 8%, which actually is in a slightly better accordwith the results from Smartt (2009, 2015), as comparedto the previous study. Since the envelope masses havenot changed appreciably and the explosion energies arealso expected to be unchanged, the light curves will beunaltered. The fractions given in their Table 4 for su-pernovae above 12, 20, and 30 M (cid:12) will also probably notchange within the variations already seen for the differentcentral engine characteristics. CONCLUSIONS
The full evolution of over 4,000 massive stars of so-lar metallicity in the mass range 12 to 60 M (cid:12) has beenstudied using unprecedented resolution, both in zonesper star and number of stars within a given mass range.The mass loss rate was varied, and an important bugin the neutrino loss routine was repaired that caused asignificant variation in the presupernova properties fromthose calculated by Sukhbold & Woosley (2014) and priorworks. Our chief conclusions are: (cid:12) ] s c a l e d P D F Su k hb o l d e t a l . ( ) T h i s W o r k He − core preSN . starimplosion implosion Fig. 16.—
The IMF-weighted black hole mass distributions forthe new ˙ M N set is shown in comparison with the results fromSukhbold et al. (2016) over the same initial mass range between 12and 27 M (cid:12) . The differences are small and well within the rangeof variation seen for different explosion models. The average blackhole masses for the older study is 15.3 and 6.52 M (cid:12) for two sce-narios respectively, while the new models yield 14.2 and 6.60 M (cid:12) . • The pattern of core compactness seen in previousstudies (e.g., Sukhbold & Woosley 2014; M¨uller etal. 2016) is robust. The “noise” in these studieswas not a consequence of inadequate zoning, butreflects real variability. The range of variation andthe location of peaks (Fig. 8) in the new studyare virtually identical to that seen by M¨uller etal. (2016), even though the new models use 4 to 10times the zoning and a much smaller surface bound-ary pressure. The results are also qualitatively sim-ilar to Sukhbold & Woosley (2014), but importantpeaks in the compactness plot are shifted down-wards by about 10 % in initial mass (Fig. 9) dueto the corrected neutrino loss rate. For the stellarphysics used and spherically symmetric nature ofthe calculation, the variations and peaks seen inFig. 8 are now well determined. • The large variation in core structure seen for starsbetween 14 and 19 M (cid:12) by Sukhbold & Woosley(2014) is not completely random. For a larger setof models with finer spacing in initial mass, severalbranches of solutions emerge. This behavior wasalso seen by M¨uller et al. (2016). The branchesapparently result from variations in the locationof the oxygen burning shells in the presupernovastar, which in turn “remember” the location of sev-eral carbon burning shells of variable extent. Somenoise is introduced by the fact that the carbonabundance and size of the carbon-oxygen core arenot precisely monotonic with mass (Fig. 7). This,in turn, reflects the operation of semiconvection atthe boundary of the hydrogen and helium convec-tive core. In the mass ranges 14 to 19 M (cid:12) and 22to 24 M (cid:12) , the presupernova core structure is moresensitive to events in the last year of a star’s life(and sometimes the last hours) than to the star’sinitial mass. Presupernova structure in these massranges does not necessarily follow the Vogt-Russell2 . . . . . . M [M (cid:12) ]1 . . . . . b a r y o n i c N S m a ss [ M (cid:12) ]
12 14 16 18 20 22 24 26 M ZAMS [M (cid:12) ]1 . . . . . . . . g r a v i t a t i o n a l N S m a ss [ M (cid:12) ] Fig. 17.—
Neutron star mass distribution resulting from analyz-ing the new model set with normal mass loss ( ˙ M N ) for stars inthe mass range 12 to 27 M (cid:12) using the approach of M¨uller et al.(2016). (Top:) The baryonic mass of the expected neutron staris plotted as a function of M , the mass where the entropy perbaryon reaches 4 . k B in the presupernova star ( § M . Many points are found here because M is usually the“mass cut” in a successful explosion, unless the final entropy pro-file is significantly more shallow. (Bottom:) The gravitational mass(adjusted from the baryonic mass in the same way as in M¨ulleret al. 2016) of the expected neutron star as a function of mainsequence mass. More massive neutron stars are typically madeby more massive main sequence stars with a greater ξ . and µ ,but the most massive points on this plot come from models with M ZAMS ∼ −
15 M (cid:12) , the most massive models with a significantsecond oxygen shell burning. The existence of multiple branchesfor these quantities (Fig. 12 and Fig. 13) results in a large rangeof neutron star masses being acessible by stars of nearly the samemain sequence mass for some initial masses. Note also the clusterof neutron stars with mass near 1.65 M (cid:12) coming from the mostmassive supernova progenitors. theorem (Vogt 1926; Russell 1927). It may be morecomparable to weather on earth. • Even accounting for these systematics, thebranches of solutions are still noisy. The level ofnoise is reduced if one characterizes the presuper-nova star by its helium core mass, not its startingmass or final total mass. The scatter is also reducedin the Ertl two-parameter characterization of core-structure rather than the O’Connor-Ott compact- ness parameter. • Given the large number of models, it is possible togive statistically meaningful results for the radius,luminosity and effective temperature of supernovaprogenitors ( § • The mass distributions of neutron stars and blackholes resulting from supernovae in the mass rangestudied are not greatly altered from our earlier sur-veys ( § (cid:12) . The average black hole mass, if onlythe helium core implodes, is 8.61 M (cid:12) . If the en-tire presueprnova star collapses, the average blackhole mass is 13.5 M (cid:12) . The fraction of stars above9 M (cid:12) that explode rather than collapsing is esti-mated to be 65%. The fraction above 18 M (cid:12) is8%. Nucleosynthesis and light curves are basicallyunchanged.Binary mass exchange and rotation undoubtedly playkey roles, but have been neglected in this study. Tofirst order, stars that end up with the same final he-lium core mass will have similar presupernova compact-ness and fates. Similar systematic variations and multi-valued solutions are expected to persist. The statisticalaverages of compact remnant masses may vary, however,and certainly the light curves will be different. Our mod-els are publically available to those wanting to estimateoutcomes based upon their own distribution of heliumcore masses, carbon-oxygen core masses, etc. In the fu-ture, we will consider rotating models, but this was prin-cipally a study of how resolution affects the solution toa well-defined, frequently studied problem.Our results suggest a slightly different strategy to thestudy of presupernova evolution and supernova model-ing than sometimes used in the past. Given the varia-tion in outcome for stars of nearly the same mass, or thesame mass with different codes, a statistically meaning-ful sample of models must be calculated before drawingstrong conclusions about supernova mass ranges, rem-nant mass distributions, nucleosynthesis, etc. Histori-cally, researchers have sometimes focused on the calcu-lation of just a few masses, e.g., 15, 20, 25 M (cid:12) , andsought to test the senitivity to changes in physics in justthose cases. Here we see that calculating a statisticallymeaningful sample may be just as important as gettingthe physics precisely right. The size of such a sampledepends upon the need to resolve regions of rapid vari-abilility found with a given code and physics, but a min-imum initial mass grid of 0.1 M (cid:12) is reasonable withinsuch regions.Some of the models calculated here had merged car-bon and and oxygen convective shells at the end. Manydid not. Others were only separated by a single thinzone from being coupled. Similar to the numerical arti-fact seen when helium burning is calculated with too lit-tle semiconvection, these mathematical bifurcations maynot be real and might be overcome, or at least smoothedout (Alexakis et al. 2004), by increasing the convectiveovershoot or doing a multi-dimensional calculation (e.g.,Meakin & Arnett 2007; Viallet et al. 2013; Arnett &Meakin 2016; Jones et al. 2017; Cristini et al. 2017).3In general, linked convective shells give a more compactcore structure and favor explosions (see, e.g., Collins etal. 2017). Further study with other representations ofconvective overshoot mixing and multi-dimensional codesare thus encouraged.Semiconvection and overshoot mixing remain majoruncertainties in studies of this sort and also play impor-tant roles in determining the helium and carbon-oxygencore masses. How KEPLER treats semiconvection is de-scribed in Weaver & Woosley (1993).It is expected thatother studies using different treatments will find resultsthat differ in important detail from those presented here.The overall pattern, multi-peaked structure and range ofvariability of the compactness parameter and other mea-sures of core structure should persist, however. Clustersof solutions due to variable numbers and extents of thecarbon and oxygen shells should also be a common fea-ture. Further exploration is again encouraged.All of the new presupernova models presented here areavailable online. Also online are plots of the convectivehistories of 800 models between 12 and 20 M (cid:12) of the ˙ M set (see Footnote 9). Other auxiliary prespernova models(e.g., those presented in § § ACKNOWLEDGMENTS
We thank Raphael Hirschi for a careful review of thiswork and for providing many useful suggestions. We alsothank Bernhard M¨uller for helping us to process our newprogenitor results through his semi-analytical explosionmodelling and feedback on this manuscript. ThomasJanka and Thomas Ertl provided valuable insight intothe supernova explosion models and their dependence onprogenitor properties. We thank Todd Thompson andJohn Beacom for useful comments on the manuscript.We also appreciate the efforts of Gang Gao in discover-ing the neutrino bug in the
KEPLER code. All numeri-cal
KEPLER calculations presented in this work were per-formed on the
RUBY cluster at the Ohio SupercomputerCenter (Ohio Supercomputer Center 1987). TS is partlysupported by NSF PHY-1404311 to John Beacom. SWis supported by NASA NNX14AH34G. AH is supportedby an ARC Future Fellowship, FT120100363, and, inpart, by the National Science Foundation under GrantNo. PHY-1430152 (JINA Center for the Evolution ofthe Elements).4
REFERENCESAlexakis, A., Calder, A. C., Heger, A., Brown, et al. 2004, ApJ,602, 931Appenzeller, I. 1970, A&A, 5, 355Arnett, W. D., & Meakin, C. 2016, Reports on Progress inPhysics, 79, 102901Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009,ARA&A, 47, 481Barkat, Z. 1994, Supernovae, 31Beasor, E. R., & Davies, B. 2018, MNRAS, 475, 55Brown, J. M., & Woosley, S. E. 2013, ApJ, 769, 99Buchmann, L. 1996, ApJ, 468, L127Burrows, A., & Lattimer, J. M. 1987, ApJ, 318, L63Chieffi, A., & Straniero, O. 1989, ApJS, 71, 47Collins, C., M¨uller, B., Heger, A. 2017, MNRAS, submitted;arXiv/170900236Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471,279Davies, B., Kudritzki, R.-P., Plez, B., et al. 2013, ApJ, 767, 3Ertl, T., Janka, H.-T., Woosley, S. E., Sukhbold, T., & Ugliano,M. 2016, ApJ, 818, 124Farmer, R., Fields, C. E., Petermann, I., et al. 2016, ApJS, 227,22Fryer, C. L. 1999, ApJ, 522, 413Fuller, J. 2017, MNRAS, 470, 1642Heger, A. 1998, Ph.D. Thesis,Heger, A., Langer, N., Woosley, S.E. 2000, ApJ, 52, 368Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS,102, 411Jones, S., Andrassy, R., Sandalski, S., et al. 2017, MNRAS, 465,2991Lauterborn, D., Refsdal, S., & Weigert, A. 1971, A&A, 10, 97Levesque, E. M., Massey, P., Olsen, K. A. G., et al. 2005, ApJ,628, 973Lovegrove, E., & Woosley, S. E. 2013, ApJ, 769, 109Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448Meibom, A., Krot, A. N., Robert, F., et al. 2007, ApJ, 656, L33Menon, A., & Heger, A. 2017, MNRAS, 4649-4664, 469M¨uller, B., Heger, A., Liptai, D., & Cameron, J. B. 2016,MNRAS, 460, 742Nieuwenhuijzen, H., & de Jager, C. 1990, A&A, 231, 134Ohio Supercomputer Center 1987, Ohio Supercomputer Center,http://osc.edu/ark:/19495/f5s1ph73O’Connor, E., & Ott, C. D. 2011, ApJ, 730, 70 O’Connor, E., & Ott, C. D. 2013, ApJ, 762, 126Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15Pejcha, O., & Thompson, T. A. 2015, ApJ, 801, 90Quataert, E., & Shiode, J. 2012, MNRAS, 423, L92Rauscher, T., Heger, A., Hoffman, R. D., & Woosley, S. E. 2002,ApJ, 576, 323Renzo, M., Ott, C. D., Shore, S. N., & de Mink, S. E. 2017, A&A,603, A118Russell, H. N, 1917, in