The Physics of Core-Collapse Supernovae: Explosion Mechanism and Explosive Nucleosynthesis

Abstract

Recent developments in multi-dimensional simulations of core-collapse supernovae have considerably improved our understanding of this complex phenomenon. In addition to that, one-dimensional (1D) studies have been employed to study the explosion mechanism and its causal connection to the pre-collapse structure of the star, as well as to explore the vast parameter space of supernovae. Nonetheless, many uncertainties still affect the late stages of the evolution of massive stars, their collapse, and the subsequent shock propagation. In this review, we will briefly summarize the state-of-the-art of both 1D and 3D simulations and how they can be employed to study the evolution of massive stars, supernova explosions, and shock propagation, focusing on the uncertainties that affect each of these phases. Finally, we will illustrate the typical nucleosynthesis products that emerge from the explosion.

1. Introduction

Massive stars (${\mathrm{M}}_{\mathrm{ZAMS}}$ > 8–10 ${\mathrm{M}}_{\odot }$) end their lives with the formation and subsequent collapse of a Fe core. Over the years, the multiple efforts of the stellar community have produced evolutionary codes able to accurately simulate the life of these objects [1,2,3,4,5,6,7,8,9,10]. However, the calculation of the complete evolution of a massive star, from H burning up to the formation of the Fe core is still characterized by several physical and computational challenges (e.g., [11,12,13,14]). Mixing and transport processes, magnetic fields, mass-loss, and nuclear reaction rates are only some of the many uncertainties that affect the evolution of massive stars.

The subsequent phase, i.e., a core-collapse supernova (CCSN), gives rise to new sources of uncertainties. The wide range of densities (${10}^3$–${10}^{15}$ g ${\mathrm{cm}}^{-3}$), temperatures (1–100 MeV), and electron fractions ($0.01$–$0.6$) that characterize a supernova pose an extremely challenging problem. To accurately model a CCSN, one needs a detailed description of the properties of matter at extremely high densities and temperatures, i.e., the equation of state (EOS) of nuclear matter, which is however still poorly constrained (e.g., [15,16,17,18,19,20,21,22,23,24]). As a consequence of the high densities reached deep inside the proto-neutron star (PNS), general relativistic effects can play a non-negligible role [25]. Moreover, at such high densities, the mean free path of neutrinos is very small, and therefore they are trapped inside the PNS and are only able to escape once the density goes below ∼${10}^{11}$–${10}^{12}$ g ${\mathrm{cm}}^{-3}$. Since a huge amount of neutrinos is produced, and they are responsible for triggering the explosion, one needs an accurate description of the interactions and transport of those neutrinos within the star [26,27,28]. Finally, self-consistent, high-resolution simulations are required to capture multi-dimensional effects that play a crucial role in the explosion [29,30,31,32,33,34,35,36,37,38,39].

CCSNe are also responsible for creating half of the elements of the periodic table. Given the extreme and diverse thermodynamic conditions achieved in the explosion, several nucleosynthetic processes may occur, as, e.g., the explosive burning stages relative to the major hydrostatic fuels of the stars (Si, O, Ne, C, He), $\nu$p–process, (weak) r–process, or $\gamma$–process. Most of these are currently performed in post-process, i.e., the nuclear network calculations are computed on the output of hydrodynamic simulations, as opposed to being computed out on-the-fly, which is much more computationally expensive and unfeasible with current codes, although promising efforts are underway [14,40]. These nucleosynthesis calculations are therefore heavily dependent on the characteristics of the explosion, mainly the explosion energy and the mass cut, and on the uncertainties in nuclear reaction rates.

In this review, we will summarize the pre-supernova evolution and uncertainties in the stellar evolution. Then, we will review the current state-of-the-art of CCSNe, with a particular focus on the nuclear EOS and the explodability problem. Finally, we provide a general overview of explosive nucleosynthesis.

2. The Progenitors of CCSNe: Massive Stars and Their Evolution

The evolution of massive stars is governed by a sequence of nuclear burning stages, in which lighter nuclear species are progressively transformed into heavier ones. Each burning stage primarily takes place in the core of the star and then it shifts outward in mass (shell burning) as soon as the available fuel diminishes. Normally, shell burning takes place at higher temperatures and lower densities in comparison to central burning. In each successive burning phase, the main product from the previous phase becomes the primary fuel, until the formation of a dense core made predominantly of Fe-peak elements (iron core). The energy released through nuclear reactions plays a crucial role in counterbalancing the gravitational collapse, especially in the advanced stages of stellar evolution, where neutrino losses also come into play. In the following, we briefly discuss the key aspects of the evolution that lead massive stars to explode as a CCSN.

2.1. Pre-Supernova Evolution

Six major fuels identify the principal burning stages in massive stars, i.e., H, ⁴He, ¹²C, ²⁰Ne, ¹⁶O, and ²⁸Si [41].

In the core of massive stars, H fusion occurs through the CNO cycle, while the proton-proton (p-p) chain plays only a marginal role [42,43]. The conversion of H in He lasts a few million years and scales inversely with the initial mass of the star. H burning in massive stars is characterized by the presence of a convective core which recedes in mass by decreasing the amount of available H. Once the H is depleted in the core, the star contracts until the He burning reactions, i.e., the $3\alpha$ reaction and ¹²C($\alpha$,$\gamma$)¹⁶O, are activated, while the H burning shifts in a shell. In this phase, the radius expands, the surface cools down, and the star becomes a red supergiant (RSG). Moreover, mass-loss may critically increase and play a crucial role, since it may reduce the size of the H-depleted core (He core) and therefore modify the evolutionary path of the star. Core He burning occurs in a convective core as well and it is controlled by the $3\alpha$ reaction in the early stages, which produces ¹²C. As ¹²C increases, the ¹²C($\alpha$,$\gamma$)¹⁶O reaction gradually becomes more and more efficient, until it takes over the $3\alpha$ reaction. The competition between these two nuclear processes controls the C/O ratio at the time of the He exhaustion in the core.

At the end of central He burning, two main quantities govern the further evolution of the star: the ¹²C/¹⁶O ratio and the mass of the CO core left by the He depletion. The CO core mass plays the same role the initial mass has in H burning. The ¹²C abundance left by central He burning influences both the formation of the convective core during the central C burning phase and the behavior of the C burning shells during the central burning of Ne, O, and Si, contributing to determine the final compactness of the star and eventually its final fate (see Section 4.2). Neutrino energy losses begin to become relevant at a temperature of the order of ∼800 MK, i.e., the typical temperature of C ignition, and enter into the total energy balance together with gravitational and nuclear energy. From this point on, the evolution of the star can be described by a sequence of progressively faster contraction, central and shell burning stages, that lead ¹²C to be converted into ²⁰Ne, ²⁰Ne into ¹⁶O, ¹⁶O into ²⁸Si, and finally ²⁸Si into Fe-peak elements.

It is however interesting to note that the processes that lead from the exhausted O core to the formation of the iron core deviate from the classical definition of nuclear burning. Instead, they can be interpreted as a re-adjustment of the chemical composition after the nuclear reactions have reached a (partial) equilibrium. For a detailed discussion of the latest stages of massive stars, we refer the reader to more specialized works (e.g., [2,8,44,45]).

Figure 1 shows the distribution of the abundances of the major fuels in a non-rotating 20 ${\mathrm{M}}_{\odot }$ star at solar metallicity, calculated with two different stellar evolutionary codes, namely FRANEC [7,8,46], and MESA [47,48,49]. As discussed in the next section, different codes can produce significantly different results for the same set of initial conditions (e.g., mass, metallicity, rotation, magnetic field…).

