Physical mechanism of core-collapse supernovae that neutrinos drive

Abstract

The current understanding of the mechanism of core-collapse supernovae (CCSNe), one of the most energetic events in the universe associated with the death of massive stars and the main formation channel of compact objects such as neutron stars and black holes, is reviewed for broad readers from different disciplines of science who may not be familiar with the object. Therefore, we emphasize the physical aspects than the results of individual model simulations, although large-scale high-fidelity simulations have played the most important roles in the progress we have witnessed in the past few decades. It is now believed that neutrinos are the most important agent in producing the commonest type of CCSNe. The so-called neutrino-heating mechanism will be the focus of this review and its crucial ingredients in micro- and macrophysics and in numerics will be explained one by one. We will also try to elucidate the remaining issues.

Keywords: supernovae, core collapse, neutrinos, neutron stars, hadronic equation of state, radiation transfer

Introduction

It has been almost a century since Baade and Zwicky^1)^ proposed the idea that supernova explosion, an optical transient as bright as its host galaxy, is associated with gravitational collapse of massive stars and is a channel to produce a neutron star, a compact object mostly comprising neutrons. More than half a century has passed since the first serious simulation performed by Colgate and White.^2)^ In retrospect, they were foresighted. They assumed that the explosion is instigated via energy deposition by neutrinos. In fact, investigation of core-collapse supernova (CCSN) mechanism since then has been an accumulation of improvements of this idea in one way or another, although it has not been a straight path, with promising ideas fizzled with a disappointment repeatedly. The CCSN mechanism is probably one of the longest-standing problems in astrophysics.

Observationalists have presented unequivocal evidence^3)^ that the supernovae are indeed associated with the gravitational collapse of massive stars and the subsequent formation of compact objects such as neutron stars and black holes, except for those classified as type Ia,^#1^ which are originated from white dwarfs in binary. The reason why we are so attracted to these core-collapse events among various astronomical phenomena is (1) being the final part of the evolution of massive stars, the CCSNe and the formation of compact objects are one of the major issues remaining in stellar evolution theory, one of the pillars of theoretical astrophysics; we are interested in which stars produce a black hole instead of a neutron star. It is indeed an important input for stellar population synthesis calculations and for some cosmological simulations. (2) The CCSNe are the main driver of chemical evolution in the universe^4)^; they enrich the universe with various nuclides by ejecting stellar envelopes abundant with heavy elements newly synthesized in the progenitor star during its quasi-static evolution and the explosion. Some neutron-rich elements are produced by rapid neutron capture process (r-process); although they are not supposed to be the main contributor these days, the core-collapse supernovae may still not be completely ignored.^5)^ (3) In addition to these astronomical interests, they are also attractive from the physics point of view. As explained in the next section, the CCSNe are probably neutrino-driven^#2^; as such, they provide us with an invaluable opportunity to study weak interaction and the properties of neutrinos themselves under the condition that cannot be realized in terrestrial experiments.^9)^ The flavor conversion of neutrinos, also called neutrino oscillation, is one of the hot topics these days; in particular, collective neutrino oscillations^10)^ are attracting many researchers, as unlike other types of neutrino oscillations, they are nonlinear phenomena and show a rich variety of behavior. The presence of neutrinos in abundance is crucial for this mode of neutrino oscillation to occur, and that is why CCSNe are the focus of this subject. (4) Neutrinos copiously emitted from individual CCSNe have been accumulated in the universe to form a diffuse background ever since the very first supernova had emitted neutrinos^11)^; as such, it encodes the history of supernova explosions and hence that of the massive-star formation in the universe. It is an important target of terrestrial neutrino detectors, such as Super Kamiokande.^12)^ (5) Neutrinos also serve as a probe into the properties of dense and hot nuclear matter.^9)^ The luminosities and spectra of neutrinos diffusing out of a protoneutron star (PNS), which is produced postcore bounce and is still very hot and rich in protons and leptons, are affected by the equation of state (EOS). As density and temperature regime realized in PNS is inaccessible to the current terrestrial experiments, the detailed observation of those neutrinos from a nearby CCSN provides us with a big chance to inspect the nuclear matter.^13)^ (6) The CCSNe are a promising source of gravitational waves (GWs); in fact, they are generated via a couple of processes and, if observed, will provide us with information on rotation and turbulence in massive star cores and on oscillations of PNS.^14)^ These are the representative rationale for supernova research, and there are actually a lot more.

The CCSN is not a rare event in the universe. In fact, we observe multiple explosions a day on average these years. What is rare and we are all looking forward to is a nearby CCSN, in particular, the one in our Milky Way Galaxy, which would allow us to observe it with not only electromagnetic waves but also neutrinos and GWs. Unfortunately, SN1987A remains the only supernova to date, from which neutrinos were detected.^15,16)^ Even SN1987A was not a Galactic supernova but occurred in Large Magellanic Cloud, which is why the terrestrial detectors observed only a dozen of neutrinos. The lack of CCSNe in our neighborhood in the modern time has forced us to rely heavily on theoretical explorations in the study of the CCSN mechanism. In fact, we have tried to perform as realistic simulations as possible, first in one spatial dimension (1D) under the assumption of a spherical symmetry, then in 2 dimensions (2D) with axisymmetry, and finally, in 3 dimensions (3D) without any symmetry assumptions. This is indeed a long endeavor, improving the treatment of various physical processes one after another from the top of the priority list down. They include a numerical scheme for neutrino transport and general relativistic gravity, various weak interactions for neutrinos of all species, the equation of state of dense and hot matter, to mention a few. This process is still underway in fact.

In the past few years, however, we may have reached a point at which the crucial ingredients for neutrino-driven explosion have been identified and their essential workings are understood, as multiple groups succeeded in producing explosions across progenitors of different characteristics in their realistic simulations. A consensus is now developing that neutrinos can in fact generate CCSNe. Of course, we still need to be cautious, as CCSN is a highly nonlinear complex system. Over the years, we have repeatedly seen in the CCSN study that when one improves or implements a process that seems favorable for explosion at the first glance, a feedback, often unexpected, occurs, and the result turns out to be not as expected. Herein, we will examine the supposedly crucial elements one by one, emphasizing their physical aspects, and try to clarify what issues remain to be addressed further.

This review is intended for rather broad readers from different disciplines of science who are unfamiliar with CCSNe but have physics expertise. We emphasize the explanation of the physical processes themselves that occur in various stages of CCSNe rather than the presentation of simulation results. We hope that some topics will touch their own subjects somehow and get them interested in CCSNe as a whole. As such, the references in this review are highly the choice mainly depends not on the importance of their contributions to the issue of CCSNe but on the convenience for the readers in understanding the issue quickly and for the authors in explaining the subject. Therefore, they are neither the first nor the most comprehensive ones in most cases. We refer those readers who are attracted by the subject and want to understand it deeply to more comprehensive review papers for experts^6–8)^ (and also to papers cited in the references for the individual subjects).^#3^

The manuscript is organized as in the next section, we outline the neutrino-heating mechanism and then explain in detail how it works, listing important elements for its quantitative studies; in Section 3, we summarize the current state of our understanding, citing some results of recent simulations that we think are important; in Section 4, we provide some critical assessments of the achievements up to the present time to suggest what to do more in the future; in Section 5, we conclude this review.

Neutrino-heating mechanism

In this section, we will describe the scenario of the neutrino-driven explosion and explain how it works, paying attention to its physical aspects. We will then summarize the essential ingredients for high-fidelity simulations of CCSNe together with the latest observational knowledge.

Physical scenario.

Figure 1 shows how CCSN is supposed to proceed. At some point in the final stage of the massive-star (≳8 M⊙ ) evolution, the gravitational collapse of the core occurs. The contraction of the core is halted when the nuclear saturation density (∼3 × 10^14^ g/cm^3^) is reached and a shock wave is generated. This shock wave is the mediator of the supernova explosion, pushing outward the matter ahead as it runs through the outer part of the core initially and the stellar envelope later. When it breaks out of the stellar surface, optical brightening occurs.

Figure 1. (Color online) Scenario of core-collapse supernova explosion.

As indicated in the figure, the longevity of massive stars is roughly 10^7^ yrs and a core of the white-dwarf size is formed at the center toward the end of their quasi-static evolutions. The core is mainly composed of iron-group elements (or Ne and O for the stars from ∼8 M⊙ to ∼10 M⊙ ). The upper left picture in the figure that shows a meridian section of such a massive star before the core-collapse is very inaccurate as the core is actually very small, ∼1/10^5^ stellar size, and the onion-like structure may not be so orderly. In fact, as we will see in the following sections, convections in O- and Ne-burning shells are so violent that these layers are nonspherical and sometimes merge with each other.

Mechanical equilibrium of a spherical star is unstable to radial perturbations when the adiabatic index Γ ≡ (∂ ln p/∂ ln ρ)s becomes smaller than 4/3, where p, ρ and s are pressure, density, and specific entropy, respectively. The core just prior to collapse is mostly supported by relativistic electrons in Fermi degeneracy, and endothermic photodissociations of heavy nuclei and electron captures on protons inside and outside nuclei effectively reduce the adiabatic index and induce the core collapse at some point. The electron capture proceeds as the density increases and the electron chemical potential rises with core contraction. As a result, matter becomes neutron-rich; this process is called neutronization. Initially, the neutrinos produced by the electron capture escape the star freely. As the density increases (ρ ≳ 10^11^ g/cm^3^), however, the mean free path ℓmfp of neutrino for coherent scattering on heavy nuclei becomes shorter than the core radius and neutrino starts random walks in the core. Then, the neutrinosphere emerges, at which the optical depth τ≡∫Rν∞dr/ℓmfp becomes 2/3. This is a neutrino analog to the photosphere for photons. In this equation, Rν is the radius of the neutrinosphere. Neutrino is eventually trapped inside the core once the diffusion timescale gets longer than the dynamical timescale at ρ ≳ 10^12^ g/cm^3^.^17)^ Then, neutrino is thermally and chemically equilibrated with matter and the neutronization is essentially halted.

The electron fraction Ye ≡ ne/nb is the ratio of the number density of electron ne to that of baryon nb = np + nn, the sum of those of proton np and neutron nn, and is equal to the proton fraction Yp ≡ np/nb owing to charge neutrality. It is ≲0.3 at this point. This value is larger than the typical value, Yp ≲ 0.1, in the neutron star, implying that the collapsing core is not yet so neutron rich as the neutron star. It is divided into inner and outer cores, each evolving in self-similar manners. Mass of the inner core is essentially given by the Chandrasekhar mass, which is the maximum mass that Fermi-degenerate electrons can support, and is ∼0.5 M⊙ (Ye/0.3)^2^. The core contraction is decelerated when the central density exceeds the nuclear saturation density at which matter becomes very stiff (Γ ≳ 2). The inner core, which has subsonically contracted, integrally bounces and generates a shock wave at the boundary between the two cores, where they collide with each other headon. The outer core is supersonically collapsing and is causally decoupled from the inner core.

Pushed by the inner core, the shock wave rams through the outer core, which is still contracting. The more massive the inner core is, the stronger the shock wave is, as the former is a piston to push the latter. The outer core, on the other hand, is less massive, the fact that is also favorable for explosion as the outer core is an obstacle for the shock wave. The value of Ye is determined by the electron capture and neutrino transfer in the core. The value of Ye decreases as the theoretical estimate of the electron capture rates becomes sophisticated in the modern CCSN simulation (see the next section for more information).

As the shock wave propagates through the outer core, it is continuously weakened by photodissociations of the heavy elements that trespass the shock front and are heated. After the shock wave passes the neutrinosphere, the neutrino emissions in the shock-heated matter downstream also sap the strength of the shock wave. If the shock wave managed to exit the core despite these obstacles, it could make all the way up to the stellar surface, ejecting the envelope rather easily, as its gravitational binding energy ∼10^50^ erg is much smaller than the explosion energy ∼10^51^ erg.

What we have described so far is a generic picture not limited to neutrino-heating scenario. The issue now is that the shock wave cannot reach the core surface directly in most cases. The energy loss of the shock wave via the photodissociations and the neutrino emissions is so severe that the shock wave is stalled in the middle of the outer core and becomes a standing wave, which no longer propagates and remains at almost the same location; the infalling matter is no longer pushed back outward and just passes through the shock wave down onto a PNS.^#4^ If this situation continues, the PNS mass would monotonically increase and exceed a critical value at some point, thereby giving rise to a collapse of PNS to a black hole. Therefore, we need an agent to reinvigorate the stagnant shock wave somehow.

Here, neutrinos play a key role. In fact, they are copiously emitted from the neutrinosphere with the total amount of energy of ∼10^53^ ergs over the next tens of seconds. This is indeed much larger than the typical explosion energy ∼10^51^ ergs, implying simply that not all energy is available for explosion. For one thing, the shock revival by neutrinos will not be possible at late times (t ≳ 1 sec). Then, roughly 10^52^ ergs that are emitted earlier can be used. The second thing is that the shock wave is standing outside neutrino sphere, i.e., in optically thin region; this implies that most of the neutrinos are just streaming past the shock wave. In fact, the optical depth of the relevant region downstream the shock wave is ∼0.1, i.e., ∼10% of the energy will be actually absorbed by the matter. The fact that this proportion looks just right (∼10^52^ ergs) × (∼0.1) ≈ 10^51^ ergs indicates that we need a quantitative evaluation of the actual value.

Workings of neutrino-heating mechanism.

One may wonder why the neutrino heating occurs. Following the arguments in Refs. 18 and 19, we will explain its workings more in detail in this section. Figure 2 exhibits a schematic picture of the postbounce core, where the shock wave is stagnant, with matter falling through it onto PNS; neutrinos are mostly emitted from the neutrinosphere and are partially absorbed at large radii; they are also locally emitted everywhere. If this local emission dominates over absorption, a net cooling occurs there, whereas a net heating occurs if the opposite is true.

Figure 2. (Color online) Schematic picture of postbounce core.

We first consider the condition, under which the absorption could match the local emission and there should be no net heating or cooling. The local absorption and emission rates may be written as

and

respectively. In these equations, Lν is the total luminosity of electron-type neutrinos ν~e~ and antineutrinos ν¯e ^#5^; Tν and T are the matter temperatures at the neutrinosphere and at the radius of our concern, respectively. The first factor on the right-hand side of Eq. [1] originates from the local energy density of the neutrinos, whereas the second factor comes from the cross section of the neutrino absorptions, which is approximately proportional to the square of the neutrino energy. The right-hand side of Eq. [2] considers the fact that the local energy density of electrons and positrons is roughly proportional to T^4^. They locally produce neutrinos when they are captured on protons and neutrons, respectively. From these equations, we find that the net heating would be vanishing everywhere outside the neutrinosphere if

were satisfied. We also consider here the fact that the emission is roughly balanced with the absorption at the neutrinosphere.

