Essential ingredients in core-collapse supernovae

With a visual brightness that can compete with their entire host galaxy, supernovae are among the most energetic astrophysical phenomena. As “new stars” that make a sudden appearance in the sky, supernovae have been reported throughout at least 2 millennia of recorded history. The ejecta of a supernova delivers 10⁴⁴ J of kinetic energy and a rich mix of recently synthesized elements into the interstellar medium, making these explosions the principal foundry of the heavy elements that have been formed since the Big Bang, a major source of heat in the interstellar medium (ISM), and a potential trigger for star formation.^1,2 Observationally, these explosions are categorized based on the presence or absence of hydrogen in their spectra soon after the explosion into two types, Type I and Type II. Numerous sub-types are based on the shape of the light curve and the presence or shape of other spectral lines.^3,4 While Type Ia supernovae, classified by the absence of hydrogen lines but the presence of strong Si II lines in their early-time spectra, are believed to be powered by the thermonuclear ignition of white dwarves^5–8 in binary systems, the remaining supernova types and subtypes are thought to result from the deaths of individual massive stars. The wide variations between types and sub-types are the result of differences in the stellar envelope as a result of stellar evolution and mass loss.

To understand the deaths of massive stars, one must begin long before the explosions that mark these deaths and distribute these materials into the interstellar medium. These powerful explosions are the inevitable consequence of the stellar lifecycle. From its initial contraction to form a hydrogen-burning main sequence star, nuclear burning, in the stellar core and in shells surrounding the core, provides the energy to resist gravitational contraction and power the star. Once a star exhausts the hydrogen in its core, the now helium-rich core contracts. This ignites hydrogen in a shell lying atop the core, boosting the luminosity of the star and driving the stellar surface outward to become a Red Giant. Continued contraction of the core is halted by the ignition of helium burning. This evolution sets the pattern for the star's life, with the ash of each stage becoming the fuel for its successor.

In stars below a critical mass — 8 solar masses (hereafter denoted M_⊙) is commonly adopted,⁹ though some argue the number may be as small as 6 M_⊙¹⁰ — degeneracy in the Carbon-Oxygen (CO) core prevents the ignition of carbon burning. Such low-mass stars end their lives as cooling CO white dwarves with masses less than 1.1 M_⊙, surrounded by the envelope they ejected as stellar winds while Asymptotic Giant Branch (AGB) stars. For much more massive stars, from perhaps 11 M_⊙ up to roughly 100 M_⊙,¹¹ hydrostatic carbon, neon, oxygen, and silicon burning leave a core composed of iron, cobalt, nickel, and neighboring species, referred to as the iron-peak nuclei because they populate a prominent peak in the solar system abundance distribution. Collapse ensues once this core grows beyond the maximum stable mass for a system supported by electron degeneracy pressure (the so-called Chandrasaekhar mass). For stars in the lower portion of this mass range, the collapse produces a neutron star and a shock wave that disrupts the stellar envelope, leading to a supernova. In the higher mass region, a black hole forms, preventing the shock wave from being sufficiently energized to disrupt the stellar envelope, and the supernova fails. However, for a rapidly-rotating progenitor star, an alternate mechanism, the collapsar mechanism¹² — driven by an accretion disk surrounding the black hole and producing jets along the rotational axis — may produce a peculiar, hyper-energetic supernova explosion (with explosion energies as much as an order of magnitude larger) and an associated gamma-ray burst. Stars in the intermediate range, between 8 M_⊙ and perhaps 11 M_⊙,¹³ successfully ignite carbon and neon burning, but under degenerate conditions, leading to massive Oxygen-Neon-Magnesium (ONeMg) white dwarves. Near the upper border of this range, above 10.5 M_⊙ in some models,¹⁴ simulations suggest that shell burning causes the mass of the ONeMg core to increase beyond the stable Chandrasaekhar mass before the star loses its envelope. Here too, the core collapses, producing a supernova.¹⁵ Two critical distinctions differentiate the supernovae of these ONeMg cores, also known as electron-capture supernovae (ECSN), from their iron-core brethren: (1) a collapse driven by electron capture on ²³Na, ²⁰Ne, and ²⁴Mg rather than on iron-peak nuclei and (2) a low-density helium shell surrounding the core in the former case, in place of dense layers of silicon and oxygen in the latter.

The collapse of the core of the star, whether composed initially of iron or oxygen–neon–magnesium, proceeds until the stellar core reaches densities similar to those of the nucleons in a nucleus, whereupon the repulsive center of the nuclear interaction renders the stellar core incompressible, halting the collapse. The inner core, which has collapsed homologously, bounces and launches a shock wave into the supersonically falling overlying layers, driving these layers outward. However, this bounce shock is enervated by the escape of neutrinos and nuclear dissociation, stalling the shock before it can drive off the envelope of the star.¹⁶ The failure of the bounce shock to drive off the stellar envelope sets the stage for a delayed mechanism. More than 10⁵³ ergs of gravitational binding energy is released from the newly-formed proto-neutron star (PNS) at the center of the explosion, 100 times the observed kinetic energy of the supernova explosion. Past simulations^17,18 demonstrate that this energy, streaming out from the PNS in the form of neutrinos and antineutrinos of all three flavors (electron, muon, and tau), can be deposited behind the shock and may revive it. Under this neutrino reheating paradigm, the shock remains an accretion shock until sufficiently reenergized by neutrino heating to overcome the gravity of the PNS and the ram pressure of the infalling matter, whereupon it propagates outward, shock heating and transmuting the overlying layers and ejecting the stellar envelope.

From the discussion above, many of the forms of physics operating in core-collapse supernovae are immediately clear; fluid dynamics, neutrino radiation transport with neutrino-matter interactions, gravity (in its full General Relativistic form) and thermonuclear kinetics. However, the question of what of these physics ingredients are truly essential to the core-collapse supernova mechanism, and at what level of detail, is a subject of ongoing debate among the modelers of these events. This question is made important by the increased computational cost of improved physical fidelity, leading modelers to adopt the least expensive implementation of each ingredient that they believe has sufficient fidelity. In the following subsections we will highlight several of these ingredients and the required fidelity.

While a prodigious amount of neutrino energy emerges from the PNS, the neutrinos are weakly coupled to the material directly below the shock, thus the neutrinos travel over significant distances before interacting with matter. Simulating this neutrino radiation problem requires sophisticated neutrino radiation transport. Unlike in photon radiation transport, neutrinos obey Fermi-Dirac statistics and both their energy and number (lepton number) are conserved. These complications are somewhat ameliorated by the absence of lines in neutrino interactions, reducing the level of spectral resolution required. The neutrino heating, which occurs primarily by the absorption of electron neutrinos and antineutrinos on the dissociation-liberated free nucleons, is very sensitive to the spectral distribution of neutrinos^19–22 and direction of propagation (specified uniquely by two angles), at any given spatial point behind the shock.^23,24 Ultimately this requires spectral, multi-angle (Boltzmann) neutrino transport in order to fully compute the neutrino distributions in energy and angle in this region. This renders neutrino transport in core-collapse supernovae a potentially six dimensional (physical space plus neutrino phase space) problem, one that is too expensive to solve even on current supercomputer architectures with sufficient phase space resolution.

Intellectual and computational progress has led to successively better treatments of the neutrino transport in simulations, beginning from the pioneering work of Colgate and White.²⁵ Leakage schemes,^26,27 where the local neutrino emission rate and the opaqueness of the overlying material are used to estimate the cooling rate and hence the neutrino luminosity, and two-fluid schemes,^28,29 where the neutrinos are treated as a Fermi-Dirac fluid with a temperature different than the local matter temperature, are the computationally cheapest. A different technique to reduce the computational cost is to replace the neutrino transport calculation in and around the proto-neutron star with a proscribed neutrino luminosity. This neutrino “lightbulb” approximation, like the leakage and two-fluid schemes, does not include the self-consistent feedback provided by complete transport methods. A separate approach sometimes taken to reduce the computational cost is to track average neutrino energies (gray transport) instead of the discretized spectrum, however such methods are challenged by the steep energy dependence of the neutrino-matter interactions.

