A global model of the magnetorotational instability in protoneutron stars

Magnetars are highly magnetized neutron stars whose magnetic dipole ranges from $10^{14}$ to $10^{15}$ G. The MRI is considered to be a promising mechanism to amplify the magnetic field in fast-rotating protoneutron stars and form magnetars. This scenario is supported by many local studies showing that magnetic fields could be amplified by the MRI on small scales. However, the efficiency of the MRI at generating a dipole field is still unknown. To answer this question, we study the MRI dynamo in an idealized global model of a fast rotating protoneutron star with differential rotation. We perform 3D incompressible MHD simulations in spherical geometry with explicit diffusivities where the differential rotation is forced at the outer boundary. We vary the initial magnetic field and investigated different magnetic boundary conditions. These simulations were compared to local shearing box simulations. We obtain a self-sustained turbulent MRI-driven dynamo, whose saturated state is independent of the initial magnetic field. The MRI generates a strong turbulent magnetic field of $B \geq 2\times 10^{15}$ G and a non-dominant magnetic dipole, which represents systematically about $5\%$ of the averaged magnetic field strength. Interestingly, this dipole is tilted towards the equatorial plane. We find that local shearing box models can reproduce fairly well several characteristics of global MRI turbulence such as the kinetic and magnetic spectra. The turbulence is nonetheless more vigorous in the local models than in the global ones. Overall, our results support the ability of the MRI to form magnetar-like large-scale magnetic fields. They furthermore predict the presence of a stronger small-scale magnetic field. The resulting magnetic field could be important to power outstanding stellar explosions, such as superluminous supernovae and GRBs.


Introduction
Magnetars are a class of young and highly magnetized neutron stars that produce a wide variety of outstanding emission at Xray or soft γ-ray energies (Kouveliotou et al. 1998;Kaspi & Beloborodov 2017, and references therein). The period and the spin-down of these objects is measured by the long time followup of their pulsed X-ray activity and a surface dipolar magnetic field of B dip = 10 14 P 5 s 1 2 Ṗ 10 −11 s s −1 1 2 G 10 14 -10 15 G can be inferred under the assumption of magnetic dipole spindown (Olausen & Kaspi 2014) 1 . 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 have 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, e-mail: alexis.reboul-salze@cea.fr 1 http://www.physics.mcgill.ca/~pulsar/magnetar/main. html 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. 2011Metzger et al. , 2018 associated with "hypernovae" or supernovae type Ic BL (Drout et al. 2011). These rare events are characterized by a kinetic energy ten times higher than standard supernovae. Furthermore, millisecond magnetars have been invoked to explain some superluminous supernovae through a delayed energy injection due to the dipole spin down Article number, page 1 of 15 arXiv:2005.03567v1 [astro-ph.HE] 7 May 2020 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 to 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 milliseconds 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 protoneutron star, 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 MHD and with an initial vertical uniform magnetic field, the growth rate of the instability is of 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 time scales than the successful explosion time and that an efficient amplification of the magnetic field occurs at small scales (Obergaulinger et al. 2009;Rembiasz et al. 2016). The influence of the specific physical conditions of protoneutron stars 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 . Finally, the buoyancy forces driven by the entropy and lepton fraction can reduce the MRI turbulence in the case of stable stratification .
The impact of the spherical geometry of the full PNS on the MRI turbulence and the ability of the MRI to generate a largescale 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 have also shown 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 non-linear 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.

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.3M 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 choose 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 ν . The incompressible approximation is used for the sake of simplicity. We assume a uniform density ρ 0 = 4 × 10 13 g cm −3 and neglect buoyancy effects. The diffusive incompressible MHD equations describing the dynamics of the PNS in a rotating frame at an angular frequency Ω 0 = 1000 s −1 read where u is the flow velocity, B the magnetic field, p the reduced pressure (i.e. the pressure normalized by the density) and µ 0 the vacuum permeability.

