Abstract

Core-collapse supernovae are dramatic explosions marking the catastrophic end of massive stars. The only means to get direct information about the supernova engine is from observations of neutrinos emitted by the forming neutron star, and through gravitational waves which are produced when the hydrodynamic flow or the neutrino flux is not perfectly spherically symmetric. The multidimensionality of the supernova engine, which breaks the sphericity of the central core such as convection, rotation, magnetic fields, and hydrodynamic instabilities of the supernova shock, is attracting great attention as the most important ingredient to understand the long-veiled explosion mechanism. Based on our recent work, we summarize properties of gravitational waves, neutrinos, and explosive nucleosynthesis obtained in a series of our multidimensional hydrodynamic simulations and discuss how the mystery of the central engines can be unraveled by deciphering these multimessengers produced under the thick veils of massive stars.

1. Introduction

The majority of stars more massive than ~8 end their lives as core-collapse supernovae. They have long attracted the attention of astrophysicists because they have many facets playing important roles in astrophysics. They are the mother of neutron stars as well as black holes; they play an important role for acceleration of cosmic rays; they influence galactic dynamics triggering further star formation; they are gigantic emitters of neutrinos and gravitational waves. They are also a major site for nucleosynthesis, so, naturally, any attempt to address human origins may need to begin with an understanding of core-collapse supernovae (CCSNe).

Current estimates of CCSN rates in our Galaxy predict one event every ~40 ± 10 years [1]. When a CCSN event occurs in our galactic center, copious numbers of neutrinos are produced, some of which may be detected on the earth. Such “supernova neutrinos” will carry valuable information from deep inside the core. In fact, the detection of neutrinos from SN1987A (albeit in the Large Magellanic Cloud) opened up the Neutrino Astronomy, which is an alternative to conventional astronomy by electromagnetic waves [2, 3]. Even though there were just two dozen neutrino events from SN1987A, these events have been studied extensively (yielding ~500 papers) and have allowed us to have a confidence that our basic picture of the supernova physics is correct (e.g., [4], see [5, 6] for a recent review). Recently significant progress has been made in the large water Čherenkov detectors such as Super-Kamiokande [7] and IceCube [8], and also in the liquid scintillator detector as KamLAND [9]. If a supernova occurs in our Galactic center (~10 kpc), about events are estimated to be detectable by IceCube (e.g., [10] and references therein). Those successful neutrino detections are important not only to study the supernova physics but also to unveil the nature of neutrinos itself such as the neutrino oscillation parameters and the mass hierarchy (e.g., [5] for a recent review).

