Two-dimensional numerical study for magnetic field dependence of neutrino-driven core-collapse supernova models

We study the effects of the magnetic field on the dynamics of non-rotating stellar cores by performing two-dimensional (2D), magnetohydrodynamics (MHD) simulations. To this end, we have updated our neutrino-radiation-hydrodynamics supernova code to include MHD employing a divergence cleaning method with both careful treatments of finite volume and area reconstructions. By changing the initial strength of the magnetic field, the evolution of $15.0$, $18.4$ and $27.0$ $M_\odot$ presupernova progenitors is investigated. An intriguing finding in our study is that the neutrino-driven explosion occurs regardless of the strength of the initial magnetic field. For the 2D models presented in this work, the neutrino heating is the main driver for the explosion, whereas the magnetic field secondary contributes to the pre-explosion dynamics. Our results show that the strong magnetic field weakens the growth of the neutrino-driven turbulence in the small scale compared to the weak magnetic field. This results in the slower increase of the turbulent kinetic energy in the postshock region, leading to the slightly delayed onset of the shock revival for models with the stronger initial magnetic field.


INTRODUCTION
Core-collapse supernovae (CCSNe) are one of the most energetic explosions in the universe, marking the catastrophic end of massive stars. Extensive investigations over the decades (e.g. Colgate & White 1966, see also Mezzacappa 2005; Janka 2012; Kotake et al. 2012;Burrows 2013;Foglizzo et al. 2015;Müller 2020 for reviews) have shown that the most promising way to explode massive stars ( 8M ) is the neutrino mechanism (Bethe 1990). In this mechanism, neutrinos from the protoneutron star (PNS) heat the material behind the stalled shock, leading to the shock revival into explosion. The neutrino heating efficiency is significantly enhanced by non-radial flows triggered by various hydrodynamic instabilities including neutrino-driven/PNS convection and the standing accretion shock instability (SASI; Blondin et al. 2003;Foglizzo et al. 2006). This has been confirmed in a growing number of self-consistent CCSN simulations (e.g. Hanke et al. 2013;Takiwaki et al. 2014;Lentz et al. 2015;Müller 2015;Summa et al. 2016;Takiwaki et al. 2016;O'Connor & Couch 2018;Pan et al. 2018;Ott et al. 2018;Kuroda et al. 2018;Nagakura et al. 2019a;Vartanyan et al. 2019;Nakamura et al. 2019;Melson et al. 2020), some of which could closely account for canonical CCSNe with the explosion energies of the order of 10 51 erg (≡ 1 Bethe, 1 B in short) or less.
However, there are some very energetic subclasses of supernovae, which is highly unlikely to be explained by the conventional neutrino mechanism. Those events so far observed include hypernovae (e.g. Iwamoto et al. 1998;Soderberg et al. 2006) and superluminous SNe (SLSNe; see, e.g. Gal-Yam 2012;Nicholl et al. 2013;Moriya et al. 2018, for reviews). The most plausible scenario to account for these extreme events requires additional energy injection via the magnetohydrodynamically-driven (MHD, in short) explosions (e.g. Wheeler et al. 2002;Burrows et al. 2007;Dessart et al. 2008Dessart et al. , 2012. The kinetic energy of hypernovae exceeds ∼ 10 B, which is ten times larger than that of the canonical CCSNe 1 . Concerning the SLSNe, the luminosity is 10-100 times higher than that of the typical CCSNe. To explain the excess of the bright luminosity, several scenarios have been proposed including the interaction scenario between the SN ejecta and its dense circumstellar medium (Chevalier & Irwin 2011;Moriya et al. 2013;Sorokina et al. 2016) or the pair-instability SN scenario (Barkat et al. 1967;Rakavy & Shaviv 1967;Terreran et al. 2017). Besides, the injection of the additional energy into the SN ejecta by the central engine of a rapidly rotating proto-magnetar is also a promising scenario for the excess of luminosity (Kasen & Bildsten 2010;Woosley 2010;Wang et al. 2015;Chen et al. 2016). For the two classes of the extreme events mentioned above, one could speculate that a common solution requires the strong (magnetar-class), large-scale magnetic fields, which can directly couple the newly-born PNS (or proto-magnetar) to its surroundings via various MHD processes (Usov 1992;Thompson 1994;Bucciantini et al. 2009;Metzger et al. 2011).
Assuming the magnetic flux conservation, the strong surface magnetic field (∼ 1 kG) of OB-type stars (e.g. Donati et al. 2002Donati et al. , 2006Hubrig et al. 2006) is considered as a possible candidate to account for the magnetic field of magnetars. On the other hand, the majority of massive stars possess the weak magnetic field, which is supported by both the observation ( 100 G; Wade & MiMeS Collaboration 2015) and stellar evolution calculations (Heger et al. 2005). In the latter case, enough amplification of the magnetic field during and after the collapse of the massive star is necessary to facilitate the MHD mechanism. Since the magnetic field amplification due to the field wrapping as the consequence of rapid and strong differential rotation scales linearly in time and takes a long time compared to the dynamical time scale of the system, the drastic and efficient field amplification mechanism is required when the seed field of the massive star is adequately weak. The exponential growth of the magnetic field is beneficial to amplify the weak seed field of the stellar core to the dynamically relevant level. The magnetic field amplification due to the magnetorotational instability (MRI; Balbus & Hawley 1991) is a candidate for the efficient field amplification mechanism in the rotating stellar cores (Akiyama et al. 2003;Masada et al. 2006Masada et al. , 2007Masada et al. , 2012Masada et al. , 2015Obergaulinger et al. 2009;Sawai et al. 2013b;Sawai & Yamada 2014, 2016Mösta et al. 2015;Rembiasz et al. 2016;Reboul-Salze et al. 2020).
On the other hand, stellar evolution calculations pointed out that the majority of the magnetic core of the massive star is expected to be rotating slowly at the pre-collapse stage (Heger et al. 2005;Ott et al. 2006;Langer 2012). Especially, the stronger magnetic field would lead to the more efficient angular momentum loss of the stellar core by the magnetic braking even if some of the stars rapidly rotate initially (Ramírez-Agudelo et al. 2013) or the stars experience spin-up due to stellar mergers (Chatzopoulos et al. 2020). Asteroseismology of lowmass stars also suggested that an unmodeled, more efficient angular momentum transport process is necessary to explain the spin period of the cores Fuller et al. 2014). Observations of surface rotational velocities of B-type stars with strong magnetic fields also favor slow rotators (Shultz et al. 2018). Even in such slowly rotating progenitors, it has been pointed out that the precollapse magnetic field, if sufficiently strong, could affect the explosion dynamics. In the context of the slowly-and non-rotating progenitor, Endeve et al. (2010Endeve et al. ( , 2012 and Obergaulinger et al. (2014) studied the dynamics of the MHD core-collapse and the exponential amplification of the magnetic field due to the SASI and convection. More recently, Müller & Varma (2020) have addressed the role of the magnetic field in the neutrino-driven explosion by performing three-dimensional (3D) MHD simulations of a slowly rotating progenitor of a 15M star.
Joining in these efforts to study the impact of the magnetic field on both the extreme and ordinary explosions of the massive stars, we investigate the magnetic field dependence of neutrino-driven explosion in the non-rotating cores by performing two-dimensional (2D), axisymmetric, MHD core-collapse simulations for several representative progenitors. In the context of 2D simulations of the non-rotating and magnetized cores, Obergaulinger et al. (2014) were the first to point out that magnetic pressure support in the gain region (via turbulence) fosters the onset of neutrino-driven explosion. This result clearly presented evidence that implementation of appropriate neutrino transport is needed for a quantitative study of MHD CCSN modeling. However, only electron and anti-electron neutrinos were taken into account in Obergaulinger et al. (2014) at that time (see, however, Obergaulinger & Aloy 2020). To revisit the problem, we have updated our supernova code (3DnSNe) to include MHD by implementing a divergence cleaning method (Dedner et al. 2002) with both finite volume and area reconstructions based on Mignone (2014). Our base code deals with three-flavor neutrino transport (namely, ν e ,ν e , ν X with ν X denoting the heavy-lepton neutrinos) ) based on the Isotropic Diffusion Source Approximation (IDSA) scheme (Liebendörfer et al. 2009), in which a detailed code comparison was already made in spherically symmetric simulations ) and in 2D simulations ) using a widely used 20 M star of Woosley & Heger (2007). This paper is organized as follows: In Section 2, the numerical methods and models are described. Our numerical results of the MHD core-collapse of non-rotating stellar cores in axisymmetry are presented in Section 3. In Section 4, we discuss the field configuration of the proto-magnetar based on our simulation results. Finally, we summarize and discuss our findings in Section 5.