Numerical methods
In order to solve the incompressible MHD system of equations (1)-(4), we used the pseudo-spectral code MagIC (Wicht 2002;Gastine & Wicht 2012;Schaeffer 2013) 2 . MagIC solves the 3D MHD equations in a spherical shell using a poloidaltoroidal decomposition for the velocity and the magnetic field, 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 non-linear terms and the Coriolis force. The linear terms are computed in the spectral space, while the non-linear 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.

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). We use a cylindrical rotation profile inspired by their results with a central part that rotates like a solid body and an outer part that is differentially rotating with a power law dependency: where s is the cylindrical radius, q 0 corresponds to the shear rate q ≡ − r Ω dΩ dr in the outer part and Ω i is the rotation rate of the inner core. Ω i is computed so that the ratio of total angular momentum over moment of inertia is equal to the frame rotation rate Ω 0 defined in Sect. 2.1. The profile shown in Fig. 1 with q 0 = 1.25 has a smooth transition between solid body rotation when s < 0.4r o and power law differential rotation when s > 0.4r o .
Sustaining this rotation profile for several hundreds of milliseconds with a strong magnetic field is impossible without a mechanism to force the differential rotation. In a core-collapse supernova, the acceleration of the inner layers and the deceleration of the outer layers are due to the contraction of the PNS and the angular momentum transported outside the PNS by the MRI. In our simplified setup, we choose to artificially force the differential rotation by keeping the initial differential rotation profile 2 https://magic-sph.github.io on the outer boundary constant throughout the simulation. This leads to a quasi-stationary state that is simpler to study and to compare with different models.
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;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 b l,m ∝ 0 for r < 7.5 km cos(k r (r − r i )) for r > 12.5 km , where k r is the radial wavenumber and with a smooth transition between 7.5 km and 12.5 km. We select Fourier radial modes and spherical harmonic modes, whose wavelengths λ r = 2π/k r and λ l = r 2 l(l+1) lie within the range [L min , L max ]. 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 The initial root mean square magnetic field strength B 0 is varied from B 0 = 6.31 × 10 14 G to B 0 = 3.36 × 10 15 G. The fastest growing mode of MRI in ideal MHD (Balbus & Hawley 1991) has a wavelength of It ranges from λ MRI = 1.9 × 10 5 cm to λ MRI = 1.0 × 10 6 cm and is well resolved with our resolution of ∆ r ≡ D/N r = 7.32 × 10 3 cm.

Boundary conditions
We assume non-penetrating boundary conditions (u r = 0). At the inner boundary, we use 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 φ = Ω i r i and u θ = 0). At the outer boundary, we use a modified no-slip condition where we force 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), or insulating (matching a potential field outside the domain).

Physical parameters and dimensionless numbers
We choose the physical parameters of the simulations to represent a fast rotating PNS model similar to the study of . All of our models have a uniform viscosity of ν = 7.03 × 10 11 cm s −2 . Given the values for the viscosity, the rotation rate Ω 0 and the characteristic length D (Sect. 2.1), the dimensionless Ekman number (characterizing the importance of viscosity over Coriolis force) is 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 13 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 vary the resistivity from η = 4.39 × 10 10 cm s −2 to η = 1.66 × 10 11 cm s −2 , 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

