Multimessengers from core-collapse supernovae: multidimensionality as a key to bridge theory and observation

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.

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.

Introduction
The majority of stars more massive than ∼8 M ⊙ 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 year . 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 (Hirata et al. 1987;Bionta et al. 1987). 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., Sato & Suzuki (1987), see Raffelt (2012Raffelt ( , 2002 for a recent review). Recently significant progress has been made in the large waterCherenkov detectors such as Super-Kamiokande (Totsuka 1992) and IceCube (Hultqvist & IceCube collaboration 2011), and also in the liquid scintillator detector as KamLAND (Suzuki 1999). If a supernova occurs in our Galactic center (∼ 10 kpc), about 10 5ν e events are estimated to be detectable by IceCube (e.g., IceCube Collaboration et al. (2011) 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., Raffelt (2012) for a recent review). axis aligned with the rotation axis inferred from the ring (Wang et al. (2002), see however Kjaer et al. (2010) for a recent counter argument). 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., Wang et al. (2001); Leonard et al. (2006a,b)), a fact which has been used to argue that the central engine of the explosion is responsible for the non-sphericity (e.g., Kifonidis et al. (2003Kifonidis et al. ( , 2006; Scheck et al. (2006); Burrows et al. (2007c,a)). By performing a series of time-dependent, non-LTE(local thermodynamic equilibrium), radiationtransport simulations (e.g., Dessart & Hillier (2010 3 ), Dessart & Hillier (2011) recently pointed out that asymmetry in the ejecta can explain the increase in the continuum polarization observed at the nebular phases (Leonard et al. 2006c;Chornock et al. 2010). Dense knots, indications of ejecta clumpiness, and filaments seen in supernova remnants by HST in the visual (Fesen et al. 2006(Fesen et al. , 2011 and by ROSAT, Chandra, and XMM-Newton (Aschenbach et al. 1995;Hughes et al. 2000;Miceli et al. 2008) 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 Smartt (2009) for a review). Although observational measurements of the progenitor masses are currently not many and still highly uncertain (see  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 M ⊙ (Smartt 2009). A generic explosion energy of the SN II-P in the mass range is roughly on the order of 10 51 erg (Utrobin & Chugai 2011), 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 (Smartt 2009). For example, the inferred 56 Ni 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 (Woosley & Bloom 2006), all show much broader lines than SNe Ic. Due to the large kinetic energies of 2 − 5 × 10 52 erg, they have been referred to as "hypernovae" (e.g., Modjaz (2011) 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., Podsiadlowski et al. (2004Podsiadlowski et al. ( , 2010; Fryer et al. (2007); Smith et al. (2011), see however Yoon & Langer (2005a); Yoon et al. (2010)). In addition to the two branches mentioned above, Nomoto et al. (2006) 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 vs. 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 4 , secondly because we are still inaccessible to multi-D stellar evolution models, in which multi-D modeling has been a major undertaking (see Arnett & Meakin (2011a); Meakin & Arnett (2007) 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 selfinteraction that induces collective neutrino flavor oscillations predominantly in the vicinity of the neutrino sphere (see Raffelt & Smirnov 2007;Duan et al. 2008Duan et al. , 2010Dasgupta 2010, for reviews of the rapidly growing research field and collective references therein) and the Mikheyev-Smirnov-Wolfenstein (MSW) effect (Mikheev & Smirnov 1986) both in the stellar envelope and in the earth (see Raffelt (2012) for a review). Although the time profiles of neutrino signals can be potentially used like a tomography to monitor the envelope profile (e.g., Duan & Kneller (2009) 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., Janka et al. (2007) for a review), the acoustic mechanism (Burrows et al. 2007b), or the magnetohydrodynamic (MHD) mechanism (e.g., Kotake et al. (2006) see references therein). To pin down the true answer among the can-didate 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) (Antonioli et al. 2004)). In addition to optical observations using largest 8-10 m telescopes such as VLT and Subaru telescope (e.g., Patat et al. (2009); Tanaka et al. (2009)), the planned thirty-metertelescope (TMT) 5 with a refined spectropolarimetric technique is expected to detect more than 20 events of the SN polarization per year 6 . This should provide valuable information to understand the SN asymmetry with increasing statistics. Not only for understanding the origin of non-spherical ejecta morphology 7 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, Woosley & Heger (2007); Thielemann et al. (2011) for recent reviews) 8 .
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 review article, 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 Ando et al. (2012) Aso et al. (2008) 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". anism 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 article, 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 article, this one might go beyond such reviews by its new perspectives on the multimessenger astronomy.

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 proto-neutron 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 world-wide, multi-disciplinary collaboration to clarify the theory of massive stellar core-collapse and the formation mechanisms of compact objects.
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 non-uniform 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 Fryer et al. (2002);Fryer (2004a)). The neutrino heating proceeds dominantly via the charged current interactions (ν e + n ⇄ e − + p,ν e + p ⇄ e + + n) in the gain region. The neutrino heating rate in the gain region can be roughly expressed as Q + ∝ L ν µ ν −1 /r 2 where L ν 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 r and µ ν is the distance from the stellar center and the flux factor 9 , respectively (e.g., Janka (2001)). For example, L ν is treated as an input parameter in the so-called 'light-bulb" approach (e.g., Janka & Müller (1996)). This is one of the most prevailing approximations in recent 3D simulations (Iwakami et al. 2008(Iwakami et al. , 2009Wongwathanarat et al. 2010;Nordhaus et al. 2010) 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, L ν 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 multi-messengers, 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 abbreviation of "NH", "WHW", and "WW" means Nomoto & Mashimoto (1988), Woosley et al. (2002), and Woosley & Weaver (1995). 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 H.T. Janka and E.
Müller. "Princeton+" stands for the group chiefly consisting of the staffs in the Princeton University (A. Burrows, J.Murphy), Caltech (C.D.Ott), Hebrew University (E. Livne), Université de Provence (L. Dessart), and their collaborators. The SN group in the Basel university is led by M. Liebendörfer and F.K. Thielemann. "OakRidge+" stands for the SN group mainly consisting of the Florida Atlantic University (S. Bruenn) and the Oak Ridge National Laboratory (A. Mezzacappa, O.E.B. Messer) and their collaborators. Tokyo+ is the SN group chiefly consisting of the staffs in the National Astronomical Observatory of Japan (myself, T. Takiwaki), Kyoto university (Y. Suwa), Numazu College (K. Sumiyoshi), Waseda University (S. 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. t exp in the fifth column indicates an approximate typical timescale when the explosion initiates and E exp represents the explosion energy normalized by Bethe(=10 51 erg) given at the postbounce time of t pb , both of which are attempted to be sought in literatures 10 . 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 O(v/c). The definition of the "RBR", "IDSA", and "MGFLD" will be given soon in the following.
Due to the page limit of this article, we have to start the story only after 2006 (see, e.g., Janka et al. (2007) for a complete review, and also Cardall (2005) 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 M ⊙ star (Kitaura et al. 2006;Janka et al. 2008) 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 M ⊙ progenitors of WHW02 triggered by the hypothesised first-order QCD phase transition in the protoneutron star (PNS) (Sagert et al. 2009). 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 non-rotating 11.2 M ⊙ progenitor of WHW02 (Buras et al. 2006a), and then for a 15M ⊙ progenitor ) of WW95 with a relatively rapid rotation imposed 11 . They newly brought in the so-called "ray-by-ray" approach (indicated by "RBR" in the table), in 10 but if we cannot find them, we remain them as blank "-". Table 1: Selected lists of recent milestones reported by SN groups around the world ("Group"), which obtained explosions by the neutrino-heating mechanism (indicated by "ν-driven") or the acoustic mechanism ("Acoustic") (See text for more details).
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 (Rampp & Janka 2000) 12 . 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 (∼ 10 51 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 Buras et al. (2006a). Very recently, Mueller et al. (2012) reported that more energetic explosions are obtained for the 11.2 and 15M ⊙ stars in their GR 2D simulations compared to the corresponding post-Newtonian models (e.g., ). 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., Dimmelmeier et al. (2002b); Cordero-Carrión et al. (2009), indicated by "C-GR" in the table).
It is rather only recently that fully 2D multi-angle Boltzmann transport simulations become practicable by the Princeton+ group Brandt et al. 2011). 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. (2011) succeeded in following the dynamics until ∼ 400 ms after bounce for a non-rotating and a rapidly rotating 20 M ⊙ 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 Multi-Group Flux-Limited Diffusion scheme which eliminates the angular dependence of the neutrino distribution function (see, e.g., Bruenn (1985) 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., Bruenn et al. (2010)), Bruenn et al. (2009) obtained neutrino-driven explosions for non-rotating progenitors in a relatively wide range in 12, 15, 20, 25 M ⊙ of WHW02 (see table). These models tend to start exploding at around 300 m after bounce, and the explosion energy for the longest running model of the 25 M ⊙ progenitor is reaching to 1 B at 1.2 s after bounce ).
On the other hand, the 2D MGFLD simulations implemented in the VULCAN code (Burrows et al. 2006) obtained explosions for a variety of progenitors of 11, 11.2, 13, 15, 20, and 25 M ⊙ not by the neutrino-heating mechanism but by the acoustic mechanism 13 . 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. (2007) (Burrows et al. 2007c) 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 ), the Tokyo+ group reported explosions for a non-rotating and rapidly rotating 13 M ⊙ progenitor of NH88. They pointed out that a stronger explosion is obtained for the rotating model comparing to the corresponding non-rotating 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 Liebendörfer et al. (2009) 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 non-local 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 (ν x , i.e., ν µ , ν τ and their anti-particles) as well as the energy-coupling weak interactions have yet to be implemented. Extending the 2D modules in (Suwa et al. 2010) to 3D, they recently reported explosions in the 3D models for an 11.2 M ⊙ progenitor of WHW02 . 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 (Liebendörfer et al. 2005)).
In the next section, we briefly summarize the findings obtained in our 2D (Suwa et al. 2010) and 3D  simulations, paying particular attention to how multidimensionality such as SASI, convection, and rotation could affect the neutrino-driven explosions.  (Nomoto & Mashimoto 1988). Thick lines in red (for 2D) and black (for 1D) show the position of shock waves, noting for 2D that the maximum (top) and average (bottom) shock position are shown. The red dashed line represents the position of the gain radius, which is similar to the 1D case (not shown). Figure 1 depicts the difference between the time evolutions of 1D (thin gray lines) or 2D (thin orange lines) simulation of the 13 M ⊙ 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     . The explosion energy is defined as the total energy (internal plus kinetic plus gravitational), integrated over all matter where the sum of the corresponding specific energies is positive. Models with "-hr" indicates the ones with higher numerical resolution, in which the mesh numbers in the lateral direction are doubled. Middle and right panels are snapshots of the density (left half) and the entropy (right half) for 1D (left) and 2D rotating (right) model at the epoch when the shock reaches 1000 km. These figures are taken from Suwa et al. (2010). 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 τ adv /τ heat , which is the ratio of the advection to the neutrino heating timescale 14 . For the 2D model (right panel), it can be shown that the condition of τ adv /τ heat 1 is satisfied behind the aspherical shock, which is deformed predominantly by the SASI (ℓ = 2 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 (Foglizzo & Tagger 2000;Foglizzo 2001), 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.

Multidimensionality in multi-D radiation hydrodynamic simulations
To see clearly the effects of stellar rotation on the postbounce evolution, a rapidly rotating model was taken in Suwa et al. (2010), in which a constant angular frequency of Ω 0 = 2 rad/s is imposed inside the iron core with a dipolar cut off (∝ r −2 ) outside. This corresponds to β ∼ 0.18% 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 ℓ = 2 mode for the rotating model (e.g., Figure 3, although the ℓ = 1 mode can be as large as the ℓ = 2 mode when the SASI enters the non-linear regime ( 200 ms after bounce). In contrast to this rotation-induced ℓ = 2 deformation, the ℓ = 1 mode tends to be larger than the ℓ = 2 mode for the 2D models without rotation in the saturation phase.
The left panel of Figure 4 shows the comparison of the explosion energies for the 2D models with and without rotation (models with "-rot" or "-2D", respectively). 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 180 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 north-south symmetric (ℓ = 2) explosion can expel more material than for the unipolar explosion (compare the middle and right panels in Figure 4). The explosion energies when we terminated the simulation are less than 10 50 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 Burrows et al. (2007a), which we are going to study as a follow-up 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 M ⊙ star of Woosley et al. (2002) with spectral neutrino transport that is solved by the IDSA scheme in a ray-by-ray manner . We briefly summarize the results in the following.   , visualized by mass shell trajectories in thin gray lines (left panel). Thick red lines show the position of shock waves, noting that the maximum (top), average (middle), and the minimum (bottom) shock position are shown, respectively. The green line represents the shock position of the 1D model. "1.30" and "1.40" indicates the mass in unit of M ⊙ enclosed inside the mass-shell. Right panel shows the evolution of average shock radii for the 2D (green line), and 3D (red line) models. The "3D low" (pink line) corresponds to the low resolution 3D model, in which the mesh numbers are taken to be half of the standard model (see text). These figures are taken from Takiwaki et al. (2012).
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. (2006a). The right panel of Figure  5 shows a comparison of the average shock radius vs. 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. (2006a) (see their Figure 15 for model s112 128 f) 15 .
Comparing the shock evolution between our 2D (green line in the right panel of Figure  5) and 3D model (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 (n φ =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. (2010) 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 16 . In the 3D model (left panel), non-axisymmetric 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 non-axisymmetric 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 Takiwaki et al. (2012) 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 netheating 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 Takiwaki et al. (2012) was that the 15 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. 16 Note that the polar axis is tilted (about π/4) both in the left and middle panel.
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 17 . 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 is 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 (Misner & Sharp 1964), the consideration of GR was standard in the pioneering era of supernova simulations (e.g., May & White (1966)). One year after Colgate & White (1966), Schwartz (1967) 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 18 . Using GR Boltzmann equations de-17 It is of crucial importance to conduct a convergence test in which a numerical gridding is changed in a systematic way (e.g. Hanke et al. (2011)).
18 Citing from his paper, "In this calculation, the neutrino luminosity of the core is found to be 10 54 erg/s, or 1/2 a solar rest mass per second !! .... This is the mechanism which the supernova explodes". The  rived by Lindquist (1966), Wilson (1971) developed a 1D GR-radiation-hydrodynamic code including a more realistic (at the time) description of the collisional term than the one in Schwartz (1967). 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 (van Riper 1979;van Riper & Lattimer 1981;van Riper 1982). 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 Bruenn et al. (2001) for a complete list of references for the early GR studies). In the middle of the 1980s, Bruenn (1985) developed a code that coupled 1D GR hydrodynamics to the MGFLD transport of order (v/c) 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 Sumiyoshi-Yamada et al. (Yamada 1997;Yamada et al. 1999;Sumiyoshi et al. 2005bSumiyoshi et al. , 2007 19 and by Liebendörfer-Mezzacappa-Bruenn et al. (Mezzacappa & Matzner 1989;Bruenn et al. 2001;Liebendörfer et al. , 2004 (and by their collaborators).
Among them, Bruenn et al. (2001) 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. (2011) reported the update of Bruenn et al. (2001) 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 Lentz et al. (2011)).
Among the most up-to-date multi-D models with spectral neutrino transport mentioned earlier (Table 1), the GR effects are at best attempted to be modeled by using a modified gravitational potential that takes into account a 1D, post-Newtonian correction (Buras et al. 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. 2006c,a; Bruenn et al. 2010). 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 (Shibata & Sekiguchi 2005a) or 3D (Shibata & Sekiguchi 2005b;Ott et al. 2007d), and references therein) or using the conformally-flatness approximation (CFC) (e.g., Dimmelmeier et al. (2002b); Cordero-Carrión et al. (2009)). Although extensive attempts have been made to include microphysics such as by the Y e formula (Liebendörfer 2005a) or by the neutrino leakage scheme (Sekiguchi 2010), 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, Mueller et al. (2012); Kuroda et al. (2012)).
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 multi-messengers 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 article 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 20 . 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.

Gravitational waves
The paper by Müller (1982) 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 Ott (2009b) 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) 21 . 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 timedependent 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., Müller (1982); Mönchmeyer et al. (1991); Yamada & Sato (1995); Zwerger & Müller (1997); Kotake et al. (2003bKotake et al. ( , 2004b; Shibata & Sekiguchi (2004); Ott et al. (2004Ott et al. ( , 2007a; Dimmelmeier et al. (2002aDimmelmeier et al. ( , 2007Dimmelmeier et al. ( , 2008; Scheidegger et al. (2008Scheidegger et al. ( , 2010b), and references therein).
As summarized by Ott (2009b), 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 (Dimmelmeier et al. 2002a(Dimmelmeier et al. , 2008 and also by fully GR 3D simulations (Ott et al. 2004) both including realistic EOSs and a deleptonization effect based on 1D-Boltzmann simulations (Liebendörfer 2005b).
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 (Maeder & Meynet 2000;Heger et al. 2005;Ott et al. 2006b). To explain the observed rotation periods of radio pulsars, the precollapse rotation periods are estimated to be larger than ∼ 100 sec (Ott et al. 2006b). 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., Kotake et al. (2004b)).
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., Kotake (2011) 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 pre-collapse density inhomogeneities (Burrows & Hayes 1996;Müler & Janka 1997;Fryer 2004b,a), moderate rotation of the iron core ), nonaxisymmetric rotational instabilities (Rampp et al. 1998;Ott et al. 2007b), gmodes (Ott et al. 2006a) and r-modes pulsations (Andersson et al. 2010) of PNSs, and the SASI (Kotake et al. ( , 2009a; Murphy et al. (2009)) 22 . Among them, the most promising GW sources may be convection and SASI, because the degree of the initial inhomogeneities (Arnett & Meakin 2011b) 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 lightbulb scheme (which we shortly call as parametric SASI simulations hereafter), we pointed out that the GW amplitudes from anisotropic neutrino emission 23 increase almost monotonically with time, and that such signals may be visible to next-generation detectors for a Galactic source (Kotake et al. , 2009a. By performing such a parametric simulation but without the excision inside the PNS, Murphy et al. (2009) 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. (2010) 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).  analyzed the GW emission based on their long-term 2D ray-by-ray Boltzmann simulations, which seem very close to produce explosions ) (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 Kotake et al. (2007Kotake et al. ( , 2009a; Yakunin et al. (2010)), 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 signalto-noise ratio.
Regarding the GW predictions in 3D models, Müler & Janka (1997) coined the first study to analyze the GW signature of 3D non-radial 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 60 • . 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 22 see also Ott et al. (2011) for the GW signals at the black hole formation. 23 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 (1978); Turner & Wagoner (1979) (see recent progress in Favata (2010)). It is expected as a primary GW source also in gamma-ray bursts (Hiramatsu et al. 2005;Suwa & Murase 2009) and Pop III stars (Suwa et al. 2007a). models, and the wave amplitudes from neutrinos are a factor of 10 larger than those due to non-radial matter motions. With another pioneering (2D) study by Burrows & Hayes (1996) Burrows & Hayes (1996), 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 (Fryer et al. 2002, 2004 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 non-spherical 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 (Kotake et al. 2009a(Kotake et al. ,b, 2011. 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. Figure 9 shows the evolution of 3D hydrodynamic features from the onset of the nonlinear regime of SASI (top left) until the shock break-out 24 (bottom right) with the gravitational waveform from neutrinos inserted in each panel. After about 100 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  Kotake et al. (2009b)). The second and fourth quadrant of each panel shows the surface of the standing shock wave. In the first and third quadrant, the profiles of the high entropy bubbles (colored by red) inside the section cut by the ZX plane are shown. The side length of each plot is 1000km. The insets show the gravitational waveforms from anisotropic neutrino emissions, with '+' on each curves representing the time of the snapshot. Note that the colors of the curves are taken to be the same as the top panel of Figure 10. This figure is taken from Kotake et al. (2009b).   (bottom) and from the sum of neutrinos and matter motions (top), seen from the polar axis and along the equator (indicated by 'Pole' and 'Equator') with polarization (+ or × modes) for two representative 3D models of A and B (see Kotake et al. (2009b) for details). The distance to the SN is assumed to be 10 kpc.

Stochastic Nature of Gravitational Waves in 3D simulations
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 quadrant 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 (Kotake et al. , 2009a. 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 becomes 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 0.5% 25 . 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 Müller et al. (2011)).
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 neutrinoheating 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 (Ott et al. 2006a) and the non-axisymmetric instabilities (Ott et al. 2007b;Scheidegger et al. 2010b) of the PNSs, and the enhanced neutrino emissions inside the PNSs ). Bearing these caveats in mind, a piece of encouraging news is that the gravitational waveforms obtained in the 2D radiation-hydrodynamic simulations (Yakunin et al. 2010) are similar to the ones obtained in our 2D study using the light-bulb scheme (Kotake et al. 2009a). More recently, Müller et al. (2011) confirmed the stochastic nature by analyzing the GW features obtained in their 3D models (Wongwathanarat et al. 2010). 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 Murphy et al. 2009;Yakunin et al. 2010 2011)). The time is measured from the epoch when the neutrino luminosity is injected from the surface of the neutrino sphere. For this 3D model with rotation, the rotational flow is imposed to advect to the PNS surface at around t = 400 ms. The supernova is assumed to be located at the distance of 10 kpc.

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 . In 3D, the modes of SASI are divided into sloshing modes and spiral modes (e.g., Iwakami et al. (2008)). Asymmetric m = 0 modes so far studied in 2D models and axisymmetric m = 0 modes are classified into the sloshing modes, where m stands for the azimuthal index of the spherical harmonics Y m l . In the latter situation, the ±m modes degenerate so that the +m modes has the same amplitudes as the −m 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 (Blondin & Mezzacappa 2007;Iwakami et al. 2009). In this situation, the +m modes has the different amplitudes from −m modes. Such rotating non-axisymmetric m = 0 modes are called as spiral modes. These non-axisymmetric 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 , and also because they may have a potential to generate pulsar spins (Blondin & Mezzacappa (2007); Rantsiou et al. (2010), see also Yamasaki & Foglizzo (2008); Wongwathanarat et al. (2010)). These finding 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 Iwakami et al. (2009), whose specific angular momentum is assumed to agree with recent stellar evolution models (Maeder & Meynet 2000;Heger et al. 2005). 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 Kotake et al. (2011) 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 non-linear 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 (z-axis) at this epoch. It should be emphasized that although the dominance of h equ ν,+ observed in the current 3D simulations is similar to the one obtained in previous 2D studies (Kotake et al. 2009a), 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 patters behind the standing shock simultaneously change in every direction like the non-rotating 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 equation (9) in Kotake et al. (2009a)) is positive near the north and south polar caps, the dominance  Fig. 13.-Spectral distributions of GWs from matter motions ("Matter") and neutrino emission ("Neutrino") seen from the pole or the equator for a representative 3D rotating model (e.g., Kotake et al. (2011)) with the expected detection limits of TAMA300 (Ando & the TAMA Collaboration 2005), first LIGO and advanced LIGO (Thorne 1995), Large-scale Cryogenic Gravitational wave Telescope (LCGT) (Kuroda & LCGT Collaboration 2010) and Fabry-Perot type DECIGO (Kawamura 2006). It is noted that h char is the characteristic gravitational wave strain defined in Flanagan & Hughes (1998). The distance to the supernova is assumed to be 10 kpc. Note that for the matter signal, the + mode seen from the polar direction is plotted (from Kotake et al. (2011)).
of the polar neutrino luminosities leads to make the positively growing feature of h equ ν,+ in Figure 11 (blue line). From the spectral analysis of the gravitational waveform (Figure 12), it can be readily seen that it is not easy to detect these neutrino-originated GW signatures with slower temporal evolution ( O(10)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 (Kawamura (2006), black line in Figure 12). 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., ) and the acoustic mechanism (Ott et al. 2006a). Therefore the detection of such signals could be expected to provide an important probe into the explosion mechanism (e.g., Ott (2009b,a)).
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  or by the gray transport scheme (Müller et al. 2011). 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 cut-out, efficient GW emission of the oscillating neutron star (Ott et al. 2006a) and non-axisymmetric instabilities (Scheidegger et al. 2008(Scheidegger et al. , 2010b of the PNSs, and the enhanced neutrino emissions inside the PNSs ) 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.

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 26 . 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., Woosley et al. (2002); Woosley & Weaver (1995); Hashimoto (1995), and collective references in Jerkstrand (2011)). The observational determination of the masses of the three main radioactive isotopes 56 Ni, 57 Ni, and 44 Ti 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., Woosley & Hoffman (1991)). During the shock propagation, iron group elements such as 56 Ni and its daughter nucleus 56 Co are predominantly produced, which are radioactive with a life-time 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, 56 Co 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 56 Co nuclide. To explain the bolometric light curve of SN1987A in such a phase, the 56 Ni mass was determined to be 0.07M ⊙ (Bouchet et al. 1991).
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 (Fransson & Kozma 1993). 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., Fransson & Kozma (1993); Mazzali et al. (2001); Maurer et al. (2011);Jerkstrand et al. (2011) and references therein). Such a time-dependent modeling by Fransson & Kozma (1993) revealed the 57 Ni mass of ∼ 3.3×10 −3 M ⊙ of SN1987A, which agrees well with observations (e.g., Varani et al. (1990); Kurfess et al. (1992) and collective references in Filippenko (1997)).
By the similar reason to 57 Ni just mentioned above, the determination of the 44 Ti is also complicated. The most recent study by Jerkstrand et al. (2011) gives an estimate of the 44 Ti to be 1.5 +0.5 −0.5 × 10 −4 M ⊙ , which is in good agreement with the eight-year spectrum analysis of SN1987A (e.g., Chugai et al. (1997), see also Lundqvist et al. (2001)). As shown above, the amount of 44 Ti is typically one order-of-magnitude smaller than that of 57 Ni, however it is crucially important for young supernova remnants due to its long life-time (∼ 86 years). It is worth mentioning that NASA will launch the satellite NuSTAR (Nuclear Spectroscopic Telescope Array) to study 44 Ti production in CCSNe. The detector will be able to map out the 44 Ti distribution of the supernova remnant Cassiopeia A and can get velocity distributions of the 44 Ti in SN 1987A. By comparing detailed modeling of the SN nucleosynthesis in the context of 2D and 3D models (e.g., Fryer et al. (2006); Magkotsios et al. (2010)), these are expected to both provide direct probes of the explosion asymmetry.
Ever since SN1987A, challenges to the classical spherical modeling (Hashimoto 1995;Woosley & Weaver 1995;Thielemann et al. 1996;Rauscher et al. 2002) have been building 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 44 Ti synthesis along the poles as indicated in SN1987A (e.g., Nagataki et al. (1997)). More recently, 3D effects have been more elaborately studied (Hungerford et al. (2005), see also Young et al. (2006)) as well as the impacts of different explosions by employing a number of progenitors (Joggerst et al. 2010) or by assuming a jet-like explosion (Couch et al. 2009;Tominaga 2009), which is one of the possible candidates of hypernovae (e.g., Maeda & Nomoto (2003)).
In addition to the above-mentioned work, nucleosynthesis in a more realistic simulation that models multi-D neutrino-driven explosions has been also extensively studied (e.g., Kifonidis et al. (2003Kifonidis et al. ( , 2006; Gawryszczak et al. (2010) and references therein). Although a small network has ever been included in the computations, these 2D simulations employing a light-bulb scheme (Kifonidis et al. 2003) or a more accurate gray transport scheme Kifonidis et al. 2006) 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. (2006) 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. (2010) 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(p-) nuclei including 92,94 Mo and 96,98 Ru (e.g., Pruet et al. (2006);Fröhlich et al. (2006); Wanajo & Janka (2011) 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 ν e orν e , which is decisive for quantifying the νp processes. Therefore the two approaches, namely light-bulb scheme vs. 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., Fujimoto et al. (2011) for more details). By changing neutrino luminosities from PNSs systematically (L ν ), we discuss how the multidimensionality formed by the SASI and convection could impact on the explosive nucleosynthesis. We employ a non-rotating 15M ⊙ 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 (Cyburt et al. 2010). 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 review article which shows only a snapshot of the moving theoretical terrain. Figure 14 shows how explosion energies (top panel), explosion timescales as well as the PNS masses (bottom panel) change with the input neutrino luminosity (L νe ) for the 2D parametric explosion models (e.g., Fujimoto et al. (2011) for more details). As seen, higher L νe makes explosion energy greater, the onset of explosion (t exp ) earlier, the PNS masses (M PNS ) smaller. The PNS masses range in (1.54 − 1.62)M ⊙ for models with higher luminosity (L νe 4.25 × 10 52 erg s −1 , panels (b), blue line), which are comparable to the baryonic mass (around 1.5M ⊙ ) of neutron stars observed in binaries (Schwab et al. 2010). The explosion energies for these high luminosity models (panel (a), red line) are also close to the canonical supernova kinetic energy of 10 51 erg. To discuss explosive nucleosynthesis in the following, we thus choose to focus on the model with L νe = 4.5 × 10 52 erg s −1 as a reference.

Symmetry breaking in explosive nucleosynthesis
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 (θ < π/4, θ > 3π/4), where both Si and O burning proceed to produce abundantly 56 Ni and 28 Si. These aspherical  Uchida et al. (2009) 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 − 15M ⊙ . 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 M ⊙ progenitor have a close correlation with the ones in the Cygnus loop (see Fujimoto et al. (2011) 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 56 Ni and 44 Ti, they are ∼ 0.06M ⊙ and 4×10 −5 M ⊙ , respectively for the model (left panel in Figure 16). For SN1987A, masses of 56 Ni and 44 Ti are deduced to be ∼ 0.07M ⊙ (Shigeyama et al. 1988;Woosley 1988) and 1 − 2 × 10 −4 M ⊙ (Nagataki 2000, and references therein). The obtained mass of 44 Ti is, therefore, not enough for SN1987A as well as for Cas A (1.6 +0.6 −0.3 × 10 −4 M ⊙ ) (Renaud et al. 2006). The shortage of 44 Ti is a long-standing problem, which might be solved by 3D effects . For the abundance patterns, our reference model (top panel in Figure 16) is shown to reproduce a similar trend to the solar system (right panel, top), which is in sharp contrast to the model with lower neutrino luminosity (right panel, bottom).
Finally we note that there should are at least two barriers to get over in the currentgeneration 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 ∼ 10 12 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., Gawryszczak et al. (2010)). 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. A encouraging news is that several groups Ti). Right panel shows element abundance yields (X i /X ⊙,i , normalized to solar) vs. mass number for our reference model of L νe = 4.5 × 10 52 erg s −1 and for (b) L νe = 3.9 × 10 52 erg s −1 Thick horizontal-dashed lines represent a factor equals to that of 16 O, while two normal and two thin lines denote a factor equals to that of 16 O times 2, 1/2, 3, and 1/3, respectively. These figures are taken from Fujimoto et al. (2011).
are pursuing it with the use of advanced numerical techniques such as the adaptive-meshrefinement approach  or Ying-Yang gridding (Wongwathanarat et al. 2010). 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.

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.0.2 and discuss properties of the so-called magnetorotational instability ( (Wang et al. 2001(Wang et al. , 2002, and the interpretation of magnetars (Duncan & Thompson 1992;Lattimer & Prakash 2007), collapsars and their relevance to gamma-ray bursts ; Woosley & Heger (2006); Yoon & Langer (2005b)) 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 (Meier et al. 1976). 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 Balbus & Hawley (1998)). Each of them we are going to review briefly in the following.  2009)). The outer edge of the sphere colored by blue represents the radius of 7.5 × 10 7 cm. These panels highlight not only the wound up magnetic field around the rotational axis (left), but also the fallback of the matter from the jet head of the downwards to the equator, making a cocoon-like structure behind the jet (right).   Figure 17) along the polar axis (left panel) or the equatorial plane (middle panel) at 2 ms after the revival of the stalled bounce shock. The matter pressure is shown by green line as a reference. The right panel shows a velocity profile along the pole measured from the shock-stall (0.0 ms in the panel). For the equator, the magnetic pressure is much less than the ram pressure (middle panel), while the magnetic pressure amplified by the field wrapping along the pole becomes as high as the ram pressure of the infalling material at the shock front (left panel, note that the shock position can be inferred by the discontinuity at around 150 km), leading to the MHD-driven shock formation (right panel). These figures are taken from .

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 ), the precollapse magnetic field is set to be as high as 10 12 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 (B φ ) can be expressed as (e.g., Meier et al. (1976)), where X denotes the distance from the rotational axis and B X represents the X component of the poloidal fields in the cylindrical coordinates (X, Z). Given that a precollapse iron core has strong magnetic field (such as B X,0 ∼ 10 12 G with B X,0 representing the initial poloidal fields) and rapid rotation (such as P init 4 s with P init being the precollapse rotation period), the typical amplification of B φ near core bounce may be estimated as where t represents the typical timescale when the field amplification via the wrapping processes saturates, P rot denotes the rotational period near bounce. Equation (2) shows that the toroidal fields grow linearly with time by wrapping the poloidal fields (B X ), which is the so-called Ω dynamo (e.g., Cerdá-Durán et al. (2007)). The 10 15 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 ∼ 10 15 G. 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, P = 4 × 10 28 ρ 10 10 g/cm 3 ∆v 2 × 10 9 cm/s 2 erg/cm 3 , where ρ and ∆v denotes typical density and radial velocity in the vicinity of the stalled shock. When the toroidal magnetic fields are amplified as large as ∼ 10 15 G due to the field wrapping behind the shock, the resulting magnetic pressure, B 2 8π ∼ (10 15 ) 2 /8π ∼ 10 29 erg/cm 3 , 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., Burrows et al. (2007a)). 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.4c), with c 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 10 51 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., Burrows et al. (2007a); Dessart et al. (2008)) 27 .

On the MRI
We are now in a position to discuss possible impacts of the MRI in supernova cores. Akiyama et al. (2003) was 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., Equation (2)). After the MRI enters to the saturated state, the field strength might reach ∼ 10 15−16 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 Hawley et al. (1995); Balbus & Hawley (1998); Masada et al. (2006)). 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. (2005) suggested that the additional energy input by the viscous heating can help the neutrinodriven 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 ( 10 12 G), are predominantly driven by the field wrapping processes, not by the MRI. Shibata et al. (2006) 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 Etienne et al. (2006)). As well-known, the growth rate of MRI-unstable modes depends on the product of the field strength and the wave number of the mode (Balbus & Hawley 1998). When a precollapse rapidly core is strongly magnetized as 10 11 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 (Shibata et al. 2006). On the other hand, the fastest growing modes drop below several meters (Obergaulinger et al. 2009) for canonical supernova progenitors, which rotates much more slowly with weaker initial fields (∼ 10 9 G) as predicted by recent stellar evolution models (Maeder & Meynet 2000;Heger et al. 2005). At present, it is computationally too expensive to resolve those small scales in the global MHD simulations, typically more than two or three orders-ofmagnitudes 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 Balbus & Hawley (1998)). Obergaulinger et al. (2009) were the first to study the growth and saturation level of the MRI in the supernova environment 28 . 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 semi-global simulation systematically in 2D and 3D , they derived scaling laws for the termination of the MRI. As estimated in Akiyama et al. (2003), the MRI was shown to amplify the seed fields exceeding 10 15 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 non-linear 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 29 to the thermal energy occurs via the MRI turbulence, which was treated parametrically in Thompson et al. (2005). 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. (1998). The hydrodynamic part is based on the second-order Godunov scheme (van Leer 1979), 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 (Clarke 1996). 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 (Evans & Hawley 1988;Stone & Norman 1992).
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., Pessah & Goodman (2009);Longaretti & Lesur (2010)). 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 30 and discuss their possible impacts on the explosion dynamics.
29 (that the differential rotation taps) 30 which is not always familiar with supernova modellers, The vertical axis is normalized by the initial magnetic energy E mag,0 ≡ B 2 0 /8π, that is E mag ≡ E mag /E mag,0 , M rφ ≡ M rφ /E mag,0 , and R rφ ≡ R rφ /E mag,0 where single bracket denotes the volume average of the physical variables. Horizontal axis is normalized by the rotational period t rot ≡ 2π/Ω. Three typical evolutionary stages, (i) exponential growth stage, (ii) transition stage, and (iii) nonlinear turbulent stage, are denoted by white, dark gray and light gray shaded region. (b). The amplitude of power spectra of the magnetic energy |B(k)| 2 along the k x and k z axes. The dashed and dotted curves describes the k x and k z components respectively. The vertical axis is normalized by the initial magnetic energy E mag,0 . The wave-numbers shown in the horizontal axis are normalized by the 2π/L. The upper thick curve demonstrates the Kolmogorov slope for isotropic incompressible turbulence ∝ k −11/3 . The lower slope is proportional to k −6 just for the reference. In our local simulations, we prepare a 2D shearing box which is threaded by the initial vertical field B 0 = 1.4 × 10 12 G with imposing initial gas pressure P 0 = 4 × 10 28 dyn cm −2 , angular velocity Ω 0 = 2000 sec −1 and shear parameter q 0 = 1.25, respectively. Here the shear parameter represents the degree of differential rotation as q ≡ −d ln Ω/d ln r with r being the radial coordinate. The other parameters are fixed to mimic the properties in the vicinity of the PNS's surface in the postbounce phase (i.e. the density is taken to be ρ 0 = 10 12 g cm −3 and the box size is taken to be L x = 4 km (horizontal) × L z = 1 km (vertical)).  c) and (d) correspond to the time slice t = 3t rot , 4t rot , 6t rot and 20t rot respectively. The color bar indicates the size of the toroidal magnetic field normalized by the strength of the initial vertical magnetic field (B 0 = 1.4 × 10 12 G). Note that the vertical and horizontal are normalized by the typical spatial scale of the local portion in the proto-neutron stars which is selected as L = 1km in the fiducial run. Three typical evolutionary stages are observed, (i) exponential growth stage, (ii) transition stage, and (iii) nonlinear turbulent stage (compare with Figure 19). The initial field strength is varied from ∼ 10 11 to ∼ 6 × 10 12 G. The dashed lines denote reference slopes proportional to (a) E 1/2 mag and (b) E mag respectively. We take the time-average of volume-averaged quantities during 90t rot < t < 100t rot .

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 (Hawley & Balbus 1992;Hawley et al. 1995). 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 a exact solution for the nonlinear MHD equations (Goodman & Xu 1994). 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 (Goodman & Xu 1994;Obergaulinger et al. 2009). 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 19a shows that, at this stage, the turbulent Maxwell stress (M rφ = −B r B φ /4π) dominates over the Reynolds stress (R rφ = ρv r δv φ ). 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 19b shows the power spectra of the magnetic energy in the nonlinear turbulent stage. The logarithmic slope of the spectrum roughly coincides with −11/3 in the regime 10 k 100 while it steepens at higher wave-numbers where numerical effects become important (like k −6 ). This indicates that the turbulent state in the MRI can be represented by the Kolmogorov spectrum for a homogeneous, incompressible turbulence that scales as k −11/3 . These features would be consistent with the 3D properties of the MRI-driven turbulence (Hawley et al. 1995). 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 (Sano & Inutsuka 2001;Obergaulinger et al. 2009). We however adopted the simulation box with a large aspect ratio (L x /L z = 4) and quasi-isothermal equation of state (γ ≃ 1), both of which could facilitate to disrupt the channel structures (Bodo et al. 2008) 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 E mag,ini ≡ B 2 ini /8π. 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 (Hawley et al. 1995;Sano et al. 1998Sano et al. , 2004Pessah & Chan 2008) 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 E mag ∝ E 1/2 mag,init ∝ B ini (Figure 21(b))). This is also qualitatively consistent with Hawley et al. (1995). In this way, we investigated parameter dependence of the rotational shear (q parameter) as well as the initial pressure (P ) on the saturated value of the Maxwell stress and deduced a scaling relation that yields M rφ ∝ B 2 P 1/4 [q/(2−q)] 1/4 Ω 1/2 . Here the shear parameter represents the degree of differential rotation as q ≡ −d ln Ω/d ln r with r 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 f 30 represents a typical amplification factor of the stress normalized by 30 (e.g., Figure  21), q pb,1 is the shear rate normalized by unity and s pb,1 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 B 0 = 1.4 × 10 12 G, P 0 = 4 × 10 28 erg cm −3 , Ω 0 = 2000 rad s −1 , and q = 1.25 are replaced by B pb , P pb , Ω pb and q pb with some amplification factors.
A typical volume in which the MRI-driven turbulence is active, may be estimated as V MRI = 4πR 2 h ≃ 10 21 cm 3 R 2 7 h 6 , where R 7 = R/10 7 cm is typical radius of the PNS normalized by 10 7 cm, and h 6 = h/10 6 cm is the typical radial width of the MRI-active layer normalized by 10 6 cm. Then the energy releasing rate L MRI ≡ ǫ MRI V MRI becomes, Typical timescales of neutrino-driven explosions observed in recent supernova simulations are 400 ms after bounce (e.g., Suwa et al. (2010)). So the energy deposition could be as high as 10 50 erg 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 (Pessah & Goodman 2009;Longaretti & Lesur 2010). To obtain more accurate estimate of the amplification factor of f (equation (6)), 3D simulations are unquestionably required, which we are currently undertaking.

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., Woosley & Bloom (2006)). This can be really the case, albeit minor, for progenitors of rapidly rotating metalpoor stars, which experience the so-called chemically homogeneous evolution (Woosley & Heger 2006;Yoon & Langer 2005b). 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., MacFadyen & Woosley (1999);MacFadyen et al. (2001), see Dessart et al. (2008); Harikae et al. (2009Harikae et al. ( , 2010 for collective references), which are presumably linked to the formation of long-duration gamma-ray bursts (GRBs) (e.g., Meszaros (2006) 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 (Yamada & Sawai 2004;Kotake et al. 2004b;Obergaulinger et al. 2006b;Cerdá-Durán et al. 2007;Shibata et al. 2006;Scheidegger et al. 2010a). 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 B 0 10 12 G (Obergaulinger et al. 2006b;Kotake et al. 2004b). Considering that the typical magnetic-field strength of GRB progenitors is at most ∼ 10 11−12 G (Woosley & Heger 2006), this is already an extreme situation. In a more extremely case of B 0 ∼ 10 13 G, a secularly growing feature in the waveforms was observed (Obergaulinger et al. 2006a;Shibata et al. 2006;Scheidegger et al. 2010a). In the following, we summarize the GW signatures based on our 2D SRMHD simulations ) 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 (B 0 = 10 12 G) also with rapid rotation initially imposed (β parameter = 0.1 %). Such a feature cannot be observed for a weakly magnetized model (B 0 = 10 11 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 T E2,20 (e.g., Thorne (1980)). A E2 20 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, Initial rotation parameter is both set to be β = 0.1 %, while the precollapse magnetic field is taken as 10 12 G (left panel) and 10 11 G (right panel), respectively (from ). The total wave amplitudes are shown by the red line, while the contribution from the magnetic fields and from the sum of hydrodynamic and gravitational parts are shown by blue and green lines, respectively (e.g., equation (15) and equations (11,13)). Note that the bounce GW signals ((t − t b 20 ms) are not affected by the magnetic fields significantly, and they are categorized into the so-called type I or II waveforms.
The right panel of Figure 23 shows contributions to the total GW amplitudes (equation (6)) for the strongly magnetized model (left panel of Figure 22), in which the lefthand-side panels are for the sum of the hydrodynamic and gravitational part, namely Right panel shows the sum of the hydrodynamic and gravitational parts (indicated by "hyd and grav" in the left-hand side) and the magnetic part (indicated by "mag" in the right-hand side), respectively. The top and bottom panels represent the positive and negative contribution (indicated by (+) or (-)) to A E2 20 , respectively (see text for more detail). The side length of each plot is 4000(km)x8000(km). This figure is taken from . equations (11,13)), and the right-handside panels are for the magnetic part, namely log ±f E2 20 (mag) (right top(+)/bottom(−)) (e.g., equation (15)). 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., +ρ * W 2 v r 2 in equation (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., +b φ 2 in equation (15)), which dominantly contribute to drive MHD explosions. Figure 24 shows the GW spectra for a pair models of B12X5β1 and B11X5β1 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 1 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 (Kawamura 2006;Kudoh et al. 2006) (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 10 51 erg.   (Kuroda & LCGT Collaboration 2010), and Fabry-Perot type DECIGO (Kawamura 2006;Kudoh et al. 2006). It is noted that h char is the characteristic gravitational wave strain defined in Flanagan & Hughes (1998). The supernova is assumed to be located at the distance of 10 kpc. This figure is taken from .

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., Sato & Suzuki (1987), see Raffelt (2012Raffelt ( , 2002 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 (Engel et al. 2003) that is an upgrade of SNO) or liquid-Argon detectors (GLACIER (Autiero et al. 2007) that is an upgrade of ICARUS) will also play an important role (see Scholberg (2010)) 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 ν e burst. At this neutronization epoch (in the first 10-20 ms postbounce), the largest difference in the light-curves among ν e and the rest of the neutrino species can be seen (e.g., Figure 1 in Kachelrieß et al. (2005)). Since neutrino oscillations take place only when there is a difference in the neutrino fluxes among different species (see section 3 Kotake et al. (2006) 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ν e signals that are produced by the inverse beta processν e + p → e + + n. It should be noted that the detection the ν e bursts is important for studying not only neutrino oscillations but also for determining the direction and distance to the source (Kachelrieß et al. 2005). In a megaton-classCherenkov detector in which gadolinium is added to catch neutrons (Beacom & Vagins 2004), and large liquid-Argon detectors (like ICARUS and GLACIER cited above) are expected to be powerful ν e 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 M ⊙ (e.g., Kitaura et al. (2006); Thompson et al. (2003); Kachelrieß et al. (2005); Hüdepohl et al. (2010);Fischer et al. (2010) and references therein), one generally finds L νe ∼ Lν e > L νx and E νx > Eν e E νe with L ν i and E ν i representing luminosity and average neutrino energy for the neutrino species of i = ν e ,ν e , ν x (;ν µ , ν τ and their anti-particles). The accretion phase is followed by the PNS 31 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 . In the PNS cooling phase, the luminosities of all neutrino species become similar L νe ∼ Lν e ∼ L νx (with the energy hierarchy of E νx > Eν e > E νe ) and decrease monotonically with time. Due to the similarity in the luminosities and spectra amongν e and ν x 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., Mikheev & Smirnov (1986)) 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., Takahashi et al. (2003a); Kachelrieß et al. (2005)) and on the Earth matter effects (e.g., Lunardini & Smirnov (2001); Dighe et al. (2004)). 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., Schirato & Fuller (2002); Lunardini & Smirnov (2003); Tomàs et al. (2004); Dasgupta & Dighe (2007); Lunardini & Peres (2008), see Duan & Kneller (2009) for a recent review). The shock passage to the so-called high-resonance regions may be observed as a sudden decrease in ν e events in the case of normal mass hierarchy (orν e 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., Dighe & Smirnov (2000); Takahashi et al. (2003b); Tomàs et al. (2004);Fogli et al. (2005)).
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 non-linearity and the presence of multi-angle 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., Chakraborty et al. (2011);Dasgupta et al. (2011), see however Suwa et al. (2011);Pejcha et al. (2011)), they are predicted to emerge as a distinct observable feature in their energy spectra (see Raffelt & Smirnov 2007;Duan et al. 2008Duan et al. , 2010Dasgupta 2010, 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 not other energy-loss channel but for the ordinary neutrinos in the context of the Standard Model of particle physics (e.g., Raffelt (1996Raffelt ( , 2002). 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., Schechter & Valle (1981); Lim & Marciano (1988); Akhmedov & Khlopov (1988) and references therein) and numerically (e.g., Ando & Sato (2003)), which can transform some of the prompt ν e burst intoν e in highly magnetized supernova envelopes, leading to a hugē ν e 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., Totani et al. (1998);Lunardini & Smirnov (2003); Takahashi et al. (2003a); Tomàs et al. (2004); Kachelrieß et al. (2005); ; Duan et al. (2008);Hüdepohl et al. (2010); Cherry et al. (2011);Sarikas et al. (2011) and references therein). Apart from a simple parametrization to mimic anisotropy and Stochasticity of the shock (e.g., Friedland & Gruzinov (2006); Choubey et al. (2007); Reid et al. (2011)), there have been not so many studies focusing how global asymmetry of the explosion dynamics obtained in recent supernova simulations (e.g., section 2 and 3) would affect the neutrino oscillations (see, however, Duan & Kneller (2009);Dasgupta et al. (2008)). 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 (Lund et al. 2010;Brandt et al. 2011), in which the collective or MSW effects are rendered to be treated in an approximate manner. By analyzing the 2D results of a 15 M ⊙ model by  who included one of the best available neutrino transfer approximations (e.g., Table 1), Lund et al. (2010) 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 M ⊙ models by Ott et al. (2008) who reported the first 2D multi-angle transport simulations, Brandt et al. (2011) 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 (Kotake et al. 2003a) as well as in its light curves (Janka & Mönchmeyer 1989;Brandt et al. 2011).
In the following, we briefly summarize our findings (Kawagoe et al. 2009) 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 ), we calculated numerically the neutrino flavor conversion in the highly non-spherical envelope through the MSW effect. The neutrino transport was simply treated by a leakage scheme (e.g., Takiwaki et al. (2009)). 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 section 3.0.2, 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 the mixing angle of θ 13 . 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., Kawagoe et al. (2009)).

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 break-out 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 ρ res ∼ 1.4 × 10 3 g/cm 3 ∆m 2 10 −3 eV 2 where ∆m 2 is the mass squared difference, E ν is the neutrino energy, Y e 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. ∆m 2 and θ correspond to ∆m 2 13 and θ 13 at the H-resonance and to ∆m 2 12 and θ 12 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(10 3 )g/cm 3 at ∼ 0.5 s, and the L-resonance, ∼ O(1)g/cm 3 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 Lresonances, respectively (e.g., Fogli et al. (2005); Tomàs et al. (2004)). 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 (Fogli et al. (2005); Tomàs et al. (2004)). The progenitor of the MHD models, possibly linked to long-duration gamma-ray bursts, is more compact due to a chemically homogeneous evolution (Yoon & Langer 2005b;Woosley & Heger 2006), which is also the reason for the early shock-arrival to the resonance regions.  .2 (f) s, respectively. Note that a sharpness of the shock along the polar direction is modified from the one in the MHD simulations (green line showing an original density profile for 2.0 s), which is made to be inevitably blunt due to the employed numerical scheme using an artificial viscosity to capture shocks (taken from Kawagoe et al. (2009)). Figure 27 is the energy spectra ofν e 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ν e can remain nonzero when the steep decline of the density at the shock front changes the resonance into non-adiabatic (see Kawagoe et al. (2009) 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, ρ res , is proportional to E −1 ν . 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., Fig 11 in Tomàs et al. (2004)), is thanks to the mentioned early shock-arrival to the resonance layer, peculiar for the MHD explosions.  -ν e spectra at the surface of the star in the case of (a), (b) and (c). We assume the inverted mass hierarchy and sin 2 2θ 13 = 10 −3 . The solid lines (red lines) and the dotted lines (blue lines) are the spectra of polar direction and equatorial direction, respectively. Note in this figure that the neutrino luminosity and spectra are modeled to be the same as those obtained in a 1D full-scale numerical simulation by the Lawrence Livermore group (Wilson 1985). These figures are taken from Kawagoe et al. (2009). is for the polar and the equatorial direction, respectively. The supernova is assumed to be located at the distance of 10 kpc. The error bars represent the 1σ statistical errors only. Right panel is for the ratio of high-to low-energy events which can be a useful quantity to characterize the change (dip in the right panel) of neutrino signals by the shock effect (see Kawagoe et al. (2009) for details).
Finally we exploratory discuss possible impacts of the spectral swap betweenν e and ν x , which is one of a probable outcome of the collective neutrino oscillations (Fogli et al. 2007). In our model, the original neutrino flux ofν e is larger than that of ν x 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 Kawagoe et al. (2009) 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 multi-angle approach (e.g., Duan et al. (2010)), which is a major undertaking.  ) and their expected features of the emergent multimessengers (horizontal direction). Regarding the photon messengers, possible signatures associated with nucleosynthesis are only partly presented (see section 1 for a number of important findings regarding the lightcurve and spectra modeling).

