© A. Reboul-Salze et al. 2021

1. Introduction

Magnetars are a class of young and highly magnetized neutron stars that produce a wide variety of outstanding emission at X-ray or soft γ-ray energies (Kouveliotou et al. 1998; Kaspi & Beloborodov 2017, and references therein). The period and the spin-down of these objects are measured by the long time follow-up of their pulsed X-ray activity, and a surface dipolar magnetic field of $B_{\mathrm{dip}}={10}^{14}{\left(\frac{P}{5\mathrm{s}}\right)}^{\frac{1}{2}}{\left(\frac{\dot{P}}{{10}^{-11}\mathrm{s}{\mathrm{s}}^{-1}}\right)}^{\frac{1}{2}}\mathrm{G}\simeq {10}^{14}-{10}^{15}$ G can be inferred under the assumption of magnetic dipole spin-down (Olausen & Kaspi 2014)¹. Their activity can be explained by the decay of their ultrastrong magnetic field and also includes short bursts (Götz et al. 2006), large outbursts (Coti Zelati et al. 2018), giant flares (Hurley et al. 2005), and quasi-periodic oscillations (Israel et al. 2005). Absorption lines are also been detected in outbursts for two objects and interpreted as proton cyclotron lines (Tiengo et al. 2013; Rodríguez Castillo et al. 2016). This suggests the presence of a strong non-dipolar surface field, whereas a weaker dipolar component is derived from the timing parameters of these objects.

Magnetic fields, especially in the presence of fast rotation, could play an important role in the dynamics of core-collapse supernovae and have garnered considerable interest in the last decade. Indeed, a strong magnetic field would efficiently extract the large rotational energy of a protoneutron star (PNS) rotating with a period of a few milliseconds. The presence of a magnetic field can impact the explosion by converting the energy of differential rotation into thermal energy (Thompson et al. 2005) and/or into a large-scale magnetic field, which can launch jets and lead to a magnetorotational explosion (Moiseenko et al. 2006; Shibata et al. 2006; Dessart et al. 2008; Winteler et al. 2012; Mösta et al. 2014; Obergaulinger et al. 2018; Bugli et al. 2020; Kuroda et al. 2020). Moreover, a millisecond proto-magnetar is a potential central engine for long gamma-ray bursts (Duncan & Thompson 1992; Metzger et al. 2011, 2018) associated with “hypernovae” or supernovae type Ic Broad Lined (Drout et al. 2011). These rare events are characterized by a kinetic energy ten times higher than that of standard supernovae. Furthermore, millisecond magnetars have been invoked to explain some superluminous supernovae through a delayed energy injection due to the dipole spin-down luminosity (Nicholl et al. 2013; Inserra et al. 2013; Margalit et al. 2018). Millisecond magnetars may also be formed in binary neutron star mergers. This can provide a natural explanation for the extended emission and the plateau phase in X-ray sources associated with a fraction of short gamma-ray bursts (Metzger et al. 2008; Bucciantini et al. 2012; Rowlinson et al. 2013; Gompertz et al. 2014). An X-ray transient has also recently been detected and interpreted as the formation of a magnetar in the aftermath of a binary neutron-star merger (Xue et al. 2019).

Several scenarios have been invoked to explain the origin of the magnetic field in magnetars. It may stem from the field of the progenitor amplified by magnetic flux conservation. This scenario can lead to the strongest magnetic fields in the case of highly magnetized progenitors that could be formed in stellar mergers (Schneider et al. 2019), although it may not explain the formation of millisecond magnetars since highly magnetized progenitors are slow rotators (Spruit 2008; Shultz et al. 2018). An alternative process is an in situ magnetic field amplification by a turbulent dynamo in the PNS, either by a convective dynamo (Thompson & Duncan 1993; Raynaud et al. 2020) or the magnetorotational instability (MRI, see Akiyama et al. 2003; Obergaulinger et al. 2009).

The first local analytical study of the MRI in the context of Keplerian accretion disks by Balbus & Hawley (1991) showed that in presence of differential rotation small seed perturbations are amplified exponentially with time. The first local simulations by Hawley & Balbus (1992) showed that the turbulent velocity and magnetic field reach a statistically stationary state. In ideal magnetohydrodynamics (MHD) and with an initial vertical uniform magnetic field, the growth rate of the instability is on the order of the rotation rate Ω. In this case, the wavelength of the fastest growing mode is proportional to the magnetic field intensity. Therefore, the weaker the magnetic field, the shorter the MRI wavelength and the more difficult it becomes to resolve in a global model. This instability has thus widely been studied in the local approximation, either analytically or by using “shearing box” simulations representing a part of the accretion disk. The linear growth of the MRI has been studied with thermal stratification (entropy and composition gradient) and with diffusion processes (viscosity and resistivity; Balbus 1995; Menou et al. 2004; Pessah & Chan 2008).

Core-collapse supernova simulations show a strong differential rotation in the PNS (Akiyama et al. 2003; Ott et al. 2006). Numerical models in the context of supernovae have shown that the MRI can grow on shorter timescales than the successful explosion time and that an efficient amplification of the magnetic field occurs at small scales (Obergaulinger et al. 2009; Guilet & Müller 2015; Rembiasz et al. 2016). The influence of the specific physical conditions of PNSs were studied in these local models. First, the pressure gradient rather than the centrifugal force balances gravity, which can lead to a non-Keplerian rotation profile. In the case of a steeper rotation profile, a stronger MRI turbulence develops in shearing box simulations (Masada et al. 2012). Secondly, due to the high density inside the PNS, neutrinos are in the diffusive regime and their transport of momentum can be described with a high viscosity, which can limit the growth of the MRI if the initial magnetic field is too low (Guilet et al. 2015). Finally, the buoyancy forces driven by the entropy and lepton fraction can reduce the MRI turbulence in the case of stable stratification (Guilet & Müller 2015).

The impact of the spherical geometry of the full PNS on the MRI turbulence and the ability of the MRI to generate a large-scale field, similar to the inferred magnetic field of magnetars, is still unknown. The first attempts to address this question rely on semi-global models that include radial gradients of density and entropy (Obergaulinger et al. 2009; Masada et al. 2015). However, these models remain local at least vertically and therefore cannot investigate the generation of a large-scale magnetic field. Global axisymmetric simulations of the MRI also show magnetic field amplification (Sawai et al. 2013; Sawai & Yamada 2016). Mösta et al. (2015) performed the first simulations describing a quarter of the PNS with a high enough resolution to resolve the MRI wavelength and showed the development of the MRI turbulence. The model was, however, started with an initial magnetic dipole and therefore did not demonstrate the generation of a magnetar-like magnetic dipole.

This paper studies for the first time the global properties of the MRI in a full 3D spherical model, where no initial large-scale magnetic field is assumed. To resolve the MRI wavelength with a reasonable resolution, we use a sufficiently strong initial magnetic field. However, only the small scales are initialized in order to study the generation of the magnetic dipole. This intense and small-scale magnetic field can be interpreted as the result of the first amplification described in local models. With respect to previous studies, a different approach is also used for the physical setup, which is reduced to its most fundamental ingredients. This has the advantage of providing a useful reference for our physical understanding, while at the same time drastically reducing the computational cost and enabling long simulation times and the exploration of the parameter space.

The paper is organized as follows. In Sect. 2, we describe the physical and numerical setup. The results are then presented in Sect. 3 for the saturated nonlinear phase of the MRI, and in Sect. 4 for the comparison of our global model to local models of the MRI. Finally, we discuss the validity of our assumptions in Sect. 5 and draw our conclusions in Sect. 6.

2. Numerical setup

2.1. Governing equations