The actual temperature profile is drawn schematically in Fig. 3. Near the neutrinosphere (r ∼ Rν), where electrons are Fermi degenerate, the temperature is almost constant, implying that the local emission is dominant over the absorption and the net cooling occurs. This is understandable, as the neutrinosphere is the place where neutrinos start to be decoupled from matter. At larger radii (r > Reos), where the degeneracy is resolved and electron-positron pairs become abundant, the temperature is observed to drop with radius roughly as ∝1/r in realistic simulations. This means that the absorption catches up with the cooling at some point (r = Rg) again unless the shock wave (r = Rs) is encountered at a small radius before that happens. Fortunately, the shock wave generated by core bounce is strong enough to hover at a large radius and the gain radius, outside which the neutrino absorption overwhelms the local emission and the net heating of matter occurs, normally emerges. This region of the net matter heating is called the gain region and is responsible for shock revival.

Figure 3. Schematic picture of temperature profiles.

The existence of the gain region is certainly a necessary condition for the shock-heating mechanism to operate but is not a sufficient condition. This may be demonstrated^20)^ by constructing spherically symmetric, neutrino-irradiated, and steady accretion flows that pass through a standing shock wave (see Fig. 2). They are time-independent solutions of hydrodynamical equations augmented with Rankine-Hugoniot relation at the shock wave. The neutrino luminosity and mass accretion rate are the parameters that specify each flow.^#6^ In fact, this can be regarded as an eigenvalue problem, in which the radius of the shock wave is searched for a given pair of the neutrino luminosity and mass accretion rate.

There is indeed a critical neutrino luminosity^20)^ for a given mass accretion rate, above which no steady solution exists. This may be regarded as an indication of shock revival, the picture that has been confirmed in various ways over the years. For example, it was demonstrated that the steady accretion flows become unstable to radial perturbations at the critical luminosity, which suggests that what follows is a runaway expansion of the shock wave, i.e., shock revival.^21)^ The existence of the critical luminosity were demonstrated with idealized but more realistic simulations.^22)^ It has become a common practice also in realistic simulations of CCSNe to plot the neutrino luminosities and mass accretion rates at shock revival (if it occurs) for different models as a critical curve.

The attempt to analytically derive the critical curve has been also made.^18)^ It may be obtained from the comparison of heating timescale with advection timescale (see also Ref. 19). In fact, the realistic simulations of CCSNe have demonstrated over the years that the former timescale becomes shorter than the latter timescale around shock revival. It is intuitively understandable as as matter is accreting, the energy deposition by neutrino should be rapid so that matter gains enough energy before it traverses the gain region, where the net heating occurs.

The two timescales may be evaluated as follows. The advection time τadv may be defined as

where Mgain is the mass in the gain region and M˙ is the mass accretion rate. The heating timescale τheat, on the other hand, may be given as

where εbind = GMPNS/rgain is the gravitational binding energy per baryon and H is the heating rate given in Eq. [1]. The mass contained in the gain region between the shock wave at rshock and the gain radius rgain can be crudely evaluated as

where ρpost is the density just downstream the shock wave, and we have used the following empirical relations for the matter profile in the postshock region (see Figs. 2 and 3):

The postshock density ρpost and the shock radius rshock may be derived in the following way. The pre- and postshock quantities are related with each other via jump condition at the shock wave, which is written in the strong shock limit as

where ρpre and vpre are the pre-shock density and velocity, respectively; β is a constant determined by the stiffness of matter.^#7^ As matter is free-falling in the upstream of the shock wave approximately, ρpre and vpre are given as

In the downstream of the shock wave, the matter profile is given approximately as Eqs. [7]–[9] and the local neutrino cooling matches the heating at the gain radius, i.e.,

where Tgain is the matter temperature at the gain radius. Using Tgain ∝ pgain1/4 (see Eqs. [7] and [9]) as well as Eq. [11], one obtains from Eq. [14] the following relation between rshock and rgain:

Demanding now that the two timescales, which are now written as

be equal to each other, and using Eq. [15], one obtains finally the following

which reproduces the critical curve approximately. This relation is further elaborated in Ref. 19 to incorporate such effects as turbulence and rotation.

The turbulence downstream the standing shock is believed to be a crucial ingredient in the neutrino-heating mechanism. In fact, it has been repeatedly demonstrated by the experimental and realistic simulations that they tend to reduce the critical neutrino luminosity, thus facilitating shock revival.^22,23)^ The main agent to induce turbulence is hydrodynamical instabilities such as convection and the standing accretion shock instability (SASI),^21,24,25)^ to which the neutrino-irradiated and spherically symmetric steady accretion flow through the standing shock wave is susceptible. In fact, the gain region inevitably has a negative entropy gradient, which can be understood as the specific entropy s of a fluid element satisfies the following

where D/Dt is the Lagrangian derivative along the stream line and vr is the radial component of velocity. In the first equation, we use the time-independence and spherical symmetry of the accretion flow. In the absence of dissipative processes such as viscosity and heat conduction, the neutrino heating and cooling are the only process to change the specific entropy of the fluid element. As H − C ≥ 0 in the gain region by definition and vr < 0 for the accretion, we find that ds/dr ≤ 0 there. As is well known in the stellar evolution theory, the negative entropy gradient is indicative of the occurrence of convection. As matter is flowing unlike in ordinary stellar interior, it is a necessary but not sufficient condition. In fact, it was shown^26)^ that the growth time needs to be shorter than the advection time by a certain factor (∼3). The convection, which normally occurs, gives rise to the turbulence once it becomes full fledged.

The SASI is another hydrodynamical instability^24)^: the spherically symmetric and steady accretion flow through a standing shock wave is unstable to nonspherical perturbations even in the absence of the negative entropy gradient^27)^ and produces sloshing and spiraling motions in the postshock flow. It was first found in idealized simulations^24)^ but has also been observed in realistic simulations since then. The linear analysis was also conducted for the idealized accretion flow and demonstrated that SASI exists as linearly unstable modes.^21)^ The most promising mechanism of SASI is probably a cycle of the inward advection of perturbations in entropy and vortices and the outward propagation of acoustic waves, each reflected near PNS and by the shock wave, respectively, and producing perturbation going in the opposite directions anew.^25)^ The study is still ongoing, however, to understand the underlying mechanism(s) involved in SASI, particularly in the presence of rotation.^28)^ Regardless, what is important is the fact that, fully grown to a nonlinear regime, SASI also induces the turbulence in the downstream of the shock wave.

The reason why the turbulence is helpful for shock revival is multifold.^23,29)^ For example, complex nonradial trajectories of fluid elements prolong the dwell time, i.e., the time it takes the fluid element to exit the gain region, thereby heating up the matter more. Lifting the heated matter and lowering the unheated matter in exchange reduce the cooling in the former and enhance the heating in the latter. The most important of all, however, is probably Reynolds stress that the turbulence produces and is denoted by a second-rank tensor R hereafter. It is the macroscopic analog of the turbulent matter motion to the ordinary stress of microscopic origin generated by thermal random motions. It is expressed in terms of velocity deviation δv = v − 〈v〉 from velocity 〈v〉 of the imaginary mean

where ⟨⋯⟩ represents the average of some kind.^#8^

This turbulence-generated stress acts as an additional pressure to the thermodynamic one and expands the gain region, which is certainly helpful for neutrino heating and shock revival. Even the jump condition at the shock wave is effectively modified. The important thing is that the turbulence is anisotropic, Rrr ≈ Rθθ + Rϕϕ, owing to the inflow of matter on average; in this expression, the subscripts indicate the component of the Reynolds stress tensor. This anisotropy has an interesting implication for shock if one defines the specific turbulence energy as εturb ≡ 1/2|δv^2^| and writes the radial component of the stress tensor as

then one obtains γturb ≈ 2, which should be compared with γ ≈ 4/3 ⋯ 5/3 for thermodynamic pressure. This means that the turbulence energy is more efficient in producing pressure than thermal energy of matter.^23)^

Turbulence-assisted neutrino heating is currently believed to be the mechanism that gives rise to most of CCSNe.^#9^ The condition for shock revival is now well understood. All we have to do is to incorporate all relevant physics in numerical simulations and see if this condition is actually satisfied or not. In Section 3, we will look at the current state of the art by citing some of the important results.

Essential ingredients in CCSNe simulations.

The core-collapse supernovae are not rare, being observed on earth multiple times a day. Only those supernovae that allow us to probe deep into their cores via neutrinos and/or GWs are rare (∼1 in 100 yrs). Therefore, theoretical approaches, particularly with large-scale high-fidelity numerical simulations, have led the research on inner workings of CCSNe till date.

The difficulty in the simulation of CCSNe lies in the fact that many ingredients of different scales and nature are involved. We list them in Fig. 4. They can be classified into macrophysical and microphysical elements. Each of them can be further divided into numerical and physical nature, although the dichotomy is sometimes blurred. In this section, we summarize these ingredients incorporated in the recent realistic CCSNe simulations. In Section 4, we will critically look into them.

Figure 4. (Color online) Various elements in CCSNe simulations. The entries in dark color are limited to some simulations.

We begin with the macrophysical elements, which include (magneto)hydrodynamical processes, self-gravity, neutrino transfer, progenitor models, and nuclear burnings. The microphysical elements, including the nuclear equation of state and weak interaction, will follow later.

Hydrodynamics.

The numerical scheme for (magneto)hydrodynamics is required to stably capture a strong shock wave generated at core bounce. It is mostly based on some approximate solutions of Riemann problem, i.e., the initial-value problem for the initial condition with an arbitrary discontinuity. The issue of relevance is probably to provide a spatial resolution good enough to capture the important features of postshock turbulence. We should mention that hydrodynamical method is still progressing in this respect. All interested readers are referred to Refs. 33 and 34 for instance.

As we mentioned in the previous section, the turbulence is thought to be a crucial ingredient in the neutrino-heating mechanism. Its numerical treatment is also critical. It has been known for a long time^35)^ that vortices interact with one another and cascade into small ones in three-dimensional turbulence. At its saturation, inertial range, in which the velocity fluctuations satisfy Kolmogorov law, is developed in the Fourier space. It is terminated at the scale of the mean free path, which is smaller by many orders of magnitude than the macroscopic scale, and the dissipation occurs there. It has been demonstrated by many simulations (see, e.g., Ref. 23 and references therein) that this picture is essentially true of the turbulence that develops in postshock flow.

The two-dimensional turbulence under Plane symmetry is qualitatively different.^36)^ This is because of the existence of the conserved quantity called enstrophy, Ω≡∫(1/2)ω2d2x=(1/2)∑k|ω^(k)|2 , in addition to the kinetic energy, E=∫(1/2)v2d2x=∑k|ω^(k)|2/k2 in this situation. Here, we assume that matter is incompressible; ω = [rot v]⊥ is vortex perpendicular to the plane of symmetry; ω^ is the Fourier transform of ω. In 2D turbulence, the enstrophy cascades into the smaller scale, whereas the kinetic energy is transferred to the larger scale; this phenomenon is called inverse cascade. This has been also observed in 2D CCSNe simulations (see Ref. 23 and references therein), although they assumed axisymmetry instead of plane symmetry. This inverse cascade is thought to facilitate shock revival in 2D. Therefore, we need to be cautious with 2D simulations.

In the 3D simulations, on the other hand, the resolution is a concern.^37)^ As it is impossible to resolve the entire inertial range, we need to ensure to capture those features in the turbulence that are of critical relevance for shock revival with limited resolution. In some simulations, the treatment of numerical grid is elaborated, as grid geometry can be a source of spherical coordinates, for instance, are singular on Z-axis; in the cartesian coordinates, on the other hand, flows may tend to be aligned with coordinate axes. The use of Ying–Yang grids^38)^ is effective for avoiding coordinate singularity. Deployment of the dendritic coordinates^39)^ is an interesting variant of static mesh-refinement technique to avoid a fine resolution near the center and Z-axis in the spherical coordinates.

Gravity.

General relativistic treatment of self-gravity is very important.^40,41)^ This may be understandable if one compares Schwarzschild radius RPNSSch and the actual radius of PNS RPNS. In self-gravitating system, the following relation

where ϕPNS is the Newtonian gravitational potential per mass. The general relativistic gravity, which is effectively stronger than the Newtonian gravity, makes PNS more compact and hotter than the Newtonian counterpart. Then, the neutrinos emitted from PNS have higher energies and heat the matter downstream the shock wave more efficiently, as the neutrino-absorption rate is roughly proportional to the square of the neutrino energy as mentioned earlier. The gravitational energy liberated by the collapse and tapped by the neutrinos is also large. On the other hand, the area of the neutrinosphere becomes smaller and gravitational well, from which matter should be pulled up, gets deeper; these facts are unfavorable for explosion. The former effects dominate over the latter and the general relativistic gravity is favorable for shock revival.

In most of recent simulations, the general relativistic corrections to the Newtonian gravity are approximately considered. This may be justified if one recalls that they are of the order of ≲30% (see Eq. [22]). The most popular approach is to modify Newtonian potential as follows^42)^: the monopole component in multipole expansion of the Newtonian gravitational potential is replaced by a function constructed to mimic the general relativistic gravity in Tolman–Oppenheimer–Volkov equation, which describes the static and spherically symmetric star in general relativity. The procedure is easy and convenient. One should consider, however, that the effective potential approach is ad hoc and may produce qualitatively different results in some cases.^43)^ The TOV equation is meant for the spherically symmetric stars and is formulated on a particular choice of coordinates (or gauge), and its applicability to rotational cases is unclear. It was demonstrated, by contrast, that spacetime realized in CCSNe is conformally flat approximately even when the core rapidly rotates.^44)^ Although not popular these days, the conformally flat approximation of spacetime may be a good option. Last but not least, we should mention fully general relativistic simulations by Kuroda et al. (see, e.g., Ref. 45 and also references therein).

Neutrino transfer.

The neutrino transfer is another key element in the simulation of CCSNe. In fact, the essential thing in the neutrino-heating mechanism is that the neutrinos emitted from the deep interior of the core are absorbed in the outer region. The precise calculation of the neutrino transfer in between is crucially important to ascertain the success of the neutrino-heating mechanism.

The difficulty with the neutrino transfer originates from the facts that the mean-free path λmfp of the neutrino changes with radius by orders and depends on the neutrino energy εν rather 1/λmfp ∝ εν2 (see also Eq. [1]). The neutrinos have very short mean-free paths (λmfp ≪ r) near the center and their transfer is regarded as the diffusion via random walk, whereas the mean-free path becomes very long (λmfp ≫ r) near the core surface and the neutrinos are streaming essentially freely there. This corresponds to the fact that the neutrinosphere is formed somewhere in between. The location of the neutrinosphere depends on the neutrino energy.^#10^

We need to explicitly consider the energy spectrum of neutrinos in their transfer; we should also accurately treat the transition from diffusion to free-streaming. In fact, it is in this transition region that the neutrinosphere is located and the luminosity and spectrum relevant for the neutrino heating in the gain region at large radii are formed. Under such circumstances, we need to use the kinetic equation to describe the neutrino transfer.

Boltzmann equation is usually the one used for that purpose.^#11^ It is an equation for the distribution function f(r, p) in phase space, which is reduced to Fermi–Dirac

in thermal equilibrium; in this expression, p, μν are the momentum and chemical potential of neutrino, respectively. As three flavors of neutrinos are confirmed to date, we have to deal in principle with six Boltzmann equations including those for antineutrinos.^#12^ They are coupled with one another and should be simultaneously solved, as there are pair processes and flavor-changing reactions as well (see Section 2.3.7). If the possible presence of muons is ignored, no difference is observed between the muon- and tau-type neutrinos.^#13^ The difference between the neutrino and antineutrino of these flavors is also minor and is often ignored in the simulation of CCSNe. Then, we distinguish three types of ν~e, ν¯e , and νx, the last of which collectively stands for νμ~, ν¯μ , ντ, and ν¯τ .