NUMERICAL METHODS AND MODELS
We have updated our supernova code, 3DnSNe 2 , that is designed for CCSN simulations in a 3D spherical coordinate system to the latest version. In this work, the code is now extended to an MHD code from a hydrodynamic (HD) one with spectral neutrino transport that is solved by the IDSA scheme (Liebendörfer et al. 2009). We have updated the original (two-flavor, i.e. ν e ,ν e ) IDSA scheme in several manners, such that the evolution of the streaming neutrinos is self-consistently solved (Takiwaki et al. 2014) and that three-flavor neutrino transport is solved including approximate general relativistic corrections (e.g. Kotake et al. 2018 for more details). A detailed code comparison has been performed in O'Connor et al. (2018) with one-dimensional (1D) geometry. The 3DnSNe code has been used in the following works: Cherry et al. (2020); Zaizen et al. (2020); Sasaki et al. (2020); Nakamura et al. (2019); Sasaki et al. (2017); Sotani & Takiwaki (2016, 2020; Nakamura et al. (2015).
In the latest code, we solve the ideal MHD equations in the spherical coordinate system (r, θ, φ). The governing equations are ∂ρ ∂t where ρ, v, B, P t , e and Φ are the mass density, the fluid velocity vector, the magnetic field vector, the total (thermal and magnetic) pressure, the total energy density and the gravitational potential, respectively. Y l is the lepton fraction and the subscription l denotes the species of leptons: l = e, ν e ,ν e , ν X and Z m is the specific internal energy of the trapped neutrinos and m represents the species of neutrinos: m = ν e ,ν e , ν X . Q E , Q m are the change of the energy and Γ l is the change of number fraction due to the interaction with the fluid and neutrinos. I is the unit matrix. The explicit expressions of the governing equations in the spherical polar coordinates are given in Appendix A. As for the approximate Riemann solver, the HLLD scheme (Miyoshi & Kusano 2005) is newly implemented in our code to solve equations (1)-(4) in a conservative form (see Appendix A for the treatment of equations 6 and 7). To get around the Carbuncle instability, we switch from the HLLD scheme to the HLLE scheme (Einfeldt 1988) in the vicinity of strong shocks (see Kim et al. 2003 and the references therein). Equation (5) and ψ in the induction equation (4) are related to the divergence cleaning method proposed by (Dedner et al. 2002). This method reduces numerical errors of the solenoidal property of the magnetic field within minimal levels. The average relative divergence error estimated by Error3(B) ave proposed by Zhang & Feng (2016) is less than 1% in our work. c h and c p characterize the propagation speed and the damping rate of the numerical divergence of the magnetic field, respectively (see Appendix B for details). Equation (5) is solved using HLLE scheme. To retain the total energy including gravitational binding energy, we use the method of Müller et al. (2010) in solving equation (3). In equation (8), the spherically symmetric gravitational potential is taken in the form of phenomenological general relativistic potential of Case A in Marek et al. (2006) and the multi-pole components are added following Wongwathanarat et al. (2010).
We use a finite volume method (Li & Li 2003;Mignone 2014) to solve conservation equations basically. However, since the divergence cleaning method is employed in our code, both the finite volume and area methods are used to solve the induction equation. The details of this special treatment, in addition to the reconstructions of the physical variables for the second-order accuracy in space, are described in Appendix B and C. A code test is given in Appendix D.
Our setups for microphysics are similar to those of O' Connor et al. (2018). The adopted neutrino reaction rate is set5a of Kotake et al. (2018), i.e. the weak magnetism and recoil correction (Horowitz 2002) as well as nucleonâȂŞnucleon bremsstrahlung is added to the standard opacity set of Bruenn (1985). In this run, 20 energy groups that logarithmically spread from 1 to 300 MeV are employed. We use the equation of state (EOS) by Lattimer & Swesty (1991) (incompressibility K = 220 MeV).
We employ the non-rotating presupernova progenitors of 15.0, 18.4 and 27.0 M of Woosley et al. (2002). As for the initial configuration of the magnetic fields, we assume a simple topology following Suwa et al. (2007); Takiwaki et al. (2014); Obergaulinger et al. (2014). The magnetic field is given by a vector potential in the φ-direction of the form where r 0 = 1000 km characterizes the topology of the field. The magnetic field is uniform when the radius, r, is smaller than r 0 , while it is like dipole field when r is larger than r 0 . B 0 determines the strength of the magnetic field inside the core (r < r 0 ). In this study, we set B 0 = 10 10 , 10 11 or 10 12 G. The model name is labelled as 's27.0B10', which represents the 27.0 M model with B 0 = 10 10 G. We choose s27.0B10 as a fiducial model because 2D (albeit, non-magnetized) results using this progenitor are available in the literature (e.g. Hanke et al. 2013;Summa et al. 2016). We follow the dynamics up to t fin ∼ 400 − 500 ms after bounce, depending on the progenitor models. In most of the models, we terminate the simulations at the final time seeing that the diagnostic explosion energies are greater than 10 50 erg. We leave the more long-term simulation for future work. The calculations are performed in axisymmetry. Therefore, the derivatives with respect to the φ-direction (i.e. ∂ ∂φ ) are taken to be zero in the governing equations when we run 2D simulations. The grid spacing in this work is the similar to that of 2D runs in Takiwaki et al. (2014). In the radial direction, a logarithmically stretched grid is adopted for 480 zones that cover from the center up to 5000 km, whereas the polar angle in the θ-direction is uniformly divided into ∆θ = π/128. The innermost 10 km are computed in spherical symmetry to avoid excessive time-step limitations. Reflective boundary conditions are imposed on the inner radial boundary (r = 0), while fixed-boundary conditions are adopted for the outer radial boundary (r = 5000 km) except the gravitational potential that is inversely proportional to the radius at outer ghost cells. A reflecting boundary condition is imposed on the 2D symmetry axis (e.g. the z-axis in our 2D run). A numerical resolution test is given in Appendix E.