Because the neutrino energy deposition behind the shock depends sensitively not only on the neutrino luminosities and spectra, but also on the path lengths traversed by the neutrinos as they propagate outward (and hence the neutrino angular distributions), spherically symmetric simulations implementing spectral Boltzmann neutrino transport were developed.^30–32 One example of these is the Agile-BOLTZTRANcode,^33–35 some results from which we will present in later sections. However, spherically symmetric models with conventional physics have generally failed to produce explosions because they do not deliver sufficient energy to the envelope as a result of the stratification imposed by spherical symmetry.^32,36–40 One useful and less costly approximation are moments methods, where in the neutrino distributions in angle are approximately represented by the lowest order angular moments of the neutrino distribution function. The reduced angular dimensionality (a few moments rather than hundreds of discrete angles) makes today's multidimensional radiation hydrodynamic simulations possible^41–48 Flux-limited diffusion (FLD), a low order moments transport method, wherein the diffusive flux is limited to prevent superluminal neutrino propagation as the mean free path becomes infinite at low density, is widely used in its spectral (also called multigroup) version (MGFLD).^49–52 Moment methods with computed closures from simplified full Boltzmann transport solves⁴⁰ or analytic closures like M1^53,54 provide better estimates of the evolution of the neutrino angular distribution from the diffusive, isotropic interior to the free-streaming exterior than FLD but at less expense than Boltzmann.

Beyond the general treatment of the neutrino transport, seemingly more subtle issues can result in significant differences. An example is the treatment of velocity-dependent ‘observer corrections', studied among other issues by Lentz et al.⁵⁵ using Agile-BOLTZTRAN. Solving the transport equation requires a choice of both coordinates and frame. Lagrangian coordinates moving with the fluid are natural for 1D and particle-based codes, while fixed Eulerian coordinates are natural for multidimensional grid-based codes since multidimensional fluid flow would eventually stretch and twist Lagrangian grids. A frame locally co-moving with fluid may be chosen, which is the frame where the radiation-matter interactions are defined. However, this requires numerous velocity-dependent terms on the left-hand side of the transport equation (see Eq. (1)). Alternately, the external observer (or lab) frame in which there are no velocity-dependent terms on the l.h.s. may be chosen, but the interactions that compose the collision integral must be transformed from their natural co-moving frame. Including the effects of fluid motion (and General Relativity) in the transport of photons or neutrinos requires either a transformation of the collision term or the l.h.s. of the transport equation.

To better understand the observer corrections we start with the co-moving form (E₀ and μ₀ as the energy and cosine of the neutrino propagation angle seen by a co-moving observer) of the

$O(v)$ spherical transport equation for the invariant neutrino distribution function, f, in Eulerian coordinates

$(r,\tilde{t})$

$(r,t̃)$⁠:

The velocity-dependent terms with the E₀ and μ₀ derivatives include corrections to the energy and angle due to observer motion and we classify these terms as “observer corrections.” When removed (see, for example,⁵⁶), the transport equation becomes much simpler,

lacking the coupling between energy groups required by the observer corrections and more closely resembling the lab-frame formulation. However, this is not a lab-frame formulation, as lab-frame transport requires a more complex Lorentz-transformed lab-frame collision integral expressed in terms of lab-frame 4-momenta.

Agile-BOLTZTRAN³³ uses F = f/ρ for the transport variable in Lagrangian coordinates (t, m), so the “no observer corrections” transport equation (Eq. (2)) must be translated into an appropriate form:

If we multiply Eq (3) by E₀² and a constant and integrate over dE₀ and dμ₀, the ∂/∂μ₀ term will vanish and we have an expression for the evolution of neutrino number density. We can again integrate over the mass of the entire model and the ∂/∂m term will be reduced to a surface flux. Integration of the collision integral over mass gives the total lepton conversion rate (ν_e ↔ e⁻, etc.), which should balance with the integrated total change in neutrino number from the integrated time derivative term plus the surface flux. However, in this form we can see that this is not the case, there is a “compression” term, (F/ρ)(∂ρ/∂t), which does not integrate away. This term is part of the terms remaining in Eq. (2) after the observer corrections are removed with a counterpart “hidden” in the observer corrections that were eliminated. One consequence of this “compression” term is that when the fluid containing neutrinos is compressed, neutrinos will be artificially destroyed, and when the fluid expands neutrinos will be artificially created.

Lentz et al.⁵⁵ computed early phase 1D supernova models using Agile-BOLTZTRAN, starting from a full-featured, spherically symmetric supernova model and stripped away major features to approximate the variation among the multi-dimensional supernova codes. The starting point for these tests is a model (GR-FullOp, black lines, in Figures 1–3) that includes a full treatment of general relativity in the hydrodynamics and transport and up-to-date neutrino opacities that include energy-exchanging scatterings on nucleons and electrons, nucleon-nucleon bremsstrahlung, and the detailed LMSH nuclear electron capture (EC) table. First the GR effects, then many of the opacities, and finally the velocity-dependent ‘observer corrections' were removed to examine their importance. As the bounce profiles (Figure 1), shock trajectories (Figure 2), and neutrino luminosities and mean energies (Figure 3) all show, each of these broad changes has a noticeable impact on the conditions and evolution of the collapsed core. To consider the effect of ignoring the observer corrections, compare the N-ReducOp model (green) with the N-ReducOp-NOC model (blue) that uses the “no observer corrections” transport equation (Eq. (3)) in Figures 1–3. At bounce (Figure 1) the shock forms at much smaller mass without observer corrections (coincidentally matching the most complete GR-FullOp model) due to the sharp reduction in Y_e and Y_L to their lowest values at bounce. This is most starkly seen in the change from the largest core Y_ν (gap between dashed and solid lines in center lower panel of Figure 1) for the N-ReducOp model to the smallest Y_ν for N-ReducOp-NOC. This reduction in the Y_ν is a direct consequence of the “compression” term in Eq. (3), also visible in the total lepton number shown in Figure 2 (lower panel). The conserved lepton number is a integration over the entire model for net electron and electron neutrino number including the appropriate boundary flux term. The first three models show excellent total lepton conservation over the entire simulation, but the N-ReducOp-NOC model drops slowly just before bounce and sharply during shock breakout with continued loss during the post-bounce accretion phase. This non-conservation demonstrates the inconsistency and unphysical nature of Eq. (3) and likely causes the large reduction in peak

$Lνe$ and the ∼40% drop in neutrino luminosities during the accretion phase when multi-dimensional models should begin the shock revival process. The reduced luminosities likely reflect that some of the neutrinos are (non-physically) destroyed instead of escaping. There are other effects missing from the N-ReducOp-NOC model including Doppler shifts and angular aberration whose effects are more difficult to extract. In the N-ReducOp-NOC model the ⟨E_ν⟩_RMS is only slightly larger and the luminosity is significantly lower indicating that, if anything, the unfavorable nature of models without observer corrections will be even more unfavorable in multi-dimensional models where this reduction in neutrino luminosity may prevent the development of a successful explosion.

Clearly, from the failure to achieve explosions in spherically symmetric models that are as physically complete as any ever performed, we can judge that accurate neutrino transport, while an essential ingredient to models of core collapse supernovae, is not the only essential ingredient. However, the results of Lentz et al.⁵⁵ caution us that, as we seek less costly approximations to the neutrino transport suitable for multi-dimensional calculations of reasonable cost, care must be taken in many aspects of the treatment of neutrino transport.

Core-collapse supernova models that break the assumption of spherical symmetry have achieved more success than their spherically symmetric brethren. Instabilities in the PNS driven by the lepton and entropy gradients that result from radiating neutrinos enhance the neutrino luminosity in some models^{37,41,57–60} by increasing the rate that energy is transferred to the neutrinosphere, aiding the explosion. Large scale convection behind the shock originates from gradients in entropy that are born from the stalling of the shock and grows as the matter is heated from below. These fluid motions break the thermal stratification imposed by spherical symmetry and allow low entropy matter to descend close to the neutrinosphere, where it can be heated. This can enhance the efficiency of the neutrino heating, allowing simultaneous downward flows that fuel the neutrino luminosities by accretion onto the PNS and upward flows that bring energy to the shock and lead to explosions.^61,62 However, these fluid motions do not guarantee explosions.^20,21,41,63 Convection works in tandem with a computationally-discovered, multidimensional instability of the shock wave itself, the Standing Accretion Shock Instability (SASI).⁶⁴ This instability to non-radial perturbations favors low order modes, ultimately leading to gross distortions of the shock, which dramatically alters the shock and explosion dynamics.^47,64–67 Multidimensional simulations also bring with them the potential to accurately include the effects of rotation. Centrifugal effects in a rotating stellar core,^41,68 and other rotational effects, can change supernova dynamics quantitatively and perhaps qualitatively.

