Neutrino-driven massive stellar explosions in 3D fostered by magnetic fields via turbulent $\alpha$-effect

We investigate the influence of magnetic field amplification on the core-collapse supernovae in highly magnetized progenitors through three-dimensional simulations. By considering rotating models, we observe a strong correlation between the exponential growth of the magnetic field in the gain region and the initiation of shock revival, with a faster onset compared to the non-rotating model. We highlight that the mean magnetic field experiences exponential amplification as a result of $\alpha$-effect in the dynamo process, which works efficiently with the increasing kinetic helicity of the turbulence within the gain region. Our findings indicate that the significant amplification of the mean magnetic fields leads to the development of locally intense turbulent magnetic fields, particularly in the vicinity of the poles, thereby promoting the revival of the shock by neutrino heating.


INTRODUCTION
The impact of the magnetic field on the explosion mechanisms of core-collapse supernovae (CCSNe) is a matter of significant importance.The combination of the rapid rotation and magnetic field of a massive star leads to the launching of magnetically-driven jets during the collapse of the massive star (e.g., Takiwaki et al. 2009;Sawai & Yamada 2016;Kuroda et al. 2020;Obergaulinger & Aloy 2020, 2021;Bugli et al. 2021Bugli et al. , 2023;;Powell et al. 2023).Even in the context of a slowly-rotating progenitor, the amplification of the magnetic field due to the turbulence is expected to foster the onset of neutrino-driven explosions (Obergaulinger et al. 2014;Müller & Varma 2020;Matsumoto et al. 2022;Varma et al. 2023).From the seminal work by Müller & Hillebrandt (1979), one of the grand challenges in the CCSN theory has long been to unveil the role of the magnetic fields in the explosion mechanism (e.g., Burrows & Vartanyan 2021;Janka 2012;Kotake et al. 2006 for review).
In the dynamo theory, a mean kinetic helicity of the turbulence is one possible origin to amplify the mean magnetic field drastically (Brandenburg & Subramanian 2005).This is the so-called α-effect.The mean kinetic helicity is naturally generated in a rotating convection system by a swirling motion of the fluid through the interaction between the Coriolis force and the convection flow (Spruit et al. 1990;Miesch 2005;Miesch & Toomre 2009).In the context of CCSNe, the proto-neutron star (PNS) and neutrino-driven convection develop while the massive star collapses (e.g., Janka 2012; Radice et al. 2018).Although the generation mechanism of ⋆ Email:jin@rk.phys.keio.ac.jp, jin@kusastro.kyoto-u.ac.jp the strong magnetic field of PNSs is still under debate (Akiyama et al. 2003;Guilet et al. 2015;Reboul-Salze et al. 2020, 2022;Barrère et al. 2022), the convection dynamo is considered as a key process for the origin of the magnetic field of PNSs and has been extensively studied in the PNS convection (Thompson & Duncan 1993;Bonanno et al. 2003;Rheinhardt & Geppert 2005;Raynaud et al. 2020;Masada et al. 2022;White et al. 2022).In fact, recent 3D dynamo simulations report that the turbulent α-effect, the origin of which is the kinetic helicity in the turbulent electromotive force based on the mean-field theory works for the amplification mechanism of the magnetic field in the PNS convection (Raynaud et al. 2020;Masada et al. 2022).
While the turbulent α-effect has been investigated in the PNS convection, its effect based on the generation of the kinetic helicity in the neutrino-driven convection on the explosion mechanism of the rotating massive star is not fully understood yet (see also Endeve et al. 2012 for the effect of the kinetic helicity in non-rotating progenitor).In this work, we focus on the impact of the magnetic field amplification due to the turbulent α-effect in the gain region on the explosion mechanism of CCSN in the context of the slowly rotating model through 3D magnetohydrodynamic (MHD) simulations.
This Letter is organized as follows.In Section 2, we describe our model in our calculations.We present the main results in Section 3. Section 4 is devoted to conclusions and discussions.