2.2. Challenges and Uncertainties in Stellar Evolution

While significant progress has been made in modeling the evolution of massive stars, several challenges and uncertainties still persist. In the following, we mention and briefly discuss the main open problems affecting the calculation of CCSN models. Note that this list does not claim to be exhaustive and only provides a starting point for the reader to frame the open problems in the modeling of stellar evolution.

2.2.1. Mass-Loss

Accurate modeling of mass-loss in massive stars remains a challenge. The mechanisms responsible for mass-loss, such as radiation-driven winds and eruptions, are not fully understood and they can have a critical impact on the evolution, the nucleosynthesis, and on the final fate of massive stars (e.g., [50,51]). The mass-loss prescriptions in stellar models are often empirical and rely on arbitrary parameters; moreover, their dependency on the metallicity of the star is still unclear.

2.2.2. Mixing and Transport Processes

The most significant source of uncertainty in stellar models is undoubtedly represented by the treatment of convection, which is intrinsically a multi-dimensional mixing and transport process. It is impossible to determine the convective transport, the extension of the convective zones, and the velocity of the eddies of material by first principles only. In order to include its (inherently multi-dimensional) effects in 1D, spherically symmetric codes, it is commonly modeled by means of the mixing length theory (MLT). Broadly speaking, the larger the magnitude of the mixing length, the stronger the convection. The mixing length is usually expressed as:

$${\Lambda }_{\mathrm{MLT}}={\alpha }_{\mathrm{MLT}}\frac{P}{\rho g},$$

(1)

where ${\alpha }_{\mathrm{MLT}}$ is a parameter of order ∼$\mathcal{O}\left(1\right)$, and P, $\rho$, and g are the local pressure, density, and gravitational acceleration, respectively. Typically, the value of ${\alpha }_{\mathrm{MLT}}$ is calibrated to reproduce the properties of the Sun [41,52,53]. There is a huge amount of literature on this topic and definitions of the mixing length beyond the simple expression given in Equation (1) have been proposed, as well as models beyond MLT calibrated on multi-dimensional simulations [54]. However, this goes beyond the scope of this article, and we point the reader to dedicated reviews on this topic [55].

Another intrinsically multi-dimensional phenomenon in stars is rotation, which is responsible for the transport of angular momentum and chemical species within the star. Rotation is usually parameterized in 1D and calibrated to reproduce the observed properties of a sample of nearby rotating stars, such as the surface chemical enrichment and the equatorial velocity [7,8,56,57].

These two processes play a significant role in determining the distribution of elements within the star, impacting subsequent nucleosynthesis, the composition of the stellar envelope, and eventually even enhancing the mass lost from the surface of the star. In particular, they highlight the need to develop multi-dimensional stellar codes that can adequately treat the transport mechanisms inside a star.

2.2.3. Magnetic Fields

The role of magnetic fields in massive star evolution is not well-constrained. Much like mass-loss, convection, and rotation, it depends on arbitrary parameters. Modeling stellar magnetic fields involves solving magneto-hydrodynamic (MHD) equations, which describe the interactions between magnetic fields and fluid motion. Accurate numerical solutions are computationally intensive and require advanced algorithms [58]. Therefore, as it is common practice in studying the evolution of massive stars, we will ignore the effects of magnetic fields in the remainder of this paper.

2.2.4. Nuclear Reaction Rates

Nuclear reaction rates are responsible for the energy generation and the nucleosynthesis of all the nuclear species in every layer of the star. Large uncertainties affect these reaction rates, which propagate, in turn, on the calculation of massive star models (e.g., [59,60]). More precise experimental data and improved theoretical models are needed to better constrain these rates.

2.2.5. Binary Nature of Massive Stars

The majority of massive stars are in binary systems, and interactions with a companion star can have a significant impact on their evolution [50,61,62,63,64]. Understanding the effects of binary interactions, such as mass transfer, common envelope evolution, and mergers, is a current challenge that in recent years is becoming more and more approachable with modern tools [65]. For the remainder of this paper, we will ignore binary interactions and focus solely on the evolution of single stars.

2.2.6. Core-Collapse Supernova Mechanism

The core-collapse supernova mechanism, responsible for the explosive death of massive stars, is still not fully understood. The details of the shock revival, neutrino interactions, and explosion dynamics are complex and require advanced and multi-dimensional numerical simulations. We will discuss this topic in detail in Section 3.

2.2.7. Stellar Codes

Performing accurate and computationally intensive simulations of massive star evolution, including hydrodynamics, radiation transport, and nuclear reactions, remains a challenge. Moreover, the adoption of different criteria for mass-loss, convection, rotation, magnetic fields, as well as of different physics inputs and different numerical methods, leads inevitably to a degeneration of solutions for the calculation of a model using the same initial conditions (see, e.g., [66,67] and Figure 1). Regarding this last point, one key aspect is the nuclear network adopted, namely, the list of all the isotopes explicitly included in the calculation and of all the reactions that link them with each other. Recent studies (e.g., [68]) showed in fact that large nuclear networks are required especially in the latest stages of the evolution of massive stars, where multiple reactions can contribute to the energy generation that sustains the structure and regulate the behavior of burning shells. Underestimating the number of nuclear species in the nuclear network may substantially decrease the computation time and facilitate the convergence of models, but it may also lead to significantly different results relative to the calculations made using larger nuclear networks. For this reason, increasing computational capabilities is an ongoing effort to produce more refined and accurate simulations. Moreover, as mentioned above, a multi-dimensional approach is required, but still far from being realized.

3. Core-Collapse Supernova Theory

The first theoretical models of supernovae date back to the late fifties [69,70]. In these seminal papers, it was postulated that a supernova explosion could be powered by thermonuclear burning of the material right outside the core. This was later shown to be incorrect by the first numerical simulations [71,72], which instead recognized the crucial role of neutrinos as being the primary cause of the explosion.

In these simulations, neutrinos emitted during the collapse would deposit energy behind the shock, energizing the explosion right after core bounce. Later simulations [73] that employed more detailed microphysics, showed that this so-called “prompt neutrino-driven mechanism” could not provide enough energy to power an explosion. Instead, they found that a “delayed neutrino-heating”, after the initial expansion and a brief stalling of the shock, could be responsible for reviving the shock at later times and launching the supernova.

Broadly speaking, this is a well-accepted explosion mechanism even today. However, significant theoretical and computational efforts have shown that a supernova is an extremely complicated interplay between microphysics [28,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88], radiation transport, and magneto-hydrodynamics [26,27,89,90,91,92], where multi-dimensional effects play a crucial role [29,30,31,32,33,34,35,36,37,38,39,93,94].

Nonetheless, it is of significant pedagogical value to illustrate the physical mechanisms behind a simplified model of the explosion.

3.1. The Collapse Phase

After a massive star has built up its iron core, nuclear reactions turn off, since iron has the largest binding energy among all nuclear species. This removes the pressure support generated by nuclear energy that prevents the core from collapsing. Therefore, gravity becomes the dominant force and this triggers core-collapse. At this point, electrons are degenerate, and neutrino emission is the main process responsible for the cooling of the iron core. As the core collapses, densities and temperatures increase, which lifts electron degeneracy, and therefore facilitates electron capture on free protons and nuclei, which deleptonizes the core. As a consequence, the electron fraction goes from ∼0.42–0.45 down to ∼0.3. Neutrino emission tends to decrease the entropy, whereas electron captures on nuclei and photodissociation tend to increase it [95]. The net result is that entropy is roughly constant, and the collapse is therefore adiabatic. Together with the extremely low (practically zero) pressure which fails to balance gravity, this causes the collapse to be homologous, i.e., $v\propto r$, where v and r are the velocity and radius, respectively. This means that there is some radius $r_{\mathrm{s}}$ (the sonic point) at which the speed of the collapse is greater than the sound speed. Beyond that point, pressure waves cannot propagate fast enough to re-equilibrate the changing pressure gradients, and therefore matter is in approximate free fall.