Typical quasi-stationary dynamos
Let us first describe the results from one fiducial simulation where we obtain a self-sustained dynamo. For the model Standard (see Table A.1), we apply insulating boundary conditions and initialize the magnetic field with the parameters B 0 = 8.8 × 10 14 G and [L min , Figure 3 shows the temporal evolution of the volumeaveraged toroidal and poloidal magnetic energy densities and the turbulent kinetic energy density. In order to separate the MRIdriven 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 milliseconds, we obtain a statistically stationary state with a mean magnetic field intensity B = 2.5 × 10 15 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 non-axisymmetric, 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).
A representative snapshot of the quasi-stationary state with a 3D rendering of the magnetic field amplitude (Fig. 4) and a 2D meridional slice of the toroidal field B φ (top panel of Fig. 5) can capture the complex geometry of the magnetic field due to this Fig. 4. 3D snapshot of the intensity of the magnetic field at t = 600 ms for the model Standard. The colors represent the magnetic field amplitude from weak (blue) to strong (red).
turbulence. As is common in MRI-driven 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 nonturbulent 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).
To understand how the magnetic and kinetic energies are distributed over different scales, we compute the axisymmetric and non-axisymmetric components of the toroidal and poloidal spectra ( Fig. 7). To first order, the non-axisymmetric 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 non-axisymmetric component at all scales. The non-axisymmetric 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).
The non-axisymmetric 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 −1 , which is consistent with the kinetic spectrum of high Reynolds local simulations azimuthal magnetic field B φ . Bottom: Turbulent radial velocity u r , i.e. the axisymmetric contribution is subtracted from the radial velocity. (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/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 Article number, page 5 of 15  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 14 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.
In Fig. 8, the intensity of the axial dipole is ∼4 times lower than the averaged dipole intensity, implying that the dipole is tilted towards the equatorial plane. We compute the dipole tilt angle θ dip from the magnetic dipole moment µ given by the formulas where J = ∇ × B/µ 0 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.

Dynamo threshold
We show 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;. Therefore, we expect the dynamo threshold to depend on the initial magnetic field intensity B 0 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 0 . Three behaviours are observed: for the lowest B 0 (dotted line), the turbulent energies decrease sharply with time, which indicates that no dynamo is achieved. For intermediate B 0 (dashed line), the turbulent energies have a slow decreasing phase followed by a fast drop. This behaviour will be called a transient dynamo in the following sections. For higher B 0 , the turbulent energies reach a quasi-steady state, with some fluctuations around a nonevolving average. In this case, the simulations have reached a self-sustained quasi-stationary dynamo.  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 0 for B 0 ≥ 10 15 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.

Robustness of the results
The previous subsections demonstrated the ability of selfsustained 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 will show that the quasi-stationary state is independent of the initial conditions and that the boundary conditions only have a minor impact.
To characterize the early phase of magnetic amplification, we measure the early maximum of the toroidal magnetic field after a few milliseconds (see for example Fig. 11). We vary the initial magnetic field amplitude B 0 and its minimum length scale L min and display the results in Fig. 12. For a given smallest initial magnetic length scale L min (indicated by the symbol color), the maximum toroidal magnetic field is clearly increasing with the initial magnetic amplitude by a factor ∼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 0 , the early maximum toroidal magnetic field increases by a factor ∼5 when L min is increased. The parameter L min seems to have the most prominent effect on the 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 maximum magnetic field, independently of the magnetic boundary conditions. Since this maximum toroidal magnetic field is taken at the beginning of the simulation, boundary conditions may have a minor effect at this stage.
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 aver- aged magnetic field now displays much less variation than the early maximum. Indeed, all models lie between 1.8 × 10 15 G and 2.7 × 10 15 G, with a mean of B = (2.27 ± 0.23) × 10 15 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 nonvanishing 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.
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.
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.

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.

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, and references therein). The rotation rate Ω 0 and the diffusivities in the local model are equal to the values used in the global model. The initial magnetic field is 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 quasistationary 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 choose different box sizes (L s , L φ , L z ). If the properties of the MRI dynamo were independent of the domain geometry, then we may expect similar results when models have the same volume. The box (L s , L φ , L z ) = (0.5r o , 4.5r o , 1.2r o ) matches the volume of the most turbulent region V turb in the global model where differential rotation is strong, i.e. for a cylindrical radius s cyl ≥ 0.5r 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.04V turb to V = V turb (see Table 1).

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.3r o , 0.9r o , 0.9r 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.

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 are averaged for the cylindrical radius in the interval s cyl ∈ [0.525r o , 0.925r o ]. We exclude a thin layer of thickness 0.075r 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 (grey 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.
In MRI turbulence, ratios of different time averaged quantities can be more robust, showing less variation than the turbu- Magnetic field (10 15 G) Fig. 19. Averaged magnetic energy density as a function of the box azimuthal size. The symbol shape represents the magnetic Prandtl number P m . The symbol color indicates the radial length of the box (brown: L s = 0.67; orange: L s = 0.4). The grey zone represents the range of magnetic energy density averaged in the turbulent zone for our sample of global models. lent 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 use the sφ component of the Maxwell and Reynolds stress tensors, respectively defined as 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/global incompressible models, compressible stratified local/global models from an accretion disk study (Hawley et al. 2011) and stratified/unstratified local models (Shi et al. 2010(Shi et al. , 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. 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 non-axisymmetric spectra (Fig. 21). The non-radial wavenumberk is related to the spherical har- Table 1. Local quantities used as diagnostics of the turbulence in both local and global models. The Maxwell and Reynolds stresses are normalized by the magnetic energy.  monic order l by the relationk = l(l + 1)/r 2 with r taken as r = 0.75r 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 non-axisymmetric 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 ob-tained in the global models is typical of the MRI in local simulations.

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.

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.

Mechanical boundary conditions
We choose 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 protoneutron star in core-collapse supernovae. Indeed, the protoneutron star contraction accelerates its rotation. At the same time, the MRI and the MHD turbulence transport the angular momentum towards the outer cylindrical radii and therefore slow down the rotation. This behaviour may even last longer than the duration of the corecollapse 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 forma-  Fig. 21. Spectrum of the magnetic energy and kinetic energy as a function of the wavenumberk. The blue spectra is averaged on the full duration of the local model Ls03Lp12Lz12 and the red spectra is the total non-axisymmetric spectra of the global model Standard (Fig. 7). tion 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. 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 contraction is artificially taken into account on the outer boundary. 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.

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 use explicit diffusivities and make sure that the diffusive scales are 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 10 cm 2 s −1 at a radius of ∼20 km . The corresponding Ekman number is then E ∼ 6 × 10 −6 , 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 13 . Accretion discs 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 non-diffusive regime on length scales at which the MRI grows. The impact of the neutrino drag regime on the linear phase has been studied by , but has never been studied in non-linear numerical simulations. The evolution of the MRI in this regime remains an open question.

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 (v A ≡ B/ √ µ 0 ρ). With c s ∼ 5 × 10 4 km s −1 at r ∼ 20 km, we check a posteriori that this condition is satisfied in our global models, for which u 2 /c 2 s ≤ v 2 A /c 2 s ≤ 10 −4 . This is consistent with the discussion of the Boussinesq approximation in .
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 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.

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 protoneutron star. 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 15 G. The toroidal component of the magnetic field dominates over the poloidal field by a factor ∼2. A non-dominant magnetic dipole B dip ∼ 10 14 G is generated, which represents about 5% of the averaged magnetic field strength. Interestingly, this dipole is tilted to-wards 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 of 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 dipole and explain magnetar formation.
The magnetic dipole amplitude obtained from the time and volume averages ranges from 9.1 × 10 13 G to 1.35 × 10 14 G, which is within the lower end of the observed range for the dipole of galactic magnetars. Note furthermore 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 4, possibly bringing the dipolar component of the magnetic field in the middle of magnetar range 10 14 − 10 15 G. To compare our results directly with observations, the 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. (2016Makishima et al. ( , 2019 have detected phase modulations in X-ray emissions from two magnetars, that 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(Rea et al. , 2013(Rea et al. , 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 3 , the magnetic field should be scaled down by a factor 30. The final neutron star would have a magnetic dipole of 1 − 2 × 10 13 G and a total magnetic field of 3 × 10 14 , 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 3 The Ekman number used in our simulations would then correspond to the realistic value of the viscosity ν = 2 × 10 10 cm 2 /s. 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 16 G (Siegel et al. 2013;Kiuchi et al. 2018). Our study supports the ability of the MRI to generate a magnetar-like dipole in binary neutron star mergers with a long-lived remnant. A detailed study with an appropriate thermodynamical background state and a 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 et al. 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).