NUMERICAL MODELS
All numerical settings in this study are the same as our previous work (Matsumoto et al. 2022) except the initial angular velocity of the progenitor, ω ini .The employed pre-supernova progenitor model is 27 M ⊙ of Woosley et al. (2002), s27.We consider simple rotation profiles of s27 in a parametric manner because of the limited understanding of the spatial distribution of the rotation in the precollapse phase of massive stars.Following Takiwaki et al. (2004) and Takiwaki et al. (2009), we assume a cylindrical rotation as follows: where x and z indicate distance from the rotational axis and the equatorial plane and x 0 = z 0 = 1000 km.We set ω 0 = 0.3, 0.1 or 0 rad s −1 .The last one corresponds to a non-rotating model calculated and labelled as s27.0B12PPM5 in the previous work.The initial strength of the magnetic field in the z-direction of a central core of s27 is roughly 10 12 G.Matsumoto et al. (2022) provides the detailed structure of the field.The equation of state of Lattimer & Swesty (1991) with a nuclear incomprehensibility of K = 220 MeV is employed.Three models that are identified by the value of ω 0 are calculated using our MHD supernova code (3DnSNe; Takiwaki et al. 2016;Matsumoto et al. 2020Matsumoto et al. , 2022) ) with 5th-order spatial accuracy in a spherical coordinate system (r, θ, ϕ).The calculation domain covers a sphere whose radius is 5000 km with a resolution of n r × n θ × n ϕ = 480 × 64 × 128.The grid-cell size to the radial direction logarithmically stretches.The resolution of the polar angle is given by ∆(cosθ) = const.covering 0 ≤ θ ≤ π.The azimuthal angle is uniformly divided into ∆ϕ = π/64 covering 0 ≤ ϕ ≤ 2π.