Recent axisymmetric (2D) simulations using spectral neutrino transport exhibit successful neutrino-driven explosions for a range of progenitors. The most physically complete of such simulations, using multigroup flux-limited diffusion^42,45,69 or a multigroup two moment method,^41,46,70 are the only simulations that include all of the physics that recent studies by Lentz et al.^55,71 find to be important in CCSN. The essential difference between the current exploding models^69,70 and similar models of the past decade,^41,72 which produced at best weak explosions for low mass progenitors, is the opening of the computational domain from 90° in latitude to 180°. This allows the SASI to excite the lowest order modes which work to push the shock outward gradually over several hundred milliseconds. The extended heating region behind the shock in turn allows increased neutrino heating to ultimately re-invigorate the shock and drive an explosion. In Bruenn et al.,⁴⁵ still strengthening explosions approaching the observed 1 Bethe (=10⁵¹ erg) are shown, half a second after the shock formed. Müller et al.⁴⁶ demonstrate qualitatively similar behavior, though with a much lower explosion energy. Three snapshots of hydrodynamic motion are visible in Figure 4, which shows the entropy (upper half panel) and radial velocity (lower half panel) for the 12 M_⊙ model from Bruenn et al.⁴⁵ at 150 ms, 300 ms, and 600 ms after bounce. At 150 ms, roughly 100 ms before rapid shock expansion heralds the onset of an explosion, asphericity is developing as a result of vigorous neutrino-driven convection and SASI. Bruenn et al.⁴⁵ find, as do Müller et al.,⁷³ that convection dominates over SASI in the development of explosion for lower mass progenitors with SASI playing the dominant role for more massive stars. By 300 ms the explosion is surging outward. However, low entropy down-flows still connect the unshocked regions with the PNS surface, continuing to supply accretion energy to power the neutrino luminosity and accelerate the explosion. By 600 ms, these down-flows have been cut off by the expanding ejecta, but their remnants continue to accrete onto the PNS, allowing the explosion to continue to gain in strength. This ongoing energizing of the explosion even after the shock has begun to move outward is a hallmark of self-consistent multidimensional simulations, but is not captured by parameterized light bulb models. That the source of the explosion, despite the importance of multidimensional fluid flow, is ultimately neutrino heating is driven home by comparison of recent self-consistent models with multigroup neutrino transport^45,46 to older self-consistent models with gray neutrino transport.^61,62 Compared to the past multi-dimensional simulations, all modern simulations exhibit a delay between core bounce and explosion that is several hundred milliseconds longer, the result of less efficient, but more realistic, heating in the spectral models.

While two dimensional axisymmetric simulations admit all of the important hydrodynamic effects, two dimensional hydrodynamics do not always respond in the same way as three dimensional hydrodynamics. For example, the turbulent cascade that proceeds to smaller scales in 3D actually proceeds to larger scales in 2D.⁷⁴ Thus the extent to which 2D supernova models represent 3D supernovae remains an open question. Fryer and Warren^68,75 demonstrated that three-dimensional models, at least when using gray transport, can exhibit large scale convective behavior similar to these two-dimensional models. Given the difference in the behavior of convection in 2D and 3D, this is a surprising result, possibly explained by the nature of the SASI.^64,76 A series of simulations in recent year, using variations on the neutrino light bulb approximation, but higher resolution than Fryer and Warren,^68,75 have investigated the differences between 2D and 3D supernova models. While Nordhaus et al.⁷⁷ find 3D simulations to be more favorable that 2D, producing explosions at lower neutrino luminosities, Hanke et al.⁷⁸ and Couch⁷⁹ find 3D to be at best similar to 2D, and likely more difficult. The answer to this question will ultimately come from 3D simulations with self-consistent spectral neutrino transport, the 3D analogues of the multigroup flux-limited diffusion or multigroup moments models discussed above.

Through advances in supercomputer hardware and supernova simulation software, such simulations are beginning to arrive. Fig. 5 shows a sample of the 3D fluid behavior from the early stages of a 3D simulations using CHIMERA, our MGFLD code. This equatorial slice shows that multidimensional fluid flow is clearly well developed by this point. Fig. 6 compares the entropy in a polar slice of this 3D model (at a slightly later point in time) to a 2D model based on the same progenitor. Together Fig. 5 and the upper panel of Fig. 6 show a roughly spherical shock stalled at approximately 200 km. This is a similar value to the equatorial shock radius shown in the lower panel of Fig. 6 at the same time. However the 2D shock surface in the lower panel of Fig. 6 is prolate, reaching roughly 300 km in the direction of both poles. Thus there is a significantly smaller shocked volume in the 3D simulation than in 2D simulation. Simulations by Hanke et al.,⁴⁷ using their multigroup moments code VERTEX, reach even further in time and show an even longer delay between core bounce and explosion than is evidenced in the axisymetric models. The next few years will see many more 3D simulations with accurate self-consistent neutrino transport, which will ultimately tell us the value of their axisymmetric analogues.

The PNS is an extremely dense, compact object, and its gravitational field is not well described by Newtonian gravity. General Relativity (GR) is required to properly capture the structure and dimensions of a neutron star. The release of gravitational binding energy as the core collapses and the PNS forms provides the energy in neutrinos, which ultimately powers the expulsion of the stellar envelope in the neutrino reheating paradigm. Thus, GR is fundamental to the energetics of core-collapse supernovae. In spherical symmetry, the formulation required to numerically solve the equations of General Relativity is straight-forward,^81,82 however equivalent formulations in multidimensions remained elusive until the relatively recent development of the BSSN formalism.^83–85 

Lentz et al.⁵⁵ examined the difference between models which include Einsteinian gravity (GR-FullOp) or Newtonian gravity (N-FullOp, red lines in Figs. 1–3) in both the hydrodynamics and transport. At bounce (Figure 1), the deeper gravitational well of the GR-FullOp model results in a more compact core. After bounce, the GR model's deeper gravitational well and more compact PNS result in a smaller gain region that inhibits heating and results in a smaller shock radius relative to the Newtonian equivalent.^38,39,41 The more compact PNS also leads to a larger accretion luminosity³⁸ and higher mean energy, ⟨E_ν⟩_RMS for all species of neutrinos, especially the ν_e and

$ν¯e$ that drive the primary heating that could revive the shock. On the basis of 1D models, Liebendörfer et al.³⁹ suggested that once the multi-dimensional nature of the mechanism lengthened the time matter spent in the gain region, the higher luminosity and harder spectrum in GR models would result in stronger or more viable explosions. This has been recently confirmed by Müller et al.,⁴⁶ whose GR models explode when their Newtonian and approximate GR equivilents do not. Of the multi-dimensional codes only VERTEX⁴⁶ now uses a GR formulation (the conformably flat approximation), while CHIMERA⁸⁶ and previously VERTEX^40,41 use approximate GR. Many multidimensional models,^43,44,56 remain Newtonian. To date, the most physically complete models to employ the BSSN formalism are leakage models by Ott et al.,⁸⁷ though the future will certainly see the increasing use of this approach coupled to more complete physics in the other sectors.

Pre-collapse stellar cores are excellent electrical conductors and likely harbor magnetic fields, sustained by some dynamo mechanism,⁸⁸ though the magnitude and topology are uncertain. We also know that neutron stars, the remnants of core-collapse supernovae, show evidence of extremely intense magnetic fields. In particular, magnetars^89,90 have inferred fields in the 10¹⁴-10¹⁵ G range, while “normal” pulsars have dipole fields of order 10¹² G.⁹¹ Beyond understanding the origin of pulsar magnetic fields, it is also natural to investigate what role magnetic fields can play in the core-collapse supernova explosion mechanism.

Indeed, not long after the discovery of pulsars, investigations into the potential role of magnetic fields in the explosion mechanism began using magnetohydrodynamic (MHD) models.^92–95 As shown by these early investigations, amplification by compression during core collapse and by rotation via the Ω-dynamo (continuing after bounce) can, in principle, give rise to rotationally-powered, magnetically-driven explosions. More recently, a renewed interest in the role of magnetic fields in the explosion of rotating progenitors has emerged,^96–105 motivated by observations of asymmetries in explosion ejecta,^96,106 natal neutron star kick velocities,¹⁰⁷ the association of Soft Gamma Repeaters and Anomalous X-ray Pulsars with magnetars,⁹⁰ and the theoretical discovery of the magnetorotational instability (MRI).¹⁰⁸ In particular, the MRI drives exponential field growth, as opposed to linear growth by the Ω-dynamo. These studies have shown that a sufficiently strong magnetic field, amplified by rapid rotation, can result in magneto-rotationally driven, and magnetically collimated, jet-like explosions. However, the requirement for rapid rotation (i.e. rotation rates of order milliseconds) limits the magneto-rotational mechanism to a subset of progenitor stars.