The neutrinos are in neither thermal nor chemical equilibrium, and their distribution functions are time dependent in general. The Boltzmann equation, which describes their time evolutions, can be written as

where v and p˙ are the velocity and time derivative of momentum of neutrino, respectively. The left-hand side of this equation expresses advection in phase space and is essentially the directional derivative of f along the trajectory of the neutrino in the absence of interactions, which should be described by collision terms on the right-hand side of the equation.

One should recognize the different natures of the left- and right-hand side of the Boltzmann equation. The interactions of neutrinos with matter are treated most conveniently in the instantaneous rest frame of matter, whereas the advection in phase space will be simplest in laboratory frame chosen appropriately. In fact, if one uses the rest frame instead for the left-hand side of the Boltzmann equation, the third term will include the aberration effect of O(v/c) in addition to general relativistic effects such as the ray-bending and redshift of O(v2/c2) . The distribution function is a Lorentz invariant; those in different frames are nothing but different parametrizations of the same function. This may lead us to use these two distribution functions in a hybrid way^46)^: the one in the laboratory frame for the left-hand side and the other in the matter rest frame for the right-hand side, appropriately transforming one from the other.

However, it is rather common to use the matter rest frame alone, considering the relativistic effects up to the order of O(v/c) . The magnitude of these relativistic corrections estimated in Eq. [22] is an overestimate at early times, when hot and lepton-rich PNS is not so compact as cold counterpart. This may justify common approximation.

It is also a common practice these days to consider not the Boltzmann equation itself but the equations for the angular moments thereof. The k-th angular moment of neutrino distribution function f is defined as

where k is a nonnegative integer, n is the spatial unit vector that specifies the direction of neutrino momentum and is multiplied k times with f in the integrand; the moment is energy dependent as shown explicitly in the equation (the spatial dependence is omitted for notational simplicity); it is a k-th order spatial tensor (see Refs. 47 and 48 for the general relativistic treatment).

The Boltzmann equation is numerically demanding, as it depends on six arguments, r and p, in addition to time t. By integrating over two angles in momentum space, the moments Mk depend on r and εν. This apparent reduction is of course compensated for by the infinite number of moments. The angular integration of the Boltzmann equation [24] multiplied with k ns gives the time-evolution equation for Mk. Recalling that n is related with v as v = |v|n, we recognize that the resultant equation has the following

where the contributions from the third term with the momentum derivative in Eq. [24] are replaced with ellipses on the left-hand side.^#14^ As the equation for the k-th moment includes at least (k + 1)-th moment, the equation system is not closed unless it is truncated at a finite order.^#15^

In the truncated moment method, one uses some relation to close the system. It is a relation that expresses high-order moments, the equations of which are not solved, in terms of the lower-order ones to be solved. It is most popular at present to truncate the hierarchy at the first order, i.e., only the equations for the zero-th and first moments are retained and the second moment (and higher ones also) is given by the closure relation. As no such functional relation actually exists among moments, it is a phenomenological approximation and, as such, various relations have been proposed in the literature (see Ref. 49 for a recent comprehensive review and references therein). They are physically motivated and constructed so that some constraints for the moments, including the optically thick and thin limits as a prerequisite, should be satisfied. Their performances have been investigated over the years, particularly in spherically symmetric situations. It seems, however, that no such thing exists as the best closure relation, that is, one is better than the others in some cases but not always.^50)^

We have been developing the neutrino transfer code that solves the Boltzmann equation as it is (of course, after the mandatory finite differencing. See Section 3.3). It is called the discrete ordinate or SN, method. This code also adopts the hybrid approach mentioned above, in which the distribution function in the matter rest frame is used for the collision terms on the right-hand side of the Boltzmann equations, whereas the laboratory frame is used for the advection terms on the left-hand side; they are Lorentz-transformed from each other when necessary. The nice thing of this method is, of course, the formulation has no systematic issues involved in the truncated moment method.

On the other hand, it has its own demerits. Among other things, the resolution issue is probably the most serious. The finite-differenced Boltzmann equations are a very large number of coupled nonlinear algebraic equations. The available computational resource puts a strong limit on the numbers of grid points we can afford to deploy. In particular, momentum space tends to be underresolved. In our typical (spatially) 2D simulations under axisymmetry, we use spherical momentum-space coordinates (εν, θν, ϕν), the radial coordinate of which is nothing but the energy. We deploy 20 energy-grid points, extending up to ∼300 MeV. The entire solid angle is covered with 10(θν) × 6(ϕν) grid points. The latter resolution is especially concerning, as the angular distribution in momentum space becomes highly radially outward peaked in the optically thin region. The results of our simulations done so far will be shown in Section 3.3 together with the issue mentioned just now.

The neutrino transfer can also be treated in a stochastic manner using Monte Carlo method.^51)^ As this method is grid free, no resolution issue is associated with the finite differencing as in the SN method in principle. The interactions of neutrino with matter occur in a probabilistic way in reality. The solution of Monte Carlo simulation can serve as a surrogate of the exact solution and is often used for calibrating other methods.^52)^ It has its own difficulties. In fact, the random noise inherent to the stochastic method gets smaller with the number N of particles deployed rather slowly, ∝1/ N . It may not be fast enough for time-dependent problems. In fact, the simulations for the calibration mentioned above are limited to the neutrino transfer on top of steady matter distributions.

The Monte Carlo method in its original form is not suitable for use in optically thick region, as a very large number of interactions should be treated. The application of the Monte Carlo method to the diffusion equation,^51)^ which is an appropriate approximation in the optically thick regime as mentioned above, instead of the Boltzmann equation, is an interesting the two equations are switched particle wise, depending on the optical depth at the particle position. It is also an attractive idea to use the Monte Carlo simulation to obtain a closure relation in situ.^53)^ We are awaiting serious applications of these methods to realistic simulations yet.

Progenitors.

Progenitors of CCSNe are massive stars (≳8 M⊙ ). They are characterized mainly by their mass at birth called zero age main sequence (ZAMS) mass and metallicity. Their profile has a substantial influence on shock revival and subsequent evolution. Particularly important are the compactness of the core, i.e., the mass-radius ratio,^54)^ the density gradient in the envelope,^#16^ and the density jumps at the boundaries of different layers such as Si and O/Ne layers that will hit the stagnant shock wave at the time of relevance for shock revival. In fact, we often observe in the CCSNe simulations that shock revival is triggered by the passage of one of these layer boundaries. The sudden drop in the mass accretion rate associated with it is thought to be one of the important elements for shock revival (see the recent reviews^6–8)^ for references). Shock propagation is accelerated when the density profile satisfies ρ(r)r^3^ < 1. Some authors claim that these quantities can be used to judge the success or failure of a given progenitor without simulations (e.g., Ref. 55).

The structure of the progenitor just prior to the core collapse is determined by the preceding evolution from ZAMS. It is notoriously nonmonotonic with the ZAMS mass.^56)^ In the past, the progenitor model was just an initial condition, although crucially important, for the CCSNe simulation provided by experts of the stellar evolution. This is still true; however, the landscape has been drastically changed in the past 5 years or so, which is worth mentioning here.

It is almost a decade ago that the importance for shock revival of the asphericity in the envelope of progenitors induced by convective activities was first highlighted.^57)^ It was 2 years later that the final few-minute evolution of the progenitor up to the core collapse was multidimensionally computed for the first time^58)^: the simulations demonstrated that the Si burning is indeed violently convective and that the turbulence it induces has an important ramification for shock revival. As the turbulence is expected to be amplified by inward advection,^59)^ it serves as (not so small) seed perturbations for the convection and SASI in the postshock region (see also Ref. 60 for the linear analysis of this process). Other authors^61,62)^ have followed suit since then, confirming their findings (but see also Ref. 63), and showed that the convection in the O/Ne layer may be even more important and, in fact, its structure itself may be qualitatively changed.^64)^

This had been actually highlighted several years earlier.^65)^ In the conventional stellar evolution calculations, which are essentially one dimensional, convection is phenomenologically treated with mixing-length theory, the validity of which is now seriously challenged and a possible extension has been also discussed.^66)^ The convection in other layers at different evolutionary stages has also been investigated by various authors (e.g., Refs. 67 and 68). The stellar rotation was also discussed in this context,^69)^ and an interesting pattern in the resultant stellar rotation, which is reminiscent of the solar rotation, was reported.

The simulation originally meant for the study of CCSNe mechanism in Ref. 57 has a far-reaching influence on the stellar evolution theory much beyond the context of CCSNe. What is important in these progresses from the point of view of the supernova theory, is the fact that the multidimensional simulations of (at least) the late stage of the stellar evolution up to the core collapse is now a part of the CCSNe simulation.

Nuclear reactions.

Nuclear reactions certainly occur during the entire evolution of CCSNe. In fact, various nuclei are either fused or decomposed in different layers of the envelope according to their compositions and thermodynamical conditions. As they infall, following the gravitational collapse of the core, and their temperatures rise, heavier and more diverse nuclei are produced. The abundance eventually approaches the one in nuclear statistical equilibrium (NSE), in which essentially all nuclei are in chemical equilibrium and the following relations hold among their chemical

where μ~A(Z,N~~)~ is the chemical potential of a nucleus composed of Z protons and N neutrons, whereas μ~p~ and μ~n~ are the chemical potentials of unbound protons and neutrons, respectively; they all include the rest mass. The NSE abundance is a function of the temperature, baryon number density, and the proton fraction and is derived from Eq. [27] as a nuclear part of EOS (see Section 2.3.6 below for more details).

In the iron core, NSE is already established by the time of the core collapse. In the envelopes, by contrast, heavy nuclei are not yet in chemical equilibrium and their abundance changes in time toward NSE via nuclear reactions. It is accelerated by the infall of the envelopes, as mentioned above, and NSE is reached at some point either during the infall down to the stalled shock wave or just after the passage. Such evolutions should be calculated according to the following

where ρ is the mass density and NA is Avogadro’s number, Yi = ni/(ρNA) with the number density ni of nucleus i stands for its number fraction, D/Dt is the Lagrangian derivative, i.e., the derivative along a stream line; the three sums on the right-hand side concern reactions with 1–3 reactant nuclei, respectively, that produce or destroy the nucleus of species i; λs encode the reaction rates for these interactions; N s are positive for the creation and negative for the destruction of nucleus i and also play the role of the symmetry factors required to avoid double counting when identical nuclei interact with each other.

This nuclear network calculation was often neglected in the past partly because the nuclear reactions are not so much energetically for example, burning of 0.1 M⊙ of ^16^O to ^56^Fe generates ∼10^50^ erg, which is smaller by an order than the canonical explosion energy of 10^51^ erg. The nuclear burnings were later approximated in the CCSNe simulations by instantaneous conversions of nuclear species accompanied by energy generations, known as flashing approximation.^70)^ Concurrently with the progress in multidimensional treatment of progenitors (see Section 2.3.4), the nuclear network is implemented in some latest advanced codes.^38,71)^

From the Lagrangian derivative in Eq. [28], the evolution of the nuclear abundance is associated with fluid element, which may be approximated by the cell in the numerical grid. It is understandable that the size of the network, i.e., the number of nuclear species considered, is normally rather small, ≲15. They are mostly confined to α-nuclei, that is, integral multiples of the alpha particle.

It was already pointed out more than 15 years ago^72)^ that the energy generation of ∼10^50^ erg by the O-burning that occurs when the O-rich material reaches the stalled shock wave can actually give a substantial impact on shock revival.^#17^ It was based on the 2D simulations that incorporated a nuclear network with 14 alpha nuclei from ^4^He to ^60^Zn. This claim was confirmed later^73)^ in the 1D and 2D simulations, in which shock revival was parametrically induced.

The nuclear reactions are also important for postshock-revival evolution. Matter is expanding this time and the density and temperature are lowering in time. The matter in the downstream of the shock wave has a high entropy per baryon s ≳ 5 kB owing to the heating in passing through the shock wave and, as a result, consists mostly of α particles and nucleons, as they are preferred to heavy nuclei from the point of view of free energy at such high entropies. As the temperature is lowered with the expansion, however, they are reassembled into heavy nuclei again; this process is called recombination. Initially, it is again understood as the change of the nuclear abundance in NSE but at some point later NSE is no longer sustained.^74)^ In particular, the formation of ^12^C via the fusion of three α particles cannot catch up with the decline of density and temperature by matter expansion just as in primordial nucleosynthesis in the early universe, and the α-rich freeze out of nuclear reactions occurs. It has been argued that the recombination will help shock revival and contribute also to the explosion energy^7,75)^ as the recombinations of nucleons into α particles and heavy elements are exothermic. They need to be quantitatively evaluated, however. The network calculation is then required once NSE is terminated.

In fact, the recombination is not the only process of relevance. Very recently, the impact of the nuclear reactions in non-NSE regime on the explosion dynamics of CCSNe was quantitatively explored with a large network and the differences that the simplified treatments such as the flashing approximation and/or the network size can make were investigated in the 1D and 2D simulations.^76)^ It was demonstrated that they affect the elemental composition, particularly the abundance of unbound nucleons, ahead of the shock wave, thus affecting the neutrino heating there and eventually the shock propagation. The energy released by the reactions heats up matter and slows down the accretion, thereby making the shock propagation easier again. The nuclear energy generated in the postshock region, on the other hand, increased the explosion energy by ≲20%.

As the revived shock wave propagates outward further through the envelope, the shock heating triggers nuclear reactions to produce heavy elements in the matter that has just passed through the shock wave. This is called the explosive nucleosynthesis. Although these reactions are also exothermic, their contribution to the explosion energy is minor.

Nuclear equation of state.

In this section, we proceed with the microphysical elements. Among them are the nuclear EOS and the weak interactions of neutrinos. We separately consider them in this order in the following.