The simulations performed in this article are designed to represent a fast rotating PNS. We assume that the hot PNS has a mass of 1.3 M_⊙ and a radius of r_o = 25 km. As the PNS is in solid body rotation for a radius r ≤ 10 km (see Sect. 2.3 for the assumed rotation profile), the inner core of the PNS is stable to the MRI and it can be excluded from the simulation domain. We chose an inner core with a radius r_i = 6.25 km. The shell gap D ≡ r_o − r_i = 18.75 km is the characteristic length of the simulation domain. In our model, we assume that neutrinos are in the diffusive regime such that their effects on the dynamics can be appropriately described by a viscosity ν (Guilet et al. 2015). The incompressible approximation is used for the sake of simplicity. We assume a uniform density ρ₀ = 4 × 10¹³ g cm⁻³ and neglect buoyancy effects. The diffusive incompressible MHD equations describing the dynamics of the PNS in a rotating frame at an angular frequency Ω₀ = 1000 s⁻¹ read

$$\begin{array}{cc}& {\displaystyle \frac{\partial \mathit{u}}{\partial t}}+(\mathit{u}·\mathbf{\nabla })\mathit{u}=-\mathbf{\nabla }p\prime -2{\mathbf{\Omega }}_0\times \mathit{u}+{\displaystyle \frac{1}{{\mu }_0{\rho }_0}}(\mathbf{\nabla }\times \mathit{B})\times \mathit{B}+\nu \mathrm{\Delta }\mathit{u},\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\displaystyle \frac{\partial \mathit{B}}{\partial t}}=\mathbf{\nabla }\times (\mathit{u}\times \mathit{B}-\eta \mathbf{\nabla }\times \mathit{B}),\hfill \end{array}$$(2)

$$\begin{array}{cc}& \mathbf{\nabla }·\mathit{u}=0,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \mathbf{\nabla }·\mathit{B}=0,\hfill \end{array}$$(4)

where u is the flow velocity, B the magnetic field, η the magnetic diffusivity, p′ the reduced pressure (i.e., the pressure normalized by the density) and μ₀ the vacuum permeability.

2.2. Numerical methods

In order to solve the incompressible MHD system of Eqs. (1)–(4), we used the pseudo-spectral code MagIC (Wicht 2002; Gastine & Wicht 2012; Schaeffer 2013)². MagIC solves the 3D MHD equations in a spherical shell using a poloidal-toroidal decomposition for the velocity and the magnetic field,

$$\begin{array}{cc}& \mathit{u}=\mathbf{\nabla }\times \mathbf{\nabla }\times \left(W{\mathit{e}}_{\mathrm{r}}\right)+\mathbf{\nabla }\times \left(Z{\mathit{e}}_{\mathrm{r}}\right),\hfill \end{array}$$(5)

$$\begin{array}{cc}& \mathit{B}=\mathbf{\nabla }\times \mathbf{\nabla }\times \left(b{\mathit{e}}_{\mathrm{r}}\right)+\mathbf{\nabla }\times \left(a_j{\mathit{e}}_{\mathrm{r}}\right),\hfill \end{array}$$(6)

where W and Z are respectively the poloidal and toroidal kinetic scalar potentials, while b and a_j are the magnetic ones. The scalar potentials and the reduced pressure p′ are decomposed on spherical harmonics for the colatitude θ and the longitude ϕ angles, together with Chebyshev polynomials in the radial direction. The time stepping scheme is a mixed implicit-explicit scheme: A Crank-Nicolson scheme is used to advance the linear terms and a second order Adams-Bashforth scheme to advance the nonlinear terms and the Coriolis force. The linear terms are computed in the spectral space, while the nonlinear terms and the Coriolis force are computed in the physical space and transformed back to the spectral space. For more detailed descriptions of the numerical method and the associated spectral transforms, the reader is referred to Gilman & Glatzmaier (1981), Tilgner & Busse (1997), and Christensen & Wicht (2015).

All the simulations presented in this paper were performed using a standard grid resolution of (n_r,  n_θ,  n_ϕ) = (257,  512,  1024). The resolution was chosen to ensure that the dissipation scales are resolved. Indeed, the maxima of viscous and resistive dissipation are respectively at the spherical harmonic orders l_ν ≃ 40 and l_η ≃ 85, which means that around ten cells resolve the resistive maximum.

2.3. Initial conditions

Many core-collapse simulations with a fast rotating progenitor have shown that the PNS is differentially rotating for several hundreds of milliseconds (e.g, Akiyama et al. 2003; Ott et al. 2006; Obergaulinger et al. 2018; Bugli et al. 2020). The contraction of the PNS and accretion on the PNS accelerates it while the angular momentum is transported outside the PNS by the MRI. However, it would not be possible to sustain a differentially rotating profile within an isolated shell of the PNS and in presence of strong magnetic fields, unless some forcing mechanism were to be applied at the boundaries. Therefore, in our simplified setup, we chose to artificially force the differential rotation by keeping the initial differential rotation profile at the outer boundary constant during the duration of the simulation. This leads to a quasi-stationary state that is simpler to study and to compare among different models. The internal rotation profile of this quasi-stationary state can evolve and is determined by the combined effects of the backreaction of the magnetic field and the rotational forcing.

As initial and outer boundary condition, we used a cylindrical rotation profile inspired by core-collapse simulations with a central part that rotates like a solid body and an outer part that is differentially rotating with a power law dependency:

$$\begin{array}{c}\hfill \mathrm{\Omega }(s)=\frac{{\mathrm{\Omega }}_{\mathrm{i}}}{{(1+{\left(\frac{s}{0.4r_{\mathrm{o}}}\right)}^{20q_0})}^{0.05}},\end{array}$$(7)

where s is the cylindrical radius, q₀ corresponds to the shear rate $q\equiv -\frac{r}{\mathrm{\Omega }}\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}r}$ in the outer part, and Ω_i is the rotation rate of the inner core. This rotation rate Ω_i was computed so that the ratio of total angular momentum over moment of inertia is equal to the frame rotation rate Ω₀ defined in Sect. 2.1. The profile shown in Fig. 1 with q₀ = 1.25 has a smooth transition between solid body rotation when s <  0.4 r_o and power law differential rotation when s >  0.4 r_o.

Fig. 1. Rotation profile inside the PNS as a function of the cylindrical radius. The inner region (s <  10 km) is in solid body rotation and the outer region (s >  10 km) follows a power law Ω ∝ r−q0 with q0 = 1.25.

Our model assumes a weak progenitor magnetic field and is designed to study in situ magnetic field amplification by the MRI, especially the large-scale field generation. Several local studies have already shown an efficient amplification of the magnetic field on small scales (e.g, Obergaulinger et al. 2009; Masada et al. 2015; Guilet & Müller 2015; Rembiasz et al. 2016). Therefore, the magnetic field may be initialized by small-scale modes with magnetic field strength high enough to resolve the MRI fastest growing mode. The initial poloidal magnetic potential is a random superposition of modes with spherical harmonics indices (l, m) and a radial profile of the form

