# The Essential Character of the Neutrino Mechanism of Core-Collapse Supernova
Explosions

Tianshu Wang 1 David Vartanyan2, Adam Burrows1, Matthew S. B. Coleman1
1Department of Astrophysical Sciences, 4 Ivy Lane, Princeton University,
Princeton, NJ 08544, USA
2Department of Astronomy, University of California, Berkeley, CA 94720, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Calibrating with detailed 2D core-collapse supernova simulations, we derive a
simple core-collapse supernova explosion condition based solely upon the
terminal density profiles of state-of-the-art stellar evolution calculations
of the progenitor massive stars. This condition captures the vast majority of
the behavior of the one hundred 2D state-of-the-art models we performed to
gauge its usefulness. The goal is to predict, without resort to detailed
simulation, the explodability of a given massive star. We find that the simple
maximum fractional ram pressure jump discriminant we define works well
$\sim$90% of the time and we speculate on the origin of the few false
positives and false negatives we witness. The maximum ram pressure jump
generally occurs at the time of accretion of the silicon/oxygen interface, but
not always. Our results depend upon the fidelity with which the current
implementation of our code Fornax adheres to Nature and issues concerning the
neutrino-matter interaction, the nuclear equation of state, the possible
effects of neutrino oscillations, grid resolution, the possible role of
rotation and magnetic fields, and the accuracy of the numerical algorithms
employed remain to be resolved. Nevertheless, the explodability condition we
obtain is simple to implement, shows promise that it might be further
generalized while still employing data from only the unstable Chandrasekhar
progenitors, and is a more credible and robust simple explosion predictor than
can currently be found in the literature.

###### keywords:

stars - supernovae - general

††pubyear: 2022††pagerange: The Essential Character of the Neutrino Mechanism
of Core-Collapse Supernova Explosions–References

##  1 Introduction

The astrophysics of core-collapse supernova (CCSN) explosions has achieved a
significant degree of sophistication since the pioneering work of Colgate &
White (1966). In that paper, the connection between the neutrino agency of
explosion and core mantle heating by the neutrinos that are copiously produced
due to core dynamical compression upon achieving the Chandrasekhar instability
at the termination of stellar evolution was first hypothesized. However, the
neutrino physics and transport, the initial stellar model, and the nuclear
equation of state (EOS) were all quite primitive. Arnett (1966) improved upon
all these factors, but did not find a robust explosion. In the interim, the
neutral current and $\tau$ neutrinos ($\nu_{\tau}$) were discovered and the
nuclear EOS was made much more sophisticated, a process which continues to
this day.

It was not until the work of Wilson (1985) and Bethe & Wilson (1985) that the
“delayed” neutrino-driven mechanism was introduced. In this model, the stalled
shock generated upon core bounce after collapse and infall is reenergized
after a post-bounce delay of hundreds of milliseconds by neutrino heating
behind the shock in the so-called “gain region.” The latter region of low, but
importantly non-zero, neutrino absorption optical depth lies between the
electron-type ($\nu_{e}$) and anti-electron ($\bar{\nu}_{e}$) neutrinospheres
and the shock. This model in broad outline is our current paradigm, but as it
was originally conceived has many limitations. It was a spherical 1D model
that did not explode without a boost in the $\nu_{e}$ and $\bar{\nu}_{e}$
neutrino luminosities of nearly $\sim$25% due to inner core convection driven
by fiat and a mixing-length prescription by a “neutron-finger” (aka “salt-
finger”) instability. Such an instability, however, is not seen nor expected
in modern proto-neutron-star (PNS) simulations (Bruenn & Dineva, 1996).
Lepton-driven Ledoux convection in the PNS is in fact seen in modern multi-D
simulations (Dessart et al., 2006; Nagakura et al., 2020), but does not boost
the $\nu_{e}$ and $\bar{\nu}_{e}$ neutrino luminosities by the required amount
and is not a central factor in reenergizing the shock in modern theory. The
central factor in driving a stalled shock into explosion are now thought to be
neutrino energy deposition/heating in the gain region, aided most of the time
by neutrino-driven convection (Herant et al., 1994). The latter contributes a
turbulent stress (Burrows et al., 1995) behind the shock and increases the
dwell time in the gain region (Murphy & Burrows, 2008) that facilitates net
neutrino heating, with the added stress likely more important than the
increased dwell time. This overall picture is now strongly buttressed by the
modern concordance in the results of many recent state-of-the-art 3D
simulations (Lentz et al., 2015; Janka et al., 2016; Roberts et al., 2016;
Müller et al., 2017a; O’Connor & Couch, 2018; Vartanyan et al., 2019a; Ott et
al., 2018; Summa et al., 2018; Glas et al., 2019; Burrows et al., 2019;
Nagakura et al., 2019; Burrows et al., 2020; Kuroda et al., 2020; Bollig et
al., 2021; Vartanyan et al., 2021; Burrows & Vartanyan, 2021). However,
despite this palpable progress, it has proven difficult to convert this
qualitative picture into a quantitative and predictive condition without these
sophisticated 3D radiation/hydrodynamic simulations. There have, however, been
numerous attempts to do so, all of which are lacking in some important
aspects.

The first attempt at this was the luminosity ($L_{\nu}$)/mass-accretion rate
($\dot{M}$) critical condition of Burrows & Goshy (1993). These authors noted
that above a critical $L_{\nu}$/$\dot{M}$ curve no hydrostatic structure,
including the core and the mantle interior to the shock, can be found.
Increasing the luminosity above this curve led to a bifurcation, and did not
merely increase the stalled shock radius quasi-statically. Furthermore, these
authors speculated that the imaginary part of the frequency of the oscillation
of the fundamental radial mode of shock pulsation would change sign into
instability in traversing this curve. That pulsation timescale is roughly
$\frac{\Delta{M_{g}}}{\dot{M}}$, where $\Delta{M_{g}}$ is the mass in the gain
region. This is also roughly the advection time between the shock and the
inner core in accretion flow. The supernova explosion solution itself was a
dynamical transition between two steady eigenfunction solutions $-$ the quasi-
static accretion phase and a steady-state neutrino-heated wind phase 111This
bifurcation was identified as a “dual-cusp catastrophe” (Burrows, 2013), but
that suggestion has not been further developed.. In any case, though there
seems to be a core of truth in the concept of a critical curve, there was no
analytic or physics-based derivation of the curve itself. Such requires a
sophisticated model of the luminosity, and the latter requires radiative
transfer, neutrino-matter microphysics, and a time-dependent model employing
PDEs $-$ in short, a detailed numerical model. Therefore, though the critical
condition remains a robust metric for determining explosion, it requires
simulation results and is not an ab-initio predictor of explosion.
Nevertheless, the concept of the crossing of a critical curve upon reaching
the explosive state has gained traction in the literature, with numerous
papers attempting improvements.

For example, Summa et al. (2016) expanded this criterion to include multi-
dimensional turbulence, explicit reference to the proto-neutron star mass, and
an estimate of the gain radius. Summa et al. (2018) and Janka (2012) (see also
Janka et al. 2016) generalized this prescription to ostensibly include some of
the the effects of rotation. Pejcha & Thompson (2012) and Raives et al. (2018)
introduced the “antesonic condition,” deriving from a simple model a
relationship between the ratio of the isothermal sound and the escape speeds
behind the stalled shock above which no steady state solution exists. They
later extended this to include rotation and turbulence and factored in the
adiabatic sound speed (Raives et al., 2021). Murphy & Dolence (2017)
introduced the “Force Explosion” condition Gogilashvili & Murphy (2021), which
uses two radiation parameters, one comparing the neutrino power deposited in
the gain region with the accretion power and the second explicitly factoring
in the neutrino optical depth in the accreted matter near the neutron star
surface. Using this approach, they rederive the antesonic and critical
conditions and attempt a more general perspective.

Inspired by the original critical condition concept, but wanting to produce a
condition from which one can predict the outcome of collapse from the
progenitor density profile alone (i.e., determine “explodability”), Ertl et
al. (2016) derived a curve in $\mu_{4}$ and $\mu_{4}\times M_{4}$ space, where
$\mu_{4}$ is a measure of the slope of the mass density around the jump
(frequently near the silicon/oxygen [Si/O] interface) in specific entropy to a
value of four (per baryon per Boltzmann’s constant) and $M_{4}$ is the
interior mass at that jump. There is a corresponding jump in density at that
same approximate position and we will return to its importance in the body of
the paper. Only the initial progenitor mass density profile is employed by
Ertl et al. (2016) to determine explodability. Calibrating on their estimates
of the explosion energy and progenitor mass of the Crab and SN1987A
supernovae, they derive a critical curve above which models would explode. In
their construction, $\mu_{4}$ and $\mu_{4}\times M_{4}$ are related indirectly
to $L_{\nu}$ and $\dot{M}$, so the heritage in Burrows & Goshy (1993) is
manifest. However, as Müller et al. (2016) suggest and we advance below, the
Ertl et al. (2016) critical curve doesn’t in fact seem to be a robust index of
explodability.

Quantities which have proven important in characterizing the potential
explodability of the quasi-stationary shock through which matter is accreting
onto the PNS core and behind which neutrinos are heating that same matter are
the advection timescale, $t_{adv}$, and the heating timescale, $t_{heat}$
(Thompson, 2000; Murphy & Burrows, 2008). In some senses, it is the
competition between these quantities that determines instability to explosion
and the potential for a dynamical transition from the steady state. Sample
definitions of these quantities are:

