Physical Correlations and Predictions Emerging from Modern Core-collapse Supernova Theory

In this paper, we derive correlations between core-collapse supernova observables and progenitor core structures that emerge from our suite of 20 state-of-the-art 3D core-collapse supernova simulations carried to late times. This is the largest such collection of 3D supernova models ever generated and allows one to witness and derive testable patterns that might otherwise be obscured when studying one or a few models in isolation. From this panoramic perspective, we have discovered correlations between explosion energy, neutron star gravitational birth masses, 56Ni and α-rich freezeout yields, and pulsar kicks and theoretically important correlations with the compactness parameter of progenitor structure. We find a correlation between explosion energy and progenitor mantle binding energy, suggesting that such explosions are self-regulating. We also find a testable correlation between explosion energy and measures of explosion asymmetry, such as the ejecta energy and mass dipoles. While the correlations between two observables are roughly independent of the progenitor zero-age main-sequence (ZAMS) mass, the many correlations we derive with compactness cannot unambiguously be tied to a particular progenitor ZAMS mass. This relationship depends on the compactness/ZAMS mass mapping associated with the massive star progenitor models employed. Therefore, our derived correlations between compactness and observables may be more robust than with ZAMS mass but can nevertheless be used in the future once massive star modeling has converged.


INTRODUCTION
Building on the early pioneering work of Colgate & White (1966), Arnett (1967), Wilson (1971), Bowers & Wilson (1982), and Bethe & Wilson (1985), with the recognition of the importance of neutrino-driven turbulence and convection (Herant et al. 1994;Burrows et al. 1995;Janka & Mueller 1996), and paralleling the developments in neutrino physics and supercomputers over the decades (Burrows 2019), core-collapse supernova (CCSN) theory has emerged after ∼60 years of progress to become a mature, if complicated and multi-physics, international enterprise.Neutrino heating, aided by the stress of neutrino-driven turbulence, is commonly accepted as the key mechanism of explosion.However, turbulence is chaotic and this implies that a range of outcomes and observable values are likely, even for the same massive star progenitor.The explosive expansion of the initially stalled bounce shock wave seems to be a critical phenomenon (Burrows & Goshy 1993) and is associated with the shock's pulsational instability, but a robust analytic model has proven elusive.This has necessitated detailed multi-dimensional numerical simulations, whose complexity and expense had slowed conceptual and quantitative progress.The accretion before explosion of the density jump associated with the Si/O interface in the massive star progenitor, itself correlated with "compactness" (see §2.1), oftimes inaugurates explosion (Vartanyan et al. 2018;Wang et al. 2022).However, this is not always the case, and its invocation as an exclusive trigger for explosion should be nuanced.Curiously, a quick explosion does not in general lead to an energetic explosion (see §3.1).
Rotation, though likely always present at some level, seems critically important in only a small subset of core collapses (hypernovae and Type Ic broad-line CCSNe), and its role in general may be subdominant, but has yet to be determined1 .Moreover, the role and evolution of magnetic fields, likely tightly coupled with rotation, is a subject of current intense research (Müller & Varma 2020;Varma & Müller 2021;Obergaulinger & Aloy 2020;Kuroda et al. 2020;Obergaulinger & Aloy 2021), but is thought to be only rarely of central importance.
Interestingly, in recent late-time 3D simulations, Wang & Burrows (2023a) observe aspherical winds (Duncan et al. 1986;Burrows 1987a;Burrows et al. 1995;Qian & Woosley 1996) emerging from the PNS after a delay of seconds following the onset of explosion; these produce secondary shock waves upon encountering the primary ejecta.As much as ∼10% of the explosion energy can be due to these winds.However, other groups have yet to simulate to late enough times after bounce to verify the details of this phenomenon, itself also associated with the important nucleosynthetic phase of α-rich freeze-out (Woosley et al. 2002;Nomoto et al. 2013;Arcones & Thielemann 2023;Wang & Burrows 2023b).Furthermore, which stars leave neutron stars and which black holes is still not resolved, though it is known that previous mappings were too simplistic (Burrows 1987b;Heger et al. 2003;Sukhbold et al. 2016;Ertl et al. 2016).There are new ideas on this partitioning (Burrows et al. 2023a), but these have not been thoroughly vetted and depend a great deal on the still-evolving progenitor models inherited.
Despite these many remaining issues and the physical and numerical complexity of the phenomenon, and despite the fact that chaotic turbulence renders a clear one-to-one mapping of progenitor to outcome naturally blurred, there has emerged a broad consensus on the viability of the turbulence-aided neutrino-driven mechanism of core-collapse supernovae (CCSNe).Modern sophisticated 3D simulations witness neutrino-driven explosions naturally and without artifice (Lentz et al. 2015;Roberts et al. 2016;Takiwaki et al. 2016;Müller et al. 2017;Ott et al. 2018;Glas et al. 2019;Müller et al. 2019;Burrows et al. 2019;Burrows et al. 2020;Müller & Varma 2020;Powell & Müller 2020;Stockinger et al. 2020;Kuroda et al. 2020;Bollig et al. 2021;Vartanyan et al. 2022Vartanyan et al. , 2023) ) and numerous reviews have articulated the mechanism's details (Janka 2012;Burrows 2013;Müller et al. 2016aMüller et al. , 2017;;Burrows et al. 2018Burrows et al. , 2020;;Mandel & Müller 2020;Burrows & Vartanyan 2021;Wang et al. 2022).With this paper, we summarize the physical and observable correlations that have emerged from the panoramic view enabled by our new and unprecedented set of twenty state-of-the-art 3D core-collapse simulations taken to late times.
In §2, we discuss progenitor issues and their structures to provide context.In §3, the basic explosion phenomenology is provided.This is followed by a graphical summary of the variety of immediate post-explosion morphologies and their general characteristics.Then, in §5, we provide for the first time using such a comprehensive long-term suite of 3D models the physical and observable correlations that define the explosion systematics and byproducts.These include the neutron star birth masses, explosion energies, 56 Ni yields, and recoil kicks.Finally, in §6, we summarize our salient results and provide a set of caveats to consider going forward.To create the set of twenty sophisticated 3D late-time simulations that we use in this correlation study, we employed the CPU machines TACC/Frontera (Stanzione et al. 2020) and ALCF/Theta and the GPU machines ALCF/Polaris and NERSC/Perlmutter.The specific setups and run parameters are described in Burrows et al. (2023a) and Vartanyan et al. (2023) and in the Appendix.