$$\begin{array}{c}\hfill b_{l,m}\propto \{\begin{array}{cc}0\hfill & \mathrm{for}r<7.5\mathrm{km}\hfill \\ cos(k_{\mathrm{r}}(r-r_{\mathrm{i}}))\hfill & \mathrm{for}r>12.5\mathrm{km}\hfill \end{array},\end{array}$$(8)

where k_r is the radial wavenumber and has a smooth transition between 7.5 km and 12.5 km. We selected Fourier radial modes and spherical harmonic modes, whose wavelengths λ_r = 2π/k_r and ${\lambda }_{\mathrm{l}}=\sqrt{\frac{r^2}{l(l+1)}}$ lie within the range [L_min,  L_max]. We note that this initial magnetic field has a vanishing net magnetic flux over the simulation domain. The amplitude of these modes is initialized randomly with a wavelength dependence that provides a flat energy spectrum. Figure 2 shows an example of the initial poloidal magnetic field on a ϕ-slice with a typical value of [L_min,  L_max]=[0.23 r_o,  0.38 r_o].

Fig. 2. Slice of the Bθ component of the initial magnetic field for the meridional plane ϕ = 0 with [Lmin,  Lmax]=[0.23 ro,  0.38 ro] and B0 = 8.8 × 1014 G.

The initial root mean square magnetic field strength B₀ is varied from B₀ = 6.31 × 10¹⁴ G to B₀ = 3.36 × 10¹⁵ G. The fastest growing mode of MRI in ideal MHD (Balbus & Hawley 1991) has a wavelength of

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{MRI}}=\frac{8\pi }{q_0(4-q_0){\mathrm{\Omega }}_0}\frac{B_0}{\sqrt{4\pi {\rho }_0}}·\end{array}$$(9)

It ranges from λ_MRI = 1.9 × 10⁵ cm to λ_MRI = 1.0 × 10⁶ cm and is well resolved with our resolution of Δ_r ≡ D/N_r = 7.32 × 10³ cm.

2.4. Boundary conditions

We assumed non-penetrating boundary conditions (u_r = 0). At the inner boundary, we used a standard no-slip condition with an inner core in solid body rotation at a rate Ω_i that evolves with the viscous torque (i.e., u_ϕ = Ω_ir_i and u_θ = 0). At the outer boundary, we used a modified no-slip condition where we forced u_ϕ to match the initial rotation profile at all times (and u_θ = 0). The outer boundary is therefore not in solid body rotation. For the magnetic field, we compared three different boundary conditions: pseudo-vacuum (imposing a radial field: B × n = 0, where n is the normal vector of the outer surface), perfect conductor (imposing a tangential field: B ⋅ n = 0), and insulating (matching a potential field outside the domain).

2.5. Physical parameters and dimensionless numbers

We chose the physical parameters of the simulations to represent a fast rotating PNS model similar to the study of Guilet et al. (2015). All of our models have a uniform viscosity of ν = 7.03 × 10¹¹ cm s⁻². Given the values for the viscosity, the rotation rate Ω₀ and the characteristic length D (Sect. 2.1), the dimensionless Ekman number (characterizing the importance of viscosity over Coriolis force) is

$$\begin{array}{c}\hfill E\equiv \frac{\nu }{{\mathrm{\Omega }}_0{D}^2}=2\times {10}^{-4}.\end{array}$$(10)

For a PNS, the viscosity is large due to the impact of neutrinos while the resistivity is very small, which leads to a high magnetic Prandtl number P_m ≡ ν/η. A realistic value of the magnetic Prandtl number for a PNS is P_m ≈ 10¹³ according to Thompson & Duncan (1993) and Masada et al. (2007). However, numerical simulations for these values of P_m are not possible due to numerical constraints. We varied the resistivity from η = 4.39 × 10¹⁰ cm s⁻¹ to η = 1.66 × 10¹¹ cm s⁻¹, which corresponds to magnetic Prandtl numbers in the range P_m ∈ [4, 8, 16]. Our standard value P_m = 16 is a compromise between the wish to be in the high magnetic Prandtl regime and the computing time constraints. The corresponding magnetic Reynolds number that characterizes the relative importance of magnetic advection to magnetic diffusion (resistivity) is

$$\begin{array}{c}\hfill R_{\mathrm{m}}\equiv \frac{D^2{\mathrm{\Omega }}_0}{\eta }\in [2\times {10}^4,4\times {10}^4,8\times {10}^4].\end{array}$$(11)

3. Results

3.1. Typical quasi-stationary dynamos

We first describe the results from one fiducial simulation where we obtain a self-sustained dynamo. For the model Standard (see Table A.1), we applied insulating boundary conditions and initialized the magnetic field with the parameters B₀ = 8.8 × 10¹⁴ G and [L_min,  L_max]=[0.23 r_o,  0.38 r_o] (see Sect. 2.3).

Figure 3 shows the temporal evolution of the volume-averaged toroidal and poloidal magnetic energy densities and the turbulent kinetic energy density. In order to separate the MRI-driven turbulent flow from the differential rotation, the turbulent kinetic energy density is computed by subtracting the contribution of the axisymmetric azimuthal velocity from the averaged kinetic energy density. After approximately 400 ms, we obtain a statistically stationary state with a mean magnetic field intensity B = 2.5 × 10¹⁵ G. The main contribution is from the toroidal magnetic field, which is ∼2 times larger than the poloidal magnetic field. The magnetic field is predominantly nonaxisymmetric, the axisymmetric magnetic field being ∼3 times weaker than the total magnetic field. The averaged magnetic energy density is more than ten times stronger than the turbulent kinetic energy density. A similar ratio is observed in local models of MRI-driven turbulence (see Sect. 4.3).

Fig. 3. Temporal evolution of the magnetic and turbulent kinetic energy density for the model StandardB0 = 8.8 × 1014 G, [Lmin,  Lmax]=[0.23 ro,  0.38 ro]. The blue and green lines are the toroidal and poloidal contributions of the magnetic energy density, while the orange line is the turbulent kinetic energy density (axisymmetric toroidal contribution is removed). The blue dotted line is the axisymmetric contribution to the toroidal magnetic energy density.

Figures 4 and 5 are representative snapshots of the quasi-stationary state and illustrate the complex geometry of the magnetic field due to MRI-driven turbulence. As is common in this turbulence, the winding of the magnetic field by the shear produces elongated structures in the azimuthal direction, which are clearly seen in the equatorial plane (Fig. 4). On the meridional cuts, the magnetic field is distributed on smaller scales. The turbulent structure of the magnetic field has lost memory of the initial magnetic configuration. Near the rotation axis, the magnetic field is not very intense and the MRI-driven turbulence is weaker for cylindrical radii smaller than s_turb ≈ 9.4 km. In the most turbulent zone, the turbulent radial velocity field has similar small-scale structures, while weak large-scale flows develop near the rotation axis (bottom panel of Fig. 5). The weak magnetic field and the non-turbulent flow near the rotation axis are expected because the vertically averaged shear rate of the simulation is small for s ≤ 9 km and the inner part of the radial profile rotates as a solid body (blue line in Fig. 6).

Fig. 4. Top: 3D snapshot of the intensity of the magnetic field at t = 600 ms for model Standard. The colors represent the magnetic field amplitude from weak (blue) to strong (red). Bottom: 3D snapshot of the magnetic field lines at t = 600 ms for the model Standard. The colors represent the magnetic field amplitude.

Fig. 5. Slices for ϕ = 0 and t = 600 ms for the model Standard. Top: azimuthal magnetic field Bϕ. Bottom: turbulent radial velocity ur (i.e., the axisymmetric contribution is subtracted from the radial velocity).

Fig. 6. Comparison of the vertically and azimuthally averaged shear rate q inside the PNS for different MHD simulations as a function of the cylindrical radius. The values are time-averaged in the quasi-stationary phase of the dynamo. The black line represents the radius of the inner core ri.

To understand how the magnetic and kinetic energies are distributed over different scales, we computed the axisymmetric and nonaxisymmetric components of the toroidal and poloidal spectra (Fig. 7). To first order, the nonaxisymmetric toroidal magnetic spectrum as a function of the spherical harmonics order l (top panel of Fig. 7) can be decomposed in two parts: an increasing part up to order l ≃ 20 and a sharp decrease for the small scales where the dissipation occurs. The axisymmetric toroidal contribution dominates for the largest scales at l <  3, while it is clearly subdominant for the small scales. The quadrupole mode (l = 2) is particularly strong for this model but this varies between different models. Contrary to the toroidal component, the poloidal magnetic field is dominated by its nonaxisymmetric component at all scales. The nonaxisymmetric poloidal magnetic spectrum is similar to the toroidal one but peaks at larger scales for l ≃ 10. These spectra show that the main contribution to the total magnetic energy comes from the toroidal component at intermediate scales (l ∼ 40).