The matter in the massive star core is a mixture of hadrons, leptons, and photons. The EOS used in the CCSNe simulation covers a wide range of density, temperature, and electron fraction (which is normally equal to the proton fraction owing to the charge neutrality^#18^): 10^5^ ≲ ρ ≲ 10^15^ g/cm^3^, 0.1 ≲ T ≲ 100 MeV, and 0 ≲ Ye ≲ 0.6. The hadrons are normally nucleons and nuclei. Their abundances differ from layer to layer in the progenitor and are determined by the precollapse evolution of the progenitor except for the core, in which it is given by the nuclear statistical equilibrium, NSE, distribution; although the iron-group elements are most abundant initially (that is why it is called the iron core), various nuclei exist in fact. As the collapse proceeds and the density rises, the NSE distribution changes and heavier and more neutron-rich nuclei become abundant. On the other hand, the electron and positron as well as all flavors of neutrinos and antineutrinos are the main components of leptons.

The entropy per baryon is low (s ∼ 1 kB) in the core just before the collapse because of the efficient cooling by the neutrino emissions in the progenitor in its advanced evolutionary stages after He burning. The entropy remains low during the core collapse, as the entropy productions associated with the electron captures on nucleons and nuclei, which are irreversible processes, are almost canceled by the neutrino emissions. The NSE may be regarded as the mixed phase between the gas and liquid phases of nucleons. The liquid phase is nothing but the uniform nuclear matter at supranuclear densities (ρ ≳ ρ~s~ ≈ 3 × 10^14^ g/cm^3^), in which nucleons are strongly correlated with one another. At very high densities (ρ ≳ a few times ρ~s), more exotic phases, such as those with hyperons, condensed bosons, or quarks may emerge. After core bounce, the shock wave sweeps the outer core, generating entropy (s ≳ 5 kB~) by the dissipation of kinetic energy. Heavy nuclei are mostly dissolved into α particles and nucleons, as that is the NSE distribution at such high entropies and subnuclear densities.

The difficulty in deriving the nuclear EOS is multifold. First of all, nuclear forces are the residual of the strong interactions among quarks that form bound states as nucleons. How to describe them has uncertainties (but see, e.g., Ref. 77). In addition, they are too strong for perturbative treatments to be applicable. This is true not only for quantum mechanical calculations but also for statistical ones. At subnuclear densities, as mentioned above, we need to handle the ensemble of a large number of nuclei (≳1000). This is not a big issue at low densities (≲10^12^ g/cm^3^), where we know the properties of individual nuclei rather well experimentally. The virial expansion may allow us to evaluate EOS in a model-independent manner.^78)^ However, difficulties arise at high densities. They originate from the fact that the NSE distribution depends on the binding energy of nuclei, which are modified by the presence of nucleons and other nuclei through nonvanishing pressure exerted by them as well as via the Coulomb attractions by electrons and are affected by the NSE distribution itself; this means that the properties of nuclei at such densities and their abundance should be simultaneously calculated.^79)^

Thanks to the detections of GWs from mergers of a neutron star with another neutron star^80)^ or a black hole in a compact binary,^81)^ the importance of the nuclear EOS has been more widely recognized these days and the list of EOSs that are available for the CCSNe and neutron star merger simulations has been extended.^82)^ The most popular at present is probably the ones based on the relativistic mean field theory.^83)^ The EOSs based on more realistic nuclear potentials are also available.^13,84,85)^ The nice thing is that they are confirmed to be consistent with the latest observational and experimental constraints.^82,86)^ The comparison of simulation results with different EOSs has also been done.^87,88)^

Weak interactions.

Neutrinos interact with matter via various channels of the weak interaction. It is remarkable that supernova modelers have incorporated more reactions more accurately over the years. This should be understandable, as the neutrino interactions are one of the most crucial ingredients in the neutrino-heating mechanism. In Table 1, we list the reactions that the most advanced simulations incorporate and are becoming the new standard these days.^#19^ One finds that there are emissions, absorptions, and scatterings of all flavors of neutrinos and antineutrinos on leptons, i.e., electron and positron as well as neutrinos and anti-neutrinos, and on hadrons, i.e., nucleons and nuclei. The pair processes include annihilation/creation of e^∓^ pairs and the nucleon bremsstrahlung as well as the pairwise flavor-changing interactions.

These reactions are implemented in the neutrino transfer equations as the collision terms on their right-hand side. For the emission and absorption reactions, they are written in the following

In this expression, R^e/a^ is called kernel for emission or absorption reaction and is a function of the neutrino energy, εν, as well as the thermodynamical variables, ρ, T and Ye, which are omitted here. The kernels are derived from the integral of the relevant matrix elements of the weak currents over the thermal distributions of matter particles involved; the matrix elements, on the other hand, are calculated to the lowest order in Fermi coupling constant GF. In fact, the energy scale (≲100 MeV) encountered in the CCSNe simulation is low enough to justify the use of the Fermi theory of weak interaction. The factor (1 − f) in the first term on the right-hand side of Eq. [29] expresses Pauli blocking for the neutrino in the final state.

For the scatterings, the collision terms are written in general as

Here, R^sc^ is the kernel for this process and is a function of the momenta of incoming and outgoing neutrinos specified by the first and second arguments, respectively. It also depends on the thermodynamical variables as R^e/a^ does. One recognizes again the Pauli-blocking factors in this equation. The scatterings on nuclei can be regarded as isoenergetic, as the rest mass energy of the nuclei, ≫1 GeV, are much larger than the neutrino energies, εν ≲ 100 MeV. Then, R^sc^(p, p′) ∝ δ(p − p′) holds and different energies of neutrinos are decoupled from one another in Eq. [30]. By contrast, the energy transfer is large in the scatterings on electrons and positrons as well as on other neutrinos. They are just like Compton scattering for high-energy photons. Then, essentially all energies of neutrinos are coupled with one another in Eq. [30]. Its treatment is computationally more demanding.

The scatterings on nucleons are subtle. The rest mass energies of nucleons (≈1 GeV) are certainly higher than the typical neutrino energy substantially, and therefore, the energy transfer is small in a single scattering. However, it cannot be neglected because the cross sections of the scatterings on nucleons are larger by roughly an order than those of the scatterings on leptons. Consequently, the accumulation of the small energy transfers in the former is thought to be more important than the single large energy transfer in the latter.^87,89)^

Such a small energy transfer poses another numerical challenge. As mentioned earlier, the energy resolution in the CCSNe simulation is not very high. In fact, we normally deploy ≲20 grid points to represent the neutrino spectrum up to ∼300 MeV. Then, the widths of individual energy bins are larger than the energy transfer in the scatterings on nucleons. Therefore, it is a common practice to deploy a subgrid inside each cell of the original grid to accurately count the number of neutrinos that pass through the cell boundary.^90)^ Fokker–Planck approximation, in which the scattering with a small energy transfer is replaced by the advection in the energy space, may be an interesting idea that is proposed rather recently^91)^ and awaits a serious application.

As the nucleons are composite particles comprising three quarks, their weak currents are more involved than those of leptons. At low energies, this may be encoded in the form factors. For example, charged current of nucleons relevant for the absorption or emission of ν + N1

⇄

ℓ + N2, where ℓ and Ns stand for a charged lepton and nucleons, respectively, can be expressed as

whereas the neutral current of nucleon used for the neutrino ν + N

⇄ ν + N may be given as

In the two equations above, us are the spinors for nucleons, γ^α^ and γ5 are the gamma matrices, and σ^αβ^ = i/2(γ^α^γ^β^ − γ^β^γ^α^); MN is the nucleon mass defined as the arithmetic mean of the proton and neutron masses; GV, GA, F2, Gp, G1N , G2N , and GAN are the form factors parametrized by the transferred momentum q^α^ = p^α^ − p′^α^. In the limit of q → 0, GV → gV = 1, GA → gA = 1.27, and G1p→1/2(1−4sin2⁡θW) , G1n→−1/2 , GAp→gA/2 , and GAn→−gA/2 ; the other form factors are either vanishing themselves or multiplying the terms that are vanishing in this limit. We refer readers, for instance, to Ref. 92 for the concrete parametrization.

In addition to the well-known V − A contributions, one also recognizes tensor and pseudoscalar terms in these expressions, which are also originated from the composite nature of nucleons. In fact, the tensor term is analogous to the term in the electromagnetic current that expresses anomalous magnetic moment and is called the weak magnetism. Since its importance was first pointed out,^89)^ it is commonly incorporated in the CCSNe simulation. The axial vector contribution is normally assumed to be isovector. There may be an isoscalar contribution −gs(q^2^)1/2 as well,^93)^ which may be interpreted as the contributions of strange quarks in the nucleon.^#20^ The current experimental constraint is gs(0) ≈ −0.1.^94)^

It is important to recognize that neutrinos are interacting with nucleons that are themselves interacting with one another via strong nuclear forces. The coherent scattering of neutrino on heavy nuclei may be regarded as an example. In fact, a neutrino is simultaneously interacting with all nucleons localized in a nucleus via nuclear forces, and, as a result, the scattering rate is enhanced by the positive interference of scattered neutrino waves. This rather extreme case demonstrates that the correlations of nucleons can modify their interactions with neutrinos. This many-body corrections to the neutrino reaction rates have been studied and somehow implemented in the recent CCSNe simulations. For example, the dispersion relation of nucleons is substantially modified by the nuclear interactions at high densities (ρ ≳ ρ~s): their effective mass and momentum, being different from those in vacuum, are often evaluated with mean field approximation^95,96)^; on the other hand, random phase approximation (RPA) is sometimes employed to evaluate the vertex corrections.^97–99)^ The axial vector coupling constant gA~ is known to be quenched in nuclei,^100)^ which is also considered at high densities.

At low densities (ρ ≲ 10^12^ g/cm^3^), virial expansion is used^101)^ to evaluate the many-body corrections via nuclear forces. Electromagnetic forces also play roles to induce correlations among heavy nuclei. The correction from this ion–ion correlation was evaluated in Ref. 102 with the Monte Carlo simulations. Its impact on the dynamics in the collapsing phase was investigated soon after that.^103)^ The corrections are commonly incorporated in recent CCSNe simulations. The neutrino emissions and absorptions on heavy nuclei are important elements particularly in the collapsing phase, where electron-type neutrinos are copiously produced by electron captures on nuclei and nucleons, which neutronize matter at the same time. As mentioned earlier, as very heavy and neutron-rich nuclei emerge as the density increases, the theoretical estimation of the weak interaction rates for them are indispensable (see Ref. 104 and references therein for the recent review on the subject). We are still observing updates on these rates (see, for example, Ref. 105 for these recent progresses). The recent CCSNe simulations incorporate some of those rates in one way or another.

Last but not least, muons may not be forgotten. They are essentially absent in the massive star core because of their large mass (∼100 MeV). They begin to thermally populate via electromagnetic processes such as e^−^ + e^+^ → μ^−^ + μ^+^ and γ + γ → μ^−^ + μ^+^ as well as via the weak interactions involving thermally-produced νμ and ν¯μ (see Table 2) after ∼100 ms postbounce when the temperature becomes a few tens MeV and the chemical potential of electron gets higher than the rest mass of muon. In fact, although muons are much less abundant Yμ ≲ 0.05 than electrons even in the late postbounce phase, they accelerate the contraction of PNS, which results in higher luminosities and mean energies of neutrinos and leads to the enhancement of neutrino heating.^106)^ These subtle effects may be important when the success of shock revival is marginal without them.

Outputs of CCSNe simulations and observations.

In the CCSNe simulations, we are certainly interested in whether shock revival obtains and how if it does; however, we are also concerned with various physical quantities they produce, particularly those that can be compared with observations. Such comparisons are indeed mandatory to ascertain the explosion mechanism.

In Fig. 5, we list what we can derive from the CCSNe simulations. As shock revival came to be obtained rather commonly in multidimensional simulations in the past 5 years or so, the focus has shifted to whether the resultant explosion can reproduce observations, particularly the explosion energy. We have been said for many years that the CCSN has a canonical explosion energy of 10^51^ erg. Therefore, it has been a criterion for a successful simulation to yield this amount of energy asymptotically. As we see below, this is changing these days, thanks to the systematic observational surveys of CCSNe (e.g., Ref. 107).

Figure 5. (Color online) Outputs of CCSNe simulations.

After the detection of GWs by LIGO from a neutron star merger, GW170817,^80)^ followed by observations of electromagnetic-wave counterpart, or kilonova, at various wavelengths from radio to gamma rays,^108)^ multimessenger astronomy is a reality now. It is also true for the CCSNe. Among various observables, neutrinos^109)^ and GWs^14)^ will probably be the most important for unveiling the mechanism of CCSNe, as they can penetrate the thick veil of envelopes of massive stars and enable us to directly inspect what is going on deep inside the core.

Both of them will be observed by terrestrial detectors in operation^14,110)^ if a CCSN occurs in our Galaxy now. It has been demonstrated over the years that the detectable levels of GWs will be emitted by core bounce, convections, and SASI at the frequencies from a few 100 Hz to a few 1000 Hz.^14)^ With KAGRA now in collaboration with LIGO and Virgo, the four detectors could find even individual polarization components of GWs, which would in turn enable us to reveal, for instance, the rotation of the core.^111)^ In the future, satellite-borne GW detectors^112,113)^ may have a chance to detect low-frequency (∼ deci-Hz) GWs, such as those emitted by neutrinos that are radiated from PNS anisotropically.^114,115)^

Neutrinos will be also observed by multiple detectors.^116)^ The light curves and spectra for different species will be useful to test numerical models, attest the occurrences of hydrodynamical instabilities, and deepen our understanding of neutrino oscillations, particularly those induced by neutrinos themselves. If a CCSN occurs in our vicinity (≲1 kpc), a chance exists to detect the presupernova neutrinos, i.e., those neutrinos emitted from the progenitor before the core collapse, which will help us characterize the progenitor.^117)^ In fact, we could detect neutrinos from Betelgeuse a week earlier than its supernova explosion if it exploded at all.

The simulations from the core collapse through core bounce to shock revival cannot be extended much beyond several seconds as they are because of a severe limitation to the time step at the center. The dynamical timescale of PNS is τdyn ∼ (Gρ~s~)^−1/2^ ≲ ms, and the numerical time step is actually smaller than this by a few orders at least. As a result, the subsequent long-term evolution up to the shock breakout of the stellar surface and beyond is normally computed by excising the inner region including PNS and replacing it with a certain boundary condition. With successful shock-revival models at hand these days, the latter simulations can be done with the results of the former simulations used for their initial conditions (e.g., Ref. 32). They produce nonspherical explosion morphologies as observed (see, e.g., Ref. 118) as a consequence of hydrodynamical instabilities induced in the envelope by ejecta from the core that are already asymmetric, as explained earlier. The simultaneous computation of explosive nuclear burnings in the envelope provides the nucleosynthetic yields as well, which are also observables and on which what happened deeper inside may have left some imprint.

More simplified simulations, where the explosion is parametrically induced by hand, are commonly done to fit observed data (see, e.g., Ref. 119). They are useful to deduce the parameters that characterize the progenitor and explosion, that is, the mass of ejecta, the explosion energy, the mass of ^56^Ni synthesized, and so on. Thanks to the recent systematic surveys of supernovae (e.g., Refs. 119–121) and these calculations, we now know the systematics of CCSN explosions. For example, the explosion energy is pretty diverse and there is a general trend that the more massive the progenitor is, the greater the explosion energy is. The same is true for the nickel mass (see Fig. 6). It is now insufficient just to ask if the “canonical” explosion energy of 10^51^ erg is attained in the simulation or not. It is more appropriate and probably mandatory to see if these distributions are reproduced or not.

Figure 6. (Color online) Explosion energies and nickel masses deduced from observations. Reproduction with permission of Fig. 8 in L. Martinez et al., “Type II supernovae from Carnegie Supernova Project-I. III. Understanding SN II diversity through correlations between physical and observed properties”, Astronomy and Astrophysics 660, A40, (2022). ©ESO.

Observational information on the remnant, normally a neutron star, is also useful. It includes the mass, spin, proper motion, and magnetic field of neutron star^#21^ as well as the morphology of ejecta, to list the most important ones. The recent simulations are getting capable of giving these quantities quantitatively indeed.

Current state of the art of numerical modeling

Brief summary of 3D simulations.