Some notable sources of uncertainty in the role of magnetic fields in core-collapse supernovae, common to most magneto-rotational collapse studies, include the nature (e.g., strength) of the pre-collapse magnetic field and the eventual impact of MRI-amplified magnetic fields. In particular, the wavelength of the fastest-growing MRI mode scales linearly with the magnetic field strength, leading to grid resolution requirements that are currently out of reach for global simulations starting with pre-collapse fields predicted by stellar evolution.¹⁰⁹ To partially overcome the resolution requirements needed to capture the MRI, the simulations have been initiated with (artificially) strong fields (e.g., of dipolar character). Another useful approach to studying the MRI involves local shearing box models, adopting post-bounce conditions.^110,111 These simulations have confirmed that the MRI can operate and does amplify magnetic fields exponentially. The shearing box simulations also show that the system develops into a turbulent state characterized by small-scale fields.

The MRI is a local instability, operating in the strongly differentially-rotating region above the PNS. At present, it is not clear how MRI amplified magnetic fields, as they may materialize in nature, resemble the fields obtained by collapsing strongly magnetized iron cores and how they may impact the global dynamics of a core-collapse supernova. This question must ultimately be addressed with global simulations (see, for example, the approach by Sawai et al.¹¹² in axial symmetry). However, as suggested by Masada et al.,¹¹¹ the magneto-rotational mechanism can be negatively affected by MRI-driven turbulence due to dissipation of the rotational free energy at the PNS surface. On the other hand, in the case of strong fields, energy deposition by resistive dissipation of magnetic energy may aid the neutrino-heating mechanism.^99,111

The lack of sufficient rotation in pre-collapse progenitor cores, as predicted by stellar evolution models,¹⁰⁹ has sparked recent interest in MHD processes in non-rotating CCSN environments; including Alfven wave-driven explosions,¹¹³ SASI-driven magnetic field amplification,^114–117 and magnetic amplification near an Alfven surface.¹¹⁸ Effects of magnetic resistivity on the explosion dynamics has also been investigated in the non-rotating case.¹¹⁹ As can be understood from simple considerations of the energetics, these studies show no significant influence of magnetic fields, unless the initial (i.e. pre-collapse) magnetic field is already relatively strong. For example, Endeve et al.,^115,116 using high-resolution three-dimensional MHD simulations, found that in the nonlinear phase of the SASI, turbulence develops via secondary fluid instabilities and feeds on the power in the dominant low-order SASI modes. The turbulence spreads throughput the shock volume and results in magnetic field amplification via a very efficient small-scale dynamo. This generates strong small-scale fields (see Figure 7) without significant effect on the shock dynamics (but possibly contributing to ultra-strong neutron star magnetic fields). It remains to be seen to what extent turbulence driven by neutrino-driven convection (in the PNS or the gain region) can play a similar role in amplifying magnetic fields (as predicted by Thompson & Duncan¹²⁰). Indeed, Obergaulinger & Janka¹¹⁷ carried out axisymmetric simulations of magnetized non-rotating iron cores that included neutrino transport. They observed magnetic field amplification from convection and the SASI in their models. Similar 3D simulations with sufficient spatial resolution will provide a better understanding of the role of MHD processes in core-collapse supernovae.

If not crucial to the explosion mechanism, MHD processes associated with the explosion dynamics may play an important role in magnetizing the PNS. Magnetic field amplification by the Ω-dynamo, MRI, SASI, and neutrino-driven convection (e.g., via an αΩ-dynamo) may all combine to provide a non-trivial contribution. Thus, a better understanding of the nature of pre-collapse magnetic fields is also needed to determine if PNS (and eventually pulsar) fields are due to properties of the progenitor or result from explosion dynamics.

Both the MRI and the SASI result in magnetic fields that evolve in a turbulent state, which is difficult to properly resolve with the shock capturing methods most commonly used for supernova simulation. (Endeve et al.¹¹⁵ estimated that late in their simulations the effective magnetic Reynolds number was of order unity.) Hence, one must use caution when extrapolating conclusions based on low Reynolds number simulations to systems governed by extremely large Reynolds numbers (such as the post-bounce supernova environment). Current simulations merely suggest that MHD turbulence (as opposed to hydrodynamic turbulence) may be the appropriate description in core-collapse supernovae, and a better understanding of its impact is clearly needed. Moreover, studying the confluence of the MRI, SASI, and neutrino-driven convection (which can excite an αΩ-dynamo¹²⁰) is needed to better understand the role of magnetic fields. As the available computing power increases and numerical algorithms improve, continued refinement of multiphysics simulations and simulations based on simplified models will be possible.

In general, an equation of state (EoS) is a thermodynamic relation between the pressure, temperature and internal energy, providing a closure to the hydrodynamic equations. For stellar matter, contributions from electrons/positrons, photons and atomic nuclei all must be considered. The thermodynamic equilibrium for strong and electromagnetic nuclear reactions which represents the final state of stellar nuclear fusion is termed Nuclear Statistical Equilibrium (NSE). For matter in NSE, the nuclear composition is a function of the thermodynamic state and electron fraction (Y_e, equal to ratio of protons to nucleons, and hence a measure of the neutronization), allowing the composition to be evolved as part of the equation of state, unlike the composition of the concentric layers of successively lighter elements that surround the iron core, where the evolution must be followed with a nuclear reaction network. At high temperatures, NSE favors free nucleons, with declining temperatures bringing increasing fractions of alpha particles and ultimately heavy nuclei. At low temperatures, the NSE composition is dominated by the maximally-bound nuclei from the vicinity of iron, giving rise to the iron peak observed in the solar abundance distribution. During collapse, the increasing density and neutronization causes NSE to favor heavier, more neutron-rich nuclei. This leads eventually to a composition dominated by exotic columnar and planar nuclear configurations (termed the nuclear pasta)¹²¹ and ultimately, when the core reaches densities similar to those of the nucleons in a nucleus, uniform nuclear matter. Ultimately, all of these phases of baryonic matter must be included in the supernova equation of state.

The composition and thermodynamic state of the core are clearly of profound impact to the hydrodynamic evolution and neutrino radiation transport in core-collapse supernovae. For example, the matter composition interacts with the neutrino opacities, determining how efficiently energy is transferred to the matter from the neutrino radiation field, while the matter pressure, derived from the thermodynamic state, causes the bounce and drives hydrodynamic motion. For these reasons, considerable effort has been invested by the nuclear astrophysics community in supernova equations of state. In addition to the traditional thermodynamic variables like pressure, internal energy and entropy, these CCSN EoSs provide nuclear compositional data assuming NSE, usually in the form of the mass fractions of protons, neutrons, α-particles and an average heavy nucleus, along with the atomic number and mass of this average heavy nucleus.

The equation of state of Lattimer & Swesty¹²² (LS) has long been the staple EoS for supernova simulations. It is based on a compressible liquid drop model, similar to that of Lamb et al.,¹²³ but in a lightweight form suitable for calculation within a supernova model. The nuclear composition for non-uniform nuclear matter is calculated in the Wigner-Seitz approximation of a heavy nucleus in a vapor of nucleons and α-particles. The LS EoS assumes a nuclear saturation density (ρ_s) of

$2.59×1014g cm −3$⁠, a symmetry energy (E_sym) of 29.3

$\mathrm{\thinspace MeV}$

$MeV$ and allows a selection of the nuclear incompressibility (κ_s). For consistency with earlier work, the simulations described in this and the following section use

$\kappa _s=180 \mathrm{\thinspace MeV}$

$κs=180 MeV$⁠, but simulations using a more realistic

$\kappa _s=220 \mathrm{\thinspace MeV}$

$κs=220 MeV$⁠, including those in the preceding sections, are similar.^124,125

The equation of state of Shen, Toki, Oyamatsu & Sumiyoshi¹²⁶ (STOS) is based on relativistic mean field theory. It is constructed assuming