Fig. 7. Top: spectrum of the toroidal magnetic energy (blue) and the poloidal magnetic energy (red) normalized by the total magnetic energy as a function of the spherical harmonics order l. The dotted lines correspond to the axisymmetric components of these energies. Bottom: spectrum of the toroidal kinetic energy (blue) and the poloidal kinetic energy (red) normalized by the total kinetic energy as a function of the spherical harmonics order l. The dotted lines correspond to the axisymmetric components. The black dashed line corresponds to a scaling of l−1.

The nonaxisymmetric poloidal and toroidal kinetic spectra are similar (bottom panel of Fig. 7) and can also be decomposed in two parts: a power law for the large scales and a sharp decrease for the small scales. The power law at larger scales (l <  35) seems to match a scaling of l⁻¹, which is consistent with the kinetic spectrum of high Reynolds local simulations (Fromang 2010). The axisymmetric toroidal component of the kinetic spectra is composed of the turbulence contribution dominating at small scales and the differential rotation dominating at large scales. The oscillations between odd and even modes observed for l <  20 are due to the symmetry of the differential rotation with respect to the equatorial plane. In the same way, the axisymmetric poloidal contribution contains both turbulence (for l >  20) and the meridional circulation (for l <  20), which originates from the interplay between angular momentum transport and the outer boundary forcing. The total kinetic spectrum is dominated by the axisymmetric differential rotation and the meridional circulation at the large scales (l <  10), while the nonaxisymmetric kinetic turbulence dominates at the intermediate and small scales.

Since the magnetar timing parameters only constrain the dipolar component of the magnetic field, we focus in more details on this specific mode. The dipole field strength reaches a significant intensity of B_dip = 1.25 × 10¹⁴ G but it is not the dominant mode in the simulation as it is approximately 20 times weaker than the total magnetic field intensity (see blue and black lines in Fig. 8). This is consistent with a visual inspection of the snapshots, where a dipole structure of the magnetic field cannot be seen.

Fig. 8. Temporal evolution of the averaged magnetic energy density and dipolar energy density for the model Standard. The black line is the averaged magnetic energy density, the blue dash-dotted line is the averaged dipole energy density and the green dash-dotted is the axial dipole energy density.

In Fig. 8, the intensity of the axial dipole is ∼4 times lower than the averaged dipole intensity, implying that the dipole is tilted toward the equatorial plane. We computed the dipole tilt angle θ_dip from the magnetic dipole moment μ given by the formulas

$$\begin{array}{cc}& \mathit{\mu }=\frac{1}{2}{\displaystyle \int \int \int }\mathit{r}\times \mathit{J}\mathrm{d}V,\hfill \end{array}$$(12)

$$\begin{array}{cc}& {\theta }_{\mathrm{dip}}={90}^{°}-arctan\left(\frac{{\mu }_z}{\sqrt{{\mu }_x^2+{\mu }_y^2}}\right),\hfill \end{array}$$(13)

where J = ∇ × B/μ₀ is the current and the volume integral covers the entire numerical domain. Figure 9 shows the time series of the dipole tilt angle. Its time-averaged value is θ_dip = 120° for the run Standard. This result is consistent with the ratio of the energy contained in the total dipole and axial dipole. The equatorial character of the dipole may be explained in the following way: The azimuthal magnetic field does not contribute to the axial dipole but it can contribute to the equatorial dipole through its m = 1 component. The fact that the MRI generates a predominantly azimuthal magnetic field may therefore explain the dominance of the equatorial dipole. The peaks we see in the dipole tilt angle evolution are due to the weakening of the dipole moment in the equatorial plane for short periods of time.

Fig. 9. Temporal evolution of the dipolar tilt angle in degrees for the model Standard. The aligned dipole corresponds to θdip = 0° and the equatorial dipole corresponds to θdip = 90°.

3.2. Dynamo threshold

We showed in the previous subsection that a MRI-driven dynamo can reach a quasi-stationary state but the turbulence can also be damped by the diffusion processes depending on the diffusivities or the initial magnetic field strength. In fact, diffusion processes (viscosity and resistivity) tend to limit the growth of MRI modes for weak fields (Fromang et al. 2007; Guilet et al. 2015). Therefore, we expect the dynamo threshold to depend on the initial magnetic field intensity B₀ and the magnetic Prandtl number P_m (in this study, the Ekman number E is the same for all simulations).

Figure 10 shows the time evolution of the magnetic and turbulent kinetic energies for different values of B₀. Three behaviors are observed: For the lowest B₀ (dotted line), the turbulent energies decrease sharply with time, which indicates that no dynamo is achieved. For intermediate B₀ (dashed line), the turbulent energies have a slow decreasing phase followed by a fast drop. This behavior is called a transient dynamo in the following sections. For higher B₀, the turbulent energies reach a quasi-steady state, with some fluctuations around a non-evolving average. In this case, the simulations reach a self-sustained quasi-stationary dynamo.

Fig. 10. Temporal evolution of the magnetic (black) and kinetic (red) turbulent energy densities for different values of B0 with [Lmin,  Lmax]=[0.23 ro,  0.38 ro] and Pm = 16.

Figure 11 shows that a self-sustained dynamo can be obtained for P_m ≳ 12 ± 4. A simulation with P_m = 16 gives a quasi-stationary dynamo, P_m = 8 a transient dynamo, and P_m = 4 no dynamo. This result is independent of the value of B₀ for B₀ ≥ 10¹⁵ G. It should be noted that the threshold of the dynamo action is set for relatively strong initial magnetic fields because of the high values of viscous and magnetic diffusivities. For lower diffusivities relevant to a PNS, we would expect MRI dynamo to set in for lower initial magnetic fields.

Fig. 11. Same as Fig. 10 for three different values of Pm and with [Lmin,  Lmax]=[0.75,  1] and B0 = 1.8 × 1015 G.

3.3. Robustness of the results

The previous subsections demonstrated the ability of self-sustained MRI-driven dynamos to generate a significant, though subdominant, magnetic dipole. In this section, we assess the robustness of these results by studying the impact of initial and boundary conditions of the magnetic field. We show that the quasi-stationary state is independent of the initial conditions and that the boundary conditions only have a minor impact.

We varied the initial magnetic field amplitude B₀ and its minimum length scale L_min and display the results in Fig. 12. To characterize the early phase of magnetic amplification, we measured the early maximum of the toroidal magnetic field after a few milliseconds (i.e., the first local maximum value reached in the time series, see for example Fig. 11). For a given smallest initial magnetic length scale L_min (indicated by the symbol color), the early maximum toroidal magnetic field is clearly increasing with the initial magnetic amplitude by a factor of ∼2. This dependence is shallower than in Rembiasz et al. (2016), probably because our initial magnetic field is on smaller scales. On the other hand, for a given B₀, the early maximum toroidal magnetic field increases by a factor of ∼5 when L_min is increased. The parameter L_min seems to have the most prominent effect on the early maximum toroidal magnetic field. A possible interpretation is that a small-scale magnetic field is more prone to dissipation. Finally, Fig. 12 suggests that the same initial conditions give approximately the same early maximum magnetic field, independently of the magnetic boundary conditions. Since this early maximum toroidal magnetic field is taken at the beginning of the simulation, boundary conditions may have a minor effect at this stage.

Fig. 12. Early maximum toroidal magnetic field as a function of initial magnetic field strength B0. The symbol color describes the smallest initial magnetic length scale Lmin. The symbol shape represents the magnetic boundary conditions, between perfect conductor (diamond), insulating (circle) and pseudo-vacuum (star).