RESULTS
Fig. 1 shows the temporal evolution of averaged shock radii (panel a) and the diagnostic explosion energy (panel b).Red, blue and green lines correspond to the models ω 0 = 0.3, 0.1 and 0 rad s −1 , respectively in each panel.The evolution of the mass accretion rate for the non-rotating model evaluated at r = 500 km from the center of the calculation domain, Ṁ, is plotted as a reference in panel (a).Note that it is not so different between all models before the shock revival because the centrifugal force in rotating models does not prevent the gravitational collapse of the s27 progenitor in the range of ω 0 we assume in this work.Therefore, the core bounce time (∼ 200 ms after the start of the simulation) is the same in all models.
As discussed in Matsumoto et al. (2022), in the non-rotating model (green line), the shock revival occurs at the sudden drop of Ṁ.The increased gas pressure behind the shock due to the neutrino heating overcomes the ram pressure of the accreting fluid onto the shock at this time and starts to push the shock outwardly because the ram pressure decreases thanks to the reduction of the density of the accreting fluid.The accumulation of the magnetic field just behind the shock contributes to the increase of the total (gas and magnetic) pressure and secondary supports the shock revival.
The shock evolution in rotating models (red and blue lines) is obviously fast compared to that in the non-rotating model (green line).This indicates that other mechanisms would be responsible for the faster explosion of the massive star through the interaction between the rotation and the magnetic field or the efficiency of the neutrino heating is enhanced.As discussed later, the locally amplified magnetic pressure due to the turbulent α-effect whose origin is the mean kinetic helicity in rotating models contributes to the faster explosion in this work.
On the one hand, the magnetic energy is locally larger than the thermal energy near the poles in the gain region to drive the fast explosion, but on the other hand, the total thermal energy in the gain region is large compared to the magnetic energy.Therefore, the thermal energy is dominant in the diagnostic explosion energy defined in Matsumoto et al. (2022) that includes the magnetic energy.In Fig. 1(b), the temporal evolution of the diagnostic explosion energy for all models is plotted.The diagnostic explosion energy in the faster explosion model is larger.It is reasonable because as the mass accretion rate is larger at the earlier phase of the evolution (see the black line in Fig. 1a), the neutrino luminosity also becomes larger in the faster explosion model at the timing of the onset of the shock revival.This gives plenty of thermal energy to the expanding ejecta through neutrino heating.
When we consider the overburden energy of unshocked region in the whole s27 progenitor beyond the calculation domain (Bruenn et al. 2013;Bollig et al. 2021;Mori et al. 2023), it is ∼ 10 51 erg at the final time of calculation (t pb = 400 ms) in all models.Therefore, the total energy (overburden and diagnostic explosion energy) of our models is negative.However, as pointed out in Marek & Janka (2009), the energy release by nuclear burning in shock-heated material is expected to compensate for the overburden energy.
Fig. 2 shows the snapshots of the shock surface and the magnetic field lines for model ω 0 = 0.3 rad s −1 at t pb = 100 ms (panel a) and t pb = 150 ms (panel b).Note that t pb means post-bounce time in the calculation for CCSN.The spatial scale size (illustrated by a white two-headed arrow in panel a) and the viewing position are the same in both panels.The shock position is represented by a whitish and transparent sphere in each panel.The color of the magnetic field lines represents the absolute strength of the magnetic field.Since we initially assume the strong magnetic field (B 0 = 10 12 G) within r ≤ 1000 km in this work, the magnetic field near the central region of the calculation domain is simply amplified to the B = 10 15 G level after the gravitational collapse of the progenitor through the magnetic flux conservation as follows: (2) The neutrino-driven convection starts to develop due to the negative gradient of the entropy around at t pb = 100 ms.The evolution of the magnetic field is different before/after this phase.The gravitational collapse of the matter in the radial direction leads to the split-monopole-like configuration of the magnetic field lines except the central part of the star before the onset of the neutrinodriven convection.Although the central core of the progenitor initially rotates rigidly, the accretion of the fluid results in a strong differential rotation near the rotational axis.It contributes to the magnetic field winding around the pole in the spherical coordinate system and the generation of the toroidal component of the magnetic field.These structures of the magnetic field are confirmed in Fig. 2(a).After the onset of the neutrino-driven convection, it induces the non-radial motion of the fluid and the turbulence in the gain region.They interact with the magnetic field lines through the swirling motion of the fluid.In Fig. 2(a), slightly bent magnetic field lines are observed behind the shock (gain region).Fig. 2(b) shows the magnetic field configuration after the shock revival.Although the azimuthal component of the magnetic field is dominant near the rotational axis even in this phase, each component of the magnetic field in the gain region (behind the shock) is comparable due to the interaction of the magnetic field with the turbulence.
Fig. 3 shows the 2D distribution of (a) plasma β = P gas /P mag , longitudinal averages of (b) plasma β, (c) kinetic helicity and (d) dynamo number around the gain region for model ω 0 = 0.3 rad s −1 at t pb = 139 ms (after the shock revival) where P gas and P mag are gas and magnetic pressure, respectively.Here, the longitudinal average of variable X is given by In Fig. 3(a), we can find the low β region (red) near the pole where the magnetic field is fully amplified to O(10 14 ) G as shown in Fig. 2(b).On the other hand, the strength of the magnetic field at the equatorial region is O(10 13 ) G. In the low β region, the magnetic pressure is larger than the gas pressure.This indicates that the magnetic pressure is the main driver for the explosion mechanism in our rotating model.However, as shown in Fig. 3(b), the longitudinal average of plasma β in the gain region is larger than unity in the meridional plane (r, θ).Therefore, the volume-averaged magnetic pressure in the gain region is less than the volume-averaged gas pressure.Locally amplified magnetic pressure is expected to contribute to the outward shock expansion.
As mentioned in the introduction, the rotating convection system is expected to generate kinetic helicity and result in the amplification of the magnetic field by the dynamo action called α-effect.To identify the amplification mechanism of the magnetic field in our rotating model, we follow the analysis of the mean-field theory of the magnetic field (Brandenburg & Subramanian 2005).The velocity and magnetic field are decomposed by the mean and tur- bulent components as follows: B(r, θ, ϕ) = ⟨B⟩(r, θ) + B ′ (r, θ, ϕ) . (5) Here, the prime represents the turbulent component.Note that we consider only the longitudinal average of the variable for the mean component and the time average is neglected for simplicity in our analysis.The mean kinetic helicity is defined by where ω ′ stands for the turbulent vorticity vector given by Under the first-order smoothing approximation, the induction equations of the mean and turbulent magnetic field are given by respectively.Here ϵ ≡ α⟨B⟩ − η t ∇ × ⟨B⟩ is a turbulent electromotive force and where τ cor is the correlation time.When the turbulent process is dominant in the amplification of the mean magnetic field, the first term of the RHS in the equation ( 8) is neglected.In addition, as we set η = 0 in our simulations, the induction equation for the mean magnetic field (8) is reduced as follows: By inserting the perturbation of the form δ⟨B⟩ ∝ exp(ik • x + σt) into the equation ( 12), we obtain the following dispersion relation: where σ and k = |k| are a growth rate of the mean magnetic field and a wave number, respectively.The condition for the exponential growth of the mean magnetic field due to the α-effect (σ > 0) is obtained as follows: Here D is a dynamo number.From Fig. 3, we can find a good correlation between the low β region (red) and the distribution of the high absolute value of the kinetic helicity.This indicates that the kinetic helicity amplifies the magnetic field because the strength of the magnetic field in the low β region is relatively strong compared to that in the high β region (green and blue).In addition, we can also observe a good correlation between the 2D distribution of the kinetic helicity and dynamo number.In Fig. 3(d), the coloured region represents a high dynamo number (D > 1) only in the gain region.Here we set 1/k = 100 km that roughly corresponds to the length scale of the gain region, L gain .This implies that the α-effect is responsible for the amplification of the magnetic field in our rotating model.
To verify the origin of the magnetic field amplification in our rotating model quantitatively, the time evolution of the mean and turbulent kinetic and magnetic energy in the gain region around the onset of the shock revival for model ω 0 = 0.3 rad s −1 is shown in Fig. 4. The mean and turbulent magnetic energy is exponentially amplified during the period between t = 105 and 110 ms.It is closely linked to the onset of the shock revival (see red line in Fig. 1a).Assuming a large dynamo number, the growth rate of the magnetic amplification due to the α-effect is roughly estimated from the dispersion relation (13) as follows: where τ cor ∼ τ adv ≡ L gain /v adv .Here τ adv and v adv ∼ 10 9 cm s −1 are advection timescale and advection velocity in the gain region, respectively.In our rotating model, the kinetic helicity suddenly increases to the order of 10 12 cm s −2 after around t pb = 100 ms (see Fig. 3c).The anticipated growth rate of the mean magnetic energy from the mean-field theory, 2σ, is also plotted by a black line in Fig. 4. In addition to the growth rate of the mean magnetic energy, the growth rate of the turbulent magnetic energy is roughly equal to 2σ.Since, in the first-order smoothing approximation of the mean-field theory, the time evolution of the turbulent magnetic field is given by the induction equation ( 9), it is reasonable that the growth rate of it is almost same as the growth rate of the mean magnetic field.In this work, we consider the slowly rotating model of the progenitor and the rotational energy of the fluid is smaller  than the turbulent kinetic energy in the gain region during the exponential amplification phase of the magnetic field.This implies that the convective dynamo effectively works in the magnetic amplification compared to magnetorotational instability (Balbus & Hawley 1998).It is clear from Fig. 4 that the turbulent magnetic energy (red line) in the gain region is large compared to the mean magnetic energy (blue line).Therefore, the turbulent magnetic pressure that is amplified via α-dynamo action of the mean magnetic field is responsible for the fast explosion in our rotating model.