RESULTS
We first describe overall evolution of the magnetized and non-rotating stellar core for our fiducial model (s27.0B10) in Section 3.1. Then in the subsequent sections, we move on to present results focusing on the impact of the initial magnetic field strength on the postbounce evolution. The progenitor dependence of the shock evolution is presented in Section 3.4.
3.1 Overall evolution of non-rotating and magnetized core-collapse model of a 27M star Fig. 1 shows the temporal evolution of the spatial distribution of the entropy per baryon and magnetic field for the fiducial model (s27.0B10). The 2D color map of the entropy per baryon is illustrated in the negative region of x (x < 0). The structure of magnetic field lines is drawn by a line integral convolution method (Cabral & Leedom 1993) in the positive region of x (x > 0). The color depicts the strength of the magnetic field. Panel (a), (b), (c) and (d) correspond to the time t pb = 100, 200, 300 and 500 ms after bounce, respectively. Hereafter t pb denotes the postbounce time.
The core bounce occurs after ∼ 200 ms (i.e. t pb = 0) after the start of the simulation, leading to the shock formation at the radius of ∼ 20 km. The bounce shock stalls at r ∼ 140 km around t pb = 100 ms, and then turns into the standing shock (see also, the top left panel of Fig. 2). When the shock stalls, the structure of the magnetic field lines is like a split monopole as shown in the right-half panel of Fig. 1a. Before the shock stall (t pb 100 ms), the flow is almost restricted in radial direction. The split-monopole like configuration is made because the magnetic field is "frozen-in" with respect to the matter motion. The electric resistivity of the magnetic field is so small that it is disregarded in this work, which can be well justified in the CCSN environment (Sawai et al. 2013a). The initial vector potential (equation 9) gives magnetic loops on the equatorial region at around r ∼ 1000 km. These magnetic loops also gravitationally collapse (dragged by matter infall) and are shown on the equatorial plane (x 30 km and z = 0) in Fig. 1a. The center of loops is located at around x ∼ 45 km and seen as a small blueish region.
As the (maximum) shock radius starts to gradually shrink after t pb 100 ms (e.g. Fig. 2a), it gradually deviates from the shock trajectory of the corresponding 1D model (black solid line in Fig. 2a). This marks the growth of non-spherical motions in the postshock region. One can clearly observe the deformation of the shock in the left-half panel of Fig. 1b at t pb = 200 ms. In Fig. 1b, one can also see the penetration of the magnetic field lines (thin red curves in the right-half panel) into the postshock region (high entropy region in the left-half panel), which makes the field configuration much more complicated than that outside the shock. In our ideal MHD simulations, the field amplification in the postshock region occurs due to compression and stretching of the magnetic field, which is governed by the non-radial matter motions. Note in our 2D models that we do not attempt to differentiate the origin of the "non-radial" motions either originating predominantly from the SASI or neutrino-driven convection because the SASI is liable to be overestimated in 2D compared to 3D simulations (e.g. Hanke et al. 2012Hanke et al. , 2013Fernández et al. 2014). Fig.1c shows a snapshot after the shock revival (t pb = 300 ms, see also Fig. 2a). The low-mode deformation of the shock and the formation of the high entropy region (colored by red in the entropy plot) is a common feature of 2D neutrino-driven explosion models. The highly aspherical shock (Fig. 1d) continues to propagate strongly toward the south pole up to a radius of ∼ 1000 km at the final calculation time (t fin ∼ 500 ms) for this model. The magnetic field configuration is similar to the blast morphology as in the previous snapshots.  In Fig. 2a, a black line, as already mentioned, denotes the shock evolution of the 1D HD model of the 27.0 M star, which is shown as a reference. The shock radius in 1D maximally reaches to ∼ 140 km at t pb ∼ 100 ms, then continues to contract during the simulation time. On the other hand, the shock revival occurs in the 2D MHD models regardless of the strength of the initial magnetic field. Looking at the panel very carefully, one might notice an interesting tendency about the slight delay of the onset of the shock revival between the three MHD models.