$ρs=2.43×1014g cm −3$⁠,

$E_{sym}= 36.9 \mathrm{\thinspace MeV}$

$Esym=36.9 MeV$ and

$\kappa _s=281 \mathrm{\thinspace MeV}$

$κs=281 MeV$⁠. Thomas-Fermi calculations¹²⁷ provide the nuclear composition for non-uniform nuclear matter, also using the Wigner-Seitz approximation.

The equation of state of Wilson and his collaborators is described most fully in Wilson and Mathews¹²⁸ (see also^51,129,130). For supranuclear matter, the empirical prescription of Baron et al.²⁶ is used with

$ρs=2.66×1014g cm −3$⁠,

$\kappa _s= 200 \mathrm{\thinspace MeV}$

$κs=200 MeV$ and the supranuclear adiabatic index, γ_s, = 2.75. These choices for ρ_s, κ_s and γ_s, as well as the form of the symmetry energy

$E_{sym} = 16 (1-2 \mathrm{Y}_{e})^2 (\rho /\rho _s )(1 + 4.5/(1 + 4\rho /\rho _s)) \mathrm{\thinspace MeV}$

$Esym=16(1−2Ye)2(ρ/ρs)(1+4.5/(1+4ρ/ρs)) MeV$⁠, were informed by relativistic Brueckner-Hartree-Fock calculations.¹³¹ The Wilson EoS also includes the effects of pion production at high density, using the model of Friedman et al.¹³² This model is constrained by comparision between experimental measurements and simulations of pion production in heavy ion collisions.¹³⁰ The identity of the average heavy nucleus is calculated from a density and electron fraction dependent analytic formula.¹²⁸ 

In addition to these stalwarts, recent years have seen a renaissance in nuclear equations of state valid for core-collapse supernovae, which, unlike the cold neutron star case, require the effects of non-zero temperature. These include new EoS by Hempel and Schaffner-Bielich,¹³³ based on an ensemble of nuclei in nuclear statistical equilibrium (NSE), and Shen et al.¹³⁴ and Shen et al.,¹³⁵ based on Hartree approaches in which the density profiles are self-consistently calculated. All of these approaches to the nuclear physics of CCSN matter have different benefits and drawbacks.

In order to illustrate the impact that the EoS can have on a core-collapse supernova simulation, we compared three EoS; LS, STOS and Wilson. Full simulations were run with each of the three EoSs using Agile-BOLTZTRAN,³⁵ as part of the M.S. thesis of M. Baird. These simulations include General Relativity, which is important for a nuclear EoS comparison since the inclusion of GR results in roughly 30% higher central densities. GR simulations therefore exercise the supranuclear portion of the EoS to a greater extent than their Newtonian counterparts. To facilitate comparison with the results of Hix et al.,¹³⁶ neutrino opacities consistent with those simulation were used here, those from Bruenn³⁶ with the Langanke-Martìnez-Pindeo-Sampaio-Hix (LMSH) nuclear electron prescription.^136,137 The conditions at bounce for the three EoSs as a function of enclosed mass are shown in Figure 8. The Wilson EoS results in the largest homologous core (0.58 M_⊙), much larger than either the STOS (0.52 M_⊙) or LS (0.48 M_⊙) cases. For much of the collapse, Y_e is higher for the Wilson cases than the others, producing this larger core. At bounce, this effect is still apparent in the outer part of the inner core (between 0.2 and 0.5 M_⊙), though a sharp gradient in Y_e develops in the Wilson EoS case over the last few ms before bounce in the regions with density above

$2×1014g cm −3$⁠. Here the pion contribution becomes significant in this EoS. This decrease in the Y_e is caused by the impact of the pions on the free nucleon chemical potentials and is consistent with previous results using this EoS.¹²⁹ The differences in Y_e at lower density between the three EoSs result from differences in the composition. For example, the generally lower heavy nucleus mass fraction of the Wilson EoS results in less electron capture on heavy nuclei consistently throughout the star. Similarly, the largest difference in heavy nucleus mass fractions between the STOS and LS EoS results in a noticeable rise in Y_e around 0.8 M_⊙ in the STOS case.

The time history of the shock over the first 200 ms after bounce is shown in Figure 9. The small scale noise seen beyond 10 ms results from the determination of the location of the shock in Agile-BOLTZTRAN. In simulations with the LS EoS (⁠

$κs=180 MeV$⁠), the shock is launched, and reaches a radius of approximately 189 km (enclosed mass of 1.38 M_⊙) before stalling 130 ms after bounce. In contrast, the shock in simulations using the STOS only reaches a distance of approximately 165 km (1.32 M_⊙) from the center before stalling. This maximum is reached much sooner (60 ms after bounce) than in the LS EoS model. With

$\kappa _s=180 \mathrm{\thinspace MeV}$

$κs=180 MeV$⁠, LS EoS is a relatively soft EoS and, as a consequence, stores a larger amount of energy in the PNS than a stiffer EoS, allowing the shock to travel to larger radius before stalling due to energy loses from neutrino emission and the disassociation of heavy nuclei. The STOS EoS, with

$\kappa _s=281 \mathrm{\thinspace MeV}$

$κs=281 MeV$⁠, is stiffer than the LS EoS, so while the initial pressures are higher and thus the early progress of the shock more rapid, less energy is stored in the compression of the PNS, despite the initial PNS mass in the STOS case being 8% larger than in the LS case. As a result the shock can not travel as far out from the center before stalling. The shock in simulations with the Wilson EoS travels out to a distance of approximately 209 km (1.44 M_⊙) before stalling,

$20 \mathrm{\thinspace km}$

$20 km$ further than the LS EoS case, and

$44 \mathrm{\thinspace km}$

$44 km$ further than the STOS EoS. Even though the Wilson EoS, with

$\kappa _s=200 \mathrm{\thinspace MeV}$

$κs=200 MeV$⁠, is nominally stiffer than the LS EoS, it provides more energy to the shock than does the LS EoS, likely as a consequence of both the high density pion contribution and the larger initial size of the PNS that results from reduced electron capture. Even 200 ms after bounce, by which time convective instabilities will have developed in the neutrino-heated region of the supernova (or multi-dimensional simulations thereof), the shock's location, which marks the outer boundary of this unstable region, exhibits a significant (∼20%) variation due to the EoS. This is a clear indication that the EoS is important to the physics of core collapse and bounce and its impact can extend even to times well after bounce.

The impact of nuclear physics on core-collapse supernovae after the formation of the supernova shock is limited by the shock's dissociation of heavy nuclei into free nucleons. Beyond the core collapse phase, only deep in the core, where the nucleons form nuclear matter, and in the outer regions where the shock is not strong enough to fully dissociate nuclei, does nuclear physics play a large role. The composition of the iron core of the star, as dictated by Nuclear Statistical Equilibrium, as well as the pressure, entropy and other thermodynamic quantities calculated within the Equation of State, are based on the temperature, baryon density and electron fraction.

Evolution of the electron fraction is determined not just by the nuclear composition, as we discussed in the previous section, but also by the reaction rates for processes like electron and neutrino capture that can alter the neutronization of the matter. Bethe et al.¹³⁸ pointed out that, due to the low entropy of the stellar core and resulting dominance of heavy nuclei over free nucleons, electron capture processes on heavy nuclei would dominate the evolution of the electron fraction during the late stages of stellar evolution and stellar core collapse. In the iron core, this predominantly occurs via Gamow-Teller (GT) transitions changing protons in the f_7/2 level into neutrons in the f_5/2 level. However, it was soon realized that as core collapse proceeds, increasing neutronization and density lead to average neutron numbers >40, filling the neutron f_5/2 orbital and quenching further electron capture on heavy nuclei. Calculations using the independent particle model (IPM) showed that neither thermal excitations nor forbidden transitions substantially alleviated this blocking.^139,140 However, it is well known that the residual nuclear interaction (beyond the IPM) mixes the fp and gds shells, for example, making the closed g_9/2 shell a magic number in stable nuclei (N = 50) rather than the closed fp shell (N = 40). Full shell model diagonalization calculations remain impossible in this regime due to the large number of available levels in the combined fp + gds system.¹⁴¹ Langanke et al.¹⁴² developed a “hybrid” scheme, employing Shell Model Monte Carlo (SMMC) calculations of the temperature-dependent occupation of the various single-particle orbitals to serve as input to Random Phase Approximation (RPA) calculations for allowed and forbidden transitions to calculate the capture rate. This approach enabled Langanke et al.¹³⁷ to demonstrate that electron capture on heavy nuclei dominates the capture on protons throughout core collapse. While direct measurement of electron capture on these radioactive species is prohibitive, charge exchange reactions are being used to test the theoretical Gamow-Teller strength distribution underlying these rate calculations.¹⁴³ Juodagalvis et al.¹⁴⁴ have recently updated the compilation of theoretical nuclear electron capture rates needed for core-collapse supernovae.