$t_{adv}(\sim t_{dwell})\approx\frac{\Delta{M_{g}}}{\dot{M}}\,,$ (1)

where $\dot{M}$ is the accretion rate and $\Delta{M_{g}}$ is the gain region
mass, and

$t_{heat}=\frac{\lvert{E_{g}}\rvert}{\dot{Q}_{heat}}\,,$ (2)

where $E_{g}$ could be either the total energy, the gravitational energy, or
thermal energy in the gain region and $\dot{Q}_{heat}$ is the neutrino power
deposition rate in that same region. Equations for these quantities have
varied in the literature and have not been made rigorous. However, if in a
general sense $t_{adv}$ is larger than $t_{heat}$, then neutrino power
deposition can heat the gain region before advection through this region
brings the accreting matter into the inner core interior to the gain radius
where there is net cooling. It is this cooling region just above the
approximate positions of the $\nu_{e}$ and $\bar{\nu}_{e}$ neutrinospheres
that undermines pressure support for the outer mantle matter and imperils
explosion. A quick traverse through the gain region undermines explosion,
while a long dwell time might facilitate it.

However, the astrophysics of this is a bit more subtle. Even if $t_{adv}$ were
infinite, continued neutrino heating in the gain region would quickly saturate
to an equilibrium between heating and cooling, the latter due to neutrino
emission via $e^{-}$ and $e^{+}$ capture on free nucleons. This cooling rate
is roughly proportional to the temperature ($T$) to the sixth power, and given
the roughly inverse-square neutrino flux dilution, the temperature would
equilibrate to a $\frac{1}{r^{1/3}}$ radial profile. In trying to
predict/explain the solar wind, Parker (1958) showed that an ideal-gas
atmosphere with a temperature profile shallower than $\frac{1}{r}$ could not
maintain hydrostatic equilibrium unless there was a finite pressure at
infinity. The absence of such a bounding pressure implied to Parker (1958)
instability to the solar wind. In our case, there is a bounding pressure due
to the accretion ram. Then, the question becomes: given the neutrino heating
rate (which leads to a temperature profile) is that ram pressure low enough to
lead to instability $-$ to an explosion? Since the heating rate
$\dot{Q}_{heat}$ depends upon the $\nu_{e}$ and $\bar{\nu}_{e}$ luminosities
(times roughly the mean neutrino optical depth, $\tau_{\nu}$, in the gain
region, of order $\sim$0.1) and the ram pressure depends upon $\dot{M}$ at the
shock, comparing $t_{adv}$ and $t_{heat}$ to derive a condition for explosion
is roughly equivalent to the classic critical condition pitting $L_{\nu}$ and
$\dot{M}$ (Burrows, 2013).

These considerations are qualitatively useful and provide insight into the
nature of the neutrino-driven mechanism of CCSN explosion, but they are not
predictive. To determine which massive-star progenitors explode and to derive
a working and credible explosion criterion one must have a detailed model for
the luminosities and heating, for the equation of state, for actual neutrino-
matter microphysics, and for the turbulence. Moreover, one must be able to
gauge the feedback influence of the accretion rate on the luminosities and
neutrino spectra themselves. Hence, short of detailed radiation/hydrodynamic
supernova simulations, even informed analytic treatments are likely to fall
short222Though, one may be able to find interesting potential correlations
between observables (Müller et al., 2016).. Therefore, while serving to
highlight important aspects that factor into explodability, all published
approximate criteria for explodability fail to be usefully predictive.

The compactness parameter of the initial progenitor density profile (O’Connor
& Ott, 2011, 2013) has perhaps been the most misapplied explosion metric. It
is defined as (O’Connor & Ott, 2011):

$\xi_{M}=\frac{M/M_{\odot}}{(R(M)/1000\,\mathrm{km})}\,,$ (3)

where $M$ denotes the interior mass at which the compactness parameter is
evaluated 333originally defined at 2.5 M⊙, but our preferred value is 1.75
M⊙). $\xi_{M}$ has been used as a parameter below which explosions would occur
and above which they would not. However, this dichotomy has been shown to be
false, with both low and high compactness models exploding (Burrows et al.,
2020; Burrows & Vartanyan, 2021)444though, perhaps the highest $\xi_{M}$
models abort explosion and reimplode to black holes.

Moreover, within the context of the neutrino-driven paradigm, higher
compactness seems necessary to achieve explosion energies in the $\sim
10^{51}$ ergs (1 Bethe) range (Burrows & Vartanyan, 2021). It is only such
models that have enough neutrino absorption optical depth behind the shock for
a sufficiently long period of time to absorb sufficient energy during launch.
Such models also have more core ejecta mass whose recombination energy is
thereby usefully high.

The advantage of compactness is that it is an ab-initio, simple-to-calculate
parameter that would not require running detailed core-collapse simulations.
Moreover, there is a strong correlation of compactness with post-bounce
neutrino luminosity and accretion rate and, for exploding models that birth
neutron stars, there is a strong positive correlation with residue mass. If a
black hole forms, the timescale for its formation is inversely related to
compactness, since accretion rate increases with compactness. But, compactness
per se does not serve as a useful explosion discriminant.

Finally, for most progenitors the mass at which $\xi_{M}$ is calculated also
encompasses the silicon/oxygen (Si/O) interface entropy and density jumps and,
as has been shown in the past (Sukhbold et al., 2016; Ertl et al., 2016;
Burrows et al., 2018; Vartanyan et al., 2018; Burrows et al., 2019; Burrows &
Vartanyan, 2021; Vartanyan et al., 2021), this jump is an important feature
that can affect explodability. This topic will be more extensively discussed
in this paper.

Given the degree to which attempting to derive a simple explodability
criterion in the past has been fraught, it is reasonable to cast a jaundiced
eye on attempts to do so now, or in the near future. There is still much to do
with regards to the initial progenitor models, the neutrino-matter
microphysics, numerical algorithms, and computational implementations. We do
not think that current sophisticated models have yet converged to either 1)
the ultimate answer or 2) the definitive mapping between progenitor and
outcome. Nevertheless, there has of late been sufficient progress in the
multi-dimensional modeling of collapse and explosion to take a fresh look,
guided by results emerging from the more detailed computational models. In
this paper, we acknowledge the physical insights provided by the past
literature on explodability criteria, but seek to derive a means to predict
explodability given only the initial mass density profiles calibrated on a
large set of two-dimensional (axisymmetric) CCSN simulations using our code
Fornax (Skinner et al., 2019). A large suite of three-dimensional simulations
with which to benchmark success is still too expensive and time-consuming. We
have seen that, though the details of 2D and 3D evolutions can be different,
with very few exceptions when one of our Fornax models explodes in 2D it does
so in 3D, and vice versa. Under this assumption, and with the presumption that
our simulations are capturing much of the essential truth about explodability
(our major and most fragile assumption), we in this paper take up the
challenge to develop a predictive explosion condition.

First, we assembled many of our recent 3D and 2D models and explored with them
various possible explodability conditions and metrics using only the
associated progenitor mass density profiles with which those models were
inaugurated. Settling on one measure (and allowing some exceptions), we then
picked 100 massive-star progenitor models from the compilation of Sukhbold et
al. (2016) and Sukhbold et al. (2018) and performed 2D CCSN Fornax simulations
for each of them for many seconds post-bounce to determine whether they
exploded by the criterion we arrived at during the earlier exploratory phase.
Indeed, we found that the explodability criterion we “predicted” comported
with the results for the 100 models, with few exceptions. This gave us
confidence in the relative usefulness of the simple approach to determining
explodability we found.

In §2 of this paper, we first set the stage by summarizing various physical
aspects of the neutrino mechanism rarely found written in one place. We feel
that articulating, even in bulletized form, the essential elements of the
modern neutrino mechanism of CCSN explosion will provide context and might
prove useful. Then, in §3 we summarize the results of the exploratory set of
3D and 2D models and the systematic behavior of this “training set” with
progenitor characteristics. The associated simulations are taken from Burrows
et al. (2020) and Burrows & Vartanyan (2021). In §4, we motivate our choice of
central explosion condition using the exploratory models, and in §6 we present
our results for the 100 2D models, discuss caveats, and advance our
approximate approach to determining explodability. We conclude in §7 with a
summary of what we have found and a listing of the strengths and potential
weaknesses of our findings.

##  2 Physics of the Neutrino Mechanism

We provide here some of the basic physical characteristics of the neutrino
mechanism of core-collapse supernova explosions. There has been no attempt to
fully flesh-out these features, nor to justify our choices. Rather, we hope
that a list of clear declarative points, without the baggage of detail, might
help orient the reader in advance of the more focussed discussion that
follows. We admit that our compilation might be incomplete and is a bit
idiosyncratic. However, we know that this decades-old subject can at times
seem confusing and less coherent than perceived by long-time practitioners. In
fact, there is significant agreement on many aspects of the theory that is
emerging.

  * •

The onset of the core-collapse supernova explosion is a bifurcation between a
quasi-steady PNS solution with a bounding accretion shock wave and a mantle
ejection instability. There are no intermediate stable solutions. The
asymptotic state of the core is one with a progressively weakening neutrino-
driven wind, whose onset depends upon the evolving pressures and ram pressures
of the material of the inner explosion (Burrows, 1987; Burrows & Goshy, 1993).
Should the latter be large, the wind phase will be delayed. The supernova is a
dynamic transition solution between the quasi-steady accreting PNS and wind
eigenvalue solutions. The wind-phase mechanical power is lower than that
experienced during the supernova phase due to the much larger absorbing mass
and degree of neutrino power deposition during the latter.

  * •