Neutrino-driven explosion of MHD core collapse
To show this clearly, we make a comparison in Fig. 2b of the shock evolution only close to the shock revival time (between t pb = 150 ms and t pb = 250 ms). The shock revival time is indeed delayed for the strong initial field model (blue curve) comparing to the weak initial field model (red curve). Fig. 2c shows the evolution of the timescale ratio, τ adv /τ heat . Following Summa et al. (2016), we estimate the advection timescale as   where M g is the mass enclosed in the gain layer (gain mass) andṀ is the mass-accretion rate through the shock. The neutrino-heating timescale is defined by where |E tot,g | is the total energy of the material in the gain layer andQ heat is the neutrino-heating rate in this region. Since the residency time of matter in the gain region is related to the exposure of the material to neutrino heating, τ adv /τ heat 1 is a necessary condition for the onset of the shock revival (e.g. Buras et al. 2006). As shown in Fig. 2c, τ adv /τ heat rises toward unity rapidly at around t pb ∼ 200 ms. It is noted that this shock revival timescale is consistent with that by Hanke et al. (2013) who conducted the 2D model using the same progenitor (27M star of Woosley et al. 2002) with more elaborate neutrino transport scheme. In Fig. 2c, it is important to point out that the growth rate of the timescale ratio (before exceeding unity) is highest for the weakly magnetized model (red line), which is followed in order by the moderately magnetized model (green line) and the strongly magnetized model (blue line). This feature is closely linked to the shock evolution after t pb ∼ 200 ms, namely, the onset of the shock revival and the subsequent runaway shock expansion is delayed for the strongly magnetized models. This indicates that the neutrino heating mainly contributes to the runaway shock expansion as shown in Fig. 2a. And the magnetic field secondary affects the shock evolution. The same tendency is also observed in Fig. 2d, which compares the evolution of the gain mass at around t pb ∼ 200 ms. Fig. 3 compares the time evolution of (a) neutrino luminosity and (b) mean energy of neutrinos between the model s27.0 series. Solid, dashed and dash-dot lines correspond to ν e ,ν e and ν X , respectively. Red, green and blue lines are the cases for B 0 = 10 10 , 10 11 and 10 12 G, respectively. Given the results mentioned above, it may not be so surprising that the initial magnetic fields have little impact on the luminosities and mean energies. It is also noted that they are in good agreement with those in Summa et al. (2016)   systematic 2D simulations covering a wide range of the (non-magnetized) progenitor models. Regardless of the difference in the neutrino transport scheme, the peak of the ν e /ν e luminosity (for the same progenitor) is ∼ 6 × 10 52 erg/s at around t pb ∼ 100 ms, which is consistent with our results. After the onset of shock revival (t pb 200 ms), the mean neutrino energy is in the range of 12 ∼ 14 MeV for ν e and 15 ∼ 17 MeV forν e in our model, which nicely matches with Summa et al. (2016), although ourν e mean energy is ∼ 8% higher than Summa et al. (2016) at our final simulation time.
Although there is no significant impact of the initial magnetic field strength on the neutrino properties (Fig. 3), we did witness the difference in the evolution of τ adv /τ heat (Fig. 2c) for models with stronger initial magnetic fields. This suggests that the stronger initial field affects the advection timescale predominantly than the neutrino heating timescale. In what follows, we explore how the strong initial field could affect the development of the non-radial matter motions in the postshock region, leading to the delayed onset of the shock revival.
One may wonder whether the amplified magnetic field in the postshock region could assist the explosion as reported in Obergaulinger et al. (2014). The plasma β, the ratio of the thermal pressure to the magnetic pressure, is a good indicator to check this possibility. Fig. 4 is a comparison of 2D spatial distribution of the plasma β at t pb = 200 ms (e.g. close to the shock revival time, see Fig. 2b) for the 27.0 M models with different initial magnetic fields. It is shown that the plasma β behind the shock is much larger than unity ( 10 1 ∼ 10 5 ), meaning that the shock revival in our models is predominantly driven by neutrino heating. This is absolutely not the case in the MHD explosion in the context of rapidly rotating and strongly magnetized cores (e.g. Kuroda et al. 2020), where the plasma β of ∼ unity is often achieved. Furthermore it is important to note that the neutrino-driven shock revival is obtained at t pb ∼ 200 ms for our weakly magnetized model (s27.0B10). This is in stark contrast with Obergaulinger et al. (2014)   . In order to focus on the convective motion, only θ-component of velocity is taken into account to compute the kinetic energy. The energy spectrum is defined by equation (12). Red, green and blue lines in each panel are the cases for B 0 = 10 10 , 10 11 and 10 12 G, respectively.
in the work, though assuming the same magnetic field strength (10 10 G). More detailed comparison between our models and the previous studies is presented in Section 3.4.

