Long-Term Multidimensional Models of Core-Collapse Supernovae: Progress and Challenges

TL;DR

The review surveys progress and open challenges in long-term, multidimensional core-collapse supernova simulations, emphasizing the neutrino-driven mechanism as the leading explanation for most explosions while acknowledging modeling differences across codes. It traces the field from early 2D efforts to mature 3D modeling with multigroup neutrino transport (Vertex/RbR+ and Nemesis) that now tracks the explosion through shock breakout and into the remnant phase, including detailed neutrino and gravitational-wave predictions. Key insights include the impact of GR gravity and many-body corrections on explodability, the role of progenitor perturbations, and the emergence of diverse nucleosynthesis and remnant properties (kicks, spins, Ni/titanium production) in long-term 3D runs. The article also highlights ongoing uncertainties in neutrino flavor conversion, the conditions leading to black-hole formation, and the need for broader progenitor and EoS studies, along with the importance of comparing simulations with SN observations (e.g., SN 1987A, Cas A) to constrain the theory of neutrino-driven explosions.

Abstract

Self-consistent, multidimensional core-collapse supernova (SN) simulations, especially in 3D, have achieved tremendous progress over the past 10 years. They are now able to follow the entire evolution from core collapse through bounce, neutrino-triggered shock revival, shock breakout at the stellar surface to the electromagnetic SN outburst and the subsequent SN remnant phase. Thus they provide general support for the neutrino-driven explosion mechanism by reproducing observed SN energies, neutron-star (NS) kicks, and diagnostically relevant radioactive isotope yields; they allow to predict neutrino and gravitational-wave signals for many seconds of proto-NS cooling; they confirm correlations between explosion and progenitor or remnant properties already expected from previous spherically symmetric (1D) and 2D models; and they carve out various scenarios for stellar-mass black-hole (BH) formation. Despite these successes it is currently unclear which stars explode or form BHs, because different modeling approaches disagree and suggest the possible importance of the 3D nature of the progenitors and of magnetic fields. The role of neutrino flavor conversion in SN cores still needs to be better understood, the nuclear equation of state including potential phase transitions implies major uncertainties, the SN 1987A neutrino measurements raise new puzzles, and tracing a possible correlation of NS spins and kicks requires still more refined SN simulations.

Figures (8)