As a result of the continuation of nuclear electron capture throughout collapse, Hix et al.¹³⁶ demonstrated that modern electron capture rates markedly reduced (∼10%) the electron fraction in the interior of the PNS, resulting in an ∼20% reduction in the initial mass of the PNS, for a 15 solar mass progenitor using the LS EoS. This effect is complementary to those discussed in the preceding section, as the total electron capture rate depends on both the electron capture rates per nucleus and the nuclear abundances provided by the EoS. Figure 10 demonstrates that the effect of the transition from IPM nuclear electron capture rates³⁶ to modern rates is much smaller for the Wilson EoS than for the LS EoS. While in the LS EoS case the initial mass of the PNS is reduced by nearly 0.09 M_⊙, in the Wilson EoS case, it is reduced by slightly more than 0.01 M_⊙. This is an excellent example of the interplay between different pieces of microscopic physics in the core-collapse supernova problem. Were the compositions provided by the as-yet-unknown true core-collapse supernova EoS to be more similar to those of Wilson than LS or STOS, the effects of modern nuclear electron capture could be more easily neglected. Lentz et al.⁷¹ have recently shown that the softer neutrino spectrum from modern nuclear electron capture also fills the low-energy neutrino spectrum directly without requiring neutrino–electron scattering, with we will discuss further in the following section.

As we discussed earlier, Lentz et al.⁵⁵ computed early phase spherically symmetric supernova models using Agile-BOLTZTRAN, by starting from a full-featured supernova model and stripped away major features to approximate the variation among the multi-dimensional supernova codes. Starting with a model with full GR hydrodynamics and transport and a modern neutrino opacity set, first the GR effects, then many of the opacities, and finally the velocity-dependent ‘observer corrections' were removed. As we can see from the bounce profiles (Figure 1), shock trajectories (Figure 2), and neutrino luminosities and mean energies (Figure 3) each of these broad changes has a noticeable impact on the conditions and evolution of the collapsed core. Electron (positron) capture on nuclei in the collapsing core and on free protons (neutrons) in the shocked core are the dominant sources of cooling and deleptonization by emission and heating through the inverse absorption processes. Though scattering on the nucleons is the dominant source of total opacity and therefore controls the trapping of neutrinos, scattering on electrons was long thought to also be important as it preferentially scatters neutrinos to lower energies where the total opacity is lower and allows more neutrinos to escape during collapse.¹⁴⁵ 

To demonstrate the broad effects of modern opacities, the FullOp opacities used in the first two models (emission, absorption, and scattering on free nucleons;¹⁴⁶ isoenergetic scattering on α-particles and heavy nuclei;³⁶ scattering of neutrinos on electrons (NES) and positrons (NPS);¹⁴⁷ production of neutrino pairs from e⁺e⁻ annihilation¹⁴⁷ and nucleon-nucleon bremsstrahlung;¹⁴⁸ and electron capture (EC) on nuclei using the LMSH EC table^137,141) are simplified to ReducOp. In ReducOp the non-isoenergetic scattering (NIS) on nucleons is replaced by isoenergetic scattering (IS)³⁶ as are nucleon emission and scattering rates. To remove the rest of the NIS opacities, NES and NPS were both omitted entirely as their total opacity is small and primarily contribute through energy downscattering. Finally, the LMSH EC table was replaced by IPM nuclear EC³⁶ to complete the simplified ReducOp opacity set. The broad effect of this change in opacity can be seen in the difference between the N-FullOp (red) model and N-ReducOp (green) models in Figures 1–3. Reduced electron capture and deleptonization from the simplified EC and the lack of NIS energy downscattering leads to higher core Y_e and Y_L = Y_e + Y_ν. The larger number of trapped neutrinos Y_ν in N-ReducOp is evident in the center lower panel of Figure 1 as the large gap between Y_e(solid lines) and Y_L (dashed lines). The higher core Y_L leads to the formation of the shock at a larger mass coordinate (shallower in the gravitational well) and a vigorous initial shock with “ringing”³² (Figure 2). The IPM EC not only effectively removes all nuclear electron capture at high densities, but overestimates the EC at lower densities. This leads to more rapid collapse and higher density in the outer core (see near 0.9 M_⊙ in Figure 1) but also a steeper density gradient that slows down shock expansion. These reduced opacities together result in a lower shock trajectory that is less favorable to explosion.

To follow-up the broad effects of opacity simplification, Lentz et al.⁷¹ reexamined most of the opacities individually using the full set as a reference point. When individual opacities were removed from a full set the effect was often different, frequently weaker, than seen in adding the same opacity improvement to a less complete set such as Bruenn.³⁶ The most striking of these was the effect of neutrino-electron scattering (NES) on the deleptonization during collapse. Previously, neutrinos down scattered in energy by NES had entered a lower opacity ‘window' through which they could more easily escape and enhance the loss of lepton number from the collapsing core.¹⁴⁵ The earlier test¹⁴⁵ had been done using the IPM nuclear EC that provides most of the neutrinos during collapse, but when NES was removed using the LMSH nuclear EC, downscattering was not necessary as neutrinos are emitted directly at lower energies and there was no discernible difference in the simulations at bounce. Despite the decreased impact of NES on the dynamics of the simulation, omitting NES still altered the ν_e-luminosity and mean energy that would affect detection and nucleosynthesis. All of the opacities tested had effects on the luminosities, spectra, and/or dynamics with the exception of neutrino-positron scattering. Overall, this study⁷¹ illustrated that for testing the effects of included or missing neutrino opacities that the context of the other included opacities is important and demonstrates the importance of most of the opacities commonly used in full interaction sets.

A fundamental reason we study core-collapse supernovae is that large overabundances of elements in the periodic table spanning from oxygen through nickel are observed in core-collapse supernovae and their remnants. In the innermost regions of the ejecta, the passage of the shock heats matter to temperatures where Nuclear Statistical Equilibrium (NSE) is dominated by free nucleons and α particles.^149,150 As a result, most of the iron-peak species synthesized in core-collapse supernovae result from α-rich freezeout.¹⁵¹ As matter expands outward, it cools, allowing the light nuclei to recombine into iron, nickel, and neighboring nuclei. In the case of α-rich freezeout, this recombination is incomplete, leaving a significant fraction of the matter still in the form of free nucleons and α particles. The detailed composition of this ejecta depends on its neutronization, which is set in the inner regions of the ejecta by neutrino interactions. Parameterized spherically symmetric models^152–154 find that above the innermost α, iron, and nickel dominated regions, passage of the shock leaves a layer rich in the α isotopes: ⁴⁰Ca, ³⁶Ar, ³²S, and ²⁸Si — the products of incomplete silicon burning. Above this is a layer of ¹⁶O, in the outer portions of which significant fractions of ²⁰Ne, ²⁴Mg, and ¹²C are found. Finally, above the ¹⁶O layer we find the helium layer and hydrogen envelope, if they were not driven off as part of a stellar wind.

Observations of elemental abundances and the distribution of these newly-formed species allow nucleosynthesis calculations to place powerful constraints on conditions deep in the interior of supernovae and their progenitors, places hidden from direct observation. The temperatures to which the supernova shock heats the ejecta make ultraviolet observations a very useful way of studying the enriched composition of recent ejecta.^155–162 With missions like Chandra, Swift and NuSTAR, such elemental observations can be extended to include the X-ray region.^163–165 Furthermore, X-ray observations of supernova remnants, like the observations of Cassiopeia A^166,167 and SN1987a,¹⁶⁸ also provide detailed data on the elemental composition and distribution of heavy elements in supernova ejecta, data which must be explained by explosion models. Infrared observations provide another view of heavy element abundances and their distribution in supernovae^169,170 and supernova remnants^171–173 and future spectroscopic missions like the James Webb Space Telescope will report on the temporal history of element production from supernovae.