It is thought that the explosion condition might be derivable from a suitable
pulsational analysis which would lead to eigenfrequencies whose fundamental
(and dipole? Dolence et al. 2013) modes would experience at explosion a change
in the sign of their imaginary parts. This has not been definitively explored
in a realistic context, but seems compelling.

  * •

Gravitational energy released due to collapse from a Chandrasekhar white dwarf
to a proto-neutron star is the ultimate source of the explosion energy 555Here
and below, when we refer to the explosion energy, we always mean the sum of
all its contributions: the internal, recombination, kinetic, gravitational,
and total binding energy of the off-grid stellar envelope (negative). The use
by some of what is called the “diagnostic energy” is misleading and can be
confusing, and we recommend it never be invoked.. However, though the total
binding energy of a neutron star is some multiple of 1053 ergs, most of this
energy is radiated by neutrinos over tens of seconds. Only a fraction is
radiated during the launch and powering of the explosion and only about half
of that is radiated in the $\nu_{e}$ and $\bar{\nu}_{e}$ neutrinos that
dominate neutrino heating. Therefore, the neutrino mechanism is not merely
$\sim$1% (or less) efficient, but closer to $\sim$5$-$10% efficient (Bruenn et
al., 2016; Bollig et al., 2021; Burrows & Vartanyan, 2021).

  * •

Nuclear burning in the inner mantle shells of the progenitor contributes at
most $\sim$5$-$10% to the total explosion energy 666The energy available from
the burning of oxygen to 56Ni is $\sim$0.1 Bethes per 0.1 M⊙ of oxygen. The
amount of oxygen that can be burned to, for instance, 56Ni in the extended
core generally scales with the compactness. The fraction of this oxygen that
will explosively be burned scales with the neutrino-driven explosion energy.
Hence, since we argue that the latter scales roughly with compactness, the
total nuclear energy available very roughly scales with compactness.
Therefore, the fractional contribution of nuclear energy to the total is very
roughly constant. Since the measured amount of ejected 56Ni varies from $\sim$
a few$\times$10-3 to $\sim$0.1 M⊙, and the inferred CCSN explosion energies
vary from $\sim$0.1 to $\sim$2 Bethes, one obtains the roughly $\sim$5$-$10%
quoted (Müller et al., 2017b)..

  * •

Unlike Type Ia thermonuclear explosions, the mechanism of core-collapse
supernovae need not disassemble the residual core to leave nothing behind. The
bound neutron star early created in collapse quickly becomes too bound to
disassemble. Rather, the ongoing and prodigious neutrino losses from the PNS
“neutrino star” that leaves a progressively more bound PNS transfers energy to
its mantle and it is the mantle that explodes. It is this energy transfer from
core to mantle and the latter’s ejection that is the essence of the neutrino
mechanism. During the launching of the shock wave much of the deposited
neutrino energy is used to do work against gravity.

  * •

Hence, the CCSN mechanism is an energy transducer from gravitational to core-
thermal to neutrino to mantle-thermal to (mostly) kinetic energy.

  * •

Furthermore, the gravitational energy of infall that is used to dissociate the
nuclei accreted through the stalled shock is retrieved upon the recombination
of the ejected nucleons. This source (recall that the specific binding energy
of nuclei is near $\sim$10 MeV/baryon) constitutes a large fraction of the
total asymptotic supernova explosion energy.

  * •

Radiation pressure is not a significant factor in the neutrino mechanism of
explosion. The neutrino Eddington luminosity is near $\sim$1055 erg s-1, far
in excess of even the peak $\nu_{e}$ luminosity at breakout (which itself
lasts only $\sim$10 milliseconds).

  * •

The cumulative energy deposited by neutrinos in the gain region and mantle
interior to the stalled shock bears no relation to the explosion energetics.
In fact, the total mantle energy is negative at the onset of explosion and has
been secularly decreasing during the delay to explosion. Afterwards, neutrino
heating and recombination together combine over time after the delayed
instability to achieve the positive asymptotic energy of the supernova.

  * •

It is thought that, in broad outline, the energy scale of CCSN explosions is
ultimately and approximately set by the total binding energy of the stellar
mantle exterior to the Chandrasekhar core. This is analogous to the self
regulation seen in stellar winds, for which the escape velocity from the star
is roughly the asymptotic velocity of the wind at infinity (Vink, 2011). This
implies that the greater this binding energy, the greater the explosion
energy, unless the progenitor does not explode (perhaps due to an over-large
binding energy that can’t be overcome by the available neutrino heating). This
feedback concept is distantly analogous to that for a stellar wind, for which
asymptotic speeds scale with surface escape speeds.

  * •

Neutrino heating of the mantle drives turbulent convection between the stalled
shock and the PNS core. This is akin to boiling water on a stove. The
turbulent stress contributes an effective additional pressure that enlarges
the gain region and increases the stalled shock radius. The larger gain region
absorbs slightly more neutrino energy and the larger shock radius puts some of
the gain region mass lower in the gravitational potential well. Thus less
bound, this mass is more easily unbound.

  * •

The stress tensor of the turbulence is anisotropic (Murphy & Burrows, 2008;
Murphy et al., 2013), with the radius-radius component in the radial direction
larger than the other components. Moreover, the effective $\gamma$ connecting
turbulent energy with turbulent pressure is near $\sim$2 (compare with
$\sim$4/3), making turbulence a more efficient means of realizing pressure
with energy density than thermal gas.

  * •

$t_{adv}$ is also relevant to the strength of the turbulence, since a large
value facilitates the growth of the seed perturbations in the pre-shock
infalling matter to the non-linear amplitude regime before settling into the
inner core by increasing the number of e-foldings (Foglizzo et al., 2007). The
magnitude and spatial scales of the seed perturbations in the progenitor
convective zones are not currently well-known, but can factor into the
explosion details. Current CCSN simulations witness turbulent Mach numbers
between $\sim$0.1 and $\sim$0.5.

  * •

$\nu_{\mu}$ and $\nu_{\tau}$ neutrino emission constitutes $\sim$50% of the
total binding energy losses. However, their contribution to heating in the
gain region is subdominant. Their major positive effect on explodability is
indirectly by helping to drive the Kelvin-Helmholtz shrinkage of the core
through their radiation, which leads to an increase in the temperatures of the
$\nu_{e}$ and $\bar{\nu}_{e}$ neutrinospheres and thereby to a hardening of
the emergent $\nu_{e}$ and $\bar{\nu}_{e}$ spectra. Such a hardening leads to
higher neutrino absorption and heating cross sections in the gain region.

  * •

The breaking of spherical symmetry enabled in 2D and 3D allows simultaneous
accretion and explosion, the former maintaining neutrino accretion power that
can continue to energize the latter. This can’t happen in 1D (spherical) and
is an important reason, along with the modeling of turbulence, CCSN
simulations need to be done in multi-D.

##  3 Training Set of 3D and 2D Core-Collapse Models

As a prelude to this study, we investigated a set of legacy 3D (Burrows et
al., 2020) and 2D (Burrows & Vartanyan, 2021) Fornax simulations777These 3D
and 2D studies were each done for a different set of progenitors and the
associated papers describe the numerical and computational setups. The early
evolution of the maximum shock radius after bounce for these models is given
in Figure 1. The solid lines are for models that exploded and the dashed lines
are for models that did not. The filled circles indicate the time when the
Si/O interface (or the interface with a large density drop) accretes through
this maximum shock radius. As is clear, for many models the time of explosion
is roughly when (or near when) this occurs, and if a progenitor from this set
explodes it generally does so after $\sim$70 to $\sim$350 milliseconds (ms)
after bounce and the inauguration (at “$t=0$”) of the quasi-steady delay
phase. Note that the models that don’t explode (at least among these simulated
models) have a low to intermediate ZAMS mass and that this model suite does
not include ZAMS masses above 26 M⊙. We leave more massive progenitors to a
future study. 888Discussion on why progenitors in a certain mass range are
less explosive can be found in Section 6.

The associated initial baryon mass density profiles of this model set are
given in Figure 2, with the dots identifying the approximate positions of the
Si/O interfaces. It is our suggestion, shared by others, that, under the
assumption of adequate microphysics and algorithms, this density profile is
the major determinant of explosion for non-rotating progenitors. This
perspective emphasizes the central role of stellar evolution prior to the core
Chandrasekhar instability at the terminal phase of a massive star. However,
there exist in these profiles numerous features that will in principle affect
the competition between “$L_{\nu}$” and $\dot{M}$ or between $t_{adv}$ and
$t_{heat}$. Among these are the depths and positions of the various shell
interfaces, the central densities at collapse, and the slopes of the density
profiles along the smooth segments throughout the inner region interior to
$\sim$2.5 M⊙. There are in principle numerous relevant parametrizations for
these profiles. What might be the correlation between any one of these
quantities and explodability? This is what we set out to research, though we
again emphasize that our conclusions are constrained by whatever limitations
our current Fornax simulations may embody. We also note that the wide range in
general structures, from steep outer profiles for the less massive progenitors
to much more shallow outer profiles for the more massive progenitors, is not
rigorously monotonic with ZAMS mass, but only crudely so (Sukhbold et al.,
2018). The mapping between ZAMS mass and terminal density profile may be as
subtle and “chaotic” as the current literature suggests, but the general range
of density profiles depicted in Figure 2 may reflect reality. However, one
should keep in mind that the detailed mapping between ZAMS mass and terminal
profile is still a work in progress (Chatzopoulos et al., 2016; Müller, 2016;
Müller et al., 2017a; Müller et al., 2019; Fields & Couch, 2020, 2021).