SUMMARY AND DISCUSSION
The impact of magnetic field amplification on the explosion mechanism of a strongly magnetized 27 M ⊙ pre-supernova progenitor is studied by performing 3D MHD supernova simulations in the context of the slowly-rotating progenitor models.We find that the exponential growth of the magnetic field in the gain region in our rotating models has a good correlation with the onset of the shock revival, which is faster than that of the non-rotating model.In addition, the explosion energy in the more rapidly rotating model is larger.We point out that the mean magnetic field is exponentially amplified due to the α-effect in the dynamo process, which works efficiently with the growing kinetic helicity of the turbulence in the rotating gain region.Our results show that the drastic amplification of the mean magnetic fields leads to the generation of localized and strong turbulent magnetic fields near poles in the gain region, promoting the shock revival by neutrino heating.
Obviously, we need more sophisticated numerical methods to obtain more accurate predictions to model shock revival through neutrino-heating in multi-dimension, as has been established by various studies (e.g., Buras et al. 2006;Suwa et al. 2010;Takiwaki et al. 2012;Nakamura et al. 2015;Pan et al. 2016;Summa et al. 2016;O'Connor & Couch 2018a;Burrows et al. 2020;Nagakura et al. 2020;Bollig et al. 2021;Vartanyan et al. 2022).Despite the existence of some general relativistic simulations (Müller et al. 2012;O'Connor & Couch 2018b;Kuroda et al. 2022;Rahman et al. 2022), the majority of the studies use phenomenological GR potential.The ab-initio neutrino transport requires huge computational resources (Iwakami et al. 2020) and approximate methods are usually employed.Furthermore, the effect of neutrino oscillation should be taken into account (Nagakura 2023;Ehring et al. 2023) and the impact of neutrino reactions remains to be examined (Kotake et al. 2018;Sugiura et al. 2022).Moreover, properties of high density equation of state are still under debate (Fischer et al. 2014(Fischer et al. , 2020)).
In this Letter, we have proposed the so-far-unidentified mechanism to forge the explosion onset via the turbulent α-effect for a strongly magnetized and slowly-rotating progenitor.The diversity of supernova explosions (Kulkarni 2012;Taddia et al. 2018;Martinez et al. 2022) might originate from the property of the progenitors (e.g., Smartt 2009;Aguilera-Dena et al. 2020;Takahashi & Langer 2021).Especially, the binary evolution of stars is important because it may give complicated features and a wide range of spin parameters for the progenitor, which are not established in a single star evolution through mass transfer from a companion, tidal locking between binaries or stellar mergers (Cantiello et al. 2007;Chatzopoulos et al. 2020).Recent advancements in 3D MHD simulations of the stellar evolution enable us to unveil the dynamics and evolution of the system where the rotation of the progenitor interacts with the magnetic field (Varma & Müller 2021;Yoshida et al. 2021;McNeill & Müller 2022;Fields 2022).In order to comprehend the origin of the diversity of CCSNe, sophisticated and a wide range of parameter studies for the magnetic field and spin period of the progenitor are needed, though computationally expensive, towards which we, including CCSN modellers in the globe are making the steady step.