To characterize our self-sustained dynamos, we performed time and volume averages on the magnetic field once the simulation has reached a quasi-stationary state (Fig. 13). The averaged magnetic field now displays much less variation than the early maximum. Indeed, all models lie between 1.8 × 10¹⁵ G and 2.7 × 10¹⁵ G, with a mean of B = (2.27 ± 0.23) × 10¹⁵ G. Contrary to Fig. 12, no systematic impact of the initial conditions can be observed as different simulations give similar magnetic fields. For a given boundary condition, the small variations can be explained by the stochasticity of MRI-driven dynamos. On the other hand, the magnetic field is slightly stronger in the case of perfect conductor boundary conditions. This trend may be explained by a stronger toroidal magnetic field close to the outer boundary, since the perfect conductor condition allows for a non-vanishing toroidal magnetic field. Overall, these results indicate that the initial conditions and boundary conditions have a small impact on the global properties of the final quasi-stationary state.

Fig. 13. Magnetic field strength averaged in the quasi-stationary state as a function of initial magnetic field strength B0. The symbol color describes the smallest initial magnetic length scale Lmin. The symbol shape represents the magnetic boundary conditions. This figure shows that the quasi-stationary state is roughly independent of the initial and boundary conditions, contrary to the early maximum.

It is important to assess to which extent our setup impacts the magnetic field morphology, and particularly the dipole component that is constrained by observations. The normalized magnetic and kinetic spectra depending on the Legendre polynomial order l (Fig. 14) and the azimuthal degree m (Fig. 15) show the differences at small and large scales for different models. The magnetic and kinetic spectra coincide well for the intermediate and small scales (top and bottom panel of Fig. 14), while the magnetic spectrum displays significant stochasticity at the large scales. For example, some models have a higher magnetic quadrupole l = 2, while other models have a maximum for a higher order l. It seems also that the models with perfect conductor boundary conditions (PC and PCBIS) have a lower kinetic l = 1 mode. As described in Sect. 3.1, the oscillations in the kinetic spectra are due to the meridional circulation. The magnetic spectra as a function of the azimuthal degree m (Fig. 15) is dominated by the large scales with a constant energy for m ≲ 10 and decreases sharply for m ≳ 10. If we compare the l and m spectra, we note that the dominating scales are larger for the azimuthal spectrum, which is consistent with the elongated structures we can see in the equatorial plane of Fig. 4. Similarly, the kinetic energy decreases more steeply as a function of the azimuthal degree m than as a function of l. All models show a very good agreement for these kinetic and magnetic spectrum.

Fig. 14. Top: spectrum of the magnetic energy as a function of the spherical harmonics order l. Bottom: spectrum of the poloidal kinetic energy as a function of the spherical harmonics order l. Each line is from a different model. The dotted line corresponds to a scaling of l−1.

Fig. 15. Magnetic (top) and kinetic (bottom) normalized energy spectra as a function of the azimuthal degree m for different models.

Even though the magnetic spectra of different simulations show significant dispersion at large scales, the dipole amplitude scales linearly with the averaged magnetic field (Fig. 16). The generation of a dipolar magnetic field is therefore a robust feature of the MRI and both initial and boundary conditions have a small impact on our qualitative and quantitative results, such as the dipole and the total magnetic field.

Fig. 16. Dipole magnetic field strength as a function of the averaged magnetic field strength. Symbols are defined in the caption of Fig. 12.

4. Comparison to local simulations

In this section, we would like to understand the impact of the geometry of the domain on the properties of the MRI. As our study is one of the first global spherical models that resolves the MRI, it is important to compare our results to local model results. We intend to establish whether a local Cartesian box can reproduce faithfully some of the properties of the MRI in a global spherical shell. A schematic view of the comparison is shown on Fig. 17. The Cartesian box is meant to represent the differentially rotating zone in the equatorial region of our global model. The Cartesian coordinates (x, y, z) in the local model corresponds to the cylindrical coordinates (s, ϕ, z) at the equator of the global model. The center of the box is at a radius of s₀ ∼ 0.7 r_o. For the cylindrical radii s_cyl ∈ [0.525 r_o,  0.9 r_o], the shear rate q is approximately constant and takes varying values in the range [0.7,  1.0] for different models (see Fig. 6).

Fig. 17. Toroidal magnetic field Bϕ of the local model Ls03Lp09Lz09 with a box size (Ls,  Lϕ,  Lz) = (0.3 ro,  0.9 ro,  0.9 ro) compared to the spherical domain of our global model.

4.1. Setup for local models

The local models were computed with the code SNOOPY, which has been widely used to study the MRI (Lesur & Longaretti 2005, 2007; Guilet & Müller 2015, and references therein). The rotation rate Ω₀ and the diffusivities in the local model are equal to the values used in the global model. The initial magnetic field was initialized with a sinusoidal radial profile of vertical field, which has been chosen such that the magnetic flux vanishes similarly to the global model. We have checked that the quasi-stationary turbulence is independent of the initial condition, provided that the magnetic field is strong enough to initiate an MRI dynamo.

To investigate the effect of the size of the domain and its geometry, we chose different box sizes (L_s,  L_ϕ,  L_z). If the properties of the MRI dynamo were independent of the domain geometry, then we would expect similar results when models have the same volume. The box (L_s,  L_ϕ,  L_z) = (0.5 r_o,  4.5 r_o,  1.2 r_o) matches the volume of the most turbulent region V_turb in the global model where differential rotation is strong, which corresponds to a cylindrical radius s_cyl ≥ 0.5 r_o. Models with reduced box dimensions, especially the azimuthal length L_ϕ, have also been considered in an attempt to obtain results more similar to the global models. The volume of the box for the different box dimensions considered ranges from V = 0.04 V_turb to V = V_turb (see Table 1).

4.2. Dynamo threshold

For comparison with the results of Sect. 3.2, we determined the threshold in magnetic Prandtl number P_m necessary for dynamo action as a function of the azimuthal length L_ϕ. Figure 18 compares the dynamo threshold in the global model (right part) and the local model (left part). For the global model we obtain no dynamo for P_m ≤ 4, a transient dynamo for P_m = 8 and a stable dynamo for P_m ≥ 16. The local model that matches this global behavior corresponds to a box size of (L_s,  L_ϕ,  L_z) = (0.3 r_o,  0.9 r_o,  0.9 r_o) and the P_m threshold decreases for larger boxes. This shows that the dynamo threshold is higher in the global model than the local model with the same volume.

Fig. 18. Summary of the state (dynamo or not) of the model in a (Lϕ, Pm) plane for both local and global models. A green circle corresponds to a self-sustained dynamo. An orange triangle corresponds to a transient dynamo. A red cross corresponds to no dynamo.

4.3. Turbulent energies and angular momentum transport

To make a more detailed comparison between the local and the global models, we now consider the turbulent transport of angular momentum as well as the kinetic and magnetic energies. For comparison with the local models, the quantities from the global models were averaged for the cylindrical radius in the interval s_cyl ∈ [0.525 r_o,  0.925 r_o]. We excluded a thin layer of thickness 0.075 r_o near the outer boundary, where the turbulence is strongly affected by the boundary conditions.

Figure 19 shows the magnetic energy density as a function of the azimuthal length L_ϕ for different magnetic Prandtl number P_m, and different box lengths L_s and L_z. All local models have a significantly higher energy density than the global models (gray area), with the smaller boxes having a magnetic energy density closer to the global models. The impact of the box dimensions can also be deduced from these results. First, for a same radial length L_s (same color) and P_m (same marker), the magnetic energy density increases with the azimuthal length L_ϕ. Moreover, the magnetic energy density increases with the magnetic Prandtl number P_m, which was expected from previous local MRI studies (Lesur & Longaretti 2007; Fromang et al. 2007; Longaretti & Lesur 2010). Figure 19 shows that the global simulations tend to have a lower averaged magnetic energy density than the local simulations.