Given this, we note some of the interesting systematic dependences we do see
in the 2D model set between residual PNS baryon mass and asymptotic explosion
energy (y-axes) and the compactness and the ZAMS mass (x-axes). These are
portrayed in Figure 3. We see that the final PNS mass is strongly correlated
with compactness ($\xi_{1.75}$) and, with a bit more noise, with ZAMS mass,
while the explosion energies are less correlated with either. However, there
is a hint that low-mass, low-compactness progenitors are less energetic than
high-mass, high-compactness progenitors. The former trend is to be expected,
since the post-bounce mass accretion rates increase reliably with compactness,
and the latter is reasonable, since high compactness boosts both the driving
accretion neutrino luminosities and the gain region mass upon which a fraction
of the emerging neutrinos are absorbed.

Figure 4 renders the dependence of the explosion energy on the mantle binding
energy and suggests that this binding energy does indeed help to set the
energy scale of explosion (when models do explode). This figure simultaneously
recapitulates the trends with compactness and ZAMS mass in Figure 3,
emphasizing again the systematics, loose or tight, between these quantities.
Figure 4 is an important possible finding of modern CCSN theory, but should be
viewed as provisional until verified more broadly. The binding energy versus
interior mass for the host of 3D and 2D models of our “training set” are given
in Figure 5. This figure suggests that there is a general expectation
connecting the “mass cut” (interior to which all materials will be left behind
after the explosion) and the explosion energy in the context of the possible
correlations depicted by Figures 3, 4, and 5.

Using these heritage 2D and 3D models, we attempted to find indices of
explodability that might be used to map supernova explosion outcome (a binary
“yes/no”) with initial density profile (Figure 2). We explored numerous
approaches, but settled on one which we now describe in §4. This condition,
while not perfect, seems to have predictive value, though there are a few
exceptions that need explaining.

##  4 Motivation for Our Explosion Criterion

Motivated by the general concept of the critical condition between luminosity
and mass accretion rate, we extracted from our reference suite of 3D and 2D
simulations $L_{\nu_{e}}$, $L_{\bar{\nu}_{e}}$, and $\dot{M}$ and plotted in
Figure 6 $L_{\nu_{e}}$ \+ $L_{\bar{\nu}_{e}}$ versus $M\dot{M}$ for these 3D
(top, left) and 2D (top, right) models. It is the $\nu_{e}$ and
$\bar{\nu}_{e}$ neutrinos that are responsible for the lion’s share of the
heating in the gain region. Since $\dot{M}$ is a monotonically-decreasing
function of time, time is increasing to the left. Scrutiny of the top panels
in Figure 6 reveals that the exploding models jump up and to the left from the
general diagonal trend, the latter defined more clearly by the trajectories
executed by the non-exploding models. This evolution is reminiscent of the
expected behavior of an exploding model upon exceeding a critical curve.
Interestingly, the thickened portions of the curves trace the phase of the
accretion by the shock of the Si/O interface. This is seen more clearly for
each of these models in the lower panels of Figure 6, which depict $M\dot{M}$
versus time after bounce. The exploding models show a marked declivity, while
for those that don’t explode this feature is more muted. The corresponding
plot versus time for $L_{\nu_{e}}$ \+ $L_{\bar{\nu}_{e}}$ is given at the top
of Figure 7.

This difference between exploding and non-exploding trajectories suggests a
common role for the accretion through the shock of the density drop of this
(or a similar) interface in often igniting a neutrino-driven explosion. Upon
encountering the shock, this rather abrupt drop in density translates into an
immediate drop in ram pressure. However, the corresponding drop in the
accumulation of mass in the deep core is delayed by a time related to
$t_{adv}$ and it is in the deeper core that the accretion component of the
luminosity is realized. The bottom of Figure 7 portrays the correlation
function between the summed luminosity and $M\dot{M}$ versus the time delay
($\Delta t$). This correlation function is calculated using
$\langle(f(t)-at-b)(g(t+\Delta t)-c(t+\Delta t)-d)\rangle$, where $a$, $b$,
$c$, $d$ are the parameters of the best linear fits for the functions $f$ and
$g$. There is a bit of noise in these curves that can compromise
interpretation, but overall the delay times are positive. Hence, the decrease
in the bounding ram pressure slightly leads the associated decrease in the
driving luminosity. During this interval, the trajectories depicted in the top
panels of Figure 6 can jump leftward999In the top panels, the trajectories
proceed from top right to bottom left.. If there is a proximate critical
curve, then such a jump can lead to explosion. Once explosion commences and
the matter starts to expand, the cooling component in the gain region
subsides, while the heating rate is more stable. Therefore, the instability
can run away. The accretion of the interface and the relative delay in the
response of the neutrino luminosity facilitate explosion. Note that in multi-
dimensional simulations the radial position of the shock wave is angle-
dependent. As a result, the accretion time of an interface is not
simultaneous. This complicates its role and this interpretation, but we
surmise not by much.

In Figure 8, we plot for the “training” 3D and 2D simulations the luminosity
and average neutrino energy evolution after bounce during the first 300
milliseconds for the $\nu_{e}$ and $\bar{\nu}_{e}$ neutrinos. These panels
show not only that the plateau luminosity is approximately monotonic with
progenitor mass (and compactness), but that the $\bar{\nu}_{e}$ luminosity
takes almost $\sim$3$-$50 milliseconds to achieve its level. Before this time
it is not substantially contributing to neutrino heating in the gain region.
This suggests that, unless $\nu_{e}$-type neutrinos alone can ignite
explosion, if the Si/O (or any other) interface is accreted too early
explosion by this route may be inhibited. Figure 8 also shows that during the
crucial early hundreds of milliseconds after bounce the average neutrino
energies are increasing. Since the neutrino-matter absorption cross sections
increase roughly as the square of the neutrino energy, due to this factor
alone the mantle power deposition grows with time and the mantle is getting
more explodable. This trend and the early lethargy in the growth of
$L_{\bar{\nu}_{e}}$ together delay explosion for many models. It is only, we
suggest, for those models with very steep outer mass density profiles and the
lowest compactness ratios that these considerations might be less important.

Importantly, these considerations are made in the context of our Fornax
simulations. Though we have incorporated a complicated suite of microphysical
processes (including inelastic scattering on electrons and nucleons); weak
magnetism and many-body corrections; velocity-dependent and multi-group
transport; and approximate general relativity (Skinner et al., 2019; Vartanyan
et al., 2019b; Burrows et al., 2020), no doubt there are limitations to our
calculations, as there are limitations in the initial models. In addition, the
outcome of the simulations will depend upon spatial and energy-group
resolution (Nagakura et al., 2019). The degree to which these limitations
undermine our findings is unknown, but we believe our simulations and the
results that have emerged from them are competitive and informative.

##  5 Test Explodability Condition

Given the earlier discussions, we now motivate a possible explodability
condition; it involves the evolution of the ram pressure at the shock and the
effect of the accretion of the Si/O (or other) interface and its associated
density drop. Exceptions to the direct predictability of this condition may be
a very low compactness parameter (leading to early explosion) or the very
early accretion of an interfacial density drop (forestalling explosion). It
may also be that the region in the space in which we suggest one look to
determine explodability possesses a fuzzier band than we suggest separating
the exploding models from the non-exploding models. This condition is derived
from the initial density profile of the progenitor alone and does not require
a detailed simulation capability. However, it is calibrated on our more
detailed Fornax results (§3;§4).

We assume that two of the main conditions for explosion are (1) a strong
density discontinuity at the Si/O interface (as long as it is not accreted
through the stalled shock too early after bounce) and/or (2) a steep density
profile (e.g., that for the 3D 9$M_{\odot}$ model). We suggest that
progenitors that meet either of these two conditions will probably explode.