INITIAL PROGENITOR STRUCTURES
For this study we employ solar-metallicity, initially non-rotating, spherical massive star progenitor models from Sukhbold et al. (2016) and Sukhbold et al. (2018) and we do not include magnetic fields.The Fornax code and calculational details are described in the Appendix and in Skinner et al. (2019).As discussed in Wang et al. (2022), it is predominantly the density profile at collapse that determines explodability and this profile is not rigorously monotonic with zero-age main-sequence (ZAMS) mass.Figure 1 depicts these mass density profiles for the models we include in this paper.As one sees, while the steepness of the envelope in the Chandrasekhar-like core is roughly an increasing function of ZAMS mass, it is not rigorously so, nor is the "compactness" parameter (O'Connor & Ott 2011) (see §2.1); we will find in this paper a tighter relation of the observables with compactness than with ZAMS mass ( §5).
Importantly, there remain many ambiguities in the initial progenitor models and, as important as they are to corecollapse supernova theory, the progenitor models can not yet be considered definitive.Issues such as the 12 C(α,γ) rate, binarity, overshoot, wind mass loss, shell mixing, semi-convection, magnetic field effects, and rotation are still to be retired.Moreover, a new generation of 3D progenitor models (Couch et al. 2015;Jones et al. 2016;Chatzopoulos et al. 2016;Müller et al. 2016b;Jones et al. 2017;Yoshida et al. 2019;Fields & Couch 2020;Bollig et al. 2021;McNeill & Müller 2021;Fields & Couch 2021;Varma & Müller 2021;Yoshida et al. 2021) is emerging.As a result of these many remaining questions and uncertainties, models for the stellar evolution of massive stars to the point of collapse can not be considered to have converged.Nevertheless, we use the Sukhbold et al. (2016) and Sukhbold et al. (2018) model suite as the best and most comprehensive available.Moreover, we suggest that our derived correlations between observables and compactness may be usefully robust, providing ongoing insight as the progenitor model studies evolve into the future.
Figure 2 shows the total stellar mass (left) that remains at collapse and the corresponding C/O and He core masses (right) as a function of ZAMS mass for the models we have chosen for this comprehensive study.We note that despite the unprecedented number of models we have included in this study, there are many more models from the Sukhbold et al. (2016) and Sukhbold et al. (2018) left unsimulated.Nevertheless, our goal was to well sample the collection, if not optimally.

Compactness
The shallowness of the mass density profiles depicted in Figure 1 is related to the post-bounce mass accretion rate history and the accretion component of the neutrino luminosities.It is through these quantities that the complicated dynamics of the supernova mainly rests.One crude, but useful, metric of this shallowness is the compactness.Compactness (O'Connor & Ott 2011) 2 is not a predictor of explodability (Wang et al. 2022) in our simulations.Nevertheless, it is a useful, if gross, measure of progenitor structure and correlates well with numerous observables and physical characteristics of supernova dynamics.
Figure 3 shows the relation between the compactness parameters, ZAMS masses, and Si/O interface interior masses.Echoic in some sense of the plots in Figure 2, the differences are important.The accretion of the Si/O interface frequently kicks the core's mantle into explosion (Wang et al. 2022) and the near monotonicity of compactness with Si/O mass is suggestive.Figures 1, 2, and 3 set the initial structural context of our results.
2 Introduced to explore black hole formation, it is defined as , for which we prefer to set M equal to 1.75 M ⊙ when studying supernova explosions, not 2.5M ⊙ , as did O 'Connor & Ott (2011).

BASIC EXPLOSION RESULTS
Table 1 provides various final-state quantities for each of the twenty 3D progenitor evolutions that comprise this study.Included are the baryon mass, gravitational mass, explosion energy, 56 Ni mass, and kick speed at the end of each simulation (at post-bounce time t max )3 .Also given is whether we witness the formation of a neutron star or black Before shock revival, we set the "explosion energy" arbitrarily to zero.Upon shock revival, the matter behind the shock is still bound and early on most of the deposited neutrino energy goes into overcoming this negative gravitational binding energy and the (negative) contribution of the overburden exterior to the simulation grid (which is included in the accounting).This is why the explosion energy starts negative and only gradually grows to become positive over time.See text in §3.1 for a discussion.
hole, the ZAMS mass, and the compactness at 1.75 M ⊙ .It is this collection of numbers that informs our discussions to follow of the systematics found between stellar property, explosion characteristics, and observables.
To set the stage for the subsequent comparisons, we display in Figure 4 the temporal evolution of the averaged shock radius (top left), PNS baryon mass (top right), explosion energy (bottom left), and kick speed (bottom right).Note that explosion occurs over a range of post-bounce times from ∼100 milliseconds (ms) to ∼800 ms, with the lower mass progenitors generally exploding more quickly and the more massive progenitors taking longer.In fact, the shock waves for the 9 to 9.5 M ⊙ progenitors barely stall, while that for the 20 M ⊙ progenitor recedes significantly before exploding.We note that among those models for which the shock initially recedes it is only those for which the compactness and associated mass accretion rate and accretion luminosity are high that explode.The 12.25 and 14 M ⊙ models have too low a compactness to reverse the recession of their stalled accretion shocks.As Figure 4 indicates, the baryon mass asymptotes earlier than the energy and kick speeds.The latter both take many seconds to "flatten," except in the case of the lowest-mass and lowest compactness progenitors.We note that explosion energies can ramp up quickly in some MHD explosions (Obergaulinger et al. 2018;Müller et al. 2018;Obergaulinger & Aloy 2020).The explosion energy (here including the stellar mantle overburden binding energy) starts negative and can take more than a second to become positive.Most of the neutrino energy deposited early on by absorption goes into overcoming the gravitational term and only later is the positive value necessary for explosion achieved.Those numerical simulations in papers that stop within ∼1 second of bounce can not possibly provide an informed estimate of the explosion energy.The same can be said for the kick speed (Coleman & Burrows 2022;Burrows et al. 2023b) − as Figure 4 reveals, models halted within a second of bounce can not provide a useful final kick magnitude.

Time of Explosion
Figure 5 displays the relation between the compactness parameter and the explosion time.The explosion time in this work is defined as the time after which the shock never again recedes in radius.We see that the relation between the two is very roughly monotonic − models with larger compactness take longer to explode, but there is a lot of scatter.This scatter may in part reflect the chaos of the dynamics, but the degree to which this is true has yet to be determined.We emphasize that this theme of chaos, stochasticity, and scatter in the correlations suffuses our results and may reflect the very real scatter in the outcomes, even for the same progenitor.In fact, one expects a distribution function in the observables for the same progenitor, arising from the turbulent chaos in the dynamics.The width and form of such distribution functions has yet to be determined and is an important topic for future investigation.One of the motivations for this study of many models, as opposed to only a few, is the idea that in so doing one can glimpse the various correlations that would otherwise be obscured.We think we have achieved this in this paper, but clearly an even more expanded (and expensive!)investigation would be better.
A more interesting plot is Figure 6, which depicts the relationship between explosion energy (left) and cold, finalstate neutron star gravitational mass (right) and time to explosion.The behavior of the neutron star mass with post-explosion time is roughly monotonic and expected.The longer the delay to explosion, the more mass can be accumulated for a given accretion rate.Moreover, as Figure 5 reflects, the explosion time is greater for greater  accretion rate (compactness), making the trend all the more robust.However, the roughly increasing explosion energy with explosion time might seem surprising and this behavior is opposite to naive expectation.Nevertheless, the greater delay to explosion is generally (if not universally) associated with higher accretion-powered luminosities, which after the onset of explosion drive it more vigorously.This trend is hinted at in the bottom left panel of Figure 4.The aspherical morphology of the explosion debris reflects the turbulence of the pre-and post-explosion dynamics that has been a central feature of CCSN theory for decades (Burrows et al. 1995) 4 .Bubble structures are universally in evidence (see Sato et al. (2021)), with larger bubbles generally driving the explosions.Importantly, the longer the delay to explosion, the more vigorous the turbulence achieved, with the Mach number of the convective overturning motions prior to explosion getting as high as ∼0.5.However, if the explosion is more prompt, the post-shock turbulence does not have sufficient time to grow to such vigor and large convective structures are not produced.The result in such cases is a more spherical explosion, albeit with many small bubble structures.Given the trends of compactness and progenitor ZAMS mass with explosion time exposed in Figure 5, we would expect lower compactness, lower ZAMS mass progenitors would explode more spherically, while those with larger compactness would explode more aspherically.
This is qualitatively what we see.Figures 7 and 8 depict 56 Ni isosurfaces, colored by Mach number (red positive and blue negative), rendered in order of ZAMS mass.We see that the lower compactness progenitors explode more spherically, while those with higher compactness and longer explosion times explode more axially, with larger driving bubbles.In fact, the latter generally explode along an axis with a net dipolar and/or quadrupolar structure.Interestingly, we generically see simultaneous explosion in one set of directions, accompanied by accretion in the others (Burrows et al. 2020).This breaking of spherical symmetry is a core behavior of 3D CCSN models.The accretion in one set of directions continues to feed neutrino emissions that drive the explosion in the others, something that can't happen in 1D.Moreover, explosive nucleosynthesis is predominantly along the directions of the driving bubbles.Complementary to those directions, the shock is weaker and nucleosynthesis by the primary shock is reduced.The  upshot is the expectation that unburned oxygen will be concentrated very roughly perpendicular to the regions where most of the 56 Ni and explosive burning products reside.
Such structures have a direct bearing on the relative line profiles and the spatial distribution of the elements in supernova debris (Taubenberger et al. 2009;Sato et al. 2021;Fang et al. 2023) and the associated (anti)correlation of 56 Ni and progenitor oxygen shell material is a qualitative prediction of these 3D models.Moreover, one would expect The approximately monotonic behavior of the neutron star mass with all three of these quantities is a robust prediction of the theory.We note that Nakamura et al. (2015) derived a similar correlation between neutron star mass and compactness, but calculated in 2D, used the all-too-approximate "IDSA" neutrino scheme, and did not calculate beyond ∼1 second after bounce.A related set of correlations was anticipated by Ertl et al. (2016).
a correlation with the polarization of the inner debris (Leonard et al. 2000(Leonard et al. , 2006;;Tanaka et al. 2017;Nagao et al. 2023).Figure 9 portrays the relationship between the explosion energy and the dipole anisotropies of the ejecta mass (left) and the ejecta energies (right) among our set of twenty models.The dipoles are the amplitudes of the l = 1 spherical harmonic components of their angular distributions.These panels suggest a roughly monotonic correlation of two indices of anisotropy with explosion energy that is a new prediction that emerges from our study.We expect that, though with the common caveat of some scatter, more energetic explosions to evince greater anisotropy, but this needs to be tested.

Neutron Star Mass
The panels in Figure 10 show the dependence of the cold neutron star gravitational mass upon compactness parameter and upon initial central density (left) and the dependence of the neutron star gravitational mass upon the interior mass of the Si/O interface (right).The correlation between the compactness and the neutron star mass succinctly expresses the direct mapping between initial progenitor structure and an observable.
This correlation is stronger than with ZAMS mass, which reflects the non-monotonicity of ZAMS mass with collapsing structure in the Sukhbold et al. (2016) and Sukhbold et al. (2018) progenitor compilation.Since compactness and mass accretion rate are tightly aligned, and the mass accretion rate is the primary factor in the accumulated mass of the residue of collapse, this relationship is a natural prediction of the neutrino-driven mechanism of CCSNe.With this derived relationship between compactness and gravitational mass, one can use the compactness/ZAMS mass correlation (itself a function of the provenance of the progenitor models) to derive a birth mass function for neutron stars ( Özel & Freire 2016;Fan et al. 2023;Vartanyan et al. 2023).This should be a program for future work.

Explosion Energy
Two of the most important observables of CCSN theory are the explosion energy and the gravitational mass of the residual neutron star.Our theoretical results reveal that the two are tightly coupled.Figure 11 portrays the relation between the two and indicates that, with the caveat that some of our models have not quite asymptoted to their final explosion energies (but are generally within ∼10% or better), the explosion energy is an almost monotonic function of neutron star gravitational mass.This relationship between two observables is one of the central results of this paper and speaks to the importance of our philosophy to explore many long-term 3D state-of-the-art simulations to extract meaningful predictions.This strong mapping between observables would not have been clearly discerned without a 3D investigation of many progenitors to late times and should inform future observational campaigns to extract physically useful parameters of core-collapse supernova explosions.The study of but one single model would not be so illuminating.
Figure 12 plots the connection between binding energy and explosion energy (left) and compactness and explosion energy (right).Both relationships are roughly monotonic and robust, but the explosion energy/binding energy correlation is physically one of the most important findings from the theoretical perspective (Janka 2017a).It implies that the supernova energy is a "self-regulating" phenomenon, akin in broad outline to winds.Such self-regulation, mediated through the shared implicit dependence upon compactness, may become one of the unifying principles of CCSN theory going forward.
One of the central observables of core-collapse supernova explosions is the 56 Ni yield (in solar masses), since the 56 Ni decay chain to 56 Fe powers some of the supernova light curve and can easily be measured.Hence, the predicted 56 Ni yields are central physical quantities of the CCSN phenomenon.Figure 13 provides the relation that emerges in our 3D simulation suite between explosion energy (in Bethes ≡10 51 ergs) and 56 Ni production (in M ⊙ ). 56Ni production, and in general the ejecta isotopic abundances, are calculated using SkyNet (Lippuner & Roberts 2017) with a 1530-isotope network including elements up to A = 100, together with about 340,000 post-processed backwardly-integrated tracers in each simulation.More details concerning the formalism can be found in Wang & Burrows (2023b).This too is one of the central results of this paper, but the correlation between the two has already been anticipated theoretically (e.g.Sukhbold et al. (2016)) and observed astronomically (Rodríguez et al. 2023).The left hand side provides the 56 Ni correlation with explosion energy, with the range of predicted energies due to the fact that not all 3D models have yet asymptoted to their final values indicated.The right hand side shows the same numerical data, but superposed on the Rodríguez et al. (2023) measurements.The overlap between the two suggests a reasonable concordance, though both trends are subject to natural variation.Note that the strong theoretical trend of 56 Ni mass with explosion energy includes the exploding black hole formers (triangles).
Figure 14 shows the relationship between 56 Ni production, ZAMS masses, and compactness.While there is some scatter, the general monotonic trend is good, with two exceptions versus ZAMS mass.The first in the 19.56 M ⊙ exploding black hole former, which produces a lot of 56 Ni and the other is the 60 M ⊙ model, which is underenergetic and leaves behind a neutron star.It also experienced pronounced mass loss prior to explosion, which may not comport with the emerging understanding of stellar mass loss (Smith 2014).The deviations for both models from the trend is less pronounced versus compactness.Though we still witness some degree of stochasticity, the better (though not perfect) correlation with compactness than with ZAMS mass reflects again the non-monotonicity of the two, one with the other in the Sukhbold et al. (2016) and Sukhbold et al. (2018) progenitor collection.What other progenitor collections would provide is yet to be determined and is an important future question in supernova and stellar evolution theory.
Figure 15 shows the final ejecta electron fraction (Y e ) distributions of all our exploding models.The competition between ν e and νe absorption in the ejecta when near the inner core pushes much of the Y e above 0.5 into the protonrich domain.Most of the mass, however, is near 0.5.However, for some models (such as model 9(a)), which experience faster ejecta speeds during the wind and α-rich freeze-out phases, some neutron-rich starting material does not have enough time to be pushed above Y e = 0.5 before being ejected.The result is a component of the explosion debris  that emerges neutron-rich, with consequences for the nucleosynthesis.In particular, isotopes near the first peak of the r-process can be produced.There are many other consequences of this (Pruet et al. 2006;Wang & Burrows 2023a,b), but we prefer here to highlight only a few.
Figure 16 shows the dependence of the freeze-out mass5 upon compactness (top left), ejecta mass dipole (top right), and explosion energy (bottom).We see here too the rough monotonicity, in particular versus the explosion energy, of the mass in the freeze-out component.This follows the corresponding correlation of the explosion energy with 56 Ni mass depicted in Figure 13 and is not surprising.It is nevertheless a strong prediction of the theory of CCSN explosions, as we see it, and is testable.et al. (2023b) have recently crafted and had accepted a paper on the physics and systematics of the recoil kicks and induced spins, and their potential correlations, for our long-term 3D simulations.Though that paper contains the same models, the time between this paper and that has allowed us to carry those models even further.However, the general conclusions in Burrows et al. (2023b) are unaltered and we won't repeat its discussion and conclusions in detail here.We found an important correlation between the kick speed and both the energy of explosion and the ejecta mass dipole.We identified a weakly inverse relation between induced spin and kick speed, with a lot of scatter, and found a weak induced-spin/kick correlation.If the initial rotational vector of the collapsing core correlates somewhat with the direction of explosion, a spin/kick correlation would naturally emerge (Holland-Ashford et al. 2017;Katsuda et al. 2018;Ng & Romani 2007).We also found a compelling explanation for the lower average gravitational masses of binary neutron stars in the neutron star mass/kick connection that emerges from the simulations, viewed collectively.The lowest mass progenitors give birth to the lowest mass neutron stars with the slowest kicks, and this would make binary disruption harder.

Burrows
Therefore, in this paper we add only two more plots that reveal two additional correlations.Figure 17 indicates that the kick speed and the energy dipole of the ejecta (left) and the combination ( Eexp M baryon ) 1/2 (right) (where E exp is the explosion energy and M baryon is the baryon mass of the neutron star) are tightly correlated.The more asymmetrical the explosion the greater the kick speed.Moreover, the kick speed is related to the explosion energy itself, divided by the residual neutron star mass.In Burrows et al. (2023b), we found similar relations between kick speed and explosion energy and ejecta mass dipole, as did Janka (2017b)6 .However, the correlation between kick speed and ejecta energy dipole is actually a bit tighter.In any case, the connections between kick speed and various measures of anisotropy, while not unexpected, are testable.

CONCLUSIONS
In this paper, we have derived correlations between core-collapse supernova observables and progenitor core structures that naturally emerge from our extensive suite of twenty state-of-the-art 3D CCSN simulations carried to late times.This is the largest such collection of 3D supernova models ever generated and allows one to witness and derive patterns that might otherwise be obscured when studying one or a few models in isolation.Moreover, such long-term simulations are necessary to capture the asymptotic values of many of the most important observable quantities, and such simulations have until recently been too expensive to generate.However, with the efficiency of our new supernova code Fornax and recent generous allocations of supercomputer time, this limitation has now been overcome.Furthermore, a large collection of such models allows one partially to mitigate against the fact that the physical turbulence and chaos that attends supernova dynamics will naturally translate into a distribution of outcomes and explosion parameters in deriving trends.7 .Looking at single models one at a time tends to obscure the overall correlations, with both progenitor structure and between observables.
From this new panoramic perspective, we have identified 1) observationally useful correlations between explosion energy, neutron star gravitational birth masses, 56 Ni and α-rich freeze-out yields, and pulsar kicks (Burrows et al. 2023b) and 2) physically and theoretically important correlations with the compactness parameter of progenitor structure.The theoretical correlation between 56 Ni and explosion energy has been seen before (e.g.Sukhbold et al. (2016)), but not in the context of such a large 3D model suite carried to, or near to, the asymptotic state.Such a correlation is also strongly suggested by observations (Rodríguez et al. 2023).The intriguing correlation between explosion energy and neutron star gravitational mass is new in the detailed 3D modeling context, but was anticipated in the 2D results of Nakamura et al. (2015) and in the semi-analytic modeling of Müller et al. (2016a).In addition, we find a correlation between explosion energy and progenitor mantle binding energy, suggesting that neutrino-driven explosions are selfregulating.We also find a testable correlation between explosion energy and measures of explosion asymmetry, such as the ejecta energy and mass dipoles.
While the correlations between two observables are roughly independent of the progenitor ZAMS mass, the many correlations we have derived with compactness can not unambiguously be tied to a particular progenitor ZAMS mass.The compactness/ZAMS mass mapping depends upon the massive star progenitor models employed and the mapping between compactness and/or initial density structure and progenitor mass is still in flux.Therefore, our derived correlations between compactness and observables may be more robust than with ZAMS mass.There remain may important issues, such as binarity, shell mixing, overshoot, wind mass loss, semi-convection, the 12 C(α,γ) nuclear rate, and Roche-lobe overflow, that remain to be resolved related to the endpoint structures of massive stars.Hence, it is likely that the mapping between compactness and ZAMS mass found in the Sukhbold et al. (2016) and Sukhbold et al. (2018) compendia is not the final word.Using a simple analytic model for core-collapse supernovae and estimates of some explosion observables, this point is clearly made in Temaj et al. (2023).
In particular, we find that there is a range of progenitors in the Sukhbold et al. compilation that don't explode, but this gap in progenitor masses may be displaced from Nature's and from those in other progenitor model sets.Hence, a gap may indeed exist to explain one channel of black hole formation (Burrows et al. 2023a) and an associated measured luminosity gap (Smartt et al. 2009;Smartt 2015;Adams et al. 2017), but it may reside at a different range of progenitor masses than we have found.Nevertheless, the theoretical derivation of a range of ZAMS masses in which explosions are difficult may survive, but what that range is is likely still in play.
Though Fornax is a sophisticated, multi-physics multi-group radiation/hydrodynamics code that has been matured over the last eight years to address the complicated multi-dimensional phenomena encountered in core collapse, there are numerous upgrades and improvements that suggest themselves.We are bundling the ν µ , νµ ν τ , ντ neutrinos into one pseudo-species.This is done because their neutrino-matter interactions are quite similar and to trim the computational load.However, their interactions are not exactly the same (Bollig et al. 2017).Also, we are using 12 (×3) energy groups and though we have in the past determined that this number is adequate, more would be better.In the same spirit, an augmentation of the number of zones and enhanced spatial resolution are always good, if affordable.Currently, we are using a 1024(r)×128(θ)×256(ϕ) grid and this is better than most 3D production simulations.However, we have shown that resolution can make a difference, perhaps qualitatively (Nagakura et al. 2019).Importantly, we are using approximate general relativity (Marek et al. 2006), which captures the basic differences in the strength of gravity between the Newtonian and Einsteinian realizations, provides an accurate lapse function, and allows for the incorporation of gravitational redshifts.However, full general relativity is the correct approach.Also, we have not attempted to capture the effects of neutrino oscillations and we are wedded in this simulation suite to the Sukhbold et al. (2016) and Sukhbold et al. (2018) initial models.The latter are certainly not the last word, are not 3D, and do not include the many effects of binarity.Furthermore, as we have demonstrated, the initial structures are crucial to determining the mapping between star and outcome and we have not incorporated the effects of magnetic fields (Müller & Varma 2020;Varma & Müller 2021;Obergaulinger & Aloy 2020;Kuroda et al. 2020;Obergaulinger & Aloy 2021), nor addressed the effects of initial core rotation.Finally, the chaos of the turbulent context of post-bounce evolution introduces a degree of randomness and stochasticity that will translate into distribution functions for the observables, even for the same progenitor, whose widths are currently unknown.Hence, one can envision many potential upgrades to current computational platforms, and many necessary and more detailed investigations to pursue, in order to further improve the 3D theory of core-collapse supernova explosions.Nevertheless, our current collection of twenty state-of-the-art 3D models provides a synoptic view unavailable in the past upon which future work can profitably build.
and NERSC/Perlmutter.To this list of machines and architectures can be added local (though smaller) CPU and GPU clusters at Princeton.Both the CPU and GPU variants scale well on all machines out to ∼150,000 CPU cores and ∼480 GPU nodes.
Fornax is written in C. We use an MPI programming model for the CPU variant and an MPI/OpenMP-4.5/5.0 hybrid paralellism model for the GPU variant.The code employs spherical coordinates in one, two, and three spatial dimensions, solves the comoving-frame, multi-group, two-moment, velocity-dependent transport equations, and uses the M1 tensor closure for the second and third moments of the radiation fields Vaytet et al. (2011) 8 Three species of neutrino (ν e , νe , and "ν µ " [ν µ , νµ , ν τ , and ντ lumped together]) are followed using an explicit Godunov characteristic method applied to the radiation transport operators, but an implicit solver for the radiation source terms.This solver for the transport operators can be explicit since the speed of light and the speed of sound in the inner core are close and we are already timestep-limited by the hydro Courant condition (Just et al. 2015;O'Connor 2015;Skinner et al. 2019).By addressing the transport operator with an explicit method, we significantly reduce the computational complexity and communication overhead of traditional multi-dimensional radiative transfer solutions by bypassing the need for global iterative solvers.Radiation quantities are reconstructed with linear profiles and the calculated edge states are used to determine fluxes via an HLLE solver.In the non-hyperbolic regime, the HLLE fluxes are corrected to reduce numerical diffusion (O'Connor & Ott 2013).The momentum and energy transfer between the radiation and the gas are operator-split and addressed implicitly.
The hydrodynamics in Fornax is based on a directionally unsplit Godunov-type finite-volume method.Fluxes at cell faces are computed with a fast and accurate HLLC approximate Riemann solver based on left and right states reconstructed from the underlying volume-averaged states.Without gravity, the coupled set of radiation/hydrodynamic equations conserves energy and momentum to machine accuracy.With gravity, energy and momentum conservation, handled with source terms in the energy and momentum equations, is excellent before and after core bounce.
Gravity is handled in 2D and 3D with a multipole solver (Müller & Steinmetz 1995), where we generally set the maximum spherical harmonic order necessary equal to twelve.The monopole gravitational term is altered to approximately accommodate general-relativistic gravity (Marek et al. 2006), and we employ the metric terms, g rr and g tt , derived from this potential in the neutrino transport equations to incorporate general relativistic redshift effects.
In the interior, to circumvent Courant limits due to converging angular zones, the code deresolves in both angles (θ and ϕ) independently with decreasing radius, conserving hydrodynamic and radiative fluxes in a manner similar to the method employed in AMR codes at refinement boundaries.The use of such static-mesh refinement, or "dendritic grid," allows us to avoid angular Courant limits at the center, while maintaining accuracy and enabling us to employ the useful spherical coordinate system natural for the supernova problem.We have set up eleven nuclear equations of state in Fornax form, but focused on the SFHo EOS (Steiner et al. 2013) for this study.
A comprehensive set of neutrino-matter interactions are followed in Fornax and these are described in Burrows et al. (2006).They include weak magnetism and recoil corrections to neutrino-nucleon scattering and absorption (Horowitz 2002) and ion-ion-correlations, weak screening, and form-factor corrections for neutrino-nucleus scattering.For inelastic neutrino-electron scattering, we use the scheme of Thompson et al. (2003) and the relativistic formalism summarized in Reddy et al. (1999).Inelastic neutrino-nucleon scattering is now handled using the much more efficient, faster, and more accurate Kompaneets approach we pioneered (Wang & Burrows 2020).We note that most other supernova codes, with the exceptions of those fielded by the Garching (Buras et al. 2006) and Oak Ridge (Bruenn et al. 2020) groups, do not generally include inelastic redistribution.Neutrino sources and sinks due to nucleon-nucleon bremsstrahlung and electron-positron annihilation are included, as described in Thompson et al. (2000).We currently include a many-body structure factor (S A ) correction to the axial-vector term in the neutrino-nucleon scattering rate due to the neutrino response to nuclear matter at (low) densities below ∼10 13 gm cm −3 derived by Horowitz et al. (2017) using a virial approach.Horowitz et al. (2017) attach their fit to S A to that of Burrows & Sawyer (1998) at higher densities.
The code checkpoints periodically, and the I/O is fully parallel.In fact, the code uses parallel HDF5 writes to a single file with underlying MPIIO coordination.Each processor writes at a different offset into a single HDF5 file, which we have found useful for restarting checkpoints with different processor numbers and topologies.We have also optimized the striping on simulation runs to get further speed in the I/O and can dump data on local SSDs (if available) using 48 stripes.Twelve energy groups arrayed logarithmically per neutrino species works well for 3D simulations.As stated, the current main advantages of the Fornax code are its efficiency due to its explicit nature, its truly multi dimensional transport solution, and the interior static mesh derefinement in the core.For the long-term 3D simulations employed for this paper, we use 12 (×3) energy groups and a 1024(r)×128(θ)×256(ϕ) spatial grid.
We realized early in 2022 a factor of ∼2 speed-up in the CPU version of Fornax.This was the result of more aggressive inlining and vectorization, and of the incorporation of the Kompaneets scheme to handle inelastic scattering on nucleons (Wang & Burrows 2020).Also in 2022, we created, tested, and fielded the MPI/OpenMP-4.5GPU version of Fornax.To achieve this, we converted the entire code to C (it was originally a hybrid of C and Fortran).Our dendritic grid is now mapped to an 1D array which is linearly allocated in memory.The whole grid is now divided into three parts: a bulk segment (not receiving MPI messages) and a sending and a receiving segment, now all continuous in memory.First focussing on the hydro only, we added GPU pragmas, writing some functions to fit the GPU format and redefining variables to minimize local memory usage.All calculations are done on the GPUs, while communication is done on the CPUs.Only the sending and receiving components are mapped between the CPUs and GPUs, greatly reducing data transfer.The grid bulk is transferred to the CPU only to write dump and restart/checkpoint files.We executed the same porting steps for the radiation transport modules.This entailed a lot more work.Noting that statically allocated arrays can not be used as pointers on GPUs using OpenMP, we redefined them dynamically.We also rewrote the flux calculation.Previously, this was based on a 1D pencil layout and parallelization could only be done along the perpendicular two dimensions (in 3D).This has been altered so that all three dimensions can be parallelized.Though this entails a bit more calculation, the overall speed-up on the GPUs is improved.We also decreased the memory footprint of the gravitational potential calculation, and substituted "atomic add" for "array reduction."This enabled the use of more threads and, hence, is much faster.We changed the index order of all physical variables so that the fastest changing index is the grid cell index.This reduced the uncoalesced memory access on the GPUs.
Finally, in 2022 we optimized the grid-to-node decomposition mapping and the load balance figure-of-merit in anticipation of the much higher node counts anticipated in the future.Importantly, we have already achieved very good load balancing efficiency for all decompositions for both the GPU and CPU variants of Fornax.

Figure 1 .
Figure 1.Progenitor mass density profiles as a function of interior mass coordinate for the twenty models we employ for this 3D supernova explosion investigation from the Sukhbold et al. (2016) and Sukhbold et al. (2018) collection of unstable massive star cores.The colors correlate with ZAMS mass and the dashed curves are for models we see leave black holes, either explosively (the 19.56-and 40-M⊙ models) or quiescently.

Figure 2 .
Figure 2. Total stellar mass (left) that remains at collapse and the corresponding C/O and He core masses (right) as a function of the log10 of the ZAMS mass.On the right panel, the filled symbols are the C/O core masses (left label), while the unfilled symbols are the He core masses (right label).The 60 M⊙ progenitor has the same C/O and He core mass because it has ejected its outer envelope exterior to the C/O core in strong winds.The black dots in both panels are all the models in the Sukhbold et al. (2016) and Sukhbold et al. (2018) solar-metallicity set.The triangles are the black hole formers.

Figure 3 .
Figure 3.The dependence of the compactness on the log10 of the ZAMS mass (left) and the Si/O interface interior mass (right).The black dots identify the full Sukhbold et al. (2016) and Sukhbold et al. (2018) collection of solar-metallicity models, while the colored symbols indicate the models in this study.The colored triangles identify the models that we see form black holes.

Figure 4 .
Figure4.The temporal evolution of averaged shock radius (top left), PNS baryon mass (top right), the minimum explosion energy (bottom left, in Bethes ≡10 51 ergs), and kick speed (bottom right).Before shock revival, we set the "explosion energy" arbitrarily to zero.Upon shock revival, the matter behind the shock is still bound and early on most of the deposited neutrino energy goes into overcoming this negative gravitational binding energy and the (negative) contribution of the overburden exterior to the simulation grid (which is included in the accounting).This is why the explosion energy starts negative and only gradually grows to become positive over time.See text in §3.1 for a discussion.

Figure 5 .
Figure 5. Relation between the explosion time and compactness parameter for the set of twenty models in this study.Note the roughly monotonic behavior, with lower ZAMS mass and lower compactness models exploding more quickly, though with some scatter.See the text in §3.1 for a discussion.

Figure 6 .
Figure6.The behavior of the explosion energy and final-state neutron star mass with explosion time.For some models that haven't quite asymptoted, a range of explosion energies is given, from minimum to projected (likely).Note the roughly monotonic behavior of both the energy and neutron star mass with explosion time.See the discussion in §3.1.

Figure 7 .
Figure 7. Ejecta morphology of the 9(a), 9(b), 9.25, 9.5, 11, 15.01, 16, 17, and 18 M⊙ models.For most models, the snapshot is taken at the end of the simulation.The 10% 56 Ni isosurface is shown, colored by the Mach number of the radial velocity.Red means positive Mach number (outflow), while blue means negative Mach number (inflow).The light blue envelope shows the shock position.The shock is rendered as a bluish bounding veil.Unbound 16 O is distributed between the nickel isosurface and the shock.The only difference between models 9(a) and 9(b) is the imposition of slight velocity perturbations upon infall in the former (Wang & Burrows 2023a).

Figure 9 .
Figure9.Relation between explosion energy and 1) the ejecta mass dipole (left) and 2) explosion energy dipole (right) for the exploding models in this compilation.The mass and energy dipoles are the amplitudes of the l = 1 spherical harmonic components of their angular distributions.The roughly monotonic behavior of these indices of asymmetry with explosion energy is a prediction of the neutrino-driven model of core-collapse supernovae.