Fig. 19. Averaged magnetic energy density as a function of the box azimuthal size. The symbol shape represents the magnetic Prandtl number Pm. The symbol color indicates the radial length of the box (brown: Ls = 0.67; orange: Ls = 0.4). The gray zone represents the range of magnetic energy density averaged in the turbulent zone for our sample of global models.

In MRI turbulence, ratios of different time averaged quantities can be more robust, showing less variation than the turbulent energies. For example, the ratio of the magnetic to kinetic turbulent energies is around 10 for both local and global models (see Table 1). One can also look at the turbulent stresses and the different component contributions to the turbulent energies to compare in more details the local and global models (Table 1). For this diagnostic, we used the sϕ component of the Maxwell and Reynolds stress tensors, respectively defined as

$$\begin{array}{cc}& M_{\mathrm{s}\varphi }=-\frac{1}{{\rho }_0{\mu }_0}\langle \frac{B_{\mathrm{s}}{B}_{\varphi }}{D^2{\mathrm{\Omega }}_0^2}\rangle ,\hfill \end{array}$$(14)

$$\begin{array}{cc}& R_{\mathrm{s}\varphi }=\langle \frac{u_{\mathrm{s}}{u}_{\varphi }}{D^2{\mathrm{\Omega }}_0^2}\rangle ,\hfill \end{array}$$(15)

where the brackets ⟨ ⋅ ⟩ corresponds to the volume averaged value.

The Maxwell and Reynolds stresses normalized by the magnetic energy are weaker in global models, which shows that the angular momentum transport is less efficient. This difference can be explained for the Maxwell stress by the lower value of the magnetic radial contribution b_s, while the other contributions are similar in proportion. The difference is also reduced for the local models that have a low magnetic energy.

Figure 20 shows that there is a correlation between the radial and vertical components of the magnetic field for physically different models, such as our local and global incompressible models, compressible stratified local and global models from an accretion disk study (Hawley et al. 2011) and stratified and unstratified local models (Shi et al. 2010, 2016). Our local models seem to match well this trend, while our global models have a slightly lower radial contribution. It is unclear whether this small difference hints at an underlying systematic effect or is instead consistent with an inherent spread in the distribution of the data.

Fig. 20. Comparison for different models of the radial component $B_{\text{s}}^2$ as a function of vertical component $B_z^2$. Both components are normalized by the azimuthal component $B_{\varphi }^2$. The different values are taken from Hawley et al. (2011) and Shi et al. (2010, 2016).

Among local models, Ls03Lp09Lz09 is the one showing features most similar to the global models (see the 3D snapshot in Fig. 17 and the one shown in Fig. 4 for a comparison). In addition to sharing the same threshold in P_m for dynamo action (Sect. 4.2), it has the smallest magnetic energy among the local models at P_m = 16. Its kinetic and magnetic spectra can therefore be compared to the global nonaxisymmetric spectra (Fig. 21). The non-radial wavenumber $\stackrel{\sim }{k}$ is related to the spherical harmonic order l by the relation $\stackrel{\sim }{k}=\sqrt{l(l+1)/r^2}$ with r taken as r = 0.75 r_o for the global model in order to simplify the comparison. The magnetic spectrum of the local model is similar to the global spectrum. It can also be decomposed in an approximately flat part at large scales and an exponential decrease due to dissipation at small scales. The local kinetic spectrum is also in good agreement with the nonaxisymmetric spectrum global model. The local magnetic and kinetic spectra are stronger than the global ones at small scales, which is consistent with the stronger turbulence in our shearing boxes leading to slightly smaller dissipative scales. The main difference between the global and local spectra is therefore the presence of scales larger than the box, especially the axisymmetric component of the kinetic spectrum. Overall, these results suggest that the self-sustained dynamo obtained in the global models is typical of the MRI in local simulations.

Fig. 21. Spectrum of the magnetic energy and kinetic energy as a function of the wavenumber $\stackrel{\sim }{k}$. The blue spectrum is averaged on the full duration of the local model Ls03Lp12Lz12 and the red spectrum is the total nonaxisymmetric spectrum of the global model Standard (Fig. 7).

5. Discussion

In this section, we discuss the influence of the domain geometry and the different limits of the simplifying assumptions made in our study: the impact of the curvature, the mechanical boundary conditions, the diffusive processes and the incompressible approximation.

5.1. Impact of the curvature

As the local and global models are designed to have a similar differential rotation, one might have naively expected to obtain quantitatively similar results. However, the different geometry of the domain is sufficient to obtain different results between the local and global models. Our interpretation is that the curvature of the sphere and non-periodic boundaries reduce the coherence of the magnetic field and velocity field. Thus, local models and large boxes in particular tend to overestimate the field amplification, since they allow for larger coherence structures to develop and thus favor the MRI action over a broader range of scales.

5.2. Forcing of the differential rotation

We chose to rely on the outer mechanical boundary to force the differential rotation in the global simulations, which is justified by the dynamical evolution of the PNS in core-collapse supernovae. Indeed, the PNS contraction accelerates its rotation. At the same time, the MRI and the MHD turbulence transport the angular momentum toward the outer cylindrical radii and therefore slow down the rotation. This behavior may even last longer than the duration of the core-collapse supernova because the contraction of the PNS can still continue after the explosion and the rotation can slow down due to winds or magnetic braking. These processes lead to the formation of two distinct regions within the PNS: The inner cylindrical part of the PNS is accelerated, while the outer cylindrical part is slowed down by the angular momentum transport. Our boundary conditions try to mimic this effect. We note that, despite being fixed at the boundary, the rotation profile inside the domain is free to evolve as a consequence of the backreaction of the magnetic field. As illustrated in Fig. 6, we indeed observe that a stronger magnetic field leads to a lower degree of differential rotation.

However, one of the main limitations of this boundary condition is that the evolution of the global rotation profile is not well described as the effects of contraction, accretion, magnetic braking and transport of angular momentum is artificially taken into account on the outer boundary. Indeed, the timescale over which the rotation profile is significantly changed can be estimated by the amount of shear energy extracted by the combined effect of Reynolds and Maxwell stresses in the turbulent region. We estimate a shear energy of 3.9 × 10⁵⁰ erg by comparing the rotational energy contained in the differentially rotating PNS to that of a PNS in solid body rotation with the same total angular momentum. The typical mean rate at which the shear energy is extracted is ≈5.2 × 10⁵¹ erg s⁻¹. The typical timescale at which the energy is extracted from the shear is therefore ≈75 ms. This means that differential rotation would decrease rapidly without contraction and accretion. In MHD core collapse simulations, it remains unclear whether the speed-up due to contraction and accretion or rather the action of stresses and magnetic braking have a stronger impact on the rotation profile.

The energy injection by this outer boundary condition also depends on the viscous layer, and therefore on the diffusion processes. In addition, the energy injection of the outer boundary is different from the constant shear used in local models, which might be an other way to explain the differences between local and global models.

5.3. Diffusive processes

One could be tempted to solve our set of equations in the ideal MHD limit. However, resolution studies have shown that the turbulent transport is then sensitive to the numerical resolution (Pessah et al. 2007). Therefore, in order to avoid the pitfall of resolution-dependent results, we used explicit diffusivities and made sure that the diffusive scales were resolved. The low numerical diffusion of the pseudo-spectral codes ensure that only physical diffusivities impact the dynamics.