In an approximate sense, the density profiles become progressively less steep
and the density jumps at the Si/O interfaces become more pronounced when
transitioning from low-mass progenitors to high-mass progenitors. Progenitors
that have both a steep density profile and a strong Si/O discontinuity (like
the 3D-11 model) will explode very easily and usually quite early. Progenitors
with a steep density profile, but without a strong discontinuity (or the
discontinuity falls in too early), also explode, but there won’t be a sudden
shock jolt caused by the discontinuity. Progenitors without a steep density
profile, but with a strong Si/O discontinuity, will explode when the interface
reaches the shock, even if the shock is already shrinking in average radius
(as with the 3D-25). But if the interface falls into the shock too early, the
luminosity and the averaged neutrino energy are still growing (particularly
the $\bar{\nu}_{e}$ luminosity, so the neutrino heating rate may not yet be
strong enough for an explosion (witness the case of 3D-15). Progenitors that
have neither a steep density profile nor a strong discontinuity seem hard to
explode. Examples of such progenitors are the 3D-13 and 3D-14.

A measure of the strength of the density discontinuity is the ram pressure
jump when the discontinuity reaches the shock. If the pressure jump is large
enough, the shock is kicked into explosion. A simple estimate of the ram
pressure at the shock radius can be obtained using the following method: We
assume that the shock radius is at $R_{s}=200$ km and that the infall velocity
is given by the free-fall velocity
$v(M)=\sqrt{\frac{2GM}{R_{s}}-\frac{2GM}{r(M)}}$. We derive an approximate
fitted relation to estimate the density just exterior to the shock
$\rho_{s}(M)=\frac{r(M)}{R_{s}}\rho(M)$, where $r(M)$ and $\rho(M)$ are the
initial radius and density of a mass shell. Then, the ram pressure is given by
$P_{\rm ram}=\rho_{s}(M)v(M)^{2}$. The infall time is estimated by a fitted
fractional free-fall time $t(M)=\sqrt{\frac{\pi}{4G\bar{\rho}(M)}}$ (Müller et
al., 2016), where $\bar{\rho}(M)=\frac{3M}{4\pi r(M)^{3}}$ is the average
density interior to a mass coordinate. The fractional change in ram pressure
is then calculated by $\frac{P_{\rm ram}(t)-P_{\rm ram}(t+\Delta t)}{P_{\rm
ram}(t)}$, where $\Delta t$ is estimated as $\frac{7R_{s}}{v(M)}$. The factor
of 7 is motivated by the shock jump for a $\gamma=4/3$ gas.

There is another expression we can use to estimate the ram pressure: $P_{\rm
ram}=\frac{v(M)}{4\pi R_{s}^{2}}\frac{dM}{dt}$. Although the ram pressure
prescriptions may differ, the fractional change across the interface is almost
the same, since we are looking at the fractional change in a very short time
period $\Delta t<10$ ms, and all slowly changing quantities (such as $v$) do
not by such changes influence the result.

Figure 9 shows the calculation of this explodability condition for the 3D
(left) and 2D (right) “training” models using the actual numbers of the
detailed simulations. Max($\frac{\Delta P_{\rm ram}}{P_{\rm ram}}$) versus
$t-t_{0}$ is plotted for each model. The straight dashed lines are those that
might be suggested to separate models that explode (upper) from those that
don’t (lower) and the circles and triangles depict actual exploding and non-
exploding models, respectively. The condition seems not bad, but not perfect,
predicting explodability better than non-explodability (particularly for the
2D models) 101010We include $t-t_{0}$ because when calculating the ram
pressure jump we search over the entire progenitor profile. Thus, all
interfaces (not only the Si/O) will be considered, and some interfaces can be
significantly more pronounced than the Si/O interface. However, those
interfaces occur either too early ($t-t_{0}<0$) or too late
($t-t_{0}>1.5\,\,{\rm seconds}$), which makes them less useful for triggering
explosion. Those that are accreted too early don’t ignite explosion because
the driving luminosities have not yet built up. Those accreted too late may
not lead to explosion because the driving luminosities have already begun to
abate..

However, in Figure 10 we make the same plot, but use the analytic
approximations described above for the same quantities. Figures 9 and 10 look
much the same. Given this, we then ask the question. If we now calculate using
Fornax with its full physics capabilities 100 2D models and then compare the
outcomes (exploding or not) using the analytic approach to “$\frac{\Delta
P_{\rm ram}}{P_{\rm ram}}$” versus time after bounce employing only the
initial model mass density profiles how does this explosion condition fare?

We recapitulate and summarize here our simple prescription for determining the
explodability of a model given only the initial density profile:

  1. 1.

Assume that the shock radius is $R_{s}=200$ km which is a reasonable estimate
for most progenitors. This assumed radius is too large for interfaces that
fall in very early, so we ignore matter interior to the Si/O
interface111111The Si/O interface in this paper is defined as the density
discontinuity closest to the inner boundary of the oxygen shell in which the
oxygen abundance is above 15%. Sometimes the Si core size given by this
definition is underestimated, but this definition does not compromise the
accuracy of our criterion.. More specifically, we ignore data with mass
coordinate $M<M_{\rm Si}-0.05M_{\odot}$.

  2. 2.

Infall velocity and infall time are given by the free-fall velocity and the
fitted fractional free-fall time:
$v(M)=\sqrt{\frac{2GM}{R_{s}}-\frac{2GM}{r(M)}}$ and
$t(M)=\sqrt{\frac{\pi}{4G\bar{\rho}(M)}}$, where $r(M)$ is the initial radius
of the mass shell and $\bar{\rho}(M)=\frac{3M}{4\pi r(M)^{3}}$ is the average
density interior to it.

  3. 3.

The density just exterior to the shock is given by an approximate fitted
relation $\rho_{s}(M)=\frac{r(M)}{R_{s}}\rho(M)$, where $\rho(M)$ is the
initial density of the mass shell. This initial density is extracted when the
core collapse speed exceeds a maximum of 900 km s-1 (Sukhbold et al., 2018).

  4. 4.

The ram pressure is assumed given by $P_{\rm ram}=\rho_{s}(M)v(M)^{2}$.

  5. 5.

The time delay is estimated as $\Delta t=\frac{7R_{s}}{v(M)}$.

  6. 6.

The fractional change in ram pressure is then calculated by $\frac{\Delta
P_{\rm ram}}{P_{\rm ram}}=\frac{P_{\rm ram}(t)-P_{\rm ram}(t+\Delta t)}{P_{\rm
ram}(t)}$. (The infall time $t$ is a monotonically increasing function of mass
coordinate $M$, so we can use its inverse function to replace all $M$
dependencies with $t$ dependencies.)

  7. 7.

Based on our simulation results, we make the following prediction: the
progenitor will explode if Max$(\frac{\Delta P_{\rm ram}}{P_{\rm ram}})>0.28$,
and vise versa.

  8. 8.

The bounce time is estimated using a fitted function
$t_{0}=0.218\left(\frac{\rho_{c}}{10^{10}\text{g cm}^{-3}}\right)^{-0.354}$ s,
where $\rho_{c}$ is the initial central mass density. This is not a very good
fit, but it gives us a sense on whether an interface reaches the shock too
early or too late.

##  6 The 2D Model Suite

To build this model suite, we employ the SFHo nuclear equation of state
(Steiner et al., 2013) and the same simulation code, Fornax, we have employed
in previous recent investigations. The numerical and physical details
incorporated into Fornax can be found in previous papers (Vartanyan et al.,
2019a; Burrows et al., 2019; Skinner et al., 2019; Burrows et al., 2020;
Burrows & Vartanyan, 2021)

The density and specific entropy profiles of the 100 initial progenitor models
taken from (Sukhbold et al., 2016) and Sukhbold et al. (2018) are given in
Figures 11 and 12. This collection of models clearly spans a wide range of
model structures, encompassing 9 to 27 M⊙. However, we do not in this study go
beyond 27 M⊙, leaving that range and (perhaps) black-hole formation to another
study.

Figure 13 depicts the a priori expectations for the complete set of models in
Sukhbold et al. (2018); we use only a fraction of those for our 100-model
study. The color codes for the ZAMS mass, indicating where in progenitor-mass
space we might expect duds.

The 100 2D models were simulated out to about 1 seconds after bounce or until
explosion, with outer radii from 30,000 to 70,000 kilometers (km). The default
number of radial zones was 1024 and of angular ($\theta$) zones was 128.
However, we redid numerous models using 256 angular zones. This was guided by
the results in Nagakura et al. (2019), who suggested models were more
explodable at higher resolution. For a subset of these models this was indeed
what we found. In fact, including the higher resolution cohort increased the
fidelity of the “Max$\frac{\Delta P}{P}$” test.

The results for our 100-model cohort using only initial profile data, but
determining explodability using the detailed 2D simulations (“circles” or
“triangles”), is given in Figure 14. This plot is to be compared with Figures
10 and 13. The number of false positives is five in $\sim$40, with six false
negatives. However, the percentage success rate is high. Moreover, the false
positives are all close to the putative demarcation line, suggesting a
slightly more nuanced definition may be called for. We note that the
preponderance of the failed models are in the low to intermediate ZAMS mass
range and that all the more massive models in the upper progenitor mass range
of this study explode. This is in stark contrast with the expectation of the
naive compactness condition. Figure 15 makes this point more succinctly.

We emphasize that detailed stellar evolution calculations (Sukhbold et al.,
2016; Sukhbold et al., 2018) find that the initial structures of collapse
depend in a “chaotic” way on, among other things, shell mergers and the
transition from radiative to convective core carbon burning (near $\sim$18
M⊙?). The different branches that emerge are seen clearly in Figure 15, where
almost all the failed explosions reside on one branch. The larger significance
of this clustering remains to be determined, but it seems to coincide with the
model space in the Sukhbold et al. (2018) collection that experiences two
simultaneous oxygen burning shells. In their model suite, above the associated
ZAMS mass range ($\geq$15 M⊙) one of these shells diminishes or disappears,
leaving the other to dominate. This seems to result in a more substantial
density jump at the “residual” Si/O interface for those more massive
progenitors. A similar branch dependence is seen in Figure 16, which renders
the initial central mass density versus ZAMS mass. This plot too suggests an
interesting relationship between branches in density structure that result
from the consequences of the chaotic complications that attend core and shell
burning of massive stars during their terminal stages.

Finally, in Figure 17 we test the Ertl et al. (2016) condition against our
100-model suite. We see that though it finds most of the non-exploding models,
it fails to find many of the exploding ones.

##  7 Conclusions

Calibrating with detailed 2D core-collapse supernova simulations, we derived a
simple explosion condition based solely on the density profiles of the
terminal cores provided by state-of-the-art stellar evolution calculations of
progenitor massive stars. This condition captures the vast majority of the
behavior of the 100 2D models we performed subsequently to gauge its
usefulness. The goal was to be able to predict, without resort to detailed
simulation, the explodability of a given massive star. Under the assumption
that our 2D multi-group radiation/hydrodynamic simulations using the code
Fornax yield reliable results, we find that the simple maximum fractional ram
pressure jump discriminant we defined works well ($\sim 90$% of the time). We
emphasize that this method employs simple analytic prescriptions for the
various inputs needed that are based solely upon a given initial mass density
profile inherited from detailed stellar model calculations. The maximum ram
pressure jump generally occurs at the time of accretion of the Si/O interface,
but not always. We note that the false positives are all close to the
horizontal demarcation line (see Figure 14), suggesting that this line can be
altered in shape only slightly to include them or that there may be other
“second-order” criteria that come into play. These could be the initial
thermal and/or Ye profiles. The latter have a slight influence over the timing
of accretion not factored into our crude analytic model for the mass accretion
rate evolution through the shock, incorporating as it does only the initial
density profile.

There seem to be two types of false negatives. The first is the models on
Figure 14 with $\frac{\Delta P_{\rm ram}}{P_{\rm ram}}\sim 0.2$. These models
have masses below $\sim$12.5 M⊙ and have steep initial density profiles and
low compactness. Such models don’t need a density discontinuity (such as at a
pronounced Si/O interface) to explode. The second class of false negatives
involves models that do not have steep density profiles and explode later
after bounce. They are much closer to the horizontal line. Though in drawing
our sample of progenitors we were initially guided by the Salpeter initial
mass function, we can not claim to have done so rigorously. Hence, the
fraction of models we see explode in this investigation should not be confused
with the fraction of massive stars in Nature above $\sim$8 M⊙ that do explode.
This fraction, and the fraction that result in black holes, has yet to be
determined. However, should our explodability condition have some validity,
this fraction can easily be calculated in the context of a definitive suite of
progenitor models spanning the requisite range of massive-star masses.

It should be mentioned that we drew our initial models from both the Sukhbold
et al. (2016) and Sukhbold et al. (2018) compilations. All the models with
$\leq$12 M⊙ came from the former, which sampled the associated mass range far
more sparsely than did Sukhbold et al. (2018) for masses above 12 M⊙. This is
reflected in the greater number of studied models above 12 M⊙. Also, we did
not include progenitors with ZAMS masses above 27 M⊙. We suspect that many
models in this higher mass range do not explode, or may explode initially, but
experience fallback later upon encountering an outer mantle with significant
binding energy that the blast can not fully overcome. The result would be
late-time fallback, but there might still be some mass ejection to infinity
with a reduced asymptotic explosion energy (Chan et al., 2018). This
intriguing possibility remains to be demonstrated, but would suggest that even
black hole formation might be associated with some sort of (supernova)
explosion.

There are many caveats that deserve mention. Our results are conditioned on
the fidelity with which our current implementation of Fornax adheres to
Nature. There remain issues concerning the neutrino-matter interaction, the
nuclear equation of state, the possible effects of neutrino oscillations, and
the accuracy of the numerical algorithms employed. Furthermore, we have noted
in the past that the spatial resolution of the computational grid can be a
factor in the outcome, with higher resolution models exploding when lower
resolution models don’t (Nagakura et al., 2019). We have endeavored in this
investigation to employ adequate resolution, but this remains to be definitely
demonstrated. Moreover, we are presuming that models that explode in 2D also
explode in 3D, and vice versa. This has generally been our experience. We note
that in 3D the horizontal explosion/no-explosion demarcation line may need to
be slightly adjusted (see Figures 9 and 10). However, determining whether this
is necessary must await a large study using a corresponding number of detailed
3D simulations. In addition, we have not addressed the effect of the initial
seed perturbations (Couch & Ott, 2013; Müller et al., 2017a; Burrows et al.,
2018) in kick-starting the turbulence behind the stalled shock, nor the fact
that chaotic flow will by its nature lead to stochastic outcomes 121212though
whether only in detail, or qualitatively, has not been determined. It is also
the case that even state-of-the-art progenitor modeling is still in flux; the
mapping between progenitor structure and ZAMS mass, for instance, has not
converged. Finally, we have not addressed in this study the possible roles of
the equation of state, rotation or magnetic fields (Burrows et al., 2007; Suwa
et al., 2013; Couch, 2013; Mösta et al., 2014; Yasin et al., 2020; Kuroda et
al., 2020; Obergaulinger & Aloy, 2020; Aloy & Obergaulinger, 2021;
Obergaulinger & Aloy, 2021; Kuroda, 2021).

Nevertheless, the straightforward explodability condition we have obtained has
been calibrated with a large collection of state-of-the-art 2D simulations and
is simple to implement. It demonstrates fidelity to our more detailed
simulations $\sim$90% of the time and shows promise that it might be further
generalized while still employing data from only the unstable Chandrasekhar
progenitors core-collapse modelers inherit. We suggest it is a more credible
and robust predictor than can currently be found in the literature, but plan
to refine it in numerous ways in the future.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## Acknowledgements

We thank Chris White, David Radice, and Hiroki Nagakura for insights and
advice during the germination of this project. We acknowledge support from the
U. S. Department of Energy Office of Science and the Office of Advanced
Scientific Computing Research via the Scientific Discovery through Advanced
Computing (SciDAC4) program and Grant DE-SC0018297 (subaward 00009650) and
support from the U. S. National Science Foundation (NSF) under Grants
AST-1714267 and PHY-1804048 (the latter via the Max-Planck/Princeton Center
(MPPC) for Plasma Physics). The three-dimensional simulations were performed
on Blue Waters under the sustained-petascale computing project, which was
supported by the National Science Foundation (awards OCI-0725070 and
ACI-1238993) and the state of Illinois. Blue Waters was a joint effort of the
University of Illinois at Urbana–Champaign and its National Center for
Supercomputing Applications. We also acknowledge access to the Frontera
cluster (under awards AST20020 and AST21003), and this research is part of the
Frontera computing project at the Texas Advanced Computing Center (Stanzione
et al., 2020). Frontera is made possible by NSF award OAC-1818253.
Additionally, a generous award of computer time was provided by the INCITE
program, enabling this research to use resources of the Argonne Leadership
Computing Facility, a DOE Office of Science User Facility supported under
Contract DE-AC02-06CH11357. Finally, the authors acknowledge computational
resources provided by the high-performance computer center at Princeton
University, which is jointly supported by the Princeton Institute for
Computational Science and Engineering (PICSciE) and the Princeton University
Office of Information Technology, and our continuing allocation at the
National Energy Research Scientific Computing Center (NERSC), which is
supported by the Office of Science of the U. S. Department of Energy under
contract DE-AC03-76SF00098.

Figure 1: Evolution of the maximum shock radius. Left: 3D simulations from
Burrows et al. (2020). Right: 2D simulations from Burrows & Vartanyan (2021).
The numbers after the “3D” or “2D” indicate the ZAMS mass of the progenitor in
units of solar masses (M⊙). Solid lines are exploding models, while dashed
lines are non-exploding ones. The dot on each line marks the time when the
maximum shock position encounters the Si/O interface. From the 3D simulations,
we see that the maximum shock radius in each model evolves in a very similar
fashion until the Si/O interface hits the shock. Despite experiencing stronger
oscillations, the 2D simulations evince the same behavior. The dashed
horizontal lines denote 250 km on left panels (re 3D) and 400 km on right
panels (re 2D).

Figure 2: The $\rho$ (in gm cm-3) versus interior mass ($M$, in M⊙) profiles
of the progenitors. Left: 3D simulations published in Burrows et al. (2020).
All these progenitors are taken from Sukhbold et al. (2016), except for the 25
$M_{\odot}$ model, which is taken from Sukhbold et al. (2018). Right: 2D
simulations published in Burrows & Vartanyan (2021). The solid lines are
exploding models, while the dashed lines are non-exploding models. The dot on
each line marks the position of the Si/O interface defined via an abundance
transition, which may slightly deviate from the start point of the associated
density jump. Most non-exploding models generally have a weak density jump at
the Si/O interface.

Figure 3: Residual baryonic PNS mass in M⊙ (top) and explosion energy in
Bethes (bottom) as a function of progenitor compactness (at 1.75 M⊙) (left)
and progenitor mass (right) from the set of 2D late-time simulations published
in Burrows & Vartanyan (2021) (“BV 21”) and using the Sukhbold et al. (2018)
progenitors. Models 12- and 15-M⊙ do not explode (marked with “X”s). We
identify a strong trend with compactness of PNS mass (in particular) and
explosion energy. There is a corresponding trend with progenitor ZAMS mass,
but it is weaker (Burrows & Vartanyan, 2021; O’Connor & Ott, 2013). See text
for a discussion. Figure 4: The explosion energy (in Bethes) versus the
initial binding energy of the progenitor exterior to the residual PNS baryon
mass (in Bethes) for the 2D models published in Burrows & Vartanyan (2021)
(“BV 21”), colored by compactness. Each model is identified by its ZAMS mass
and colored by its compactness. Non-exploding models are indicated with “X”s.
A clear trend, with a large spread, is in evidence.

Figure 5: The exterior binding energy of the massive-star progenitors in
Bethes (1051 ergs) versus interior mass (in M⊙). Left: 3D simulations; Right:
2D simulations. The solid lines are exploding models, while the dashed lines
are the non-exploding models. The circular dot on each line shows the position
of the Si/O interface, while the square shows the PNS mass at the end of each
simulation.

Figure 6: The sum of the electron and anti-electron neutrino luminosities (at
100 km) versus $M\dot{M}$ (top) and $M\dot{M}$ versus t (time after bounce)
(bottom). The left plots are for the 3D models and the right plots are for the
2D models. The accretion rates ($\dot{M}$) and interior masses ($M$) are
measured at 250 km (for 3D) and 400 km (for 2D), i.e., at the radii shown by
the dashed horizontal curves in Figure 1. As before, the 3D simulations are
from Burrows et al. (2020), while the 2D simulations are from Burrows &
Vartanyan (2021). Solid lines are exploding models, while dashed lines are
non-exploding models. The thick portion of each curve marks the density
discontinuity (and, hence, the discontinuity in the accretion rate) at the
Si/O interface, which is more clearly shown in the bottom panels. In these
plots, we show only the phase before the maximum shock radius reaches 250 km
or 400 km so that the accretion rate shown is not influenced by the explosion.
From the top panels, we see that the exploding and non-exploding models behave
differently. There is a roughly diagonal relation that all simulations (both
2D and 3D) adhere to at later times, and only the exploding models deviate (to
the left on the top plots) from that diagonal relation. See text for a
discussion.

Figure 7: Top panels show the luminosity and accretion rate (times $M$) as a
function of time for the 3D (left) and 2D (right) model sets. The solid lines
are the sum of the $\nu_{e}$ and $\bar{\nu}_{e}$ luminosities (left ordinate),
while the dashed lines are $M\dot{M}$ (right ordinate). The bottom curves
render the luminosity-$M\dot{M}$ correlation functions. We use the thick parts
of the top curves to calculate this correlation. Although the uncertainty of
this time delay measurement is largely due to noisy luminosities (particularly
in 2D), we still see that the luminosity lags the accretion rate near when the
Si/O interface (or any significant interface) is accreted. Such a time delay
causes excursions from the diagonal relation in Figure 6 that are often
germane to the onset of explosion. Figure 8: Electron neutrino (left) and
anti-electron (right) neutrino luminosity and average neutrino energy versus
time after bounce for the reference 3D simulations. The top panels are
luminosities measured at 100 km and the bottom panels are the average neutrino
energies. The solid lines trace the exploding models, while the dashed lines
trace the non-exploding models. The dot on each line marks the time when the
Si/O interface is accreted through the shock. We see that the average electron
neutrino energy first decreases, and then increases, with time after bounce.
For the anti-electron-types, the average energy only increases during the
first many hundreds of milliseconds after bounce. Since the cross sections are
proportional to the square of the neutrino energy, the neutrino heating rates
track these behaviors for a given luminosity. Hence, if the Si/O interface
meets the shock too early, the heating rate could be too low to support an
explosion. See the text for a discussion.

Figure 9: Maximum fractional ram pressure change over a time shift $\Delta t$
measured for the 3D (left) and 2D (right) simulations. The circular dots
identify exploding models, while the triangles indicate the non-exploding
models. $\Delta t$ is taken to be equal to $\frac{7R_{s,max}}{v}$, where
$R_{s,max}$ is the maximum shock radius and $v$ is the infall velocity
measured at 250 km (for 3D) and 400 km (for 2D). For this plot, we ignore the
pressure change caused by the ram pressure change at bounce. The maximum
fractional change generally occurs at the Si/O interface, and the amplitude of
this fractional change shows how strongly the discontinuity influences the
explosion. From this plot, we can see that the ram pressure changes in non-
exploding models are smaller, while for most of the exploding ones they are
larger. We draw a horizontal dashed line at 0.28 to distinguish the strong-
discontinuity and weak-discontinuity clusters which can be more clearly seen
in Figure 10. Most exploding models satisfy this criterion, with three
exceptions (3D-15, 2D-13 and 2D-14). However, when we look at the shock
evolution of these three exceptions in Figure 1, we find that the shocks of
2D-13 and 2D-14 expand much more slowly compared to the other exploding models
at the time the Si/O interface falls in. This indicates that either they are
marginally explosive or that the explosion is not directly triggered by the
ram pressure jump. The radius of the shock of model 3D-15 expands earlier and
faster compared to other non-exploding models, which indicates that it might
be marginally explosive or that the interface falls in too early to have
enough neutrino heating rate for an explosion. See text for a discussion.

Figure 10: The same as Figure 9, but estimated from the initial progenitor
models themselves without reference to the simulations. Circular dots are
exploding models, while the triangles are the non-exploding models. $\Delta t$
is given by $\frac{7R_{s}}{v(M)}$, where $R_{s}=200$ km. The factor of 7 comes
simply from the density jump at a shock for a $\gamma=4/3$ gas. The infall
velocity is estimated using the free-fall velocity
$v(M)=\sqrt{\frac{2GM}{R_{s}}-\frac{2GM}{r(M)}}$. The infall time is
calculated using a fraction of this free-fall time fit to the simulations:
$t(M)=\sqrt{\frac{\pi}{4G\bar{\rho}(M)}}$, where $\bar{\rho}(M)=\frac{3M}{4\pi
r(M)^{3}}$ is the average density interior to a certain mass coordinate. The
bounce time is estimated with the following approximate fitted function:
$t_{0}=0.218\left(\frac{\rho_{c}}{10^{10}\text{g cm}^{-3}}\right)^{-0.354}$ s.
Comparing these plots to Figure 9, we see that this approach predicts
explodability quite well. The general trend persists: exploding models tend to
have stronger maximum ram pressure jumps, while the non-exploding ones have
weaker maximum ram pressure jumps.  Figure 11: The mass density (in gm cm-3)
versus interior mass (in M⊙) for the 100 progenitors explored in our large 2D
study. The models were selected to span the distribution in density space,
including both steep and shallow profiles as well strong and weak
compositional interfaces. Figure 12: The entropy profiles that correspond to
the models depicted in Figure 11. Figure 13: Maximum fractional ram pressure
change given a time shift $\Delta t$, estimated from all the progenitor models
we downloaded as a prelude to the large 2D explosion study. Progenitors with
$M<12M_{\odot}$ are from Sukhbold et al. (2016), while progenitors above 12 M⊙
are from Sukhbold et al. (2018). We see that most progenitors with $M>$ 16 M⊙
are above the approximate explosion line, while many low-mass progenitors are
below it. However, we would expect that some low-mass progenitors below the
line will still explode because of their steep density profiles. Figure 14:
Same as Figure 13, but for the 100 models we actually simulated in 2D for this
study. The dashed line shows the simple explosion line from the right panel of
Figure 10. Circular dots are the exploding models, while the triangles are
non-exploding ones. Among 100 simulations, there are only $\sim$10 exceptions
to our pressure jump criterion. Five of them do not explode, but are above the
line (false positives), and six of them explode, but are below the line (false
negatives). Three false positive dots and two false negative dots are all very
close to the criterion line, which means that 1) there might be a marginally
exploding region, 2) the demarcation line can be slightly shifted and/or
needn’t be straight, or 3) there are secondary parameters of relevance when
things are marginal. Other false negative dots can be explained by steep
density profiles, in which the abrupt density jump is not strong, but the
overall steepness of the density profile triggers the explosion. See text for
a discussion. Figure 15: Compactness at 1.75 M⊙ versus oxygen core mass. The
black dots are a comprehensive set of progenitor models found in the Sukhbold
et al. (2016) and Sukhbold et al. (2018). Blue and red triangles are the 100
2D simulations in this work. Blue triangles are the exploding models, while
red triangles are the non-exploding models. From this figure we see that most
of the non-exploding models reside along a specific progenitor branch. We note
that in Sukhbold et al. (2018) the carbon-oxygen core mass is proportional to
the ZAMS mass. This plot is therefore similar to their Figure 15, but
contradicts its conclusions concerning which models explode. Figure 16:
Initial central density versus ZAMS mass. The black dots are a comprehensive
set of progenitor models found in the Sukhbold et al. (2016) and Sukhbold et
al. (2018). Blue and red triangles are the 100 2D simulations in this work.
Blue triangles are the exploding ones, while red triangles are non-exploding
models. Most non-exploding models have a similar central density around
$6\times 10^{9}$ g cm-3 Figure 17: The Ertl condition plot. $M_{4}$ is the
mass coordinate where the specific entropy first exceeds 4$k_{b}/$baryon,
which roughly marks the position of the Si/O interface, and
$\mu_{4}=\frac{0.3M_{\odot}\times
1000\text{km}}{R(M_{4}+0.3M_{\odot})-R(M_{4})}$. The black dots are a
comprehensive set of progenitor models found in the Sukhbold et al. (2016) and
Sukhbold et al. (2018). Blue and red triangles are the 100 2D simulations in
this work. Blue triangles are the exploding ones, while red triangles are non-
exploding models. The dashed line is the w18.0 condition in Ertl et al.
(2016). Almost all non-exploding models are below the Ertl et al. line, but
the exploding models are both above and below the line.