In this section, we are now ready to look at the current status of the CCSN modeling. In the past decade, we have seen numerous 3D simulations without any preimposed symmetry published in the literature. They produced shock revival (and hence explosions later probably) more often than not, the fact that led many to think rather confidently that the neutrino heating is the mechanism behind most of CCSNe.

Here, we show, as an example of the recent 3D simulations, the results of Princeton University group. The purpose of this section is to provide readers a quick grasp of the sophistication of the latest 3D simulations, and this choice is rather arbitrary. In fact, a couple of other groups exist that have published their own 3D simulations more or less on the same level as those presented here, some being even better on some of the ingredients incorporated (see, e.g., Refs. 39, 123–127 and references therein). They would serve the same purpose equally well. These simulations commonly used a high-resolution shock-capturing scheme based on an approximate Riemann solver for hydrodynamics; the neutrino transfer equation approximated to O(v/c) was solved with the moment method truncated at the first moment; the Newtonian gravity was modified in its monopole component for general relativity; some of the modern nuclear EOSs mainly based on the relativistic mean field theory was adopted; as for the neutrino-matter interactions, the new standard set was implemented (see Section 2.3.7); the turbulence in the progenitor (see Section 2.3.4) was also considered in some cases. Those who are interested in the latest results of these other simulations should consult the original papers or the reviews for experts^6–8)^ and references therein.

Figure 7 shows the temporal evolutions of the average shock radius, diagnostic explosion energy, PNS baryonic mass, and radius for a suite of progenitor models specified by their ZAMS mass, as shown in the panels.^39)^ Although these progenitors are nonrotating, the dynamics are not spherically symmetric because of hydrodynamical instabilities instigated by the perturbations added by hand to the progenitors initially (see Sections 2.2 and 2.3.1). In 3D, even rotation is induced differentially by the spiral modes of SASI. The total angular momentum is, of course, vanishing owing to its conservation. Consequently, the shock wave is not spherically symmetric either. That is why its average is shown in the top panel of the figure. In these simulations, the M1-closure relation was adopted; only the monopole gravity with the general relativistic correction was incorporated; SFHo EOS^83)^ was used for the nuclear EOS.

Figure 7. (Color online) Results of recent 3D simulations. From top to bottom, the average shock radii, diagnostic explosion energies, PNS masses, and radii are shown as functions of time from the bounce for various progenitors with different ZAMS masses. Reproduction with permission of Figs. 2, 6, and 7 in A. Burrows et al., “The overarching framework of core-collapse supernova explosions as revealed by 3D fornax simulations”, Monthly Notice of Royal Astronomical Society 491, 2715–2735, (2020).

In the top panel of Fig. 7, shock revival occurs for most models after relatively short periods (≲300 ms) of stagnation at r ∼ 200 km, triggered by the hit of the boundary of the Si and O layers. The exceptions are those models with the ZAMS masses of 13, 14 and 15 M⊙ . These unsuccessful models have commonly a smaller density jump at the Si/O interface than others. One recognizes the trend that lower-mass progenitors experience shock revival earlier when they do have one. On the other hand, they achieve lower energies by the time when the simulations are terminated, as exhibited in the middle panel of the figure, which shows diagnostic explosion energy. The explosion energy as we observe it cannot be obtained with these short simulations,^#22^ as demonstrated shortly in the next section.

The diagnostic energy is defined at each time normally as the sum of kinetic, internal, and gravitational energies of the matter that have positive values of the sum of these three energies. It changes in time and asymptotes either the explosion energy if an explosion occurs indeed or zero if it fails. Therefore, it is encouraging that the diagnostic energy is positive and increasing toward the end of the simulation as it is really the case for the models with shock revival. However, the nonvanishing diagnostic energies obtained at the end of the simulations are O (10^50^) erg, i.e., one order of magnitude smaller than the canonical explosion energy 10^51^ erg. This is an issue common to most of the latest 3D simulations conducted up to t ∼ 1 sec. As mentioned in Section 2.4 and will be discussed in the next section, this may not be a problem.

The bottom panel of Fig. 7 shows the time evolution of the baryonic mass (and radius) of the region with ρ ≥ 10^11^ g/cm^3^, which is a convenient surrogate of PNS. The former increases in time as matter accretes but are supposed to be settled to the PNS mass after it is isolated from the expanding ejecta. It seems that the baryonic mass approaches the asymptotic value faster than the diagnostic energy does. This PNS mass is defined also for unsuccessful models, and it continues to increase toward the end of their simulations and shows no hint of leveling off as expected for the continuous accretion. The asymptotic PNS masses of the successful models, on the other hand, seem to range from 1.3 M⊙ (for the model with the ZAMS mass of 9 M⊙ ) to 2.0 M⊙ (for the model with the ZAMS mass of 25 M⊙ ).

The ejecta mass was estimated to range from 0.05 M⊙ to 0.15 M⊙ . Of course, one needs to add the mass of the envelope not included in the computational region of these simulations to find the actual ejecta masses. As the nuclear network was not implemented in these simulations, unfortunately, the mass of ^56^Ni was not derived (but see Ref. 128 for the postprocess estimate). However, they found that the electron fraction Ye of the ejecta, an important indicator for the possibility of the rapid neutron-capture process (r-process), one of the two major channels to produce elements heavier than the iron-group elements, is close to 0.5 or higher. This is consistent with the observations over the years that the CCSNe ejecta are not neutron rich and will not be suitable for the r process (see Ref. 129 and reference therein).

In the early phase of the progress in the 3D CCSNe simulations, there were some confusions on the issue of whether the 3D dynamics is more favorable for shock revival than the 2D dynamics in axisymmetry. The nonradial motions associated with the convection and SASI are restricted by the symmetry in 2D and hence qualitatively different in nature from the 3D counterpart. As mentioned earlier (see Section 2.3.1), the most remarkable is the opposite directions of the cascade in the inertial range of in 2D the inverse cascade tends to produce large structures. It is naively supposed then that this property of the 2D dynamics tends to overestimate the role of the turbulence and, as a result, the neutrino heating, and that shock revival is easier in 2D. That was actually not the case in some simulations.^130)^ Although the subsequent resolution studies made the situation more complicated initially, it seems that rather limited numerical resolutions in the early 3D simulations combined with different coordinates used in those computations are partially responsible for the controversies (see, e.g., Ref. 131). For the moment, it is certain that the 2D and 3D dynamics are both advantageous for shock revival than the 1D dynamics in spherical symmetry. It may be interesting that Burrows et al.^39)^ found in most cases that their 3D models explode when the 2D models of the same progenitors do.

In addition to the resolution dependence, the sensitivities of the simulation results to various ingredients, such as the microphysical inputs,^87,106)^ have been also investigated.^6–8)^ It seems that there is no single decisive factor for the successful shock revival. In fact, a turnoff of one of seemingly small corrections can sometimes turn an otherwise successful model into failure (see Section 4). This implies probably that the multi-D dynamics is so close to the criticality (see Section 2.2) that the accumulation of small effects enables to go beyond the threshold. It also means that everything has to be modeled well.

The effect of rapid rotation is not simple. The response of the centrifugal force against compression is similar to that of a gas with the adiabatic index of 5/3. Its extra support tends to expand the shock wave and hence the gain region, which will be helpful for shock revival. This is vindicated by the toy model analysis of the critical luminosity based on steady neutrino-irradiated rotational accretion flows through a standing shock wave.^132)^ Counteracting effects always exist the centrifugal force also supports PNS, thereby lowering its density and temperature, which results in the reduction of the neutrino luminosity and energy.^133)^ The SASI is also affected by rapid the nonaxisymmetric spiral mode with m = 1 in the spherical harmonic expansion becomes dominant over other modes such as the sloshing one with m = 0, which is prominent in 2D with axisymmetry.^134)^ It was pointed out^135)^ that low-T/|W| instability can be induced. The consequence of these instabilities is often the formation of one or two arms, producing shock revival eventually in the equatorial plane^133,135)^ (but see also Ref. 136 where the explosion occurred without help from these instabilities and the resultant bipolar geometry is roughly aligned with the rotation axis). We pointed out its possible imprints in the gravitational wave signals^111)^ already in Section 2.4. It is also worth a mention that the rotation of the arms will induce a characteristic temporal modulation in the neutrino luminosity that may be observed on earth.^137)^

Long-term multi-D simulations.

After 3D-simulation uproar, one of the important progresses is the long-term multi-D simulations.^7,38,138)^ By the long-term simulation, we mean here the one well beyond a second postbounce, i.e., up to the time when the explosion energy is finally determined. One should be able to understand how tough this will be if one recalls that the typical time step is a few orders shorter than the dynamical timescale of PNS, i.e., ≲1 ms. The 3D simulations up to a second postbounce are already a remarkable achievement. Those much beyond it should be highly appreciated.

Shock revival is certainly a necessary condition in the neutrino-heating scenario to obtain CCSN as we observe it but not a sufficient one. In fact, as the top panel of Fig. 7 shows clearly, the success or failure of shock revival can be reliably judged by the postbounce time of ≲1 sec; the middle panel demonstrates, on the other hand, the (diagnostic) explosion energy at this time is still much smaller than the canonical explosion energy of 10^51^ erg. It has been often argued that if it is extrapolated to, say, 10 seconds, it will reach the canonical value. Provided that the CCSN dynamics is highly nonlinear, this argument cannot be easily accepted. To settle this issue, much longer simulations were needed.

In 2D, such results were reported rather recently.^7)^ They are simulations under axisymmetry up to ∼5 sec postbounce for a suite of progenitor models with the ZAMS mass of 9–27 M⊙ , with essentially the same physics incorporated as for their 3D simulations presented in the previous section. Again, most models give shock revival successfully. This time only those models with 12 M⊙ and 15 M⊙ are unsuccessful. There is a general trend for the successful models that the more massive the progenitor is, the later shock revival occurs. The comparison with the 3D counterparts demonstrates that the 2D models tend to give earlier shock revivals.

More important is the temporal evolutions of the diagnostic explosion energy for these models. The diagnostic explosion energy increases almost linearly with time indeed and that some models have reached or even exceeded the canonical explosion energy in a few seconds. For other models the diagnostic explosion energy is leveled off at smaller values. It actually ranges from 0.09 × 10^51^ erg to 2.3 × 10^51^ erg at the end of the simulations. Here, there is a trend that the more massive the progenitor is, the higher energy it yields. The evolutions of PNS mass for the same models, on the other hand, confirm the earlier observation (Section 3.1) that it reaches the asymptotic value more quickly than the diagnostic explosion energy in the successful models.

The most impressive of all is probably the explosion energies as a function of the ejecta mass compared with the corresponding observational data^#23^ (see Fig. 6 in their paper). The numerical results nicely mesh with the observations. As the ejecta mass is the total mass minus the PNS mass and the former overwhelms the latter normally, the trend is essentially the same as the one mentioned the more massive the progenitor is, the higher the explosion energy is. The explosion energy inferred from the observations covers a wide range with under- and over-energetic explosions existing. The canonical energy of 10^51^ erg seems rather useless in judging the validity of numerical models. It is probably these kinds of observational data that should be reproduced by successful CCSNe simulations.

Why does the diagnostic explosion energy continue to increase in time for such a long time? The naive expectation is that once shock revival occurs, the accretion of matter through the shock wave down to PNS is inverted to expansion; then, the neutrino luminosities will drop, as the contribution from the matter accretion is lost, and the neutrino-heating is essentially turned off. This is not what happens in the multi-D dynamics. In fact, the accretion continues from certain directions after the shock wave revives and the matter starts to expand in other directions. Such cohabitation of accretion and expansion is possible in multi-D and is responsible for the continued neutrino heating and the lasting gain of the explosion energy.

Essentially the same picture was obtained in the 3D simulation in Ref. 38 (see also Ref. 128). This is a computation for a single progenitor model with the ZAMS mass of 19 M⊙ . The simulation was initiated at 7 minutes before the core collapse and continued up to 7 sec postbounce. The neutrino transfer is calculated for all six neutrino species with ray-by-ray plus approximation, in which spatially one-dimensional truncated moment equations are solved for each radial direction individually, with the nonradial neutrino advection neglected except for the contribution from lateral matter motions; the closure relation is derived by solving simplified Boltzmann equations. They adopted rather conventional Lattimer–Swesty’s (LS) EOS, which was derived from a parametrized free energy based on Skyrme interaction^140)^; the incompressibility at the saturation density is K = 220 MeV; the ensemble of nuclei is replaced by a single representative heavy nucleus. More importantly, a nuclear network calculation, albeit small with 23 nuclei, was incorporated so that the recombination of nucleons into nuclei is handled, which was lacked in Ref. 7. They replaced these detailed treatments with something simpler after ∼2 sec postbounce, provided that the 3D simulations are much computationally demanding.

In the top panel of Fig. 8, the angle-averaged temporal evolutions of some representative mass shells are shown with black lines, whereas those of the maximum, minimum, and average shock radii are indicated with white dashed, dotted, and solid lines, respectively; the white dash-dotted line marks the position of the gain radius; the color represents the entropy per baryon; the two red lines are the trajectories of the Si/O interface and the outermost mass shell that will be settled eventually on the surface of PNS. Shock revival is instigated when the Si/O interface hits the stagnated (in fact, receding) shock wave.

Figure 8. (Color online) 3D long-term simulation. (Top panel) the time evolutions of some representative mass shells (black lines); the maximum, minimum, and average shock radii (white dashed, dotted, and solid lines, respectively); and the gain radius (white dash-dotted line); the color represents the entropy per baryon. (Middle panel) The diagnostic energy with (red solid line) and without (green solid line) overburden and their time derivatives (red- and green-dashed lines) as well as the neutrino-heating rates (blue and orange solid lines) multiplied with factors indicated as functions of time. (Bottom panel) The time evolution of the abundance of various nuclei. Reproduction with permission of Figs. 1, 2, and 7 in Bollig et al., “Self-consistent 3D supernova models from −7 minutes to +7 A 1-bethe explosion of a ∼19 M⊙ progenitor”, Astrophysical Journal 915, 28 (23pp), (2021), DOI: https://doi.org/10.3847/1538-4357/abf82e. ©AAS.

The middle panel presents the evolutions of the diagnostic energy and the energy deposition rate. It is clear that the former continues to rise for ∼5 sec. The asymptotic value with the overburden corrected exceeds 0.9 × 10^51^ erg.^#24^ The continuous rise is supported by the almost constant energy deposition of ∼0.2 × 10^51^ erg/s from ∼2 to ∼5 sec. They found that the recombinations also contribute to this long-term increase in the diagnostic explosion energy. They emphasize that the long-lasting accretion of matter after shock revival feeds them all. The gain region is replenished with cold accreting matter, which is then heated up by neutrinos and joins the ejecta later. This is also confirmed from the evolution of the PNS mass. Although not shown in the figure, it is already settled at ∼1 sec to the asymptotic value of ∼1.65 M⊙ .

