Magnetorotational Explosion of A Massive Star Supported by Neutrino Heating in General Relativistic Three Dimensional Simulations

We present results of three-dimensional (3D), radiation-magnetohydrodynamics (MHD) simulations of core-collapse supernovae in full general relativity (GR), for the first time, with spectral neutrino transport. In order to study the effects of progenitor's rotation and magnetic fields on the dynamics and the emergent neutrino signals, we compute three models, where the precollapse rotation rate and magnetic fields are included parametrically to a 20 M$_{\odot}$ star. While we find no shock revival in our two non-magnetized models during our simulation times ($\sim500$ ms after bounce), the magnetorotationally driven shock expansion immediately initiates after bounce in our rapidly rotating and strongly magnetized model. We show that the expansion of the magnetorotationally-driven flows toward the polar directions is predominantly driven by the magnetic pressure, whereas the shock expansion toward the equatorial direction is supported by neutrino heating. After $\sim 150$ ms postbounce, the magnetorotationally-driven bipolar flows eventually cascades to an one-sided, unipolar flow. Our detailed analysis indicates that the growth of the so-called kink instability may hinder the collimation of jets, resulting in the formation of broader outflows. Furthermore we find a clear signature of the lepton-number emission self-sustained asymmetry (LESA), though only in the MR explosion model, whose asymmetry is consistent with the unipolar explosion morphology. We also report several unique neutrino signatures, which are significantly dependent on both the time and the viewing angle, if observed, possibly providing a rich information regarding the onset of the magnetorotationally-driven explosion.