Figure 1 .
Figure 1.Panel (a): Evolution of the averaged shock radii, r sh,ave , for all models and the mass accretion rate at r = 500 km, Ṁ, (black line) for the non-rotating model.Panel (b): Comparison of the diagnostic explosion energy including magnetic energy defined in Matsumoto et al. (2022).

Figure 2 .
Figure 2. Snapshots of the shock surface (whitish sphere) and magnetic field lines for model ω 0 = 0.3 rad s −1 .Panels (a) and (b) correspond to t pb = 100 and 150 ms, respectively.The spatial scale is represented by a white two-headed arrow in panel (a) that is parallel to the z-axis.The spatial scale and the viewing angle of the central object are fixed for both panels.

Figure 3 .
Figure 3. 2D distribution of (a) plasma β, longitudinal averages of (b) plasma β, (c) kinetic helicity and (d) dynamo number around the gain region for model ω 0 = 0.3 at t pb = 139 ms.Note that x ′ in panel (a) means one axis, whose direction is along ϕ = 0.025 on x-y plane and origin matches that of x.Data in the negative region of x are copied from data in the positive region in panels (b), (c) and (d).

Figure 4 .
Figure 4. Temporal evolution of the kinetic and magnetic turbulent and mean energy in the gain region around the onset of the shock revival for model ω 0 = 0.3 rad s −1 .