The bottom panel of Fig. 8 displays the temporal change in the nuclear abundance in the ejecta. As mentioned above, it is given by the network calculation up to ∼2 sec and estimated under the assumption of NSE thereafter. Among other things, the abundance of ^56^Ni is the most important. One finds that it increases in time and levels off after ∼5 sec. There is a small jump discernible at the transition from the network calculation to NSE around ∼2 sec post bounce. It is conservative to say that the asymptotic value is MNi56∼0.05M⊙ . Although this is slightly smaller than the value of 0.07 M⊙ estimated for SN1987A, which occurred in the Large Magellanic Cloud in 1987 and is the only one so far, for which neutrinos were detected by the terrestrial detectors,^15,16)^ it is well within the range of diverse values inferred from observations of the light curve and spectrum (see Fig. 6).^#25^

The two simulations cited here (see also Refs. 128 and 138) demonstrated clearly that this long computations are indeed needed to predict unambiguously some observable quantities that are critically important to judge if the numerical model is truly successful. So far, the results look very encouraging.

Boltzmann ν–radiation–hydrodynamics simulations.

Although this paper is mainly meant for review on the recent progress in the CCSN theory, this section contains an original piece concerning authors’ we have been working on the multi-D CCSNe simulations that solve the Boltzmann equations as they are.^#26^ We believe that this is an important contribution to the community, as the others are commonly using truncated moment approximation these days and therefore may share the same systematic issues. In the following section, we report what we have found, with some unpublished results included.

Figure 9 displays the time evolution of shock radii for some of our axisymmetric 2D CCSN simulations with full Boltzmann neutrino transport. The color-shaded regions show the ranges of the shock radii. The EOS dependence can be seen by comparing models 11.2LS, 11.2FS, and 11.2VM, in which we use the LS EOS with the incompressibility K = 220 MeV,^140)^ Furusawa (FS) EOS,^141)^ and a multinuclear EOS based on the variational method (VM),^142)^ respectively. For this comparison, we use the 11.2 M⊙ progenitor of Woosley and Heger.^143)^

Figure 9. (Color online) Time trajectories of shock radii for some selected CCSN models obtained by our axisymmetric Boltzmann ν–Radiation–Hydrodynamics simulations. We display the results for 11.2– and 15 M⊙ progenitors with different EOSs (FS, LS, and VM) and a stellar rotation (rot); the color distinguishes these models (see the text for more details). The thick line represents the angular-averaged shock radii, whereas the angular variations are expressed as the shaded regions.

As shown in the figure, the LS model achieves shock revival, whereas the FS model fails. For the VM model, the outcome at the end of the simulation is difficult to judge, as the shock wave has not yet entered a runaway expansion but seems to be still in a large-amplitude oscillation. Provided that the shock radius is quite similar to that in the LS model, the VM model has a similar level of explodability to the LS model. We find that the baryon mass in the gain region and the strength of turbulence vary with models, the fact that accounts for the different shock evolutions.^#27^ The strength of turbulence is associated with the vigor of the prompt convection occurring at ∼10 ms after bounce,^144)^ which depends in turn on EOS through the nuclear composition of accreting matter. In fact, the nuclear binding energy per nucleon is higher for the LS model, implying that the energy loss by photodissociations is large. As a result, the entropy gradient becomes steeper downstream the shock wave and triggers the stronger prompt convection.

The progenitor dependence of shock evolution can be found by the comparison between two VM 11.2VM and 15VM in Fig. 9. The time evolutions of shock radius are almost identical initially but they start to deviate from each other at ∼60 ms, which should be associated with some difference in the structures of the progenitors such as the density jump at the Si/Si–O interface. In fact, an earlier arrival of the interface accounts for the more rapid shock expansion in the VM11.2 model. The density jump at this interface reduces the ram pressure of the accreting matter, facilitating the shock expansion (see Section 2.3.4). Our result is in line with recent studies of the progenitor dependence in the CCSNe simulation by other groups.^145–147)^

The diagnostic explosion energy of the 11.2LS model is less than 10^50^ erg at the end of our simulation.^148)^ However, the value is provisional and a long-term simulation is certainly needed to determine unambiguously the explosion energy that can be compared with observations (see Section 3.2). An earlier shock revival may not always lead to a high explosion energy^138)^ as the mass accretion onto PNS after shock revival is a major source of explosion energy.

Let us give a brief discussion on a possible source of neutron star kick. It is well known that neutron stars are moving at a few 100 km/s on average.^149)^ This is called the neutron star kick, as the linear momentum is supposed to be imparted to the neutron star at birth. In Fig. 10, we show the color map of Ye together with the vectors of fluid velocity in the meridian section at 200 ms for the 11.2VM (left panel) and 15VM (right panel) models; the inner region covering −50 km < x < 50 km and −50 km < z < 50 km is exhibited. It is clear that PNS is displaced in the Z-direction for the 11.2VM model, whereas it remains at the center for the 15VM model. As discussed by Nagakura et al.,^150)^ the PNS proper motion in the 11.2VM model is not driven by the hydrodynamical or gravitational force; in fact, they almost balance each other. Instead, it is the asymmetric neutrino emission, which is associated with a dipolar distribution of Ye near the surface of PNS, that induces this motion. On the other hand, the asymmetric Ye profile is sustained in turn by the PNS proper motion itself. It should be obvious that the proper handling of the feedback loop between the PNS proper motion and the asymmetric neutrino emission is crucially important to describe this proper motion of PNS correctly. Our Boltzmann-neutrino-transfer code is fully general relativistic, based on 3 + 1 formalism, and can treat the PNS proper motion using coordinate transformations.^150)^

Figure 10. (Color online) Color maps of Ye with the vectors of matter velocity in the meridian section. The left and right panels show the results of the 11.2VM and 15VM models, respectively, at 200 ms after bounce. The PNS is clearly displaced from the origin in the 11.2VM model (left panel), which is owing to the asymmetric neutrino emissions (see the text for more details).

The conservation of the linear momentum is a challenge in the Boltzmann–radiation–hydrodynamic simulations, particularly in the diffusion regime. We resolve this issue using a novel method, in which the momentum transfer between neutrinos and matter is treated not with the collision term of the Boltzmann equation, in which severe cancelations occur, but with the advection term in the momentum conservation equations in the diffusion regime.^151)^ This improves the linear momentum conservation in our CCSN simulations, enabling us to compute the proper motion of PNS reliably. The recent study of PNS kick based on their 3D CCSN models supports our claim that the impact of asymmetric neutrino emissions on the PNS kick can be comparable with or even larger than the hydrodynamic contribution.^152)^ These results may change the current theory of the NS kick, although more systematic and longer-term simulations are required to draw a robust conclusion.

It is an interesting question how stellar rotation affects the results and may be inferred from the comparison between 15VM and 15VM-rot in Fig. 9. In the latter model, we add rotation by hand to the nonrotating progenitor according to the following angular velocity

when the central density reaches 3 × 10^10^ g/cm^3^. In this expression, Ω and r denote the angular velocity and radius, respectively. As shown in the figure, the rotation facilitates the shock expansion. This is owing to the centrifugal force, which works to reduce the ram pressure of accreting matter in the preshock region. In fact, the shock morphology tends to be oblate in the 15VM-rot model, which is shown in Fig. 11. The neutrino luminosities and mean energies are substantially lower in the rotational model than those in the nonrotational model, implying that the neutrino heating becomes less efficient by rotation.^133)^ This may be true for rotations slower than those assumed in our simulation.

Figure 11. (Color online) 2D color maps of the entropy per baryon in the meridian section at 200 ms after bounce for the nonrotating (left panel) and rotating (right panel) models of 15 M⊙ progenitor with the VM EOS.

The stellar rotation can provide a preferable condition for the occurrence of collective neutrino oscillations^153)^ (see Section 4.4 for more details). This is mainly owing to the expansion of the low-Ye region near equatorial plane, which tends to make the ν~e~ distribution more isotropic in momentum space through the charged-current reactions than the ν¯e distribution. Therefore, the electron neutrino lepton number (ELN) crossing occurs in the energy-integrated angular distributions of ν~e~ and ν¯e , and the fast-pairwise flavor conversions (FFC) are triggered. This is different from other mechanisms observed in nonrotating models (see, e.g., Refs. 154–157) and may be intrinsic to rapidly rotating CCSNe. The nonlinear evolution of FFC coupled with matter dynamics should be investigated in the future.

Using the neutrino data obtained in the Boltzmann-neutrino-transfer simulations, we can make an assessment of the approximate transfer methods used in the literature (see Sections 2.3.3 and 4.2). This inspection also highlights the strengths and weaknesses of each method (including the discrete ordinate method for our Boltzmann solver). As we have discussed in the previous Section 2.3.3, the truncated moment method is a modern approach that is most popular in the recent CCSNe simulations and needs scrutiny of their fidelity. Here, we pay a particular attention to the closure relations. As the angular distribution in momentum space is calculated in the Boltzmann neutrino transfer, it can provide all components of the second angular moment (or the Eddington tensor), which can be used for comparison with the ones given by the closure relations. They may be used to develop better closure prescriptions (see Ref. 158 for an application to the deep learning of the closure relation).

The sensitivity to the closure relation used of the resultant neutrino distributions has been studied by various authors (see, e.g., Refs. 49, 50, 52, 159 and Section 4.2). Although some of these studies suggest that the result is insensitive to the choice of closure relation, we had better be cautious, as the accuracy of closure relations depends on the background matter profile.^160)^ In strongly nonspherical matter distributions, nonradial neutrino fluxes can be comparable with the radial flux but may not be accurately handled by the analytical closure relations. As such an example, we present here a comparison of the r − θ component of Eddington tensor directly derived in a Boltzmann-radiation-hydrodynamics simulation and that given by the M1 closure relation for the 15VM-rot model (see Fig. 12).

Figure 12. (Color online) Color maps of the r − θ component of Eddington tensor for the 15VM-rot model in the meridian section at 12 ms after bounce. They cover the spatial region of 0 km < x < 150 km and −150 km < z < 150 km. The left panel represents the result directly computed from the neutrino distribution functions obtained in the neutrino-radiation-hydrodynamics simulation with full Boltzmann transfer; the middle panel shows the result obtained with the M1 closure relation from the zero-th and first angular moments that are calculated from the neutrino distribution functions obtained in the same neutrino-radiation-hydrodynamics simulation; the right panel displays the difference between the two results.

As shown in this plot, the difference is large indeed in the optically thin region. Importantly, the r − θ component of the Eddington tensor is related with the time evolution of lateral neutrino advection (see, e.g., Ref. 144). The result indicates that the distribution computed by the truncated moment method with the M1 closure may have systematic errors. However, the zero-th and first moments used in the closure relation in this study are derived from the same Boltzmann-radiation-hydrodynamics simulation by appropriately integrating the angular distribution. More quantitative assessments of the difference between the Boltzmann and truncated moment methods are currently underway (see also the principal-axis analysis below). These efforts will hopefully provide hints to eventually improve the closure relation.

The second moment discussed above (not to mention the zero-th and first moments) is just one of low-order moments. Let us stress again that the Boltzmann-transfer simulation grants us full access to the angular distributions of neutrinos in momentum space. In Fig. 13, we show such distributions of ν¯e obtained in the VM15-rot model at an early postbounce time of 12 ms. The angular distribution depends on the frame of reference. The figure presents the distributions for different energies of ν¯e in the laboratory frame at the spatial point of r = 57 km on the equatorial plane. They are plotted as surfaces in such a way that the distance from the origin to each point on the surface is proportional to the the distribution function for the direction to that point.

Figure 13. (Color online) Angular distributions of ν¯e in momentum space at r = 57 km on the equatorial plane at 12 ms after bounce for the 15VM-rot model. From left to right, the neutrino energy increases from 4 MeV to 19 MeV to 42 MeV. From top to bottom, we change the viewing angle.

For the low-energy (4 MeV) neutrinos displayed in the left column, the angular distribution is radially forward peaked, indicating the matter is nearly transparent to those neutrinos at that point. As the neutrino energy increases, the backward-propagating neutrinos get populated, resulting in more isotropic distributions. However, they are not perfectly isotropic but are nonaxisymmetric in fact as should be evident from the plots in the middle and right columns. This is owing to the rotation of matter. The breaking of axisymmetry in the neutrino angular distribution implies that the azimuthal component of neutrino flux is nonvanishing, which is but a natural outcome for the neutrino emission from a rotating matter. It also means that neutrinos carry away angular momentum from PNS.^161)^ These complex but characteristic features in the neutrino angular distribution will be useful and basic data to develop a better closure relation for the truncated moment approximation.

Last but not least, we describe two recent progresses in our Boltzmann-radiation-hydrodynamics code. One of them is the general relativistic extension. Although our Boltzmann-transfer code was fully general relativistic from the beginning, based on the conservative form given in Ref. 162 (see also Ref. 163), we have turned off general relativistic terms in most of our CCSNe simulations.^#28^ Recently, however, we conducted a suite of test calculations in black hole spacetimes to demonstrate the general relativistic capability of our Boltzmann solver.^164)^ More recently, we applied the code to the computation of PNS convection in a general relativistic but fixed spacetime.^165)^ We are currently working on the integration of the metric solver with our Boltzmann solver, which will enable us to perform multidimensional Boltzmann-radiation-hydrodynamics simulations in dynamical spacetimes.

The other progress worth a mention is that 3D Boltzmann-radiation-hydrodynamics simulations are now getting available. We reported the first result paying attention to the early postbounce phase in Ref. 166 (see Fig. 14). Herein, we also presented an in-depth analysis of the Eddington tensor, using a new approach, i.e., the principal-axis analysis. Put succinctly, it enables us to visualize the Eddington tensor as an ellipsoid whose principal axes are parallel to the eigenvectors of the Eddington tensor. We made an assessment of the ability of the M1 closure relation to reproduce the ellipsoid. It assumes the longest principal axis to be always aligned with the neutrino energy flux. However, this was not always true (see Ref. 166 for more details). This is an important observation, suggesting a possible direction of improvement for the closure relation, although more studies are certainly needed. Collaborations with other groups are crucially important and the entire community will benefit from such joint efforts.

Figure 14. (Color online) Vectors of the number flux of ν~e~ normalized by its number density, or the flux factor, and the color map of the number obtained by 3D ν-radiation-hydrodynamic simulations with full Boltzmann neutrino transport. The snapshot is taken at 11 ms after bounce.

Critical assessments and remaining issues

In the previous section, we described the achievements that the CCSNe modelers have accomplished so far, paying particular attention to the latest progresses. They are remarkable various ingredients in macro- and microphysics have been scrutinized and their implementations in the numerical code have been continuously improved. The results seem to be very well compatible with observations. They may be too good actually, provided there are still unsatisfactory issues remaining. In this section, we will look at them rather critically.

Turbulence.

The hydrodynamical instabilities are ubiquitous in the quasi-static evolution of progenitors and in the subsequent development of CCSN explosions. The turbulence is a common outcome in their nonlinear phase and, as mentioned in Sections 2.2, 2.3.1, and 2.3.4, it is supposed to play a crucially important role in the neutrino-heating mechanism, pushing the shock wave outward with its anisotropic Reynolds stress. Therefore, it is understandable that its proper handling in the simulation is equally crucial. In the case of 3D turbulence, eddies are cascaded into smaller ones and the numerical resolution becomes an issue. In the energy spectrum of the turbulence, the inertial range is extended down to the scale of dissipation, which is essentially of the order of the mean-free path of constituent particles and is hence smaller by many orders of magnitude than the scale that numerical simulations can normally resolve. Then, the issue is whether this underresolution has an influence on numerical results.