• Figure 1: Equation-of-state dependence of CCSN explosions of an 18.88 $M_\odot$ star (R. Bollig, private communication, and References Bollig+2021Janka+2023; Kresse et al., in preparation). All 3D SN simulations were started from a progenitor model with 3D perturbations of velocity, density, and chemical composition due to convective oxygen-shell burning Yadav+2020 and employed different nuclear EoSs widely used in SN simulations: LS220 Lattimer+1991, SFHo, SFHx Steiner+2013Hempel+2010, DD2 Typel+2010Hempel+2012, and APR Schneider+2019b. Left panel: Spherically averaged shock radius (solid lines) and PNS radius (defined at a density of $10^{11}$ g cm$^{-3}$; dashed lines) as functions of post-bounce time. Right panel: Time evolution of diagnostic explosion energy (solid lines) and explosion energy with overburden energy of the progenitor taken into account (dashed lines; for a definition of these energies, see Section \ref{['sec:lessons']}). In all cases successful neutrino-driven explosions were obtained, but there is a clear correlation between the onset of the explosion and the contraction of the PNS. The faster the PNS contracts, the earlier the explosion sets in. This signals a crucial and sensitive influence of the PNS's radius evolution on the start of the explosion. Note that the PNS radii in the simulations with the APR and SFHx EoSs nearly overlap, and the shocks and explosion energies in both cases show similar behavior. Figure courtesy of Daniel Kresse and Robert Bollig.
• Figure 2: "Explodability" as function of the ZAMS mass of progenitor stars in large sets of CC simulations using different "neutrino engines" in 1D calculations and different neutrino-hydrodynamics codes in 2D models to obtain neutrino-driven explosions (Glas et al., in preparation). Circles mark BH forming cases and crosses indicate successful explosions. Upper four panels: Onset time of the explosion (defined by the post-bounce time when the average shock radius exceeds 300 km) versus ZAMS mass with averages indicated by horizontal lines. The results of the Push 1D engine models ( top left) were taken from Ebinger+2019, where progenitors from Woosley+2002WoosleyHeger2007 were considered; the P-Hotb 1D engine models ( middle left) are from Sukhbold+2016 using progenitors from Sukhbold+2014WoosleyHeger2015; and the Fornax 2D simulations VartanyanBurrows2023 as well as the Vertex 2D simulations (R. Bollig, private communication; middle right) employed the SFHo EoS Steiner+2013Hempel+2010 and progenitors for $\ge$12 $M_\odot$ from Sukhbold+2018 and for $<$12 $M_\odot$ from WoosleyHeger2015 and, in the case of Vertex, also from Woosley+2002. Both Push and P-Hotb triggered explosions artificially in 1D with neutrino engines whose parameters were calibrated by reproducing explosion properties of SN 1987A. Push yields 79 (81%) explosions and 18 (19%) BH cases out of 97 models, P-Hotb 90 (57%) explosions and 67 (43%) BHs in 157 models, Fornax 63 (63%) explosions and 37 BH cases (37%) of 100 models, and Prometheus-Vertex 73 (41%) explosions and 104 BH cases (59%) of 177 models. Bottom panels: Core compactness $\xi_{1.75}$ versus ZAMS mass for the Fornax and Vertex 2D simulations with successful and failed explosions indicated by different symbols. Gray horizontal lines mark $\xi_{1.75} = 0$ to guide the eye. Figure courtesy of Robert Glas.
• Figure 3: PNS convection in MD SN simulations compared to MLT convection in 1D PNS cooling calculations with the Prometheus-Vertex code (Heinlein et al., in preparation). Left panels: Radial profiles of matter entropy and electron fraction at different post-bounce times (denoted in the box above the panels) for a 1D PNS model with a baryonic mass of 1.44 $M_\odot$ and a 3D simulation of a PNS with a baryonic mass growing from 1.48 $M_\odot$ at 0.5 s to 1.55 $M_\odot$ at 5.0 s. The PNS is defined by a mass-density $\rho \ge 10^{11}$ g cm$^{-3}$ and contracts from initially 30 km to 14 km. Small differences between the 1D and the 3D model near the surface are a consequence of continuous accretion in the 3D case. The convective layer can be recognized by flat entropy profiles. Right panels: Lab-frame luminosities ( top) and mean energies ( bottom) of radiated $\nu_e$, $\bar{\nu}_e$, and a single species of heavy-lepton neutrinos, $\nu_x$, as functions of time for a 1D PNS cooling model with MLT convection compared to a 2D simulation that was chosen because its constant baryonic PNS mass of 1.36 $M_\odot$ closely matches that of the 1D model. Note that the luminosities of $\nu_e$ and $\nu_x$ are shifted by constant factors to avoid cluttering of the different lines. The 2D neutrino data are angle-averaged. MLT convection is able to reproduce the results obtained by MD simulations with good accuarcy; here it is important to know that there are no appreciable differences between PNS convection in 2D and 3D. Figure courtesy of Malte Heinlein.
• Figure 4: Long-term 3D CCSN simulations of the Garching group for 9.0, 9.6, 12.28, 15, 18.88, and 20 $M_\odot$ progenitors (Stockinger+2020Bollig+2021Janka+2024; Kresse et al., in preparation; models s12.28b, m15a, m15b, s18.88b, s20a, and s20b correspond to models s12.28, m15, m15e, s18.88, s20, and s20e, respectively, already discussed in Reference Janka+2024). The calculations were performed with the Prometheus hydrodynamics code including Vertex neutrino transport for periods between 0.5 s and 5.1 s and continued with the Nemesis neutrino treatment to later times. The m15 models are based on a rotating progenitor Summa+2018, and the s12.28 and s18.88 models made use of progenitors whose convective oxygen-shell burning had been computed in 3D for the final hour before CC in the 12.28 $M_\odot$ case and the final 10 min in the 18.88 $M_\odot$ case Yadav+2020. The CCSN runs were either performed with the LS220 EoS Lattimer+1991 or the SFHo EoS Steiner+2013Hempel+2010Janka+2024. Top left: Spherically averaged shock radii as functions of post-bounce time. Bottom left: PNS baryonic masses (solid lines) and gravitational masses (dashed lines) as functions of post-bounce time. Top right: Explosion energies versus post-bounce time with overburden energies (i.e., the binding energies of the progenitor stars above the outward moving shocks) taken into account (for a definition, see Section \ref{['sec:lessons']}). The a and b versions of the models differ slightly in their explosion energies (b versions being more energetic) and roughly bracket uncertainties connected to neutrino heating and cooling as well as hydrodynamic stochasticity. Bottom right: Time evolution of the PNS kicks due to asymmetric mass ejection (dashed lines) and asymmetric mass ejection plus anisotropic neutrino radiation (solid lines). Figure courtesy of Daniel Kresse.
• Figure 5: Mass downflow and outflow properties, neutrino heating rates, and explosion energetics in the long-term 3D CCSN simulations of 9.0, 9.6, 12.28, 15, 18.88, and 20 $M_\odot$ progenitors shown in Figure \ref{['fig:3Dmodels']} (Bollig+2021Kresse2023; Kresse et al., in preparation). Upper row: Mass downflow rates ( left), mass outflow rates ( middle), and ratios of outflow rate to inflow rate as functions of post-bounce time ( right), all measured at a radius of 400 km. Inflow and outflow rates are effectively equal between $\sim$1 s and $\sim$(6--8) s, because inflowing matter absorbs energy from neutrinos close to the PNS and is re-ejected afterwards as long as neutrino heating is still strong enough. Afterwards $\dot M_\mathrm{out}/\dot M_\mathrm{in}$ begins to decrease, thus signaling the onset of fallback accretion by the PNS. Models s9.0 and z9.6 are exceptions as their mass inflows to the PNS stop early, giving way to NDW outflows, which, however, are much weaker than the inflow-outflow rates of all other models. Model m15a is also an exception with an outflow-to-downflow ratio smaller than unity, signaling continuous mass accretion by the PNS (see lower left panel of Figure \ref{['fig:3Dmodels']}). Lower row: Total (i.e., volume-integrated) neutrino-heating rates ( left), ratios of the growth rate of the explosion energy (with overburden taken into account) to the total neutrino-heating rate ( middle), and ratios of the growth rate of the explosion energy ($\dot E$ with overburden taken into account) to the mass outflow rate at 400 km ( right). The energy deposited by neutrino heating around the PNS leads to a continuous, long-lasting growth of the explosion energy (see upper right panel of Figure \ref{['fig:3Dmodels']}). Models s9.0 and z9.6 are again exceptions, because their neutrino heating achieves to gravitationally unbind near-surface PNS matter in the NDWs, however without providing any relevant net gain for the explosion energies. In phases where the growth rate of the explosion energy exceeds the total rate of neutrino energy deposition, $\dot E/\dot Q_\nu > 1$, nuclear energy generation contributes to the growth of the explosion energy. This typically happens at very late times ($\mathrel{\vcenter {\hbox{$>$}\hbox{$\sim$}}}$6 s after bounce) except in models m15a and m15b, where the heating rates drop more rapidly in comparison. Figure courtesy of Daniel Kresse.