Positive effects of rotation in favor of the onset of explosion include the larger shock radius due to the centrifugal force (Nakamura et al. 2014), vigorous spiral SASI activity (Summa et al. 2018), and energy transport via the rotational instability (Takiwaki et al. 2016). On the other hand, rotation weakens the explodability because it results in a more extended and cooler PNS, which reduces the neutrino luminosities and energies (Marek & Janka 2009). These studies show that the impact of rotation on the neutrino mechanism depends sensitively on the precollapse rotation rate. Supported by the outcomes from these state-of-the-art multi-D simulations, we are now reaching a broad consensus that the multi-D neutrino mechanism is the most promising way to account for canonical CCSNe with the explosion energies of the order of 10 51 erg (≡ 1 Bethe, 1 B in short) or less.
The neutrino mechanism, however, is likely to fail in a subclass of CCSNe with very energetic explosion of ∼ 10 B, which is termed as hypernova (HN) (Iwamoto et al. 1998). Observationally a HN is associated with collapse of very massive star typically with 30 − 40 M ⊙ in the main sequence stage (Tanaka et al. 2009).
The magnetorotational (MR) explosion mechanism relies on the extraction of the rotational free energy from the central compact objects via the magnetic fields (see also Blandford & Znajek 1977;McKinney 2006, in various contexts).
Rapid rotation of the iron core is a necessary condition for the working of the MR mechanism (see Kotake et al. 2006, for collective references of early studies therein). In the collapsing core, the magnetic fields are amplified to dynamically relevant field strengths by rotational winding and/or magnetorotational instability (MRI) (Akiyama et al. 2003;Obergaulinger et al. 2009;Masada et al. 2015;Rembiasz et al. 2016). After bounce, the strong magnetic pressure launches the jets along the rotational axis (Ardeljan et al. 2000;Burrows et al. 2007;Takiwaki et al. 2009;Scheidegger et al. 2010;Winteler et al. 2012;Mösta et al. 2014;Obergaulinger et al. 2014).
The highly aspherical explosion is also observationally supported by the analysis of the line profiles (e.g., Maeda et al. 2008). Note in the non-rotating progenitors, Obergaulinger et al. (2014) were the first to point out that magnetorotationally-driven pressure support in the gain region (via turbulence) fosters the onset of neutrino-driven explosion. This result clearly presents evidence that implementation of sophisticated neutrino transport is needed for a quantitative study of magnetorotationally-driven CCSN modeling.
In the context of purely neutrino-driven models (without magnetic fields), it becomes certain that twodimensional (2D) simulations overestimate the explodability for a wide variety of progenitor (Hanke et al. (2012; Couch (2013); Dolence et al. (2013); Takiwaki et al. (2014)). In order to correctly capture the evolution and dynamics of the postshock turbulence, three-dimensional (3D) modeling is required. The higher explodability in 2D is also reported in MR models. Mösta et al. (2014) has shown that a full 3D model leads to the formation of the less collimated (bipolar) jets than those in the counterpart octant symmetry model, which mimics 2D. They pointed out that the less collimated outflow in 3D is an outcome of the so-called |m| = 1 kink instability (Lyubarskii 1999;Begelman 1998;Narayan et al. 2009). It has been demonstrated that the kink instability displaces the jet center from the rotational axis and prevents the magnetic fields amplification preferentially on the axis (see also Li 2000). More recently, Obergaulinger & Aloy (2020) has reported the first 3D special-relativistic MHD simulations with spectral neutrino transport. Their 3D models showed slightly longer explosion times, although the explosion energy and ejecta mass were higher and larger, respectively, compared to those in the counterpart 2D models. Any remarkable non-axisymmetic instabilities, including the kink instability, were not seen in the 3D models of Obergaulinger & Aloy (2020), which is in contrast with Mösta et al. (2014). Therefore the multi-D effects in MHD models are still controversial, due partly to the limited number of full 3D MHD CCSN simulations reported so far (Mikami et al. 2008;Scheidegger et al. 2010;Mösta et al. 2014;Obergaulinger & Aloy 2020).
In this paper, we report first results of full 3D-GR, magnetorotational core-collapse simulations of a 20 M ⊙ star with spectral neutrino transport. We calculate three models, rotating magnetized, rotating non-magnetized, and non-rotating non-magnetized models. Our results show that the MR explosion occurs in the rotating magnetized model shortly after core bounce, whereas the shock revival is not obtained in both non-magnetized models during our simulation time (∼ 500 ms after bounce). While our results basically confirm the previous results (Mösta et al. 2014), our findings include detailed analysis of the kink instability, the LESA activity in the MR explosion, and the neutrino signals from the 3D-GR MHD models with self-consistent neutrino transport.
This paper is organized as follows. Section 2 starts with a concise summary of our GR-MHD neutrino transport scheme, which is followed by the initial setup of the simulation. The main results and detailed comparison with previous studies are presented in Section 3. We summarize our results and conclusions in Section 4. Note that the geometrized unit is used in Section 2, i.e., the speed of light, the gravitational constant and the Planck constant are set to unity: c = G = h = 1, and cgs unit is used in Section 3. The metric signature is (−, +, +, +).
Greek indices run from 0 to 3 and Latin indices from 1 to 3, except ν and ε denoting neutrino species and energy, respectively. We also use a conventional expression for spatial coordinates (x 1 , x 2 , x 3 ) = (x, y, z).

SETUP
In our full GR radiation-MHD simulations, we solve the evolution equations of metric, magnetohydrodynamics, and energy-dependent neutrino radiation. Each of the evolution equations is solved in an operator-splitting manner, while the system evolves self-consistently as a whole satisfying the Hamiltonian and momentum constraints (Kuroda et al. 2012(Kuroda et al. , 2016b.

Basic ν-GRMHD Equations
Regarding the metric evolution, we evolve the standard BSSN variablesγ ij , w(= e −2φ ) (Marronetti et al. 2008),Ã ij , K, andΓ i (Shibata & Nakamura 1995;Baumgarte & Shapiro 1999). Here φ ≡ ln(γ)/12 with γ = det(γ ij ). The gauge is specified by the "1+log" lapse and by the Gamma-driver-shift condition. Evolution equation of these variables are solved with a fourthorder finite difference scheme in space (Zlochower et al. 2005) and with a fourth-order Runge-Kutta time integration. In appendix 4, we show results of the polarized Gowdy wave test (Alcubierre et al. 2004) to show the fourth-order convergence of our metric solver.
In the radiation-magnetohydrodynamic part, the total stress-energy tensor T αβ (total) is expressed as where T αβ (matter) , T αβ (EM) , and T αβ (ν,ε) are the stress-energy tensor of matter, electro-magnetic, and energy (ε) dependent neutrino radiation field of specie ν, respectively. We consider all three flavors of neutrinos (ν e ,ν e , ν x ) with ν x representing heavy-lepton neutrinos (i.e. ν µ , ν τ and their anti-particles). ε represents the neutrino energy measured in the comoving frame. In this paper, although we omit to describe detailed evolution equations of the neutrino radiation field (we refer the reader to Kuroda et al. 2016b), we solve spectral neutrino transport of the zeroth and first order radiation momenta, based on the truncated moment formalism (Shibata et al. 2011) employing an M1 analytical closure scheme.
In the following, we briefly describe our GRMHD formulation. The stress-energy tensor of electro-magnetic field T αβ (EM) is expressed as where F αβ is the electromagnetic field tensor. Since we currently consider the ideal MHD case, Maxwell's equations are written in terms of the dual tensor F * αβ = 1 2 ǫ αβγδ F γδ as We define the magnetic field four vector b α as below with ǫ αβγδ and u α being the Levi-Civita tensor and matter four velocity, respectively. In addition, for later convenience, the magnetic field three vector B i should also be introduced as where W = −u µ n µ is the Lorentz factor (do not confuse with w = e −2φ of geometrical variables) and n µ = (−α, 0, 0, 0) is a unit vector normal to the spacelike hypersurface foliated into the spacetime. Then, using the orthogonality condition B α n α = 0, the time and spacial components of Eq. (3) can be rewritten as i.e., the solenoidal constraint of B i , and respectively, where v i ≡ u i /u t . Additionally to the evolution equation (7) of the magnetic field, we solve the following ideal hydrodynamic equations (see, e.g. Shibata & Sekiguchi 2005;Gammie et al. 2003) including electron number conser- and where On the right hand side of Eq.(10), D i represents the covariant derivative with respect to the three metric γ ij . ρ is the rest mass density and h = 1+e mat +P mat /ρ is the specific enthalpy of matter (composed of baryons, electrons, and photons) with e mat and P mat being the specific internal energy and pressure of matter, respectively. b 2 = b α b α , P tot = P mat + P mag is the total pressure, P mag = b 2 /2 is the magnetic pressure, Y e ≡ n e /n b is the electron fraction (n e and n b are the number densities of electrons and baryons, respectively), and m u is the atomic mass unit. P mat (ρ, s, Y e ) and e mat (ρ, s, Y e ) are given by an equation of state (EOS) with s denoting the entropy per baryon.
We thus evolve the following magnetohydrodynamic and radiation conservative variables where (E (ν,ε) , F (ν,ε) i ) are the zeroth and first order moments of neutrino radiation (Shibata et al. 2011;Kuroda et al. 2016b). Every time we update the conservative variables Q, we obtain the following primitive variables by Newton's method.

Constrained Transport
We solve the conservation equations (8)-(11) using the HLL scheme (Harten et al. 1983). Meanwhile the induction equation (7) is solved by a constrained transport (CT) method (Evans & Hawley 1988) to satisfy the solenoidal condition Eq. (6). For the CT method, we also utilize the HLL scheme when we reconstruct the electric field that will be mentioned shortly here. To solve the (HLL) Riemann problem, we need to evaluate the left and right states at cell surface. The left and right states are interpolated from cell centered primitive variables P and some of the metric terms (w, α, β i , γ ij ), which are needed to evaluate the full conservative variables Q, by a spatial reconstruction. We use the Piecewise Parabolic Method (PPM) for the spatial reconstruction (Colella & Woodward (1984) or Hawke et al. (2005) for more suitable upwind reconstruction method in GR.). After the spatial reconstruction step, we calculate the fastest left-and right-going wave speeds (e.g. Antón et al. 2006) and the HLL fluxes.
We also introduce the electric field E i defined by for the CT method. Then the equation (7) can be rewritten as Employing a usual staggered mesh algorithm, we define the magnetic field B i and the electric field E i at cell surface and edge, respectively, while the rest of variables are defined at cell center. For instance, B x and E x are defined at (i+1/2, j, k) and (i, j+1/2, k+1/2), respectively, where (i, j, k) denotes the cell center and, e.g. j + 1/2 represents a displaced position from cell center along y axis with a half cell width. As in the usual manner, the electric field E i defined on the cell edge is obtained from the HLL flux for B i , corresponding to the advection term in Eq. (7). We use the nearest four electric fields defined on the cell surface, i.e., corresponding terms in the HLL flux, and take their simple arithmetic average 1 . Our numerical grid employs a fixed nested structure and there is a boundary between different numerical resolutions. Therefore we apply a refluxing procedure both for the HLL fluxes and the electric field E i (Kuroda & Umeda 2010) before solving Eqs.(8)-(11) and (15) to satisfy the conservation law and solenoidal constraint in the whole computational domain.

Initial Setup
We study the frequently used solar-metallicity model of the 20 M ⊙ star "s20a28n" from Woosley & Heger (2007). Although one of our final aims is to understand the hypernova explosion mechanism of very massive stars ( 30 M ⊙ ), this progenitor star is widely used in previous studies (e.g., Melson et al. 2015;Ott et al. 2018;O'Connor & Couch 2018;Burrows et al. 2019) and our non-rotating, non-magnetized model (see below) could thus be a reference model to calibrate our 3D ν-GRMHD code. For the nuclear EOS, we use SFHo of Steiner et al. (2013). The 3D computational domain is a cubic box with 3 × 10 4 km width in which nested boxes with 10 refinement levels are embedded. Each box contains 128 3 cells and the minimum grid size near the origin is ∆x = 458 m. In the vicinity of the stalled shock at a radius of ∼ 100 km, our resolution achieves ∆x ∼ 1.9 km, i.e., the effective angular resolution becomes ∼ 1 • . The neutrino energy space ε logarithmically covers from 1 to 300 MeV with 12 energy bins. Regarding neutrino opacities, the standard weak interaction set in Bruenn (1985) plus nucleon-nucleon Bremsstrahlung (Hannestad & Raffelt 1998) is taken into account. The original progenitor model "s20a28n" assumes neither rotation nor magnetic fields. We thus artificially add them to the non-rotating progenitor model. We employ a widely used cylindrical rotational profile where u φ ≡ u 2 x + u 2 y . ̟ 0 and Ω 0 indicate the size and angular frequency of a rigidly rotating central cylinder, respectively. This rotational profile gives the angular frequency falling off with For the initial magnetic fields that should satisfy the solenoidal constraint, we use the following purely toroidal vector potential Here (r, θ, φ) denote the usual coordinates in the spherical polar coordinate system. By defining these vector potentials on the cell edge and taking their curl B = ∇×A, the magnetic field defined on the numerical cell surface automatically satisfies the solenoidal constraint. This vector potential gives nearly uniform magnetic field parallel to the rotational axis (i.e. z-axis) for r R 0 and dipolar magnetic field for r R 0 .
We set ̟ 0 = R 0 = 10 8 cm corresponding roughly to the iron core size at the precollapse stage. We calculate three models: rotating magnetized, rotating nonmagnetized, and non-rotating non-magnetized. For the rotating models, we set Ω 0 = 1 rad s −1 . This value is very reasonable compared to the one of a rotating 20 M ⊙ model in Heger & Langer (2000) that gives Ω 0 ∼ 3 rad s −1 . Regarding the magnetic field strength at origin, we set B 0 = 10 12 G that can be amplified strongly enough to affect the dynamics through simple linear amplification mechanisms, i.e., compression and rotational wrapping, during collapse and is also widely used in most of previous MHD simulations (Burrows et al. 2007;Takiwaki et al. 2009;Scheidegger et al. 2010;Mösta et al. 2014;Obergaulinger & Aloy 2020). Three models are labeled as R0B00, R1B00, and R1B12, where the integer after R denotes Ω 0 . B00 and B12 represent B 0 = 0 and 10 12 G, respectively.

RESULTS
In this Section, we explain our main results. Sections 3.1 and 3.2 are devoted to explaining general hydrody-namic properties in the post-bounce evolution. In Section 3.3, we discuss non-axisymmetric instabilities in the PNS and MHD outflow, which is relevant to the neutrino signals in Section 3.4. The role of neutrino heating in the MR mechanism is addressed in Section 3.5. We explain the LESA activity in our MR explosion model in Section 3.6.

Postbounce Evolution
We begin with a brief description of the postbounce evolution of all the three models in this work. After the start of calculation (t = 0), the core bounce occurs at t = 0.261 s, 0.264 s, and 0.265 s for model R0B00, R1B00, and R1B12, respectively. The central (maximum) rest mass density ρ max reaches 4.42 × 10 14 g cm −3 (model R0B00), 4.37 × 10 14 g cm −3 (R1B00), and 4.35 × 10 14 g cm −3 (R1B12). A monotonic feature that rapid rotation and high initial magnetic field delay the bounce time and decrease ρ max is due to the stronger centrifugal force and magnetic pressure at bounce. The lapse function at the center also shows the similar trend, where it takes the smallest and highest value for model R0B00 and R1B12, respectively. For the computed three models, ρ max and the (minimum) lapse function evolve with time after bounce, keeping the above trend at bounce (for example, smallest ρ max for model R2B12 relative to other models).
To visualize the postbounce evolution, Fig. 1 shows the volume rendered entropy (upper three panels) and inverse of the plasma β for model R1B12 in the logarithmic scale (lower three panels) at selected postbounce times (t pb ). Here the plasma β is defined by the ratio of the gas to the magnetic pressure, i.e., β ≡ P gas /P mag . After bounce, the formation of the bipolar flow can be clearly seen in the left panels. Inside the expanding blobs, the magnetic pressure dominates over the gas pressure as shown by the yellowish region (log 10 β −1 0.5) in the lower panels. This is a clear evidence of the magnetorotationally-driven shock revial for model R1B12. As an important 3D effect, we see that the shock morphology is less collimated compared to the previous 2D axisymmetric studies, although similar initial rotation and magnetic fields were adopted (Burrows et al. 2007;Takiwaki et al. 2009;Mösta et al. 2014;Obergaulinger & Aloy 2017). The middle panels show that the jet head is not aligned with the rotational axis at t pb ∼ 100 ms, but is displaced from the axis (indicated by the deviation from the white line). In Sec. 3.3, we will discuss the reason of this in more detail. Fig. 2 shows the volume rendered entropy structure for models R0B00 (upper row panels) and R1B00 (lower row panels) from t pb ∼ 245 ms to ∼ 500 ms. Comparing with the uni-/bipolar like structure seen in the magnetized model R1B12, the shock morphology of these two nonmagnetized models is obviously different. Models R0B00 and R1B00 show roundish and oblate shock morphology, respectively. During our simulation time up to t pb ∼ 500 ms, we do not find a shock revival in these two nonmagnetized models. Fig. 3 shows the maximum (thick lines) and averaged (thin) shock radii in top panel and the time evolution of the diagnostic explosion energy E exp in the bottom panel for model R0B00 (black line), R1B00 (blue line), and R1B12 (red line), respectively. Here, E exp is defined by

Shock Wave Evolution
which is analogous to Eq.
(2) of Müller et al. (2012), but takes into account the additional contribution from magnetic fields. From the top panel, one can see that the shock revival is not obtained for the non-magnetized models R0B00 (black line) and R1B00 (blue line) for the simulation time, whereas the shock propagates outwards in the magnetized model R1B12 (red line). The shock is slightly energetized at t pb ∼ 180 ms for model R1B12 and at t pb ∼ 200 ms for models R0B00/R1B00, when the Si/O interface accretes onto the shock. This leads to the runaway shock expansion for model R1B12, whereas it only results in the slight shock expansion maximally up to the radius of ∼ 400 km for model R1B00, gradually shifting to the standing shock later on (see, blue and black lines). The time when the Si/O interface accretes onto the shock differs about ∼ 20 ms between model R1B12 and the other two models. The time lag is because of the difference in the (maximum) shock position (∼ 4 × 10 7 cm) at t pb ∼ 180 ms. Since the typical accretion velocity is ∼ 2 × 10 9 cm s −1 there, this can be translated into the crossing time of ∼ 20 ms, which is consistent with the time difference.
The diagnostic explosion energy in the bottom panel basically correlates with the shock evolution. In the successful explosion model R1B12, the diagnostic explosion energy increases significantly faster than the other two non-explosion models already at ∼ 20 ms after bounce. It reaches ∼ 10 50 erg around t pb ∼ 100 ms. The value E exp ∼ 10 50 erg at the time when the shock reaches R shock ∼ 1000 km is very similar to the ones in previous 2D (Takiwaki et al. 2009;Obergaulinger & Aloy 2017) and 3D (Obergaulinger & Aloy 2020) studies with the similar initial rotation and magnetic fields strength. In the non-magnetized models R0B00 and R1B00, E exp temporally reaches ∼ 10 49 erg at T pb ∼ 220 ms when the Si/O interface accretes and a temporal shock expansion occurs, though it soon decreases.
We can also find a typical signature of SASI in the evolution of shock radii. From top panel in Fig. 3, a time modulation is visible in the maximum shock radii, particularly in the model R0B00 (thick black line) for t pb 100 ms. Such a modulation reflects the appearance of SASI Foglizzo et al. 2006). To see more quantitatively the shock morphology and also the dominant SASI mode, we plot time evolution of normalized mode amplitudes A ℓm ≡ c ℓm /c 00 of spherical polar expansion of the shock surface R shock (θ, φ) for several dominant modes in Fig. 4. Here we adopt the same definition for c ℓm as in Burrows et al. (2012) with ℓ and m representing the quantum number with respect to the real spherical harmonics of Y m ℓ , respectively. In the top panel, the dominant mode is (ℓ, m) = (2, 0) (black line) for the first ∼ 120 ms after bounce. Since its sign is positive, the shock morphology is prolate as also shown in the left and center columns in Fig. 1. -Snapshots of the volume rendered entropy (upper panels) and inverse of the plasma β in the logarithmic scale (log β −1 , lower panels) for model R1B12. From left to right panels, the post-bounce time of t pb ∼ 56 ms, 100 ms, and 250 ms are depicted, respectively. In the upper panels, the central bluish spherical/spheroidal object roughly corresponds to the unshocked PNS core. Note that the inclination angle of the coordinates is not fixed in each time snapshot to visualize the expansion morphology more clearly. The white line denotes the length scale that is parallel to the rotational axis (z-axis).  Fig. 1, but for only entropy of non-magnetized models at different time slices t pb ∼ 245, 370, and 500 ms. The upper and lower panels are for model R0B00 and R1B00, respectively. Note again that the white line denotes the length scale that is parallel to the rotational axis (z-axis).  However, for t pb 120 ms in the same model R1B12, (ℓ, m) = (1, 0) (red line) gradually takes over as the dominant term with its sign being negative. Therefore the shock morphology at the end of simulation time is unipolar toward the negative z-axis, which is again consistent with the right column in Fig. 1. In the middle panel, R1B00 shows that A 20 becomes negative for t pb 50 ms which reflects a rotating oblate spheroid (see, bottom panels in Fig. 2). At the same time, (ℓ, m) = (1, ±1) (blue and green lines) also show comparable amplitudes with that of (2,0), but with clear quasi-periodic oscillations. Between these two |m| = 1 modes, i.e., (ℓ, m) = (1, 1) and (1, −1), a phase shift seemingly with ∼ π/2 exists which indicates that the spiral SASI motion appears (Blondin & Mezzacappa 2007). In the non-rotating model R0B00, all the three modes with ℓ = 1 and m = 0, ±1 show basically the same amplitude with almost no phase shift up to t pb ∼ 120 ms. Therefore, the dominant SASI mode is the sloshing mode firstly after bounce. Afterward the (1, 0) mode gradually decouples from the other two different azimuthal modes. There seems to be a phase shift of ∼ π/2 between (1,0) (red line) and the other two with (1, ±1) (green and blue). This can be explained by the dominant SASI motion changing from the sloshing motion to the spiral one around t pb ∼ 120 ms. Note that the growth of the spiral SASI in the non-rotating progenitors (Blondin & Mezzacappa 2007) is consistent with the outcomes of previous 3D core-collapse models Kuroda et al. 2016a;Ott et al. 2018).
3.3. Non-axisymmetric instabilities inside the MHD outflow In this subsection, we discuss non-axisymmetric instabilities inside the MHD outflow and their potential impact on the shock evolution. In a 3D-GR model using the similar precollapse rotation rate and magnetic fields to our model R1B12, Mösta et al. (2014) observed the appearance of the kink instability (Lyubarskii 1999;Begelman 1998;Narayan et al. 2009). According to their analysis, the linear growth of the kink instability shortly starts after bounce, which is followed by the non-linear phase already at t pb ∼ 20 ms. At that moment, the jet barycenter showed a significant displacement from the rotational axis, which is one of the main features of the growth of the kink instability, leading to a broader and less energetic outflow compared to the counterpart axisymmetric case. We also check if this instability appears and affects the dynamics of outflow in model R1B12.
The condition |b φ /b z | > ̟/L, i.e., the well-known Kruskal-Shafranov criterion, is the major factor that determines whether the system is unstable to the most dominant screw mode, i.e., for |m| = 1 mode with a condition mb φ < 0. Here L and ̟ denote the inverse of minimal wave number of the unstable mode propagating parallel to the rotational axis and distance from the rotational axis, respectively. In a sufficiently rapidly rotating case, one should also take into account the rotational stabilizing effect that relaxes the Kruskal-Shafranov criterion to |b φ /b z | > Ω̟ (Tomimatsu et al. 2001), where Ω is the angular frequency in geometric units. In our magnetized model R1B12, the toroidal magnetic field dominates over the poloidal one |b φ /b z | > 1 just above the PNS core (z ∼ 10 − 50 km). Such a configuration is usually seen in the magnetized collapse model as the initial poloidal field can be very efficiently converted into the toroidal one mainly through the field-wrapping. As a consequence, the value |b φ /b z |/(Ω̟) inside the MHD outflow reaches O(10 2 ∼ 10 3 ) in our model. We therefore consider that the MHD outflow appeared in our model R1B12 can also be subject to the kink instability.
Following Mösta et al. (2014), we monitor how the barycenter of MHD outflow is displaced from the rotational/magnetic field axis, i.e., z-axis. We take the same definition for the barycenter x i c written by (Mösta et al. 2014) x for i = 1 and 2, where we perform the surface integral ds over the domain with |x, y| ≤ 50 km at z = ±50 km. In addition, to see the mode propagation direction properly in a rotating system, we map the original Cartesian coordinates x i to a rotating framex i bȳ with Q i j being the usual rotation operator with respect to z-axis Θ(t, z) measures the cumulative rotation angle of the system at a given slice z(=z) after core bounce and is given by a following rough estimation whereω(t, z) is the mean angular frequency of the plane. Since the PNS differentially rotates, the rotational angle Θ(t, z) is just a rough measurement. We evaluate the mean angular frequencyω(t, z) simply bȳ where ω z = v z / x 2 + y 2 is the angular frequency measured in the Eulerian frame and we use the rest mass density as a weight. After mapping, we plot the barycenter x i c on the rotating planexȳ atz = ±50 km. In top and bottom panels of Fig. 5, we show the trajectory of barycenter of MHD outflow (solid lines) on the original xy and rotatingxȳ planes at z = ±50 km. To highlight the initial linear growth phase, we show only from the bounce time up to t pb = 30 ms that is indicated by the arrow. In addition, we show direction of b φ averaged over ̟ 40 km by dash-dotted line for reference. Because of our initial purely poloidal magnetic field with dipole-like structure orienting toward positive z-axis, direction of the toroidal component generated after core-collapse mainly through the field wrapping is basically clockwise (b φ < 0) and counterclockwise (b φ > 0) for z > 0 and z < 0, respectively, on these planes. Note that the positive z-axis points toward us and from the condition mb φ < 0 that selects the leading mode to develop, the propagation direction of the most unstable 40 km, which is clockwise (b φ < 0) and counterclockwise (b φ > 0) for z > 0 and z < 0, respectively, on these planes. mode in a comoving frame is expected to be counterclockwise (m = 1) and clockwise (m = −1) for z > 0 and z < 0, respectively. From top panel in Fig. 5, both of the solid lines show basically counterclockwise propagation direction, i.e., m = 1 mode. For the region with z > 0 (blue line in the top panel), the mode propagation direction is opposite to that of b φ and it is consistent with a linear analysis. On the other hand, both solid and dash-dotted lines are pointing toward the same counterclockwise direction on the plane at z = −50 km (red lines), which is not in accordance with the theoretical expectation mb φ < 0. We think that this inconsistency is an apparent as the red solid line in the rotating frame (bottom panel) is showing a clockwise propagation direction opposite to that of b φ . These facts support that the kink instability likely appears, displaces the shock center, and consequently makes the shock morphology broader compared to the corresponding 2D model.
We should also mention another relevant nonaxisymmetric instability that might influence on the growth of the above kink instability. As we have already mentioned, the ratio of rotational to gravitational potential energy after bounce in both of our rotating models reaches several percent, which makes the PNS core being subject to the low-T /W instability (Watts et al. 2005;Saijo & Yoshida 2006). Once the instability appears, it produces an instability mode that propagates in the same direction as the fluid motion, i.e., this time with m = 1 mode in both the northern and southern hemispheres. Therefore, it means that the two different instabilities, namely the low-T /W and kink instabilities, could simultaneously exist possibly with the same m = +1 mode for z > 0 and with two opposite m = ±1 modes for z < 0 breaking the parity between northern and southern hemispheres.
It is beyond the scope of this paper to quantify how the two instabilities coexist, how they affect the PNS core dynamics and the disrution of the bipolar flows as seen in model R1B12. Once the bipolar flows are disrupted, the mass accretion rate becomes higher on the weaker explosion side as a consequence of deflection of mass accretion on the stronger explosion side. This could explain a relatively weak explosion (E exp ∼ 10 50 erg) of model R1B12.
Apparently we need more studies with varying the initial magnetic fields and rotational profiles systematically in order to clarify the disruption mechanism of the MHD outflows.

Rotational Effects on Neutrino Profiles
The time modulation of CCSN neutrino signals reflects the hydrodynamics evolution of the postbounce core (e.g., Tamborra  In this section, we describe how we can make the connection between the core dynamics and neutrino signals. In Fig. 6, we plot the neutrino luminosity L ν (top row) and mean neutrino energy ε ν (bottom row) for specific observer angles for electron type (left), anti-electron type (middle), and heavy lepton type neutrinos (right). Here we evaluate these signals by averaging the neutrino's energy flux at r = 400 km following Tamborra et al. (2014). We choose three observer angles relative to the rotational axis that are denoted by N (north pole), E (equatorial plane, here represented by the positive x-axis), and S (south pole). To prevent too many lines, we plot only spherical averaged values for the non-rotating model R0B00 (solid black line) as it shows basically no significant asymmetry.
Common features among all models are as follows. The neutrino luminosities of all flavours plateau at t pb ∼ 50 − 100 ms. At that moment, the luminosities reach L ν ∼ 6 × 10 52 erg s −1 for ν e andν e and L ν ∼ 3.5 × 10 52 erg s −1 for ν x . Although such values depend on the progenitor star, EOS, and neutrino matter interactions employed, the peak luminosities are in good agreement with those in recent studies with detailed neutrino transport (Müller et al. 2017;O'Connor & Couch 2018;Summa et al. 2018;Vartanyan et al. 2019). The luminosities become nearly constant at t pb ∼ 220 ms when the mass accretion decreases suddenly. We can also see how the progenitor rotation and magnetic field affects the neutrino signals. The non-rotating model R0B00 shows basically the highest luminosity and mean energy in all flavor of neutrinos (see black lines). Meanwhile, the rotating magnetized model R1B12, which explodes shortly after bounce, shows lowest values in both luminosities and mean energies, though there is a slight observer angle dependence. The model R1B00 appears in between them. Such features stem from that the most compact PNS of R0B00 without being subject to the rotational flattening emits higher neutrino luminosities and energies due to its hotter core temperature. On the other hand, the rotating magnetized model R1B12, which shows a lower mass accretion rate due to the centrifugal force and also experiences the mass ejection through bipolar outflow, has a less compact PNS leading to lower neutrino energies and luminosities.
Next, we focus on the viewing angle dependence of the neutrino signals. In Fig. 7, we show a magnified view of Fig. 6 from t pb = 120 ms to 180 ms. In both the rotating models R1B00 and R1B12, the neutrino luminosity and energy observed along the equatorial plane (solid red and blue lines) show the lowest value compared to those along the rotational axis (dash-dotted and dotted lines labeled by N and S). This is because of the rotational flattening of the PNS, where the neutrino-sphere radius along the equatorial plane is located outward than that of the rotational axis, making the neutrino temperature seen along the equator lower than that from the rotational axis (e.g., Kotake et al. 2003;Ott et al. 2008;Harada et al. 2019).
There is yet another neutrino signature for model R1B00. Seen from the equatorial plane (blue solid line), a clear periodic time modulation can be seen. On the other hand, the modulation is hard to be seen from the rotational axis (blue dash-dotted and dotted lines). Furthermore, the degree of the rotational effect differs depending on the neutrino flavour. It is particularly strong in ν e and becomes weaker in order ofν e and ν x . Fig. 7 shows that ν e signals have a time modulation with amplitudes of ∼ 5 × 10 51 erg s −1 and ∼ 0.5 MeV for the luminosity and mean energy, respectively, while those values decrease to ∼ 1 × 10 51 erg s −1 and ∼ 0.2 MeV for ν x . Such a modulation was first discussed in  and is associated with the growth of the so-called low-T /W instability (Ott et al. 2005;Saijo & Yoshida 2006;Watts et al. 2005) and the (one-armed) spiral flows. In both of our rotating models, the ratio of rotational to gravitational potential energy after bounce reaches several percent, which is close to the onset of the low-T /W instability. The neutrino spheres of all flavors are located above the PNS core surface at R ∼ 10 km, where the low-T /W instability starts to (typically) develop, and also below the shock which is deformed by the spiral SASI (for model R1B00). Once the two instabilities appear, they can deform the neutrino spheres and potentially be the origin of the neutrino time modulation (see Kazeroni et al. (2017) for the possible connection of the two instabilities). However, we note that the smaller modulation in the ν x signals seem to favor that the outermost ν e sphere is more strongly affected by the spiral SASI.
Indeed there is a quantitative evidence that the deformation of neutrino sphere creates the time modulated neutrino signals. In Fig. 8, we plot spectra of the (angledependent) neutrino luminosity Lν e corresponding to the  at R = 400 km and θ = 90 • with N ν being the number luminosity estimated in the same way as the luminosity L ν (Tamborra et al. 2014). Although we here use the number luminosity N ν , we can do the same discussion using the luminosity L ν . R 11 is the isodensity surface extracted at the rest mass density of ρ = 10 11 g cm −3 corresponding roughly to the neutrino sphere. The normalized mode amplitude of spherical polar expansion of R 11 is evaluated exactly in the same manner as that of the shock surface. Here we focus on theν e signals (120 ≤ t pb ≤ 180 ms) bearing in mind the detectability (Abe 2016;Abbasi et al. 2011) (which will be reported elsewhere).
The black line in Fig. 8 shows that the time modulation seen in Fig. 7 peaks at F ∼ 125 Hz. This component is mainly composed of m = 1 neutrino number-flux as the two peaks of red and black lines are appearing nearly the same frequency. The peak of Nν e with m = 2, which is a daughter mode of m = 1, appears closely at a double frequency F ∼ 240 Hz of that of m = 1 as expected, but the m = 2 mode seems to contribute less to the total  Fig. 7, of normalized mode amplitudes of the number luminosity Nν e,m with the lower index m being the azimuthal mode either m = 1 (red) or m = 2 (blue), and of the normalized mode amplitude of spherical polar expansion of the isodensity surface R 11 , extracted at the rest mass density of ρ = 10 11 g cm −3 , for (ℓ, m) = (1, 1) (green). The vertical axis is in arbitrary unit. Here R 11 is roughly representing the neutrino sphere. The spectra are obtained by the Fourier transformation for the time interval of 120 ≤ t pb ≤ 180 ms.
neutrino signals than the m = 1 mode. Finally, as it is obvious from the peak at F ∼ 125 Hz in green line, the origin of these time modulations of the neutrino signals is m = 1 deformation of neutrino sphere represented by R 11 . We thus conclude that the strong spiral SASI appearing in R1B00 deforms the neutrino sphere with the same m = 1 mode and leads to the characteristic neutrino signals.
We also mention that we observe a clear north-south asymmetry in neutrino signals in model R1B12 for t pb 120 ms, i.e., between dash-dotted and dotted red lines, which cannot be seen in the corresponding lines of R1B00. In this model R1B12, the neutrino emission toward the north pole is significantly stronger than the one toward south. The excess toward north is consistent with the one-sided explosion to the south pole (see the red line in the top panel of Fig. 4 for A 10 mode). Due to the shock expansion mainly toward the south, the mass accretion is stronger in the northern hemisphere, which results in higher accretion luminosities and neutrino energies in the north pole.
3.5. The Role of Neutrino Heating Next we make a comparison of the energetics and discuss the role of neutrino heating among the models, particularly how the neutrinos contribute to the shock expansion. In Fig. 9, we plot the mass in the gain region M gain (top left), heating rateQ (top middle), gain and shock radii, R gain and R shock , respectively (top right), specific heating rate ζ =Q/M gain (bottom left), heating efficiency η =Q/(L νe +Lν e ) which measures how much of the emergent ν e andν e 's contribute to the matter heating (bottom middle), and ratio of advection to heating time scale τ adv /τ heat (bottom right) for each model. To obtain these values, we first define the gain radius R gain (θ, φ) at each radial direction (θ, φ). R gain is defined at the first point where the net energy deposition rateq becomes zero behind the shock, withq being defined bẏ Then each value is defined by and where the surface integral ds appearing in the denominator of Eq . (29) is performed in front of the shock surface and v r is the radial component of the three velocity v i . In the top right panel, we show spherical averaged shock (solid lines) and gain radii (dashed). While in the rest of panels, to illustrate how the values vary relative to the rotational axis, we divide the space into two equal volume regions, polar and equator, and show the values evaluated in each region. Here, we define the polar region (labeled by "Pol") by the cone angle of 60 • around the rotational axis, i.e., θ ≤ 60 • or θ ≥ 120 • , and the equatorial region (labeled by "Equ") by 60 • < θ < 120 • . These ranges are used in the volume and surface integrals in Eqs. (27)-(30). When we evaluate η(=Q/(L νe + Lν e )), ζ, and τ adv /τ heat , we first evaluate every quantity, e.g.Q and L νe + Lν e , in each region and then take their ratio. Regarding the model R0B00, we show its values integrated over all solid angles (labeled by "4π") since it has basically no significant angle dependence. Note that we show half values for extensive variables, i.e., M gain anḋ Q, for model R0B00 to compare with other models. Fig. 9 clearly shows how the rotational and magnetic field effects appear in general and also how they change the values relative to the rotational axis. The (spherically averaged) gain radius locates more inward in rotating models R1B12 (red dashed line in the top right panel) and R1B00 (blue dashed) than the non-rotating model R0B00 (black dashed). As can be seen in the top-left and -middle panels, the more inward R gain and larger R shock produce a more extended gain region and consequently a larger mass and total heating rate integrated over that region. The non-rotating model R0B00 shows smallest M gain andQ, typically several times smaller than the other two. The specific heating rate ζ (bottom left panel) also shows a rotational dependence. In general, R0B00 presents higher ζ(=Q/M gain ), although M gain andQ themselves are smaller than the other two. On the other hand, from the perspective of neutrino heating efficiency, η(=Q/(L νe + Lν e )) in R0B00 shows the least efficiency (bottom middle). Therefore rotation works to lower the specific heating rate ζ but raise the heating efficiency η. Such a trend is consistent with previous rotating models with detailed neutrino transport in Summa et al. (2018). The solid and dashed lines represent that the volume/surface integral is performed around the polar axis (labeled by "Pol") and equatorial plane (labeled by "Equ"), respectively. Regarding the non-rotating model R0B00, we integrate over all solid angles (solid black line labeled by "4π"). Note that we show half values for extensive variables, i.e., M gain andQ, for model R0B00 for comparison with other models.
In the bottom right panel, all these features mentioned above are aggregated in a value τ adv /τ heat . Higher τ adv /τ heat represents that the dwell time of matter in the gain region is relatively long in terms of heating time scale. It thus leads to a more favorable condition for the explosion. Particularly τ adv /τ heat larger than one can be a measurement of the onset of runaway shock expansion due to neutrino-heating (see, Müller et al. (2017); Summa et al. (2018); Ott et al. (2018) for the latest 3D successful explosion models and also O'Connor & Couch (2018) for the 3D non-explosion models). In the bottom right panel, model R1B12 which has the largest gain region shows highest τ adv /τ heat (red lines), while model R0B00 shows the lowest value (black line). Therefore our result also shows that rotation makes τ adv /τ heat higher. This tendency is again consistent with Summa et al. (2018). In addition, the magnetic fields also assist the expansion of the shock surface and produce higher τ adv /τ heat than the corresponding non-magnetized model R1B00.
Next we discuss how rotation affects the energetics in each region relative to the rotational axis. First, in the model R1B00, both M gain andQ show significantly higher values along the equator (blue dashed lines) than those in the polar region (blue solid). The blue dashed and solid lines start to diverge when the second shock expansion takes place at t pb ∼ 220 ms. The higher values seen in the equatorial region are again due to the rotational shock expansion. These rotational effects were already discussed by Nakamura et al. (2014), though with a very simplified neutrino light bulb method, and we obtain a consistent result in our self-consistent M1 neutrino transport simulations. The heating efficiency η in the equatorial region is also nearly twice as high as that in the polar region. As a consequence, τ adv /τ heat exceeds one, only in the equator (blue dashed line), and not in the polar region (blue solid line). If the neutrino heat-ing were more efficient and could actually aid the second shock expansion, it would directly lead to the shock runaway phase. The model R1B00, however, deflates and does not enter the runaway phase during our simulation time up to t pb ∼ 500 ms.
We see an interesting feature in model R1B12. In this rotating magnetized model, as we have explained in Sec. 3.2, it exhibits a rapid shock expansion toward the rotational axis soon after core bounce. Therefore τ adv /τ heat in the polar region (red solid line) shows slightly higher value than the equatorial one (red dashed). However, the higher value in the red solid line only persists during the first ∼ 100 ms after bounce and afterward the red dashed line takes over the solid one with largely exceeding one. Interestingly, τ adv /τ heat in the polar region shows basically less than unity till t pb ∼ 200 ms, although the shock runaway already occurs mainly toward the polar region. The trend is thus completely opposite to that of R1B00 in which the region with larger shock expansion exhibits larger τ adv /τ heat . We interpret these behaviors as that the neutrino heating is not the main mechanism of the bipolar shock expansion in R1B12, but the magnetic fields play the leading role to aid the shock expansion. On the other hand, as the red dashed line is exceeding unity, the shock expansion along the equator is mainly supported by neutrino heating.
3.6. The LESA Tamborra et al. (2014) reported the existence of the lepton-number emission self-sustained asymmetry termed LESA. This phenomenon is characterized by a spherical symmetry breaking of the lepton number emission, basically dominated by a dipole mode. In their paper, they explained the fundamental mechanism of LESA by a partial distribution of Y e on the PNS core surface r ∼ 25 km. In addition, such a partial distri-bution of Y e is an outcome of non-spherical mass accretion, basically with low mode ℓ = 1, on to the PNS core surface that replenishes the lepton-rich matter on a specific point. According to the schematic figure 15 in Tamborra et al. (2014), the non-spherical mass accretion is caused by shock deflection when the stellar mantle with rich Y e passes through the globally deformed shock surface. O' Connor & Couch (2018) and Vartanyan et al. (2019) also reported the appearance of LESA using M1 neu-trino transport method, although the relative magnitudes of the dipole mode are significantly weaker than the values in Tamborra et al. (2014) employing the ray-by-ray transport in the 3D models. Furthermore O'Connor & Couch (2018) pointed out the importance of velocity dependent terms in the neutrino transport as the models without that term do not show any conclusive evidence for LESA. Vartanyan et al. (2019) also showed that the dipole mode can be comparable to the monopole one in the late post bounce phase t pb ∼ 650 ms. Therefore, although the significance of dipole magnitude may actually depend on the detailed neutrino transport scheme, the LESA seems to be a common phenomenon in CCSNe.
While in the lower three panels, we plot the value 3L 1m /L 0 for each quantum number m in each model to discuss the correlation with the shock morphology. From top panel, we see that the absolute magnitude of normalized dipole mode in model R1B12 shows significantly larger value than the other two non-explosion models. Although the highest value ∼ 0.4 (red line in top panel) is significantly smaller than the values of Tamborra et al. (2014), in which they find the excess of dipole mode, interestingly in all models irrespective of the explosion, our value is comparable to those of O'Connor & Couch (2018) and Vartanyan et al. (2019). They reported the relative magnitudes of the dipole mode of ∼ 0.5 for t pb 400 ms. In this model R1B12, the dominant contribution to the total dipole mode is mainly coming from m = 0 mode (blue line in the second panel). Since it basically exhibits the positive value for t pb 100 ms, the relativeν e 's number flux is less toward positive z-axis and higher toward negative z-axis. From Fig. 4, the shock morphology with (ℓ, m) = (1, 0) mode becomes stronger for t pb 100 ms with negative value that reflects that the shock expansion takes place relatively stronger toward the negative z-axis (also see, the final snapshot of the shock morphology in Fig. 1). It is thus opposite to the dipole mode of the lepton number flux.
As first pointed out by Tamborra et al. (2014), the LESA is driven by a replenishment of rich Y e material onto one side of the PNS core that can be established once the one-sided explosion morphology (note in our case, which is magnetorotationally-driven) becomes almost steady. Our results first show that the LESA originally proposed in the non-rotating progenitor also develops in our 3D MHD exploding model.

DISCUSSION AND CONCLUSIONS
We have presented the first 3D-GR MHD simulations of a 20 M ⊙ star with spectral neutrino transport. For the nuclear EOS and neutrino opacities, we used SFHo of Steiner et al. (2013) and a baseline set of weak interactions (Bruenn 1985;Rampp & Janka 2002), where nucleon-nucleon bremsstrahlung is additionally taken into account, respectively. Neutrino transport is handled by M1 closure scheme with the red and Doppler shift terms being fully considered.
We calculated three models, non-rotating nonmagnetized, rotating non-magnetized, and rotating magnetized models to explore the effects of progenitor's rotation and magnetic field both on the dynamics and neutrino profiles. Regarding the dynamics, while no shock revival was observed in two non-magnetized models during our simulation times, the shock expansion initiated shortly after bounce in a rotating magnetized model. Initially the shock morphology takes a bipolar structure, which was eventually taken over by a unipolar one. The shock front reached 1000 km at t pb ∼ 220 ms and still continued expansion at the end of our simulation time. From our analysis for the rotating magnetized model, we interpreted that the polar expansion is driven mainly by the magnetic pressure, while the equatorial expansion is facilitated by the neutrino heating. Although we did not see the shock revival in two non-magnetized models, the standing shock locates further outward in the rotating model, which expands the gain region and increases the mass in the region. Therefore, we obtained a consistent result with previous studies that the (moderate) rotation makes the condition more favorable for the explosion than the non-rotating case.
Using the same (or very similar) non-rotating 20 M ⊙ progenitor star as in this study, some previous 3D studies have shown a successful explosion (Melson et al. 2015;Ott et al. 2018;Burrows et al. 2019), while the others have not (Tamborra et al. 2014;Melson et al. 2015;O'Connor & Couch 2018). It is thus worth comparing our non-rotating and non-exploding model R0B00 with these previous studies. One of major limitations in this work is its relatively lower numerical resolution compared to the previous ones. It has been thoroughly examined that insufficient resolution can potentially inhibit the shock revival due to less turbulent pressure (e.g., Couch & Ott 2015;Müller & Janka 2015;Roberts et al. 2016;Takiwaki et al. 2016;Burrows et al. 2019;Nagakura et al. 2019). For instance, Ott et al. (2018) performed full relativistic 3D calculations with M1 neutrino transport and obtained the shock revival. This might be possibly due to their higher numerical resolution within the shock surface that achieves a factor of ∼ 2 − 4 higher than ours. The higher numerical resolution allows the growth of turbulence leading to an additional pressure support. It should be also noted that more up-to-date neutrino opacities, e.g., a strangenessdependent contribution to the axial-vector coupling constant or many-body corrections to neutrino-nucleon scattering (Burrows & Sawyer 1998;Horowitz et al. 2017), generally benefit to facilitate the shock revival (e.g., Kotake et al. 2018;Burrows et al. 2019). We are currently conducting CCSN simulations with better neutrino opacities following Kotake et al. (2018), which would be reported elsewhere in the near future.
We investigated the effect of the precollapse rotation and magnetic fields on the neutrino signals. In general, both of the rotation and magnetic field decrease the neutrino luminosity and energy as they make the PNS core less compact due to the centrifugal force and/or mass ejection. In addition, the rotation produces angle dependent neutrino signals relative to the rotational axis. The neutrino luminosity and energy along the equator are significantly lower than those along the rotational axis. We observed a quasi-periodic time modulation in the neutrino signals especially in model R1B00 toward the equator that is greatly suppressed along the rotational axis. From our spectral analysis, the peak frequencies of the time modulated signals and of the m = 1 deformation of neutrino sphere(s) have nearly the same value. Therefore, together with the less modulation in heavier type neutrino signals, we consider that the spiral SASI mode deforms the neutrino spheres leading to the quasi-periodic signals. As previously identified (Tamborra et al. 2014;O'Connor & Couch 2018;Vartanyan et al. 2019), we also witnessed the LESA activity (albeit weak), which is found to be strongest for our MR-explosion model. Our results showed clear dependencies of neutrino signals on progenitor's rotation, magnetic field, and the observation angle. A more systematic study (such as changing the progenitor model, the initial magnetorotational strength, and the inclination between the rotation and magnetic axis) is needed for clarifying the multi-messenger signals from magnetorotationallydriven CCSNe.
As an important 3D effect, we showed that the kink instability is most likely to appear in the magnetized model that can potentially broaden the expanding blob, leading to weaker bipolar jets. However, the PNS core may also be subject to the low-T /W instability, we could not disentangle the outcomes of these two possibly coexisting instabilities. Further numerical simulations by other independent groups, preferably with finer numerical resolutions, are definitely required to clarify the interplay between the two instabilities.
In the end of our discussion, we briefly mention the possible role of MRI. Although the stellar magnetic field configuration and its strength at pre-collapse phase are poorly understood, the initial value employed in this study 10 12 G might be too strong according to magnetized stellar evolution calculations by Heger et al. (2005), which gives 10 9 G. To see how the MRI amplifies such plausibly weak magnetic fields, Obergaulinger et al. (2009) conducted local shearing disk simulations. Their results showed that the initial seed magnetic fields inside the PNS O(10 12 ) G can be amplified to dynamically relevant strengths O(10 15 ) G within several ms. Since the main magnetic field amplification mechanism during core-collapse is compression, their initial seed magnetic fields inside the PNS O(10 12 ) G could originate from the pre-collapse phase O(10 9 ) G, which seems compatible with the stellar evolution calculation. Sawai & Yamada (2014) has shown in their global 2D axisymmetric simulations that the MRI can not only amplify the initial seed magnetic fields but also produce a global magnetic field in the postshock region. Later, a globally ordered field amplification in the PNS was found in full 3D-GR MHD simulations by Mösta et al. (2014). All these facts indicate that model R1B12 in this study is not too extreme but might be plausible, although the typical length scale of the MRI O(10) m is far too small to resolve by our current numerical grid size (simply limited by our currently available computational resources). Other than the MR explosion scenario, the turbulence in the MRI could enhance the neutrino heating efficiency, which could impact the neutrino mechanism (Sawai & Yamada 2014;Masada et al. 2015). All these subjects remain to be studied. As such, we can see a vast untouched (research) territory lying in front of us, into which we have just made a first jump with a newly developed tool (our 3D-GR MHD code) in hand.
We thank Shota Shibagaki, Martin Obergaulinger, and Federico Maria Guercilena for helpful discussions and useful comments. This research was supported by the ERC Starting Grant EUROPIUM-677912 (TK and AA), JSPS KAKENHI Grant Number (JP17H05206, JP17K14306, and JP17H01130, JP17H06364, JP18H01212 (KK and TT)), and JICFuS as a priority issue to be tackled by using the Post 'K' Computer. Numerical computations were carried out on Cray XC50 at CfCA, NAOJ.

APPENDIX
In this appendix, we show that our metric evolution implementation has a fourth-order convergence in space by checking the well-known polarized Gowdy wave test (Alcubierre et al. 2004). We omit to write the Gowdy wave metric and initial condition that can be found elsewhere (e.g., Alcubierre et al. 2004). We evolve the collapsing Gowdy-wave metric backwards in time using the harmonic slicing condition with zero shift vector β i = 0 as for the gauge condition. Although the Gowdy-wave is a plane wave, we perform the test both in full 1D and 2D space. In the latter 2D case, we tilt the propagation direction of the plane wave at 45 • in the xy-plane. We employ two different grid spacing dx = 1/N with N = 64 or 128 to check for the numerical convergence. Fig. 11 shows the L 2 norm of violation of the Hamiltonian constraint |H| 2 for coarser spacing model with N = 64 (black line) and finer one with 128 (red). For finer resolution models (red lines), we multiply |H| 2 by 2 4 , since we use fourth-order spatial finite differencing. From the figure, we see that there is almost a perfect overlap during the first ∼ 180 and ∼ 40 crossing times in 1D and 2D test, respectively, which shows that our metric evolution scheme actually achieves a fourth-order convergence in space. Fig. 11.-We plot the L 2 norm of violation of the Hamiltonian constraint |H| 2 for coarser spacing model with N = 64 (black lines) and finer ones with 128 (red). For finer resolution models (red lines), we multiply |H| 2 by 2 4 , since we use fourth-order spatial finite differencing. We evolve the metric backward in time starting at T ∼ 9.875.