Magnetic field dependence on non-radial motion and development of turbulence
In the previous subsections, we have shown that the neutrino heating plays a dominant role in triggering the explosion of our 27 M models, whereas the magnetic field plays a secondary role, namely, to delay the explosion onset. We now proceed to clarify the reason by focusing on the role of the magnetic field on the non-radial matter motions and turbulence in the postshock region. Similar to Fig. 2, but Fig. 5a, 5b and 5c show the time evolution of the lateral kinetic energy (top panels), and the advection timescale in the gain region for the model series of s27.0, respectively. In Fig. 5a, the contribution from the lateral (θ)-component of the velocity (v θ ) is taken into account as a measure to quantify the vigor of the non-radial motions and turbulence in the gain region. The entire evolution of the lateral kinetic energy is shown in Fig. 5a (log-scale in the y-axis), whereas Fig. 5b focuses on the time around the shock revival (t pb ∼ 200 ms) (linear scale in the y-axis). From Fig. 5a, one can see that the lateral kinetic energies firstly increase exponentially before the shock revival (t pb ∼ 200 ms), and then reach asymptotically to ∼ 10 50 erg toward the final simulation time regardless of the different initial field strength.
Looking more closely at the linear phase (t pb 200 ms), Fig. 5b depicts that the growth of the kinetic energy is fastest (biggest) for model s27.0B10 (red line) compared to the more strongly magnetized models (green and blue lines). This feature, as previously identified in the 3D MHD simulations of Endeve et al. (2012) (but with more idealized setting), is also consistent with the earlier onset of the shock revival as seen in Fig. 2b.
Since the neutrino heating timescale is similar among the s27.0 models (e.g. Section 3.2), the difference of the advection timescale in the gain region should play a key role of the explosion onset, i.e. the longer the better. As shown in Fig. 5c, the advection timescale of the weakly  Fig. 5c. Likewise, these results are in favour of explaining the delayed onset of the shock revival of the strongly magnetized models. In order to investigate the role of the initial magnetic field on the turbulent activity in the postshock region, we compute the turbulent energy spectrum, e kine (l). Following Hanke et al. (2012), it is defined as, where Y lm is the spherical harmonics of degree l and m, and Ω is a solid angle. Note in our 2D simulations, m = 0 is only considered in equation (12). The turbulent energy spectra in Fig. 5d are evaluated at a fixed radius (r = 100 km) in the postshock region and at around the explosion onset. They are time averages at the center of t pb = 195 ms. Fig. 5d clearly shows that the energy density of higher-order modes (l 10) is bigger for the weaker magnetic field case (red line) than that for the stronger field case (blue line), although the energy density of the smaller-order modes is comparable in all the three cases. This implies that the strong magnetic field, most likely due to the magnetic tension, prevents the growth of the turbulent motions down to small scales (at larger l). The suppression of the turbulent energy at larger l is reconciled with the slow increase of the non-radial kinetic energy for the strongly magnetized model as shown in Fig. 5b. Again this is consistent with the delayed onset of the shock revival. Similar to Fig. 5a and 5d, but Fig. 6 shows the evolution of the lateral magnetic energy in the gain region (left panel) and the corresponding time-averaged energy spectra of the magnetic turbulence at t pb = 195 ms (right panel). Fig. 6a shows that the lateral magnetic energy is exponentially amplified up to the explosion onset (t pb ∼ 200 ms) in all the three models. The exponential growth terminates when the shock revival initiates (t pb 200 ms). In the non-linear phase (t pb 200 ms), the lateral magnetic energies are shown to be almost kept constant with time, whose strength is bigger for the strongly magnetized model. It is noted that the saturated values of model s27.0B12 (∼ 10 48 erg, blue line), s27.0B11 (∼ 10 46 erg, green line), and s27.0B10 (∼ 10 44 erg, green line) differ by the two orders-of-magnitudes, respectively, which is proportional to the square of the initial magnetic energy (i.e. ∝ B 2 0 ). Since the magnetic pressure/force does not play a crucial role in the shock revival (see Fig. 4 and the explanation in the last paragraph of Section 3.2), the magnetic field can be amplified by the compression and stretching due to the non-radial fluid motions. Actually, by comparing Fig. 5a and Fig. 6a, one can see that the magnetic energy inside the gain region is less than the kinetic energy in all the three models. The ratio of the magnetic energy to the kinetic energy, E mag (B θ )/E kine (v θ ), in models s27.0B10, s27.0B11 and s27.0B12 at around t pb ∼ 200 ms are 10 −5 , 10 −3 and 10 −1 , respectively, whereas the kinetic energy at the time is almost 10 49 erg similar to all the three models. The shock revival occurs before the magnetic energy in the gain layer reaches to the equipartition with the kinetic energy. We speculate if the explosion were much more delayed like in Obergaulinger et al. (2014), the magnetic field amplification could continue until it becomes sufficiently high (such as the equipartition level) enough to assist the shock revival.
To quantify the vigor of turbulence of the magnetic field, we estimate the spectrum of the lateral magnetic field as in equation (12). Plotted in Fig. 6b is In order to normalize the spectra with the different initial field strength, the spectra for the case with B 0 = 10 10 and 10 11 G are multiplied by a factor of 10 4 and 10 2 , respectively. This is reasonable because the magnetic energy is proportional to B 2 and that of B 0 = 10 10 and 10 11 G are 10 4 and 10 2 times smaller than that of B 0 = 10 12 G, respectively. From Fig. 6b, one can see that the energy spectrum of the magnetic turbulence typically decreases with l, i.e. the magnetic energy is mainly stored in a large-scale field. The slope is, however, shallower for the weak field model (s27.0B10, red line) comparing with the strong field models (green and blue lines). This indicates that a small scale (turbulent) magnetic field can be more preferentially developed for the weak field model. The relative excess of the lateral (turbulent) magnetic energy at larger l is consistent with the excess of the lateral (turbulent) kinetic energy as shown in Fig. 5d. These results suggest that the strong initial field could act to suppress the magnetic turbulence in the small scale compared to that for the weak initial field model.

Impact of different progenitor models
In order to investigate the impact of the initial magnetic field on the different progenitor models, we present results for the 15.0 and 18.4 M models. Similar to Figs. 2a and 2b, but Fig. 7 shows the evolution of the shock for model 15.0 M (top panels) and model 18.4 M (bottom panels), respectively. In both of the models, the shock revival occurs at t pb ∼ 200 ms (Fig. 7) regardless of the difference in the initial magnetic field strength. Remarkably, the slight delay of the shock revival for the strongly magnetized models is also obtained (s15.0B12 vs. s15.0B10 and s18.4B12 vs. s18.4B10). These features are common to those obtained in the 27M models as already mentioned in the last section.
The 2D HD simulation of the 18.4 M progenitor model was reported in Summa et al. (2016). They obtained the shock revival at t pb ∼ 520 ms, which is ∼ 320 ms later compared to our counterpart model (s18.4B10). Given the big difference of neutrino transport scheme between the two codes, we cannot unambiguously specify the reason of the discrepancy. But already in the 1D comparison work of O'Connor et al. 2018 (e.g. green line of their Fig. 2), we can see that our HD code (3DnSNe) leads to a larger shock radius especially later than t pb ∼ 200 ms relative to the other codes, leading to the enhanced heating rate in the gain region t pb 250 ms (see, green line of their Fig. 5). This could be one of the reasons of the early shock revival seen in our 18.4 M run. Be that as it may, the shock revival time and the neutrino properties (Figs. 2 and 3) show a good agreement with those of the 27M model of Summa et al. (2016). The match may be simply by chance, and a detailed comparison of 2D CCSN models from different groups is apparently needed to clarify these features. By running the ALCAR code for a given progenitor model (s20 of Woosley & Heger 2007) but with varying the neutrino opacities and the transport schemes, Just et al. (2018) showed in their detailed, 2D systematic simulations that a simplified treatment in the neutrino-pair processes makes the onset time of the shock revival significantly earlier (∼ 100 − 400 ms) compared to the case without such simplification (compare the explosion time of models s20-rbr and s20-rbr-pp{1/2/3} in their Table 1). Since we employ the simplified prescription (e.g. the assumption of the isotropy and local-thermal-equilibrium with respect to the heavy-lepton neutrinos), this might be also the reason of our earlier onset of the explosion.
As already mentioned before, the 2D MHD simulations of the non-rotating cores of the 15.0 M star (Woosley et al. 2002) were reported in Obergaulinger et al. (2014) with varying the initial field strength (like in this study). In their strong initial field model (10 12 G), they obtained a magnetically assisted explosion at t pb ∼ 600 ms, which is about ∼ 200 ms earlier than their weak field model (10 10 G). If the neutrino-driven shock revival is so much delayed, they pointed out that an equipartition between the turbulent kinetic and magnetic energy was archived due to the growing turbulence over the long period of time. This may seem to contradict with our results, i.e. the slight delay of the explosion onset for the strongly magnetized model (s15.0B12) as shown in Fig. 7a and Fig. 7b. However, in our models, the shock revival occurs much earlier until the equipartition could be achieved. The difference of the 2D hydrodynamics, which could stem from the omission of heavy-lepton neutrinos in Obergaulinger et al. (2014) and also from the difference of the neutrino transport scheme (M1 vs. ray-by-ray, e.g. Just et al. 2018), could cause these discrepancies. However, we think that both of the results are important in the sense that they illuminate two faces of the magnetic fields in the explosion dynamics of non-rotating CCSNe, namely, either in the case of relatively earlier onset of explosion as pointed out in our study or in the later onset of explosion as studied by Obergaulinger et al. (2014). 3