## References

  * Aloy & Obergaulinger (2021) Aloy M. Á., Obergaulinger M., 2021, MNRAS, 500, 4365
  * Arnett (1966) Arnett W. D., 1966, Canadian Journal of Physics, 44, 2553
  * Bethe & Wilson (1985) Bethe H. A., Wilson J. R., 1985, ApJ, 295, 14
  * Bollig et al. (2021) Bollig R., Yadav N., Kresse D., Janka H.-T., Müller B., Heger A., 2021, ApJ, 915, 28
  * Bruenn & Dineva (1996) Bruenn S. W., Dineva T., 1996, ApJ, 458, L71
  * Bruenn et al. (2016) Bruenn S. W., et al., 2016, ApJ, 818, 123
  * Burrows (1987) Burrows A., 1987, ApJ, 318, L57
  * Burrows (2013) Burrows A., 2013, Reviews of Modern Physics, 85, 245
  * Burrows & Goshy (1993) Burrows A., Goshy J., 1993, ApJ, 416, L75
  * Burrows & Vartanyan (2021) Burrows A., Vartanyan D., 2021, Nature, 589, 29
  * Burrows et al. (1995) Burrows A., Hayes J., Fryxell B. A., 1995, ApJ, 450, 830
  * Burrows et al. (2007) Burrows A., Dessart L., Livne E., Ott C. D., Murphy J., 2007, ApJ, 664, 416
  * Burrows et al. (2018) Burrows A., Vartanyan D., Dolence J. C., Skinner M. A., Radice D., 2018, Space Science Reviews, 214, 33
  * Burrows et al. (2019) Burrows A., Radice D., Vartanyan D., 2019, MNRAS, 485, 3153
  * Burrows et al. (2020) Burrows A., Radice D., Vartanyan D., Nagakura H., Skinner M. A., Dolence J. C., 2020, MNRAS, 491, 2715
  * Chan et al. (2018) Chan C., Müller B., Heger A., Pakmor R., Springel V., 2018, ApJ, 852, L19
  * Chatzopoulos et al. (2016) Chatzopoulos E., Couch S. M., Arnett W. D., Timmes F. X., 2016, The Astrophysical Journal, 822, 61
  * Colgate & White (1966) Colgate S. A., White R. H., 1966, ApJ, 143, 626
  * Couch (2013) Couch S. M., 2013, ApJ, 765, 29
  * Couch & Ott (2013) Couch S. M., Ott C. D., 2013, ApJ, 778, L7
  * Dessart et al. (2006) Dessart L., Burrows A., Livne E., Ott C. D., 2006, ApJ, 645, 534
  * Dolence et al. (2013) Dolence J. C., Burrows A., Murphy J. W., Nordhaus J., 2013, ApJ, 765, 110
  * Ertl et al. (2016) Ertl T., Janka H.-T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124
  * Fields & Couch (2020) Fields C. E., Couch S. M., 2020, ApJ, 901, 33
  * Fields & Couch (2021) Fields C. E., Couch S. M., 2021, arXiv e-prints, p. arXiv:2107.04617
  * Foglizzo et al. (2007) Foglizzo T., Galletti P., Scheck L., Janka H.-T., 2007, ApJ, 654, 1006
  * Glas et al. (2019) Glas R., Just O., Janka H. T., Obergaulinger M., 2019, ApJ, 873, 45
  * Gogilashvili & Murphy (2021) Gogilashvili M., Murphy J. W., 2021, arXiv e-prints, p. arXiv:2110.10173
  * Herant et al. (1994) Herant M., Benz W., Hix W. R., Fryer C. L., Colgate S. A., 1994, ApJ, 435, 339
  * Janka (2012) Janka H.-T., 2012, Annual Review of Nuclear and Particle Science, 62, 407
  * Janka et al. (2016) Janka H.-T., Melson T., Summa A., 2016, Annual Review of Nuclear and Particle Science, 66, 341
  * Kuroda (2021) Kuroda T., 2021, ApJ, 906, 128
  * Kuroda et al. (2020) Kuroda T., Arcones A., Takiwaki T., Kotake K., 2020, ApJ, 896, 102
  * Lentz et al. (2015) Lentz E. J., et al., 2015, ApJ, 807, L31
  * Mösta et al. (2014) Mösta P., et al., 2014, The Astrophysical Journal, 785, L29
  * Müller (2016) Müller B., 2016, Publ. Astron. Soc. Australia, 33, e048
  * Müller et al. (2016) Müller B., Heger A., Liptai D., Cameron J. B., 2016, MNRAS, 460, 742
  * Müller et al. (2017a) Müller B., Melson T., Heger A., Janka H.-T., 2017a, MNRAS, 472, 491
  * Müller et al. (2017b) Müller T., Prieto J. L., Pejcha O., Clocchiatti A., 2017b, ApJ, 841, 127
  * Müller et al. (2019) Müller B., et al., 2019, Monthly Notices of the Royal Astronomical Society, 484, 3307–3324
  * Murphy & Burrows (2008) Murphy J. W., Burrows A., 2008, ApJ, 688, 1159
  * Murphy & Dolence (2017) Murphy J. W., Dolence J. C., 2017, ApJ, 834, 183
  * Murphy et al. (2013) Murphy J. W., Dolence J. C., Burrows A., 2013, ApJ, 771, 52
  * Nagakura et al. (2019) Nagakura H., Burrows A., Radice D., Vartanyan D., 2019, MNRAS, 490, 4622
  * Nagakura et al. (2020) Nagakura H., Burrows A., Radice D., Vartanyan D., 2020, MNRAS, 492, 5764
  * O’Connor & Couch (2018) O’Connor E. P., Couch S. M., 2018, ApJ, 865, 81
  * O’Connor & Ott (2011) O’Connor E., Ott C. D., 2011, ApJ, 730, 70
  * O’Connor & Ott (2013) O’Connor E., Ott C. D., 2013, ApJ, 762, 126
  * Obergaulinger & Aloy (2020) Obergaulinger M., Aloy M. Á., 2020, MNRAS, 492, 4613
  * Obergaulinger & Aloy (2021) Obergaulinger M., Aloy M. A., 2021, Monthly Notices of the Royal Astronomical Society, 503, 4942–4963
  * Ott et al. (2018) Ott C. D., Roberts L. F., da Silva Schneider A., Fedrow J. M., Haas R., Schnetter E., 2018, ApJ, 855, L3
  * Parker (1958) Parker E. N., 1958, ApJ, 128, 664
  * Pejcha & Thompson (2012) Pejcha O., Thompson T. A., 2012, ApJ, 746, 106
  * Raives et al. (2018) Raives M. J., Couch S. M., Greco J. P., Pejcha O., Thompson T. A., 2018, MNRAS, 481, 3293
  * Raives et al. (2021) Raives M. J., Thompson T. A., Couch S. M., 2021, MNRAS, 502, 4125
  * Roberts et al. (2016) Roberts L. F., Ott C. D., Haas R., O’Connor E. P., Diener P., Schnetter E., 2016, ApJ, 831, 98
  * Skinner et al. (2019) Skinner M. A., Dolence J. C., Burrows A., Radice D., Vartanyan D., 2019, ApJS, 241, 7
  * Stanzione et al. (2020) Stanzione D., West J., Evans R. T., Minyard T., Ghattas O., Panda D. K., 2020, in PEARC ’20. Practice and Experience in Advanced Research Computing. Portland, OR, pp 106–111
  * Steiner et al. (2013) Steiner A. W., Hempel M., Fischer T., 2013, ApJ, 774, 17
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H.-T., 2016, ApJ, 821, 38
  * Sukhbold et al. (2018) Sukhbold T., Woosley S. E., Heger A., 2018, ApJ, 860, 93
  * Summa et al. (2016) Summa A., Hanke F., Janka H.-T., Melson T., Marek A., Müller B., 2016, ApJ, 825, 6
  * Summa et al. (2018) Summa A., Janka H.-T., Melson T., Marek A., 2018, ApJ, 852, 28
  * Suwa et al. (2013) Suwa Y., Takiwaki T., Kotake K., Fischer T., Liebendörfer M., Sato K., 2013, ApJ, 764, 99
  * Thompson (2000) Thompson C., 2000, The Astrophysical Journal, 534, 915
  * Vartanyan et al. (2018) Vartanyan D., Burrows A., Radice D., Skinner M. A., Dolence J., 2018, MNRAS, 477, 3091
  * Vartanyan et al. (2019a) Vartanyan D., Burrows A., Radice D., Skinner M. A., Dolence J., 2019a, MNRAS, 482, 351
  * Vartanyan et al. (2019b) Vartanyan D., Burrows A., Radice D., 2019b, MNRAS, 489, 2227
  * Vartanyan et al. (2021) Vartanyan D., Laplace E., Renzo M., Götberg Y., Burrows A., de Mink S. E., 2021, ApJ, 916, L5
  * Vink (2011) Vink J. S., 2011, Ap&SS, 336, 163
  * Wilson (1985) Wilson J. R., 1985, in Centrella J. M., Leblanc J. M., Bowers R. L., eds, Numerical Astrophysics. p. 422
  * Yasin et al. (2020) Yasin H., Schäfer S., Arcones A., Schwenk A., 2020, Phys. Rev. Lett., 124, 092701