In fact, in the early development of the 3D simulations, the identification of the inertial range in the simulations was controversial and may have been a source of confusions on the role of the turbulence in the neutrino-heating mechanism.^23)^ In this respect, the bottleneck effect, that is, the slowdown of the energy cascade by viscosity that gives rise to the accumulation of the turbulence kinetic energy and to the flattening of the turbulence spectrum toward the end point of the inertial range, is particularly important.^167)^ A more thorough investigation on this issue in Ref. 168 demonstrated that the bottleneck appears indeed in the simulation if the resolution is high enough; by contrast, if the simulation is underresolved, the bottleneck becomes unrecognizable, being merged with the inertial range under development, and may even change the slope of the spectrum (see Fig. 15). Interestingly, the bottleneck effect occurs on the dissipation scale of the simulation, i.e., roughly the cell size in the numerical mesh, which is much greater than the real dissipation scale.

Figure 15. (Color online) The energy spectra of turbulence obtained in the simulations with different spatial resolutions. The bottleneck effect emerges as a hump at the right (high wave number) end of the inertial range, i.e., the flat portion of the spectrum. Reproduction with permission of Fig. 13 in D. Radice et al., “Neutrino-driven convection in core-collapse high-resolution simulations”, Astrophysical Journal 820, 76 (19pp), (2016). DOI: 10.3847/0004-637X/820/1/76. ©AAS.

The emergence of the bottleneck indicates that the formation of the inertial range is captured in the simulation, which is certainly great but actually a necessary condition. It does not guarantee that the resolution is sufficient. Unfortunately, even this necessary condition is rarely demonstrated in recent 3D simulations because of the severe limitations in the numerical resources available. The goodness of numerical resolution is normally justified by the robustness of the results, e.g., whether shock revival occurs, against the change in resolution. One may argue that the global dynamics such as shock revival converges faster than the turbulence features with the increasing number of grid points. However, the resolution test is rather limited in those simulations and the same argument was also made in previous papers. One may find a qualitative difference if the numerical resolution is changed by orders.

Provided that it is all but impossible for the numerical simulation to resolve all features in turbulence down to the dissipation scale, it may be an alternative idea to construct a subgrid model of the turbulence effects and implement it in the numerical code for the CCSNe simulations. The model should be not only physically motivated but also well calibrated, for example, by local high-resolution simulations. The investigation in that direction is much awaited.

Neutrino transport.

The numerical treatment of neutrino transfer is one of the key ingredients in the CCSNe modeling. This could not be overemphasized and indeed has been vindicated over the years as the numerical method is improved continuously and the spatial dimensions are extended from 1D to 2D, and then to 3D. In fact, the very first simulation by Colgate and White^2)^ that envisaged the neutrino heating as the agent of explosion replaced neutrino transport with a local energy deposition and easily produced an explosion. The necessity of the neutrino transfer calculation was soon realized.^169)^ After the CCSNe scenario based on the electroweak theory was established,^17,170)^ Livermore group successfully produced explosions via the neutrino-heating mechanism in their 1D simulations with the flux-limited diffusion approximation for neutrino transfer^171)^ during 1980s. The success is now thought to be owing to their assumption on the doubly diffusive instability (also called the neutron-finger instability),^#29^ which is thought at present unlikely to occur in the supernova core.^172)^ In fact, by the early 1990s, unsuccessful explosions had been consistently reported in spherically symmetric simulations with the flux-limited diffusion approximation for neutrino transfer.^173)^ This conclusion was finally confirmed with the Boltzmann neutrino transport^174–177)^ by the early 2000s.

Back in the middle of 1990s, some 2D simulations that used light-bulb approximation for neutrino transfer in the optically thin region^#30^ produced explosions^178–180)^; the optimism fizzled soon, however, when more elaborate neutrino transport based on the ray-by-ray extension of the sophisticated transport method used in the 1D simulations^181)^ or the flux-limited diffusion^182)^ yielded no explosion by neutrino heating again in the first decade of 2000s. Into the second decade, however, successful 2D simulations started to be reported (see, e.g., Refs. 139 and 183 but also, e.g., Ref. 184 for the discrepancies among different groups at that time). The subsequent improvements in the treatment of various physics such as the turbulence, general relativistic gravity, neutrino physics (see Section 4.3) and, among other things, the spatial dimension in the simulation eventually led to successful explosions in the past 5 years or so (see, e.g., Refs. 6–8).

The current de facto standard for neutrino transfer is the two-moment method with some closure relation prescribed (see Section 2.3.3). It is truly a multidimensional transport method unlike the ray-by-ray treatment. There are various closure relations proposed in the literature (see Ref. 49 and references therein) and comparative studies have been done.^50,159)^ However, most of these comparative studies were done for the stationary spherical background (but see Ref. 52). It may be particularly reassuring that it was demonstrated very recently^159)^ that the outcomes of the multi-D simulations with different closure relations produced similar results, with the differences smaller than those from the uncertainties in EOS (see Fig. 16). However, it is a comparison among different closures, and they are all approximations. What is really needed is the comparison with the exact result, which is unfortunately inaccessible for the moment. In the same vein, it is a bit concerning that most of the recent simulations have used the two-moment method, although the closure relation is different, as they may have the same systematic issues. Simulations based on other formulations are much awaited.

Figure 16. (Color online) Comparison of closure relations. From top to bottom, the shock radii, heating rates, neutrino luminosities, and mean energies are shown. Reproduction with permission of Fig. 16 in T. Wang et al., “Effects of different closure choices in core-collapse supernova simulations”, Astrophysical Journal 943, 78 (18pp), (2023). DOI: https://doi.org/10.3847/1538-4357/aca75c.

The CCSNe simulation with the Boltzmann neutrino transport should play that role. Direct comparisons of our simulation outcomes are currently hampered by the facts that the neutrino processes incorporated in our Boltzmann solver are rather limited compared with the most advanced codes with the two-moment method and that the general relativistic corrections to gravity are still not considered.^#31^ We hence investigated the performance of some closures, using snapshots from our 2D and 3D simulations (see Section 3.3) and demonstrated that the neutrino distribution in momentum space has the longest principal axis neither aligned with nor perpendicular to the flux as in the M1 closure relation. The implications of this finding for the CCSNe simulation are remaining to be studied.

In addition, our current multi-D Boltzmann simulations have a resolution issue. They are computationally demanding even under axisymmetry, in which the total dimensions are 5 (= 2 (in space) + 3 (in momentum space)). Consequently, the number of the grid points deployed in momentum space is rather small (typically 10 in θν and 6 in ϕν). Although the neutrino distribution is rather featureless at small angular scales in momentum space, the globally forward-peaked distribution expected in the optically thin region is underresolved in general. In fact, a comparison with a Monte-Carlo code^52)^ demonstrated that much larger number (∼40 in θν) of grid points are needed to fully resolve such distributions (see Fig. 17). Although the angular distribution is not so forward peaked in the gain region, this issue is certainly a concern, which needs to be addressed.

Figure 17. (Color online) Comparison of Boltzmann (DO) and Monte-Carlo codes. The neutrino angular distributions in momentum are compared. The colors indicate the spatial positions, whereas the line thickness shows the resolutions, with the thin lines corresponding to the standard resolution and the medium and thick lines to twice and four times higher resolutions. Reproduction with permission of Fig. 5 in S. Richers et al., “A detailed comparison of multidimensional Boltzmann neutrino transport methods in core-collapse supernovae”, Astrophysical Journal 847, 133 (21pp), (2017). DOI: https://doi.org/10.3847/1538-4357/aa8bb2. ©AAS.

Many-body effects.

We saw in Section 2.3.7 that the currently most advanced simulations incorporate various neutrino processes in detail, which should be highly appreciated. However, not all implementations are of the same quality (see the followings for details). This is understandable if one considers difficulties to theoretically estimate some effects reliably. The nuclear many-body effects are probably one of such things.

Unbound protons and neutrons are the most abundant constituent particles in the postshock matter. They are the most important reactants for neutrinos. As we touched in Section 2.3.7, nucleons are themselves interacting with one another via strong nuclear forces. These interactions modify the dispersion relations of nucleons and their responses to neutrinos via weak interactions. If the nucleons are spatially and/or temporally correlated or anticorrelated in their density and/or spin density owing to the mutual interactions, their interactions with neutrinos via the vector and/or axial vector current will be enhanced or suppressed, respectively.^9,185)^ An extreme example of such a spatial correlation is, as mentioned earlier, the coherent scattering of neutrinos on the nucleons bound in heavy nuclei, in which the cross section is enhanced if the wavelength of neutrino is longer than the nuclear radius. At supranuclear densities in PNS, by contrast, the nuclear interactions are dominantly repulsive, producing anticorrelations among nucleons, which tend to reduce their weak interaction rates in general. The scatterings via the axial vector current are normally dominant over those via the vector current; the spin correlation may be more important.^101)^ It should be reminded that not only the spatial correlations but also the temporal correlations are important.^9)^ Intuitively, if a nucleon has its spin flipped by another nucleon during its interaction with a neutrino, the latter reaction will be affected.

These correlations are difficult to theoretically calculate. The nuclear interaction is strong and does not allow perturbative treatments. As explained in Section 2.3.7, the recent advanced CCSNe simulations consider the many-body effects at high densities (ρ ∼ ρ~s) based on the RPA with the Fermi liquid theory, in which the nuclear interactions are described by Landau–Migdal parameters and the effective mass of nucleon is considered^185)^; at low densities ρ ≲ 10^12^ g/cm^3^, the Virial expansion is adopted, and EOS as well as the response to weak interactions of nucleons in the long-wavelength limit are calculated on the same basis.^101)^ This is nice because it is model-independent. As the Virial expansion is an expansion with respect to the fugacity z=eμN/kBT (or the density more intuitively), it is valid at low densities. In this expression, μN~ is the chemical potential of nucleon. The RPA calculation, on the other hand, is supposed to be valid at high densities as the parameters used are suitable at ρ ∼ ρ~s~.

In fact, these two theoretical estimates, one from the Virial expansion and the other from RPA, are incorporated in the actual CCSNe simulations, somehow being extrapolated to ρ ∼ a few × 10^12^ g/cm^3^ and merged (see Fig. 18).^101)^ An artificial nature of this prescription is particularly evident at low values of the proton fraction Yp. In fact, the validity of the extrapolations themselves is suspicious in the first place. The Landau–Migdal parameters used in Ref. 185 are calibrated at the saturation density^101)^ and will not be appropriate at such “low” densities; in addition, as shown in Ref. 186, the postshock matter is more like a classical liquid rather than the Fermi liquid in the relevant region. It is also an issue whether RPA properly accounts for the relevant correlations at these densities. On the other hand, the fugacity, the expansion parameter in the Virial expansion, gets larger with density and is thought to be very close to the validity limit at these “high” densities, and therefore, the accuracy of the expansion may be a concern; moreover, the momentum transfer dependence, which is ignored in the long-wavelength limit considered in the calculation of the weak interaction rates in the Virial expansion,^101)^ may be actually important.^187)^ These intermediate densities ρ ∼ 10^11^–10^12^ g/cm^3^ are crucially important for the neutrino-heating mechanism, as they correspond to the region, where the neutrino luminosity and spectrum is formed.

Figure 18. (Color online) Interpolated many-body effects. The vertical axis represents the axial response function, i.e., the static spin-structure function evaluated in the long-wavelength limit and normalized by the total number density of nucleons so that it should be dimensionless. See Ref. 185 for more details. Reproduction with permission of Fig. 2(b) of C. J. Horowitz et al., “Neutrino-nucleon scattering in supernova matter from the virial expansion”, Physical Review C 95, 025801, (2017), DOI: https://doi.org/10.1103/PhysRevC.95.025801. ©APS.

There is an indication from the CCSNe simulations (see Ref. 87 and references therein) that the nuclear many-body effects have a nonnegligible influence on shock revival. In Fig. 19, an example of the comparison between the runs with and without the many-body effects is presented. It is found in the models with a 16 M⊙ progenitor that shock revival does not obtain without the many-body effects while it does with them. As mentioned in Section 3, there seems to be no single decisive factor for shock revival in the neutrino-heating mechanism, but the accumulation of various effects amounts to the successful shock revival in the end. The nuclear many-body effects are one of such small effects. In fact, the model for a 10 M⊙ progenitor fails to produce shock revival, even when the many-body effects are incorporated, whereas the model for an 11 M⊙ progenitor leads to shock revival irrespective of the many-body effects.^188)^ Therefore, they are certainly important but not decisive on their own.

Figure 19. (Color online) Contribution of the many-body effects. A comparison of the simulations with different microphysics incorporated are presented. Attention should be paid to the difference between the black (with the many-body effects) and green (without the many-body effects) lines. Reproduction with permission of Fig. 1 in A. Burrow et al., “Crucial physical dependencies of the core-collapse supernova mechanism”, Space Science Reviews 214, 33, (2018).

Last but not least, the weak interactions with nuclei may be regarded as one of the nuclear many-body issues. No doubt the electron capture on nuclei is the most important, dictating the electron fraction evolution in the prebounce phase. Although a remarkable progress has been made by nuclear physicists on this issue over the years (see Ref. 104 and references therein), the current rates will not be the final ones.^105)^ The neutrino absorption by heavy nuclei is normally incorporated in the CCSNe simulation as the inverse process of the electron capture but the absorption of the electron-type antineutrinos is ignored then. There is also a room for improvement in the treatment of neutrino reactions on light nuclei.^189)^ In fact, the light elements, such as deuteron, triton, hellion, and α particle, appear, although not so abundantly, in the postshock region. They also absorb and scatter (inelastically) neutrinos, thus contributing to the neutrino heating. We had better not forget the uncertainties in these reactions.

Neutrino oscillations.

Neutrino oscillations are flavor-changing propagations of neutrinos induced by their dispersion relations at low energies, which are not diagonal to the flavor eigenstates. The neutrino masses have nonvanishing off-diagonal components in flavor and induce flavor conversions during the propagation in vacuum. In matter, neutrinos interact with other particles. Such interactions modify the dispersion relations of neutrinos in general. For example, the forward scattering of electron-type neutrinos on electrons (the MSW effect) in the sun is believed to be the solution of the solar neutrino problem. The possible implications of these types of neutrino oscillations for CCSNe and its neutrino signals observed on earth have been extensively studied (see Ref. 190 for a recent review and references therein). In fact, the MSW effect occurs in the envelope of the progenitor and does not affect the explosion mechanism. It is mandatory, however, to consider this effect in the study of signals at terrestrial detectors.

The current focus is on the neutrino oscillations induced by the interactions among neutrinos themselves (see Ref. 10 for a recent review and references therein). Such interactions also give nonnegligible modifications to the dispersion relations of neutrinos if they are abundant. It is then obvious that the neighborhood of the PNS produced in CCSNe provides such an environment. What differentiate neutrino-induced neutrino oscillations from the other variants mentioned above is the fact that they are nonlinear phenomena.

Neglecting the complications from the spinor structure of neutrinos, we may conveniently describe the evolution of the flavor contents of neutrinos by the following kinetic equations (see, e.g., Ref. 191):