The central density increases extremely rapidly and at some point the neutrinos that are being emitted become trapped, at density ${\rho }_{\mathrm{trap}}\gtrsim {10}^{11}\mathrm{g}{\mathrm{cm}}^{-3}$. Therefore, a neutrinosphere is formed, i.e., the radius at which neutrinos are not trapped anymore, and therefore start free-streaming outwards. Notice that, technically speaking, the radius of the neutrinosphere is a function of neutrino energy (higher-energy neutrinos decouple at larger radii), although oftentimes in the literature the term neutrinosphere is used to indicate the energy-averaged neutrinosphere, for simplicity. Moreover, different species will decouple at different radii since they experience different interactions with matter. Specifically, for a given energy, the heavy-lepton neutrinos and antineutrinos decouple at smaller radii than electron antineutrinos, which in turn decouple at smaller radii than electron neutrinos. Given the key role that neutrinos play in the explosion, the transition between the completely trapped regime and the free-streaming regime has to be accurately simulated. This requires solving the Boltzmann transport equation, and several numerical approaches can be chosen, as we will briefly mention in the following sections.

The collapse is halted only when the inner core reaches nuclear saturation densities (∼$2.5\times {10}^{14}{\mathrm{g}\mathrm{cm}}^{-3}$), at which point the nuclear equation of state stiffens due to the strong force becoming repulsive. The collapse of the inner core stops, and a pressure wave starts propagating outwards. Once this pressure wave reaches the sonic point, it steepens into a shock wave near the edge of the homologous core, and this is what is called the “bounce”.

3.2. The Bounce Phase

Once the shock wave is launched, the inner core becomes causally disconnected from the outer layers of the star, and its electron fraction remains practically frozen at a value of ∼0.3.

The shock propagates outwards and photodissociates the infalling Fe-core material into free nucleons and alpha particles. This quickly drains the kinetic energy of the shock and, after a brief period of positive post-shock velocity [96], the shock stagnates, and moves outwards in radius simply due to the huge accretion rate that settles matter onto the newly-born Proto-Neutron Star (PNS), thus increasing the post-shock pressure in the form of neutrino heating. The shock reaches a maximum stalling radius of ∼150–200 km after 100–200 ms from the bounce. Now, it will either slowly recede back following the contraction of the PNS, producing a failed SN, or it will be re-energized by some heating mechanism, producing a successful explosion.

3.3. Neutrino Delayed Heating

Once the shock has reached its maximum stalling radius of ∼150–200 km, matter in the post-shock region is emitting neutrinos as a cooling mechanism but is also absorbing neutrinos escaping the PNS. The net effect of these two competing processes is that matter is cooling in the region closer to the PNS, due to neutrino emission being stronger. However, in the region closer to the shock, matter is being heated, since neutrino absorption is stronger than emission. This is the so-called gain region. Therefore, the shock loses energy due to photodissociation of the infalling material, and it gains energy due to neutrino heating. However, all modern spherically symmetric simulations [97] have shown that, without any multidimensional effects, the stalling shock slowly recedes quasi-statically. Eventually, a large fraction of the star will fall back on the central compact object, leaving behind a black hole. Exceptions to this are extremely low massive stars with steep density profiles [98], or the addition of some exotic physics, such as a phase transition to quark matter at high densities [99].

After the neutrino-delayed heating mechanism was established, several studies investigated whether an increased amount of neutrinos emitted from the PNS (mostly due to convective motions below the neutrinosphere) could increase the heating in the gain region enough to trigger an explosion [100,101,102]. However, this was later found to be incorrect [103,104], although PNS convection is still important in several aspects of CCSNe [105,106,107,108,109].

3.4. Multi-Dimensional Effects

It is now well-established that, to get a successful explosion, neutrino heating must be complemented by another source of energy that is inherently multi-dimensional. Rotation and magneto-hydrodynamic instabilities have been shown to aid the explosion, although only in a small fraction of CCSNe, due to the large angular momentum and magnetic fields required [38,110,111,112]. Another mechanism that has been proposed is the “jittering jets” mechanism [113], which has however not been observed so far in self-consistent simulations. For most CCSNe, the processes that can aid the explosion are mainly the Standing Accretion Shock Instability (SASI) [114,115,116], and neutrino-driven turbulent convection [31,93,94,117,118].

SASI is a large-scale instability (characterized by $l=1,2$) operating as an advective-acoustic cycle that can amplify shock expansion. Neutrino-driven turbulent convection refers instead to the convective instabilities that originate as a consequence of the negative entropy gradient in the gain region. These two mechanisms are sometimes hard to distinguish numerically (especially in 2D) since they both lead to similar effects on large scales [119], but have nonetheless very distinct features, the main one being that SASI is a quasi-periodic phenomenon, leading to quasi-periodic oscillations of the shock.

Since these multi-dimensional phenomena have such an important role in the explosion, several groups have studied how asphericities in the pre-SN star can affect either SASI or neutrino-driven turbulent convection (or both). Generally, they found that density perturbations caused by turbulent motions in the latest stages (i.e., a few minutes) before core-collapse can be accreted through the shock and act as seeds for turbulent motions in the gain region, and therefore facilitate the explosion [120,121,122,123,124].

3.5. Microphysics: The Nuclear EOS and Neutrino Physics

The bounce and post-bounce phases of the supernova, as well as the deleptonization history during the collapse phase, can dramatically change the following dynamics. Therefore, uncertainties in the nuclear equation of state and neutrino interactions can have a significant impact on the outcomes of a supernova.

As it turns out, the Equation of State (EOS) for nuclear matter is not yet well constrained (e.g., [15,16,17,18,22,23,24]) since the very large densities reached in the PNS are extremely hard to probe. Therefore, several different EOSs satisfy current experimental, observational, and theoretical constraints, despite predicting very different outcomes when it comes to supernova explosions. Given the huge impact that the EOS has on CCSNe (and any other astrophysical phenomena involving compact objects), countless studies have been carried out in the last few decades [74,78,96,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140].

In general, soft EOSs tend to facilitate the explosion, whereas stiff EOSs tend to disfavor the explosion. Qualitatively, a soft EOS is characterized by a slow increase of pressure as density increases, whereas for a stiff EOS the pressure increases more rapidly with density. Therefore, for a soft EOS, the PNS contracts to smaller radii since the pressure response to the gravitational pull is small. Therefore, it can reach higher temperatures, which produce more energetic neutrinos, increasing the neutrino heating in the gain region. However, due to the highly non-linear nature of CCSNe, as well as the complicated dependence of neutrino opacities on the EOS, this simple picture is not necessarily realized in simulations. In other words, stiffness is not necessarily correlated with stronger shock revivals and explosions. Nonetheless, there have been some recent papers [137,139] that were able to point out clear correlations between some specific property of the EOS (specifically the effective nucleon mass) and the strength of the explosion. However, this correlation is only valid for EOSs calculated using Skyrme interactions theories. Since the definition (and calculation) of the effective nucleon mass depends on the theoretical framework used to perform the calculation (e.g., Skyrme interaction theories versus relativistic mean-field theories), this correlation breaks down when considering a broader range of EOSs. Interestingly, another correlation can still be observed: the central entropy in the PNS right after bounce correlates with the strength of the subsequent explosion [140]. No other thermodynamic quantity or neutrino property (i.e., luminosity, energy, neutrinosphere properties, etc.) shows any correlation whatsoever.