The regime of the diffusive parameters that can be numerically reached is a limit to our model: A realistic value of the neutrino viscosity is ν ∼ 2 × 10¹⁰ cm² s⁻¹ at a radius of ∼20 km (Guilet et al. 2015). The corresponding Ekman number is then E ∼ 6 × 10⁻⁶, which is lower than in our simulations by a factor of about 30. Using a more realistic value of E may lower the thresholds to obtain a dynamo (see Sect. 3.2) and may lead to a stronger amplification of the magnetic field. Moreover, the resistivity used in our simulations is much larger than realistic values. Due to numerical constraints, the magnetic Prandtl number is limited to 16 in our study, much smaller than realistic values of P_m ∼ 10¹³. Accretion disk studies have shown that the MRI-driven turbulence is highly sensitive to the magnetic Prandtl number, and generally it increases with P_m (Lesur & Longaretti 2007; Fromang et al. 2007; Longaretti & Lesur 2010). Our local numerical simulations confirm this trend (see Table 1). Our quantitative results should therefore be considered as lower bounds on the magnetic field strength that would be reached at higher P_m. Investigating the regime of higher P_m in a global model is computationally demanding and it would certainly be more affordable to use local models. Further comparisons between local and global simulations would then allow us to extrapolate the results to a global model.

Finally, we assume that neutrinos are in the diffusive regime, which should be valid only in the inner part of the PNS. In the outer parts of the PNS, the neutrinos are in the nondiffusive regime on length scales at which the MRI grows. The impact of the neutrino drag regime on the linear phase has been studied by Guilet et al. (2015), but has never been studied in nonlinear numerical simulations. The evolution of the MRI in this regime remains an open question.

5.4. Incompressible approximation

The incompressible approximation was used to provide an idealized reference model and to reduce the cost of our numerical simulations by filtering out sound waves. First, filtering sound waves is justified when the sound speed c_s is much larger than both the fluid and Alfvén velocities (${\mathit{v}}_{\mathrm{A}}\equiv B/\sqrt{{\mu }_0\rho }$). With c_s ∼ 5 × 10⁴ km s⁻¹ at r ∼ 20 km, we check a posteriori that this condition is satisfied in our global models, for which $u^2/c_{\text{s}}^2\le v_{\text{A}}^2/c_{\text{s}}^2\le {10}^{-4}$. This is consistent with the discussion of the Boussinesq approximation in Guilet et al. (2015). Second, we neglect the composition and entropy gradients. The incompressible approximation could be extended to take into account the buoyancy within the Boussinesq approximation. However, we do not include it for the sake of simplicity and in order to be able to compare our results to most of the numerical studies of the MRI in the literature. Finally, density stratification is neglected and its study is postponed to a further work.

6. Conclusions

For the first time, we have investigated the generation of large-scale magnetic fields by the MRI in an idealized 3D spherical model of PNS. We performed a parameter exploration of the initial and boundary conditions and we compared in details the results of our global simulations with local shearing box simulations. The main findings of our study can be summarized as follows:

• The MRI leads to a quasi-stationary state with a strong turbulent magnetic field of B ≥ 2 × 10¹⁵ G. The toroidal component of the magnetic field dominates over the poloidal field by a factor of ∼2. A nondominant magnetic dipole B_dip ∼ 10¹⁴ G is generated, which represents about 5% of the averaged magnetic field strength. Interestingly, this dipole is tilted toward the equatorial plane with a tilt angle in the interval θ_dip ∈ [60 ° ,120 ° ].

• The magnetic field amplification and dipole generation by the MRI is a robust mechanism that operates for a large number of different initial setups and magnetic boundary conditions. Indeed, the turbulent energies and the angular momentum transport in the quasi-stationary state do not depend much on the initial magnetic field and boundary conditions.

• The comparison between global and local studies suggests that the geometry of the domain leads to a decrease in the turbulent energies. Indeed, global models have lower turbulent energies than local models, although the ratio of kinetic and magnetic energies and the relative contributions to these energies are comparable. The results from the small boxes of the local simulations have a better agreement with those from the global models, which may be interpreted by the fact that the curvature in a spherical shell limits the coherence length scale of the turbulence.

These findings support the ability of the MRI to generate a strong magnetic dipole and explain magnetar formation. One should of course keep in mind that our results still need to be confirmed by setups with a more realistic PNS structure as discussed below.

The magnetic dipole amplitude obtained from the time and volume averages ranges from 9.1 × 10¹³ G to 1.35 × 10¹⁴ G, which is within the lower end of the observed range for the dipole of galactic magnetars. We further note that these results are obtained before the full contraction to the final size of a cold neutron star, which has a radius close to 12 km. If the magnetic flux is conserved during this contraction, we expect the magnetic field to be amplified by a factor of ≃4, possibly bringing the dipolar component of the magnetic field in the middle of magnetar range 10¹⁴ − 10¹⁵ G.

To compare our results directly with observations, the turbulent magnetic field should also be relaxed to a stable state without differential rotation and be evolved on longer timescales. Under the joint effects of magnetospheric dissipation and internal dissipation, the dipole from our model could become a stable equatorial dipole (Lander & Jones 2018).

Similarly to our results, some observations suggest that the dipole may not be the strongest component of the magnetic field of magnetars. Indeed, Makishima et al. (2016, 2019) have detected phase modulations in X-ray emissions from two magnetars, which are interpreted as a free precession due to prolate deformations of the neutron star under the torque of a very intense internal toroidal field. The MRI-driven magnetic field is mainly toroidal and could be important to explain these observations.

As discussed above, our results suggest that a fast rotating PNS can become a typical magnetar through MRI action. A slower rotation, on the other hand, is expected to lead to a lower magnetic field and may explain the formation of low field magnetars (Rea et al. 2012, 2013, 2014). In fact, owing to the simplicity of our setup, our simulations can be rescaled to a different rotation frequency. For example, for a 30 times slower rotation rate³, the magnetic field should be scaled down by a factor of 30. The final neutron star would have a magnetic dipole of ≃1 − 2 × 10¹³ G and a total magnetic field of ≃3 × 10¹⁴, which corresponds to a low field magnetar. The ratio of the magnetic dipole to the total magnetic field can also be compared to the cases of SGR 0418+5729 (Tiengo et al. 2013) and J1822.3−16066 (Rodríguez Castillo et al. 2016), for which the total magnetic field can be measured using the proton cyclotron line. The observed ratios are respectively B_dip/B_tot ∼ 0.0012 − 0.024 and B_dip/B_tot ∼ 0.01 − 0.05. These are roughly in agreement with our results B_dip/B_tot ∼ 0.05. The slightly higher value found in our models may be due to the high resistivity used in our simulations, since we may expect that a lower resistivity would lead to a smaller-scale magnetic field.

The magnetic field we obtain is strong enough to potentially impact the core-collapse dynamics and launch jet-driven explosions. However, most of the studies of jet driven explosions assume a strong axial dipole field, while our study displays an equatorial dipole. Our results open the question of the impact of a tilted dipole on the dynamics of a core-collapse supernova. The dipolar tilt angle is also a key parameter to determine the amount of thermalized energy in the magnetosphere and the energy available to launch a jet (Margalit et al. 2018). Interestingly, Bugli et al. (2020) showed that higher order multipoles tend to decrease the efficiency of the magnetorotational launching mechanism but can still launch a jet, which suggests that the multipoles of our magnetic field could contribute to the explosion.

We stress that the limitations of this work need to be assessed. Our quantitative results (e.g., strength of the magnetic field and its dipolar component, dipolar tilt angle) may change with the addition of the density and entropy background gradients. A convective dynamo with fast rotation can occur in a realistic interior model of the full PNS, a mechanism that has recently been shown to be able to form magnetars (Raynaud et al. 2020). This raises the open question of the interaction of the MRI in a stably stratified zone with an inner convective dynamo in core-collapse supernovae. The stably stratified region unstable to the MRI is important to link the convective dynamo buried under the PNS surface and the explosion that is launched from the surface.