Figure 10 .
Figure 10.Left: Relation between the final-state neutron star gravitational mass and both the compactness parameter (filled symbol) and the initial central density (open symbol), and Right: the relation between the final-state neutron star gravitational mass and the Si/O interface mass.The approximately monotonic behavior of the neutron star mass with all three of these quantities is a robust prediction of the theory.We note thatNakamura et al. (2015) derived a similar correlation between neutron star mass and compactness, but calculated in 2D, used the all-too-approximate "IDSA" neutrino scheme, and did not calculate beyond ∼1 second after bounce.A related set of correlations was anticipated byErtl et al. (2016).

Figure 11 .
Figure 11.The explosion energy versus the cold, catalyzed neutron star gravitational mass.The range for some models reflects the remaining uncertainty of the final state energies.The monotonic, roughly linear, relationship emerges as one of the most interesting predictions of the neutrino-driven core-collapse model of supernova explosions.See §5.2 for a discussion.

Figure 12 .
Figure 12.The relation between explosion energy and binding energy exterior to 2.0 M⊙ (left) and to the compactness (right).The connection between binding energy and explosion energy is physically one of the most important findings from the theoretical perspective.The only major outlier is the 19.56 M⊙ exploding black hole former.See text in §5.2 for a discussion.

Figure 13 .
Figure 13.Relations between 56 Ni (in M⊙) and explosion energy (in Bethes), with the range in the former indicated for some models that have yet to completely asymptote.The triangles are for the exploding black hole formers, yet they still follow the nearly monotonic curve.The right hand side includes the data of Rodríguez et al. (2023) and suggests by the overlap of the two a reasonable concordance.Similar observational correlations are found in Müller et al. (2017) and Pejcha (2020).