Summary and Discussion
The aim of writing this article 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 multi-messengers, which is far beyond the scope of this review. We boldly present the incomplete table 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 (Maeder & Meynet 2000;Heger et al. 2005), 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 non-linear 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 non-linear epoch also change stochastically with time, which is expected to be detected by currently running detector such as by IceCube and SK (Brandt et al. 2011;Lund et al. 2010) 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., Hüdepohl et al. (2010);Fischer et al. (2010)). 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 (Pruet et al. 2006;Fröhlich et al. 2006;Wanajo & Janka 2011). When the revived supernova shock (successfully) attains the kinetic energy as energetic as 10 51 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 non-uniform 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 (∼ 10 51 erg, e.g., Table 1). Moreover, the softer EOS, such as of the Lattimer & Swesty (1991) (LS) EOS with an incompressibility at nuclear densities, K, 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 (K = 240 ± 20 MeV, Shlomo et al. (2006)), the soft EOS may not account for the recently observed massive neutron star of ∼ 2M ⊙ (Demorest et al. 2010) (see the maximum mass for the LS180 EOS in O'Connor & Ott 2011;Kiuchi & Kotake 2008). With a stiffer EOS, the explosion energy may be even lower as inferred from  who did not obtain the neutrinodriven explosion for their model with K = 263 MeV 32 . 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 (Sagert et al. 2009;Fischer et al. 2011), viscous heating by the magnetorotational instability (section 3.0.3), or energy dissipation via Alfvén waves (Suzuki et al. 2008). In addition to the 3D effects, GR is expected to help the onset of multi-D neutrino-driven explosions Kuroda et al. 2012) 33 . At present, the most up-to-date neutrino transport code can treat the multi-energy and multi-angle transport in 2D ) and even in 3D simulations (Sumiyoshi & Yamada 2012) 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 life-time!) presumably by using the exascale 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 (Guilet et al. 2010). What determines mode couplings between small-scale convection eddies and large-scale SASI modes ? To apply the theory of turbulence (Murphy & Meakin 2011) 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 (Hempel & Schaffner-Bielich 2010;Furusawa et al. 2011;Shen et al. 2011). 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 34 (indicated by "excess for 32 On the other hand, they obtained 2D explosions for Shen EOS (K = 281MeV) (H-T. Janka, private communication). 33 e.g., Mueller et al. (2012) in 2D simulations but with detailed neutrino transport or Kuroda et al. (2012) in 3D simulations but with a very approximate neutrino transport 34 It should be noted that the feature can be seen even in a slowly rotating progenitor that agrees with recent stellar evolution calculation (Maeder & Meynet 2000;Heger et al. 2005) (see section 2.2.2 for more details). 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 (Ott et al. 2006a), 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 (Janka & Mönchmeyer 1989;Kotake et al. 2003a;Brandt et al. 2011) (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 (rotationinduced) 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).

equator" in
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. (2011) 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. (2006), the disappearing neutrino signals that mark the epoch of BH formations also can be a target of SK and IceCube (IceCube Collaboration et al. 2011). When quarks and hyperons appear in the postbounce core, a neutrino burst produced by the sudden EOS softening and by the subsequent rebounce (Sagert et al. 2009) is likely to be detected by IceCube for a galactic event (IceCube Collaboration et al. 2011). The interval between core-bounce and the BH formation depends on the details of exotic physics in the super-dense PNS cores (Nakazato et al. 2010a,b). All of these observational signatures should provide an important probe into the so-called dense QCD region in the QCD phase diagram (e.g., Fukushima & Hatsuda (2011) 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 ∼ 10 11 G), the MHD mechanism can produce bipolar explosions along the rotational axis predominantly driven by the field wrapping processes (section 3.0.2). 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 (section 3.0.3 and 3.0.4). 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 (Leonor et al. 2010) 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 ν e converts to that ofν e , which is thus expected to be a probe into the magnetic moment of neutrinos (e.g., Ando & Sato (2003)). Finally the MHD explosions are considered to be a possible r-process cite (Nishimura et al. 2006;Fujimoto et al. 2007) 35 . 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 10 15 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., Wheeler et al. (2000); Thompson et al. (2004); Uzdensky & MacFadyen (2007) and see references therein). This picture, often referred as collapsar (e.g., MacFadyen & Woosley (1999);MacFadyen et al. (2001)), has been the working hypothesis as a central engine of 35 See collective references in Wanajo & Janka (2011) for other plausible r-process cites.
long-duration GRBs for these 10 years (see references in Dai & Lu (1998); Thompson et al. (2004); Uzdensky & MacFadyen (2007); Bucciantini et al. (2007); Metzger et al. (2011) for other candidate mechanism including magnetar models). However it is still controversial whether the generation of the relativistic outflows predominantly proceeds via the neutrinoheating 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 multi-angle neutrino transfer for treating highly anisotropic neutrino radiation from the accretion disk. 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 article 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 article can progress in a stepby-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 (Antonioli et al. 2004). 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 multi-messenger astronomy. The interplay between the detailed numerical modeling, the advancing supercomputing resources, and the multi-messenger 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 endeavours 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 world-wide, multi-disciplinary 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 article.