CCSNe are now about to start even another astronomy, Gravitational-Wave Astronomy. Currently long-baseline laser interferometers such as LIGO (USA, e.g., [11]), VIRGO (Italy) (http://www.ego-gw.it/) GEO600 (Germany) (http://geo600.aei.mpg.de/), and TAMA300 (Japan, e.g., [12]) are currently operational and preparing for the first observation (see, e.g., [13] for a recent review), by which the prediction by Einstein’s theory of General Relativity (GR) can be confirmed. These instruments are being updated to their Advanced status and may start taking data, possibly detecting GWs (gravitational waves) for the first time, as soon as 2015 (see [14] for a recent review). In fact, Advanced LIGO/VIRGO, which is an upgrade of the initial LIGO and VIGRO, are expected to be completed by 2015 and will increase the observable detection volume by a factor of ~1000 [15]. The Large-scale Cryogenic Gravitational-wave Telescope (LCGT [16], renamed recently as KAGRA) was funded in late 2010, which is being built under the Kamioka mine and is expected to take its first data in 2016. At such a high level of precision, those GW detectors are sensitive to many different sources, including chirp, ring-down, and merger phases of black-hole and neutron star binaries (e.g., [1719]), neutron star normal mode oscillations (e.g., [20]), rotating neutron star mountains (e.g., [21]), and CCSN explosions (e.g., [2224] for recent reviews), on the final of which we focus in this paper.

According to the Einstein’s theory of GR (e.g., [25]), no GWs can be emitted if gravitational collapse of the supernova core proceeds perfectly spherically symmetric. To produce GWs, the gravitational collapse should proceed aspherically and dynamically. Observational evidence gathered over the last few decades has pointed towards CCSNe indeed being generally aspherical (see [26] for a recent review). The most unequivocal example is SN1987A. To explain the light-curve, a large amount of mixing of Ni outward to the H-He interface and of H inward into the He-core was required (e.g., Woosley [27]; Blinnikov et al. [28]; Utrobin [29], and Woosley & Heger [30] for a review). Such mixing processes coming from the Rayleigh-Taylor instability at the interface with different compositions after shock passage have been examined extensively so far by 2D (e.g., [3133]) and 3D simulations [34]. The asymmetry of iron and nickel lines in SN1987A was proposed to be explained, if the explosion occurs in a jet-like [35, 36] or in a unipolar manner [37]. More directly, the HST images of SN1987A are showing that the expanding envelope is elliptical with the long axis aligned with the rotation axis inferred from the ring ([38], see however [39] for a recent counterargument). The aspect ratio and position angle of the symmetry axis are consistent with those predicted earlier from the observations of speckle and linear polarization. What is more, the linear polarization became greater as time passed (e.g., [4042]), a fact which has been used to argue that the central engine of the explosion is responsible for the nonsphericity (e.g., [33, 4346]). By performing a series of time-dependent, non-LTE (local thermodynamic equilibrium), radiation-transport simulations (e.g., [4749] , see [50, 51] for an alternative approach for the light-curve modeling), Dessart and Hillier [52] recently pointed out that asymmetry in the ejecta can explain the increase in the continuum polarization observed at the nebular phases [53, 54]. Dense knots, indications of ejecta clumpiness, and filaments seen in supernova remnants by HST in the visual [55, 56] and by ROSAT, Chandra, and XMM-Newton [5759] in X-rays also provide evidence that small- and large-scale inhomogeneities (and perhaps even fragmentation) are a common feature in supernova explosions.

Advancing ability of the HST has enabled the direct observation of the progenitors of nearby CCSNe from pre-explosion images (see [60] for a review). Although observational measurements of the progenitor masses are currently not many and still highly uncertain (see [61] for collective references therein), evidence has accumulated that SN type II plateau (II-P) comes predominantly from stars in the range about of 8–16 [60]. A generic explosion energy of the SN II-P in the mass range is roughly on the order of  erg [37], however, a large diversity of the explosion signatures (i.e., explosion energy and the synthesized Ni mass) have been so far observed from quite similar progenitors [60]. For example, the inferred mass and the kinetic energies differ by a factor of five between SNe 2005cs and 2003gd, both of them are among the golden events in which enough information was obtained to give an accurate estimate on a color or spectral type of the progenitor and the initial mass. More massive stars are expected to lose much of their mass and explode as hydrogen-stripped SNe (Ib/c and IIb). Among them, the type Ic-BL SNe, which are associated with long gamma-ray bursts [62], all show much broader lines than SNe Ic. Due to the large kinetic energies of 2–5 erg, they have been referred to as “hypernovae” (e.g., [63] for a recent review). These events are likely to come from interacting binaries in which the primary exploding star has a mass lower than what is usually associated with evolution to the massive WR phase (e.g., [6467], see however [68, 69]). In addition to the two branches mentioned above, Nomoto et al. [70] predicted yet another branch, in which the SNe II-P with higher progenitor mass result in fainter explosions. By connecting to candidate SN explosion mechanisms or to progenitor structures, it is indeed best if one could obtain a unified picture to understand the mentioned wide diversity which is not only dependent on the progenitor masses but also on the evolution scenarios (a single versus binary evolution). But, this is an area for future study firstly because it is too computationally expensive to perform a long-term simulation that follows multidimensional (multi-D) dynamics consistently from the onset of core-collapse, through explosion, up to the nebular phase (for example see [71] for the most up-to-date 1D modeling of spectra and light curves for binary models.), secondly because we are still inaccessible to multi-D stellar evolution models, in which multi-D modeling has been a major undertaking (see [72, 73] for recent developments).

From a theoretical point of view, clarifying what makes the dynamics of the supernova engine deviate from spherical symmetry is essential in understanding the GW emission mechanism. Here it is worth mentioning that GWs are primary observables, which imprint a live information of the central engine, because they carry the information directly to us without being affected in propagating from the stellar center to the earth. On the other hand, SN neutrinos are exposed to a number of (external) environmental effects, including self-interaction that induces collective neutrino flavor oscillations predominantly in the vicinity of the neutrino sphere (see [7477] for reviews of the rapidly growing research field and collective references therein) and the Mikheyev-Smirnov-Wolfenstein (MSW) effect [78] both in the stellar envelope and in the earth (see [5] for a review). Although the time profiles of neutrino signals can be potentially used like a tomography to monitor the envelope profile (e.g., [79] for a recent review), SN neutrinos generally could provide a rather indirect information about the central engine compared to GWs.

The breaking of the sphericity in the supernova engine has been considered as the most important ingredient to understand the explosion mechanism, for which supernova theorists have been continuously keeping their efforts for the past ~40 years. Currently multi-D simulations based on refined numerical models have shown several promising scenarios. Among the candidates are the neutrino heating mechanism aided by convection and hydrodynamic instabilities of the supernova shock (e.g., [80] for a review), the acoustic mechanism [81], or the magnetohydrodynamic (MHD) mechanism (e.g., [22] see references therein). To pin down the true answer among the candidate mechanisms, GW signatures, albeit being a primary observable, will not be solely enough and a careful analysis including neutrinos and photons should be indispensable. The current neutrino detectors is ready to broadcast the alert to astronomers to let them know the arrival of neutrinos (e.g., Supernova Early Warning Systems (SNEWS) [82]). In addition to optical observations using largest 8–10 m telescopes such as VLT and Subaru telescope (e.g., [83, 84]), the planned thirty-meter telescope (TMT) (http://www.tmt.org/) with a refined spectropolarimetric technique is expected to detect more than 20 events of the SN polarization per year (Tanaka in private communication). This should provide valuable information to understand the SN asymmetry with increasing statistics. Not only for understanding the origin of nonspherical ejecta morphology (especially for nearby event) but more importantly for understanding the origin of heavy elements, it is of crucial importance to accurately determine nucleosynthesis in the SN ejecta, which naturally requires a multi-D numerical modeling (see, [30, 85] for as recent reviews, it is worth mentioning that radioactive decay can affect the shape of the light curve for some peculiar SN 1987A-like events [37], in which explosive nucleosynthesis plays an important role for the light-curve modeling).

Putting things together, the multidimensionality determines the explosion mechanism, in turn we may extract the information that traces the multidimensionality by the SN multimessengers, which would be only possible by a careful analysis on GWs, neutrinos, and photons. In this paper, we hope to bring together various of our published and unpublished findings from our recent multi-D supernova simulations and the obtained predictions of the SN multimessengers so far (for other high-energy astrophysical sources such as magnetars, gamma-ray bursts, and coalescing binaries, see [8689] for recent reviews). Before we go into details from the next sections, we first have to draw a caution that the current generation of simulation results that we report in this article should depend on the next generation calculations by which more sophistication can be made not only in determining the efficiency of neutrino-matter coupling (the so-called neutrino transport calculation), but also in the treatment of general relativity. Therefore we provide here only a snapshot of the moving (long-run) documentary film whose headline we (boldly) chose to entitle as “multimessengers from CCSNe to bridge theory and observation.”

Among the mentioned candidate mechanisms, we focus on the neutrino-heating mechanism in Section 2 and the MHD mechanism in Section 3, respectively. In each section, we first briefly summarize the properties of the explosion dynamics and then move on to discuss possible properties of the SN multimessengers paying particular attention to their detectability. It may be best if we can cover these SN messengers once for all in this paper, but unfortunately not. What we have studied so far is limited to GWs and explosive nucleosynthesis in the neutrino-heating mechanism, and to GWs and neutrino signals in the MHD mechanism. To compensate the uncovered fields, the related references will be given. Although a number of excellent reviews already exist on various topics in this paper, this one might go beyond such reviews by its new perspectives on the multimessenger astronomy.

2. Neutrino-Heating Mechanism

CCSN simulations have been counted as one of the most challenging subjects in computational astrophysics. The four fundamental forces of nature are all at play; the collapsing iron core bounces due to strong interactions; weak interactions determine the energy and lepton number loss in the core via the transport of neutrinos; electromagnetic interactions determine the properties of the stellar gas; GR plays an important role due to the compactness of the protoneutron star and also due to high velocities of the collapsing material outside. Naturally, such physical richness ranging from a microphysical scale (i.e., femto-meter scale) of strong/weak interactions to a macrophysical scale of stellar explosions has long attracted the interest of researchers, necessitating a worldwide multidisciplinary collaboration to clarify the theory of massive stellar core-collapse and the formation mechanisms of compact objects.

Ever since the first numerical simulation of such events [90], the neutrino-heating mechanism [9193], in which a stalled bounce shock could be revived via neutrino absorption on a timescale of several hundred milliseconds after bounce, has been the working hypothesis of supernova theorists for these ~45 years. However, the simplest spherically-symmetric (1D) form of this mechanism fails to blow up canonical massive stars [9497]. Pushed by mounting observations of the blast morphology (e.g., [26]) mentioned above, it is now almost certain that the breaking of the spherical symmetry is the key to solve the supernova problem. So far a number of multi-D hydrodynamic simulations have been reported, which demonstrated that hydrodynamic motions associated with convective overturn (e.g., [98102]) as well as the Standing-Accretion-Shock-Instability (SASI, e.g., [44, 103113] see references therein) can help the onset of the neutrino-driven explosion.

To test the neutrino-heating mechanism in the multi-D context, it is of crucial importance to solve accurately the neutrino-matter coupling in spatially nonuniform hydrodynamic environments. For the purpose, one ultimately needs to solve the six-dimensional (6D) neutrino radiation transport problem (three in space, three in the momentum space of neutrinos), which is a main reason why the supernova simulations stand out from other astrophysical simulations due to their complexity. In the final sentence of the last paragraph, we wrote “demonstrated” because the neutrino heating was given by hand as an input parameter in most of the simulations cited above (see, however [101, 102]). The neutrino heating proceeds dominantly via the charged current interactions () in the gain region. The neutrino heating rate in the gain region can be roughly expressed as , where is the neutrino luminosity emitted from the surface of neutrino sphere and it determines the amplitude of the neutrino heating as well as cooling, and and are the distance from the stellar center and the flux factor (this quantity represents the degree of anisotropy in neutrino emission; near at the neutrino sphere, in the free-streaming limit ()), respectively (e.g., [114]). For example, is treated as an input parameter in the so-called “light-bulb” approach (e.g., [100]). This is one of the most prevailing approximations in recent 3D simulations [111, 112, 115, 116] because it is handy to study multi-D effects on the neutrino heating mechanism (albeit on the qualitative grounds). To go beyond the light-bulb scheme, should be determined in a self-consistent manner. For the purpose, one needs to tackle with neutrino transport problem, only by which energy as well as angle dependence of the neutrino distribution function can be determined without any assumptions. Since the focus of this review is on the SN multimessengers, a detailed discussion of various approximations and numerical techniques taken in the recent radiation-hydrodynamic SN simulations cannot be provided. Table 1 is not intended as a comprehensive compilation, but we just want to summarize milestones that have recently reported the neutrino-driven exploding models so far.

In Table 1, the first column (“Progenitor”) shows the progenitor model employed in each simulation. The abbreviations of “NH,” “WHW,” and “WW” mean [117119]. The second column shows SN groups with the published or submitted year of the corresponding work. “MPA” stands for the CCSN group in the Max Planck Institute for Astrophysics led by Janka and Müller. “Princeton+” stands for the group chiefly consisting of the staffs in the Princeton University (Burrows, Murphy), Caltech (Ott), Hebrew University (Livne), Université de Provence (Dessart), and their collaborators. The SN group in the Basel university is led by Liebendörfer and Thielemann. “OakRidge+” stands for the SN group mainly consisting of the Florida Atlantic University (Bruenn) and the Oak Ridge National Laboratory (Mezzacappa, Messer) and their collaborators. Tokyo+ is the SN group chiefly consisting of the staffs in the National Astronomical Observatory of Japan (myself, Takiwaki), Kyoto university (Suwa), Numazu College (Sumiyoshi), Waseda University (Yamada), and their collaborators. The third column represents the mechanism of explosions which are basically categorized into two (to date), namely by the neutrino-heating mechanism (indicated by “-driven”) or by the acoustic mechanism (“Acoustic”). “Dim.” in the fourth column is the fluid space dimensions which is one-, two-, or three-dimension (1-,2-,3D). The abbreviation “N” stands for “Newtonian,” while “PN”—for “Post-Newtonian”—stands for some attempt at inclusion of general relativistic effects, and “GR” denotes full relativity. in the fifth column indicates an approximate typical timescale when the explosion initiates and represents the explosion energy normalized by Bethe (= erg) given at the postbounce time of , both of which are attempted to be sought in literatures (but if we cannot find them, we remain them as blank “—”). In the final column of “ transport,” “Dim” represents momentum dimensions and the treatment of the velocity-dependent term in the transport equations is symbolized by . The definition of the “RBR,” “IDSA,” and “MGFLD” will be given soon in the following.

Due to the page limit of this paper, we have to start the story only after 2006 (see, e.g., [80] for a complete review, and also [120] for a similar table before 2006). The first news of the exploding model was reported by the MPA group. By performing radiation-neutrino-hydrodynamic simulations which includes one of the best available neutrino transfer approximations, they reported 1D and 2D explosions for the 8.8 star [121, 122] whose progenitor has a very tenuous outer envelope with steep density gradient (a characteristic property of AGB stars). Also in 1D, the Basel+ group reported explosions for 10 and 15 progenitors of WHW02 triggered by the hypothesised first-order QCD phase transition in the protoneutron star (PNS) [123]. To date, these two are the only modern numerical results where the neutrino-driven mechanism succeeded in 1D. In the 2D MPA simulations, they obtained explosions for a nonrotating 11.2 progenitor of WHW02 [124], and then for a progenitor [125] of WW95 with a relatively rapid rotation imposed (by comparing the precollapse angular velocity to the one predicted in a recent stellar evolution calculation [126, 127]). They newly brought in the so-called “ray-by-ray” approach (indicated by “RBR” in the table), in which the neutrino transport is solved along a given radial direction assuming that the hydrodynamic medium for the direction is spherically symmetric. This method, which reduces the 2D problem partly to 1D, fits well with their original 1D Boltzmann solver [94] (note that the ray-by-ray approach has an advantage compared to other approximation schemes, such that it can fully take into account the available neutrino reactions (e.g., [128] for references therein) and also give us the most accurate solution for a given angular direction). For 2D hydrodynamic simulations with the ray-by-ray transport, one needs to solve the 4D radiation transport problem (two in space and two in the neutrino momentum space). Regarding the explosion energies obtained in the MPA simulations, their values at their final simulation time are typically underpowered by one or two orders of magnitudes to explain the canonical supernova kinetic energy (~1051 erg). But the explosion energies presented in their figures are still growing with time, and they could be as high as 1 B if they were able to follow a much longer evolution as discussed in [124]. Very recently, Mueller et al. [129] reported that more energetic explosions are obtained for the 11.2 and stars in their GR 2D simulations compared to the corresponding post-Newtonian models (e.g., [125]). In the table, the data in italic character represent their most up-to-date results based on the 2D GR simulations using conformally-flatness approximation (e.g., [130, 131], indicated by “C-GR” in the table).

It is rather only recently that fully 2D multiangle Boltzmann transport simulations become practicable by the Princeton+ group [132, 133]. In this case, one needs to handle the 5D problem for 2D simulations (two in space, and three in the neutrino momentum space). However this scheme is very computationally expensive currently to perform long-term supernova simulations. In fact, the most recent 2D work by Brandt et al. [133] succeeded in following the dynamics until ~400 ms after bounce for a nonrotating and a rapidly rotating 20 model of WHW02, but explosions seemingly have not been obtained in such an earlier phase either by the neutrino-heating or the acoustic mechanism.

In the table, “MGFLD” stands for the multigroup Flux-limited diffusion scheme which eliminates the angular dependence of the neutrino distribution function (see, e.g., [134] for more details). For 2D simulations, one needs to solve the 3D problem, namely two in space, and one in the neutrino momentum space. By implementing the MGFLD algorithm to the CHIMERA code in a ray-by-ray fashion (e.g., [135]), Bruenn et al. (2009) obtained neutrino-driven explosions for nonrotating progenitors in a relatively wide range in 12, 15, 20, 25 of WHW02 (see Table 1).These models tend to start exploding at around 300 m after bounce, and the explosion energy for the longest running model of the 25 progenitor is reaching to 1 B at 1.2 s after bounce [135].

On the other hand, the 2D MGFLD simulations implemented in the VULCAN code [136] obtained explosions for a variety of progenitors of 11, 11.2, 13, 15, 20, and 25 not by the neutrino-heating mechanism but by the acoustic mechanism (for the 8.8 progenitor, they obtained neutrino-driven explosions (see Table 1)). The acoustic mechanism relies on the revival of the stalled bounce shock by the energy deposition via the acoustic waves that the oscillating protoneutron stars (PNSs) would emit in a much delayed phase (~1 second) compared to the conventional neutrino-heating mechanism (~300–600 milliseconds). If the core pulsation energy given in Burrows et al. [45] could be used to measure the explosion energy in the acoustic mechanism, they reach to 1 B after 1000 ms after bounce.

By performing 2D simulations in which the spectral neutrino transport was solved by the isotropic diffusion source approximation (IDSA) scheme [137], the Tokyo+ group reported explosions for a nonrotating and rapidly rotating 13 progenitor of NH88. They pointed out that a stronger explosion is obtained for the rotating model comparing to the corresponding nonrotating model. The IDSA scheme splits the neutrino distribution into two components (namely the streaming and trapped neutrinos), both of which are solved using separate numerical techniques (see [137] for more details). The approximation level of the IDSA scheme is basically the same as the one of the MGFLD. The main advantage of the IDSA scheme is that the fluxes in the transparent region can be determined by the nonlocal distribution of sources rather than the gradient of the local intensity like in MGFLD. A drawback in the current version of the IDSA scheme is that heavy lepton neutrinos (, i.e., , and their antiparticles) as well as the energy-coupling weak interactions have yet to be implemented. Extending the 2D modules in [138] to 3D, they recently reported explosions in the 3D models for an 11.2 progenitor of WHW02 [139]. By comparing the convective motions as well as neutrino luminosities and energies between their 2D and 3D models, they pointed out whether 3D effects would help explosions or not is sensitive to the employed numerical resolutions. They argued that next-generation supercomputers are at least needed to draw a robust conclusion of the 3D effects.

Having summarized a status of the current supernova simulations, one might easily see a number of issues that remain to be clarified. First of all, the employed progenitors usually rather scatter (e.g., Table 1). Different SN groups seem to have a tendency to employ different progenitors, providing different results. By climbing over a wall which may have rather separated exchanges among the groups, a detailed comparison for a given progenitor needs to be done seriously in the multi-D results (as have been conducted in the Boltzmann 1D simulations between the MPA, Basel+, and Oak Ridge+ groups [140]).

In the next section, we briefly summarize the findings obtained in our 2D [138] and 3D [139] simulations, paying particular attention to how multidimensionality such as SASI, convection, and rotation could affect the neutrino-driven explosions.

2.1. Multidimensionality in Multi-D Radiation Hydrodynamic Simulations

Figure 1 depicts the difference between the time evolutions of 1D (thin gray lines) or 2D (thin orange lines) simulation of the 13 progenitor model. Until ~100 ms after bounce, the shock position of the 2D model (thick red line) is similar to the 1D model (thick black line). Later on, however, the shock for 2D does not recede as for 1D, but gradually expands and reaches 1000 km at about 470 ms after bounce. Comparing the position of the gain radius (red-dashed line) to the shock position for 1D (thick black line) and 2D (thick red line), one can see that the advection time of the accreting material in the gain region can be longer in 2D than in 1D. This longer exposure of cool matter in the heating region to the irradiation of hot outstreaming neutrinos from the PNS is essential for the increased efficiency of the neutrino heating in multi-D models. This can be also depicted in Figure 2. The right half shows , which is the ratio of the advection to the neutrino heating timescale. This quantity is known as a useful quantity to diagnose the success (, i.e., the neutrino-heating timescale is shorter than the advection timescale of material in the gain region) or failure () of the neutrino-driven explosion (e.g., [141]). For the 2D model (right panel), it can be shown that the condition of is satisfied behind the aspherical shock, which is deformed predominantly by the SASI ( mode at the snapshot), while the ratio is shown to be smaller than unity in the whole region behind the spherical standing accretion shock (left panel: 1D). Here SASI that develops due to the advection-acoustic cycle in the supernova core [142, 143], is a uni- and bipolar sloshing of the stalled bounce shock with pulsational strong expansion and contraction (seen as oscillations in the red curves in Figure 1). Comparing the left half of each panel, the entropy for the 2D model is shown to be larger than for the 1D model. This is also the evidence that the neutrino heating works more efficiently in 2D.

To see clearly the effects of stellar rotation on the postbounce evolution, a rapidly rotating model was taken in [138], in which a constant angular frequency of /s is imposed inside the iron core with a dipolar cut off (r−2) outside. This corresponds to with being the ratio of the rotational to the gravitational energy when the simulation starts. For the rotating model, the dominant mode of the shock deformation after bounce is almost always the mode for the rotating model (e.g., Figure 3, although the mode can be as large as the mode when the SASI enters the nonlinear regime (200 ms after bounce). In contrast to this rotation-induced deformation, the mode tends to be larger than the mode for the 2D models without rotation in the saturation phase.

The top panel of Figure 4 shows the comparison of the explosion energies for the 2D models with and without rotation (models with “-rot” or “-2D,” resp.). Although the explosion energies depend on the numerical resolutions quantitatively, they show a continuous increase for the rotating models. The explosion energies for the models without rotation, on the other hand, peak around ms when the neutrino-driven explosion sets in (see also Figure 1), and show a decrease later on. The reason for the greater explosion energy for models with rotation is due to the bigger mass of the exploding material. This is because the northsouth symmetric () explosion can expel more material than for the unipolar explosion (compare the left and right of the bottom panel in Figure 4). The explosion energies when we terminated the simulation are less than 1050 erg for all the computed models. For the rotating models, we are tempted to speculate that the explosion energies could increase later on by a linear extrapolation. However, in order to identify the robust feature of an explosion, a longer-term simulation with improved input physics is needed. In combination with the assumed rapid rotation, magnetic fields should be also taken into account as in [46], which we are going to study as a followup of our rotating 2D models (Suwa et al. in preparation).

Extending our 2D modules mentioned above, we are currently running 3D simulations for an 11.2 star of Woosley et al. [118] with spectral neutrino transport that is solved by the IDSA scheme in a ray-by-ray manner [139]. We briefly summarize the results in the following.

The left panel of Figure 5 shows mass-shell trajectories for the 3D (red lines) and 1D model (green line), respectively. At around 300 ms after bounce, the average shock radius for the 3D model exceeds 1000 km in radius. On the other hand, an explosion is not obtained for the 1D model, which is in agreement with Buras et al. [124]. The right panel of Figure 5 shows a comparison of the average shock radius versus postbounce time. In the 2D model, the shock expands rather continuously after bounce. This trend is qualitatively consistent with the 2D result by Buras et al. [124] (see their Figure 15 for model s112_128_f). The reason that the average shock of our 2D model expands much faster than theirs would come from the neglected effects in this work including general relativistic effects, inelastic neutrino-electron scattering, and cooling by heavy-lepton neutrinos. All of them could give a more optimistic condition to produce explosions. Apparently these ingredients should be appropriately implemented, which we hope to be practicable in the next-generation 3D simulations.

Comparing the shock evolution between our 2D (green line in the right panel of Figure 5) and 3D models (red line), the shock is shown to expand much faster for 2D. The pink line labeled by “3D low” is for the low resolution 3D model, in which the mesh numbers are taken to be half of the standard model. Note that the 3D computational grid consists of 300 logarithmically spaced, radial zones to cover from the center up to 5000 km and 64 polar () and 32 azimuthal () uniform mesh points, which are used to cover the whole solid angle. The low resolution 3D model has one-half of the mesh numbers in the direction ( = 16), while fixing the mesh numbers in other directions. Comparing with our standard 3D model (red line), the shock expansion becomes less energetic for the low resolution model (later than ~150 ms). Above results indicate that explosions are easiest to obtain in 2D, followed in order by 3D, and 3D (low). At first sight, this may look contradicted with the finding by Nordhaus et al. [116] who pointed out that explosions could be more easily obtained in 3D than in 2D. The reason of the discrepancy is summarized shortly as it follows.

Figure 6 compares the blast morphology for our 3D (left panel) and 2D (right) model (note that the polar axis is tilted (about ) both in the left and middle panel). In the 3D model (left panel), nonaxisymmetric structures are clearly seen. By performing a tracer-particle analysis, the maximum residency time of material in the gain region is shown to be longer for 3D than 2D due to the nonaxisymmetric flow motions (see Figure 7). This is one of advantageous aspects of 3D models to obtain the neutrino-driven explosions. On the other hand, our detailed analysis showed that convective matter motions below the gain radius become much more violent in 3D than in 2D, making the neutrino luminosity larger for 3D (see [139] for more details). Nevertheless the emitted neutrino energies are made smaller due to the enhanced cooling. Due to these competing ingredients, the neutrino-heating timescale becomes shorter for the 3D model, leading to a smaller net-heating rate compared to the corresponding 2D model (Figure 8). Note here that the spectral IDSA scheme, by which the feedback between the mass accretion and the neutrino luminosity can be treated in a self-consistent manner (not like the light-bulb scheme assuming a constant luminosity), sounds quite efficient in the first-generation 3D simulations.

As seen from Figure 5, an encouraging finding in [139] was that the shock expansion tends to become more energetic for models with finer resolutions. These results would indicate whether these advantages for driving 3D explosions can or cannot overwhelm the disadvantages is sensitive to the employed numerical resolutions (It is of crucial importance to conduct a convergence test in which a numerical gridding is changed in a systematic way (e.g., [144]).) To draw a robust conclusion, 3D simulations with much more higher numerical resolutions and also with more advanced treatment of neutrino transport as well as of gravity are needed, which could be hopefully practicable by utilizing forthcoming petaflops-class supercomputers.

In addition to the 3D effects, impacts of GR on the neutrino-driven mechanism stand out among the biggest open questions in the supernova theory. It should be remembered that using newly derived Einstein equations [146], the consideration of GR was standard in the pioneering era of supernova simulations (e.g., [147]). One year after Colgate and White [90], Schwartz [148] reported the first fully GR simulation of stellar collapse to study the supernova mechanism, who implemented a gray transport of neutrino diffusion in the 1D GR hydrodynamics (citing from his paper, “In this calculation, the neutrino luminosity of the core is found to be 1054 erg/s, or 1/2 a solar rest mass per second !! This is the mechanism which the supernova explodes.” The neutrino luminosity rarely becomes so high in the modern simulations, but it is surprising that the potential impact of GR on the neutrino-heating mechanism was already indicated in the very first GR simulation). Using GR Boltzmann equations derived by Lindquist [149], Wilson [150] developed a 1D GR-radiation-hydrodynamic code including a more realistic (at the time) description of the collisional term than the one in [148]. By performing 1D GR hydrodynamic simulations that included a leakage scheme for neutrino cooling, hydrodynamical properties up to the prompt shock stagnation were studied in detail [151153]. These pioneering studies, albeit using a much simplified neutrino physics than today, did provide a bottom-line of our current understanding of the supernova mechanism (see [154] for a complete list of references for the early GR studies). In the middle of the 1980s, Bruenn [134] developed a code that coupled 1D GR hydrodynamics to the MGFLD transport of order including the so-called standard set of neutrino interactions. Since the late 1990s, the ultimate 1D simulations, in which the GR Boltzmann transport is coupled to 1D GR hydrodynamics, have been made feasible by Yamada et al. [155158] (very recently, they reported their success to develop the first multiangle, multienergy neutrino transport code in 3D [159]) and by Mezzacappa et al. [154, 160162] (and by their collaborators).

Among them, Bruenn et al. [154] firstly showed that the neutrino luminosity and the average neutrino energy of any neutrino flavor during the shock reheating phase increase when switching from Newtonian to GR hydrodynamics. They also pointed out that the increase is larger in magnitude compared to the decrease due to redshift effects and gravitational time dilation. By employing the current best available weak interactions, Lentz et al. [163] reported the update of Bruenn et al. [154] very recently. They showed that the omission of observer corrections in the transport equation particularly does harm to drive the neutrino-driven explosions. A disadvantageous trend of GR to produce the neutrino-driven explosions has been commonly observed in these full-fledged 1D simulations; the residency time of material in the gain region becomes shorter due to the stronger gravitational pull. Due to this competition between the gain and loss effects in the end, GR works disadvantageously to facilitate the neutrino-driven explosions in 1D. In fact, the maximum shock extent in the postbounce phase is shown to be 20% smaller when switching from Newtonian to GR hydrodynamics (e.g., Figure 2 in [163]).

Among the most up-to-date multi-D models with spectral neutrino transport mentioned earlier (Table 1), the GR effects are best attempted to be modeled by using a modified gravitational potential that takes into account a 1D, post-Newtonian correction [124, 125, 135, 164]. A possible drawback of this prescription is that a conservation law for the total energy cannot be guaranteed by adding an artificial term in the Poisson equation. Since the energy reservoir of the supernova engines is the gravitational binding energy, any potential inaccuracies in the argument of gravity should be avoided. There are a number of relativistic simulations of massive stellar collapse in full GR (e.g., 2D [165] or 3D [166, 167], and references therein) or using the conformally-flatness approximation (CFC) (e.g., [130, 131]). Although extensive attempts have been made to include microphysics such as by the formula [168] or by the neutrino leakage scheme [169], the effects of neutrino heating have yet to be included in them, which has been a main hindrance to study the GR effects on the multi-D neutrino-driven mechanism (see, however, [129, 170]).

Putting things together, a complete “realistic” supernova model should naturally be done in full GR (MHD) with multi-D GR Boltzmann neutrino transport, in which a microphysical treatment of equation of state (EOS) and nuclear-neutrino interactions are implemented as realistically as possible. Unfortunately none of the currently published SN simulations satisfy the “ultimate” requirement (e.g., Table 1). In this sense, all the mentioned studies employ some approximations (with different levels of sophistication) towards the final goal.

In the same way, theoretical predictions of the SN multimessengers that one can obtain by analyzing the currently available numerical results, cannot unambiguously give us the final answer yet. Again we hereby note the feature of this paper which shows only a snapshot of the moving theoretical terrain. Keeping this caveat in mind, it is also true that a number of surprising GW features of CCSNe have been reported recently both by the first-principle simulations (e.g., in Table 1) and also by idealized simulations in which explosions are parametrically initiated mostly by the light-bulb scheme. As will be mentioned in the next section, the latter approach is also useful to get a better physical understanding of the GW signatures obtained in the first-principle simulations (in this sense, these two approaches are complimentary in understanding the GW signatures). Having shortly summarized a current status of CCSN simulations, we are now ready to move on to focus on the GW signatures from the next section.

2.2. Gravitational Waves

The paper by Müller [171] entitled as “Gravitational Radiation from Collapsing Rotating Stellar Cores” unquestionably opened our eyes to the importance of making the GW prediction based on realistic SN numerical modeling (see e.g., Section 2 in [23] for a summary of more earlier work which had mainly focused on the GW emission in very idealized systems such as in homogeneous spheroids and ellipsoids) (needless to say, this kind of approach is still very important to extract the physics of the GW emission mechanism). As one may expect from the title of his paper, rapid rotation, if it would exist in the precollapse iron core, leads to significant rotational flattening of the collapsing and bouncing core, which produces a time-dependent quadrupole (or higher) GW emission. Following the first study by Müller, most studies of the past thirty have focused on the so-called bounce signals (e.g., [171185]), and references therein).

As summarized by Ott [23], a number of important progresses have been recently made to understand features of the bounce signals by extensive 2D GR studies using the CFC approximation [181183] and also by fully GR 3D simulations [178] both including realistic EOSs and a deleptonization effect based on 1D-Boltzmann simulations [186].

For the bounce signals having a strong and characteristic signature, the iron core must rotate enough rapidly. However recent stellar evolution calculations suggest that rapid rotation assumed in most of the previous studies is not canonical for progenitors with neutron star formations [126, 127, 187]. To explain the observed rotation periods of radio pulsars, the precollapse rotation periods are estimated to be larger than ~100 sec [187]. In such a slowly rotating case, the detection of the bounce signals becomes very hard even by next-generation laser interferometers for a Galactic supernova (e.g., [176]).

Besides the rapid rotation, convective matter motions and anisotropic neutrino emission in the postbounce phase are expected to be primary GW sources with comparable amplitudes to the bounce signals (e.g., [188] for a review). Thus far, various physical ingredients for producing asphericities and the resulting GWs in the postbounce phase have been studied, such as the roles of precollapse density inhomogeneities [102, 189, 190], moderate rotation of the iron core [191], nonaxisymmetric rotational instabilities [192, 193], g-modes [194] and r-modes pulsations [20] of PNSs, and the SASI [145, 195198] (see also [199] for the GW signals at the black hole formation). Among them, the most promising GW sources may be convection and SASI, because the degree of the initial inhomogeneities [200] and the growth of rotational instabilities as well as the r-modes and g-modes pulsations are rather uncertain.

Based on 2D simulations that parametrize the neutrino heating and cooling by the light-bulb scheme (which we shortly call as parametric SASI simulations hereafter), we pointed out that the GW amplitudes from anisotropic neutrino emission (as anisotropic matter motions generate GWs, anisotropic neutrino emission also gives rise to GWs, which has been originally pointed out in late 1970’s by Epstein [201] and Turner and Wagoner [202] (see recent progress in [203]). It is expected as a primary GW source also in gamma-ray bursts [204206] and Pop III stars [207]) increase almost monotonically with time, and that such signals may be visible to next-generation detectors for a Galactic source [195, 196]. By performing such a parametric simulation but without the excision inside the PNS, Murphy et al. [198] showed that the GW signals from matter motions can be a good indicator of the explosion geometry. These features qualitatively agree with the ones obtained by Yakunin et al. [208] who reported exploding 2D simulations in which the ray-by-ray MGFLD neutrino transport is solved with the hydrodynamics (e.g., Oak-Ridge+ simulations in Table 1). Marek et al. [197] analyzed the GW emission based on their long-term 2D ray-by-ray Boltzmann simulations, which seem very close to produce explosions [125] (e.g., MPA simulations in Table 1). They also confirmed that the GWs from neutrinos with continuously growing amplitudes (but with the different sign of the amplitudes in [195, 196, 208]), are dominant over the ones from matter motions. They proposed that the third-generation class detectors such as the Einstein Telescope are required for detecting the GW signals with a good signal-to-noise ratio.

Regarding the GW predictions in 3D models, Müller and Janka [190] coined the first study to analyze the GW signature of 3D nonradial matter motion and anisotropic neutrino emission from prompt convection in the outer layers of a PNS during the first 30 ms after bounce. Their first 3D calculations using the light-bulb recipe were forced to be performed in a wedge of opening angle of . Albeit with this limitation (probably coming from the computer power at that time), they obtained important findings that because of smaller convective activities inside the cone with slower overturn velocities, the GW amplitudes of their 3D models are more than a factor of 10 smaller than those of the corresponding 2D models, and the wave amplitudes from neutrinos are a factor of 10 larger than those due to nonradial matter motions. With another pioneering (2D) study by Burrows and Hayes [189] it is worth mentioning that those early studies had brought new blood into the conventional GW predictions, which illuminated the importance of the theoretical prediction of the neutrino GWs.

A series of findings obtained by Fryer et al. in early 2000s [101, 209] have illuminated also the importance of the 3D modeling. By running their 3D Newtonian Smoothed-Particle-Hydrodynamic (SPH) code coupled to a gray flux-limited neutrino transport scheme, they studied the GW emission due to the inhomogeneous core-collapse, core rotation, low-modes convection, and anisotropic neutrino emission. Although the early shock-revival and the subsequent powerful explosions obtained in these SPH simulations have yet to be confirmed by other groups, their approach paying particular attention to the multiple interplay between the explosion dynamics, the GW signatures, the kick and spins of pulsars, and also the nonspherical explosive nucleosynthesis, blazed a new path on which current supernova studies are progressing.

We also studied the GW signals from 3D models that mimic neutrino-driven explosions aided by the SASI [145, 196, 210]. In the series of our 3D experimental simulations, the light-bulb scheme was used to obtain explosions and the initial conditions were derived from a steady-state approximation of the postshock structure and the dynamics only outside an inner boundary at 50 km was solved. Based on the results, we show in the following that features of the gravitational waveforms obtained in 3D models are significantly different than those in 2D, which tells us a necessity of 3D modeling for a reliable prediction of the GW signals.

2.2.1. Stochastic Nature of Gravitational Waves in 3D Simulations

Figure 9 shows the evolution of 3D hydrodynamic features from the onset of the nonlinear regime of SASI (top left) until the shock breakout (this corresponds to the shock emergence at the outer boundary of the computational domain (~2000 km in radius)) (bottom right) with the gravitational waveform from neutrinos inserted in each panel. After about ms, the deformation of the standing shock becomes remarkable marking the epoch when the SASI enters the nonlinear regime (top left of Figure 9). At the same time, the gravitational amplitudes begin to deviate from zero. Comparing the top two panels in Figure 10, which shows the total amplitudes (top panel, neutrino + matter) and the neutrino contribution only (bottom), it can be seen that the overall structures of the waveforms are predominantly determined by the neutrino-originated GWs with the slower temporal variations (30–50 ms), to which the GWs from matter motions with rapid temporal variations (10 ms) are superimposed.

As seen from the top right through bottom left to right panels of Figure 9, the major axis of the growth of SASI is shown to be not aligned with the symmetric axis (z-axis in the figure) and the flow inside the standing shock wave is not symmetric with respect to this major axis (see the first and third quadrants in Figure 9). This is a generic feature in the computed 3D models, which is in contrast to the axisymmetric case. The GW amplitudes from SASI in 2D showed an increasing trend with time due to the symmetry axis, along which SASI can develop preferentially [195, 196]. Free from such a restriction, a variety of the waveforms is shown to appear (see waveforms inserted in Figure 5). Furthermore, the 3D standing shock can also oscillate in all directions, which leads to the smaller explosion anisotropy than 2D. With these two factors, the maximum amplitudes seen either from the equator or the pole become smaller than 2D. On the other hand, their sum in terms of the total radiated energy are found to be almost comparable between 2D and 3D models, which is likely to imply the energy equipartition with respect to the spatial dimensions.

The (two) pair panels of Figure 10 show the gravitational waveforms for models with different neutrino luminosity. The input luminosity for the pair panels (models A (top two) and B (bottom two)) differs only (although any seed perturbation would demonstrate the chaotic behavior, we simply chose the value of 0.5% as a reference). Despite the slight difference in the luminosities, the waveforms of each polarization are shown to exhibit no systematic similarity when seen from the pole or equator. This also reflects a chaotic nature of the SASI influenced by small differences (see also [211]).

It should be noted that the approximations taken in the simulation, such as the excision inside the PNS with its fixed inner boundary and the light bulb approach with the isotropic luminosity constant with time, are the very first step to model the dynamics of the neutrino-heating explosion aided by SASI and study the resulting GWs. As already mentioned, the excision of the central regions inside PNSs may hinder the efficient gravitational emission of the oscillating neutron star [194] and the nonaxisymmetric instabilities [185, 193] of the PNSs, and the enhanced neutrino emissions inside the PNSs [197]. Bearing these caveats in mind, a piece of encouraging news is that the gravitational waveforms obtained in the 2D radiation-hydrodynamic simulations [208] are similar to the ones obtained in our 2D study using the light-bulb scheme [196]. More recently, Müller et al. [211] confirmed the stochastic nature by analyzing the GW features obtained in their 3D models [115]. The obtained GW amplitudes as well as the degree of anisotropic neutrino emission are almost comparable to our results, which are considerably smaller than 2D models [191, 197, 198, 208].

2.2.2. Breaking of the Stochasticity due to Stellar Rotation

More recently, we studied the effects of stellar rotation on the stochastic nature of the GWs mentioned above [210]. In 3D, the modes of SASI are divided into sloshing modes and spiral modes (e.g., [111]). Asymmetric modes so far studied in 2D models and axisymmetric modes are classified into the sloshing modes, where stands for the azimuthal index of the spherical harmonics . In the latter situation, the modes degenerate so that the modes has the same amplitudes as the modes. If random perturbations or uniformly rotating flows are imposed on these axisymmetric flows in the postbounce phase, the degeneracy is broken and the rotational modes emerge [112, 212]. In this situation, the modes has the different amplitudes from modes. Such rotating nonaxisymmetric modes are called as spiral modes. These nonaxisymmetric modes are expected to bring about a breakthrough in our supernova theory, because they can help to produce explosions more easily compared to 2D due to its extra degree of freedom [116], and also because they may have a potential to generate pulsar spins ([212, 213], see also [115, 214]). These findings illuminating the importance of stellar rotation motivated us to clarify how rotation that gives a special direction to the system (i.e., the spin axis), could affect the stochastic GW features which we observed in the absence of rotation (Section 2.2.1).

Figure 11 shows the gravitational waveform for a typical 3D model with rotation (left: total amplitudes, right: neutrino only). To construct a model with rotation, we give a uniform rotation on the flow advecting from the outer boundary of the iron core as in [112], whose specific angular momentum is assumed to agree with recent stellar evolution models [126, 127]. Comparing to Figure 10 (for models without rotation), one can clearly see a sudden rise in the GW amplitude after around 500 ms for the rotating model (blue line in Figure 11), which is the plus mode of the neutrino GW seen from the equator. By systematically changing the initial angular momentum and the input neutrino luminosities from the PNS, we computed fifteen 3D models in [210] and found that the GW features mentioned above were common among the models.

Figure 12 illustrates a typical snapshot of the flow fields for the rotating model (corresponding to the one in Figure 11) when the spiral SASI modes have already entered the nonlinear regime, seen from the pole (right panel) or from the equator (left panel), respectively. From the left panel, one may guess the presence of the sloshing modes that happen to develop along the rotational axis (-axis) at this epoch. It should be emphasized that although the dominance of observed in the current 3D simulations is similar to the one obtained in previous 2D studies [196], its origin has nothing to do with the coordinate symmetry axis. The preferred direction here is determined by the spin axis. Free from the 2D axis effects, the major axis of the sloshing SASI mode changes stochastically with time, and the flow patterns behind the standing shock simultaneously change in every direction like the nonrotating models. As a result, the sloshing modes can make only a small contribution to the GW emission. The remaining possibility is that the spiral flows seen in the right panel should be a key importance to understand the GW feature mentioned above. In fact, by analyzing the matter distribution on the equatorial plane, we find that the compression of matter is more enhanced in the vicinity of the equatorial plane due to the growth of the spiral SASI modes, leading to the formation of the spiral flows circulating around the spin axis with higher temperatures. As a result, the neutrino emission seen parallel to the spin axis becomes higher than the ones seen from the other direction. Remembering that the lateral-angle () dependent function of the GW formulae (e.g., in in [196]) is positive near the north and south polar caps, the dominance of the polar neutrino luminosities leads to make the positively growing feature of in Figure 11 (blue line). From the spectral analysis of the gravitational waveform (Figure 13), it can be readily seen that it is not easy to detect these neutrino-originated GW signatures with slower temporal evolution (O  ms) by ground-based detectors whose sensitivity is limited mainly by the seismic noises at such lower frequencies. However these signals may be detectable by the recently proposed future space interferometers like Fabry-Perot type DECIGO ([215], black line in Figure 13). Contributed by the neutrino GWs in the lower frequency domains, the total GW spectrum tends to become rather flat over a broad frequency range below ~100 Hz. These GW features obtained in the context of the SASI-aided neutrino-driven mechanism are different from the ones expected in the other candidate mechanisms, such as the MHD mechanism (e.g., [216]) and the acoustic mechanism [194]. Therefore the detection of such signals could be expected to provide an important probe into the explosion mechanism (e.g., [23, 217]).

We like to draw a caution that most of the 3D models cut out the PNS and the neutrino transport is approximated by a simple light-bulb scheme [210] or by the gray transport scheme [211]. Needless to say, these exploratory approaches are but the very first step to model the neutrino-heating explosion and to study the resulting GWs. As already mentioned, the excision of the central regions inside PNSs truncates the feedback between the mass accretion to the PNS and the resulting neutrino luminosity, which should affect the features of the neutrino GWs. By the cutout, efficient GW emission of the oscillating neutron star [194] and nonaxisymmetric instabilities [184, 185] of the PNSs, and the enhanced neutrino emissions inside the PNSs [197] cannot be treated in principle. To elucidate the GW signatures in a more quantitative manner, full 3D simulations with spectral neutrino transport are apparently needed (e.g., Section 2.1). This is unquestionably a vast virgin territory awaited to be explored for the future.

2.3. Explosive Nucleosynthesis

In this section, we proceed to discuss possible signatures of supernova nucleosynthesis. The study of nucleosynthesis is of primary importance to unveil the origins of heavy elements. It could also provide a valuable information of the ejecta morphology by observing the aspherical distributions of the synthesized elements especially for a nearby CCSN event (note that nucleosynthesis is not critical for the modeling of the light-curve and spectra for the most frequent types of SNe II-P). In the following, we first present a short overview paying particular attention to explosive nucleosynthesis, and then discuss possible observational signatures that would imprint information of multidimensionalities of the supernova engine.

When in a successful explosion the shock passes through the outer shells, its high temperature induces an explosive nucleosynthesis on short timescales (e.g., [118, 119, 221], and collective references in [222]). The observational determination of the masses of the three main radioactive isotopes , , and sets one of the main constraints on the explosion dynamics, because the production of these elements is sensitive to the track of density and temperature that the expanding material traces (e.g., [223]). During the shock propagation, iron group elements such as and its daughter nucleus are predominantly produced, which are radioactive with a lifetime of 8.8 days and 111.5 days, respectively. Most CCSNe enter the so-called nebular phase after the first few months when the expanding ejecta becomes optically thin in the continuum. In the early nebular phase, is the major nuclear power source. As long as the decay particles are trapped by the ejecta, the radiation energy supplied by radioactivity is emitted instantaneously, so that the light curve can be described by an exponential decay with time, simply tracing the decay of the nuclide. To explain the bolometric light curve of SN1987A in such a phase, the mass was determined to be [224].

After several years of explosion, the radioactive output from the ejecta no longer balances with the instantaneous input by radioactivity, because the reprocessing timescale is going to be longer [225]. The bolometric light curve is affected by the delayed release of the ionization energy. After that, a self-consistent modeling is needed, in which one should include a detailed calculation of the gamma-ray/positron thermalization and a determination of the time-dependent temperature, ionization, and excitation (e.g., [225228] and references therein). Such a time-dependent modeling by Fransson and Kozma [225] revealed the mass of ~3.3 of SN1987A, which agrees well with observations (e.g., [229, 230] and collective references in [231]).

By the similar reason to just mentioned above, the determination of the is also complicated. The most recent study by Jerkstrand et al. [228] gives an estimate of the to be , which is in good agreement with the eight-year spectrum analysis of SN1987A (e.g., [232], see also [233]). As shown above, the amount of is typically one order-of-magnitude smaller than that of , however it is crucially important for young supernova remnants due to its long lifetime (~86 years). It is worth mentioning that NASA will launch the satellite NuSTAR (Nuclear Spectroscopic Telescope Array) to study production in CCSNe. The detector will be able to map out the distribution of the supernova remnant Cassiopeia A and can get velocity distributions of the in SN 1987A. By comparing detailed modeling of the SN nucleosynthesis in the context of 2D and 3D models (e.g., [234, 235]), these are expected to provide both direct probes of the explosion asymmetry.

Ever since SN1987A, challenges to the classical spherical modeling [119, 221, 236, 237] have been built also in the SN nucleosynthesis (likewise in the explosion theory and GWs mentioned so far). For many years it has been customary to simulate explosions and the effects of the shock wave on the explosive nucleosynthesis by igniting a thermal bomb in the star’s interior or by initiating the explosion by a strong push with a piston. 2D simulations with manually imparted asymmetries showed that bipolar explosion scenarios could account for enhanced synthesis along the poles as indicated in SN1987A (e.g., [35]). More recently, 3D effects have been more elaborately studied ([238], see also [239]) as well as the impacts of different explosions by employing a number of progenitors [240] or by assuming a jet-like explosion [241, 242], which is one of the possible candidates of hypernovae (e.g., [243]).

In addition to the above-mentioned work, nucleosynthesis in a more realistic simulation that model multi-D neutrino-driven explosions has been also extensively studied (e.g., [33, 43, 244] and references therein). Although a small network has ever been included in the computations, these 2D simulations employing a light-bulb scheme [33] or a more accurate gray transport scheme [43, 44] have made it possible to elucidate nucleosynthesis inside from the iron core after the shock-revival up to explosion in a more consistent manner. Kifonidis et al. [43] demonstrated that the SASI-aided low-mode explosions can naturally explain masses and distributions of the synthesized elements observed in SN1987A. Their recent 3D results by Hammer et al. [34] show that the 3D effects change the velocity profiles as well as the growth of the Rayleigh-Taylor instability, which affects properties of the SN ejecta. In simulations with spectral neutrino transport, a new nucleosynthesis process, the so-called p process, is reported to successfully explain some light proton-rich (-)nuclei including and (e.g., [245247] and references therein). It should be noted that self-consistent simulations are currently too computationally expensive to follow the dynamics of the expanding shock far outside the central iron core (≫1 s after bounce) where explosive nucleosynthesis takes place. On the other hand, the light-bulb scheme is not accurate to determine a sensitive balance between neutrino captures proceeding via or , which is decisive for quantifying the p processes. Therefore the two approaches, namely, light-bulb scheme versus self-consistent neutrino transport, may be regarded as playing a complimentary role at present.

In the next section, we briefly summarize our findings on explosive nucleosynthesis in our 2D models that utilize the light-bulb scheme to trigger explosions (e.g., [248] for more details). By changing neutrino luminosities from PNSs systematically (), we discuss how the multidimensionality formed by the SASI and convection could impact on the explosive nucleosynthesis. We employ a nonrotating 15 star with solar metalicity, which has been a best-studied supernova progenitor. The abundance pattern of the synthesized elements is estimated by a post-processing procedure, in which we have adopted a large nuclear network continuously maintained by the Joint Institute for Nuclear Astrophysics (JINA) REACLIB project [249]. Our readers might wonder why we are returning to 2D results here again after referring to 3D studies in Section 2.2. This is simply because our 3D project for the nucleosynthesis is currently in progress, and we like to note again the feature of this paper which shows only a snapshot of the moving theoretical terrain.

2.3.1. Symmetry Breaking in Explosive Nucleosynthesis

Figure 14 shows how explosion energies (left panel), explosion timescales, as well as the PNS masses (right panel) change with the input neutrino luminosity () for the 2D parametric explosion models (e.g., [248] for more details). As seen, higher makes explosion energy greater, the onset of explosion () earlier, the PNS masses () smaller. The PNS masses range in 1.54–1.62 for models with higher luminosity (, panel (b), blue line), which are comparable to the baryonic mass (around ) of neutron stars observed in binaries [250]. The explosion energies for these high luminosity models (panel (a), red line) are also close to the canonical supernova kinetic energy of  erg. To discuss explosive nucleosynthesis in the following, we thus choose to focus on the model with as a reference.

For the model, Figure 15 shows a snapshot of internal energy (left panel) and abundances of ejecta (O, Si, Fe) when the SASI-aided low-mode shock is propagating in the O-rich layers at around ~20000 km along the pole. Ejecta with higher internal energies (left panel, colored by red) concentrate on the shock front particularly in polar regions (), where both Si and O burning proceed to produce abundantly and . These aspherical element distributions might be responsible for anisotropies of the SN ejecta as observed in SN1987A. More recently, Kimura et al. [251, 252] have analyzed the metal distribution of the Cygnus loop and pointed out that the progenitor of the Cygnus loop is a CCSN explosion whose progenitor mass ranges in ~12–15. Since the material in this middle-aged supernova remnant has not been completely mixed yet, the observed asymmetry is considered to still remain a trace of inhomogeneity produced at the moment of explosion. The asymmetries of the SN ejecta obtained in our 2D explosion model of a 15 progenitor have a close correlation with the ones in the Cygnus loop (see [248] for more detailed comparison). This might allow us to speculate that globally asymmetric explosions induced by the SASI-aided neutrino-heating mechanism could account for the origin.

Concerning the synthesized amount of and , they are ~0.06 and , respectively, for the model (bottom panel in Figure 16). For SN1987A, masses of and are deduced to be ~0.07 [27, 253] and 1−2 ([36], and references therein). The obtained mass of is, therefore, not enough for SN1987A as well as for Cas A () [254]. The shortage of is a long-standing problem, which might be solved by 3D effects [34]. For the abundance patterns, our reference model (top and middle panels in Figure 16) is shown to reproduce a similar trend to the solar system (top panel, Figure 16), which is in sharp contrast to the model with lower neutrino luminosity (middle panel, Figure 16).

Finally, we note that there should be at least two barriers to get over in the current-generation studies of the CCSN nucleosynthesis. As mentioned already, the first one is the need of radiation hydrodynamic simulations that follow the dynamics from gravitational collapse, shock-revival, shock propagation in the stellar mantle, to stellar explosions in a self-consistent manner. The second one is the need of computing light curves, spectra, and the degree of polarizations in 2D or 3D hydrodynamic models, for which one needs to solve the multi-D photon transport coupled with multi-D hydrodynamics. Unfortunately however, the marriage of the two items, albeit ultimately needed for clarifying the photon messengers, may not be done immediately due to their computational expensiveness. For example, one needs to follow the dynamics more than ~1 day after the onset of gravitational core-collapse for computing the explosive nucleosynthesis, because the envelop of a typical supernova progenitor extends up to a radius of ~1012 cm. It takes furthermore more than ~1 week before the shock propagation enters to the so-called homologous phase. At present, the best available numerical simulation to this end is limited to 2D (e.g., [244]). For discussing anisotropies in emission-line profiles, one needs to follow the dynamics later than ~1 year after explosions when the remnant becomes transparent. It is a very challenging (but very important) task for the first principle simulations to overcome the big gaps regarding the very different timescales. An encouraging news is that several groups are pursuing it with the use of advanced numerical techniques such as the adaptive-mesh-refinement approach [116] or Ying-Yang gridding [115]. With an accelerating power (like peta- or exa-scale) of supercomputers, it would be unsurprising that the next (or the next after next!) generation supernova modellers can execute the ultimate ab initio simulations, whose findings would be immediately used to make a more precise prediction of the photon messengers.

3. MHD Mechanism

In this section, we move on to focus on the MHD mechanism. After we briefly summarize the current status of this topic below, we mention the explosion dynamics that proceeds by the field-wrapping mechanism in Section 3 and discuss properties of the so-called magnetorotational instability (MRI) in Sections 3 and 3 based on our preliminary results. Possible signatures of GWs and neutrinos are discussed in Sections 3.1 and 3.2, respectively.

Numerical simulations of MHD stellar explosions have started already in early 1970s shortly after the discovery of pulsars [255258]. However, it is rather only recently that the MHD studies come back to the front-end topics in the supernova research followed by a number of extensive MHD simulations (e.g., [46, 176, 207, 259270] for references therein). Main reasons for this activity are observations indicating very asymmetric explosions [38, 40], and the interpretation of magnetars [271, 272], collapsars and their relevance to gamma-ray bursts [273275] as a possible outcome of the magnetorotational core-collapse of massive stars.

The MHD mechanism of stellar explosions relies on the extraction of rotational free energy of collapsing progenitor core via magnetic fields. Hence a high angular momentum of the core is preconditioned for facilitating the mechanism [276]. Given (a rapid) rotation of the precollapse core, there are at least two ways to amplify the initial magnetic fields to a dynamically important strength, namely, by the field wrapping by means of differential rotation that naturally develops in the collapsing core, and by the MRI (MRI, see [16]). Each of them we are going to review briefly in the following.

Field-Wrapping Mechanism
Three dimensional plots of Figure 17 are useful to see how the field wrapping occurs. For the 2D model taken from our 2D special relativistic (SR) MHD simulation [268], the precollapse magnetic field is set to be as high as  G that is uniform and parallel to the rotational axis and the initial parameter (ratio of initial rotational energy to the absolute value of the initial gravitational energy) is taken to be 0.1% with a uniform rotation imposed in the iron core. From the left panel, it can be seen that the field lines are strongly twisted around the rotational axis. From the induction equation of ideal MHD equations, the time evolution of the toroidal fields () can be expressed as (e.g., [276]),
where denotes the distance from the rotational axis and represents the component of the poloidal fields in the cylindrical coordinates (). Given that a precollapse iron core has strong magnetic field (such as with representing the initial poloidal fields) and rapid rotation (such as with being the precollapse rotation period), the typical amplification of near core bounce may be estimated as
where represents the typical timescale when the field amplification via the wrapping processes saturates, denotes the rotational period near bounce. Equation (2) shows that the toroidal fields grow linearly with time by wrapping the poloidal fields (), which is the so-called dynamo (e.g., [265]). The  G class magnetic fields are already dynamically important strength to affect the postbounce hydrodynamics as will be mentioned below.
The white lines in the right panel in Figure 17 show the streamlines of matter. A fallback of material just outside of the jet head advecting downwards to the equator (like a cocoon) is seen. In this jet with a cocoon-like structure, the magnetic pressure is always dominant over the matter pressure. As estimated above, the typical field strength behind the shock is ~1015. The critical strength of the toroidal magnetic field to induce the magnetic shock revival (i.e., the revival of the stalled bounce shock due to the MHD mechanism) is estimated as follows. The matter in the regions behind the stalled shock is pushed inwards by the ram pressure of the accreting matter. This ram pressure is estimated as,
where and denote typical density and radial velocity in the vicinity of the stalled shock. When the toroidal magnetic fields are amplified as large as ~1015 due to the field wrapping behind the shock, the resulting magnetic pressure, , can overwhelm the ram pressure, leading to the magnetic shock revival (e.g., Figure 18). The onset timescale of the magnetic shock revival is sharply dependent on the precollapse magnetic fields and rotation (e.g., [46]). As the initial field strength is larger with more rapid rotation imposed, the interval from the shock stall to the MHD-driven shock revival becomes shorter. The speed of the jet head is typically mildly relativistic (at most ~0.4), with being the speed of light (see also the right panel of Figure 18). At the shock breakout of the iron core, the explosion energies generally exceed  erg for models that produce MHD explosions in a shorter delay after the stall of the bounce shock. These features are in accord with those obtained in more detailed MHD simulations with spectral neutrino transport (e.g., [46, 273]) (see, e.g., [216] for a more detailed comparison).

On the MRI
We are now in a position to discuss possible impacts of the MRI in supernova cores. Akiyama et al. [279] were the first to point out that the interfaces surrounding the nascent PNSs quite generally satisfy the MRI instability criteria. Therefore any seed magnetic fields can be amplified exponentially in the differentially rotating layers, much faster than the linear amplification due to the field wrapping (e.g., (2)). After the MRI enters to the saturated state, the field strength might reach  G, which is high enough to affect the supernova dynamics. Not only in the field amplification, the MRI plays a crucial role also in operating the MHD turbulence (see [16, 280, 281]). The turbulent viscosity sustained by the MRI can convert a fraction of the shear rotational energy to the thermal energy of the system. Thompson et al. [141] suggested that the additional energy input by the viscous heating can help the neutrino-driven supernova explosion. Followed by the exponential field amplification and additional heating, a natural outcome of the magnetorotational core-collapse is expected to be the formation of energetic bipolar explosions, which might be observed as hypernovae (e.g., Section 2.3).
Note again that bipolar explosions obtained in the previous section with the assumption of a strong precollapse magnetic field (1012 G), are predominantly driven by the field wrapping processes, not by the MRI. Shibata et al. [282] pointed out in their fully 2D GRMHD core-collapse simulations that more than 10 grids are at least needed to capture the fastest growing mode of the MRI (see also [283]). As well-known, the growth rate of MRI-unstable modes depend on the product of the field strength and the wave number of the mode [16]. When a rapidly precollapse core is strongly magnetized as  G (like in the case of collapsar progenitor), the wavelength is on the order of km in the postbounce core. In this case, the exponential field growth of the MRI was successfully captured by high-resolution simulations [282]. On the other hand, the fastest growing modes drop below several meters [284] for canonical supernova progenitors, which rotates much more slowly with weaker initial fields (~109 G) as predicted by recent stellar evolution models [126, 127]. At present, it is computationally too expensive to resolve those small scales in the global MHD simulations, typically more than two or three orders-of-magnitudes smaller than their typical finest grid size. To reveal the nature of the MRI, local simulations focusing on a small part of the MRI-unstable regions are expected to be quite useful as traditionally studied in the context of accretion disks (see [16]).
Obergaulinger et al. [284] were the first to study the growth and saturation level of the MRI in the supernova environment (see, [285] for extensive simulations that focus on the MRI in relativistic outflows in the context of gamma-ray bursts and active galactic nuclei). To ease a drawback of the local shearing box simulation, they employed the shearing disk boundary conditions by which global radial density stratification can be taken into account. By performing such a semiglobal simulation systematically in 2D and 3D, they derived scaling laws for the termination of the MRI. As estimated in [279], the MRI was shown to amplify the seed fields exceeding  G. These important findings may open several questions that motivate us to join in this effort, such as how the nonlinear properties as well as the scaling laws could be in the subsequent nonlinear state and whether the viscous and resistive heating driven by the MRI turbulence could or could not affect the supernova mechanism.
In the following, we briefly summarize the current status of our MRI project, in which we perform a local shearing box simulation to study the MRI-driven turbulence and their nonlinear properties bearing in mind the application to the supernova core (Masada, Takiwaki, and Kotake in preparation). The final goal of this project is to determine how much the conversion from the rotational energy (that the differential rotation taps) to the thermal energy occurs via the MRI turbulence, which was treated parametrically in [141]. As will be shown later, our preliminary results suggest that the additional energy supply would affect the neutrino-driven explosion in the case of rapidly rotating cores at late postbounce phase.
To study the nonlinear properties of the MRI, we solve viscous MHD equations by a finite-difference code that was originally developed by Sano et al. [286]. The hydrodynamic part is based on the second-order Godunov scheme [287], which consists of Lagrangian and remap steps. The exact Riemann solver is modified to account for the effect of tangential magnetic fields. The field evolution is calculated with Consistent MoC-CT method [288]. The energy equation is solved in the conservative form and the viscous terms are consistently calculated in the Lagrangian step. The advantages of our scheme are its robustness for strong shocks and the satisfaction of the divergence-free constraint of magnetic fields [289, 290].
Rigorously speaking, 3D simulations are absolutely required for the study of the MRI because only by them the disruption of channel structures, mode couplings, and the growth of parasitic instabilities can be accurately determined (e.g., [291, 292]). But in the following, we are only able to show our 2D results assuming axisymmetry, which are currently being updated to 3D. Bearing these caveats in mind, we briefly illustrate fundamental properties of the MRI (which is not always familiar with supernova modellers,) and discuss their possible impacts on the explosion dynamics.

Local Simulations in Axisymmetry
As it is well known, there are three typical evolutionary stages observed in our local simulations. They are (i) exponential growth stage, (ii) transition stage, and (iii) nonlinear turbulent stage [280, 293]. Each stage is denoted by white, dark gray, and light gray shaded regions in Figure 19.
In the stage (i), the channel structure of the magnetic field exponentially evolves and inversely cascades to the larger spatial scale with small structures merging as the magnetic field is amplified. This is because the channel modes of the MRI are an exact solution for the nonlinear MHD equations [294]. The temporal evolution of channel structures for the toroidal magnetic field and their inversely cascading nature are shown in panels (a) and (b) of Figure 20. Then in the stage (ii), the channel structure of the magnetic field is disrupted via the parasitic instability and magnetic reconnection at the transition stage [284, 294]. The channel disruption induces a drastic phase-shift from the coherent structure to the turbulent tangled structure of the field as is illustrated by panel (c) of Figure 20. The magnetic energy stored in the amplified magnetic field is then converted to the thermal energy of the system. Finally in the stage (iii), the nonlinear stage emerges after the channel structures are disrupted to be a turbulence state driven by the MRI (see panel (d) in Figure 20). Figure 19(a) shows that, at this stage, the turbulent Maxwell stress () dominates over the Reynolds stress (). Therefore the turbulent Maxwell stress is the main player to convert the free energy stored in the differential rotation into the thermal energy of the system.
Figure 19(b) shows the power spectra of the magnetic energy in the nonlinear turbulent stage. The logarithmic slope of the spectrum roughly coincides with in the regime while it steepens at higher wave-numbers where numerical effects become important (like ). This indicates that the turbulent state in the MRI can be represented by the Kolmogorov spectrum for a homogeneous, incompressible turbulence that scales as . These features would be consistent with the 3D properties of the MRI-driven turbulence [280]. It has been pointed out in the previous axisymmetric local shearing box simulation that the channel structures do neither disrupt nor transit to the turbulent state without a relatively strong resistive effect [284, 295]. We however adopted the simulation box with a large aspect ratio () and quasi-isothermal equation of state (), both of which could facilitate to disrupt the channel structures [296] as in 3D simulations.
Figure 21 shows the time- and volume-averaged Maxwell (green squares) and Reynolds stresses (orange diamonds, panel (a)), and magnetic energy (blue circles, panel (b)) in the turbulent stage as a function of the initial magnetic energy. It is shown that the Maxwell stress depends on the initial magnetic flux and it provides the following scaling relation:
with . This suggests that the magnetic flux initially penetrating the system controls the nonlinear transport properties of the MRI-driven turbulence. This is qualitatively consistent with the previous results in the accretion disk [280, 286, 297, 298] while the power-law index obtained here is slightly larger. Note that the magnetic energy at the nonlinear turbulent stage has a linear dependence on the initial magnetic flux (because (Figure 21(b))). This is also qualitatively consistent with Hawley et al. [280]. In this way, we investigated parameter dependence of the rotational shear ( parameter) as well as the initial pressure () on the saturated value of the Maxwell stress and deduced a scaling relation that yields . Here the shear parameter represents the degree of differential rotation as with being the radial coordinate.
Using the scaling relation obtained in our local simulations, the heating rate maintained by the MRI turbulence can be estimated as:
where represents a typical amplification factor of the stress normalized by 30 (e.g., Figure 21), is the shear rate normalized by unity, and is the ratio of the vorticity to the shear normalized by unity. Since we are interested in the growth of the MRI in the postbounce phase, the original magnetic flux, gas pressure angular velocity, and the shear rate , , , and are replaced by , , and with some amplification factors.
A typical volume in which the MRI-driven turbulence is active, may be estimated as , where  cm is typical radius of the PNS normalized by cm, and  cm is the typical radial width of the MRI-active layer normalized by  cm. Then the energy releasing rate becomes:
Typical timescales of neutrino-driven explosions observed in recent supernova simulations are 400 ms after bounce (e.g., [125, 138]). So the energy deposition could be as high as if the core rotates rather rapidly as taken in the above estimation. If this is the case, the turbulent heating due to the MRI is expected to play an important role for assisting the neutrino-driven explosion.
It should be noted that nonaxisymmetric modes of the parasitic instability could disrupt the channel structures more effectively than in 2D [291, 292]. To obtain more accurate estimate of the amplification factor of (6), 3D simulations are unquestionably required, which we are currently undertaking [299].

3.1. Gravitational Waves

As already mentioned in Section 2.2, rapid rotation, necessary for producing strong bounce GW signals, is likely to obtain ~1% of massive star population (e.g., [62]). This can be really the case, albeit minor, for progenitors of rapidly rotating metal-poor stars, which experience the so-called chemically homogeneous evolution [274, 275]. In such a case, the MHD mechanism could work to produce energetic explosions, which are receiving great attention recently as a possible relevance to magnetars and collapsars (e.g., [300, 301], see [273, 302, 303] for collective references), which are presumably linked to the formation of long-duration gamma-ray bursts (GRBs) (e.g., [304] for review).

Among the previous studies focusing on the bounce signals (e.g., references in Section 2.2), only a small portion of papers has been spent on determining the GW signals in the MHD mechanism [176, 260, 264, 265, 269, 282]. This may be because the MHD effects on the dynamics as well as their influence over the GW signals can be visible only for cores with precollapse magnetic fields over  G [176, 264]. Considering that the typical magnetic-field strength of GRB progenitors is at most  G [274], this is already an extreme situation. In a more extremely case of  G, a secularly growing feature in the waveforms was observed [263, 269, 282]. In the following, we summarize the GW signatures based on our 2D SRMHD simulations [216] and discuss how a peculiarity of the MHD mechanism could be imprinted in the GW signals.

The left panel of Figure 22 shows the gravitational waveform with the increasing trend, which is obtained for a model with strong precollapse magnetic field ( G) also with rapid rotation initially imposed ( parameter = 0.1%). Such a feature cannot be observed for a weakly magnetized model ( G, right panel). To understand the origin of the increasing trend, it is straightforward to look into the quadrupole GW formula, which can be expressed as where is (e.g., [306]). can be expressed as On the right hand side, the first term is related to anisotropic kinetic energies which we call as hydrodynamic part, the second term is related to anisotropy in gravitational potentials which we call as the gravitational part, and finally the third term is related to anisotropy in magnetic energies that we refer to as magnetic part,

The right panel of Figure 23 shows contributions to the total GW amplitudes (6) for the strongly magnetized model (left panel of Figure 22), in which the left-hand-side panels are for the sum of the hydrodynamic and gravitational part, namely (left top(+)/bottom(−) (11), (13), and the right-hand-side panels are for the magnetic part, namely (right top(+)/bottom(−)) (e.g., (14)). By comparing the top two panels in Figure 19, it can be seen that the positive contribution is overlapped with the regions where the MHD outflows exist. The major positive contribution is from the kinetic term of the MHD outflows with large radial velocities (e.g., in (11)). The magnetic part also contributes to the positive trend (see top right-half in the right panel (labeled by mag(+))). This comes from the toroidal magnetic fields (e.g., in (14)), which dominantly contribute to drive MHD explosions.

Figure 24 shows the GW spectra for a pair models of B and B like in Figure 22 that with or without the increasing trend. Regardless of the presence (left panel) or absence (right panel) of the increasing trend, the peak amplitudes in the spectra are around kHz. This is because they come from the GWs near at bounce when the MHD effects are minor. On the other hand, the spectra for lower frequency domains (below ~100 Hz) are much larger for the model with the increasing trend (left panel) than without (right panel). This reflects a slower temporal variation of the secular drift inherent to the increase-type waveforms (e.g., Figure 22). It is true that the GWs in the low frequency domains mentioned above are relatively difficult to detect due to seismic noises, but a recently proposed future space interferometers like Fabry-Perot type DECIGO is designed to be sensitive in the frequency regimes [215, 278] (e.g., the black line in Figure 24). Our results suggest that these low-frequency signals, if observed, could be one important messenger of the increase-type waveforms that are likely to be associated with MHD explosions exceeding  erg.

3.2. Neutrino Signals

As of now, SN1987A remains the only astrophysical neutrino source outside of our solar system. Even though the detected events were just two dozen, these events have been studied extensively (yielding ~500 papers) and have allowed us to have a confidence that our basic picture of the supernova physics is correct (e.g., [4], see [5, 6] for a recent review). For a next nearby event, it is almost certain that the flagship detectors like Super-Kamiokande and IceCube are able to detect a SN “neutrino light curve” with high statistics. For the detectors, the horizon to the sources now extends out to about 100 kpc and thus covers our galaxy and its satellites. A future megaton-class detector reaches as far as Andromeda galaxy at a distance of 780 kpc (e.g., Hyper-Kamiokande, Memphys, and LBNE) and large-scale scintillator (e.g., HALO [307] that is an upgrade of SNO) or liquid-Argon detectors (GLACIER [308] that is an upgrade of ICARUS) will also play an important role (see [309]) for a complete list and details).

Before core bounce, neutrinos of different flavors are initially tapped in their relative neutrino spheres with increasing their Fermi energies as the central density becomes higher. After bounce, the shock dissociates nuclei into free nucleon, which drastically increases the number of protons, leading to a sharp increase in electron capture and production of the prompt burst. At this neutronization epoch (in the first 10–20 ms postbounce), the largest difference in the light-curves among and the rest of the neutrino species can be seen (e.g., Figure 1 in [310]). Since neutrino oscillations take place only when there is a difference in the neutrino fluxes among different species (see Section 3 [22] for more details on the neutrino oscillations), the oscillation effects in the neutronization epoch would be strongest. Unfortunately, however, they are very hard to measure because the currently-running detectors are primarily sensitive to signals that are produced by the inverse beta process . It should be noted that the detection the bursts is important for studying not only neutrino oscillations but also for determining the direction and distance to the source [310]. In a megaton-class Čherenkov detector in which gadolinium is added to catch neutrons [311], and large liquid-Argon detectors (like ICARUS and GLACIER cited above) are expected to be powerful detectors.

During the subsequent accretion phase (approximately before 1 s postbounce), the bounce shock stagnates and matter falls in through the stalled shock to the center, releasing gravitational energy that powers neutrino emission. For relatively well-studied progenitors in the mass range of 10.8 to 25 (e.g., [96, 121, 310, 312, 313] and references therein), one generally finds and with and representing luminosity and average neutrino energy for the neutrino species of (; and their antiparticles). The accretion phase is followed by the PNS (protoneutron star) cooling in which the accretion stops and the PNS settles to become a neutron star after a successful revival of the stalled bounce shock into explosions. Since the explosion mechanism is still under debate (see discussion in Section 2), it is at present difficult to tell the accurate turn-over time. But it should be less than 1-2 s postbounce, otherwise the central PNS would collapse into a black hole [314]. In the PNS cooling phase, the luminosities of all neutrino species become similar (with the energy hierarchy of ) and decrease monotonically with time. Due to the similarity in the luminosities and spectra among and and also due to the darkening with time, it is much harder to see flavor oscillation effects by the currently running detectors.

Neutrinos streaming out of the iron core interact with matter via the Mikheyev-Smirnov-Wolfenstein (MSW) effect (e.g., [78]) firstly in propagating through progenitor envelope and then through the Earth before reaching detectors. Such effects have been extensively investigated so far from various points of view, with a focus such as on the progenitor dependence of the early neutrino burst (e.g., [310, 315]) and on the Earth matter effects (e.g., [316, 317]). The conversion efficiency via the MSW effect depends sharply on density, equivalently electron fraction gradients, thus sensitive to discontinuities produced by the passage of supernova shocks. Such shock effects have been also extensively studied (e.g., [318322], see [79] for a recent review). The shock passage to the so-called high-resonance regions may be observed as a sudden decrease in events in the case of normal mass hierarchy (or in the case of inverted mass hierarchy). Such features monitoring time evolution of the density profile like a tomography, thus could provide a powerful test of the mixing angle and the mass hierarchy (e.g., [320, 323325]).

It is rather recently that the importance of collective flavor oscillations was widely recognized. Neutrinos streaming out of the neutrino spheres are so dense that they provide a large matter effect for each other. The collective effect usually takes place between the neutrino sphere and the mentioned MSW region. Regardless of the inherent nonlinearity and the presence of multiangle effects, the final outcome for the emergent neutrino flux seems to be converging after a series of extensive study over the past years at least for the 1D models. Although the collective effects will not be significant to assist the neutrino-driven explosions (e.g., [326, 327], see however [328, 329]), they are predicted to emerge as a distinct observable feature in their energy spectra (see [7477] for reviews of the rapidly growing research field and collective references therein).

An important lesson from SN1987A is that for explaining the duration of the events, there was no other energy-loss channel but for the ordinary neutrinos in the context of the Standard Model of particle physics (e.g., [6, 330]). The next SN event could provide an opportunity to study also a nontrivial property of neutrinos, such as the magnetic dipole moments. The resonant spin-flavor conversion has been also studied both analytically (e.g., [331333] and references therein) and numerically (e.g., [334]), which can transform some of the prompt burst into in highly magnetized supernova envelopes, leading to a huge burst. The mentioned important ingredients related to the flavor conversions due to the MSW effects, collective effects, and the electromagnetic effects in the supernova environment have been studied often one by one in each study without putting all the effects together, possibly in order to highlight the new ingredient. In this sense, all the studies mentioned above should be regarded as complimentary towards the precise predictions of SN neutrinos.

Here it should be mentioned that most of those rich phenomenology of SN neutrinos have been based on 1D spherically symmetric simulations (e.g., [75, 310, 312, 315, 319, 320, 335338] and references therein). Apart from a simple parametrization to mimic anisotropy and Stochasticity of the shock (e.g., [339341]), there have been not so many studies focusing how global asymmetry of the explosion dynamics obtained in recent supernova simulations (e.g., Sections 2 and 3) would affect the neutrino oscillations (see, however, [79, 342]). This is mainly due to the lack of multi-D supernova models, which are very computationally expensive to continue the simulations till the shock waves propagate outward until they affect neutrino transformations. It is only recently that several studies along this line have been reported [133, 343], in which the collective or MSW effects are rendered to be treated in an approximate manner. By analyzing the 2D results of a 15 model by Marek and Janka [125] who included one of the best available neutrino transfer approximations (e.g., Table 1), Lund et al. [343] pointed out that fast time variations caused by convection and SASI lead to significant modulations around a few hundred Hz, which can be visible in IceCube or future megaton-class detectors for the galactic SN source. Based on the 2D results of 20 models by Ott et al. [132] who reported the first 2D multiangle transport simulations, Brandt et al. [133] also obtained the similar results. From their rapidly rotating model, they also pointed out that a rapid rotation of precollapse SN cores imprints strong asymmetries in the neutrino flux [344] as well as in its light curves [133, 345].

In the following, we briefly summarize our findings [305] in which we studied exploratory the neutrino oscillations in the context of the MHD mechanism. As mentioned in Section 3, it is recently possible for special relativistic simulations to follow the dynamics of MHD explosions continuously, starting from the onset of gravitational collapse, through core-bounce, the magnetic shock-revival, till the shock-propagation to the stellar surface. Based on our models [268], we calculated numerically the neutrino flavor conversion in the highly nonspherical envelope through the MSW effect. The neutrino transport was simply treated by a leakage scheme (e.g., [268]). The emergent neutrino spectra was assumed to take a Fermi-Dirac type distribution function and the neutrino temperature was estimated to take the average matter temperature on the neutrino sphere. As explained in Section 3, the density profile along the polar and equatorial direction is very different due to the strong explosion anisotropy inherent to the MHD explosions. In the following, we pay attention to the anisotropic shock effect on the MSW effect. As will be discussed, we could observe a sharp dip in the neutrino event only seen from a polar direction, albeit depending on the mass hierarchy and the mixing angle of . An advantage of the MHD models is that the shock revival can occur much faster after bounce compared to the other proposed mechanism. Therefore, the neutrino luminosity could remain higher than those in the other mechanisms, which could potentially enhance the detectability due to the early shock arrival to the resonance region. For simplicity, the effects of neutrino self-interactions are treated very phenomenologically and the resonant spin-flavor is not considered. Even though far from comprehensive in this respect, we presented the first discussion how the magneto-driven explosion anisotropy has impacts on the emergent neutrino spectra and the resulting event number observed by the SK for a future Galactic supernova (e.g., [305]).

3.2.1. Possible Neutrino Signatures from MHD Explosions

One prominent feature of the MHD models is a high degree of the explosion asphericity. Figure 25 shows several snapshots featuring typical hydrodynamics of the model, from near core-bounce (0.8 s, left), during the shock-propagation (2.0 s, middle), till near the shock breakout from the star (4.3 s, right), in which time is measured from the epoch of core bounce. It can be seen that the strong shock propagates outwards with time along the rotational axis. On the other hand, the density profile hardly changes in the equatorial direction, which is a generic feature of the MHD explosions.

As well known, the flavor conversion through the pure-matter MSW effect occurs in the resonance layer, where the density is where is the mass squared difference, is the neutrino energy, is the number of electrons per baryon, and is the mixing angle. Since the inner supernova core is too dense to allow MSW resonance conversion, we focus on two resonance points in the outer supernova envelope. One that occurs at higher density is called the H-resonance, and the other, which occurs at lower density, is called the L-resonance. and correspond to and at the H-resonance and to and at the L-resonance.

Figure 26 are evolutions of the density profiles in the polar direction every 0.4 s as a function of radius. Above and below horizontal lines show approximately density of the resonance for different neutrino energies which are 5 and 60 MeV, respectively. Blue lines (left panel) show the range of the density of the H-resonance, and sky-blue lines (right panel) show that of the L-resonance. Along the polar axis, the shock wave reaches to the H-resonance, ~O )g/ at ~0.5 s, and the L-resonance, ~O  g/cm3 at ~1.2 s. It should be noted that those timescales are very early in comparison with the ones predicted in the neutrino-driven explosion models, typically ~5 s and ~15 s for the H- and L-resonances, respectively (e.g., [320, 325]). This arises from the fact that the MHD explosion is triggered promptly after core bounce, which is in sharp contrast to the neutrino-driven delayed explosion models [320, 325]. The progenitor of the MHD models, possibly linked to long-duration gamma-ray bursts, is more compact due to a chemically homogeneous evolution [274, 275], which is also the reason for the early shock-arrival to the resonance regions.

Figure 27 is the energy spectra of for (a), (b), and (c) in the case of inverted mass hierarchy. The solid lines (red lines) and the dotted lines (blue lines) are the spectra of the polar direction and the equatorial direction, respectively. At the panel (a), an enhancement of the low-energy side of the polar direction is seen, which is as a result of the supernova shock reaching to the H-resonance region. The survival probability of can remain non-zero when the steep decline of the density at the shock front changes the resonance into non-adiabatic (see [305] for more details). As the shock propagates from the high density region to low energy region, the shock effect transits from low-energy side to high-energy side in spectra, because the density at the resonance point, , is proportional to . In fact, the energy spectrum of the polar direction becomes softer in the high-energy side than for the equatorial direction as shown in Figure 27(b).

Figure 28 shows the time evolution of the event number in SK, where the solid line (red line) and the dotted line (blue line) are for the polar and the equatorial direction, respectively. The shock effect is clearly seen. The event number of the polar direction shows a steep decrease, marking the shock passage to the H-resonance layer. The change of the event number by the shock passage is about 36% of the event number without the shock. Since the expected events are ~2500 at the sudden decrease (Figure 28), it seems to be quite possible to identify such a feature by the SK class detectors. Such a large number imprinting the shock effect, possibly up to two-orders-of magnitudes larger than the ones predicted in the neutrino-driven explosions (e.g., Figure  11 in [320]), is thanks to the mentioned early shock-arrival to the resonance layer, peculiar for the MHD explosions.

Finally we exploratory discuss possible impacts of the spectral swap between and , which is one of a probable outcome of the collective neutrino oscillations [346]. In our model, the original neutrino flux of is larger than that of for the low-energy side, which is vice versa for the high-energy side. Due to the spectral swap, the above inequality sign reverses (see [305] for more detail). As a result, the event number in the polar direction is expected to increase by the influence of the shock. To draw a robust conclusion to the self-interaction effects, one needs to perform a more sophisticated analysis such as the multiangle approach (e.g., [76]), which is a major undertaking.

4. Summary and Discussion

The aim of writing this paper was to provide an overview of what we currently know about the explosion mechanism, neutrinos, gravitational waves, and explosive nucleosynthesis in CCSNe to bridge theory and observation through multimessenger astronomy. Not to inflate the volume, we had to limit ourselves to focus primarily on the findings based on our numerical studies, which are thus limited to GWs (Section 2.2) and explosive nucleosynthesis (Section 2.3) in the neutrino-heating mechanism (Section 2), and to GWs (Section 3.1) and neutrino signals (Section 3.2) in the MHD mechanism (Section 3). Apparently our current presentation is weak on the photon side mainly because there is a big gap to connect the outcomes obtained in the multi-D neutrino-radiation-hydrodynamic simulations covering only the very centre of the star, to the light-curve and spectra modeling that requires multi-D photon-radiation-hydrodynamic modeling (far) after the shock-breakout of the star. The gap is not only the matter of physical scale, but also the matter of numerical difficulty. Note that Table 2 is not intended to show an overview covering everything concerning the SN multimessengers, which is far beyond the scope of this review. We boldly present the incomplete figure here, hoping that the shortage might give momentum to theorists for fixing the problems and for gaining much more comprehensive multimessenger perspectives for the future.

In the neutrino-heating mechanism, stellar rotation holds a key importance to characterize the SN multimessengers. The whole story may be summarized as follows. If a precollapse iron core has a “canonical” rotation rate as predicted by recent stellar evolution calculations [126, 127], collapse-dynamics before bounce proceeds spherically and the postbounce structures interior to the PNS are essentially spherical. Typically later than ~100 ms after bounce, the bounce shock turns to a standing accretion shock, and the activities of convective overturns and SASI become vigorous with time behind the standing shock. Since anisotropies of the neutrino flux and matter motions in the postbounce phase are governed by the nonlinear hydrodynamics, GWs emitted in this epoch also change stochastically with time (indicated by “Stochastic” in the table, e.g., Section 3.1 for more details). For detecting these signals for a Galactic CCSN event with a good signal-to-noise ratio, we need next generation detectors such as the advanced LIGO, LCGT, and the Fabry-Perot (FP-) type DECIGO. Neutrino signals at the nonlinear epoch also change stochastically with time, which is expected to be detected by currently running detector such as by IceCube and SK [133, 197, 343] for the galactic event. It should be also mentioned that detailed neutrino transport including the most up-to-date interaction rates should be accurately implemented in the multi-D simulations to make a reliable prediction of the neutrino signals (e.g., [312, 313]). When material behind the stalled shock successfully absorbs enough neutrino energy to be gravitationally unbound from the iron core, the p process is expected to successfully explain some light proton-rich nuclei in the neutrino-driven wind phase [245247]. When the revived supernova shock (successfully) attains the kinetic energy as energetic as  erg, explosive nucleosynthesis that proceeds in the SASI-aided low-mode explosions is expected to partly explain the solar abundance yields as well as the observed nonuniform morphology of the synthesized elements (e.g., Section 2.3.1).

One important notice here is that explosion energies obtained in the state-of-the-art multi-D simulations are typically underpowered by one or two orders-of-magnitudes to explain the canonical supernova kinetic energy (~1051 erg, e.g., Table 1). Moreover, the softer EOS, such as of Lattimer and Swesty [347] (LS) EOS with an incompressibility at nuclear densities, , of 180 MeV, is employed in those simulations (e.g., in multi-D simulations of the MPA, Oak Ridge+, and Tokyo+ groups in Table 1). On top of a striking evidence that favors a stiffer EOS based on the nuclear experimental data (MeV, [348]), the soft EOS may not account for the recently observed massive neutron star of ~2 [349] (see the maximum mass for the LS180 EOS in [314, 350]). With a stiffer EOS, the explosion energy may be even lower as inferred from Marek and Janka [125] who did not obtain the neutrino-driven explosion for their model with  MeV (on the other hand, they obtained 2D explosions for Shen EOS ( MeV) (Janka, private communication)). What is then missing furthermore? We may get the answer by going to 3D simulations (Section 2.1) or by taking into account new ingredients, such as exotic physics in the core of the protoneutron star [123, 351], viscous heating by the magnetorotational instability (Section 3), or energy dissipation via Alfvén waves [352]. In addition to the 3D effects, GR is expected to help the onset of multi-D neutrino-driven explosions [170, 353] (e.g., [129] in 2D simulations but with detailed neutrino transport or [170] in 3D simulations but with a very approximate neutrino transport). At present, the most up-to-date neutrino transport code can treat the multienergy and multiangle transport in 2D [132] and even in 3D simulations [159] but only limited to the Newtonian simulations. Extrapolating a current growth rate of supercomputing power, the dream simulation (i.e., 6D simulations with full Boltzmann neutrino transport in full GR simulations) is likely to be practicable not in the distant future (in our lifetime!) presumably by using the exa-scale platforms.

In addition to these numerical advancements, the physical understanding of the supernova theory via analytical studies is also in a steady progress. What determines the saturation levels of SASI? A careful analysis on the parasitic instabilities has been reported to answer this question [354]. What determines mode couplings between small-scale convection eddies and large-scale SASI modes? To apply the theory of turbulence [355] should put a milestone to address this question. Besides the two representative EOSs of Lattimer-Swesty and Shen, new sets of EOSs have been recently reported [356358]. Some of them will enable SN modellers to go beyond a single nucleus approximation, which is an important improvement for an accurate description of neutrino-nucleus interaction.

Rotation, albeit depending on its strength, can give a special direction (i.e., spin axis) in the supernova cores. As discussed in Section 2.2.2, stochastic nature of the GWs becomes weak in the presence of rotation as a result of the spiral SASI (It should be noted that the feature can be seen even in a slowly rotating progenitor that agrees with recent stellar evolution calculation [126, 127] (see Section 2.2.2 for more details)) (indicated by “excess for equator” in Table 2). Although it is not easy to detect these low-frequency GWs from anisotropic neutrino emission by currently running laser interferometers, a recently proposed space-based interferometers (like FP-type DECIGO) would permit detection for a galactic event. Contributed by the neutrino GWs in the lower frequency domain, the total GW spectrum is expected to become rather flat over a broad frequency range below ~1 kHz. This GW feature obtained in the context of the SASI-aided neutrino-driven mechanism is different from the one in the other candidate supernova mechanisms, such as the MHD mechanism (Section 3.1) and the acoustic mechanism [194], thus could provide an important probe into the explosion mechanism. Rapid rotation induces a polar excess of neutrino emission with a harder spectrum than on the equator [133, 344, 345] (indicated by “polar excess” in the table). This should be also the case of the MHD mechanism, which deserves further investigation for a quantitative discussion (symbolized by Polar excess (?) in the table). These (rotation-induced) neutrino signatures, if observed, will carry an important information about the angular momentum profiles hidden deep inside massive stellar cores. Concerning the p process and explosive synthesis in rapidly rotating cores (as well as in the MHD explosions), there still remains a vast virgin territory for further investigation (symbolized by “p process (?)” in the table).

If the neutrino-heating mechanism (or some other mechanisms) fails to blow up massive stars, central PNSs collapse to BHs. Recent GR simulations by Ott et al. [199] show that the significant GW emission is associated at the moment of the BH formation, which can be a promising target of the advanced LIGO for a galactic source. As pointed out by Sumiyoshi et al. [359], the disappearing neutrino signals that mark the epoch of BH formations also can be a target of SK and IceCube [10]. When quarks and hyperons appear in the postbounce core, a neutrino burst produced by the sudden EOS softening and by the subsequent rebounce [123] is likely to be detected by IceCube for a galactic event [10]. The interval between core-bounce and the BH formation depends on the details of exotic physics in the super-dense PNS cores [360, 361]. All of these observational signatures should provide an important probe into the so-called dense QCD region in the QCD phase diagram (e.g., [362] for recent review) to which lattice calculations are hardly accessible at present. Concerning photons, no optical outbursts are expected in the BH-forming SNe (indicated by “no photon (?)” in Table 2), if not for rapid rotation and strong magnetic fields prior to core-collapse.

If a precollapse core rotates enough rapidly (typically initial rotation period less than 4 s) with strong magnetic fields (higher than ~1011 G), the MHD mechanism can produce bipolar explosions along the rotational axis predominantly driven by the field wrapping processes (Section 3(a)). If the MRI can be sufficiently resolved in global simulations, the MRI would exponentially amplify the initial magnetic fields to a dynamically important strength within several rotational periods (Sections 3(b) and 3(c)). The next-generation supercomputing resources are needed again to see the outcome. If such simulations would be executable, the MHD outflows will be produced even for more weakly magnetized cores than currently predicted. The GW signals in the MHD explosions are characterized by a burst-like bounce signal plus a secularly growing tail (Section 3.1, indicated by “burst and tail” in the table). Similar to the neutrino GWs in the neutrino-driven mechanism, a future detector (like FP-type DECIGO) is again necessary for detecting the low-frequency tail component. As discussed in Section 3.2.1, the pole to equator anisotropy of the shock propagation in the MHD explosions affect the neutrino signals through the MSW effects. The anisotropic shock passage to the H-resonance regions leads to a sudden decrease in the SK events (~2500 events for a galactic source), which is quite possible to identify by the SK-class detectors (indicated by “anisotropy in SK events” in Table 2). Since the MHD mechanism has such distinct signatures, a planned joint analysis of neutrino and GW data [363] would be potentially very powerful to tell the MHD mechanism from the other candidate mechanisms. If the resonant spin-flavor conversion (indicated as RSF in Table 2) occurs in the highly magnetized core, the neutronization burst of converts to that of , which is thus expected to be a probe into the magnetic moment of neutrinos (e.g., [334]). Finally the MHD explosions are considered to be a possible -process cite [364, 365] (see collective references in [247] for other plausible -process cites). This is because the mass ejection from the iron core occurs in much shorter timescales compared to the delayed neutrino-heating mechanism, so that the ejecta can stay neutron rich before it becomes proton-rich via neutrino capture reactions. Most of these possibilities regarding the MHD mechanism have been proposed so far in simulations with a crude treatment of neutrino transport, which is now rather easily reexamined by the state-of-the-art multi-D simulations (if the code is MHD).

Finally, what is about the story if the MHD mechanism fails? (indicated by "path to hypernovae in Table 2). The rapidly rotating PNS collapses into a BH, probably leading to the formation of accretion disk around the BH. Neutrinos emitted from the accretion disk heat matter in the polar funnel region to launch outflows or strong magnetic fields in the cores of the order of  G also play an active role both in driving the magneto-driven jets and in extracting a significant amount of energy from the BH (e.g., [366368] and see references therein). This picture, often referred as collapsar (e.g., [300, 301]), has been the working hypothesis as a central engine of long-duration GRBs for these 10 years (see references in [367371] for other candidate mechanism including magnetar models). However, it is still controversial whether the generation of the relativistic outflows predominantly proceeds via the neutrino-heating mechanism or the MHD mechanism. In contrast to a number of findings illustrated in Table 2, much little things are known about the BH-forming supernovae and also about the collapsar. This is mainly because the requirement for making a realistic numerical modeling is very computationally expensive, which (at least) necessitates the multi-D MHD simulations not only with GR for handling the BH formation (thus full GR), but also with the multiangle neutrino transfer for treating highly anisotropic neutrino radiation from the accretion disk [372]. To get a unified picture of massive stellar death, we need to draw a schematic picture (like Table 2) also in the case of the BH-forming supernovae. We anticipate that our long-lasting focus (and experience) on the explosion dynamics of canonical CCSNe will be readily applicable to clarifying the origin of the gigantic explosion energy of hypernovae and also unraveling (ultimately) the central engine of long-duration GRBs. This should provide yet another grand challenge in computational astrophysics in the next decades.

As repeatedly mentioned so far, all the numerical results in this paper should be tested by the next-generation calculations by which more sophistication is made not only in the treatment of multi-D radiation transport (both of neutrino and photon), but also in multi-D hydrodynamics including stellar rotation and magnetic fields in full GR. From an optimistic point of view, our understanding on every issue raised in this paper can progress in a step-by-step manner at the same pace as our available computational resources will be growing bigger and bigger from now on. Since 2009, several neutrino detectors form the Supernova Early Warning Systems (SNEWS) to broadcast the alert to astronomers to let them know the arrival of neutrinos [82]. Currently, Super-Kamiokande, LVD, Borexino, and IceCube contributes to the SNEWS, with a number of other neutrino and GW detectors planning to join in the near future. This is a very encouraging news towards the high-precision multimessenger astronomy. The interplay between the detailed numerical modeling, the advancing supercomputing resources, and the multimessenger astronomy, will remain a central issue for advancing our understanding of the theory of massive stellar core-collapse for the future. A documentary film recording our endeavors to make practicable the “multimessenger astronomy of CCSNe” seems not to show “fin” immediately and is becoming even longer as far as the three components evolve with time. To raise the edifice, it is becoming increasingly important to bring forward a worldwide, multidisciplinary collaboration among different research groups. We believe that such an approach would provide the shortest cut to get a deeper understanding of a number of unsettled and exciting issues that we were only able to touch in this paper.

Acknowledgments

K. Kotake is thankful to K. Sato for continuing encouragements. He also wishes to thank his collaborators, S. Yamada, M. Liebendörfer, N. Ohnishi, K. Sumiyoshi, K. Nakamura, T. Kuroda, M. Hashimoto, S. Harikae, N. Yasutake, N. Nishimura, K. Shaku, and H. Suzuki. The authors would like to acknowledge helpful exchanges and stimulating discussions with E. Müller, T. Foglizzo, H.-T. Janka, C. D. Ott, C. Fryer, J. Novak, S. Ando, P. Cerdá-Durán, M. Obergaulinger, J. Murphy, R. Fernández, E. O’Connor, M. Shibata, T. Font, and J. M. Ibañez. They are also grateful to our anonymous referees who gave us a number of valuable comments to enhance the quality of this paper. Numerical computations were carried out in part on XT4 and general common use computer system at the center for Computational Astrophysics, CfCA, the National Astronomical Observatory of Japan. This study was supported in part by the Grants-in-Aid for the Scientific Research from the Ministry of Education, Science and Culture of Japan (nos. 19540309, 20740150, 23540323, 23340069, and 24244036) and by HPCI Strategic Program of Japanese MEXT.