Figure 14 .
Figure 14.The relationship between the ejected 56 Ni mass and 1) the log10 of the ZAMS mass (left) and 2) the compactness parameter (right).See text in §5.3 for a discussion.

Figure 15 .
Figure15.Ejecta Ye distributions for our exploding models.Most of the ejecta are at Ye = 0.5 or higher, in the protonrich regime.However, a few of the models, in particular the 9(a) M⊙ model that boasted imposed initial perturbations in the progenitor, have an interesting neutron-rich tail.Such a tail can give rise to the first peak on the r-process (seeWang & Burrows  (2023a,b)).

Figure 16 .
Figure16.Relations between the freeze-out ejecta mass and the compactness parameter (top left), ejecta mass dipole (top right), and explosion energy (bottom).The reasonably robust correlations between the freeze-out mass and both the explosion energy and the ejecta mass dipole are predictions of the theory that can be tested.See the text at §5.3 for a discussion.

Figure 17 .
Figure 17.Relations between the neutron star kick speed and the ejecta energy dipole (left) and the quantity ( Eexp M bar ) 1/2 (right), where Eexp is the explosion energy and M bar is the baryon mass of the neutron star.See Burrows et al. (2023b) for a more thorough discussion.

Table 1 .
Summary of simulation outcomes and characteristics.The range of explosion energies provided is from minimum to projected (likely) and reflects the fact that not all models have completely asymptoted.See text for details.