Two-dimensional numerical study for magnetic field dependence of neutrino-driven core-collapse supernova models
Jin Matsumoto, Tomoya Takiwaki, Kei Kotake, Yuta Asahina, Hiroyuki R. Takahashi
MMNRAS , 1–22 (2020) Preprint 7 October 2020 Compiled using MNRAS L A TEX style file v3.0
Two-dimensional numerical study for magnetic field dependence ofneutrino-driven core-collapse supernova models
J. Matsumoto (cid:63) , T. Takiwaki , K. Kotake , , Y. Asahina and H. R. Takahashi Research Institute of Stellar Explosive Phenomena, Fukuoka University, Fukuoka 814-0180, Japan Division of Science, National Astronomical Observatory of Japan, Tokyo 181-8588, Japan Department of Applied Physics, Fukuoka University, Fukuoka 814-0180, Japan Center for Computational Sciences, Tsukuba University, Ibaraki 305-8577, Japan Faculty of Arts and Sciences, Department of Natural Sciences, Komazawa University, Tokyo 154-8525, Japan
Accepted 2020 October 06. Received 2020 October 05; in original form 2020 August 20
ABSTRACT
We study the e ff ects of the magnetic field on the dynamics of non-rotating stellar cores byperforming two-dimensional (2D), magnetohydrodynamics (MHD) simulations. To this end,we have updated our neutrino-radiation-hydrodynamics supernova code to include MHD em-ploying a divergence cleaning method with both careful treatments of finite volume and areareconstructions. By changing the initial strength of the magnetic field, the evolution of 15 . . . M (cid:12) presupernova progenitors is investigated. An intriguing finding in our studyis that the neutrino-driven explosion occurs regardless of the strength of the initial magneticfield. For the 2D models presented in this work, the neutrino heating is the main driver for theexplosion, 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-driventurbulence in the small scale compared to the weak magnetic field. This results in the slowerincrease of the turbulent kinetic energy in the postshock region, leading to the slightly delayedonset of the shock revival for models with the stronger initial magnetic field. Key words: stars: massive – stars: magnetic field – supernovae: general
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 Liebendörfer et al. 2005; 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 explodemassive stars ( (cid:38) M (cid:12) ) is the neutrino mechanism (Bethe 1990). In this mechanism, neutrinos from the protoneutron star (PNS) heat thematerial behind the stalled shock, leading to the shock revival into explosion. The neutrino heating e ffi ciency is significantly enhancedby non-radial flows triggered by various hydrodynamic instabilities including neutrino-driven / PNS convection and the standing accretionshock instability (SASI; Blondin et al. 2003; Foglizzo et al. 2006). This has been confirmed in a growing number of self-consistent CCSNsimulations (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 erg ( ≡ ∼
10 B, which is (cid:63)
Email:[email protected], [email protected] © a r X i v : . [ a s t r o - ph . H E ] O c t J. Matsumoto et al. ten times larger than that of the canonical CCSNe . Concerning the SLSNe, the luminosity is 10–100 times higher than that of the typicalCCSNe. To explain the excess of the bright luminosity, several scenarios have been proposed including the interaction scenario betweenthe SN ejecta and its dense circumstellar medium (Chevalier & Irwin 2011; Moriya et al. 2013; Sorokina et al. 2016) or the pair-instabilitySN scenario (Barkat et al. 1967; Rakavy & Shaviv 1967; Terreran et al. 2017). Besides, the injection of the additional energy into the SNejecta by the central engine of a rapidly rotating proto-magnetar is also a promising scenario for the excess of luminosity (Kasen & Bildsten2010; Woosley 2010; Wang et al. 2015; Chen et al. 2016). For the two classes of the extreme events mentioned above, one could speculatethat a common solution requires the strong (magnetar-class), large-scale magnetic fields, which can directly couple the newly-born PNS (orproto-magnetar) to its surroundings via various MHD processes (Usov 1992; Thompson 1994; Bucciantini et al. 2009; Metzger et al. 2011).The MHD explosion mechanism originally proposed in the 1970s (Bisnovatyi-Kogan 1970; LeBlanc & Wilson 1970; Meier et al.1976; Müller & Hillebrandt 1979) has recently received considerable attention. To study the MHD mechanism in various contexts, extensivenumerical simulations have been conducted so far (e.g. Symbalisty 1984; Ardeljan et al. 2000; Kotake et al. 2004; Sawai et al. 2005; Shibataet al. 2006; Moiseenko et al. 2006; Suwa et al. 2007; Takiwaki et al. 2009; Takiwaki & Kotake 2011; Obergaulinger et al. 2006a,b; Winteleret al. 2012; Sawai & Yamada 2014, 2016; Obergaulinger & Aloy 2017, 2020; Obergaulinger et al. 2018; Mösta et al. 2014, 2015; Bugli et al.2020; Kuroda et al. 2020). However, the extremely huge explosion energy ( ∼
10 B) has yet to be obtained in these simulations. And therestill remains a big issue whether a combination of strong magnetic field and rapid rotation can be achieved in the precollapse iron core.Assuming the magnetic flux conservation, the strong surface magnetic field ( ∼ (cid:46)
100 G; Wade & MiMeS Collaboration2015) and stellar evolution calculations (Heger et al. 2005). In the latter case, enough amplification of the magnetic field during and after thecollapse of the massive star is necessary to facilitate the MHD mechanism. Since the magnetic field amplification due to the field wrappingas the consequence of rapid and strong di ff erential rotation scales linearly in time and takes a long time compared to the dynamical timescale of the system, the drastic and e ffi cient field amplification mechanism is required when the seed field of the massive star is adequatelyweak. The exponential growth of the magnetic field is beneficial to amplify the weak seed field of the stellar core to the dynamically relevantlevel. The magnetic field amplification due to the magnetorotational instability (MRI; Balbus & Hawley 1991) is a candidate for the e ffi cientfield amplification mechanism in the rotating stellar cores (Akiyama et al. 2003; Masada et al. 2006, 2007, 2012, 2015; Obergaulinger et al.2009; Sawai et al. 2013b; Sawai & Yamada 2014, 2016; Guilet et al. 2015; Guilet & Müller 2015; Mö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 tobe rotating slowly at the pre-collapse stage (Heger et al. 2005; Ott et al. 2006; Langer 2012). Especially, the stronger magnetic field wouldlead to the more e ffi cient 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 low-mass stars also suggested that an unmodeled, more e ffi cient angular momentum transport process is necessary to explain the spin period ofthe cores (Cantiello et al. 2014; Fuller et al. 2014). Observations of surface rotational velocities of B-type stars with strong magnetic fieldsalso favor slow rotators (Shultz et al. 2018). Even in such slowly rotating progenitors, it has been pointed out that the precollapse magneticfield, if su ffi ciently strong, could a ff ect the explosion dynamics. In the context of the slowly- and non-rotating progenitor, Endeve et al. (2010,2012) and Obergaulinger et al. (2014) studied the dynamics of the MHD core-collapse and the exponential amplification of the magnetic fielddue to the SASI and convection. More recently, Müller & Varma (2020) have addressed the role of the magnetic field in the neutrino-drivenexplosion by performing three-dimensional (3D) MHD simulations of a slowly rotating progenitor of a 15 M (cid:12) star.Joining in these e ff orts 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-rotatingand 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 transportis needed for a quantitative study of MHD CCSN modeling. However, only electron and anti-electron neutrinos were taken into account inObergaulinger et al. (2014) at that time (see, however, Obergaulinger & Aloy 2020). To revisit the problem, we have updated our supernovacode (3DnSNe) to include MHD by implementing a divergence cleaning method (Dedner et al. 2002) with both finite volume and areareconstructions based on Mignone (2014). Our base code deals with three-flavor neutrino transport (namely, ν e , ¯ ν e , ν X with ν X denotingthe heavy-lepton neutrinos) (Kotake et al. 2018) based on the Isotropic Di ff usion Source Approximation (IDSA) scheme (Liebendörfer et al.2009), in which a detailed code comparison was already made in spherically symmetric simulations (O’Connor et al. 2018), in 2D simulations(Kotake et al. 2018) using a widely used 20 M (cid:12) star of Woosley & Heger (2007) and in 3D simulations (Cabezón et al. 2018) using 15 M (cid:12) stars of Woosley & Weaver (1995) and 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 In order to explain observational features of HNe, an extensive study has been carried out so far in various contexts (e.g. Maeda & Nomoto 2003; Woosley& Bloom 2006; Tominaga 2009; Mazzali et al. 2014; Greiner et al. 2015; Metzger et al. 2015 for references therein). MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe core-collapse of non-rotating stellar cores in axisymmetry are presented in Section 3. In Section 4, we discuss the field configuration of theproto-magnetar based on our simulation results. Finally, we summarize and discuss our findings in Section 5. We have updated our supernova code, 3DnSNe (Takiwaki et al. 2016), that is designed for CCSN simulations in a 3D spherical coordinatesystem to the latest version. In this work, the code is now extended to an MHD code from a hydrodynamic (HD) one with spectral neutrinotransport that is solved by the IDSA scheme (Liebendörfer et al. 2009). We have updated the original (two-flavor, i.e. ν e , ¯ ν e ) IDSA schemein several manners, such that the evolution of the streaming neutrinos is self-consistently solved (Takiwaki et al. 2014) and that three-flavorneutrino transport is solved including approximate general relativistic corrections (e.g. Kotake et al. 2018 for more details). A detailed codecomparison has been performed in O’Connor et al. (2018) with one-dimensional (1D) geometry. The 3DnSNe code has been used in thefollowing 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 + ∇ · ( ρ v ) = , (1) ∂ ( ρ v ) ∂ t + ∇ · ( ρ vv + P t I − BB ) = − ρ ∇ Φ , (2) ∂ e ∂ t + ∇ · [( e + P t ) v − B ( v · B )] = − ρ v · ∇ Φ + Q E , (3) ∂ B ∂ t + ∇ · ( vB − Bv + ψ I ) = , (4) ∂ψ∂ t + c h ∇ · B = − c h c p ψ , (5) ∂ρ Y l ∂ t + ∇ · ( ρ Y l v ) = Γ l , (6) ∂ρ Z m ∂ t + ∇ · ( ρ Z m v ) + ρ Z m ∇ · v = Q m , (7) ∆Φ = π G ρ , (8)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 ofleptons: 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 unitmatrix. 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 solveequations (1)–(4) in a conservative form (see Appendix A for the treatment of equations 6 and 7). To get around the Carbuncle instability, weswitch from the HLLD scheme to the HLLE scheme (Einfeldt 1988) in the vicinity of strong shocks (see Kim et al. 2003 and the referencestherein). 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 divergenceerror 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 speedand the damping rate of the numerical divergence of the magnetic field, respectively (see Appendix B for details). Equation (5) is solvedusing HLLE scheme. To retain the total energy including gravitational binding energy, we use the method of Müller et al. (2010) in solvingequation (3). In equation (8), the spherically symmetric gravitational potential is taken in the form of phenomenological general relativisticpotential 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 divergencecleaning method is employed in our code, both the finite volume and area methods are used to solve the induction equation. The details ofthis special treatment, in addition to the reconstructions of the physical variables for the second-order accuracy in space, are described inAppendix 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 Previously, Takiwaki et al. (2012) and Takiwaki et al. (2014) employed the ZEUS-MP code (Hayes et al. 2006) where the (tensor-type) artificial viscositywas used to capture the shock (see also Iwakami et al. 2008). Suwa et al. (2010) and Suwa et al. (2016) utilized the ZEUS-2D code of Stone & Norman (1992).MNRAS000
100 G; Wade & MiMeS Collaboration2015) and stellar evolution calculations (Heger et al. 2005). In the latter case, enough amplification of the magnetic field during and after thecollapse of the massive star is necessary to facilitate the MHD mechanism. Since the magnetic field amplification due to the field wrappingas the consequence of rapid and strong di ff erential rotation scales linearly in time and takes a long time compared to the dynamical timescale of the system, the drastic and e ffi cient field amplification mechanism is required when the seed field of the massive star is adequatelyweak. The exponential growth of the magnetic field is beneficial to amplify the weak seed field of the stellar core to the dynamically relevantlevel. The magnetic field amplification due to the magnetorotational instability (MRI; Balbus & Hawley 1991) is a candidate for the e ffi cientfield amplification mechanism in the rotating stellar cores (Akiyama et al. 2003; Masada et al. 2006, 2007, 2012, 2015; Obergaulinger et al.2009; Sawai et al. 2013b; Sawai & Yamada 2014, 2016; Guilet et al. 2015; Guilet & Müller 2015; Mö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 tobe rotating slowly at the pre-collapse stage (Heger et al. 2005; Ott et al. 2006; Langer 2012). Especially, the stronger magnetic field wouldlead to the more e ffi cient 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 low-mass stars also suggested that an unmodeled, more e ffi cient angular momentum transport process is necessary to explain the spin period ofthe cores (Cantiello et al. 2014; Fuller et al. 2014). Observations of surface rotational velocities of B-type stars with strong magnetic fieldsalso favor slow rotators (Shultz et al. 2018). Even in such slowly rotating progenitors, it has been pointed out that the precollapse magneticfield, if su ffi ciently strong, could a ff ect the explosion dynamics. In the context of the slowly- and non-rotating progenitor, Endeve et al. (2010,2012) and Obergaulinger et al. (2014) studied the dynamics of the MHD core-collapse and the exponential amplification of the magnetic fielddue to the SASI and convection. More recently, Müller & Varma (2020) have addressed the role of the magnetic field in the neutrino-drivenexplosion by performing three-dimensional (3D) MHD simulations of a slowly rotating progenitor of a 15 M (cid:12) star.Joining in these e ff orts 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-rotatingand 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 transportis needed for a quantitative study of MHD CCSN modeling. However, only electron and anti-electron neutrinos were taken into account inObergaulinger et al. (2014) at that time (see, however, Obergaulinger & Aloy 2020). To revisit the problem, we have updated our supernovacode (3DnSNe) to include MHD by implementing a divergence cleaning method (Dedner et al. 2002) with both finite volume and areareconstructions based on Mignone (2014). Our base code deals with three-flavor neutrino transport (namely, ν e , ¯ ν e , ν X with ν X denotingthe heavy-lepton neutrinos) (Kotake et al. 2018) based on the Isotropic Di ff usion Source Approximation (IDSA) scheme (Liebendörfer et al.2009), in which a detailed code comparison was already made in spherically symmetric simulations (O’Connor et al. 2018), in 2D simulations(Kotake et al. 2018) using a widely used 20 M (cid:12) star of Woosley & Heger (2007) and in 3D simulations (Cabezón et al. 2018) using 15 M (cid:12) stars of Woosley & Weaver (1995) and 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 In order to explain observational features of HNe, an extensive study has been carried out so far in various contexts (e.g. Maeda & Nomoto 2003; Woosley& Bloom 2006; Tominaga 2009; Mazzali et al. 2014; Greiner et al. 2015; Metzger et al. 2015 for references therein). MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe core-collapse of non-rotating stellar cores in axisymmetry are presented in Section 3. In Section 4, we discuss the field configuration of theproto-magnetar based on our simulation results. Finally, we summarize and discuss our findings in Section 5. We have updated our supernova code, 3DnSNe (Takiwaki et al. 2016), that is designed for CCSN simulations in a 3D spherical coordinatesystem to the latest version. In this work, the code is now extended to an MHD code from a hydrodynamic (HD) one with spectral neutrinotransport that is solved by the IDSA scheme (Liebendörfer et al. 2009). We have updated the original (two-flavor, i.e. ν e , ¯ ν e ) IDSA schemein several manners, such that the evolution of the streaming neutrinos is self-consistently solved (Takiwaki et al. 2014) and that three-flavorneutrino transport is solved including approximate general relativistic corrections (e.g. Kotake et al. 2018 for more details). A detailed codecomparison has been performed in O’Connor et al. (2018) with one-dimensional (1D) geometry. The 3DnSNe code has been used in thefollowing 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 + ∇ · ( ρ v ) = , (1) ∂ ( ρ v ) ∂ t + ∇ · ( ρ vv + P t I − BB ) = − ρ ∇ Φ , (2) ∂ e ∂ t + ∇ · [( e + P t ) v − B ( v · B )] = − ρ v · ∇ Φ + Q E , (3) ∂ B ∂ t + ∇ · ( vB − Bv + ψ I ) = , (4) ∂ψ∂ t + c h ∇ · B = − c h c p ψ , (5) ∂ρ Y l ∂ t + ∇ · ( ρ Y l v ) = Γ l , (6) ∂ρ Z m ∂ t + ∇ · ( ρ Z m v ) + ρ Z m ∇ · v = Q m , (7) ∆Φ = π G ρ , (8)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 ofleptons: 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 unitmatrix. 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 solveequations (1)–(4) in a conservative form (see Appendix A for the treatment of equations 6 and 7). To get around the Carbuncle instability, weswitch from the HLLD scheme to the HLLE scheme (Einfeldt 1988) in the vicinity of strong shocks (see Kim et al. 2003 and the referencestherein). 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 divergenceerror 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 speedand the damping rate of the numerical divergence of the magnetic field, respectively (see Appendix B for details). Equation (5) is solvedusing HLLE scheme. To retain the total energy including gravitational binding energy, we use the method of Müller et al. (2010) in solvingequation (3). In equation (8), the spherically symmetric gravitational potential is taken in the form of phenomenological general relativisticpotential 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 divergencecleaning method is employed in our code, both the finite volume and area methods are used to solve the induction equation. The details ofthis special treatment, in addition to the reconstructions of the physical variables for the second-order accuracy in space, are described inAppendix 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 Previously, Takiwaki et al. (2012) and Takiwaki et al. (2014) employed the ZEUS-MP code (Hayes et al. 2006) where the (tensor-type) artificial viscositywas used to capture the shock (see also Iwakami et al. 2008). Suwa et al. (2010) and Suwa et al. (2016) utilized the ZEUS-2D code of Stone & Norman (1992).MNRAS000 , 1–22 (2020)
J. Matsumoto et al. opacity set of Bruenn (1985). In this run, 20 energy groups that logarithmically spread from 1 to 300 MeV are employed. We use the equationof state (EOS) by Lattimer & Swesty (1991) (incompressibility K =
220 MeV).We employ the non-rotating presupernova progenitors of 15 .
0, 18 . . M (cid:12) of Woosley et al. (2002). As for the initial configurationof the magnetic fields, we assume a simple topology following Suwa et al. (2007); Takiwaki et al. (2014); Obergaulinger et al. (2014). Themagnetic field is given by a vector potential in the φ -direction of the form A φ = B r r + r r sin θ , (9)where r = r , is smaller than r , while itis like dipole field when r is larger than r . B determines the strength of the magnetic field inside the core ( r < r ). In this study, we set B = , 10 or 10 G. The model name is labelled as ‘s27.0B10’, which represents the 27 . M (cid:12) model with B = G. We chooses27.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 ∼ −
500 ms after bounce, depending on the progenitor models. In most ofthe models, we terminate the simulations at the final time seeing that the diagnostic explosion energies are greater than 10 erg. We leavethe 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 inthe 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 polarangle in the θ -direction is uniformly divided into ∆ θ = π/ r = r = z -axis in our 2D run). A numerical resolution testis given in Appendix E. We first describe overall evolution of the magnetized and non-rotating stellar core for our fiducial model (s27.0B10) in Section 3.1. Thenin the subsequent sections, we move on to present results focusing on the impact of the initial magnetic field strength on the postbounceevolution. The progenitor dependence of the shock evolution is presented in Section 3.4. M (cid:12) 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 < x ( x > t pb = t pb denotes thepostbounce 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 ∼ 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. Beforethe shock stall ( t pb (cid:46)
100 ms), the flow is almost restricted in radial direction. The split-monopole like configuration is made because themagnetic field is "frozen-in" with respect to the matter motion. The electric resistivity of the magnetic field is so small that it is disregarded inthis work, which can be well justified in the CCSN environment (Sawai et al. 2013a). The initial vector potential (equation 9) gives magneticloops on the equatorial region at around r ∼ x (cid:38)
30 km and z =
0) in Fig. 1a. The center of loops is located at around x ∼
45 km and seen as a smallblueish region.As the (maximum) shock radius starts to gradually shrink after t pb (cid:38)
100 ms (e.g. Fig. 2a), it gradually deviates from the shock trajectoryof the corresponding 1D model (black solid line in Fig. 2a). This marks the growth of non-spherical motions in the postshock region. Onecan 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 penetrationof the magnetic field lines (thin red curves in the right-half panel) into the postshock region (high entropy region in the left-half panel), whichmakes the field configuration much more complicated than that outside the shock. In our ideal MHD simulations, the field amplification inthe 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 di ff erentiate the origin of the "non-radial" motions either originating predominantly fromthe SASI or neutrino-driven convection because the SASI is liable to be overestimated in 2D compared to 3D simulations (e.g. Hanke et al.2012, 2013; Ferná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 theformation of the high entropy region (colored by red in the entropy plot) is a common feature of 2D neutrino-driven explosion models. The
MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe Figure 1.
Time evolution of entropy per baryon ( x <
0) and magnetic field ( x >
0) for fiducial model (s27.0B10). Panels (a), (b), (c) and (d) correspond to t pb = t pb denotes the postbounce time. highly aspherical shock (Fig. 1d) continues to propagate strongly toward the south pole up to a radius of ∼ t fin ∼
500 ms) for this model. The magnetic field configuration is similar to the blast morphology as in the previous snapshots.
Fig. 2 summarizes the time evolution of the maximum shock radius (top panels), the ratio of the advection timescale to the neutrino-heatingtimescale, τ adv /τ heat , (bottom left panel) and gain mass (bottom right panel) for the 27 . M (cid:12) model. Red, green and blue lines in each panelare the cases for B = , 10 and 10 G, respectively.In Fig. 2a, a black line, as already mentioned, denotes the shock evolution of the 1D HD model of the 27 . M (cid:12) star, which is shown asa 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 thepanel very carefully, one might notice an interesting tendency about the slight delay of the onset of the shock revival between the three MHDmodels.
MNRAS000
MNRAS000 , 1–22 (2020)
J. Matsumoto et al. r s h [ k m ] Time after bounce [ms](a) s27.0 B = 10 G B = 10 G B = 10 G1D HD 6080100120140160180200220 160 180 200 220 240 r s h [ k m ] Time after bounce [ms](b)0.00.20.40.60.81.01.21.4 50 100 150 200 250 τ a d v / τ h e a t Time after bounce [ms](c) s27.0B10s27.0B11s27.0B12 10 − − −
50 100 150 200 250 G a i n m a ss [ M ⊙ ] (d) Time after bounce [ms]s27.0B10s27.0B11s27.0B12 Figure 2.
Time evolution of shock radius (top panels), τ adv /τ heat (bottom left panel) and gain mass (bottom right panel) for s27 . M (cid:12) progenitor. Red, greenand blue lines in each panel are the cases for B = , 10 and 10 G, respectively. Black lines in top panels show the evolution of the 1D HD simulationfor the 27 . M (cid:12) model. 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 = t pb =
250 ms). The shock revival time is indeed delayed for the strong initial field model (blue curve) comparing to the weak initialfield 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 τ adv = M g ˙ M , (10)where M g is the mass enclosed in the gain layer (gain mass) and ˙ M is the mass-accretion rate through the shock. The neutrino-heatingtimescale is defined by τ heat = | E tot , g | ˙ Q heat , (11)where | E tot , g | is the total energy of the material in the gain layer and ˙ Q heat is the neutrino-heating rate in this region. Since the residency timeof matter in the gain region is related to the exposure of the material to neutrino heating, τ adv /τ heat (cid:38) τ adv /τ heat rises toward unity rapidly at around t pb ∼
200 ms. It is noted thatthis shock revival timescale is consistent with that by Hanke et al. (2013) who conducted the 2D model using the same progenitor (27 M (cid:12) starof 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 weaklymagnetized model (red line), which is followed in order by the moderately magnetized model (green line) and the strongly magnetizedmodel (blue line). This feature is closely linked to the shock evolution after t pb ∼
200 ms, namely, the onset of the shock revival andthe subsequent runaway shock expansion is delayed for the strongly magnetized models. This indicates that the neutrino heating mainly
MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe L u m i n o s i t y [ e r g/ s ] Time after bounce [ms](a) s27.0 red B = 10 Ggreen B = 10 Gblue B = 10 G ν e ¯ ν e ν X M e a n e n e r g y [ M e V ] Time after bounce [ms](b)
Figure 3.
Panel (a): temporal evolution of neutrino luminosity for fiducial progenitor model (s27.0). Panel (b): time evolution of mean energy of neutrinos forfiducial progenitor model (s27.0). Line types and colors have the same meanings as those in panel (a).
Figure 4.
2D distribution of plasma β for fiducial progenitor model (s27.0) at t pb =
200 ms. Panels (a), (b) and (c) correspond to the cases with B = ,10 and 10 G, respectively. contributes to the runaway shock expansion as shown in Fig. 2a. And the magnetic field secondary a ff ects the shock evolution. The sametendency 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 = , 10 and 10 G, respectively. The evolution of ν e and ¯ ν e with di ff erent magnetic fields in the luminosity and mean energy are almost identical up to 200ms and represented as blue solid and dashed lines, respectively. On the other hand, the evolution of ν X with di ff erent magnetic fields almostoverlaps during the whole calculation time. Given the results mentioned above, it may not be so surprising that the initial magnetic fieldshave 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)who have done a systematic 2D simulations covering a wide range of the (non-magnetized) progenitor models. Regardless of the di ff erencein the neutrino transport scheme, the peak of the ν e / ¯ ν e luminosity (for the same progenitor) is ∼ × erg / s at around t pb ∼
100 ms, whichis consistent with our results. After the onset of shock revival ( t pb (cid:38)
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 thanSumma 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 thedi ff erence in the evolution of τ adv /τ heat (Fig. 2c) for models with stronger initial magnetic fields. This suggests that the stronger initial fielda ff ects the advection timescale predominantly than the neutrino heating timescale. In what follows, we explore how the strong initial fieldcould a ff ect 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 MNRAS000
8% higher thanSumma 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 thedi ff erence in the evolution of τ adv /τ heat (Fig. 2c) for models with stronger initial magnetic fields. This suggests that the stronger initial fielda ff ects the advection timescale predominantly than the neutrino heating timescale. In what follows, we explore how the strong initial fieldcould a ff ect 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 MNRAS000 , 1–22 (2020)
J. Matsumoto et al. E k i n e ( v θ ) [ e r g ] Time after bounce [ms](a) s27.0B10s27.0B11s27.0B12 0.00.51.01.52.02.53.03.54.04.55.0 160 180 200 220 240 E k i n e ( v θ ) / [ e r g ] Time after bounce [ms](b) s27.0B10s27.0B11s27.0B12051015202530 50 100 150 200 250 A d v ec t i o n t i m e s c a l e [ m s ] Time after bounce [ms](c) s27.0B10s27.0B11s27.0B12 10 t = 195 ms e k i n e ( l ) [ e r g/ c m ] l s27.0B10s27.0B11s27.0B12 Figure 5.
Temporal evolution of kinetic energy (top panels), time evolution of advection timescale (bottom left panel) and spectrum of kinetic energy at t = θ -component of velocityis 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 = , 10 and 10 G, respectively. 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 isa 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 . M (cid:12) models with di ff erent initial magnetic fields. It is shown that the plasma β behind the shock is much larger than unity ( (cid:38) ∼ ), meaningthat the shock revival in our models is predominantly driven by neutrino heating. This is absolutely not the case in the MHD explosion inthe context of rapidly rotating and strongly magnetized cores (e.g. Kuroda et al. 2020), where the plasma β of ∼ unity is often achieved. Thisis because the rapidly rotating PNS with strong di ff erential rotation at the surface results in the magnetic field winding and the subsequentincrease of the magnetic pressure. The increased magnetic pressure derives the shock expansion toward the polar directions. Furthermore it isimportant to note that the neutrino-driven shock revival is obtained at t pb ∼
200 ms for our weakly magnetized model (s27.0B10). This is instark contrast with Obergaulinger et al. (2014) who obtained a very late onset of explosion at t pb ∼
700 ms for a 15 M (cid:12) progenitor employedin the work, though assuming the same magnetic field strength (10 G). More detailed comparison between our models and the previousstudies is presented in Section 3.4.
In the previous subsections, we have shown that the neutrino heating plays a dominant role in triggering the explosion of our 27 M (cid:12) models,whereas the magnetic field plays a secondary role, namely, to delay the explosion onset. We now proceed to clarify the reason by focusingon 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 timescalein the gain region for the model series of s27.0, respectively. In Fig. 5a, the contribution from the lateral ( θ )-component of the velocity ( v θ ) istaken 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 MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe
100 150 200 250 300 350 400(a) E m ag ( B θ ) [ e r g ] Time after bounce [ms]s27.0B10s27.0B11s27.0B12 10 t = 195 ms e m ag ( l ) [ e r g/ c m ] l s27.0B10 ( × )s27.0B11 ( × )s27.0B12 Figure 6.
Panel (a): temporal evolution of magnetic energy (only θ -componet) in gain region for fiducial progenitor model (s27.0). Red, green and blue lines arethe cases for B = , 10 and 10 G, respectively. Panel (b): spectrum of magnetic energy in gain region defined by equation (13) for fiducial progenitormodel at t =
197 ms. The meanings of colors for lines are the same as those in panel (a). In order to normalize the spectra with di ff erent field strength, thosefor the case with B = and 10 G are multiplied by a factor of 10 and 10 , respectively. 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 ∼ erg toward the final simulation time regardless of the di ff erent initial field strength.Looking more closely at the linear phase ( t pb (cid:46)
200 ms), Fig. 5b depicts that the growth of the kinetic energy is fastest (biggest) formodel s27.0B10 (red line) compared to the more strongly magnetized models (green and blue lines). This feature, as previously identifiedin the 3D MHD simulations of Endeve et al. (2012) (but with more idealized setting), is also consistent with the earlier onset of the shockrevival as seen in Fig. 2b.Since the neutrino heating timescale is similar among the s27.0 models (e.g. Section 3.2), the di ff erence of the advection timescale in thegain 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 weaklymagnetized model (red line) is actually longer than that of the strongly magnetized model (blue line) around the explosion onset ( t pb ∼ e kine ( l ). Following Hanke et al. (2012), it is defined as, e kine ( l ) = (cid:88) m = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω Y ∗ lm ( θ, φ ) √ ρ v θ ( r , θ, φ ) d Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (12)where Y lm is the spherical harmonics of degree l and m , and Ω is a solid angle. Note in our 2D simulations, m = r =
100 km) in the postshock region and at around theexplosion 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 (cid:38)
10) is bigger for the weaker magnetic field case (red line) thanthat for the stronger field case (blue line), although the energy density of the smaller-order modes is comparable in all the three cases. Thisimplies that the strong magnetic field, most likely due to the magnetic tension, prevents the growth of the turbulent motions down to smallscales (at larger l ). The suppression of the turbulent energy at larger l is reconciled with the slow increase of the non-radial kinetic energy forthe 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 correspond-ing time-averaged energy spectra of the magnetic turbulence at t pb =
195 ms (right panel). Fig. 6a shows that the lateral magnetic energy isexponentially amplified up to the explosion onset ( t pb ∼
200 ms) in all the three models. The exponential growth terminates when the shockrevival initiates ( t pb (cid:38)
200 ms). In the non-linear phase ( t pb (cid:38)
200 ms), the lateral magnetic energies are shown to be almost kept constantwith time, whose strength is bigger for the strongly magnetized model. It is noted that the saturated values of model s27.0B12 ( ∼ erg,blue line), s27.0B11 ( ∼ erg, green line), and s27.0B10 ( ∼ erg, green line) di ff er by the two orders-of-magnitudes, respectively,which is proportional to the square of the initial magnetic energy (i.e. ∝ B ).Since the magnetic pressure / force does not play a crucial role in the shock revival (see Fig. 4 and the explanation in the last paragraphof Section 3.2), the magnetic field can be amplified by the compression and stretching due to the non-radial fluid motions. Actually, bycomparing 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 MNRAS000
200 ms), the lateral magnetic energies are shown to be almost kept constantwith time, whose strength is bigger for the strongly magnetized model. It is noted that the saturated values of model s27.0B12 ( ∼ erg,blue line), s27.0B11 ( ∼ erg, green line), and s27.0B10 ( ∼ erg, green line) di ff er by the two orders-of-magnitudes, respectively,which is proportional to the square of the initial magnetic energy (i.e. ∝ B ).Since the magnetic pressure / force does not play a crucial role in the shock revival (see Fig. 4 and the explanation in the last paragraphof Section 3.2), the magnetic field can be amplified by the compression and stretching due to the non-radial fluid motions. Actually, bycomparing 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 MNRAS000 , 1–22 (2020) J. Matsumoto et al. 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 − , 10 − and 10 − , respectively, whereas the kinetic energy at the time is almost 10 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 ifthe explosion were much more delayed like in Obergaulinger et al. (2014), the magnetic field amplification could continue until it becomessu ffi ciently 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 e mag ( l ) = (cid:88) m = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω Y ∗ lm ( θ, φ ) B θ ( r , θ, φ ) d Ω (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (13)In order to normalize the spectra with the di ff erent initial field strength, the spectra for the case with B = and 10 G are multiplied bya factor of 10 and 10 , respectively. This is reasonable because the magnetic energy is proportional to B and that of B = and 10 Gare 10 and 10 times smaller than that of B = G, respectively.From Fig. 6b, one can see that the energy spectrum of the magnetic turbulence typically decreases with l , i.e. the magnetic energyis mainly stored in a large-scale field. The slope is, however, shallower for the weak field model (s27.0B10, red line) comparing with thestrong field models (green and blue lines). This indicates that a small scale (turbulent) magnetic field can be more preferentially developedfor 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 turbulencein the small scale compared to that for the weak initial field model. ff erent progenitor models In order to investigate the impact of the initial magnetic field on the di ff erent progenitor models, we present results for the 15 . . M (cid:12) models. Similar to Figs. 2a and 2b, but Fig. 7 shows the evolution of the shock for model 15 . M (cid:12) (top panels) and model 18 . M (cid:12) (bottom panels), respectively. In both of the models, the shock revival occurs at t pb ∼
200 ms (Fig. 7) regardless of the di ff erence in the initialmagnetic 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 27 M (cid:12) models as already mentioned in the lastsection.The 2D HD simulation of the 18 . M (cid:12) 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 di ff erence of neutrino transportscheme between the two codes, we cannot unambiguously specify the reason of the discrepancy. But already in the 1D comparison workof 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 especiallylater than t pb ∼
200 ms relative to the other codes, leading to the enhanced heating rate in the gain region t pb (cid:38)
250 ms (see, green line oftheir Fig. 5). This could be one of the reasons of the early shock revival seen in our 18 . M (cid:12) run. Be that as it may, the shock revival timeand the neutrino properties (Figs. 2 and 3) show a good agreement with those of the 27 M (cid:12) model of Summa et al. (2016). The match maybe simply by chance, and a detailed comparison of 2D CCSN models from di ff erent 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 andthe transport schemes, Just et al. (2018) showed in their detailed, 2D systematic simulations that a simplified treatment in the neutrino-pairprocesses makes the onset time of the shock revival significantly earlier ( ∼ −
400 ms) compared to the case without such simplification(compare the explosion time of models s20-rbr and s20-rbr-pp{1 / /
3} in their Table 1). Since we employ the simplified prescription (e.g. theassumption of the isotropy and local-thermal-equilibrium with respect to the heavy-lepton neutrinos), this might be also the reason of ourearlier onset of the explosion.As already mentioned before, the 2D MHD simulations of the non-rotating cores of the 15 . M (cid:12) star (Woosley et al. 2002) werereported in Obergaulinger et al. (2014) with varying the initial field strength (like in this study). In their strong initial field model (10 G),they obtained a magnetically assisted explosion at t pb ∼
600 ms, which is about ∼
200 ms earlier than their weak field model (10 G). If theneutrino-driven shock revival is so much delayed, they pointed out that an equipartition between the turbulent kinetic and magnetic energywas archived due to the growing turbulence over the long period of time. This may seem to contradict with our results, i.e. the slight delayof the explosion onset for the strongly magnetized model (s15.0B12) as shown in Fig. 7a and Fig. 7b. However, in our models, the shockrevival occurs much earlier until the equipartition could be achieved. The di ff erence of the 2D hydrodynamics, which could stem from theomission of heavy-lepton neutrinos in Obergaulinger et al. (2014) and also from the di ff erence 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 thatthey illuminate two faces of the magnetic fields in the explosion dynamics of non-rotating CCSNe, namely, either in the case of relativelyearlier onset of explosion as pointed out in our study or in the later onset of explosion as studied by Obergaulinger et al. (2014). In either of the cases, 3D simulation is apparently needed to draw a robust conclusion to clarify the roles of magnetic fields on the explosion onset in thenon- or slowing-rotating cores (e.g. Müller & Varma 2020). MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe r s h [ k m ] Time after bounce [ms](a) s15.0 B = 10 G B = 10 G B = 10 G1D HD 6080100120140160180200220240260 160 180 200 220 240 r s h [ k m ] Time after bounce [ms](b) s15.0050100150200250300350400 0 100 200 300 400 500 r s h [ k m ] Time after bounce [ms](c) s18.4 B = 10 G B = 10 G B = 10 G1D HD 6080100120140160180200220 160 180 200 220 240 r s h [ k m ] Time after bounce [ms](d) s18.4
Figure 7.
Evolution of shock radii for models s15.0 and s18.4. Red, green and blue lines in each panel are the cases for B = , 10 and 10 G,respectively. Black lines show the evolution of the 1D calculation (HD) for each model.
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 magneticfield of the PNS is also an interesting topic. Especially, the origin of the strong magnetic field of the magnetar is still under debate overthe three decades (e.g., Duncan & Thompson 1992; Paczynski 1992; Thompson & Duncan 1993; Kouveliotou et al. 1998, Lai 2001; Kaspi& Beloborodov 2017 for a review). Two possible formation scenarios have been proposed to account for the strong magnetic field of themagnetars. 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 mainlydue 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 estimatethe strength of the magnetic field inside the PNS. Since the initial magnetic field is given by equation (9) and uniform ( B = B ) inside thecore ( r < km), the magnetic field strength of the PNS is estimated as follows: B PNS ∼ G (cid:32) B G (cid:33)(cid:32)
30 km r PNS (cid:33) , (14)where r PNS denotes the radius of the PNS. From this rough estimate, one can anticipate that a magnetar-class magnetic field could be formedin the vicinity of the PNS. Note that the radius of the PNS is defined by the iso-density surface of 10 g / cm .Fig. 8 compares the 2D distribution of the magnetic field strength for our fiducial runs (s27.0) with B = G (panel a), B = (panel b) and B = G (panel c), respectively. One can see that the magnetic field configuration in the central region ( r (cid:46)
50 km) lookssimilar between the three panels. All the three panels reveal a common shortcoming of our 2D models, namely, the magnetic field is strongestalong the symmetry axis (the z -axis) because of the reflective boundary condition there. In addition, an artificial structure of the magnetic field MNRAS000
50 km) lookssimilar between the three panels. All the three panels reveal a common shortcoming of our 2D models, namely, the magnetic field is strongestalong the symmetry axis (the z -axis) because of the reflective boundary condition there. In addition, an artificial structure of the magnetic field MNRAS000 , 1–22 (2020) J. Matsumoto et al.
Figure 8.
2D distribution of magnetic field strength for fiducial progenitor model (s27.0) near at the final simulation time ( t pb ∼
400 ms). Panel (a), (b) and (c)correspond to the cases with B = , 10 and 10 G, respectively. 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 (cid:46)
30 km is typically ∼ G (shown in red), which isconsistent with the rough estimate of equation (14) in the PNS interior. The field strength becomes smaller below the PNS surface (30 ∼ (cid:38) r (cid:38)
12 km), which is shown in green or yellowish region in the panel. This region is convectively unstable due to the negativelepton gradient, whereas the region above is convective stable due to the positive entropy gradient. In the convective region, it is known in2D simulations that the magnetic field is expelled due to the so-called convective flux expulsion (e.g. Galloway & Weiss 1981; Davidson2001). This not only explains the reduction of the field strength there, but also the accumulation of the magnetic field above the PNS surface( r ∼ −
30 km), which can be seen as the radially ordered field in the range of 25 km (cid:46) r (cid:46)
30 km (shown as a red circular band in thepanel). This phenomenon has been already identified in the seminal work by Obergaulinger et al. (2014). Going more deeper inside ( r (cid:46) ∼ G), which corresponds to the unshocked core, where the density is close tothe nuclear density ( ∼ × g / cm ).Finally we briefly state some speculations based on our results. The magnetars are typically observed as anomalous X-ray pulsarsor soft gamma repeaters, and some of them are found in supernova remnants (e.g., Kaspi & Beloborodov 2017). The X-ray observationsof supernova remnants having an association with magnetars (e.g. Kes 73, N49, and CTB 109) show that the explosion energies of theseobjects are comparable to those of canonical CCSNe (e.g., 10 erg, Vink & Kuiper 2006; Nakano et al. 2017). That might suggest thatthese magnetars may not requre rapid rotators with highly aspherical and energetic jets, but simply the normal neutrino-driven explosion asthe central engine. Our strong initial field models (B12) would lead to the magnetar-class fields at the surface of the PNS (Fig. 8), but theexplosion 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 ( ∼ . 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 . . M (cid:12) presupernova progenitors are threaded by only the poloidal component of the magnetic field. Since theazimuthal components of the velocity and magnetic field are zero initially with the 2D assumption, the evolution of the velocity and magneticfield is restricted in the poloidal components. We have performed numerical runs for the evolution of the stellar cores by varying the strengthof the magnetic field inside the core between B = and 10 G.The stellar core collapses, and the neutrino-driven explosion occurs in all the computed models. For the 2D models explored in thiswork, an intriguing finding is that the main driver of the explosion is the neutrino heating regardless of the strength of the initial magneticfield. The magnetic field secondary contributes to the evolution of the stellar core. The strong magnetic field prevents the development of theneutrino-driven turbulence in the small scale compared to the weak magnetic field. This leads to the slow increase of the turbulent kineticenergy, 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
MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe 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 a ff ected 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). Especially,3D purely HD simulations of non-rotating progenitors indicate negative impacts on the explosion compared to 2D HD simulations. Themagnetic field may assist the explosion in the 3D simulation of CCSNe (e.g. Müller & Varma (2020)). As already mentioned above, we areplanning to update our code to make 3D-MHD CCSN modeling possible as recently reported by Obergaulinger & Aloy (2020) 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. 2015, 2018; Kuroda et al. 2016; Skinner et al. 2016), the variable Eddingtonfactor method (Buras et al. 2006; 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. 2014, 2017; Harada et al. 2019; Nagakura et al. 2019c; Iwakami et al. 2020). Another potentialingredient that can a ff ect the postbounce dynamics is the treatment of gravity: general relativistic (GR) simulations (e.g., Kuroda et al. 2012,2016, 2018, 2020; Roberts et al. 2016; Mösta et al. 2014, 2015) would be particularly important especially for the collapse of massive starswith high progenitor’s compactness (e.g. O’Connor & Ott 2013; Sukhbold et al. 2016; Ertl et al. 2016; Ebinger et al. 2020).We have investigated the dependence of the magnetic field strength on the dynamics of the core-collapse of the massive star in thecontext of the slowly-rotating core. In order to independently clarify the role of the magnetic field and rotation of the star, we assume no-rotating cores in this study. The rotational dependence of the core on the dynamics of the MHD core-collapse is within the scope of our workand will be reported in our subsequent paper.Definitively much more elaborate study is needed to unravel the formation mechanism of magnetars, although we presented a verybold speculation in this study that our 2D non-rotating and strongly mangetized models could deserve further investigation toward a futurecomparison 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, whichshowed the generation of the strong magnetic fields in the vicinity of the PNS (see, also Masada et al. 2020). The chiral magnetic e ff ects havebeen reported to account for the origin of the strong magnetic fields in magnetars (Yamamoto 2016; Masada et al. 2018). All of these studiesrequire a dedicated 3D MHD modeling, toward which we have made our first sail in this work. ACKNOWLEDGEMENTS
We thank K. Nakamura, M. Bugli, Y. Masada, Y. Matsumoto, K. Tomisaka, and H.-Th. Janka for useful and stimulating discussions. Numeri-cal computations were carried out on Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japanand on Cray XC40 at YITP in Kyoto University. This work was supported by Research Institute of Stellar Explosive Phenomena at FukuokaUniversity and the associated project (No. 207002), and also by JSPS KAKENHI Grant Number (JP17K14260, JP17H05206, JP17K14306,JP17H01130, JP17H06364, JP18H01212, JP18K13591, JP19K23443, JP20K14473, JP20K11851, JP20H01941 and JP20H00156). This re-search was also supported by MEXT as “Program for Promoting researches on the Supercomputer Fugaku” (Toward a unified view of heuniverse: from large scale structures to planets) and JICFuS.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.
REFERENCES
Akiyama S., Wheeler J. C., Meier D. L., Lichtenstadt I., 2003, ApJ, 584, 954Ardeljan N. V., Bisnovatyi-Kogan G. S., Moiseenko S. G., 2000, A&A, 355, 1181Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Barkat Z., Rakavy G., Sack N., 1967, Phys. Rev. Lett., 18, 379Bethe H. A., 1990, Reviews of Modern Physics, 62, 801Bisnovatyi-Kogan G. S., 1970, Azh, 47, 813Blondin J. M., Mezzacappa A., DeMarino C., 2003, ApJ, 584, 971Bruenn S. W., 1985, ApJS, 58, 771Bucciantini N., Quataert E., Metzger B. D., Thompson T. A., Arons J., Del Zanna L., 2009, MNRAS, 396, 2038Bugli M., Guilet J., Obergaulinger M., Cerdá-Durán P., Aloy M. A., 2020, MNRAS, 492, 58Buras R., Rampp M., Janka H. T., Kifonidis K., 2006, A&A, 447, 1049Burrows A., 2013, Reviews of Modern Physics, 85, 245Burrows A., Dessart L., Livne E., Ott C. D., Murphy J., 2007, Astrophys. J., 664, 416Cabezón R. M., Pan K.-C., Liebendörfer M., Kuroda T., Ebinger K., Heinimann O., Perego A., Thielemann F.-K., 2018, A&A, 619, A118Cabral B., Leedom L., 1993, Computer Graphics (SIGGRAPH ’93 Proceedings), 27, 263–272Cantiello M., Mankovich C., Bildsten L., Christensen-Dalsgaard J., Paxton B., 2014, ApJ, 788, 93Chatzopoulos E., Frank J., Marcello D. C., Clayton G. C., 2020, ApJ, 896, 50MNRAS000
Akiyama S., Wheeler J. C., Meier D. L., Lichtenstadt I., 2003, ApJ, 584, 954Ardeljan N. V., Bisnovatyi-Kogan G. S., Moiseenko S. G., 2000, A&A, 355, 1181Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Barkat Z., Rakavy G., Sack N., 1967, Phys. Rev. Lett., 18, 379Bethe H. A., 1990, Reviews of Modern Physics, 62, 801Bisnovatyi-Kogan G. S., 1970, Azh, 47, 813Blondin J. M., Mezzacappa A., DeMarino C., 2003, ApJ, 584, 971Bruenn S. W., 1985, ApJS, 58, 771Bucciantini N., Quataert E., Metzger B. D., Thompson T. A., Arons J., Del Zanna L., 2009, MNRAS, 396, 2038Bugli M., Guilet J., Obergaulinger M., Cerdá-Durán P., Aloy M. A., 2020, MNRAS, 492, 58Buras R., Rampp M., Janka H. T., Kifonidis K., 2006, A&A, 447, 1049Burrows A., 2013, Reviews of Modern Physics, 85, 245Burrows A., Dessart L., Livne E., Ott C. D., Murphy J., 2007, Astrophys. J., 664, 416Cabezón R. M., Pan K.-C., Liebendörfer M., Kuroda T., Ebinger K., Heinimann O., Perego A., Thielemann F.-K., 2018, A&A, 619, A118Cabral B., Leedom L., 1993, Computer Graphics (SIGGRAPH ’93 Proceedings), 27, 263–272Cantiello M., Mankovich C., Bildsten L., Christensen-Dalsgaard J., Paxton B., 2014, ApJ, 788, 93Chatzopoulos E., Frank J., Marcello D. C., Clayton G. C., 2020, ApJ, 896, 50MNRAS000 , 1–22 (2020) J. Matsumoto et al.
Chen K.-J., Woosley S. E., Sukhbold T., 2016, ApJ, 832, 73Cherry J. F., Fuller G. M., Horiuchi S., Kotake K., Takiwaki T., Fischer T., 2020, Phys. Rev. D, 102, 023022Chevalier R. A., Irwin C. M., 2011, ApJ, 729, L6Colgate S. A., White R. H., 1966, ApJ, 143, 626Davidson P. A., 2001, An Introduction to MagnetohydrodynamicsDedner A., Kemm F., Kröner D., Munz C. D., Schnitzer T., Wesenberg M., 2002, Journal of Computational Physics, 175, 645Dessart L., Burrows A., Livne E., Ott C. D., 2008, ApJ, 673, L43Dessart L., Hillier D. J., Waldman R., Livne E., Blondin S., 2012, MNRAS, 426, L76Dolence J. C., Burrows A., Murphy J. W., Nordhaus J., 2013, ApJ, 765, 110Donati J. F., Babel J., Harries T. J., Howarth I. D., Petit P., Semel M., 2002, MNRAS, 333, 55Donati J. F., Howarth I. D., Bouret J. C., Petit P., Catala C., Landstreet J., 2006, MNRAS, 365, L6Duncan R. C., Thompson C., 1992, ApJ, 392, L9Ebinger K., Curtis S., Ghosh S., Fröhlich C., Hempel M., Perego A., Liebendörfer M., Thielemann F.-K., 2020, ApJ, 888, 91Einfeldt B., 1988, SIAM Journal on Numerical Analysis, 25, 294Endeve E., Cardall C. Y., Budiardja R. D., Mezzacappa A., 2010, ApJ, 713, 1219Endeve E., Cardall C. Y., Budiardja R. D., Beck S. W., Bejnood A., Toedte R. J., Mezzacappa A., Blondin J. M., 2012, ApJ, 751, 26Ertl T., Janka H. T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124Fernández R., Müller B., Foglizzo T., Janka H.-T., 2014, MNRAS, 440, 2763Ferrario L., Wickramasinghe D., 2006, MNRAS, 367, 1323Foglizzo T., Scheck L., Janka H.-T., 2006, ApJ, 652, 1436Foglizzo T., et al., 2015, Publ. Astron. Soc. Australia, 32, e009Fuller J., Lecoanet D., Cantiello M., Brown B., 2014, ApJ, 796, 17Gal-Yam A., 2012, Science, 337, 927Galloway D. J., Weiss N. O., 1981, ApJ, 243, 945Greiner J., et al., 2015, Nature, 523, 189Guilet J., Müller E., 2015, MNRAS, 450, 2153Guilet J., Müller E., Janka H.-T., 2015, MNRAS, 447, 3992Hanke F., Marek A., Müller B., Janka H.-T., 2012, ApJ, 755, 138Hanke F., Müller B., Wongwathanarat A., Marek A., Janka H.-T., 2013, ApJ, 770, 66Harada A., Nagakura H., Iwakami W., Okawa H., Furusawa S., Matsufuru H., Sumiyoshi K., Yamada S., 2019, ApJ, 872, 181Hayes J. C., Norman M. L., Fiedler R. A., Bordner J. O., Li P. S., Clark S. E., ud-Doula A., Mac Low M.-M., 2006, ApJS, 165, 188Heger A., Woosley S. E., Spruit H. C., 2005, ApJ, 626, 350Horowitz C. J., 2002, Phys. Rev. D, 65, 043001Hubrig S., Briquet M., Schöller M., De Cat P., Mathys G., Aerts C., 2006, MNRAS, 369, L61Iwakami W., Kotake K., Ohnishi N., Yamada S., Sawada K., 2008, ApJ, 678, 1207Iwakami W., Okawa H., Nagakura H., Harada A., Furusawa S., Sumiyoshi K., Matsufuru H., Yamada S., 2020, arXiv e-prints, p. arXiv:2004.02091Iwamoto K., et al., 1998, Nature, 395, 672Janka H.-T., 2012, Annual Review of Nuclear and Particle Science, 62, 407Just O., Obergaulinger M., Janka H. T., 2015, MNRAS, 453, 3386Just O., Bollig R., Janka H. T., Obergaulinger M., Glas R., Nagataki S., 2018, MNRAS, 481, 4786Kasen D., Bildsten L., 2010, ApJ, 717, 245Kaspi V. M., Beloborodov A. M., 2017, ARA&A, 55, 261Kim S.-s., Kim C., Rho O.-H., Kyu Hong S., 2003, Journal of Computational Physics, 185, 342Kotake K., Sawai H., Yamada S., Sato K., 2004, ApJ, 608, 391Kotake K., Takiwaki T., Suwa Y., Iwakami Nakano W., Kawagoe S., Masada Y., Fujimoto S.-i., 2012, Advances in Astronomy, 2012, 428757Kotake K., Takiwaki T., Fischer T., Nakamura K., Martínez-Pinedo G., 2018, ApJ, 853, 170Kouveliotou C., et al., 1998, Nature, 393, 235Kuroda T., Kotake K., Takiwaki T., 2012, ApJ, 755, 11Kuroda T., Takiwaki T., Kotake K., 2016, ApJS, 222, 20Kuroda T., Kotake K., Takiwaki T., Thielemann F.-K., 2018, MNRAS, 477, L80Kuroda T., Arcones A., Takiwaki T., Kotake K., 2020, ApJ, 896, 102Lai D., 2001, Reviews of Modern Physics, 73, 629Langer N., 2012, ARA&A, 50, 107Lattimer J. M., Swesty D. F., 1991, Nuclear Phys. A, 535, 331LeBlanc J. M., Wilson J. R., 1970, ApJ, 161, 541Lentz E. J., et al., 2015, ApJ, 807, L31Li S., Li H., 2003, Technical Report LA-UR-03-8925, Los Alamos National LabLiebendörfer M., Messer O. E. B., Mezzacappa A., Bruenn S. W., Cardall C. Y., Thielemann F. K., 2004, ApJS, 150, 263Liebendörfer M., Rampp M., Janka H. T., Mezzacappa A., 2005, ApJ, 620, 840Liebendörfer M., Whitehouse S. C., Fischer T., 2009, ApJ, 698, 1174Londrillo P., Del Zanna L., 2000, ApJ, 530, 508Maeda K., Nomoto K., 2003, ApJ, 598, 1163Marek A., Dimmelmeier H., Janka H. T., Müller E., Buras R., 2006, A&A, 445, 273Masada Y., Sano T., Takabe H., 2006, ApJ, 641, 447Masada Y., Sano T., Shibata K., 2007, ApJ, 655, 447Masada Y., Takiwaki T., Kotake K., Sano T., 2012, ApJ, 759, 110Masada Y., Takiwaki T., Kotake K., 2015, ApJ, 798, L22Masada Y., Kotake K., Takiwaki T., Yamamoto N., 2018, Phys. Rev. D, 98, 083018 MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe Masada Y., Takiwaki T., Kotake K., 2020, arXiv e-prints, p. arXiv:2001.08452Matsumoto Y., et al., 2019, PASJ, 71, 83Mazzali P. A., McFadyen A. I., Woosley S. E., Pian E., Tanaka M., 2014, MNRAS, 443, 67Meier D. L., Epstein R. I., Arnett W. D., Schramm D. N., 1976, ApJ, 204, 869Melson T., Kresse D., Janka H.-T., 2020, ApJ, 891, 27Metzger B. D., Giannios D., Thompson T. A., Bucciantini N., Quataert E., 2011, MNRAS, 413, 2031Metzger B. D., Margalit B., Kasen D., Quataert E., 2015, MNRAS, 454, 3311Mezzacappa A., 2005, Annual Review of Nuclear and Particle Science, 55, 467Mignone A., 2014, Journal of Computational Physics, 270, 784Miyoshi T., Kusano K., 2005, in AGU Fall Meeting Abstracts. pp SM51B–1295Moiseenko S. G., Bisnovatyi-Kogan G. S., Ardeljan N. V., 2006, MNRAS, 370, 501Moriya T. J., Blinnikov S. I., Tominaga N., Yoshida N., Tanaka M., Maeda K., Nomoto K., 2013, MNRAS, 428, 1020Moriya T. J., Sorokina E. I., Chevalier R. A., 2018, Space Sci. Rev., 214, 59Mösta P., et al., 2014, ApJ, 785, L29Mösta P., Ott C. D., Radice D., Roberts L. F., Schnetter E., Haas R., 2015, Nature, 528, 376Müller B., 2015, MNRAS, 453, 287Müller B., 2020, arXiv e-prints, p. arXiv:2006.05083Müller E., Hillebrandt W., 1979, A&A, 80, 147Müller B., Varma V., 2020, arXiv e-prints, p. arXiv:2007.04775Müller B., Janka H.-T., Dimmelmeier H., 2010, ApJS, 189, 104Nagakura H., Sumiyoshi K., Yamada S., 2014, ApJS, 214, 16Nagakura H., Iwakami W., Furusawa S., Sumiyoshi K., Yamada S., Matsufuru H., Imakura A., 2017, ApJS, 229, 42Nagakura H., Burrows A., Radice D., Vartanyan D., 2019a, MNRAS, 490, 4622Nagakura H., Burrows A., Radice D., Vartanyan D., 2019b, MNRAS, 490, 4622Nagakura H., Sumiyoshi K., Yamada S., 2019c, ApJ, 878, 160Nakamura K., Takiwaki T., Kuroda T., Kotake K., 2015, PASJ, 67, 107Nakamura K., Takiwaki T., Kotake K., 2019, PASJ, 71, 98Nakano T., Murakami H., Furuta Y., Enoto T., Masuyama M., Shigeyama T., Makishima K., 2017, PASJ, 69, 40Nicholl M., et al., 2013, Nature, 502, 346Nomoto K., Tanaka M., Tominaga N., Maeda K., 2010, New Astron. Rev., 54, 191Nordhaus J., Burrows A., Almgren A., Bell J., 2010, ApJ, 720, 694O’Connor E., 2015, ApJS, 219, 24O’Connor E. P., Couch S. M., 2018, ApJ, 865, 81O’Connor E., Ott C. D., 2013, ApJ, 762, 126O’Connor E., et al., 2018, Journal of Physics G Nuclear Physics, 45, 104001Obergaulinger M., Aloy M. Á., 2017, MNRAS, 469, L43Obergaulinger M., Aloy M. Á., 2020, MNRAS, 492, 4613Obergaulinger M., Aloy M. A., Müller E., 2006a, A&A, 450, 1107Obergaulinger M., Aloy M. A., Dimmelmeier H., Müller E., 2006b, A&A, 457, 209Obergaulinger M., Cerdá-Durán P., Müller E., Aloy M. A., 2009, A&A, 498, 241Obergaulinger M., Janka H. T., Aloy M. A., 2014, MNRAS, 445, 3169Obergaulinger M., Just O., Aloy M. A., 2018, Journal of Physics G Nuclear Physics, 45, 084001Ott C. D., Burrows A., Thompson T. A., Livne E., Walder R., 2006, ApJS, 164, 130Ott C. D., Roberts L. F., da Silva Schneider A., Fedrow J. M., Haas R., Schnetter E., 2018, ApJ, 855, L3Paczynski B., 1992, Acta Astron., 42, 145Pan K.-C., Liebendörfer M., Couch S. M., Thielemann F.-K., 2018, ApJ, 857, 13Rakavy G., Shaviv G., 1967, ApJ, 148, 803Ramírez-Agudelo O. H., et al., 2013, A&A, 560, A29Raynaud R., Guilet J., Janka H.-T., Gastine T., 2020, Science Advances, 6, eaay2732Reboul-Salze A., Guilet J., Raynaud R., Bugli M., 2020, arXiv e-prints, p. arXiv:2005.03567Rembiasz T., Guilet J., Obergaulinger M., Cerdá-Durán P., Aloy M. A., Müller E., 2016, MNRAS, 460, 3316Roberts L. F., Ott C. D., Haas R., O’Connor E. P., Diener P., Schnetter E., 2016, ApJ, 831, 98Sasaki H., Kajino T., Takiwaki T., Hayakawa T., Balantekin A. B., Pehlivan Y., 2017, Phys. Rev. D, 96, 043013Sasaki H., Takiwaki T., Kawagoe S., Horiuchi S., Ishidoshiro K., 2020, Phys. Rev. D, 101, 063027Sawai H., Yamada S., 2014, ApJ, 784, L10Sawai H., Yamada S., 2016, ApJ, 817, 153Sawai H., Kotake K., Yamada S., 2005, ApJ, 631, 446Sawai H., Yamada S., Kotake K., Suzuki H., 2013a, ApJ, 764, 10Sawai H., Yamada S., Suzuki H., 2013b, ApJ, 770, L19Shibata M., Liu Y. T., Shapiro S. L., Stephens B. C., 2006, Phys. Rev. D, 74, 104026Shibata M., Kiuchi K., Sekiguchi Y., Suwa Y., 2011, Progress of Theoretical Physics, 125, 1255Shigeyama T., Kashiyama K., 2018, PASJ, 70, 107Shultz M. E., et al., 2018, MNRAS, 475, 5144Skinner M. A., Ostriker E. C., 2010, ApJS, 188, 290Skinner M. A., Burrows A., Dolence J. C., 2016, ApJ, 831, 81Soderberg A. M., et al., 2006, Nature, 442, 1014Sorokina E., Blinnikov S., Nomoto K., Quimby R., Tolstov A., 2016, ApJ, 829, 17Sotani H., Takiwaki T., 2016, Phys. Rev. D, 94, 044043MNRAS000
Chen K.-J., Woosley S. E., Sukhbold T., 2016, ApJ, 832, 73Cherry J. F., Fuller G. M., Horiuchi S., Kotake K., Takiwaki T., Fischer T., 2020, Phys. Rev. D, 102, 023022Chevalier R. A., Irwin C. M., 2011, ApJ, 729, L6Colgate S. A., White R. H., 1966, ApJ, 143, 626Davidson P. A., 2001, An Introduction to MagnetohydrodynamicsDedner A., Kemm F., Kröner D., Munz C. D., Schnitzer T., Wesenberg M., 2002, Journal of Computational Physics, 175, 645Dessart L., Burrows A., Livne E., Ott C. D., 2008, ApJ, 673, L43Dessart L., Hillier D. J., Waldman R., Livne E., Blondin S., 2012, MNRAS, 426, L76Dolence J. C., Burrows A., Murphy J. W., Nordhaus J., 2013, ApJ, 765, 110Donati J. F., Babel J., Harries T. J., Howarth I. D., Petit P., Semel M., 2002, MNRAS, 333, 55Donati J. F., Howarth I. D., Bouret J. C., Petit P., Catala C., Landstreet J., 2006, MNRAS, 365, L6Duncan R. C., Thompson C., 1992, ApJ, 392, L9Ebinger K., Curtis S., Ghosh S., Fröhlich C., Hempel M., Perego A., Liebendörfer M., Thielemann F.-K., 2020, ApJ, 888, 91Einfeldt B., 1988, SIAM Journal on Numerical Analysis, 25, 294Endeve E., Cardall C. Y., Budiardja R. D., Mezzacappa A., 2010, ApJ, 713, 1219Endeve E., Cardall C. Y., Budiardja R. D., Beck S. W., Bejnood A., Toedte R. J., Mezzacappa A., Blondin J. M., 2012, ApJ, 751, 26Ertl T., Janka H. T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124Fernández R., Müller B., Foglizzo T., Janka H.-T., 2014, MNRAS, 440, 2763Ferrario L., Wickramasinghe D., 2006, MNRAS, 367, 1323Foglizzo T., Scheck L., Janka H.-T., 2006, ApJ, 652, 1436Foglizzo T., et al., 2015, Publ. Astron. Soc. Australia, 32, e009Fuller J., Lecoanet D., Cantiello M., Brown B., 2014, ApJ, 796, 17Gal-Yam A., 2012, Science, 337, 927Galloway D. J., Weiss N. O., 1981, ApJ, 243, 945Greiner J., et al., 2015, Nature, 523, 189Guilet J., Müller E., 2015, MNRAS, 450, 2153Guilet J., Müller E., Janka H.-T., 2015, MNRAS, 447, 3992Hanke F., Marek A., Müller B., Janka H.-T., 2012, ApJ, 755, 138Hanke F., Müller B., Wongwathanarat A., Marek A., Janka H.-T., 2013, ApJ, 770, 66Harada A., Nagakura H., Iwakami W., Okawa H., Furusawa S., Matsufuru H., Sumiyoshi K., Yamada S., 2019, ApJ, 872, 181Hayes J. C., Norman M. L., Fiedler R. A., Bordner J. O., Li P. S., Clark S. E., ud-Doula A., Mac Low M.-M., 2006, ApJS, 165, 188Heger A., Woosley S. E., Spruit H. C., 2005, ApJ, 626, 350Horowitz C. J., 2002, Phys. Rev. D, 65, 043001Hubrig S., Briquet M., Schöller M., De Cat P., Mathys G., Aerts C., 2006, MNRAS, 369, L61Iwakami W., Kotake K., Ohnishi N., Yamada S., Sawada K., 2008, ApJ, 678, 1207Iwakami W., Okawa H., Nagakura H., Harada A., Furusawa S., Sumiyoshi K., Matsufuru H., Yamada S., 2020, arXiv e-prints, p. arXiv:2004.02091Iwamoto K., et al., 1998, Nature, 395, 672Janka H.-T., 2012, Annual Review of Nuclear and Particle Science, 62, 407Just O., Obergaulinger M., Janka H. T., 2015, MNRAS, 453, 3386Just O., Bollig R., Janka H. T., Obergaulinger M., Glas R., Nagataki S., 2018, MNRAS, 481, 4786Kasen D., Bildsten L., 2010, ApJ, 717, 245Kaspi V. M., Beloborodov A. M., 2017, ARA&A, 55, 261Kim S.-s., Kim C., Rho O.-H., Kyu Hong S., 2003, Journal of Computational Physics, 185, 342Kotake K., Sawai H., Yamada S., Sato K., 2004, ApJ, 608, 391Kotake K., Takiwaki T., Suwa Y., Iwakami Nakano W., Kawagoe S., Masada Y., Fujimoto S.-i., 2012, Advances in Astronomy, 2012, 428757Kotake K., Takiwaki T., Fischer T., Nakamura K., Martínez-Pinedo G., 2018, ApJ, 853, 170Kouveliotou C., et al., 1998, Nature, 393, 235Kuroda T., Kotake K., Takiwaki T., 2012, ApJ, 755, 11Kuroda T., Takiwaki T., Kotake K., 2016, ApJS, 222, 20Kuroda T., Kotake K., Takiwaki T., Thielemann F.-K., 2018, MNRAS, 477, L80Kuroda T., Arcones A., Takiwaki T., Kotake K., 2020, ApJ, 896, 102Lai D., 2001, Reviews of Modern Physics, 73, 629Langer N., 2012, ARA&A, 50, 107Lattimer J. M., Swesty D. F., 1991, Nuclear Phys. A, 535, 331LeBlanc J. M., Wilson J. R., 1970, ApJ, 161, 541Lentz E. J., et al., 2015, ApJ, 807, L31Li S., Li H., 2003, Technical Report LA-UR-03-8925, Los Alamos National LabLiebendörfer M., Messer O. E. B., Mezzacappa A., Bruenn S. W., Cardall C. Y., Thielemann F. K., 2004, ApJS, 150, 263Liebendörfer M., Rampp M., Janka H. T., Mezzacappa A., 2005, ApJ, 620, 840Liebendörfer M., Whitehouse S. C., Fischer T., 2009, ApJ, 698, 1174Londrillo P., Del Zanna L., 2000, ApJ, 530, 508Maeda K., Nomoto K., 2003, ApJ, 598, 1163Marek A., Dimmelmeier H., Janka H. T., Müller E., Buras R., 2006, A&A, 445, 273Masada Y., Sano T., Takabe H., 2006, ApJ, 641, 447Masada Y., Sano T., Shibata K., 2007, ApJ, 655, 447Masada Y., Takiwaki T., Kotake K., Sano T., 2012, ApJ, 759, 110Masada Y., Takiwaki T., Kotake K., 2015, ApJ, 798, L22Masada Y., Kotake K., Takiwaki T., Yamamoto N., 2018, Phys. Rev. D, 98, 083018 MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe Masada Y., Takiwaki T., Kotake K., 2020, arXiv e-prints, p. arXiv:2001.08452Matsumoto Y., et al., 2019, PASJ, 71, 83Mazzali P. A., McFadyen A. I., Woosley S. E., Pian E., Tanaka M., 2014, MNRAS, 443, 67Meier D. L., Epstein R. I., Arnett W. D., Schramm D. N., 1976, ApJ, 204, 869Melson T., Kresse D., Janka H.-T., 2020, ApJ, 891, 27Metzger B. D., Giannios D., Thompson T. A., Bucciantini N., Quataert E., 2011, MNRAS, 413, 2031Metzger B. D., Margalit B., Kasen D., Quataert E., 2015, MNRAS, 454, 3311Mezzacappa A., 2005, Annual Review of Nuclear and Particle Science, 55, 467Mignone A., 2014, Journal of Computational Physics, 270, 784Miyoshi T., Kusano K., 2005, in AGU Fall Meeting Abstracts. pp SM51B–1295Moiseenko S. G., Bisnovatyi-Kogan G. S., Ardeljan N. V., 2006, MNRAS, 370, 501Moriya T. J., Blinnikov S. I., Tominaga N., Yoshida N., Tanaka M., Maeda K., Nomoto K., 2013, MNRAS, 428, 1020Moriya T. J., Sorokina E. I., Chevalier R. A., 2018, Space Sci. Rev., 214, 59Mösta P., et al., 2014, ApJ, 785, L29Mösta P., Ott C. D., Radice D., Roberts L. F., Schnetter E., Haas R., 2015, Nature, 528, 376Müller B., 2015, MNRAS, 453, 287Müller B., 2020, arXiv e-prints, p. arXiv:2006.05083Müller E., Hillebrandt W., 1979, A&A, 80, 147Müller B., Varma V., 2020, arXiv e-prints, p. arXiv:2007.04775Müller B., Janka H.-T., Dimmelmeier H., 2010, ApJS, 189, 104Nagakura H., Sumiyoshi K., Yamada S., 2014, ApJS, 214, 16Nagakura H., Iwakami W., Furusawa S., Sumiyoshi K., Yamada S., Matsufuru H., Imakura A., 2017, ApJS, 229, 42Nagakura H., Burrows A., Radice D., Vartanyan D., 2019a, MNRAS, 490, 4622Nagakura H., Burrows A., Radice D., Vartanyan D., 2019b, MNRAS, 490, 4622Nagakura H., Sumiyoshi K., Yamada S., 2019c, ApJ, 878, 160Nakamura K., Takiwaki T., Kuroda T., Kotake K., 2015, PASJ, 67, 107Nakamura K., Takiwaki T., Kotake K., 2019, PASJ, 71, 98Nakano T., Murakami H., Furuta Y., Enoto T., Masuyama M., Shigeyama T., Makishima K., 2017, PASJ, 69, 40Nicholl M., et al., 2013, Nature, 502, 346Nomoto K., Tanaka M., Tominaga N., Maeda K., 2010, New Astron. Rev., 54, 191Nordhaus J., Burrows A., Almgren A., Bell J., 2010, ApJ, 720, 694O’Connor E., 2015, ApJS, 219, 24O’Connor E. P., Couch S. M., 2018, ApJ, 865, 81O’Connor E., Ott C. D., 2013, ApJ, 762, 126O’Connor E., et al., 2018, Journal of Physics G Nuclear Physics, 45, 104001Obergaulinger M., Aloy M. Á., 2017, MNRAS, 469, L43Obergaulinger M., Aloy M. Á., 2020, MNRAS, 492, 4613Obergaulinger M., Aloy M. A., Müller E., 2006a, A&A, 450, 1107Obergaulinger M., Aloy M. A., Dimmelmeier H., Müller E., 2006b, A&A, 457, 209Obergaulinger M., Cerdá-Durán P., Müller E., Aloy M. A., 2009, A&A, 498, 241Obergaulinger M., Janka H. T., Aloy M. A., 2014, MNRAS, 445, 3169Obergaulinger M., Just O., Aloy M. A., 2018, Journal of Physics G Nuclear Physics, 45, 084001Ott C. D., Burrows A., Thompson T. A., Livne E., Walder R., 2006, ApJS, 164, 130Ott C. D., Roberts L. F., da Silva Schneider A., Fedrow J. M., Haas R., Schnetter E., 2018, ApJ, 855, L3Paczynski B., 1992, Acta Astron., 42, 145Pan K.-C., Liebendörfer M., Couch S. M., Thielemann F.-K., 2018, ApJ, 857, 13Rakavy G., Shaviv G., 1967, ApJ, 148, 803Ramírez-Agudelo O. H., et al., 2013, A&A, 560, A29Raynaud R., Guilet J., Janka H.-T., Gastine T., 2020, Science Advances, 6, eaay2732Reboul-Salze A., Guilet J., Raynaud R., Bugli M., 2020, arXiv e-prints, p. arXiv:2005.03567Rembiasz T., Guilet J., Obergaulinger M., Cerdá-Durán P., Aloy M. A., Müller E., 2016, MNRAS, 460, 3316Roberts L. F., Ott C. D., Haas R., O’Connor E. P., Diener P., Schnetter E., 2016, ApJ, 831, 98Sasaki H., Kajino T., Takiwaki T., Hayakawa T., Balantekin A. B., Pehlivan Y., 2017, Phys. Rev. D, 96, 043013Sasaki H., Takiwaki T., Kawagoe S., Horiuchi S., Ishidoshiro K., 2020, Phys. Rev. D, 101, 063027Sawai H., Yamada S., 2014, ApJ, 784, L10Sawai H., Yamada S., 2016, ApJ, 817, 153Sawai H., Kotake K., Yamada S., 2005, ApJ, 631, 446Sawai H., Yamada S., Kotake K., Suzuki H., 2013a, ApJ, 764, 10Sawai H., Yamada S., Suzuki H., 2013b, ApJ, 770, L19Shibata M., Liu Y. T., Shapiro S. L., Stephens B. C., 2006, Phys. Rev. D, 74, 104026Shibata M., Kiuchi K., Sekiguchi Y., Suwa Y., 2011, Progress of Theoretical Physics, 125, 1255Shigeyama T., Kashiyama K., 2018, PASJ, 70, 107Shultz M. E., et al., 2018, MNRAS, 475, 5144Skinner M. A., Ostriker E. C., 2010, ApJS, 188, 290Skinner M. A., Burrows A., Dolence J. C., 2016, ApJ, 831, 81Soderberg A. M., et al., 2006, Nature, 442, 1014Sorokina E., Blinnikov S., Nomoto K., Quimby R., Tolstov A., 2016, ApJ, 829, 17Sotani H., Takiwaki T., 2016, Phys. Rev. D, 94, 044043MNRAS000 , 1–22 (2020) J. Matsumoto et al.
Sotani H., Takiwaki T., 2020, Phys. Rev. D, 102, 023028Stone J. M., Gardiner T., 2009, New Astron., 14, 139Stone J. M., Norman M. L., 1992, ApJS, 80, 753Stone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B., 2008, ApJS, 178, 137Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38Sumiyoshi K., Yamada S., 2012, ApJS, 199, 17Sumiyoshi K., Yamada S., Suzuki H., Shen H., Chiba S., Toki H., 2005, ApJ, 629, 922Summa A., Hanke F., Janka H.-T., Melson T., Marek A., Müller B., 2016, ApJ, 825, 6Suwa Y., Takiwaki T., Kotake K., Sato K., 2007, PASJ, 59, 771Suwa Y., Kotake K., Takiwaki T., Whitehouse S. C., LiebendÃ-rfer M., Sato K., 2010, PASJ, 62, L49Suwa Y., Yamada S., Takiwaki T., Kotake K., 2016, ApJ, 816, 43Symbalisty E. M. D., 1984, ApJ, 285, 729Takiwaki T., Kotake K., 2011, ApJ, 743, 30Takiwaki T., Kotake K., Sato K., 2009, ApJ, 691, 1360Takiwaki T., Kotake K., Suwa Y., 2012, ApJ, 749, 98Takiwaki T., Kotake K., Suwa Y., 2014, ApJ, 786, 83Takiwaki T., Kotake K., Suwa Y., 2016, MNRAS, 461, L112Terreran G., et al., 2017, Nature Astronomy, 1, 713Thompson C., 1994, MNRAS, 270, 480Thompson C., Duncan R. C., 1993, ApJ, 408, 194Tominaga N., 2009, ApJ, 690, 526Usov V. V., 1992, Nature, 357, 472Vartanyan D., Burrows A., Radice D., Skinner M. A., Dolence J., 2019, MNRAS, 482, 351Vink J., Kuiper L., 2006, MNRAS, 370, L14Wade G. A., MiMeS Collaboration 2015, in Balega Y. Y., Romanyuk I. I., Kudryavtsev D. O., eds, Astronomical Society of the Pacific Conference Series Vol.494, Physics and Evolution of Magnetic and Related Stars. p. 30 ( arXiv:1411.3604 )Wang S. Q., Wang L. J., Dai Z. G., Wu X. F., 2015, ApJ, 799, 107Wheeler J. C., Meier D. L., Wilson J. R., 2002, ApJ, 568, 807Winteler C., Käppeli R., Perego A., Arcones A., Vasset N., Nishimura N., Liebendörfer M., Thielemann F. K., 2012, ApJ, 750, L22Wongwathanarat A., Hammer N. J., Müller E., 2010, A&A, 514, A48Woosley S. E., 2010, ApJ, 719, L204Woosley S. E., Bloom J. S., 2006, ARA&A, 44, 507Woosley S. E., Heger A., 2007, Phys. Rep., 442, 269Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015Yamamoto N., 2016, Phys. Rev. D, 93, 065017Zaizen M., Cherry J. F., Takiwaki T., Horiuchi S., Kotake K., Umeda H., Yoshida T., 2020, J. Cosmology Astropart. Phys., 2020, 011Zhang M., Feng X., 2016, Frontiers in Astronomy and Space Sciences, 3, 6
APPENDIX A: FULL 3D MHD EQUATIONS IN SPHERICAL COORDINATES
The explicit formulae of the ideal MHD equations in spherical coordinates, neglecting the contribution of the gravity and the interaction withneutrinos, are ∂ρ∂ t + r ∂ ( r ρ v r ) ∂ r + r sin θ ∂ ( ρ v θ sin θ ) ∂θ + r sin θ ∂ ( ρ v φ ) ∂φ = , (A1) ∂ ( ρ v r ) ∂ t + r ∂ [ r ( ρ v r v r + P t − B r B r )] ∂ r + r sin θ ∂ [( ρ v r v θ − B r B θ ) sin θ ] ∂θ + r sin θ ∂ ( ρ v r v φ − B r B φ ) ∂φ = ρ ( v θ v θ + v φ v φ ) + P t − B θ B θ − B φ B φ r , (A2) ∂ ( ρ v θ ) ∂ t + r ∂ [ r ( ρ v θ v r − B θ B r )] ∂ r + r sin θ ∂ [( ρ v θ v θ + P t − B θ B θ ) sin θ ] ∂θ + r sin θ ∂ ( ρ v θ v φ − B θ B φ ) ∂φ = ( ρ v φ v φ + P t − B φ B φ ) cot θ r − ρ v θ v r − B θ B r r , (A3) ∂ ( ρ v φ ) ∂ t + r ∂ [ r ( ρ v φ v r − B φ B r ) ∂ r + r sin θ ∂ [( ρ v φ v θ − B φ B θ ) sin θ ] ∂θ + r sin θ ∂ ( ρ v φ v φ + P t − B φ B φ ) ∂φ = − ( ρ v φ v θ − B φ B θ ) cot θ r − ρ v φ v r − B φ B r r , (A4) ∂ e ∂ t + r ∂ [ r ( e + P t ) v r − r ( v · B ) B r ]] ∂ r + r sin θ ∂ [( e + P t ) v θ sin θ − ( v · B ) B θ sin θ ] ∂θ + r sin θ ∂ [( e + P t ) v φ − ( v · B ) B φ ] ∂φ = , (A5) ∂ B r ∂ t + r ∂ ( r ψ ) ∂ r + r sin θ ∂ [( B r v θ − v r B θ ) sin θ ] ∂θ + r sin θ ∂ ( B r v φ − v r B φ ) ∂φ = ψ r , (A6) ∂ B θ ∂ t + r ∂ [ r ( B θ v r − v θ B r )] ∂ r + r sin θ ∂ ( ψ sin θ ) ∂θ + r sin θ ∂ ( B θ v φ − v θ B φ ) ∂φ = ψ cot θ r , (A7) ∂ B φ ∂ t + r ∂ [ r ( B φ v r − v φ B r )] ∂ r + r ∂ ( B φ v θ − v φ B θ ) ∂θ + r sin θ ∂ψ∂φ = , (A8) MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe Table A1.
Summary of the discretization method and the reconstruction of physical variables at the cell surface for conservation equations.discretization method equation position of term discretization formula reconstruction of physical variablesfinite volume (A1), (A2), (A3), (A4) (A5), (A9) 1st term (B8)finite volume (A1), (A2), (A3), (A4) (A5), (A6), (A9) 2nd term (B3) (C2), (C3) with (C4), (C10)finite volume (A1), (A2), (A3), (A4) (A5), (A7), (A9) 3rd term (B4) (C15), (C16)finite volume (A1), (A2), (A3), (A4) (A5), (A8), (A9) 4th term (B5) (C25), (C26)finite area (A6), (A7), (A8) 1st term (B21)finite area (A6) 3rd term (B15) (C15), (C16)finite area (A6) 4th term (B16) (C25), (C26)finite area (A7) 2nd term (B17) (C2), (C3) with (C30), (C31)finite area (A7) 4th term (B18) (C25), (C26)finite area (A8) 2nd term (B19) (C2), (C3) with (C30), (C31)finite area (A8) 3rd term (B20) (C25) a , (C26) aa Spacial index k is replaced by j . ∂ψ∂ t + r ∂ ( r c h B r ) ∂ r + r sin θ ∂ [ c h B θ sin θ ] ∂θ + r sin θ ∂ ( c h B φ ) ∂φ = − c h c p ψ . (A9)These are the continuity equation (A1), the momentum conservation equations (A2)–(A4), the energy conservation equation (A5) and theinduction equations (A6)–(A8). The divergence cleaning method (Dedner et al. 2002) is implemented in our code to maintain the numericalerrors of solenoidal property of the magnetic field within minimal levels. The average relative divergence error estimated by Error3( B ) ave of Zhang & Feng (2016) is less than 1% in this work. Equation (A9) and ψ in induction equations (A6)–(A8) are related to the divergencecleaning method. The HLLD scheme (Miyoshi & Kusano 2005) is used to solve equations (A1)–(A8) in a conservative form. To solveequation (A9), the HLLE scheme (Einfeldt 1988) is used.In addition to the baryon number conservation (equation A1), we evolve the equation regarding the lepton number conservation (seeequation 6), whose discretization and reconstruction method is the same as that of equation (A1). In the IDSA scheme, the energy equationfor the trapped neutrinos needs to be evolved (see equation 7). Though the treatment of the advection term of equation (7) is the same as thatof equation (A1), the left-hand-side of equation (7) includes the term of p ν ∇ · v with p ν = ρ Z m the neutrino pressure, which should be treatedas a source term. Then the explicit form of equation (6) in spherical coordinates reads ∂ρ Z m ∂ t + r ∂ ( r ρ Z m v r ) ∂ r + r sin θ ∂ ( ρ Z m v θ sin θ ) ∂θ + r sin θ ∂ ( ρ Z m v φ ) ∂φ = − ρ Z m r ∂ ( r v r ) ∂ r − ρ Z m r sin θ ∂ ( v θ sin θ ) ∂θ − ρ Z m r sin θ ∂ v φ ∂φ , (A10)noting again that in this Appendix the neutrino-matter interaction term is dropped for the sake of brevity. APPENDIX B: FINITE VOLUME AND AREA DISCRETIZATION IN SPHERICAL COORDINATES
The governing equations for the conservative variables are evolved by the finite volume and area methods in our code. Both the finite volumeand area methods are used for the induction equations, while only the finite volume method is used for other equations.In the spherical coordinates, the conservation equations except the induction equations are simply described as follows; ∂ Q ∂ t + r ∂ ( r F r ) ∂ r + r sin θ ∂ ( F θ sin θ ) ∂θ + r sin θ ∂ F φ ∂φ = S , (B1)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. Followingthe 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 thevolume average of Q over the cell volume,¯ Q i , j , k = (cid:82) i + / i − / (cid:82) j + / j − / (cid:82) k + / k − / Q ( r , θ, φ ) r sin θ drd θ d φ (cid:82) i + / i − / (cid:82) j + / j − / (cid:82) k + / k − / r sin θ drd θ d φ , (B2)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-sideof equation (B1) are discretized as follows; 1 r ∂ ( r F r ) ∂ r ∼ r i + / F r , i + / , j , k − r i − / F r , i − / , j , k r i + / / − r i − / / , (B3)1 r sin θ ∂ ( F θ sin θ ) ∂θ ∼ r i + / / − r i − / / r i + / / − r i − / / F θ, i , j + / , k sin θ j + / − F θ, i , j − / , k sin θ j − / cos θ j − / − cos θ j + / , (B4) MNRAS000
The governing equations for the conservative variables are evolved by the finite volume and area methods in our code. Both the finite volumeand area methods are used for the induction equations, while only the finite volume method is used for other equations.In the spherical coordinates, the conservation equations except the induction equations are simply described as follows; ∂ Q ∂ t + r ∂ ( r F r ) ∂ r + r sin θ ∂ ( F θ sin θ ) ∂θ + r sin θ ∂ F φ ∂φ = S , (B1)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. Followingthe 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 thevolume average of Q over the cell volume,¯ Q i , j , k = (cid:82) i + / i − / (cid:82) j + / j − / (cid:82) k + / k − / Q ( r , θ, φ ) r sin θ drd θ d φ (cid:82) i + / i − / (cid:82) j + / j − / (cid:82) k + / k − / r sin θ drd θ d φ , (B2)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-sideof equation (B1) are discretized as follows; 1 r ∂ ( r F r ) ∂ r ∼ r i + / F r , i + / , j , k − r i − / F r , i − / , j , k r i + / / − r i − / / , (B3)1 r sin θ ∂ ( F θ sin θ ) ∂θ ∼ r i + / / − r i − / / r i + / / − r i − / / F θ, i , j + / , k sin θ j + / − F θ, i , j − / , k sin θ j − / cos θ j − / − cos θ j + / , (B4) MNRAS000 , 1–22 (2020) J. Matsumoto et al. r sin θ ∂ F φ ∂φ ∼ r i + / / − r i − / / r i + / / − r i − / / θ j + / − θ j − / cos θ j − / − cos θ j + / F φ, i , j , k + / − F φ, i , j , k − / φ k + / − φ k + / , (B5)where F r , i ± / , j , k , F θ, i , j ± / , k and F φ, i , j , k ± / are the numerical flues at the cell surface in each direction (or the interpolated velocity for the sourceterms 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 ofthe discrete conservation equations are replaced like below;cot θ ∼ sin θ j + / − sin θ j − / cos θ j − / − cos θ j + / , (B6)1 r ∼ r i + / / − r i − / / r i + / / − r i − / / . (B7)Physical variables in the source terms are also evaluated by the cell-volume-averaged values. The time derivative of equation (B1) is dis-cretized forward in time: ∂ Q ∂ t ∼ ¯ Q n + i , j , k − ¯ Q ni , j , k ∆ t , (B8)where the superscript n stands for the number of time steps and ∆ t is a step size. This formulation leads to the first-order accuracy in time,which will be extended in higher-order accuracy (Matsumoto et al. in preparation).In the induction equations (A6)–(A8), the terms including ψ are also discretized based on the finite volume method. However, consid-ering the magnetic flux conservation, other terms should be better discretized based on the finite area method, namely, by integrating theseterms over the cell area using the Stokes’ theorem. The terms discretized by the finite area method are listed in Table A1. The first termsin the left-hand-side of the induction equations are time derivatives of each component of the magnetic field. Here, we introduce the areaaverage of the variables over the cell areas defined at the cell center,¯ B r , i , j , k = (cid:82) j + / j − / (cid:82) k + / k − / B r ( θ, φ ) r sin θ d θ d φ (cid:82) j + / j − / (cid:82) k + / k − / r sin θ d θ d φ , (B9)¯ B θ, i , j , k = (cid:82) i + / i − / (cid:82) k + / k − / B θ ( r , φ ) r sin θ drd φ (cid:82) i + / i − / (cid:82) k + / k − / r sin θ drd φ , (B10)¯ B φ, i , j , k = (cid:82) i + / i − / (cid:82) j + / j − / B φ ( r , θ ) rdrd θ (cid:82) i + / i − / (cid:82) j + / j − / rdrd θ , (B11)respectively. The flux in the conservation form of the induction equation is related to the electric field defined as: E r = B θ v φ − v θ B φ c , (B12) E θ = B φ v r − v φ B r c , (B13) E φ = B r v θ − v r B θ c , (B14)where c is the speed of light. The flux terms of the induction equations are then discretized as follows;1 r sin θ ∂ ( cE φ sin θ ) ∂θ ∼ r i + / − r i − / r i + / / − r i − / / cE φ, i , j + / , k sin θ j + / − cE φ, i , j − / , k sin θ j − / cos θ j − / − cos θ j + / , (B15)1 r sin θ ∂ ( − cE θ ) ∂φ ∼ r i + / − r i − / r i + / / − r i − / / θ j + / − θ j − / cos θ j − / − cos θ j + / − cE θ, i , j , k + / + cE θ, i , j , k − / φ k + / − φ k − / , (B16)1 r ∂ ( − rcE φ ) ∂ r ∼ − r i + / cE φ, i + / , j , k + r i − / cE φ, i − / , j , k r i + / / − r i − / / , (B17)1 r sin θ ∂ cE r ∂φ ∼ r i + / − r i − / r i + / / − r i − / / θ j + / − θ j − / cos θ j − / − cos θ j + / cE r , i , j , k + / − cE r , i , j , k − / φ k + / − φ k − / , (B18)1 r ∂ ( rcE θ ) ∂ r ∼ r i + / cE θ, i + / , j , k − r i − / cE θ, i − / , j , k r i + / / − r i − / / , (B19)1 r ∂ ( − cE r ) ∂θ ∼ r i + / − r i − / r i + / / − r i − / / − cE r , i , j + / , k + cE r , i , j − / , k θ j + / − θ j + / . (B20)The time derivative of the induction equation is also discretized forward in time: ∂ B s ∂ t ∼ ¯ B n + s , i , j , k − ¯ B ns , i , j , k ∆ t , (B21) MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe where s represents the direction ( r , θ, φ ).The numerical fluxes at the cell surface are normally computed by solving a Riemann problem between the discontinuous states acrossthe cell interfaces. In our code, an approximate Riemann solver (HLLD; Miyoshi & Kusano 2005) is used to estimate the numerical fluxes.On the choice of the parameters in equation (A9), c h represents the propagation speed of the numerical error of the deviation from ∇ · B =
0. For the value, the largest eigenvalue in the global computational domain is often taken (e.g., Zhang & Feng 2016): c h = max i , j , k (cid:16) c s + | v r | , c s + | v θ | , c s + | v φ | (cid:17) , (B22)where c s is the total sound speed that includes the e ff ect of the magnetic fields. Note that we can also evaluate it from the Courant–Friedrichs–Lewy (CFL) condition (Matsumoto et al. 2019). In order to determine c p , we have to take care the e ff ect of non-uniform grid size thatemployed in the simulations. Following Zhang & Feng (2016), we set parameter α = ∆ h i , j , k c h / c p = . c p is chosen to satisfy thisrelation. Here ∆ h i , j , k is the minimum grid size that is taken locally: ∆ h i , j , k = min (cid:16) ∆ r i , r i ∆ θ j , r i sin θ j ∆ φ k (cid:17) , (B23)where ∆ r i = r i + / − r i − / , ∆ θ j = θ j + / − θ j − / , ∆ φ k = φ k + / − φ k − / . Practically we need only the damping rate (the coe ffi cient of ψ in RHSof equation A9), namely, c h / c p = α c h / ∆ h i , j , k , which indeed has a unit of 1 / s, and do not need to evaluate c p itself. APPENDIX C: RECONSTRUCTION OF VOLUME AND AREA AVERAGES
In order to achieve the high-order accuracy in space, the reconstruction of physical variables at the cell surface is necessary in the procedureof the estimation for the numerical fluxes. The piecewise linear method (PLM) of reconstruction (second-order) is implemented in our code.Note that since both finite volume and area methods are used in our code, the reconstructions of volume and area averages are necessary.Following the method of Skinner & Ostriker (2010) and Mignone (2014), the reconstruction of the volume average is obtained. First, weconsider 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, a i = (cid:82) r i + / r i − / a ( r ) r dr (cid:82) r i + / r i − / r dr . (C1)The physical variables at the left and right interfaces inside the zone i are then given by a L , i = a i − ∆ a i (1 + γ i ) , (C2) a R , i = a i + ∆ a i (1 − γ i ) , (C3)where γ i = ∆ r i r i + ∆ r i / r i (C4)is a correction factor for curvature in the spherical coordinates and ∆ r i = r i + / − r i − / . (C5)Here, ∆ a i is the di ff erence slope of the physical variable. In our code, the modified van Leer (VL) limiter proposed in Mignone (2014) isimplemented for the slope limiter to achieve the total variation diminishing (TVD); ∆ a i = ∆ a Fi ϕ VL ( υ ) , (C6)where υ = ∆ a Bi ∆ a Fi . (C7)The forward- and backward-di ff erence slopes are given by ∆ a Fi = ( a i + − a i ) ∆ r i (cid:104) r (cid:105) i + − (cid:104) r (cid:105) i , (C8) ∆ a Bi = ( a i − a i − ) ∆ r i (cid:104) r (cid:105) i − (cid:104) r (cid:105) i − , (C9)where (cid:104) r (cid:105) i is the cell-volume-averaged value inside zone i ; (cid:104) r (cid:105) i = (cid:82) r i + / r i − / r dr (cid:82) r i + / r i − / r dr = r i + r i ∆ r r i + ∆ r . (C10)The modified VL limiter is defined as follows (see Mignone 2014, for details); ϕ VL ( υ ) = υ ( c Fi υ + c Bi ) υ + ( c Fi + c Bi − υ + for υ ≥ , υ = , (C11) MNRAS000
In order to achieve the high-order accuracy in space, the reconstruction of physical variables at the cell surface is necessary in the procedureof the estimation for the numerical fluxes. The piecewise linear method (PLM) of reconstruction (second-order) is implemented in our code.Note that since both finite volume and area methods are used in our code, the reconstructions of volume and area averages are necessary.Following the method of Skinner & Ostriker (2010) and Mignone (2014), the reconstruction of the volume average is obtained. First, weconsider 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, a i = (cid:82) r i + / r i − / a ( r ) r dr (cid:82) r i + / r i − / r dr . (C1)The physical variables at the left and right interfaces inside the zone i are then given by a L , i = a i − ∆ a i (1 + γ i ) , (C2) a R , i = a i + ∆ a i (1 − γ i ) , (C3)where γ i = ∆ r i r i + ∆ r i / r i (C4)is a correction factor for curvature in the spherical coordinates and ∆ r i = r i + / − r i − / . (C5)Here, ∆ a i is the di ff erence slope of the physical variable. In our code, the modified van Leer (VL) limiter proposed in Mignone (2014) isimplemented for the slope limiter to achieve the total variation diminishing (TVD); ∆ a i = ∆ a Fi ϕ VL ( υ ) , (C6)where υ = ∆ a Bi ∆ a Fi . (C7)The forward- and backward-di ff erence slopes are given by ∆ a Fi = ( a i + − a i ) ∆ r i (cid:104) r (cid:105) i + − (cid:104) r (cid:105) i , (C8) ∆ a Bi = ( a i − a i − ) ∆ r i (cid:104) r (cid:105) i − (cid:104) r (cid:105) i − , (C9)where (cid:104) r (cid:105) i is the cell-volume-averaged value inside zone i ; (cid:104) r (cid:105) i = (cid:82) r i + / r i − / r dr (cid:82) r i + / r i − / r dr = r i + r i ∆ r r i + ∆ r . (C10)The modified VL limiter is defined as follows (see Mignone 2014, for details); ϕ VL ( υ ) = υ ( c Fi υ + c Bi ) υ + ( c Fi + c Bi − υ + for υ ≥ , υ = , (C11) MNRAS000 , 1–22 (2020) J. Matsumoto et al. where c Fi = (cid:104) r (cid:105) i + − (cid:104) r (cid:105) i r i + / − (cid:104) r (cid:105) i , (C12) c Bi = (cid:104) r (cid:105) i − (cid:104) r (cid:105) i − (cid:104) r (cid:105) i − r i − / . (C13)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 a j = (cid:82) θ j + / θ j − / a ( θ ) sin θ d θ (cid:82) θ j + / θ j − / sin θ d θ . (C14)The physical variables at the left and right interfaces inside the zone j are then given by a L , j = a j − ∆ a j (1 + γ j ) , (C15) a R , j = a j + ∆ a j (1 − γ j ) , (C16)where γ j = cos θ j sin θ j (cid:20) ∆ θ j − cos( ∆ θ j / ∆ θ j / (cid:21) (C17)is also the correction factor for curvature in the spherical coordinates and ∆ θ j = θ j + / − θ j − / . (C18)Here, ∆ a j = ∆ a Fj ϕ VL is the di ff erence slope of the physical variable in the θ -direction. The forward- and backward-di ff erence slopes used inthe VL limiter are given by ∆ a Fj = ( a j + − a j ) ∆ θ j (cid:104) θ (cid:105) j + − (cid:104) θ (cid:105) j , (C19) ∆ a Bj = ( a j − a j − ) ∆ θ j (cid:104) θ (cid:105) j − (cid:104) θ (cid:105) j − , (C20)where (cid:104) θ (cid:105) j is the cell-volume-averaged value inside zone j ; (cid:104) θ (cid:105) j = (cid:82) θ j + / θ j − / θ sin θ d θ (cid:82) θ j + / θ j − / sin θ d θ = θ j + cos θ j sin θ j (cid:20) − ∆ θ j ∆ θ j / ∆ θ j / (cid:21) . (C21)The coe ffi cients c Fi and c Bi in equation (C11) are replaced by c Fj = (cid:104) θ (cid:105) j + − (cid:104) θ (cid:105) j θ j + / − (cid:104) θ (cid:105) j (C22)and c Bj = (cid:104) θ (cid:105) j − (cid:104) θ (cid:105) j − (cid:104) θ (cid:105) j − θ j − / , (C23)respectively.In the φ -direction, the reconstruction of the physical variable is given by the same formula as the Cartesian coordinates. The cell-volume-averaged physical value inside zone z defined at the cell center is given by a k = (cid:82) φ k + / φ k − / a ( φ ) d φ (cid:82) φ k + / φ k − / d φ . (C24)The physical variables at the left and right interfaces inside the zone k are reconstructed as follows; a L , k = a k − ∆ a k , (C25) a R , k = a k + ∆ a k , (C26)where ∆ a k = ∆ a Fk ϕ VL is the di ff erence slope of the physical variable in the φ -direction. These are simply linear interpolations withoutcorrection factors for curvature. The forward- and backward-di ff erence slopes used in the VL limiter are defined by ∆ a Fk = a k + − a k , (C27) ∆ a Bk = a k − a k − . (C28) MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe Note that the coe ffi cients c Fk and c Bk of the VL limiter are 2 in this case.Next, we consider the reconstruction of area averaged values. The cell-area-averaged physical value inside zone i along the r -directiondefined at the cell center is given by a i = (cid:82) r i + / r i − / a ( r ) rdr (cid:82) r i + / r i − / rdr . (C29)This is the same formula for the case with the reconstruction of volume averaged value along the radial direction in the cylindrical coordinates(see Skinner & Ostriker 2010). The physical variables at the left and right interfaces inside the zone i are given by equations (C2) and (C3),respectively. However, the correction factor for curvature is γ i = ∆ r i r i (C30)in this case. In the reconstruction of area averaged values, we also employ the VL limiter (C6) and (C11). The cell-area-averaged value of r inside zone i used in the VL limiter is given by (cid:104) r (cid:105) i = r i + ( ∆ r i ) r i . (C31)In the θ -direction, there are two types of reconstructions. It depends on the direction of the area. When the area is perpendicular to the r -direction, the cell-area-averaged physical value inside zone j defined at the cell center is given by a j = (cid:82) θ j + / θ j − / a ( θ ) sin θ d θ (cid:82) θ j + / θ j − / sin θ d θ . (C32)This is the same formula as the cell-volume-averaged physical value inside zone j (C14). Therefore, the reconstructions of the physical valueat the cell surfaces are given by equations (C15) and (C16). On the other hand, when the area is perpendicular to the φ -direction, a j = (cid:82) θ j + / θ j − / a ( θ ) d θ (cid:82) θ j + / θ j − / d θ . (C33)This simply results in the linear interpolation at the cell surface like equation (C25) or (C26) although the spacial index k in the equations isreplaced by j .In the φ -direction, the reconstruction of the physical variable is the same as the case with the volume-averaged reconstruction. It issimply 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 previousworks with similar setups (Londrillo & Del Zanna 2000; Stone et al. 2008; Stone & Gardiner 2009; Skinner & Ostriker 2010). The initialcondition that we adopt is as follows. We set a static background of ρ = . , p = .
1. The background is magnetized as B x = . B is already normalized by √ π ). In the background, we put a hot gas of p =
10 in the region of r < .
1. The polytropic indexof EOS, Γ , is set to . The domain of [0 , . × [0 , π ] × [0 , π ] of spherical polar coordinate is uniformly covered by 200 × ×
400 grids,respectively.The time snapshot of t = . ff ect of the magnetic field. It is alsoimportant to check the density distribution. The top right, bottom left and bottom right panels are the 2D slice of y = x = z = x -direction is typically seen in thistest. The evolution of the shock does not depend on the coordinate system. Though the structure of the mesh is quite di ff erent 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 acoordinate independent feature would support the correctness of our code implementation. MNRAS000
400 grids,respectively.The time snapshot of t = . ff ect of the magnetic field. It is alsoimportant to check the density distribution. The top right, bottom left and bottom right panels are the 2D slice of y = x = z = x -direction is typically seen in thistest. The evolution of the shock does not depend on the coordinate system. Though the structure of the mesh is quite di ff erent 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 acoordinate independent feature would support the correctness of our code implementation. MNRAS000 , 1–22 (2020) J. Matsumoto et al.
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 showthe density structure in 2D slice of y = x = z =
0, respectively. The uniform magnetic field is initially imposed in x -direction. 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.0B12models with the resolution of ∆ θ = π/ 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 inSection 3, this tendency has been also observed in the fiducial resolution runs. Our resolution study demonstrates that the delay of the shockrevival occurs regardless of the angular resolution in the θ -direction although more careful studies covering the wide range of resolution arenecessary in order to draw a robust conclusion. MNRAS , 1–22 (2020) agnetic field dependence of neutrino-driven SNe r s h [ k m ] Time after bounce [ms]s27.0B10s27.0B12s27.0B10 highs27.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 bluelines represent the weak ( B = G) and strong ( B = G) magnetic field model, respectively. For the reference, the results of fiducial runs aredemonstrated by doted lines.MNRAS000