MAGNETIC FIELD OF PROTO-NEUTRON STAR
Shedding light on the dynamics of the MHD core collapse is a main purpose of this paper. On the other hand, the origin of the magnetic field of the PNS is also an interesting topic. Especially, the origin of the strong magnetic field of the magnetar is still under debate over the three decades (e.g., Duncan & Thompson 1992;Paczynski 1992;Thompson & Duncan 1993;Kouveliotou et al. 1998, Lai 2001Kaspi & Beloborodov 2017 for a review). Two possible formation scenarios have been proposed to account for the strong magnetic field of the magnetars. These include a turbulent dynamo amplification in a rapidly rotating PNS (Thompson & Duncan 1993) and a fossil field hypothesis (Ferrario & Wickramasinghe 2006;Vink & Kuiper 2006). In the latter scenario, the magnetic field of the progenitor star is amplified mainly due to the magnetic flux conservation as a result of the gravitational collapse of the massive star.
In the present study, we have assumed the non-rotating stellar cores. Based on the fossil field hypothesis, one can exploratory estimate the strength of the magnetic field inside the PNS. Since the initial magnetic field is given by equation (9) and uniform (B = B 0 ) inside the core (r < 10 3 km), the magnetic field strength of the PNS is estimated as follows: where r PNS denotes the radius of the PNS. From this rough estimate, one can anticipate that a magnetar-class magnetic field could be formed in the vicinity of the PNS. Note that the radius of the PNS is defined by the iso-density surface of 10 11 g/cm 3 . Fig. 8 compares the 2D distribution of the magnetic field strength for our fiducial runs (s27.0) with B 0 = 10 10 G (panel a), B 0 = 10 11 (panel b) and B 0 = 10 12 G (panel c), respectively. One can see that the magnetic field configuration in the central region (r 50 km) looks similar between the three panels. All the three panels reveal a common shortcoming of our 2D models, namely, the magnetic field is strongest along the symmetry axis (the z-axis) because of the reflective boundary condition there. In addition, an artificial structure of the magnetic field is also observed around the symmetry axis in the innermost region (r < 10 km) where the calculations are performed in spherical symmetry (see the last paragraph of Section 2). Having mentioned the shortcoming, let us focus on model s27.0B12 (panel c).
From panel (c), one can see that the field strength inside the region of r 30 km is typically ∼ 10 15 G (shown in red), which is consistent with the rough estimate of equation (14) in the PNS interior. The field strength becomes smaller below the PNS surface (30 ∼ 25 km r 12 km), which is shown in green or yellowish region in the panel. This region is convectively unstable due to the negative lepton gradient, whereas the region above is convective stable due to the positive entropy gradient. In the convective region, it is known in 2D simulations that the magnetic field is expelled due to the so-called convective flux expulsion (e.g. Galloway & Weiss 1981;Davidson 2001). This not only explains the reduction of the field strength there, but also the accumulation of the magnetic field above the PNS surface (r ∼ 25 − 30 km), which can be seen as the radially ordered field in the range of 25 km r 30 km (shown as a red circular band in the  (c) correspond to the cases with B 0 = 10 10 , 10 11 and 10 12 G, respectively. panel). This phenomenon has been already identified in the seminal work by Obergaulinger et al. (2014). Going more deeper inside (r 12 km), one can see a rather uniformally magnetized region (∼ 10 15 G), which corresponds to the unshocked core, where the density is close to the nuclear density (∼ 3 × 10 14 g/cm 3 ).
Finally we briefly state some speculations based on our results. The magnetars are typically observed as anomalous X-ray pulsars or soft gamma repeaters, and some of them are found in supernova remnants (e.g., Kaspi & Beloborodov 2017). The X-ray observations of supernova remnants having an association with magnetars (e.g. Kes 73, N49, and CTB 109) show that the explosion energies of these objects are comparable to those of canonical CCSNe (e.g., 10 51 erg, Vink & Kuiper 2006;Nakano et al. 2017). That might suggest that these magnetars may not requre rapid rotators with highly aspherical and energetic jets, but simply the normal neutrino-driven explosion as the central engine. Our strong initial field models (B12) would lead to the magnetar-class fields at the surface of the PNS (Fig. 8), but the explosion occurs by the neutrino heating. This would be consistent with the finding in the X-ray observations. To prove this bold speculation, we need to firstly follow a long-term evolution of our models because the diagnostic explosion energies (∼ 0.1 B) still fall short. Besides, the mass accretion rate of the fallback matter should be evaluated in the long-term simulation since the strong flow could bury the surface magnetic field into the crust (Shigeyama & Kashiyama 2018). These should be also verified in the 3D-MHD core-collapse simulations (e.g., Winteler et al. 2012;Mösta et al. 2014;Kuroda et al. 2020), which could be our next step to be taken.