To illustrate the wide range of possible outcomes that can be obtained by varying the EOS, we show in Figure 2 the evolution of the shock radius for a 20 ${\mathrm{M}}_{\odot }$ progenitor for several different EOSs, alongside the mass-radius relations of cold neutron stars allowed by those EOSs. Some of them satisfy all of the current constraints, whereas others have been recently ruled out, but have been widely used in simulations up to the very recent past. It is clear that, especially given the bifurcation nature of the CCSN problem, changes in the EOS can give rise to completely different supernova dynamics and, in certain cases, outcomes.

Another significant source of uncertainty is neutrino physics. Specifically, neutrino transport and neutrino opacities can both largely affect the pre and post-bounce dynamics. Significant progress has recently been made on neutrino opacities [81,86,87,149,150,151,152] as well as neutrino transport [27,81,89,90,153,154,155,156,157,158]. More recent, and significantly harder to interpret in the context of CCSNe, is the role of collective neutrino oscillations [159,160,161].

Finally, it is worth mentioning the important role that General Relativity plays in CCSNe. Due to the high densities achieved in the PNS, general relativistic effects can be quite significant. Until recently, the vast majority of simulations would solve the Newtonian equations for radiation hydrodynamics, with multipole corrections to the Newtonian potential à la Marek et al. [162], and only a few multi-dimensional simulations would solve the full GR equations [29,109]. However, current code capabilities and improvements of physical prescriptions implemented, have reached a high enough accuracy that general relativistic effects must be fully accounted for. Some parametric 1D simulations [25] have already shown that GR might have non-linear effects not reproducible by approximate methods, and therefore full GR simulations will become more and more relevant in the near future.

Due to the complexity of all of these subjects, we point the reader to dedicated reviews [19,20,26,28,159,160,161] and references therein.

3.6. State-of-the-Art CCSNe Codes

The field of CCSNe has recently seen a promising increase in the number of codes sophisticated enough to perform self-consistent simulations. Some of the most common codes employed to perform these kinds of simulations (although the list is likely longer) are:

• FLASH [163]: a multi-physics code that, when employed for CCSNe simulations, typically uses Adaptive Mesh Refinement (AMR) for the hydrodynamics, an approximate treatment of general relativistic effects following Case A of Marek et al. [162], and a two-moment M1 scheme [156] for multidimensional neutrino transport. Most versions of the code do not include effects of nuclear burning, although a more modern version of the code, i.e., FLASH-X, is currently under development [164];
• 3DnSNeIDSA [34]: it uses an approximate treatment of general relativistic effects following Case A of Marek et al. [162], and neutrino transport is solved using a spectral method based on the isotropic source diffusion approximation (ISDA) [165], where multidimensional effects are treated using a ray-by-ray approximation. The most recent version of the code is also able to solve the magnetohydrodynamic (MHD) equations [166]. The effects of nuclear burning are usually not taken into account, although some versions of the code include small nuclear networks at low densities [167];
• PROMETHEUS-VERTEX [168,169]: it uses an approximate treatment of general relativistic effects following Case A of Marek et al. [162], and neutrino transport is solved using a variable Eddington factor closure. The multi-dimensional aspects are treated using the so-called ray-by-ray plus approximation. For regions not in nuclear statistical equilibrium a “flash-ing” approximation is used to quantify the heating due to nuclear burning;
• CHIMERA [170]: it uses an approximate treatment of general relativistic effects following Case A of Marek et al. [162], and neutrino transport is solved using a multi-group flux-limited diffusion (MGFLD) approximation following [81]. Multi-dimensional neutrino transport effects are treated using the so-called ray-by-ray approximation. Moreover, it is able to handle regions not in nuclear statistical equilibrium (NSE) using a large nuclear reaction network with up to 160 species [171];
• COCONUT-FMT [172]: it solves the general relativistic metric equations using the extended conformal flatness condition. Neutrinos are treated using a simpler, fast multigroup transport (FMT). Recently, it has been extended to solve the MHD equations [173] and Newtonian hydrodynamics. The general relativistic hydrodynamics has also been coupled with other neutrino transport codes, like VERTEX [29]. The “flash-ing” approximation from Rampp and Janka [168] is used to quantify the heating due to nuclear burning;
• Aenus-Alcar [174]: it uses an approximate treatment of general relativistic effects following Case A of Marek et al. [162]. The neutrino transport is solved using a multi-dimensional two-moment M1 scheme. For regions not in nuclear statistical equilibrium a “flash-ing” approximation is used to quantify the heating due to nuclear burning [168]. Versions of the code that include a full nuclear burning network are currently being developed [175];
• SPHYNX [176]: the only smoothed particle hydrodynamics (SPH) code (to the authors’ knowledge) employed to simulate CCSNe, capable of reaching similar resolutions to the codes mentioned above. Neutrinos are treated using the relatively simple spectral leakage scheme described in Perego et al. [177]. Nuclear burning can be incorporated to the EOS routines modularly, although no simulations including it have been tested to the author’s knowledge;
• GRMHD codes. Recently, modern GRMHD codes have been developed to perform CCSNe simulations. Kuroda et al. [157] developed a code that solves the spacetime equations using a BSSN formalism [178,179], solves the MHD equations, and uses a multi-dimensional two-moment based treatment of multidimensional neutrino transport. Other codes that follow a similar approach are currently under development, and several new multi-physics, GRMHD codes are likely to be employed for production runs of CCSNe in the next few years.

4. CCSNe in Spherical Symmetry

Due to their low computational cost and (relatively) simple implementation, spherically symmetric simulations are still widely used, and comparisons of several different codes [97] show good agreement (much better compared to 2D and 3D [92,180]) in predicting shock evolution and neutrino properties. Therefore, given the prohibitive computational challenge of running hundreds of high-fidelity CCSNe simulations, several studies employ 1D simulations, where the explosion is triggered using a parametric model (although limited 2D studies with up to $\sim 100$ simulations have recently become feasible [181,182]). There is a plethora of parametric simulations developed over the last thirty years [25,44,183,184,185,186,187,188,189,190], as well as semi-analytic models [191,192,193,194,195], aimed at exploring different aspects of supernovae.

In particular, there is a category of recently developed 1D simulations that incorporate neutrino-driven convection through a parametric model based on Reynolds decomposition of the Euler equations and calibrated using multi-dimensional simulations. Such models [25,187,188,190] are oftentimes referred to as 1D+ models, due to their spherically symmetric nature supplemented by a parametric model for convection, an inherently multi-dimensional phenomenon. As we will show later in this section, some of these models tend to be quite different from previous 1D models, where the explosion is triggered using other methods. They also seem to be in excellent agreement with some 3D simulations across different progenitors, although more detailed studies need to be carried out.

Many of these studies were designed to answer an important question: what causes one star to explode and another one to fail? Two related problems, addressed by some of the same studies, are: (i) whether one can identify a condition that can differentiate between successful shock revival and failed explosion; (ii) whether it is possible to predict the explodability (and explosion properties) based on the structure of the pre-SN progenitor.