This study has been focused on the core-collapse supernova context but it is also relevant to binary neutron star mergers, which display similar conditions in terms of differential rotation and neutrino radiation (Guilet et al. 2017). The MRI may grow in the hypermassive neutron star resulting from the merger and the magnetic field may be amplified up to 10¹⁶ G (Siegel et al. 2013; Kiuchi et al. 2018; Ciolfi et al. 2019). These studies start with a strong axial magnetic dipole (≥3 × 10¹⁵ G), while in our case we start with a small-scale field and the MRI generates a strong tilted dipole. Although it is difficult to compare the results obtained with such different initial conditions, our work tends to support the ability of the MRI to generate an magnetar-like dipole in binary neutron star mergers. A detailed study with appropriate thermodynamic background state and differential rotation profile would further assess this possibility. This would support the invoked scenario of magnetars powering an X-ray transient as the aftermath of a binary neutron-star merger (Xue et al. 2019).

Investigating how the magnetic field of magnetars is generated in realistic conditions is a promising avenue of research as new observations will come in the multi-messenger era: The Large Synoptic Survey Telescope (LSST) will be able to observe one hundred times more supernovae (LSST Science Collaboration 2017) and the future spatial telescope SVOM will help us to localize and observe transient events in the γ-ray sky, like electromagnetic counterparts of binary neutron-star mergers (Wei et al. 2016).

Acknowledgments

This research was supported by the European Research Council through the ERC starting grant MagBURST No. 715368. The computations were performed on the supercomputer Occigen of the CINES (Applications A0030410317 and A0050410317). We thank Thomas Gastine, Thierry Foglizzo and Sébastien Fromang for their valuable insight and expertise.

References

Appendix A: Parameters of our global simulations

All Tables

All Figures

	Fig. 1. Rotation profile inside the PNS as a function of the cylindrical radius. The inner region (s <  10 km) is in solid body rotation and the outer region (s >  10 km) follows a power law Ω ∝ r−q0 with q0 = 1.25.
In the text

	Fig. 2. Slice of the Bθ component of the initial magnetic field for the meridional plane ϕ = 0 with [Lmin,  Lmax]=[0.23 ro,  0.38 ro] and B0 = 8.8 × 1014 G.
In the text

Fig. 3. Temporal evolution of the magnetic and turbulent kinetic energy density for the model StandardB0 = 8.8 × 1014 G, [Lmin,  Lmax]=[0.23 ro,  0.38 ro]. The blue and green lines are the toroidal and poloidal contributions of the magnetic energy density, while the orange line is the turbulent kinetic energy density (axisymmetric toroidal contribution is removed). The blue dotted line is the axisymmetric contribution to the toroidal magnetic energy density.

In the text

	Fig. 4. Top: 3D snapshot of the intensity of the magnetic field at t = 600 ms for model Standard. The colors represent the magnetic field amplitude from weak (blue) to strong (red). Bottom: 3D snapshot of the magnetic field lines at t = 600 ms for the model Standard. The colors represent the magnetic field amplitude.
In the text

	Fig. 5. Slices for ϕ = 0 and t = 600 ms for the model Standard. Top: azimuthal magnetic field Bϕ. Bottom: turbulent radial velocity ur (i.e., the axisymmetric contribution is subtracted from the radial velocity).
In the text

	Fig. 6. Comparison of the vertically and azimuthally averaged shear rate q inside the PNS for different MHD simulations as a function of the cylindrical radius. The values are time-averaged in the quasi-stationary phase of the dynamo. The black line represents the radius of the inner core ri.
In the text

Fig. 7. Top: spectrum of the toroidal magnetic energy (blue) and the poloidal magnetic energy (red) normalized by the total magnetic energy as a function of the spherical harmonics order l. The dotted lines correspond to the axisymmetric components of these energies. Bottom: spectrum of the toroidal kinetic energy (blue) and the poloidal kinetic energy (red) normalized by the total kinetic energy as a function of the spherical harmonics order l. The dotted lines correspond to the axisymmetric components. The black dashed line corresponds to a scaling of l−1.

In the text

	Fig. 8. Temporal evolution of the averaged magnetic energy density and dipolar energy density for the model Standard. The black line is the averaged magnetic energy density, the blue dash-dotted line is the averaged dipole energy density and the green dash-dotted is the axial dipole energy density.
In the text

	Fig. 9. Temporal evolution of the dipolar tilt angle in degrees for the model Standard. The aligned dipole corresponds to θdip = 0° and the equatorial dipole corresponds to θdip = 90°.
In the text

	Fig. 10. Temporal evolution of the magnetic (black) and kinetic (red) turbulent energy densities for different values of B0 with [Lmin,  Lmax]=[0.23 ro,  0.38 ro] and Pm = 16.
In the text

	Fig. 11. Same as Fig. 10 for three different values of Pm and with [Lmin,  Lmax]=[0.75,  1] and B0 = 1.8 × 1015 G.
In the text

	Fig. 12. Early maximum toroidal magnetic field as a function of initial magnetic field strength B0. The symbol color describes the smallest initial magnetic length scale Lmin. The symbol shape represents the magnetic boundary conditions, between perfect conductor (diamond), insulating (circle) and pseudo-vacuum (star).
In the text

Fig. 13. Magnetic field strength averaged in the quasi-stationary state as a function of initial magnetic field strength B0. The symbol color describes the smallest initial magnetic length scale Lmin. The symbol shape represents the magnetic boundary conditions. This figure shows that the quasi-stationary state is roughly independent of the initial and boundary conditions, contrary to the early maximum.

In the text

	Fig. 14. Top: spectrum of the magnetic energy as a function of the spherical harmonics order l. Bottom: spectrum of the poloidal kinetic energy as a function of the spherical harmonics order l. Each line is from a different model. The dotted line corresponds to a scaling of l−1.
In the text

	Fig. 15. Magnetic (top) and kinetic (bottom) normalized energy spectra as a function of the azimuthal degree m for different models.
In the text

	Fig. 16. Dipole magnetic field strength as a function of the averaged magnetic field strength. Symbols are defined in the caption of Fig. 12.
In the text

	Fig. 17. Toroidal magnetic field Bϕ of the local model Ls03Lp09Lz09 with a box size (Ls,  Lϕ,  Lz) = (0.3 ro,  0.9 ro,  0.9 ro) compared to the spherical domain of our global model.
In the text

	Fig. 18. Summary of the state (dynamo or not) of the model in a (Lϕ, Pm) plane for both local and global models. A green circle corresponds to a self-sustained dynamo. An orange triangle corresponds to a transient dynamo. A red cross corresponds to no dynamo.
In the text

	Fig. 19. Averaged magnetic energy density as a function of the box azimuthal size. The symbol shape represents the magnetic Prandtl number Pm. The symbol color indicates the radial length of the box (brown: Ls = 0.67; orange: Ls = 0.4). The gray zone represents the range of magnetic energy density averaged in the turbulent zone for our sample of global models.
In the text

	Fig. 20. Comparison for different models of the radial component $B_{\text{s}}^2$ as a function of vertical component $B_z^2$. Both components are normalized by the azimuthal component $B_{\varphi }^2$. The different values are taken from Hawley et al. (2011) and Shi et al. (2010, 2016).
In the text

	Fig. 21. Spectrum of the magnetic energy and kinetic energy as a function of the wavenumber $\stackrel{\sim }{k}$. The blue spectrum is averaged on the full duration of the local model Ls03Lp12Lz12 and the red spectrum is the total nonaxisymmetric spectrum of the global model Standard (Fig. 7).
In the text