SUMMARY AND DISCUSSION
We have investigated the impact of the magnetic field on the collapse of non-rotating stellar cores through 2D axisymmetric MHD simulations. Initially, 15.0, 18.4 or 27.0 M presupernova progenitors are threaded by only the poloidal component of the magnetic field. Since the azimuthal components of the velocity and magnetic field are zero initially with the 2D assumption, the evolution of the velocity and magnetic field is restricted in the poloidal components. We have performed numerical runs for the evolution of the stellar cores by varying the strength of the magnetic field inside the core between B = 10 10 and 10 12 G.
The stellar core collapses, and the neutrino-driven explosion occurs in all the computed models. For the 2D models explored in this work, an intriguing finding is that the main driver of the explosion is the neutrino heating regardless of the strength of the initial magnetic field. The magnetic field secondary contributes to the evolution of the stellar core. The strong magnetic field prevents the development of the neutrino-driven turbulence in the small scale compared to the weak magnetic field. This leads to the slow increase of the turbulent kinetic energy, leading to the slight delay of the neutrino-driven shock revival.
Finally we shall refer to the limitation of this study. The major limitation of this work is apparently the 2D assumption, which we shall take mainly for the sake of our code development which is doable as a first step without using huge computational resources. However, the explodability of CCSN models has been shown to be significantly affected by the dimension of the simulation (Nordhaus et al. 2010;Hanke et al. 2012;Dolence et al. 2013;Takiwaki et al. 2014;Nakamura et al. 2019;Nagakura et al. 2019b;Melson et al. 2020). As already mentioned above, we are planning to update our code to make 3D-MHD CCSN modeling possible as recently reported by  and Müller & Varma (2020). Our method of the ray-by-ray neutrino transport has room to be updated to the more sophisticated schemes such as the M1-closure scheme (Shibata et al. 2011;O'Connor 2015;Just et al. 2015Just et al. , 2018Kuroda et al. 2016;Skinner et al. 2016), the variable Eddington factor method Müller et al. 2010), and the discrete-angle (S n ) method (Liebendörfer et al. 2004;Sumiyoshi et al. 2005;Sumiyoshi & Yamada 2012;Nagakura et al. 2014Nagakura et al. , 2017Harada et al. 2019;Nagakura et al. 2019c;Iwakami et al. 2020). Another potential ingredient that can affect the postbounce dynamics is the treatment of gravity: general relativistic (GR) simulations (e.g., Kuroda et al. 2012Kuroda et al. , 2016Kuroda et al. , 2018Kuroda et al. , 2020Roberts et al. 2016;Mösta et al. 2014Mösta et al. , 2015 would be particularly important especially for the collapse of massive stars with high progenitor's compactness (e.g. O'Connor & Ott 2013;Sukhbold et al. 2016;Ertl et al. 2016).
Definitively much more elaborate study is needed to unravel the formation mechanism of magnetars, although we presented a very bold speculation in this study that our 2D non-rotating and strongly mangetized models could deserve further investigation toward a future comparison with the observations. The generation of the magnetar-class field in the context of the PNS dynamo process (Thompson & Duncan 1993) have gained considerable attention recently. In fact, Raynaud et al. (2020) presented the first 3D dynamo simulations, which showed the generation of the strong magnetic fields in the vicinity of the PNS (see, also Masada et al. 2020). The chiral magnetic effects have been reported to account for the origin of the strong magnetic fields in magnetars (Yamamoto 2016;Masada et al. 2018). All of these studies require a dedicated 3D MHD modeling, toward which we have made our first sail in this work. Table A1. Summary of the discretization method and the reconstruction of physical variables at the cell surface for conservation equations.  In the spherical coordinates, the conservation equations except the induction equations are simply described as follows; where Q is a conservative variable, F r , F θ , and F φ are the corresponding flues in the r-, θand φdirections and S is a source term. Following the methods proposed in Li & Li (2003) and Mignone (2014), we discretize the equation (B1) based on the finite volume method. Equation (B1) is integrated over the cell volume using the Gauss's theorem. The conservative variable defined at the cell center,Q i, j,k , is given by the volume average of Q over the cell volume, where i, j and k stand for the spacial index of the cell center in the r-θand φdirections, respectively. The flux terms in the left-hand-side of equation (B1) are discretized as follows; 1 r 2 ∂(r 2 F r ) ∂r ∼ r i+1/2 2 F r,i+1/2, j,k − r i−1/2 2 F r,i−1/2, j,k r i+1/2 3 /3 − r i−1/2 3 /3 , 1 r sin θ ∂(F θ sin θ) ∂θ ∼ r i+1/2 2 /2 − r i−1/2 2 /2 r i+1/2 3 /3 − r i−1/2 3 /3 F θ,i, j+1/2,k sin θ j+1/2 − F θ,i, j−1/2,k sin θ j−1/2 cos θ j−1/2 − cos θ j+1/2 , 1 r sin θ ∂F φ ∂φ ∼ r i+1/2 2 /2 − r i−1/2 2 /2 r i+1/2 3 /3 − r i−1/2 3 /3 θ j+1/2 − θ j−1/2 where F r,i±1/2, j,k , F θ,i, j±1/2,k and F φ,i, j,k±1/2 are the numerical flues at the cell surface in each direction (or the interpolated velocity for the source terms in equation (A10). In this case, the neutrino pressure is evaluated at the cell center). In addition, cot θ and 1/r in the source terms of the discrete conservation equations are replaced like below; Physical variables in the source terms are also evaluated by the cell-volume-averaged values. The time derivative of equation (B1) is discretized forward in time: where the superscript n stands for the number of time steps and ∆t is a step size. A Runge-Kutta time integration in the second-order accuracy is employed in the present code, which will be extended in higher-order accuracy (Asahina et al. in preparation).
In the induction equations (A6)-(A8), the terms including ψ are also discretized based on the finite volume method. However, considering the magnetic flux conservation, other terms should be better discretized based on the finite area method, namely, by integrating these terms over the cell area using the Stokes' theorem. The terms discretized by the finite area method are listed in Table A1. The first terms in the left-hand-side of the induction equations are time derivatives of each component of the magnetic field. Here, we introduce the area average of the variables over the cell areas defined at the cell center, respectively. The flux in the conservation form of the induction equation is related to the electric field defined as: Following the method of Skinner & Ostriker (2010) and Mignone (2014), the reconstruction of the volume average is obtained. First, we consider the reconstruction in the r-direction. Let a i be a cell-volume-averaged physical value inside zone i defined at the cell center, that is, The physical variables at the left and right interfaces inside the zone i are then given by where is a correction factor for curvature in the spherical coordinates and Here, ∆a i is the difference slope of the physical variable. In our code, the modified van Leer (VL) limiter proposed in Mignone (2014) is implemented for the slope limiter to achieve the total variation diminishing (TVD); The forward-and backward-difference slopes are given by where r i is the cell-volume-averaged value inside zone i; r i = r i+1/2 r i−1/2 r 3 dr r i+1/2 r i−1/2 r 2 dr = r i + 2r i ∆r 2 12r 2 i + ∆r 2 . (C10) The modified VL limiter is defined as follows (see Mignone 2014, for details); where We employ the modified VL limiter in all the directions for the reconstructions of both volume and area averages.
In the θ-direction, the cell-volume-averaged physical value inside zone j defined at the cell center is given by The physical variables at the left and right interfaces inside the zone j are then given by where is also the correction factor for curvature in the spherical coordinates and Here, ∆a j = ∆a F Figure D1. Snapshot of the 3D blast wave. Top left: The volume-rendered image of the magnetic pressure. Top right, bottom left and bottom right panels show the density structure in 2D slice of y = 0, x = 0 and z = 0, respectively. The uniform magnetic field is initially imposed in x-direction.
This is the same formula as the cell-volume-averaged physical value inside zone j (C14). Therefore, the reconstructions of the physical value at the cell surfaces are given by equations (C15) and (C16). On the other hand, when the area is perpendicular to the φ-direction, This simply results in the linear interpolation at the cell surface like equation (C25) or (C26) although the spacial index k in the equations is replaced by j.
In the φ-direction, the reconstruction of the physical variable is the same as the case with the volume-averaged reconstruction. It is simply the linear interpolation of the physical variable and given by equation (C25) or (C26).
For a concise summary, the physical variables reconstructed by the volume-and area-averaged methods are listed in Table A1.

APPENDIX D: BLAST WAVE IN A STRONGLY MAGNETIZED MEDIUM
For a code verification, we demonstrate a test of an MHD blast wave problem in 3D. This test is famous and performed in many previous works with similar setups (Londrillo & Del Zanna 2000;Stone et al. 2008;Stone & Gardiner 2009;Skinner & Ostriker 2010). The initial condition that we adopt is as follows. We set a static background of ρ = 1.0, p = 0.1. The background is magnetized as B x = 1.0 (in our definition, B is already normalized by √ 4π). In the background, we put a hot gas of p = 10 in the region of r < 0.1. The polytropic index s27.0B10 s27.0B12 s27.0B10 high s27.0B10 high Figure E1. Temporal evolution of the shock radii for our fiducial progenitor model (s27.0) with high numerical resolution (see the text). Solid, red and blue lines represent the weak (B 0 = 10 10 G) and strong (B 0 = 10 12 G) magnetic field model, respectively. For the reference, the results of fiducial runs are demonstrated by doted lines. of EOS, Γ, is set to 5 3 . The domain of [0, 0.5] × [0, π] × [0, 2π] of spherical polar coordinate is uniformly covered by 200 × 200 × 400 grids, respectively.
The time snapshot of t = 0.2 is shown in Fig. D1. The top left panel is the volume-rendered image of the magnetic pressure. The blast wave spherically propagates outward. The shock front is located with red sphere and green envelop in the figure. The central region has higher pressure and lower magnetic pressure. The region is significantly collimated due to the effect of the magnetic field. It is also important to check the density distribution. The top right, bottom left and bottom right panels are the 2D slice of y = 0, x = 0 and z = 0 planes, respectively. The structure of the density is quite similar to that in the other works (see Fig. 36 of Stone et al. 2008). The shock front corresponds to the high-density outer ring (shown as green) and that is almost spherical. The central low density region is the consequence of the rarefaction which propagates inward. The crescent density structure at the front of the shock in the x-direction is typically seen in this test.
The evolution of the shock does not depend on the coordinate system. Though the structure of the mesh is quite different in the x-z (r-θ) plane in the top right panel and the x-y (r-φ) plane in the bottom right panel, the density structures of them are quite similar. Such a coordinate independent feature would support the correctness of our code implementation.

APPENDIX E: RESOLUTION STUDY OF THE SHOCK REVIVAL WITH MAGNETIC FIELDS
The influence of the numerical resolution on the evolution of the MHD core collapse is investigated by changing the grid spacing in the θ-direction. The number of grid points in the θ-direction for fiducial runs is 128 as described in Section 2. We run s27.0B10 and s27.0B12 models with the resolution of ∆θ = π/256. The number of grid points in high resolution calculations is twice as larger as the fiducial runs. Fig. E1 shows the temporal evolution of the shock radius in the high resolution runs. Solid, red and blue lines correspond to the models s27.0B10 and s27.0B12, respectively. For the reference, the results of fiducial runs for both models are shown by thin dotted lines.
In both weak (s27.0B10) and strong (s27.0B12) magnetic field models, the shock radius in the high resolution run expands fast compared to that in the fiducial run. Solid lines reach 400 km faster than doted lines as shown in Fig E1. This result is consistent with a recent resolution study for the supernova simulation (Melson et al. 2020). The higher angular resolution provides more favorable explosion conditions. Also in the high resolution runs, the delay of the onset of the shock re-expansion with the stronger magnetic field models is observed at t pb ∼ 200 ms. The onset of the shock expansion in the weak magnetic field model is faster than that in the strong field model. As mentioned in Section 3, this tendency has been also observed in the fiducial resolution runs. Our resolution study demonstrates that the delay of the shock revival occurs regardless of the angular resolution in the θ-direction although more careful studies covering the wide range of resolution are necessary in order to draw a robust conclusion.