This wealth of observations reveals a consistent pattern of extensive, large-scale mixing of the chemical elements.¹⁷⁴ For example, Mazzali et al.¹⁷⁵ infer from observations that the Type Ic supernova SN 2002ap had an oxygen-rich inner core with ⁵⁶Ni at higher velocity, a strong indication of large-scale overturn. The case is particularly persuasive for supernova SN 1987A. Observed asymmetries in iron lines are most easily explained by the concentration of iron-peak elements into high-velocity “bullets”.¹⁷⁶ The Bochum event, the rapid development of fine structure in the H_α line from SN1987a roughly 2 weeks after the explosion,¹⁷⁷ was interpreted by Utrobin et al.¹⁷⁸ as an indication that a large (∼10⁻³M_⊙) clump of nickel was ejected at high velocity

$(∼4700 kms −1)$ into the far hemisphere of the supernova. Similarly, near-IR observations of He I lines arising roughly 2 months after explosion were interpreted by Fassia and Meikle¹⁷⁹ as indications of dense clumps of nickel mixed into the hydrogen envelope of SN 1987a. Indications of large-scale but clumpy elemental mixing extend to optical,¹⁵⁸ infrared^172,173 and X-ray¹⁶⁶ observations of supernova remnants. Rather than expanding chemically distinct spheres of ejecta, the picture of a supernova that emerges is one of metal-rich shrapnel tearing though the envelope of the now deceased star.

Observations of γ-ray lines are an excellent way to determine not just the elemental but isotopic production of supernovae. They are therefore of particular interest, both to supernova theory, in general, and to predictions of supernova nucleosynthesis, in particular. Observations like the Compton Gamma-Ray Observatory (CGRO) measurement of the ratio of ⁵⁷Ni to ⁵⁶Ni in SN 1987A¹⁸⁰ or the detection of ⁴⁴Ti in the Vela¹⁸¹ and Cassiopeia A¹⁸² SN remnants probe the deepest layers of the supernova ejecta, just above the mass cut.¹⁸³ A further source of isotopic data comes from presolar grains that come to Earth as part of primitive meteorites. (See, e.g.,^184–186 for reviews of the properties of presolar grains and the methods of their analysis). Grains with an excess in the abundance ratio of ¹⁸O/¹⁶O are ascribed a supernova origin¹⁸⁷ and many of these supernova grains show evidence for extinct radioactivities.^188,189 Only the neutrino and gravitational wave signals, which follow the evolution of the proto-neutron star, provide information from deeper within the explosion than γ-ray lines and premolar grains. Except for the few cases where the progenitor has been observed prior to the supernova,¹⁹⁰ elemental and especially isotopic composition of the ejecta allow the most direct determination of the mass and metallicity of the progenitor.^191,192

This nucleosynthesis involves a large number of nuclear reactions, strong and electromagnetic reactions as well as weak nuclear reactions involving capture of electrons, positrons and neutrinos occurring under conditions where nuclear statistical equilibrium does not apply. To examine the non-equilibrium evolution of the chemical composition that is a vital part of core-collapse supernova requires some form of nuclear reaction network within any simulation of core-collapse supernovae. For many models which include neutrino transport, whose goal is to understand the CCSN mechanism, this is often a small α-network (linking the 13 or 14 nuclear species of even and equal neutron and proton numbers from helium to nickel or zinc). In some such cases, an even simpler 3 level “flashing” scheme is used, which instantaneously transmutes oxygen to silicon to iron (or NSE) at appropriate temperatures.

Unfortunately, much of our current understanding of core-collapse supernova nucleosynthesis is based on parameterized models which replace the inner workings of the supernova with a kinetic energy piston^152,154,193 or a thermal energy bomb.^153,194 These parametrized models ignore much that we have learned in the past two decades about the nature of these neutrino-driven, hydrodynamically-unstable explosions. In these bomb or piston simulations, the explosion's energy, its delay time and the mass cut, which separates the ejecta from matter destined to become part of the neutron star, are externally supplied parameters. Even many recent simulations^80,195 employ such methods. Aufderheide et al.¹⁹⁶ demonstrated that the bomb and piston parameterizations are largely compatible, with the largest differences coming in the inner regions of the ejecta. It is this inner region that is also most strongly affected by the details of the explosion mechanism.¹⁹⁷ In the case of the neutrino reheating mechanism, these effects include interaction with the tremendous flux of neutrinos, and the temporal delay in achieving the explosion. This provides strong motivation to merge modeling of the nucleosynthesis with modeling of the supernova central engine, a task now being undertaken by several groups.

One common property exhibited by all recent one- and two-dimensional simulations utilizing spectral neutrino transport^{39–41,69,70} is a decrease in the neutronization in the outer part of the neutrino heating region, due to neutrino interactions. This is a feature that the parameterized bomb/piston nucleosynthesis models discussed above cannot replicate because they ignore the neutrino reheating mechanism. The neutronization is important because Galactic Chemical Evolution (GCE) calculations and the relative neutron-poverty of terrestrial iron and neighboring elements strongly limit the amount of neutronized material that may be ejected into the interstellar medium by core-collapse supernovae.¹⁹⁸ Hoffman et al.¹⁹⁹ placed a limit of 10⁻⁴M_⊙ on the typical amount of neutron-rich (Y_e < 0.47) ejecta allowed from each core-collapse supernova. Past multidimensional models of core-collapse supernovae using gray (energy-integrated or -averaged, i.e., not spectral) neutrino transport that did produce explosions tended to greatly exceed this limit.^20,61 To compensate, modelers have been forced to invoke the fallback of a considerable amount of matter onto the neutron star, occurring on a time scale longer than was simulated.²⁰⁰ 

Work by Fröhlich et al.^201,202 and Pruet et al.^203,204 indicates that neutrino-powered explosions using spectral neutrino transport result in nucleosynthesis products qualitatively different in composition from either the parameterized bomb/piston nucleosynthesis models or gray transport models of the core-collapse mechanism. In the models used by both Fröhlich et al. and Pruett et al., the neutrino physics was artificially tuned in order to drive explosions. Thus, these models are also parameterized, but the explosion is explicitly neutrino-driven and the ejecta are determined self-consistently by the radiation hydrodynamics. These models have shown that the inclusion of neutrino captures on the ejecta, which decrease the neutronization, remedies some defects in the predictions made by previous models of nucleosynthesis, all of which neglected this important piece of physics. These interactions remove the over-production of neutron-rich iron and nickel isotopes that have long plagued parameterized bomb and piston models. These simulations^201,203 also show enhanced production of Sc, Cu, and Zn; elements which observations of metal-poor stars²⁰⁵ suggest are 3–10 times more abundant than previous models predicted. Moreover, Fröhlich et al.²⁰² revealed a significant neutrino-driven flow to proton-rich nuclei above A = 64, which was also confirmed by Pruet et al.,²⁰⁴ leading to the discovery of the νp-process and suggesting that the innermost ejecta of core-collapse supernovae may be the production site of the light p-process nuclei. The discovery of a potential new nucleosynthesis process requires significant amounts of additional nuclear data be measured or calculated²⁰⁶ for accurate simulations. Given the importance of these findings, it is essential these studies now be conducted in the context of self-consistent, multidimensional models using spectral neutrino transport.

The separation between models of the central engine and the nucleosynthesis can in part be traced to two factors. First, the frequent failure of self-consistent models to produce explosions in the past limited their application to the study of supernova nucleosynthesis. Second, while the central engine plays its role in initiating the explosion in perhaps 1-2 seconds and nucleosynthesis is largely complete within a similar timeframe, the spatial and velocity distribution of the newly-made isotopes continues to develop over many minutes as the shock propagates through the rest of the star and matter falls back onto the PNS. Now that self-consistent models are producing realistic explosions, efforts are underway to remove this division. This requires including the tools to compute realistic nucleosynthesis within self-consistent models of CCSN which include neutrino transport and extending the models, in some fashion, to much later epochs than the 1-2 seconds such models typically cover. These tools include progressively more realistic nuclear reaction networks, growing from 3 species flashing schemes to small α-networks toward reaction networks which include the hundreds of species at play in CCSN, and tracer particles to allow post-processing analysis of nucleosynthesis with nuclear reaction networks larger than those included within the self-consistent simulation. In spite of the challenges it poses, truly understanding the nucleosynthesis in CCSN must start from the central engine and model the complete destruction of the star.