4.1. The Critical Luminosity Condition

The first attempts at describing an explosion condition date back to the work of Burrows and Goshy [196]. There, the supernova was framed as a bifurcation problem. The only two drastically different solutions are a successful shock revival or a failed explosion. No intermediate solution is possible, and therefore one can imagine that there exists a critical condition that, if met, would cause an explosion.

As summarized in Section 3, the stalling phase of the shock can be described as a balance between neutrino heating and energy lost by the shock due to the photodissociation of the infalling material. Therefore, the two key quantities that describe the shock evolution are the mass accretion rate $\dot{M}$ and the net energy deposited by neutrinos in the gain region in the unit time $\dot{Q}$ (or, equivalently, the neutrino luminosity $L_{\nu }$). The idea delineated for the first time by Burrows and Goshy [196] is that there is a critical curve in the $\dot{M}$-$\dot{Q}$ plane that, if crossed during the post-bounce phase, would lead to a successful explosion.

Revamped using different formalisms [169,191,197,198,199], the concept of a critical condition has been quite successful, and it was confirmed by several multi-dimensional simulations [200,201,202]. However, the existence of a critical condition was usually only qualitatively shown, and when quantitative estimates were given, they were usually only order of magnitude accurate. Recently, the more quantitative Force Explosion Condition (FEC) was derived [203] with the assumption of spherical symmetry, and extensions to multi-dimensions are currently underway [204].

4.2. The Explodability Problem

Historically, the traditional idea about explodability was that less massive stars $M\lesssim 20$ ${\mathrm{M}}_{\odot }$ successfully explode and form neutron stars. More massive stars $M\gtrsim 20$ ${\mathrm{M}}_{\odot }$ lead to failed supernovae and form black holes. This picture has however been challenged by recent studies [148,181,186,188,194,205] that showed the existence of “islands of explodability” throughout the mass range 9–120 ${\mathrm{M}}_{\odot }$. However, there is a fundamental disagreement among these studies regarding where these islands are located. The one common trait highlighted in all those papers is the importance of the Si/Si-O interface in determining the outcome of the explosion.

The first study that provided a quantitative criterion to predict the explodability of massive stars solely based on the pre-SN progenitor, was the one by Ertl et al. [186]. Their criterion was based on two quantities: (i) $M_4$, the mass coordinate at which the specific entropy equals 4 (i.e., a proxy for the Si/O interface) and (ii) ${\mu }_4$ the density gradient at that mass coordinate. That criterion was able to correctly predict the outcome of their 1D simulations. However, as mentioned in Section 3, 1D simulations do not self-consistently explode. Therefore, to achieve an explosion, Ertl et al. [186] artificially increased the number of neutrinos coming from the PNS following the parametric model of Ugliano et al. [184], calibrated on observations of SN1987A [206,207,208]. Therefore, their criterion correctly predicted the outcome of those artificially exploded simulations, and it was not clear if that would translate to high-fidelity, self-consistent explosions. From their simulations, they found that progenitors with masses 22–25 ${\mathrm{M}}_{\odot }$ and 28–80 ${\mathrm{M}}_{\odot }$ lead to failed SN and form black holes.

Later studies [181,188,205,209] also found islands of explodability, but at very different masses. Two of those studies [181,205], independently derived a criterion to predict the explosion based on properties of the Si/Si-O interface, with a similar approach to Ertl et al. [186]. However, only one parameter was used to formulate the criterion. This parameter, in both cases, could be expressed as $\delta {\rho }_{*}^2/{\rho }_{*}^2$, where $\delta {\rho }_{*}$ is the magnitude of the density drop occurring at the Si/Si-O interface, and ${\rho }_{*}$ is the density at which the jump occurs. Boccioli et al. [205] found explosions when $\delta {\rho }_{*}^2/{\rho }_{*}^2>0.08$, whereas Wang et al. [181] found explosions when $\delta {\rho }_{*}^2/{\rho }_{*}^2>0.078$, which is a discrepancy smaller than 3%.

The main difference between those two studies is that Wang et al. [181] used ∼100 2D simulations, whereas Boccioli et al. [205] used ∼300 1D+ simulations where the explosion was achieved using STIR, a parametric model for $\nu$-driven convection based on Reynolds decomposition and time-dependent mixing-length theory. The main parameter of that model is ${\alpha }_{\mathrm{MLT}}$, as defined in Equation (1). The excellent agreement between those two independent studies is a testimony to how 1D+ models can still be employed in some cases without having to compromise on the physics. Moreover, it shows that the Si/Si-O interface represents an important feature of the SN progenitor that has a huge impact on determining the outcome of the explosion. An example of the Si/Si-O interface is shown in Figure 3 for a 20 ${\mathrm{M}}_{\odot }$ progenitor from two different stellar evolution codes. The general idea is that, once the Si/Si-O interface is accreted through the shock, the mass accretion rate suddenly drops due to the drop in pre-shock density. This decreases the ram pressure and the shock experiences a temporary expansion that, if supported by sufficient neutrino heating (oftentimes in the form of neutrino-driven convection), can turn into a runaway explosion.

The explodability prediction of the criterion derived by Boccioli et al. [205] (but in the case of [181] it would be quite similar) is shown in Figure 4 for 200 1D+ simulations from [148]. The criterion correctly predicts the outcome of the simulations in more than $90\%$ of progenitors and recovers the islands of explodability extremely well. One can see that the progenitors forming black holes are mostly the low-mass ones with $12<M<15$ ${\mathrm{M}}_{\odot }$, as well as a few around 18, 28, and ≳100 ${\mathrm{M}}_{\odot }$. The progenitors used are the same as the ones employed by [148,186]. However, the explodability derived is quite different, as shown in Figure 5. The black-hole-forming regions are completely different, which highlights how large the uncertainties in the explosion mechanism are, and it shows that multi-dimensional simulations are needed to determine whether parametric 1D and 1D+ simulations correctly model the explosion. As mentioned in the previous sections, large sets of multi-dimensional simulations have now become feasible, and the explodability found by Boccioli et al. [205] and Wang et al. [181] generally agrees with the results of a dozen high-fidelity 3D simulations performed by Burrows et al. [37]. This shows how 1D+ simulations can now well reproduce the qualitative general trend and evolution of high-fidelity multi-dimensional simulations, and can therefore still prove to be useful tools to explore the parameter space of supernovae.

Impact of Stellar Evolution Uncertainties on the Explodability

It is important to point out that the supernova is only sensitive to the stellar profiles right before collapse. That is, it is not directly related to the Zero-Age Main Sequence (ZAMS) mass and metallicity of the star. Therefore, as seen in Section 2.2.7, differences in stellar evolution calculations will yield different pre-collapse profiles, and therefore different explosion outcomes. An example of this can be seen by comparing the upper and lower panels of Figure 6, adapted from [205]. The upper panel shows the explodability of stars computed using the stellar evolution code FRANEC; the lower panel shows the explodability of stars computed using the stellar evolution code KEPLER [148] (lower panel).

A detailed comparison between codes is extremely challenging given the differences in algorithms, assumptions, and overall physical processes considered. It is however apparent that there are many uncertainties in the stellar evolution of massive stars (see Section 2.2), and that these uncertainties hugely affect the final structure of a star, which in turn changes the explodability. Moreover, multi-dimensional effects during the evolution of massive stars (especially in the last burning stages) are not well understood, as mentioned in Section 3.4. Finally, the explodability discussed here assumes single-star evolution, but we know that about two-thirds of massive stars are in binary systems, and that can also significantly change the final structure of the star, and therefore the outcome of the explosion [61,210].

