Crucial Physical Dependencies of the Core-Collapse Supernova Mechanism

We explore with self-consistent 2D F{\sc{ornax}} simulations the dependence of the outcome of collapse on many-body corrections to neutrino-nucleon cross sections, the nucleon-nucleon bremsstrahlung rate, electron capture on heavy nuclei, pre-collapse seed perturbations, and inelastic neutrino-electron and neutrino-nucleon scattering. Importantly, proximity to criticality amplifies the role of even small changes in the neutrino-matter couplings, and such changes can together add to produce outsized effects. When close to the critical condition the cumulative result of a few small effects (including seeds) that individually have only modest consequence can convert an anemic into a robust explosion, or even a dud into a blast. Such sensitivity is not seen in one dimension and may explain the apparent heterogeneity in the outcomes of detailed simulations performed internationally. A natural conclusion is that the different groups collectively are closer to a realistic understanding of the mechanism of core-collapse supernovae than might have seemed apparent.


Introduction
A goal of core-collapse supernova theory is to explain the mechanism of explosion. It is an accepted truism of the field that a necessary condition for this is that complicated multi-dimensional simulation codes incorporating the requisite neutrino, nuclear, and gravitational physics reproduce such explosions robustly, yielding at the very least the requisite asymptotic energies, residual neutron star masses, and nucleosynthesis. A simple analytic explanation has not been forthcoming and, given the manifest complexities of the process, would not be deemed credible. It is thought that one litmus test of success would be such a demonstration for a non-rotating model between ∼8 and ∼20 M ⊙ , the progenitor ZAMS (Zero-Age Main-Sequence) mass regime that must provide the lion's share of core-collapse supernova events.
However, to date the various groups engaged in such efforts around the world have failed to agree on the outcomes of the collapse of otherwise similar progenitor massive star cores. This is despite claims to have embedded the necessary physics and microphysics into the simulations. In fact, the numerical algorithms, resolutions, and input physics all differ, and even when the physics is deemed similar, the implementations and approximations surely differ. The ORNL group (Bruenn et al. 2013(Bruenn et al. ,2016, with their CHIMERA code, uses multi-group flux-limited diffusion (MGFLD) neutrino transport (Bruenn 1985), the VH-1 Newtonian hydrodynamics package, a monopole correction for general relativity (GR) (Marek et al. 2006), but multiple one-dimensional solves for multi-dimensional transport using the so-called "ray-by-ray+" approach (Buras et al. 2003;Burrows, Hayes, & Fryxell 1995). Such a dimensional reduction for the transport, particularly manifest in 2D, ignores lateral, non-radial radiative transport, which has been shown to be of quantitative (Ott et al. 2008;Brandt et al. 2011;Burrows 2013;Dolence, Burrows, & Zhang 2015;Sumiyoshi et al. 2015) and qualitative (Skinner, Burrows, & Dolence 2016) importance. To the point, Skinner et al. (2016) have shown that the ray-by-ray anomalies in the angular distribution of the radiation field and corresponding neutrino heating rates can reinforce the axial sloshing motions in 2D and push the shock into explosion 1 . Bruenn et al. (2013Bruenn et al. ( ,2016 find explosions in 2D (axially symmetry) for all the progenitors they studied (the 12-, 15-, 20-, and 25-M ⊙ models of Woosley & Heger 2007), but the models all explode at about the same post-bounce time (∼100 milliseconds) and the shock radii never decrease in value. When performed in 3D for the 15-M ⊙ model of Woosley & Heger (2007), Lentz et al. (2015) obtain a weaker explosion delayed in its onset by an extra ∼100 milliseconds. Müller et al. (2012ab), using the conformally-flat CoCoNuT hydrodynamics code in combination with the VERTEX transport solver and the ray-by-ray+ approach, find that models explode in 2D, but explode with lower energies and at post-bounce explosion times that differ significantly from those found by the ORNL group. In addition, their shock radii generally decrease from a peak occurring near a time ∼200 ms after bounce, before increasing just prior to explosion hundreds of milliseconds later (and all at different post-bounce times). Using the VERTEX-PROMETHEUS code, the Garching group obtained weak explosions in 2D for a 11.2-M ⊙ progenitor ) and for a rotating 15-M ⊙ progenitor (Marek & Janka 2009), though the outcomes when proper correction is made for the corresponding mantle binding energies were not clear. More recently, the Garching group ) obtained explosions in 2D for a broad range of Woosley & Heger (2007) progenitors from 11 to 28 M ⊙ , but these calculations as well did not comport with those of the ORNL group in the timescales and energies seen. The Garching group has also used for all their multi-D VERTEX supernova simulations a variant of the problematic ray-by-ray+ dimensionallyreduced transport approach.
However,  and Couch & Ott (2015), and originally Burrows, Hayes, & Fryxell (1995), highlight the importance of turbulent pressure behind the shock as an aid to explosion, but note that such a pressure is larger (artificially) in 2D than in 3D. Turbulent pressure is also (possibly) more anisotropic in 2D, favoring the radial stress component, further enhancing the prospects (again, artificially) for overcoming the accretion ram pressure. Hence, the current explosions in 2D may in part, or at times, be numerical artifacts.
Importantly, though Hanke et al. (2012) have speculated that the sloshing motion oftimes associated in 2D with the SASI (standing accretion shock instability) may be crucial to explosion, the recent non-rotating default VERTEX-PROMETHEUS calculations in 3D by that same Garching group (that at times manifest the SASI) do not explode (Hanke et al. 2013;Tamborra et al. 2014), even though the corresponding 2D simulations did. As an aside, Burrows et al. (2012) note that such a pronounced axial sloshing motion is rarely in evidence in 3D. However, with altered microphysics, in particular with a change in the axialvector coupling constant (g A ) due to a speculative enhanced effect on the nucleon spin of the strange quark, Melson et al. (2015) do obtain a weak explosion in 3D of the Woosley & Heger (2007) 20-M ⊙ in 3D. The difference in outcome is due to the consequent decrease in the neutrino-neutron neutral-current scattering cross section in the neutron-rich envelope of the proto-neutron star bounded by the stalled shock and the resultant increase in the electron-neutrino luminosity and average energy that are instrumental in heating this envelope. However, the magnitude of the strangeness correction employed by Melson et al. may be larger than experiment allows (Ahmed et al. 2012;Green et al. 2017). Furthermore, it is not clear that all these currently published under-powered 3D explosions will actually succeed as an explosion after tranversing the entire star and after the mass cut between ejecta and residual core is determined hydrodynamically. An exception is the recent work of Müller et al. (2017), who perform 3D multi-group simulations using their simpler FMT method (Müller & Janka 2015) for an 18-M ⊙ progenitor model evolved in 3D to collapse. Such a 3D progenitor naturally manifests seed perturbations that Müller et al. (2017) show are instrumental in leading to an explosion with a respectable energy.
The calculations of Burrows et al. ( ,2007a and Ott et al. (2008), using the VULCAN/2D code, and Dolence, Burrows, & Zhang (2015), using the CAS-TRO code (Zhang et al. 2011(Zhang et al. ,2013 employed multi-dimensional transport, and not ray-by-ray+, but in neither of these 2D studies did the authors see explosions by the neutrino mechanism. VULCAN/2D did not have all the terms to order v/c in the transport, but CASTRO did, and the results were similar (Dolence, Burrows & Zhang 2015). Importantly, neither VULCAN/2D nor CASTRO made corrections for the effects of general relativity, and this could explain in part the different outcome. However, since other self-consistent calculations (summarized previously) that obtained explosions in 2D used the ray-by-ray scheme, while the VULCAN/2D and CASTRO studies did not, one is tempted to suggest that the ray-by-ray approach may in 2D be yielding qualitatively incorrect results. Suwa et al. (2010) obtain an explosion in 2D of a 13-M ⊙ progenitor, while Takiwaki et al. (2012) obtain an explosion for an 11.2-M ⊙ progenitor in both 2D and 3D. Both these efforts, however, neglect ν µ and ν τ neutrinos, which constitute roughly 50% of the total neutrino losses after bounce, and use the IDSA (Liebendörfer et al. 2009) and ray-by-ray approximations for the transport.  continued their explorations in 2D with a large suite of progenitor simulations, but continue to neglect ν µ and ν τ neutrinos and use the IDSA plus ray-by-ray+ transport methods. They highlight a possible role in the explosion systematics of transitions in the mass accretion rate.  perform fully general-relativistic 3D simulations, employing a leakage scheme for all the neutrinos. Their goal was to explore the relative role of the SASI and neutrinodriven convection, and they found for their simulation of the 27-M ⊙ progenitor of Woosley, Heger, & Weaver (2002) that neutrino-driven convection dominated once it started. Roberts et al. (2016), using full GR hydrodynamics and an M1 transport scheme ( §2) roughly similar to ours (but without the velocity-dependent terms in the transport, inelastic scattering, or many-body effects), emphasize the importance of spatial resolution in determining whether and how the same 27-M ⊙ model explodes, as well as the different outcomes for octant and full 4π steradian simulations. Kuroda et al. (2016) have achieved a fully general-relativistic 3D code, also using the M1 transport closure, but have yet to simulate progenitor models beyond a few tens of milliseconds post-bounce. Pan et al. (2016) have constructed a 2D non-relativistic (Newtonian) neutrino radiation/hydrodynamic scheme using the IDSA transport approach for the ν e andν e neutrinos and a leakage scheme for the ν µ neutrinos. Importantly, they do not use the ray-by-ray+ approach, but follow the transport of ν e andν e neutrinos multi-dimensionally. They obtain explosions for all the progenitor models studied. Suwa et al. (2010) and Nakamura et al. (2014) found that rotation aided explosion, mostly by rotationally expanding the size of the gain region and increasing the mass it contained. Iwakami, Nagakura, & Yamada (2014) and  highlight the rotational excitation of m = 1 spiral-arm modes and find a role for non-axisymmetric rotational instabilities. Iwakami, Nagakura, & Yamada (2014) used a light-bulb neutrino scheme, and neglected "ν µ " neutrinos. Both Nakamura et al. (2014) and  observed equatorial explosions in the rapidly-rotating context. Recently,  and Summa et al. (2018) published a rapidly-rotating 3D model that exploded and highlighted the potential role of an "m = 1" structure in the gain region. Earlier, Fryer & Heger (2000) used SPH for the hydrodynamics and a simplified gray scheme for the neutrino transport to explore the role of rapid rotation. Moreover, in all of these studies the initial rotation rate was not only high in the mantle, but was high in the core. Such rapid initial spins (periods of a few seconds) and final spins (periods of ∼2−10 milliseconds) seem not to be consistent with inferred pulsar spin periods at birth (crudely, ∼300 ± 200 milliseconds; Emmering & Chevalier 1989;Faucher-Giguere & Kaspi 2006;Popov & Turolla 2012;Noutsos et al. 2013) and may be associated only with hypernovae (Burrows et al. 2007c) and/or gamma-ray bursts (MacFadyen & Woosley 1999).
In this paper, we explore, with self-consistent 2D Fornax ( §2) simulations, the dependence of the outcome of collapse (most notably whether the model explodes) on neutrino-nucleon scattering rates (via modifications of in-medium response corrections due to many-body effects), pre-collapse convective perturbations, inelastic neutrino-electron scattering, and inelastic neutrino-nucleon scattering. We also continue our study, started in Skinner et al. (2016), of the issues raised by the use of the ray-by-ray+ method. What we find is that when the proto-neutron star bounded by a stalled shock is close to the critical condition for explosion (Burrows & Goshy 1993), as it easily can be in the turbulent multi-D context, the sensitivity to explosion of small changes in the physical inputs is amplified. The magnitude of such changes might be only ∼20%, but the result can be qualitatively different, in particular whether the model explodes. In the 1D (spherical) case, the core is not sensitive to comparable changes ("Mazurek's Law"), but in the multi-D turbulent (and chaotic) context, small changes in the physics can have a qualitative effect on explodability. This may explain why the various groups around the world simulating core-collapse supernovae can witness very different outcomes, despite the fact that they ostensibly are incorporating almost the same microphysics and similar computational approaches. Small differences are amplified near criticality in this chaotic context. The availability of a new generation of fast, but accurate, simulation codes, such as Fornax ( §2), and significant supercomputer resources enable rapid multi-parameter investigations in the multi-dimensional (in particular 2D), multi-physics context. Such wide-ranging explorations reveal patterns not easily discerned when one (or only a few) simulations are the focus of a paper and its (or their) results solely are mined.

Numerical Method and Computational Setup
We have developed an entirely new multi-dimensional, multi-group radiation/hydrodynamic code, Fornax, for the study of core-collapse supernovae. This code is described in detail in an upcoming paper (in preparation). For the purposes of this paper, we note that it employs spherical coordinates in one and two spatial dimensions, solves the comoving-frame, multi-group, two-moment, velocity-dependent transport equations to O(v/c), and uses the M1 tensor closure for the second and third moments of the radiation fields (Dubroca & Feugeas 1999;Vaytet et al. 2011). 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.
The hydrodynamics in Fornax is based on a directionally-unsplit Godunovtype finite-volume method. Fluxes at cell faces are computed with the fast and accurate HLLC approximate Riemann solver based on left and right states reconstructed from the underlying volume-averaged states. In the interior, to alleviate Courant limits due to converging angular zones, the code can deresolve in angle with decreasing radius, conserving hydrodynamic and radiative fluxes in a manner similar to the method employed in AMR codes at refinement boundaries. Gravity is handled in 2D and 3D with a monopole or a multipole solver (Müller & Steimetz 1995). When using the latter, we generally set the maximum spherical harmonic order necessary equal to twelve. The monopole gravitational term is altered to accommodate approximate general-relativistic gravity (Marek et al. 2006), and we employ the metric terms, g rr and g tt , derived from this potential to incorporate general relativistic redshift effects in the neutrino transport equations (in the manner of Rampp & Janka 2002). In 2D, rotation and a third component of the velocity vector can be included in the hydrodynamics. We use the SFHo equation of state (EOS) by default (Steiner et al. 2013), but in this paper do compare with results using the DD2 EOS (Banik et al. 2014) For these simulations, we follow twenty energy (ε ν ) groups for each of the ν e , ν e , and "ν µ " species (again, the four species, ν µ ,ν µ , ν τ , andν τ lumped together.) For the ν e types, the neutrino energy ε ν varies logarithmically from 1 MeV to 300 MeV, while it varies from 1 MeV to 100 MeV for theν e s and ν µ s. We have performed calculations with forty energy groups and found little difference in the results. The radial coordinate, r, runs from 0 to 20,000 kilometers (km) in 608 zones. The radial grid smoothly transitions from uniform spacing with ∆r = 0.5 km in the interior to logarithmic spacing, with a transition radius near ∼100 km. The polar angular grid spacing covers the full 180 • and varies smoothly from ≈ 0.95 • at the poles to ≈ 0.65 • at the equator in 256 zones.
A comprehensive set of neutrino-matter interactions are followed in Fornax, and these are described in Burrows, Reddy, & Thompson (2006). Inelastic neutrinonucleon scattering is handled using a modified version of the Thompson, Burrows, & Pinto (2003) approach ( §4). Our neutrino microphysics is more comprehensively listed in our upcoming code paper.

Many-Body Neutrino Response Corrections
Melson et al. (2015) invoked a modification in the axial-vector coupling constant (g A ) due to a possible strangeness contribution to the nucleon spin of g s A = −0.2. This results in an approximate decrease in the neutrino-nucleon scattering rate of ∼20% and in Melson et al. the consequence was an explosion in 3D, even though they did not witness an explosion in 3D when using their default microphysical suite. Curiously, without their strangeness correction, the same model exploded in 2D. However, the value of the correction, g s A , they employed in their 3D model is likely too large and g s A may be closer to zero (Ahmed et al. 2012;Green et al. 2017).
Many-body corrections to neutral-current and charged-current neutrino-nucleon interactions have been discussed in the past (Burrows & Sawyer 1998Hannestad & Raffelt 1998;Reddy et al. 1999;Roberts et al. 2012) in the context of protoneutron stars and supernovae, and have been known to affect the neutrino-matter reaction rates. Burrows & Sawyer (1998) in particular suggested that many-body corrections to the axial-vector and vector structure factors for neutrino-nucleon scattering could be of a magnitude sufficent to be of relevance to the viability of the neutrino-driven mechanism of core-collapse supernovae, but did not provide a robust estimate of the magnitude of this diminution much below nuclear densities.
The effect of the many-body correction to neutral-current neutrino-nucleon scattering in the supernova context is mostly due to the increase in the ν µ luminosity occasioned by the decrease in the associated ν µ + (n, p) → ν µ + (n, p) scattering cross section for densities above ∼10 12 g cm −3 ; this causes a further compression in the core. Such a compression, similar to the effect of GR, raises the temperatures near the ν e andν e neutrinospheres, thereby raising their associated luminosities and average emergent neutrino energies. These changes increase the neutrino-matter heating rates in the gain region and, hence, facilitate explosion. Since the super-allowed charged-current absorption reactions still dominate the ν e /ν e -matter interaction rates, the direct effect of this structure factor correction to the axial-vector term in the neutrino-nucleon scattering rate on the ν e andν e luminosities is slight.
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 g cm −3 was recently derived by Horowitz et al. (2017) using a virial approach. These authors derive a fit: where In these equations, T is the temperature in MeV, Y e is the electron fraction, n is the baryon density in fm −3 , A 0 = 920, B 0 = 3.05, C 0 = 6140, and D 0 = 1.5×10 13 . Horowitz et al. (2017) join their fit to Burrows & Sawyer (1998) for the higher densities, but were most careful fitting their formula for temperatures between 5 and 10 MeV. Nevertheless, though Horowitz et al. (2017) intended their fit to apply at all densities, temperatures, and Y e s, (and we use it in this paper at all thermodynamic points), one should be aware that no current approach to this physics is likely to be correct at the highest densities above ∼10 14 g cm −3 . Therefore, use of this formula in supernova and proto-neutron-star cores should be considered provisional. Fortunately, since such densities are too high to affect the evolution during the first second post bounce, and are most relevant during the later proto-neutron-star cooling phase, the values of the Horowitz correction at the highest densities are not germane to the conclusions of this paper. However, many-body effects in the charged-current sector and on absorption may be comparable (Burrows & Sawyer 1999;Roberts et al. 2012;Fischer 2016), but have not yet been factored in. Therefore, neglected is the effect of final-state nucleon blocking for charged-current absorption reactions. Blocking is important only at high densities, where many-body interaction effects are likely to be even larger. Therefore, we have postponed their inclusion until such corrections, which must for self-consistency be done with the same interaction model that informed the equation of state employed, are available.
The fit derived by Horowitz et al. (2017) to the structure factor applied to g 2 A translates into a decrease in the neutrino-nucleon scattering cross section in the crucial region at and deeper than the various neutrinospheres of ∼5% to ∼35%. This effect is a function of density, temperature, and electron fraction (Y e ). Horowitz et al. (2017) state that the corresponding structure factor for the vector current is likely greater than one, but since the vector contribution to neutrinonucleon scattering is small, this is likely to be subdominant. The upshot is a potentially important, and physically plausible, decrease in the neutrino-matter scattering rates that translates into an increase in the driving ν e andν e luminosities and average energies, thereby increasing the heating rate in the gain region.
We note that the corresponding many-body correction for charged-current interactions may be in the same direction (Fischer 2016;Burrows & Sawyer 1999, although see Roberts et al. 2012) and could also be important ( §5.3). The correction would be small at low densities in the gain region, but higher deeper inside, where the neutrinospheres reside. If the rates are suppressed, this could increase the ν e andν e luminosities, while simultaneously not decreasing the heating in the gain region and is not like a uniform correction at all radii. Such behavior, if it obtains, would be near optimal for aiding the explosion.

Inelastic Scattering
Neutrino-electron scattering rates and cross sections are much smaller than those for neutrino-nucleon scattering, which themselves are smaller still than those for super-allowed charged-current reactions such as ν e +n → p+e − . For ε ν = 10 MeV, this deficit is approximately two orders of magnitude. However, due to the small mass of the electron, the energy transfer to the matter during a "Compton-like" neutrino-electron scattering is on average quite large, while, due to the large mass of the nucleon, the corresponding average energy transfer during neutrino-nucleon scattering is rather small. Therefore, the large cross section for neutrino-nucleon scattering multiplied by the small associated energy transfer can be comparable to the product of the small neutrino-electron cross section with the large perinteraction energy transfer and depend upon temperature, density, and neutrino energy (Janka et al. 1996;Thompson, Burrows, & Horvath 2000). The upshot is that inelastic scattering off both electrons and nucleons can modify thermal profiles and heating rates exterior to the neutrinospheres and in the gain region and contribute to explosion by the neutrino heating mechanism. Moreover, Müller et al. (2012b) make the point that heating by ν µ -nucleon inelastic energy transfer can boost the temperatures near the ν e andν e neutrinospheres and results in slightly higher ν e andν e energy luminosities, which in turn enhance heating behind the shock correspondingly.
The detailed theory of the structure functions and redistribution kernels for such inelastic scattering, including the effects of final-state blocking and the thermal spectrum of the targets, can be found in Burrows, Reddy, & Thompson (2006), Reddy et al. (1999), Thompson, Burrows, & Pinto (2003), and Thompson, Burrows, & Horvath (2000). Pioneering work on inelastic scattering off electrons in the core-collapse context can be found in Bruenn (1985) and Mezzacappa & Bruenn (1993). However, those latter papers were focussed on the downscattering effect due to inelastic ν e -e − scattering on the electron neutrinos produced during infall and the consequent decrease in the trapped lepton fraction. Since a larger trapped lepton fraction could help facilitate immediate post-bounce explosions (Burrows & Lattimer 1983), this quantity was more relevant when the prompt hydrodynamic supernova mechanism still seemed viable. However, with the emergence of the delayed, neutrino-driven mechanism (Bethe & Wilson 1985), and the conclusion that the prompt mechanism could not work due to catastrophic neutrino losses at and around shock breakout, the value of the trapped lepton fraction, and its precise value, receded in significance.
Nevertheless, heating behind the stalled shock due to inelastic energy transfer from neutrinos to both electrons and nucleons, or boosting the ν e andν e luminosities by ν µ downscattering near their neutrinospheres (Müller et al. 2012b), may help achieve the critical condition for explosion, particularly when in concert with the inclusion of the in-medium neutrino response (Horowitz et al. 2017; §3) and the slightly net positive influence of GR. Heating by inelastic neutrino-electron scattering is still sub-dominant with respect to that due to charged-current ν e andν e absorption, and in 1D (spherical) simulations there is almost no hydrodynamic consequence of its inclusion (Thompson, Burrows, & Pinto 2003). The same can be said of inelastic neutrino-nucleon scattering (Janka et al. 1996;Burrows & Sawyer 1998Thompson, Burrows, & Horvath 2000). However, in the realistic multi-D context of the core-collapse phenomenon, the cumulative effect of the addition of a few sub-dominant heating mechanisms, compression due to enhanced ν µ neutrino losses and consequent ν e andν e neutrinosphere heating 3 , and greater proximity to the critical condition due to multi-D effects 4 amplifies the leverage of even small additions to the gain-region heating power over their effects individually and can convert an anemic explosion into a robust explosion, or even a dud into a blast. One recalls that the original delayed mechanism of Wilson required for explosion only a modest enhancement of ∼25% in the neutrino luminosity 5 .
From the work of Thompson, Burrows, & Horvath (2000) and Tubbs (1979), we find that the crossover energy between upscattering and downscattering is nearer 6k B T (not 3k B T , as in Müller & Janka 2015), where k B is Boltzmann's constant and T is the temperature, around densities of ∼10 11 g cm −3 to ∼10 13 g cm −3 and so our approximate redistribution rate for ν µ s is proportional to κ scat (ε µ − 6k B T )/m n c 2 , where κ scat is the ν µ scattering opacity and m n is the neutron mass.
Inelastic scattering off nucleons lowers the ν µ luminosity, while the corresponding quantities for the ν e andν e are slightly increased (Müller et al. 2012b). The latter responses explain the positive effect on explodability of the inclusion of inelastic scattering off nucleons. The direct effect of such inelasticity is in the core, and the effect in the shocked mantle is indirect.

Microphysical Dependences of the 2D Explosion Models
For a set of 4 fiducial models (with progenitor masses of 13, 15, 16, and 20 M ⊙ ), we performed a parameter study of input physics highlighting the role of small changes in promoting (or preventing) explosion. We refer the reader to the upcoming paper by Vartanyan et al. (2018, in preparation) for a more detailed description of these simulations. When models explode, the time to explosion can vary significantly with inputs. The relative time of explosion can help one gauge the role of the inputs in question 6 . While there are many combinatoric possibilities, we settled on just a few comparisons to demonstrate the effects we see universally. We employ the notation IES, INS, MB, pert, and rbrp to indicate "inelastic neutrino-electron," "inelastic neutrino-nucleon," "many-body," "perturbations,", and "ray-by-ray+," respectively.

Inelastic Scattering
Models with only inelastic scattering off electrons (IES), inelastic scattering off both electrons and nucleons (IES INS), and additionally with the many-body effect (IES INS MB) can illuminate the respective roles of each 7 . Figure 1 compares outcomes for the 16-M ⊙ progenitor. We toggled the role of inelastic scattering off electrons and nucleons as well as the many-body effect for the 16 M ⊙ progenitor of Woosley & Heger (2007), which we found to explode early at ∼300 ms using the default setup IES INS MB. However, performing the simulation with any of these components removed prevents explosion. Figure 2 depicts the role of IES, INS, and MB on the emergent spectra at two different times after bounce. One of the central results that can be gleaned is the time of explosion. This quantity can help one gauge the relative role of the inputs in question. The default model (def) here does explode. Prior to explosion, both inelastic scattering off electrons and off nucleons and the many-body correction led to higher spectral fluxes. After the default model explodes, the flux spectrum diminishes relative to that of the non-exploding model 16 def noMB. The slight increase in the heating in the gain region prior to explosion due to the inclusion of inelastic neutrino-electron and neutrino-nucleon scattering has certainly helped the core achieve the critical explosion condition. Figure 3 compares the emergent neutrino luminosities (left panel) and root-mean-square (rms) neutrino energies (right) for three models. The boosting in the ν µ luminosity and rms energy is clearly manifest, as are the corresponding boosts in those quantities for the ν e andν e neutrinos. While one may have speculated that enhanced neutrino losses, particularly due to ν µ s that are almost ineffectual in heating the shocked mantle, 6 The "time of explosion" is sometimes defined as the time the mean shock radius achieves 400 km. We think this definition somewhat arbitrary, but acknowledge some degree of arbitrariness in any definition, however useful. We here define the time of explosion as the approximate time at which the curve of the mean shock radius versus time is inflected upwards. All our models that inflect in this way explode and all these models have been continued to at least one second after bounce. The mean shock radii of these exploding models all achieve radii beyond 5000 kilometers by the end of the simulation.
7 In a previous version of this paper, we had employed an opacity file that was compromised by a compiler bug. The net result of this error was a slightly enhanced magnitude of the neutral-current many-body effect. This error has now been fixed. . This model explodes at ∼250 ms post-bounce. We then remove and add certain inputs, denoted by a subscript with "def". Removing either the many-body correction (blue, "16 def noMB") or inelastic scattering off nucleons (green, "16 def noINS") leads to a dud. However, even without the many-body correction (noMB), adding either perturbations (red) ( " pert"; §6) or modifying the opacity table to include Fischer's (2016) correction to the nucleon-nucleon bremsstrahlung rate ("bf") and only 20% of the electron capture rate (Juodagalvis et al. 2010) on heavy nuclei (orange, "0.2j"), leads to an explosion ∼50 and ∼100 ms, respectively, after our default model. This helps illustrate the sensitive dependence of the outcome − explosion or dud − on the microphysical inputs when near criticality.
would have had a negative effect on explodability, the converse is true. Greater losses lead to a further compaction of the core with an increase in the matter temperatures near the ν e andν e neutrinospheres. The result is similar in effect to that of GR, whereby such core heating enhances the driving ν e andν e luminosities and the average neutrino energies, which in turn enhance the heating power in the gain region. Since it is this power deposition that ultimately drives explosions, the net effect is quite supportive of explosion. However, the actual magnitude and form of the correction to the axial-vector coupling term in the expression for the neutrino-nucleon scattering rate may be different from that derived by Horowitz et al. (2017) and this still needs to be verified. Moreover, the effects of similar many-body corrections to the absorption rates need to be incorporated, as do those for the vector coupling strengths. One prediction of the consequence of the structure-factor correction we have employed is the enhanced ν µ luminosities and average energies seen in Figure 3. It is note- worthy that, as it stands, the effect on the outcome of core-collapse of many-body rate suppressions (Burrows & Sawyer 1998) might be large.

Equation of State
Through its control of the pressure for a given temperature, density, and Y e , the equation of state of hot, lepton-rich nuclear matter will determine the structure of the proto-neutron star and its evolution after bounce. The stiffer the EOS, the more extended the core. On the one hand, a stiff EOS will resist the quick increase in temperature near the ν e andν e neutrinospheres, and the consequent enhancement of the neutrino heating rate in the gain region, seen both in the comparison of GR and Newtonian models and in the increase in the ν µ losses due to the many-body effect. On the other, a stiffer EOS will provide a more stable inner platform that won't as easily (by its inward motion with time after bounce) send out weakening rarefactions to the outer bounce shock that could inhibit explosion. Figure 4 indicates (at least for this 16-M ⊙ model) that the first effect seems to win, since the DD2 EOS model does not explode, while the SFHo model (the softer of the two at high densities) does. On this figure we also show that, though the DD2 model didn't explode, the behavior of the shock with time was significantly less vigorous when the inelastic scattering effects were turned off, reiterating the conclusions of §4 and §5.1. We also note that since the shock in the default DD2 model achieved a large mean shock radius (∼200 km) before subsiding, this model was nevertheless very close to exploding.

Nucleon-Nucleon Bremsstrahlung and Electron-Capture on Heavies
On Figure 1, we also provide a model that drops the many-body correction to neutrino-nucleon scattering, as derived by Horowitz et al. (2017), but substitutes in the Fischer (2016) correction to the nucleon-nucleon bremsstrahlung rate ("bf") and divides the electron capture rate on heavy nuclei derived for the SFHo EOS by Juodagalvis et al. (2010) by five ("0.2j"). The former effect lowers the bremsstrahling production rate of ν µ neutrinos, as well as the inverse bremsstrahlung absorption. The latter effect slightly retards the infall rate, thereby decreasing the mass accretion rate post-bounce at a given time. Decreasing this rate can faciliate explosion. As the figure suggests, these microphysical changes can compensate for ignoring the many-body effect to result in a similarly-aided explosion. Clearly, such modest alterations in the microphysics of relevance, still within the envelope of our ignorance, could play a positive role. In Figure 5 below, we see a similar positive effect of these potential changes in the default microphysics for a 13-M ⊙ model. Together with the model in Figure 1, this suggests that a more extensive exploration of such physics would be fruitful.

Progenitor Perturbations
Performing the last stages of stellar evolution before collapse hydrodynamically in 2D and 3D has been shown to alter, perhaps in significant ways, the compositional, entropy, and density profiles of the core (Meakin et al. 2011;Chatzopoulos et al. 2016;Abdikamalov et al. 2016;Müller et al. 2017), and it has long been known that progenitor density profiles have an impact on the outcome of collapse. This is implicit in the critical curve analysis of Burrows & Goshy (1993), whereṀ and the accretion ram pressure play central roles. It is also a factor in discussions of the compactness parameter (O'Connor & Ott 2011) and its extensions (Ertl et al. 2015). In this vein, one notes that at least three groups (Kitaura et al. 2006;Burrows et al. 2007b;Radice et al. 2017) have already demonstrated that the 8.8-M ⊙ "electron-capture" supernova progenitor of Nomoto & Hashimoto (1988), with its extremely steep density ledge, can explode in 1D by the neutrino-driven wind mechanism, though the explosion energy is low (∼1 − 2× 10 50 ergs). However, when spherical progenitor models do not readily lead to explosion, the initial perturbation spectrum in the progenitor's convective silicon and oxygen zones could certainly affect the timescales for the generation of turbulence behind the stalled shock and be a factor in the onset of explosion (Couch & Ott 2013;Müller & Janka 2015;Abdikamalov et al. 2016;Müller et al. 2017). Specifically, the magnitude, character, and spectra of seed perturbations will affect how quickly turbulence reaches the non-linear regime and, perhaps, whether turbulence grows to non-linearity at all during the finite time the accreta are in the unstable gain region.
Therefore, introducing physically-motivated seed perturbations that originate from and reflect the three-dimensional character of the convective core of an actual massive star at its terminal stage of evolution is a topic of some interest. However, most calculations done to date do not start with true 3D convective structures, but with 1D models from the literature, and impose either ad hoc perturbations in density or velocity randomly at the grid level or allow grid asphericities (such as obtained when using a Cartesian grid) or truncation errors to act as seeds. Neither of these approaches is physical, and the resulting initial perturbation spectra lead to early growth rates in the linear regime that reflect not the multi-D progenitor perturbation structure, but the numerical development of convenient artificial power spectra. This will affect how quickly turbulence reaches the non-linear regime and, perhaps, whether turbulence grows to non-linearity at all during the finite time the accreta are in the unstable gain region. In addition, this may have a bearing on the post-bounce delay to a turbulence-aided explosion, with the consequent effect on the time and energy of that explosion.
Therefore, the seed perturbations that arise during oxygen and silicon burning prior to collapse might be key inputs into core-collapse supernova theory, and Couch & Ott (2013), Müller & Janka (2015), , Müller (2016), and Müller et al. (2017) have begun to explore this. However, mixing-length theory, though inadequate as a comprehensive theory, still provides a measure of the magnitude of velocity perturbations at the onset of collapse , and they are only a few hundred to ∼500 km s −1 , with Mach numbers bounded by ∼0.08 (Woosley & Heger 2007). This is not large.
In this section, we provide a quick glimpse at the possible relative role of significant perturbations on the timing and character of explosion in light of the other physics. To this end, we employ the methodology of Müller & Janka (2015). These authors impose a vector velocity perturbation map on their progenitor that is more realistic than many past approaches and renders a perturbation field that is independent of grid resolution. This (surprisingly) was rarely attempted in the past and provides a specific context for future comparison. We set the maximum perturbation speed on the grid to 1000 km s −1 , which as indicated earlier may be near or beyond the expected upper end of the range, a spherical harmonic index ℓ of 2, and a radial "quantum number" n of 5. Both ℓ and n are parameters in the Müller & Janka (2015) formulation. With this parameter set, we simulate in 2D the self-consistent multi-group evolution of 13-and 15-M ⊙ progenitors and compare the result to a default model for which the initial perturbations are much smaller and arose numerically from grid noise.  (2007) models, showcasing the effect of significant perturbations and moderate rotation. We follow Müller & Janka (2015) prescription in implementing perturbations to radial velocities on infall over three regions with the maximal radial velocity (1000 km s −1 ), n (number of radial convective cells), and l (number of angular convective cells) as parameters. The interior region spans 1000 to 2000 km, outside the nascent core, the middle region 2100 to 4000 km, and the outer region 4100 to 6000 km, truncated roughly where accretion ends after the first second for our simulations. All regions have l = 2, n = 5 and maximum radial velocity of 1000 km s −1 . See text for a discussion. Figure 5 portrays the evolution of the mean shock radii, with and without progenitor perturbations. The left panel shows that the imposed perturbations converted failure into success for the 15-M ⊙ . We also provide on this panel a model including the effect of a modest rotation rate, along with the perturbation. Including the latter resulted in an explosion, though later. Our rotational study has revealed that the effect of rotation is not monotonic with initial internal rotational speed or spatial profile. On the right panel of Figure 5, we show the corresponding behavior for the 13-M ⊙ progenitor, which does not explode in our study for the default microphysics. It does, however, explode when similar perturbations are imposed and when a similar rotational profile is assumed. In this case, rotation promotes explosion, highlighting its non-monotonic effect upon outcome.
We note that the initial perturbations we imposed have an amplitude that is somewhat larger than expected  and that their character is still rather artificial. The magnitude of the initial perturbations may well be lower, but their character and magnitude will certainly vary from progenitor to progenitor. Therefore, a much more thorough study with 3D progenitors and 3D collapse models is called for.

Ray-by-ray+ Anomalies
As shown by Skinner et al. (2016), the ray-by-ray+ approach to neutrino transport, whereby multi-D transport is replaced by multiple 1D transport calculations with corrections for matter advection, but not lateral transport, can introduce systematic errors in the heating rates along the poles in axisymmetric 2D simulations. Such enhancements, when in proximity to criticality, may be producing explosions artificially. At the very least, the time to explosion is artificially shortened, perhaps significantly. Since there is little or no accumulation of explosion energy prior Shock radius (km) versus time after bounce (s) for the 20 M ⊙ WH07 progenitor using the LS220 EOS, with and without the ray-by-ray-plus (rbrp) approximation to neutrino transport. All models include inelastic scattering off electrons and nucleons. We see that including ray-by-ray+ (rbrp) leads to an explosion when otherwise there was none. The model with ray-by-ray+ and without many-body (green) explodes as well, though 700 ms after the model with ray-by-ray+ and many-body, suggesting that the ray-by-ray+, though artificial, is more significant to explosion than the physical inclusion of the many-body correction.
to global instability (Burrows, Hayes, & Fryxell 1995), 8 an earlier explosion may make more of the emitted neutrinos available for explosive driving. Figure 6 compares three 20-M ⊙ models, using this time the LS220 EOS. Here, the default model does not explode, but the one employing the ray-by-ray+ simplification does. Note that, for both ray-by-ray+ models, the model without the many-body (MB) correction explodes much later (by 700 ms) than the model with MB. Our results suggest, for this progenitor, that the artificial ray-by-ray+ approximation is more significant to explosion than the physical many-body effect. One is left to speculate whether 2D simulations in the literature that employ rayby-ray+ would indeed explode if they used more realistic transport. This is all the more relevant in 3D, given that extant published models that do explode in 2D have more difficulty exploding in 3D, a context in which it is not clear that the ray-by-ray+ method introduces as great an artifact as in 2D. It is true that the turbulent pressure spectra in 2D and 3D are different, with the turbulent cascade in 2D resulting in enhanced stresses on larger scales, and that the turbulent-stress boost to explodability may be smaller in 3D. This could also be a factor in the more anemic outcomes in published 3D models. Nevertheless, it may be that the more problematic nature of published 3D models vis-à-vis 2D models is a conse- quence of some combination of the use of ray-by-ray+ and the reduced turbulent stress in 3D.
We conclude this section by noting that the calculations of Skinner et al. (2016) did not include various physical effects (such as inelasticities and the many-body correction) that we highlight here. In Skinner et al., we obtained explosions using ray-by-ray+ for some models that did not explode otherwise, and when they exploded they did so late.

Compactness
Finally, we conclude that although the compactness parameter (O'Connor & Ott 2011,2013Pejcha & Thompson 2015), defined at bounce for a given mass interior (M ) at a given radius (R) as ξ ≡ M/M ⊙ (1000km/R), is one measure of the important density structure of the progenitor, it is not necessarily predictive of explosion, at least during the first second after bounce. We find that, depending upon the neutrino physics employed, the temporal order in which models with different compactness parameters explode after bounce varies. If time of explosion is a fit measure of explodability, then this alone would challenge the usefulness of the compactness concept visà vis explodability. Nevertheless, though the compactness parameter is large for the 21-M ⊙ model and small for the 12-M ⊙ , we found that the former can be more explodable. One might have thought that if compactness were the sole predictor the details of the neutrino-matter interaction would have mattered little in this regard. However, we see that the time to explosion does not necessarily correlate well with compactness. Models with dense envelopes have larger accretion and total neutrino luminosities after bounce, and it is these luminosities that drive explosion. In addition, models with dense envelopes and higher compactness have a greater optical depth to neutrinos in the gain region. Since the neutrino power deposition goes approximately as the product of this depth and the luminosities, high-compactness progenitors have an advantage. Therefore, it is feasible that more massive models might be more explosive, and this trend has been seen by other modelers Bruenn et al. 2016). In fact, in the Summa et al. and Bruenn et al. papers for the Woosley & Heger (2007) models, the post-bounce explosion times increase in the sequence 20, 25, 15, and 12 M ⊙ , which is clearly not monotonic with compactness. This behavior is the reverse of what might be expected if low compactness signaled greater explodability. Finally, the work of Nakamura et al. (2015) suggests that though they see a weak correlation, it is not monotonic with compactness. In fact, many of their plots versus compactness resemble scatter plots.
A critical issue is whether these more massive models with shallower density profiles that explode early can maintain explosion during the traversal of the shock through the outer stellar mantle. The binding energy penalty of the outer envelope generally increases with progenitor mass and might be too steep a price to pay during subsequent evolution for those more massive cores that explode earlier and, perhaps, more energetically. As previously stated, for most of the relevant mass function, the envelope binding energy exterior to a given interior mass is an increasing function of progenitor mass. It is this "barrier" that may set the limit to the range of massive stars that can explode and leave behind neutron stars, although it cannot be excluded that the progenitor mass range that yields neutron stars, and not black holes, may be discontinuous (Sukhbold et al. 2016). Figure 7 portrays various exterior binding energies and Figure 8 follows with a depiction of the close correspondence between the compactness and the envelope binding energy exterior to a baryonic mass cut of 1.5 M ⊙ . Therefore, we contend that whatever significance there may be to the compactness parameter is likely due to its correspondence with the binding energy of the outer envelope exterior to a given mass cut and that compactness need not correlate with the apparent explodability during the first second after core bounce. Importantly, the ultimate outcome will depend upon the progress of the shock at post-bounce times that generally exceed those to which most core-collapse simulations currently go.

Conclusions
In this paper, we have generated and explored detailed 2D (axisymmetric) approximate GR models of core-collapse supernovae using the new code Fornax, treating all the relevant physics to determine the dependence of the mechanism of explosion and its timing on the physical and numerical inputs and assumptions. These include inelastic neutrino-electron and neutrino-nucleon scattering, many-body/structure-factor corrections to the neutrino-nucleon scattering rates, nucleon-nucleon bremsstrahlung, electron capture on heavies, and physically-motivated initial perturbations. We have also reexamined in brief the effect of using the ray-  , calculated at various interior masses, versus progenitor ZAMS mass, as well as the corresponding envelope binding energy (blue dots; in Bethes [10 51 ergs]) for a baryon mass cut of 1.5 M ⊙ (see Figure 7). The progenitor models of Sukhbold & Woosley (2014) (also in Sukhbold et al. 2016) are used. As this figure shows, whatever the position at which compactness is defined, it correlates extremely well with envelope binding energy. It is our contention that it is the latter quantity that is more germane to the outcome of core collapse.
by-ray+ simplification to neutrino transport. We found that much of the wide variation between the results obtained by different groups around the world simulating stellar collapse (as well as whether the core explodes at all) might be explained by slight variations at the ∼10-30% level in the microphysical inputs when the models are near the critical condition for explosion. In the process, we gauged the relative importance of otherwise sub-dominant neutrino physics processes to the outcome of collapse. Proximity to criticality amplifies the dependence upon small changes in the neutrino sector that translate into slight, but crucial, changes in the emergent luminosities and average neutrino energies and the consequent post-shock heating rates; such sensitivity is not manifest in 1D simulations for which the core is (most often) far from explosive.
Thus, "Mazurek's Law" of severe feedback under variations in neutrino cross sections and rates is overturned due to the proximity to instability possible in the realistic multi-D turbulent context. While alterations in the neutrino coupling rates have little effect in 1D, in multi-D the hydrodynamic response to even small changes and/or corrections to neutrino interaction rates can be more substantial due to greater proximity to the critical curve.
The upshot is that small variations between the methods, microphysics, and resolutions used by groups who ostensibly are incorporating the "same" inputs translates naturally into post-bounce explosion time differences that can range by many hundreds of milliseconds, and in some cases can turn a dud into an explosion (or vice versa). We suggest that this thereby explains in large measure the apparent heterogeneity in the outcomes of detailed simulations performed internationally. A natural conclusion is that, viewed correctly, the different groups are collectively closer to a realistic understanding of the neutrino-driven mechanism of supernova explosion than might have seemed apparent and that a push to rationalize approaches and understand microphysical details in the neutrino-matter interaction sector and the nuclear equation of state could bring a resolution to the decades-long quest for a predictive model of core-collapse supernova explosions.
We have found that many-body corrections to neutrino-matter interaction rates, even at sub-nuclear densities, can have similar effects as general relativity, progenitor perturbations, or inelastic neutrino-electron and neutrino-nucleon scattering. However, we caution that the actual magnitude and functional form of all the various many-body corrections to the neutrino-matter rates (both neutraland charged-current), however important they seem from our current simulations, still need to be explored and verified.
We performed a test using the ray-by-ray+ approximation to neutrino transport in a manner similar to that employed by Skinner et al. (2016) to gauge its effect on the outcome of collapse when a full physics suite is employed. In Skinner et al. (2016), it was shown that ray-by-ray+ artificially enhanced heating along the poles in synchrony with the axial sloshing seen in 2D simulations and thereby made models more explosive. The shift in the explosion times can be as large as the full range currently witnessed by the various supernova simulation groups for a given progenitor. One can speculate that had groups that use the ray-by-ray+ simplification used real multi-D transport their 2D models would have exploded later, perhaps much later or not at all. Speculating further, one wonders whether the fact that 3D models explode later than the corresponding 2D models (Lentz et al. 2015) or not at all (Melson et al. 2015, if without their strangeness fix) may be connected with their use of ray-by-ray+, since non-rotating 3D simulations seldom manifest the axial sloshing seen in 2D. In short, without ray-by-ray+, it is not clear that 2D would explode much earlier than 3D. However, the differences between the 2D and 3D convective cascades may be equally in play here and the associated simulations still need to be performed to assess this.
The possible role of initial perturbations as seeds to the growth of convective instability in and around the gain region has been a subject of recent focus. Clearly, allowing grid noise, truncation error, or other numerical noise to initiate linear growth to the non-linear phase, or imposing artificial initial perturbations, is less than satisfactory. This is particularly true given that seeds have a finite time to grow after accreting through the shock before leaving the unstable region and given that the time to instability and explosion is germane to which phase of the neutrino light curve is driving explosion. Inaugurating the non-linear convective phase early and maintaining it until explosion may be important, and the magnitude and timing of convective seeding is therefore of interest. Though we have deferred until later a more comprehensive study of this subject, we tested the effect of adding to the progenitor velocity perturbations whose magnitude was informed by mixinglength theory . In fact, we allowed the maximum perturbation speed to be 1000 km s −1 , with a Mach number as high as ∼0.12, which is a bit larger than found in 1D progenitor models (cf. Woosley & Heger 2007). Indeed, we found that large imposed perturbations of this sort could enable explosion, though the magnitude of the imposed perturbation seems large. Clearly, it will be important to determine their character, and the many 3D-progenitor studies now in process promise to do just that.
Many-body corrections to scattering rates, inelastic scattering, decreases in nucleon-nucleon bremstrahlung rates and in electron capture rates on heavy nuclei, and initial seed perturbations all boost "explodability" and shorten the time to explosion. In fact, these effects add synergistically and non-linearly to aide explosion, despite the fact that they individually amount to effects at the ∼10−30% level. This is due to the proximity of multi-D models to criticality, and is not seen in 1D simulations.
Furthermore, the turbulence behind the shock that has been shown to aid explosion in the realistic multi-D context naturally introduces indeterminacy in detail. Even the same progenitor star, but with different random seed perturbations and rotational structures at collapse, should yield a range of explosion energies, nucleosynthesis, 56 Ni yields, pulsar kicks, and explosion morphologies. Therefore, it is expected that Nature provides distribution functions, and not one-to-one maps, in all signatures of explosion. In the long run, theory will need to come to grips with this, but in the short run one should not expect that in the chaotic context of turbulent convection and pre-collapse structures the best models will correspond in detail. This is the physical and natural consequence of turbulence and chaos, and will be paralleled in comparison verification studies.
The next stage is to explore the same issues in three dimensions, and one expects there to be important differences Lentz et al. 2015;Müller 2015;Müller et al. 2017). It is only after performing such simulations, and their subsequent verification, that a robust resolution to the core-collapse supernova problem can be claimed. Nevertheless, there has been significant progress of late in unraveling this central mystery in astrophysics, in demonstrating the viability of the neutrino-driven mechanism of explosion, and in illuminating its component physics.
Center for Supercomputing Applications. This paper has been assigned a LANL preprint # LA-UR-16-28849.