where f(t, x, p) is an Nflavor × Nflavor matrix-valued function in phase space (x, p); here, the number of flavors is denoted by Nflavor. The diagonal components of f are the ordinary distribution functions for individual neutrino flavors and the off-diagonal components express the coherence between two flavors specified by the row and column of the matrix. The unit vector n specifies the direction of neutrino momentum. C on the right-hand side is the collision term and is also an Nflavor × Nflavor matrix. H is another matrix of the same size called Hamiltonian and is given as

in which M^2^ is the mass-squared matrix, E is the neutrino energy, and U is Pontecorvo–Maki–Nakagawa–Sakata matrix; v^μ^ = (1, n) is the neutrino four-velocity^#32^ and Λμ=2GF diag[ jαμ ] with jαμ being the lepton number current for the charged lepton of species α; f′ in the last term of Eq. [35] should be understood as a function of v′^μ^; the volume element in momentum space is denoted by dΓ′. The flavor-isospin convention is used, in which the negative energy is reserved for the antineutrino such that f(E, v) = f^T^(−E, v), with T meaning the transposition. The Hamiltonian is essentially the self-energy of neutrinos evaluated at the mean-field level. The third term on the left-hand side of Eq. [34], which is absent in the ordinary kinetic equation, produces the flavor conversions and the three terms in the Hamiltonian given in Eq. [35] are responsible, respectively, for the vacuum oscillation, the MSW effect and neutrino-induced oscillation.

One finds that the third term in the Hamiltonian depends on f. This makes the kinetic equation [34] nonlinear unlike other terms, except the collision term, which will be considered later. This means that the flavor conversion occurs according to the current flavor content, which is modified by the conversion. As a consequence of this nonlinearity, the neutrino-induced neutrino oscillation is much more complex than the other linear oscillations. Its analysis is also complex. It is understandable then to resort to the linearization of Eq. [34] at the flavor eigenstate. This is justified if one focuses on the condition for the occurrence of the neutrino-induced flavor conversion and its exponential growth in the linear regime. In fact, the flavor eigenstate is a fixed point if the vacuum oscillation term is ignored, and this linear analysis may be regarded as a familiar stability analysis of the fixed point, in which we inspect the dispersion relation for the normal modes.^192)^

There are actually two unstable modes, i.e., slow and fast modes. The former grows at the rate proportional to ωvμ , where ω~v~ is the vacuum oscillation frequency ω~v~ = δm^2^/2E ≈ 10^5^ sec^−1^ ≈ 10^−5^ cm^−1^ ≈ 10^−16^ MeV for E ∼ 10 MeV, and μ = GFnν ≈ 10^−8^ MeV ≈ 10^3^ cm^−1^ ≈ 10^13^ sec^−1^ for nν ∼ 10^35^ cm^−3^, which may be appropriate near the neutrinosphere. The growth rate of the latter, by contrast, is given by μ. If μ ≫ ω~v~ as given above, the fast mode grows faster indeed than the slow mode. This is the reason why the fast mode is attracting many researchers at present. The fact that the growth rate of the fast mode does not depend on the vacuum oscillation frequency indicates that the neutrino mass, which was the original driver of the neutrino oscillation in vacuum, plays no essential role in the fast mode. It is indeed a self-induced flavor conversion. It requires only an initial perturbation, which may be given by the mass term.

The linear stability analysis has made clear that the fast mode becomes unstable, i.e., the flavor conversion occurs very quickly when neutrino flavor lepton number (NFLN) defined as

crosses zero at some v for some flavor sector specified by the pair (a, b).^193)^ In the above equations, the spatial dependence is suppressed for notational simplicity. As Fa is the energy-integrated angular distribution function in momentum space for the neutrino species a (= e, μ, τ), the NFLN crossing means that ν~a~ is dominant over ν~b~ in some directions, whereas the opposite is true for the other directions. Such a situation may be rather ubiquitous in the postbounce core^157)^ for various reasons (but see also Ref. 156): for example, the angular distributions in momentum space of all neutrinos are outward-peaked in the preshock region, and ν~e~ is normally dominant in number; hence, Fe is positive in the outward directions; by contrast, as the inward-going neutrinos are mainly produced by backward scatterings of the outward-going neutrinos on heavy nuclei and the scattering rates are larger for ν¯e than for ν~e~ because of the higher average energies for the former, Fe may be negative in the inward direction; |Fμ| and |Fτ| are normally much smaller than |Fe|; then, there should be a crossing at some intermediate direction. In the gain region, on the other hand, the absorption of neutrinos by unbound nucleons may be responsible for the crossing. Complex matter motions are also thought to contribute to the generation of the NFLN crossings deeper inside (see Ref. 157 for more details).

At present, many researchers are attracted by the nonlinear evolution of the fast flavor conversion. In general, numerical simulations are required to study it quantitatively. According to some recent local simulations,^194–196)^ the flavor contents after the nonlinear saturation seem to depend on the angle. Suppose, for example, that Gab(v) is positive in a certain angular domain denoted by Ω~+~ and negative in the remaining domain designated as Ω~−~ and, as a result, there is a crossing on the boundary of the two domains. Let us suppose further that the angular integral of |Gab(v)| over Ω~+~ is larger than its integral over Ω~−. Then, the two flavors are completely mixed up in Ω−, whereas flavor a prevails in Ω+~ so that the angle-integrated NFLN should be conserved (see Fig. 20). The fast flavor conversion occurs pairwise νaν¯a→νbν¯b , and hence, the NFLN integrated over the angle is conserved for each flavor. If their results hold generically indeed, one needs to track the neutrino angular distributions to accurately describe the fast flavor conversion (see also Ref. 197 for more recent developments).

Figure 20. (Color online) NFLN after nonlinear saturation. The vertical axis shows the depolarization factor, which is 0.5 for the complete mixing, whereas it is smaller than 0.5 if the flavor mixing is incomplete. Reproduction of Fig. 2 in S. Bhattacharyya and B. Dasgupta, “Fast flavor depolarization of supernova neutrinos”, Physical Review Letters 126, 061302, (2021), DOI: https://doi.org/10.1103/PhysRevLett.126.061302. (CC-BY)

So far, we have excluded the ordinary collisions of neutrinos, that is, the right-hand side of Eq. [34]. This is because of the fact that the neutrino oscillation is a refractive phenomenon proportional to GF, whereas the collision term is proportional to GF2 and is normally smaller^#33^ (but see also Ref. 198). If one is interested in the neutrino propagation beyond the mean-free path (indeed we are), the effects of the collisions may not be neglected and interesting. In fact, the fast flavor conversion may be enhanced by the collisions^199,200)^ in some cases. There are also papers reporting the suppression of the conversion^201,202)^ and the relevance of this enhancement in the realistic setting remains controversial.^203,204)^ Very recently, it is found that the ordinary collisions may induce flavor conversions on their own.^205)^ This is surprising, as the collisions are normally supposed to be a source of decoherence^198)^ and, as such, are not imagined to be a source of flavor conversions. A detailed linear analysis (not to mention their nonlinear evolution) remains to be conducted to clarify the condition for this collision-induced flavor conversion.^206)^

The eventual goal is to incorporate all relevant flavor conversions in the CCSNe simulations and quantitatively investigate their implications for the neutrino-heating mechanism, or for the CCSNe mechanism more generally. This is going to be a challenge, as the length- and timescales of the flavor conversions may be quite different among different conversion modes and from those of macroscopic hydrodynamics and neutrino transfer. This will be even more so if the saturation of the fast conversion occurs anisotropically in momentum space. Nevertheless, there have been some attempts to perform CCSNe simulations with the flavor conversions considered somehow.^207–209)^ Given this vast difference in the scales, it may be a good idea to construct and implement a subgrid model in the CCSNe simulations^210)^ (see also Ref. 211 for another possibility). Considering that it will be a long way until we can perform the CCSNe simulations with all the relevant flavor conversions fully considered, we had better combine such phenomenological and/or approximate numerical approaches with efforts to understand the properties of individual conversion modes. In fact, new modes may be still waiting for uncovering.^212)^

Last but not least, the mean-field approach, which has been used in most studies so far and on which the derivation of the kinetic equation is based [34], is recently challenged seriously.^212,213)^ In fact, qualitatively different behavior in the flavor conversion is demonstrated for the quantum mechanical approach that solves the basic equations without using the mean-field approximation but for a rather small number of neutrinos. Some counter-arguments to this approach have been also presented.^214,215)^ The controversy reminds us again that we are still in the process of development.

Concluding remarks

We have overviewed the basic scenario, essential ingredients, the current status, and the remaining issues of the neutrino-heating mechanism, the currently most promising one, to drive most of the core-collapse supernovae. Lacking nearby events that allow us to study the inner workings of CCSNe in detail observationally, possibly with neutrinos and GWs in addition to the conventional electromagnetic signals of all wavelengths, we have relied heavily on theoretical investigations to understand what is going on deep inside the massive star. Those who have read through this manuscript will probably agree with us that what the CCSNe modelers have achieved over the last half century is remarkable indeed.

The supernova simulation is still evolving. Very recently, for example, the impact of the numerical treatments of the nuclear reactions in the CCSNe simulation was investigated in detail, based on the 1D and 2D simulations for 20 M⊙ model up to 1.5 sec postbounce.^216)^ The authors implemented not only a usual 16-species α network but also a much larger one with 94 nuclear species included.^#34^ They found that the nuclear burning during the collapse slows down infall and reduces ram pressure. They pointed out that the neutrino-heating rate is also affected by the change in the abundance of nucleons. They claimed that the nuclear energy generation in the postshocked matter amplified the explosion energy by ∼20% in their models. The 3D simulations possibly from 3D progenitors are much awaited.

Recent developments in the chiral effective field theory of nuclear physics^217,218)^ are impressive enough for us to feel rather optimistic about the improvements of the nuclear and weak-interaction physics incorporated in the CCSNe simulations in the near future. The evolution of our understanding of the collective neutrino oscillations is rapid and earnest efforts to include this physics in supernova models are ongoing; however, there is still much remaining to be understood before we can incorporate it in the CCSNe simulation as vindicated by the controversy on the basic equations we have been using.

We have not considered magnetic fields in this article. As mentioned earlier, combined with rapid rotation, they are supposed to play a critical role to produce hypernovae, more than an order more energetic supernovae. They have also roles to play even in the neutrino-heating mechanism for the canonical supernovae. In fact, the magnetic pressure can give an important aid to the thermal pressure and turbulent stress in sustaining the heating region.^45)^ The turbulent magnetic field may serve as an extra energy reservoir and facilitate shock revival.^219)^ We need a better understanding of the strength and configuration of the magnetic field just before the core collapse. The chiral magnetic effects may be another interesting topic in the recent development. They may be induced by a chiral imbalance in electrons and/or charged current interactions of neutrinos out of equilibrium,^220)^ resulting in an additional contribution to the electric current of electrons. The current may amplify the magnetic field via the chiral plasma instability, which may inverse cascade later to help shock revival eventually.^221)^ More quantitative assessments with self-consistent simulations are awaited.

On the other hand, another direction of efforts is going the construction of 1D models, in which an explosion is triggered by artificial methods.^222,223)^ They are meant for systematic studies of various effects parametrically. As such, they are expected to provide complementary information to realistic multidimensional simulations. Although they are calibrated either with some observed quantities or features in the multidimensional simulations, they are admittedly approximate and still need improvements (see, e.g., Ref. 88 for criticisms). Nevertheless, they will be useful to make comparisons with observations more systematically and statistically, and probably more handily.

Combining all these efforts together, we hope that we could obtain a coherent picture of neutrino-driven core-collapse supernovae before the next Galactic event finally occurs.

Acknowledgment

This research used the K and Fugaku supercomputers and the high-performance computing resources of FX100 at Nagoya University ICTS, Grand Chariot at Hokkaido University, and Oakforest-PACS at JCAHPC via the HPCI System Research Project (Project IDs: hp130025, 140211, 150225, 150262, 160071, 160211, 170031, 170230, 170304, 180111, 180179, 180239, 190100, 190160, 200102, 200124, 210050, 210051, 210164, 220047, 2201173, and 220223). The new system Supercomputer “Flow” at Nagoya University ICTS, NEC SX Aurora Tsubasa at KEK, the Research Center for Nuclear Physics (RCNP) at Osaka University, and XC50 and the general common-use computer system operated by CfCA at the National Astronomical Observatory of Japan (NAOJ) also contributed to this study. Large-scale storage of numerical data is provided by Japan Lattice Data Grid (JLDG). This work was supported in part by Grants-in-Aid for Scientific Research (26104006, 15K05093, 19K03837, 19K14723, 20H01905, 21H01083, and 23K03468) and the Grant-in-Aid for Scientific Research on Innovative areas “Gravitational Wave Physics and Astronomy: Genesis” (17H06357 and 17H06365) and “Unraveling the History of the Universe and Matter Evolution with Underground Physics” (19H05802 and 19H05811) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This work was also partly supported by research programs at the K-computer of the RIKEN R-CCS, HPCI Strategic Program of Japanese MEXT, Priority Issue on Post-Kcomputer (Elucidation of the Fundamental Laws and Evolution of the universe), and the Joint Institute for Computational Fundamental Sciences (JICFus). R.A. is supported by the JSPS Grant-in-Aid for JSPS Fellows (Grant No. 22J10298) from MEXT. A.H. is supported by the JSPS Grant-in-Aid for Research Activity Start-up (19K23435) and for Early-Career Scientists (21K13913). S.Y. is supported by the Institute for Advanced Theoretical and Experimental Physics, Waseda University, and the Waseda University Grant for Special Research Projects (project No. 2022C-140).

Profile

Shoichi Yamada, the lead author of this article, was born in Tokyo in 1964. He graduated from the University of Tokyo in 1987 and received his Ph.D. degree from the same university under the supervision of Prof. Katsuhiko Sato in 1992. After a year as a Research Fellow of the Japan Society for the Promotion of Science (JSPS), he became a research associate at the School of Science of the University of Tokyo. He was then awarded a Postdoctoral Fellowship for Research Abroad by the JSPS and visited the Max-Planck-Institute for Astrophysics in Germany and stayed there for 2 years. He extended his stay at the same institute for another 7 months as a guest scientist. After returning to the University of Tokyo and continuing as a research associate for two more years, he moved to the Institute of Laser Engineering, Osaka University, to join the laser astrophysics project led by Prof. Hideaki Takabe. After 2 years, he moved to the Department of Physics, Faculty of Science and Engineering, Waseda University, which is his current affiliation. The supernova explosion, one of the most energetic cosmic events, and the subsequent formation of compact objects such as neutron stars and black holes have been his main research topics. He was attracted to the various physical aspects of this astronomical phenomenon. Starting with hydrodynamics, he studied the emissions of gravitational waves and neutrinos. He was also fascinated with the equation of state of hot and dense hadronic matter as well as with the weak interactions therein. Currently, he is eagerly studying neutrino oscillations induced by their mutual interactions. For supernova research, large-scale numerical simulations are indispensable. The coauthors of this study are all close collaborators of the lead author, who embarked on the first-principles simulation of core-collapse supernova.

Footnotes

References