Figure 6. Figure adapted from [205]. The upper panels show the outcome of simulations run using the stellar evolution code FRANEC [211], whereas the lower panel shows the outcome of simulations run using the stellar evolution code KEPLER [148]. The details of the simulations can be found in [205]. Grey areas indicate that no pre-SN progenitor at that ZAMS mass was available, and therefore no simulation was performed.

Figure 6. Figure adapted from [205]. The upper panels show the outcome of simulations run using the stellar evolution code FRANEC [211], whereas the lower panel shows the outcome of simulations run using the stellar evolution code KEPLER [148]. The details of the simulations can be found in [205]. Grey areas indicate that no pre-SN progenitor at that ZAMS mass was available, and therefore no simulation was performed.

5. Propagation of the Shock and Explosive Nucleosynthesis

In the case of a successful explosion, the re-energized shock escapes the Fe core and is free to propagate through the mantle of the star, accelerating and pushing away the stellar material above the forming PNS. The result is that part of the envelope is ejected and pollutes the surrounding interstellar medium, while the innermost layers above the Fe core fall back onto the PNS because of its strong gravitational field, accreting on the central compact object. Due to the lack of a full self-consistent theory for the final stages of the hydrostatic evolution of a massive star, the core-collapse, and the post-bounce evolution (see previous sections), the late stages of CCSNe are usually simulated by injecting a certain amount of energy at the outer edge of the Fe core, and then following the propagation of the shock with semi-analytical prescriptions or hydrodynamic simulations.

The injected energy is often calibrated to obtain a kinetic energy of the ejecta of the order of ${10}^{51}\mathrm{erg}$ (1 foe) and it may be injected in the form of thermal energy (thermal bomb), kinetic energy (kinetic bomb), or by giving a certain velocity to the outer edge of the Fe core (piston). Usually, in a CCSN calculations, the mass-coordinate that separates the ejecta from the fallback is called mass-cut and it is an arbitrary quantity. It can be fixed at the end of the simulation to eject a certain amount of ⁵⁶Ni (e.g., [8]), by using interpolation formulae (as used, e.g., by [49]), or as the result of the hydro calculation [192,212,213,214]. Regardless of the method used to simulate the supernova explosion, the passage of the shock wave has significant effects on the stellar material above the Fe core.

5.1. Explosive Nucleosynthesis

During its propagation through the stellar structure, the shock locally induces compression and heating. The sudden variation of temperature and density in each layer of the star triggers the so-called “explosive nucleosynthesis”, i.e., a variation of the chemical composition due to nuclear reactions occurring on a timescale of a few seconds. It is possible to derive the basic properties of the explosive nucleosynthesis without the aid of hydrodynamic simulations, under the assumption that the shock is adiabatic, radiation dominated, and it is propagating in spherical symmetry. The relation between the densities in the regions behind and after the shock wave, in case of strong shock limit [215,216], is:

$${\rho }_{\mathrm{shock}}=\frac{\gamma +1}{\gamma -1}{\rho }_{\mathrm{pre}};$$

(2)

with $\gamma$ the adiabatic index. For two typical values of $\gamma$, i.e., 5/3 and 4/3, we obtain ${\rho }_{\mathrm{shock}}=4{\rho }_{\mathrm{pre}}$ and ${\rho }_{\mathrm{shock}}=7{\rho }_{\mathrm{pre}}$, respectively. A similar relation holds for the temperature:

$$T_{\mathrm{shock}}\propto {\mathrm{Ma}}^2{T}_{\mathrm{pre}};$$

(3)

where $\mathrm{Ma}\gg 1$ is the Mach number in the case of strong shock. Hence, the passage of the shock induces a compression that keeps the density almost of the same order of magnitude as in the pre-shock zones, while simultaneously generating a very strong heating of the matter crossed. Furthermore, the relation:

$$\begin{array}{cc}\hfill T_{\mathrm{shock}}& ={\left(\frac{3}{2a}\right)}^{1/4}{\left[(\gamma +1){\rho }_{\mathrm{pre}}\right]}^{1/4}{v}_{\mathrm{shock}}^{1/2}\hfill \\ \hfill & =3.8\times {10}^9·{\left(\frac{{\rho }_{\mathrm{pre}}}{{10}^7{\mathrm{g}/\mathrm{cm}}^3}\right)}^{1/4}·{\left(\frac{v_{\mathrm{shock}}}{2\times {10}^9\mathrm{cm}/\mathrm{s}}\right)}^{1/2}\mathrm{K}\hfill \end{array}$$

(4)

shows that the velocity of the fluid hit by the shock is directly related to the shock temperature. In other words, the kinetic energy is converted into internal energy and heats up the matter. Note that ${(\gamma +1)}^{1/4}\sim 1.25$, since the adiabatic index is always between 4/3 and 5/3. Figure 7 shows how different shock velocities lead to different peak temperatures in a semi-analytical simulation. It is worth mentioning that the shock wave accelerates where the quantity $\rho r^3$, i.e., the density of the star multiplied by the radius to the power of three, decreases [2]. This acceleration occurs when the shock wave encounters the interface between the CO core and the He layer. Hydrodynamic simulations show that [88,217,218], as a consequence of the density contrast between the inner CO-core and the outer He layer, a reverse shock (i.e., a shock propagating inward) is formed at this interface. However, the effect of this reverse shock on the amount of fallback is negligible.

The compression and the heating induced by the passage of the shock wave trigger a rapid nucleosynthesis that occurs within seconds, much faster than the characteristic timescale of the pre-supernova nucleosynthesis. In the post-shock zones, pressure, density, and temperature are almost constant in mass and they decrease as matter expands and cools down. In order to determine the typical explosive burning temperatures we require that the nuclear burning timescale $\tau$ for each major stellar fuel (namely H, He, C, Ne, O, and Si, see Section 2) must be comparable with the explosion timescale. Assuming $\tau \sim 1$ s and considering only the principal channel of destruction for each fuel, it is possible to find the minimum temperature required to burn a considerable fraction of a given element, since the explosive timescale depends just on temperature and density. The variation of $\tau$ as a function of $T_9=T/{10}^9$ K is shown in Figure 8. From the Figure, we note that the threshold temperatures are $4\times {10}^9$ K for ²⁸Si, $3.3\times {10}^9$ K for ¹⁶O, $2.1\times {10}^9$ K for ²⁰Ne, and $1.9\times {10}^9$ K for ¹²C. For temperatures $T_{\mathrm{shock}}<1.9\times {10}^9$ K the timescale for nuclear burning is too long compared to the explosion timescale to trigger any relevant nuclear reaction. These reference temperatures are also reported in the background of Figure 7 as horizontal grey bands.

Using these threshold temperatures it is also possible to define the geometrical volumes in which these burning stages occur. Since the shock is radiation dominated, we can write the radius of the shock as a function of the explosion energy and the peak temperature only:

$$R={\left(\frac{3E_{\exp }}{4\pi a}\right)}^{1/3}{T}_{\mathrm{peak}}^{-4/3}$$

(5)

where R is the radial coordinate, $E_{\mathrm{expl}}$ is the explosion energy, and a is the radiation density constant. Equation (5) identifies the spherical volumes within which each explosive burning stage occurs and these volumes depend just on the explosion energy. The specific amount of mass that is involved in the explosive nucleosynthesis within each volume depends on the mass-radius relation in the PSN model, while the relative proportions among the nucleosynthesis products depend both on the $Y_e$ profile and on the pre-supernova abundances.