The post-processing approach is widely used in supernova nucleosynthesis calculations. In this approach, a simplified approximation is used to calculate the compositional change and resulting rate of nuclear energy generation within the hydrodynamic model. Supernova mechanism models typically use an α-network or an even simpler 3 level “flashing” scheme. Subsequently, the temperature and density histories for individual mass elements from the supernova model are then used as input for separate, larger nucleosynthesis calculations. Performing post-processing calculations based on one-dimensional models is relatively straightforward, since most one-dimensional simulations are Lagrangian. Thus the needed temporal histories of temperature and density are simply those of the individual Lagrangian mass elements. However, in multi-dimensions, fluid motions can result in highly tangled Lagrangian grids. As a result, Eulerian hydrodynamics, where the discretization occurs in space rather than mass, is used to perform most multi-dimensional stellar astrophysics simulations. Because Eulerian codes use spatial discretization, the Lagrangian thermodynamic histories which are a natural result in a Lagrangian code are unavailable. Instead passive tracer particles must be employed.

Convection and hydrodynamic instabilities greatly complicate the picture developed from spherically symmetric models, potentially destroying strict compositional layering in the progenitor. Meakin and Arnett²⁰⁷ (see also^208,209) have demonstrated that convection occurs in the oxygen shell even before the explosion. For two decades, it has been known that Rayleigh-Taylor instabilities originate at the Si/O and (C+O)/He boundaries,^{194,210–213} however these instabilities do not mix nickel to sufficiently high velocities to account for observations of Supernova 1987A²¹⁴ (and references therein) and other core-collapse supernovae. The implication is that gross asymmetries must be present in the core and be part of the mechanism itself. However few simulations to date have directly considered the impact of this multi-dimensional behavior on the nucleosynthesis. Nagataki et al.¹⁹⁴ and Maeda et al.²¹⁵ employed parameterized thermal bomb models which demonstrated that in aspherical models a larger fraction of the ejecta experiences α-rich freezeout. Unfortunately, these simulations tracked the composition only via post processing, i.e. no evolution of the composition was included within the hydrodynamic simulation. Kifonidis et al.^216,217 considered the impact of the neutrinos and the effects of the convective overturn in 2D simulations that reached beyond the central engine. They used a “neutrino lightbulb” similar to Janka and Müller,²⁰ coupled to a 13 species α-network, while separately tracking the neutronization. These simulations showed significantly higher velocities than previous models.^210–212 However, as a result of crude neutrino transport, the simulations of Kifonidis et al.²¹⁶ still predict the ejection of much larger quantities of neutron-rich iron group elements than Galactic chemical evolution seems to allow. These hydrodynamic models have recently been expanded to three dimensions by Hammer et al.,²¹⁸ starting from the asymmetric parameterized neutrino-driven explosions modeled by Scheck.^219,220 The simulations of Hammer et al.²¹⁸ are hampered in their ability to study the nucleosynthesis in detail, since the initial models of Scheck²¹⁹ did not consider nuclear burning, but they do show that in 3D, metal-rich clumps suffer much less deceleration at the shell interfaces. This results in asymptotic velocities of metal-rich clumps that are twice or thrice those found in 2D models, better matching observations. 3D late time supernova simulations by Young et al.,^221,222 which extended a spherically symmetric initial explosion model, also reveal dense knots in the ejecta.

The complex fluid flow that results from convection and the SASI witnessed in multi-dimensional models makes separating the ejecta from material trapped in the proto-neutron star more difficult. In the parameterized spherically symmetric models that dominate our understanding of supernova nucleosynthesis, the mass cut, a Lagrangian coordinate defining the inner edge of the ejecta, is a parameter determined by comparison to limits on the ejected mass of ⁵⁶Ni and other constraints. As Figure 11, (taken from a 2D simulation using our CHIMERA supernova code⁶⁹) illustrates, convection and the SASI complicate this picture, with some parcels of material with the same initial radius being ejected while others are captured in the neutron star. Clearly in multi-dimensional models, and likely in Nature, the mass-cut loses its strictly spherical shape. It may even fail to be simply connected in the topological sense, with regions of eventual ejecta being immersed in material destined to form part of the neutron star. Limits on the neutron-richness of the iron and nickel ejecta in the neutrino-less parameterized nucleosynthesis models generally place the mass-cut in the oxygen layer. However, multi-dimensional simulations are beginning to reveal greater complexity. For example, Figure 12 (taken from the same CHIMERA simulation as Figure 11) shows that some of the matter that began deep in the silicon layer (red circles) of a 12 solar mass star seems likely to be ejected. Because of its exposure to the neutrino field, this matter is likely to be less neutron-rich than it was when the star began to collapse, easing the constraints on the neutron-richness of the ejecta.

As a result of their gray neutrino transport, the simulations of Kifonidis et al.²¹⁶ predict the ejection of much larger quantities of neutron-rich iron group elements than Galactic Chemical Evolution (GCE) seems to allow. In contrast, preliminary analysis of our spectral 2D CHIMERA results, based on analytic extension of tracer histories to late times,^223,224 document the development of the proton-rich ejecta and the νp-process, though at a somewhat reduced level compared to previous work,²⁰² because of faster expansion time scales. This highlights the critical dependence of the strength of the νp-process on the expansion time scale. Recent extended simulations^225–227 highlight the impact (or lack thereof) of the forming reverse shock on the expansion time scale and hence the strength of the νp-process, providing further justification for the merger of explosion mechanism models and nucleosynthesis models, as well as the extension of these models in time until the nucleosynthesis is complete. Such merged simulations are necessary to fully understand the quantity and distribution of intermediate-mass and iron-peak nuclei made in core-collapse supernovae, likely resolving the existing issue of too-neutron-rich ejecta, and helping us understand the role of massive stars in GCE. These simulations will also likely reveal the location and strength of the νp-process, indicating the range of p-process species whose cosmic abundance also result from the deaths of massive stars.

Advancement in our understanding of core-collapse supernovae, both the mechanism that drives the explosion and the nucleosynthesis that results, is driven simultaneously by improvements in our ability to model the macroscopic phenomena like hydrodynamics and radiation transport with greater physical fidelity and improvements in our understanding of the microscopic physics that drives these events. Over the past decade, improvements in computational technology and our ability to harness it have revealed the importance of spectral neutrino transport, reinforced our appreciation for the convective engine that drives these events and furthered the development of a numerical formalism that enables General Relativity to be included in multidimensional core-collapse supernova simulations. Matching this greatly improved physical fidelity at macroscopic scales has been a continued improvement in our understanding of the physics at microscopic scales; the nuclear equation of state, nuclear electron capture and other neutrino opacities. Here too, improvements in computational technology, and our ability to harness it, constrained where possible by experimental measurements, have improved our understanding of core-collapse supernova physics.

In the next decade, we will see three dimensional models that finally include all of the essential physics that our current models suggest is necessary to understand core-collapse supernovae; multidimensional spectral neutrino transport coupled to magnetohydrodynamics and General Relativity. These models will include ever improving nuclear physics and neutrino opacities and realistic nuclear reaction networks. These simulations will definitively answer the question of how well axisymetric simulations can reproduce three dimensional reality and, we may hope, produce explosions that mimic the topology observed in supernovae and supernova remnants. They will let us understand the development of pulsar magnetic fields and the potential these fields have to alter the dynamics of the explosion. They will show unequivocally how nucleosynthesis is impacted by the convectively unstable, neutrino-driven explosion and illuminate the contribution that core-collapse supernovae make to the p and r processes. They will also open the door to physically faithful simulations of more exotic events; failed supernovae that explode as collapsars producing γ-ray bursts, the formation of magnetars and the mergers of neutron stars. Achieving these goals will require the expenditure of considerable human effort and computational power, and thus the continued support of the nuclear and astrophysical communities.

This work has been supported by grants from the NASA Astrophysics Theory Program (NNH11AQ72I) and the NSF PetaApps, Nuclear Physics and Stellar Astrophysics Programs (OCI-0749242,AST-0653376) and by the Department of Energy Offices of Nuclear Physics and Advanced Scientific Computing Research. Computational resources were provided by the NSF's TeraGrid at the National Institute for Computational Sciences under grant number TG-MCA08X010; by the National Energy Research Scientific Computing Center, supported by the U.S. DoE Office of Science under Contract No. DE-AC02-05CH11231; and by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program at the Oak Ridge Leadership Computing Facility, supported by the U.S. DoE Office of Science under Contract No. DE-AC05-00OR22725.