5.2. Explosive Burning Stages

The explosive Si burning occurs in all the zones heated up to $T_9=4$. For temperatures $T>5\times {10}^9$ K, a full nuclear statistical equilibrium is achieved, i.e., a unique cluster of elements in equilibrium is formed and all the Si is converted into light particles (p, n, and $\alpha$) and later into Fe peak nuclei. Their relative proportions are a function of temperature, density, expansion (cooling) timescale, and local neutron excess. This stage is called “complete Si burning”. If the timescale of the expansion is shorter than the timescale required to reach NSE, the chemical composition is dominated by $\alpha$ nuclei ($\alpha$-rich freeze-out). Otherwise, the main products are Sc, Ti, Co, Zn, and Ni. In particular, ⁵⁶Ni is the most abundant nuclear specie, and it is the main driver of the supernova light curve, through its decay to ⁵⁶Co first and to ⁵⁶Fe then.

If the peak temperature is $4\times {10}^9\mathrm{K}<T<5\times {10}^9\mathrm{K}$, not all the processes are in statistical equilibrium, and two equilibrium clusters of elements are formed. They are defined by $A>46$ and $A<44$, respectively, with A being the atomic mass. Since the matter crossed by the shock at these temperatures is rich in elements with $A\ll 44$, the value $A=44$ is a bottleneck for all the reactions that produce elements with higher A. As a consequence, Si is only partially exhausted and the remaining Si abundance decreases with temperature. For this reason, this stage is called “incomplete Si burning”. In this case, the most produced isotopes are Ca, Ar, and Si. A few processes may overcome the $A=44$ barrier and produce V, Cr, Mn, Fe, and Ni.

In the temperature range $3.3\times {10}^9\mathrm{K}<T<4\times {10}^9\mathrm{K}$ only one equilibrium cluster is established at $24<A<44$. Therefore, also in this case $A=44$ constitutes a barrier for the production of heavier elements. This temperature range is usually reached in O-rich zones, and the most produced nuclear species are ²⁸Si, ³²S ³⁶Ar, and ⁴⁰Ca.

A further decrease in temperature, down to $1.9\times {10}^9\mathrm{K}<T<3.3\times {10}^9\mathrm{K}$, does not allow the formation of any cluster of equilibrium. Therefore, explosive Ne ($T>2.1\times {10}^9$ K) and C ($T>1.9\times {10}^9$ K) burning occur at the same pace as during the hydrostatic evolution of the star. The main explosive products are ²⁴Mg, ²⁷Al, ^29,30Si, ³¹P, and ³⁷Cl for explosive Ne burning and ²⁰Ne and ²³Na for explosive C burning.

Below $1.9\times {10}^9$ K the temperature is not high enough, relative to the explosion timescale, to trigger further major nuclear burning. However, in the case of highly energetic supernovae, some additional nucleosynthetic processes may occur, such as the activation of ²²Ne($\alpha$,n)²⁵Mg in the He shell, which is crucial for the interpretation of the abundances measured in pre-solar silicon carbide grains [219].

Figure 9 shows the explosive nucleosynthesis production within the stellar structure of a selected number of nuclear species for the same FRANEC model shown in Figure 1. The explosion in this case is simulated in the HYPERION code [220] with an explosion energy of 2 foe. The geometrical volumes corresponding to each explosive burning stage (Equation (5)) are represented as grey shaded areas in the figure and are also reported in

Table 1. Explosive burning stages and corresponding volumes in a 20 ${\mathrm{M}}_{\odot }$ star at solar metallicity exploding with $E_{exp}=2$ foe.

Stage	Temperature	Mass	Radius
	(GK)	(${\mathrm{M}}_{\odot }$)	(km)
complete Si	>5	1.86	4657
incomplete Si	>4	2.03	6270
explosive O	>3.3	2.25	8104
explosive Ne	>2.1	2.88	14,806
explosive C	>1.9	3.07	16,919

. This figure clearly shows that the chemical composition is deeply modified by the explosive nucleosynthesis up to a fraction of the C convective shell. The ratio between the explosive and the pre-supernova yields is instead shown in Figure 10. We remind the reader that the yield of a certain element is the integral of its abundance over the mass of the star, from the mass-cut up to the surface and including the wind contribution. The mass-cut is such that the ejecta contains 0.07 ${\mathrm{M}}_{\odot }$ of ⁵⁶Ni and corresponds to 1.85 ${\mathrm{M}}_{\odot }$. The explosion mainly produces elements between Ca and Zn, whereas, for all of the other elements, the main contribution comes from the hydrostatic evolution of the star.

5.3. The Production of Exotic Nuclei from Neutrino Winds

CCSNe may also host some peculiar nucleosynthetic processes that account for the production of a number of exotic nuclei, due to the extreme physical conditions imposed to the stellar matter by the passage of the shock wave. In particular, neutrino-driven winds drastically alter the electron fraction $Y_e$ of the innermost ejecta, reducing or enhancing the neutronization of the material in proximity of the newly formed PNS [221]. Depending on the value of $Y_e$, the ejecta can be either exposed to a neutron-rich (with $Y_e<0.5$) or proton-rich wind (characterized instead by $Y_e>0.5$). In the first case, some of the lightest r-process nuclei between Sr and Ag can be synthesized by rapid neutron captures (weak r–process, see, e.g., [222,223,224]). In the latter case, instead, the production is concentrated toward the left side of the stability valley of the nuclide chart and favors the synthesis of neutron-deficient isotopes of nuclei between Ni and Ag (the $\nu p$–process, see, e.g., [225,226,227]). In both cases, neutrino-driven wind nucleosynthesis may account for the production of rare nuclei heavier than Fe up to Ag.

5.4. The Production of p-Nuclei: The $\gamma$-Process in CCSNe

CCSNe are one of the possible production sites of 35 rare neutron-deficient isotopes of elements heavier than Fe (namely: ⁷⁴Se, ⁷⁸Kr, ⁸⁴Sr, ^92,94Mo, ^96,98Ru, ¹⁰²Pd, ^106,108Cd, ^112,114,115Sn, ¹¹³In, ¹²⁰Te, ^124,126Xe, ^130,132Ba, ^136,138Ce, ¹³⁸La, ¹⁴⁴Sm, ¹⁵²Gd, ^156,158Dy, ^162,164Er, ¹⁶⁸Yb, ¹⁷⁴Hf, ¹⁸⁰Ta, ¹⁸⁰W, ¹⁸⁴Os, ¹⁹⁰Pt, and ¹⁹⁶Hg), the so-called p-nuclei [228,229,230,231,232,233]. These peculiar nuclei can not be synthesized by neutron capture processes. Instead, they are produced via chains of photodisintegrations on trans-iron seeds in O/Ne–rich layers of exploding massive stars ($\gamma$–process), at temperatures ranging between 2 and 3.5 GK. Some of the lightest ones can also be synthesized by the $\nu p$–process (see Section 5.3). Figure 11 shows the abundances of three typical p-nuclei after the supernova explosion of the MESA model (published in [49]) already presented in Section 2 and in the upper panel of Figure 1. Note that the mass region where the $\gamma$–process occurs coincides with the explosive Ne burning zone.

Despite decades of exploration into the $\gamma$–process nucleosynthesis in CCSNe, model predictions continue to exhibit discrepancies when compared to solar abundances [233]. Specifically, two major challenges persist: (1) the average $\gamma$–process yields indicate an underproduction by approximately a factor of three relative to the amount required for explaining the solar abundances of p-nuclei [234], and (2) the p-isotopes of Mo and Ru systematically experience underproduction by more than an order of magnitude in comparison to other p-nuclei [232,235]. Historically, to estimate the contribution of CCSNe to the solar abundances of p-nuclei without relying on galactic chemical evolution (GCE) calculations, the $\gamma$–process production has been traditionally compared with the production of oxygen. Indeed, massive stars stand as the principal contributors to the galactic production of this element [236]. We therefore define the over-production factor $F_{\mathrm{i}}$ for a specific isotope i as the ratio between its mass fraction $X_{\mathrm{i}}$ in the supernova ejecta (calculated as the total integrated yield divided by the total ejected mass) and the corresponding solar mass fraction ${\mathrm{x}}_{\mathrm{i},\odot }$ (in our case from [237]). The average overproduction factor ${\mathrm{F}}_0$, representing all p-nuclei, is defined as $F_0=\left({\sum }_{\mathrm{i}=1}^{35}{F}_{\mathrm{i}}\right)/35$. This average factor serves as a proxy for analyzing p-nuclei production (see, e.g., [231,238]). Since not all the p-nuclei are solely produced by the $\gamma$–process [239,240,241,242,243,244], Roberti et al. [233] defined an additional factor, F_fl, as the average over-production of the three most produced $\gamma$–only nuclei. To evaluate the $\gamma$–process nucleosynthesis in CCSNe, we compare ${\mathrm{F}}_0$ and F_fl with the over-production factor of ¹⁶O (${\mathrm{F}}_{\mathrm{O}16}$). Figure 12 shows the comparison among ${\mathrm{F}}_{\mathrm{i}}$, ${\mathrm{F}}_0$, F_fl, and ${\mathrm{F}}_{\mathrm{O}16}$, in the case of the same model shown in Figure 11. Most of the ${\mathrm{F}}_{\mathrm{i}}$ factors, together with ${\mathrm{F}}_0$, are under-produced relative to ${\mathrm{F}}_{\mathrm{O}16}$. Instead, F_fl is a factor of 1.3 higher than ${\mathrm{F}}_{\mathrm{O}16}$. However, given that p-nuclei are secondary nuclear species, a prerequisite for deeming their synthesis significant is that their over-production must exceed a factor of 2 relative to a primary element such as oxygen [238,245].

Various alternative nucleosynthetic processes and astrophysical phenomena (see [233], for a recent overview of the production sites of p–nuclei) have been proposed to account for the solar abundances of p–nuclei. However, the challenge of explaining the distribution of p–nuclei abundances in the solar composition persists as an open problem.

6. Summary and Conclusions

The evolution, explosion, and nucleosynthesis of massive stars is an extremely challenging problem, both in terms of the computational power required to describe it, as well as in terms of the experimental, observational, and theoretical uncertainties that affect every aspect of this fascinating branch of astrophysics.

The evolution of massive stars is an inherently three-dimensional problem, given the impact that convection (and to a lesser extent also rotation and magnetic fields) has on the life cycle of stars. However, full 3D simulations of the entire life of the star are practically impossible, and therefore only spherically symmetric simulations are feasible. Experimental nuclear reaction rates, mass loss prescriptions, and convective models are the main sources of uncertainty in these models. Additionally, the size of the nuclear burning network can also have a significant impact, especially in the late phases of the evolution. Finally, different algorithms for nuclear burning, hydrodynamics, and convection can also significantly affect the evolution, and assessing the impact and correctness of these numerical methods on the results is far from trivial. Therefore, different stellar evolution codes often predict different evolutionary paths for a given star, even when adopting similar physical prescriptions. Some promising attempts at simulating small parts of late-stage evolution in 3D are currently underway, and they will be able to provide a better insight into how convection operates in different stellar environments.

An important caveat is that the above discussion is related to single-star evolution. However, about 70% of massive stars occur in binaries. Depending on the properties of the stars and their orbit, binary interactions can dramatically change the evolutionary paths of stars. Since only spherically symmetric simulations are currently feasible, correctly describing these interactions is extremely hard, but promising efforts are underway [246].

Uncertainties in single-star evolution also indirectly affect the explosion phase. A core-collapse supernova directly depends on the thermodynamic profiles of the star right before collapse. Therefore, in reality, it is not sensitive to quantities such as ZAMS mass or metallicity. This means that, since uncertainties in the stellar evolution change the final structure of a star with a given ZAMS mass and metallicity, they will also indirectly change the outcomes of a supernova. After the core collapses, there are additional uncertainties that can affect the supernova dynamics. The largest ones are arguably the equation of state of nuclear matter and neutrino opacities. Constraining the equation of state of nuclear matter is particularly challenging since direct experiments at such extreme thermodynamic conditions are not feasible in a laboratory. Therefore, several constraints are based on indirect measurements of microscopic quantities as well as astrophysical observations and theoretical calculations. Neutrino opacities, instead, rely on theoretical calculations and approximations. In recent times, self-consistent simulations of the explosion phase in 3D have become feasible and were able to shed light on several multi-dimensional effects (neutrino-driven turbulent convection in primis) which play a key role in the explosion. Finally, recent studies based on physically reliable simulations have started to quantitatively analyze the properties of the explosion, and promising efforts towards predicting the outcome of the explosion based on the pre-supernova properties have already been carried out [181,186,194,203,205]. In the near future, large suites of 3D simulations will also be able to shed some light on this paramount explodability problem.

The chemical composition resulting from CCSNe fundamentally relies on the characteristics of the explosion. At present, the computation of explosive nucleosynthesis using extensive nuclear networks (i.e., including thousands of nuclear species and tens of thousands of reaction rates) is almost exclusively based on post-processing the outcomes derived from one-dimensional simulations. Moreover, the separation between the ejecta and the remnant still depends on a mass-cut that is arbitrarily chosen. However, given that CCSNe are intrinsically multi-dimensional phenomena, relying on yields from 1D models introduces uncertainties that remain challenging to quantify, especially within the context of GCE calculations and studies on the enrichment of the interstellar medium. Relying on 1D simulations may not fully capture the intricate dynamics and nuanced processes inherent in the multi-dimensional nature of core-collapse supernovae, prompting concerns about the accuracy and comprehensiveness of the derived chemical composition. This discrepancy underscores the urgent need for advancements in simulation techniques capable of better encapsulating the intrinsic multi-dimensional nature of core-collapse supernovae, thereby enhancing our capacity to predict and interpret resulting chemical abundances. Within this context, several groups exploring the 3D CCSN problem are yielding promising results in simulating the latest stages of massive star evolution and early phases of the explosion, as well as its impact on chemical evolution [13,14,40,88,247].

The intricate life cycle of a star, from its formation up until its explosive demise, is a fascinating computational problem. The first semi-analytical and computer models in the fifties and sixties started tackling this extremely tough computational challenge that astrophysicists are still facing today. The tremendous progress in our understanding of stellar evolution, explosion, and subsequent nucleosynthesis was enabled by increasingly more sophisticated codes and increasingly faster CPUs (and, more recently, GPUs). The Exascale era is upon us, and faster and more flexible supercomputers have been recently deployed. This will facilitate large-scale, multi-dimensional, multi-physics simulations that will be able to investigate many phenomena that are currently very uncertain: mass loss, the role of magnetic fields and rotation in the late stages of stellar evolution, the explosion of a CCSN from collapse to shock breakout, nucleosynthesis networks coupled directly to hydrodynamic simulations, and many others.