# MRI-driven α – Ω dynamos in protoneutron stars

Alexis Reboul-Salze, Jérôme Guilet, Raphaël Raynaud, Matteo Bugli

► To cite this version:

Alexis Reboul-Salze, Jérôme Guilet, Raphaël Raynaud, Matteo Bugli. MRI-driven α – Ω dynamos in protoneutron stars. Astronomy and Astrophysics - A&A, 2022, 667, pp.A94. 10.1051/0004-6361/202142368. hal-03433445

HAL Id: hal-03433445

https://hal.science/hal-03433445

Submitted on 23 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MRI-driven $\alpha\Omega$ dynamos in protoneutron stars

A. Reboul-Salze¹,², J. Guilet¹, R. Raynaud³, and M. Bugli¹

¹ Laboratoire AIM, CEA/DRF-CNRS-Université Paris Cité, IRFU/Département d'Astrophysique, CEA-Saclay 91191, France
e-mail: alexis.reboul-salze@aei.mpg.de

² Max Planck Institute for Gravitational Physics (Albert Einstein Institute), 14476 Potsdam, Germany

³ Université Paris Cité, Université Paris-Saclay, CNRS, CEA, Astrophysique, Instrumentation et Modélisation, 91191 Gif-sur-Yvette, France

Received 4 October 2021 / Accepted 12 July 2022

## ABSTRACT

**Context.** Magnetars are highly magnetized neutron stars that can produce a wide diversity of X-ray and soft gamma-ray emissions that are powered by magnetic dissipation. Their magnetic dipole is constrained in the range of $10^{14}-10^{15}$ G by the measurement of their spin-down. In addition to fast rotation, these strong fields are also invoked to explain extreme stellar explosions, such as supernovae, which are associated with long gamma-ray bursts and superluminous supernovae. A promising mechanism for explaining magnetar formation is the amplification of the magnetic field by the magnetorotational instability (MRI) in fast-rotating protoneutron stars (PNS). This scenario is supported by recent global incompressible models, which showed that a dipole field with magnetar-like intensity can be generated from small-scale turbulence. However, the impact of important physical ingredients, such as buoyancy and density stratification, on the efficiency of the MRI in generating a dipole field is still unknown.

**Aims.** We assess the impact of the density and entropy profiles on the MRI dynamo in a global model of a fast-rotating PNS. The model focuses on the outer stratified region of the PNS that is stable to convection.

**Methods.** Using the pseudo-spectral code MagIC, we performed 3D Boussinesq and anelastic magnetohydrodynamics simulations in spherical geometry with explicit diffusivities and with differential rotation forced at the outer boundary. The thermodynamic background of the anelastic models was retrieved from the data of 1D core-collapse supernova simulations from the Garching group. We performed a parameter study in which we investigated the influence of different approximations and the effect of the thermal diffusion through the Prandtl number.

**Results.** We obtain a self-sustained turbulent MRI-driven dynamo. This confirms most of our previous incompressible results when they are rescaled for density. The MRI generates a strong turbulent magnetic field and a nondominant equatorial dipole, which represents about 4.3% of the averaged magnetic field strength. Interestingly, an axisymmetric magnetic field at large scales is observed to oscillate with time, which can be described as a mean-field $\alpha\Omega$ dynamo. By comparing these results with models without buoyancy or density stratification, we find that the key ingredient explaining the appearance of this mean-field behavior is the density gradient. Buoyancy due to the entropy gradient damps turbulence in the equatorial plane, but it has a relatively weak influence in the low Prandtl number regime overall, as expected from neutrino diffusion. However, the buoyancy starts to strongly impact the MRI dynamo for Prandtl numbers close to unity.

**Conclusions.** Our results support the hypothesis that the MRI is able to generate magnetar-like large-scale magnetic fields. The results furthermore predict the presence of a $\alpha\Omega$ dynamo in the protoneutron star, which could be important to model in-situ magnetic field amplification in global models of core-collapse supernovae or binary neutron star mergers.

**Key words.** stars: magnetars – supernovae: general – gamma-ray burst: general – dynamo – magnetohydrodynamics (MHD) – methods: numerical

## 1. Introduction

Magnetars are a special class of neutron stars that corresponds to isolated young neutron stars that are characterized by their variable high-energy emission, which is powered by the dissipation of enormous internal magnetic fields (Kouveliotou et al. 1998; Turolla et al. 2015; Kaspi & Beloborodov 2017; Esposito et al. 2021, and references therein). Observations of their pulsed X-ray activity constrain their rotation period $P$ and spin-down $\dot{P}$, two important timing parameters from which a magnetic field

$$B_{\text{dip}} = 10^{14} \left(\frac{P}{5 \text{ s}}\right)^{1/2} \left(\frac{\dot{P}}{10^{-11} \text{ s s}^{-1}}\right)^{1/2} \text{ G} \quad (1)$$

can be inferred under the assumption that only the dipolar field brakes the neutron star. Their activity can be explained by the decay of their ultrastrong magnetic field and includes short bursts (Götz et al. 2006), large outbursts (Coti Zelati et al.

2018, 2021), and giant flares (Hurley et al. 2005; Svinkin et al. 2021) with quasi-periodic oscillations in their signal (Israel et al. 2005; Roberts et al. 2021). In addition, absorption lines that were interpreted as proton cyclotron lines have been detected in two outbursts (Tiengo et al. 2013; Rodríguez Castillo et al. 2016), suggesting a nondipolar surface field stronger than the dipolar component derived from Eq. (1). Another class of transients that can be modeled using magnetars is represented by fast radio bursts (FRBs), that is, millisecond-duration bursts of radio waves discovered in 2007 (Lorimer et al. 2007) that are observed from all directions in the sky. While the exact origin of these events is not yet well known, the observation of the signal FRB 200428 from the Galactic magnetar SGR 1935+2154 (Bochenek et al. 2020; CHIME/FRB Collaboration 2020) confirms that the magnetar model may be able to explain FRBs.

The formation of magnetars in the presence of fast rotation might also be important to explain the most extreme supernovae:
hypernovae (also called broadline type Ic supernovae) associ-
ated with long gamma-ray bursts (Drout et al. 2011) and super-
luminous supernovae (Nicholl et al. 2013; Inserra et al. 2013;
Margalit et al. 2018). These two classes of outstanding stellar
explosions are rare events, representing 0.1% and 1% of all
supernovae, respectively. The first are characterized by a kinetic
energy that is ten times higher than that of standard supernovae,
and the second by a luminosity that is one hundred times higher
than that of standard supernovae. The central engine that can
explain the kinetic energy of hypernovae is based on a strong
large-scale magnetic field that can efficiently extract the high
rotational energy of a fast-rotating protoneutron star (PNS).
This large-scale magnetic field would then launch jets and lead
to a magnetorotational explosion (e.g., 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). The jets that are launched by a millisec-
ond protomagnetar could also lead to (ultra) long gamma-ray
bursts (Duncan & Thompson 1992; Metzger et al. 2011, 2018a).
Millisecond magnetars have also been invoked to power super-
luminous supernovae with their spin-down luminosity, which
corresponds to a delayed energy injection (Nicholl et al. 2013;
Inserra et al. 2013; Margalit et al. 2018). Other high-energy
events could be explained by millisecond magnetars, such as
short gamma-ray bursts (Metzger et al. 2008) and X-ray tran-
sients in the aftermath of binary neutron star mergers (Xue et al.
2019). A millisecond magnetar might indeed form during
neutron star mergers, thus providing an explanation for the
plateau phase and the extended emission in X-ray sources
associated with short gamma-ray bursts (Metzger et al. 2008;
Bucciantini et al. 2012; Rowlinson et al. 2013; Gompertz et al.
2014).

A scenario that might explain the magnetic field of some magnetars is magnetic flux conservation during the collapse of a highly magnetized progenitor. For example, the collapse of progenitors resulting from a stellar merger could lead to the strongest magnetic fields (Schneider et al. 2019). However, this scenario predicts slowly rotating progenitors and therefore can- not explain the formation of millisecond magnetars. To have both fast rotation and a strong magnetic field, an in-situ mag- netic field amplification must be considered. Two mechanisms have been studied: convective dynamos (Thompson & Duncan 1993; Raynaud et al. 2020, 2021; Masada et al. 2022), and the magnetorotational instability (MRI; see Akiyama et al. 2003; Obergaulinger et al. 2009; Reboul-Salze et al. 2021).

The MRI has first been studied in Keplerian accretion disks, both analytically (Balbus & Hawley 1991) and numerically in local shearing box models (Hawley & Balbus 1992). The first studies, which considered the incompressible ideal magneto-hydrodynamics (MHD) framework, showed that in the presence of differential rotation, the turbulent velocity and magnetic field reach a statistically stationary state. Important physical ingredients, such as thermal stratification due to entropy and composition gradients and diffusivities (viscosity, resistivity, and thermal diffusion), were then taken into account by semi-analytical linear analysis (Acheson & Gibbons 1978; Balbus 1995; Menou et al. 2004; Pessah & Chan 2008) and numerical studies (Fromang et al. 2007; Lesur & Longaretti 2007; Simon & Hawley 2009; Longaretti & Lesur 2010; Meheut et al. 2015; Guilet & Müller 2015). An important aspect of disk models that include vertical density stratification is the appearance of oscillating dynamo cycles shaping the structure of the axisymmetric magnetic field (Davis et al. 2010; Shi et al. 2016; Deng et al. 2019).

In the context of core-collapse supernova (CCSN), numeri-
cal studies have shown that the MRI is a viable mechanism to
efficiently amplify the small-scale magnetic field on timescales
shorter than the successful explosion time (typically a few hun-
dred milliseconds; Obergaulinger et al. 2009; Guilet & Müller
2015; Rembiasz et al. 2016). In protoneutron stars, the physi-
cal conditions differ from accretion disks in several ways. The
strong differential rotation inside the PNS (Akiyama et al. 2003;
Ott et al. 2006; Bugli et al. 2020) is non-Keplerian because of
the important role of the pressure gradient in the hydrostatic
force balance (Masada et al. 2012). Due to the high density and
temperature, neutrinos play an important role in the dissipation
of the kinetic energy, which can be modeled as a strong viscos-
ity (Masada et al. 2007; Guilet et al. 2015). Finally, the MRI tur-
bulence can be reduced by the buoyancy forces driven by the
gradients of entropy and lepton fraction in the case of stable
stratification (Guilet & Müller 2015).

The impact of spherical geometry, global entropy, and
density gradients on MRI turbulence 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 only
describe the small-scale turbulence. Mösta et al. (2015) showed
the development of MRI turbulence in the first simulations by
describing a quarter of the PNS whose resolution was high
enough to resolve the MRI wavelength. These simulations
started with a relatively strong initial large-scale magnetic field,
and due to their high computational cost, they could not last
long enough to show dynamo cycles. Reboul-Salze et al. (2021,
hereafter referred to as Paper I) showed that the MRI is able
to robustly generate a large-scale field with magnetar-like inten-
sities under the incompressible approximation, and we extend
these promising results to a more realistic setup including global
gradients here.

This article investigates the impact of a density gradient and
stable stratification on the global properties of the MRI in a 3D
spherical model. The anelastic approximation allows us to take
the effects of density and entropy gradients into account while
filtering sound waves, which drastically increases the time step
and enables long simulation times. 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 nonlin-
ear phase of the MRI and in Sect. 4 for the comparison between
different global models of the MRI. Finally, we discuss the valid-
ity of our assumptions in Sect. 5 and draw our conclusions in
Sect. 6.

## 2. Numerical setup

### 2.1. 1D protoneutron star model

The simulations performed in this article are designed to rep-
resent the outer region of a fast-rotating PNS, which is stable
to the convection and unstable to the MRI. To model this sta-
bly stratified region, we used the same methods and internal
structure as Raynaud et al. (2020), but focused on a different
part of the PNS. The PNS we considered has a baryonic (final
gravitational) mass of 1.78 (1.59) $M_{\odot}$ and was taken from a
1D core-collapse supernova simulation from Hüdepohl (2014).
This simulation used the high-density equation of state LS220
(Lattimer & Swesty 1991) and the nonrotating 27 $M_{\odot}$ progeni-
tor s27.0 by Woosley et al. (2002). The calculations were per-
formed with the code Prometheus-Vertex, which combines the
<figure><img src="image_2.png" /><figcaption>Fig. 1. 1D radial profiles of the entropy per baryon (black line) and density (blue line) of our PNS model at t = 0.2 s post-bounce. Our simulation domain corresponds to the outer stably stratified layer delimited by the thick black line.</figcaption></figure>

hydrodynamics solver Prometheus with the neutrino transport module Vertex (Rampp & Janka 2002). This module solves the energy-dependent moment equations of three species of neutrinos and antineutrinos using a variable Eddington factor closure and including a most recent set of neutrino interaction rates. The energy and lepton number transport by convection is modeled with a mixing length treatment. We chose the physical parameters to represent the PNS at t = 0.2 s after bounce according to the 1D CCSN simulation, whose profiles of density and specific entropy are shown in Fig. 1. While the models presented in Raynaud et al. (2020) focused on the convective zone, our simulation domain spans the outer stably stratified region, thus extending from $r_i = 25.5$ km to $r_o = 39.25$ km, corresponding to the PNS surface defined by the density $\rho_o = 10^{11}$ g cm$^{-3}$. All background profiles of density, temperature, entropy gradient, and gravitational acceleration ($\tilde{\rho}, \tilde{T}, \nabla \tilde{S}$, and $\tilde{g}$) were fit from the PNS model described above with fifth-order to eleventh-order polynomials, hence reproducing the profiles with a good accuracy. Their values at the outer boundary are noted with the letter ‘o’ as a subscript. To have a self-consistent anelastic model, the thermal expansion coefficient at constant pressure $\tilde{\alpha}_T = - (\frac{\partial \rho}{\partial T})_P$ was computed using the following thermodynamic relation:

$$ \frac{d \ln \tilde{T}}{d r} = \frac{1}{c_p} \frac{d}{d r} \tilde{S} - \frac{\tilde{\alpha}_T \tilde{g}}{c_p}, \qquad (2) $$

where $c_p = 3.6 \times 10^8$ erg K$^{-1}$ g$^{-1}$ is the specific heat capacity at constant pressure, assumed to be uniform. The profiles of the 1D CCSN model and the anelastic reference state agree well (see Appendix A). For the sake of simplicity, several assumptions are made in this work. In order to describe the buoyancy associated with both entropy and lepton number gradients, we used an effective entropy gradient assuming the thermal and lepton number diffusivities to be identical. In our model, we also assumed that the neutrinos are in a diffusive regime, so that the effects on the dynamics can be modeled by thermal and viscous diffusivities, $\kappa$ and $\nu$, respectively. The magnetic diffusivity $\eta$ was also explicitly included in our simulations. Finally, we assumed all the diffusivities $\kappa$, $\nu$, and $\eta$ to be constant inside the simulation domain.

## 2.2. Governing equations

We adopted the anelastic approximation to model the flow inside the PNS in order to take the density and effective entropy gradient into account while filtering out sound waves. The MHD anelastic equations describing the dynamics of the PNS in a rotating frame at an angular frequency $\Omega_0 = 958$ rad s$^{-1}$ read (Jones et al. 2011)

$$ \nabla \cdot \tilde{\rho} \mathbf{u} = 0, \qquad (3) $$

$$ \frac{D \mathbf{u}}{Dt} = -\nabla \left(\frac{P}{\tilde{\rho}}\right) - 2\Omega_0 \mathbf{e}_z \times \mathbf{u} - \tilde{\alpha}_T \tilde{T} \frac{s'}{c_p} \tilde{g} \mathbf{e}_r + \frac{1}{\mu_0 \tilde{\rho}} (\nabla \times \mathbf{B}) \times \mathbf{B} + \mathbf{F}^v, \qquad (4) $$

$$ \tilde{\rho} \tilde{T} \left( \frac{D s'}{Dt} + \mathbf{u} \cdot \nabla \tilde{S} \right) = \kappa \nabla \cdot \left( \tilde{\rho} \tilde{T} \nabla s' \right) + \Phi_v + \frac{\eta}{\mu_0} (\nabla \times \mathbf{B})^2, \qquad (5) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \eta \nabla^2 \mathbf{B}, \qquad (6) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (7) $$

where $\mathbf{u}$ is the flow velocity, $\mathbf{B}$ is the magnetic field, $P$ is the pressure, $s'$ is the entropy perturbation, and $\mu_0$ is the vacuum permeability. In this system, viscous force and viscous heating are given by $F_i^v = \tilde{\rho}^{-1} \partial_j \sigma_{ij}$ and $\Phi_v = \partial_j u_i \sigma_{ij}$, where $\sigma_{ij} = 2\tilde{\rho}\nu(e_{ij} - e_{kk}\delta_{ij}/3)$ is the rate of strain tensor and $e_{ij} = (\partial_j u_i + \partial_i u_j)/2$ is the deformation tensor. Tensors are expressed with the Einstein summation convention and the Kronecker symbol $\delta_{ij}$.

## 2.3. Numerical methods

In order to integrate the system of Eqs. (3)–(7) in time, we used the benchmarked pseudo-spectral code MagIC<fnref n="1" /> (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,

$$ \tilde{\rho} \mathbf{u} = \nabla \times \nabla \times (W \mathbf{e}_r) + \nabla \times (Z \mathbf{e}_r), \qquad (8) $$

$$ \mathbf{B} = \nabla \times \nabla \times (b \mathbf{e}_r) + \nabla \times (a_j \mathbf{e}_r), \qquad (9) $$

where $W$ and $Z$ are the poloidal and toroidal kinetic scalar potentials, respectively, while $b$ and $a_j$ are the magnetic potentials. The scalar potentials and the pressure $P$ are decomposed on spherical harmonics for the colatitude $\theta$ and the longitude $\phi$ angles, together with Chebyshev polynomials in the radial direction. 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, we refer to Gilman & Glatzmaier (1981), Tilgner & Busse (1997), and Christensen et al. (2015).

The simulations presented in this paper were performed either using a standard grid resolution of $(n_r, n_\theta, n_\phi) = (257, 512, 1024)$ or a higher resolution of $(n_r, n_\theta, n_\phi) = (385, 768, 1536)$. The resolution was chosen to ensure that the dissipation scales were resolved. For model Standard, about nine cells resolve the resistive scale as the maxima of viscous and resistive dissipation are at the spherical harmonic orders $l_\nu \approx 70$ and $l_\eta \approx 100$, respectively.

<footnote n="1">¹ https://magic-sph.github.io</footnote>
## 2.4. Initial conditions

Many core-collapse simulations with a fast-rotating progenitor have shown that the PNS rotates differentially for several hundred milliseconds (e.g., Akiyama et al. 2003; Ott et al. 2006; Obergaulinger et al. 2018; Bugli et al. 2020). In order to sustain differential rotation in our simulations for a similar duration, we forced the outer boundary to rotate according to the initial rotation profile, with a similar method as we used in our incompressible study (Paper I). The rotation profile was then evolved dynamically inside the simulation domain. The initial cylindrical rotation profile was inspired by the simulations of Bugli et al. (2020) and was composed of an inner part in solid-body rotation and an outer part in differential rotation with a cylindrical profile,

$$ \Omega(s) = \frac{\Omega_i}{\left(1 + \left(\frac{s}{0.25 r_o}\right)^{20 q_o}\right)^{0.05}}, \quad (10) $$

where s is the cylindrical radius, $q_o = 1.5$ corresponds to the shear rate

$$ q = -\frac{s}{\Omega} \frac{d\Omega}{ds} \quad (11) $$

at the outer boundary, and $\Omega_i = 3822 \text{ rad s}^{-1}$ is the rotation rate of the inner core. This rotation rate $\Omega_i$ was computed so that the ratio of total angular momentum over the moment of inertia was equal to the frame rotation rate $\Omega_0$ defined in Sect. 2.2. This was computed with the following formula:

$$ \Omega_i = \frac{I \Omega_0}{\int_V \tilde{\rho} s^2 \left( 1 + \left( \frac{s}{0.25 r_o} \right)^{20 q_o} \right)^{-0.05} dV}, \quad (12) $$

where s is the cylindrical radius, V is the volume of the domain, and I is the moment of inertia of the simulation domain. This initial rotation profile was inspired by those found in the fast-rotating and magnetized core-collapse supernova simulations by Bugli et al. (2020). The shear rate $q_o = 1.5$ might seem high for protoneutron stars, but the shear rate q is lower inside the domain than at the boundary, as already shown in Paper I.

For the initial magnetic field, the toroidal component was set to zero and the poloidal component was initialized with a random superposition of modes with spherical harmonics indices $(l, m)$ with a radial dependence based on Fourier modes as in Paper I, but with a radial modulation of their amplitude to keep a constant Alfvén speed, defined by

$$ v_A = \frac{B}{\sqrt{\mu_0 \tilde{\rho}}}. \quad (13) $$

This initialization implies that the initial magnetic field is stronger in the inner region, which has a higher density (see the left panel of Fig. 2). For all the simulations, the Fourier and spherical harmonic modes were selected so that their wavelength was between [0.3D, 0.5D], with D the shell width,

$$ D = r_o - r_i = 13.7 \text{ km}. \quad (14) $$

The initial root mean-square magnetic field strength at the outer boundary ranged from $B_o = 1.5 \times 10^{14} \text{ G}$ to $B_o = 3.3 \times 10^{14} \text{ G}$ depending on the model. This strong magnetic field allowed us to be sure that the MRI was well resolved: the wavelength of

the fastest-growing mode in ideal MHD with $B_o = 1.5 \times 10^{14} \text{ G}$ reads (Balbus & Hawley 1991)

$$ \lambda_{\text{MRI}} = \frac{8\pi}{q_o(4-q_o)\Omega_0} \frac{B_o}{\sqrt{\mu_0 \tilde{\rho}_o}} \approx 0.69D \approx 9.5 \times 10^5 \text{ cm}. \quad (15) $$

Our strong initial magnetic field may represent the magnetic field generated by the first MRI amplification on small scales found in local models of PNS (Obergaulinger et al. 2009; Guilet & Müller 2015; Rembiasz et al. 2016), as explained in Paper I.

## 2.5. Boundary conditions

We assumed nonpenetrating boundary conditions ($u_r = 0$). At the inner boundary, we used a stress-free condition, where the viscous stress vanished. For the outer boundary, we forced $u_\phi$ to match the initial profile at all times, and the other components of the velocity were set to zero, exactly like in our previous setup (Paper I). For the magnetic field, we used insulating boundary conditions (matching a potential field outside the domain), while the entropy perturbations were set to zero at both boundaries.

## 2.6. Physical parameters and dimensionless numbers

In fast-rotating fluids, the impact of stable stratification on the MRI is characterized by the ratio of the Brunt–Väisälä frequency squared $N^2$ over the rotation frequency squared $\Omega^2$, with

$$ N^2 = - \frac{\tilde{\rho}}{\tilde{\rho}} \left( \left. \frac{\partial \tilde{\rho}}{\partial \tilde{S}} \right|_{\tilde{P}, Y_c} \frac{d\tilde{S}}{dr} + \left. \frac{\partial \tilde{\rho}}{\partial Y_c} \right|_{\tilde{P}, \tilde{S}} \frac{dY_c}{dr} \right), \quad (16) $$

where $Y_c$ is the electron fraction and $\tilde{P}$ is the pressure of the reference state of the PNS model. In our model, this ratio is higher than for the model studied in Guilet & Müller (2015), especially at the equator (see Fig. 3). Our stratification remains small compared to typical values of $N^2/\Omega^2 \sim 10^3-10^4$ in the radiative zone of intermediate and massive stars (Fuller et al. 2019). The Brunt–Väisälä frequency at the outer boundary is $N_o = 4136 \text{ s}^{-1}$. Accordingly, we might expect the buoyancy to dampen MRI-driven turbulence (Guilet & Müller 2015), but the thermal diffusivity $\kappa$ also plays an important role as it reduces the impact of the stable stratification on radial scales smaller than a critical length $L_c$. One way to estimate $L_c$ is to compare the timescale for thermal diffusion to the timescale of gravity waves, from which we obtain

$$ L_c = \sqrt{\frac{\kappa}{N_o}}. \quad (17) $$

For the uniform thermal and viscous diffusivities, we used the values from the 1D CCSN model taken in the middle of the stable zone $r \approx 3.3 \times 10^6 \text{ km}$, which are $\kappa = 1.61 \times 10^{14} \text{ cm s}^{-2}$ and $\nu = 8.03 \times 10^{11} \text{ cm s}^{-2}$, respectively. While all our models share the same constant kinematic viscosity, we considered three different values for the thermal diffusivity, $\kappa = \{1.61 \times 10^{14}, 4.02 \times 10^{13}, 8.03 \times 10^{12}\} \text{ cm s}^{-2}$, which means that the thermal Prandtl number is

$$ \text{Pr} = \frac{\nu}{\kappa} = \{0.005, 0.02, 0.1\}. \quad (18) $$

Our standard value of the thermal Prandtl number $\text{Pr} = 0.005$ corresponds to the value resulting from the 1D core-collapse supernova model (Hüdepohl 2014). With these values, the critical length ranges from $L_c = 0.14L = 1.97 \times 10^5 \text{ cm}$ to
$L_c = 0.031L = 4.4 \times 10^4$ cm. In our simulations, the ratio of the Brunt–Väisälä frequency $N_0$ at the outer boundary to frame rotation $\Omega_0$ was fixed by the following formula:

$$ \frac{N_0}{\Omega_0} = \sqrt{\frac{-RaE^2}{Pr}}, \qquad (19) $$

where the Ekman number $E$ (characterizing the importance of viscosity over the Coriolis force) reads

$$ E = \frac{\nu}{\Omega_0 D^2} = 4.44 \times 10^{-4}, \qquad (20) $$

and the Rayleigh number $Ra$ is defined by

$$ Ra = - \frac{\tilde{\alpha}_{T_o} T_o g_o D^4}{c_p \nu K} \left| \frac{d\tilde{\xi}}{dr} \right|_o. \qquad (21) $$

The Ekman number was kept identical for all the models presented in this article. Therefore, to keep the Brunt–Väisälä frequency constant between our models with buoyancy effects while varying the thermal Prandtl number, the Rayleigh number was also varied from $Ra = -4.73 \times 10^5$ to $Ra = -9.46 \times 10^6$.

With these dimensionless numbers, we can estimate the regime of boundary forcing in our simulations. The viscous timescale $\tau_v = L^2/\nu$ can be compared to the Ekman spin-up timescale $\tau_E = \sqrt{L^2/(\Omega_0\nu)}$ and the Eddington-Sweet circulation timescale $\tau_{ED} = L^2 N_0^2/(\kappa\Omega_0^2)$, (Gouhier et al. 2021, 2022). This gives

$$ \frac{\tau_E}{\tau_v} = \sqrt{E} = 0.021 \ll \frac{\tau_{ED}}{\tau_v} = \text{Pr} \left(\frac{N_0}{\Omega_0}\right)^2 = 0.093 < \frac{\tau_v}{\tau_v} = 1 \qquad (22) $$

for the simulation Standard. This order (i.e., $\tau_E \ll \tau_{ED} \ll \tau_v$) corresponds to the Eddington-Sweet regime, but as Pr increases for some simulations, it approaches the viscous regime (i.e., $\tau_E \ll \tau_v \ll \tau_{ED}$).

We used a constant resistivity of $\eta = 5.0 \times 10^{10}$ cm s$^{-2}$, which gives a magnetic Prandtl number

$$ Pm = \frac{\nu}{\eta} = 16. \qquad (23) $$

As in Paper I, the choice of this value for our simulations was forced by the high cost of simulating high Pm models, and it lies well below a realistic parameter regime. The neutrino viscosity is very high compared to the typical resistivity of a degenerate electron gas inside a PNS. Realistic estimates of the magnetic Prandtl number predict $Pm \approx 10^{13}$ (Thompson & Duncan 1993; Masada et al. 2007), which cannot be reached in direct numerical simulations. The magnetic Reynolds number that characterizes the relative importance of induction to magnetic diffusion is

$$ Rm = \frac{D^2 \Omega_0}{\eta} = 3.6 \times 10^4, \qquad (24) $$

and, for the same reasons, it is quite small compared to realistic estimates. Overall, the only parameters that were varied in this study were the thermal Prandtl number Pr and the Rayleigh number $Ra$, which in two simulations was also set to zero to remove the thermal stratification and allow a clearer assessment of the impact of buoyancy in our models.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Ratio of the squared Brunt–Väisälä frequency $N^2$ (Eq. (16)) to the squared initial rotation frequency $\Omega^2$ (Eq. (10)). The buoyancy influence is strongest in the red regions.</figcaption></figure>

### 3. Typical anelastic simulation

#### 3.1. Quasi-stationary state of an MRI-driven dynamo

We start by describing the magnetic field produced by one fiducial simulation in which we obtain an MRI-driven dynamo, hereafter called model Standard. For this model, we used the standard diffusivities and an initial magnetic field intensity of $B_0 = 1.5 \times 10^{14}$ G (see Table B.1). Figures 3 (right panel) and 4 show selected snapshots of the toroidal magnetic field and the magnetic field lines that highlight the complex geometry of the MRI-driven turbulence. On the meridional cuts in Fig. 3, the spatial distribution of the small-scale turbulence is particularly striking for several reasons. First, almost no turbulence is observed in the equatorial plane, while strong turbulence develops in mid-latitude regions. The same feature can be seen for the kinetic turbulence on the radial velocity $u_r$. This phenomenon can be explained by the influence of the buoyancy force (see Sect. 4.2 for more details). Second, the high-density regions at lower radii are devoid of MRI turbulence because only wide patches of magnetic field can be observed there without small-scale turbulence. This is more unexpected and is discussed in Sect. 5.1. Moreover, the region close to the rotation axis is in solid-body rotation and is therefore stable to the MRI, which leads to no turbulence, as expected. Finally, the comparison between the initial
<figure><img src="image_4.png" /><figcaption>Fig. 3. Meridional slices at $\phi = 0$ of the initial poloidal magnetic field $B_r$ (left) and of the toroidal magnetic field $B_\phi$ at $t = 773$ ms (right) for model Standard.</figcaption></figure>

and saturated magnetic field shown in Fig. 3 suggests that the magnetic field has lost the memory of its initial configuration.

The contrast between the equatorial plane and mid-latitude regions is present in the magnetic field lines as well, where strong and elongated field lines in the azimuthal direction can be seen in the mid-latitude regions of Fig. 4. These field lines are generated by the expected winding of the magnetic field by the shear. Some small-scale loops of field lines can also be seen, with a weaker magnetic field in the mid-latitude regions.

For a more quantitative assessment, we now examine the energetics of this MRI-driven dynamo. Figure 5 shows the time evolution of the toroidal and poloidal magnetic energy densities and the turbulent kinetic energy density. The latter was computed by subtracting the contribution of the axisymmetric rotation from the averaged kinetic energy density in order to separate the differential rotation and the MRI-driven turbulent flow. After several hundred milliseconds, we obtain a statistically stationary state with a mean magnetic field of $B \approx 1.4 \times 10^{14}$ G in the full simulation volume. If we reduce the volume to take the localization of the turbulence into account, the resulting mean magnetic field is $B \approx 2.25 \times 10^{14}$ G.

In order to compare the intensity in the full volume to our previous incompressible study (Paper I), where the density $\rho_0 = 4 \times 10^{13}$ g cm$^{-3}$ was higher, we used the dimensionless magnetic

field strength defined by the Lorentz number,

$$ \mathcal{B} = \frac{B}{\sqrt{\bar{\rho}_0 \mu_0 D \Omega}} \qquad (25) $$

We have $\mathcal{B} \approx 0.096$ for model Standard and $\mathcal{B} \approx 0.064$ for Paper I, which is lower but on the same order of magnitude. This difference is discussed in further details in Sect. 4.4. Compared to Paper I, most of the ratios of the energies or magnetic fields are similar. The kinetic energy is about ten times lower than the total magnetic energy, and the toroidal magnetic field is about ~2.5 times larger than the poloidal magnetic field. The main difference with Paper I is the increased axisymmetric component of the toroidal magnetic field, which is higher than the total poloidal magnetic field and is almost equal to the nonaxisymmetric toroidal magnetic field, as we discuss in Sect. 4.1. An interesting new feature compared to Paper I are synchronous oscillations of the magnetic and kinetic energy densities. These oscillations suggest a dynamo cycle, as we show in Sect. 3.2.

To understand how the magnetic and kinetic energies are distributed over different scales, instantaneous toroidal and poloidal spectra are presented in Fig. 6. On the small scales, the nonaxisymmetric poloidal and toroidal magnetic spectra are similar to those of Paper I. For the poloidal component, the nonaxisymmetric contribution dominates at all scales, as in Paper I. By
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** 3D rendering of the magnetic field lines in simulation Standard at *t* = 773 ms. The color represents the magnetic field amplitude in Gauss.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Time evolution of the magnetic and turbulent kinetic energy densities for model Standard with $B_0 = 1.5 \times 10^{14}$ G. The black and blue lines are the toroidal and poloidal contributions of the magnetic energy density, and the red line is the turbulent kinetic energy density (i.e., the energy of the nonaxisymmetric component of the velocity). The magenta line is the axisymmetric contribution to the toroidal magnetic energy density.</figcaption></figure>

contrast, at large scales, the toroidal magnetic spectrum is dominated by the axisymmetric component, with a particularly strong quadrupole. These strong axisymmetric modes are a new feature of the toroidal spectrum and are linked to the stronger axisymmetric toroidal energy and the mean-field dynamo described in the next section. Overall, these spectra show that the total magnetic energy comes essentially from the axisymmetric contribution at large scales and the nonaxisymmetric contribution at small scales.

Axisymmetric structures that are symmetric with respect to the equatorial plane, such as the rotation and meridional circulation, translate into oscillations between odd and even harmonic order *l* visible in the kinetic spectra (Gubbins & Zhang

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Instantaneous volume-averaged spectra of magnetic and kinetic energies at *t* = 773 ms. **Top:** spectrum of the poloidal magnetic energy (red) and the toroidal magnetic energy (blue) normalized by the total magnetic energy as a function of the spherical harmonics degree *l*. The dotted (solid) lines correspond to the axisymmetric (nonaxisymmetric) contributions to these energies. **Bottom:** spectrum of the poloidal kinetic energy (red) and the toroidal kinetic energy (blue) normalized by the total kinetic energy as a function of the spherical harmonics order *l*.</figcaption></figure>

1993). The kinetic spectra are comparable to those of Paper I because small and intermediate scales are dominated by turbulence, while differential rotation dominates at large scales.

We now highlight the time evolution of the magnetic dipole because it is the component inferred by observations. The dipole energy translates into an intensity $B_{dip} \approx 6.3 \times 10^{12}$ G, which represents $\approx 4.5\%$ of the total magnetic field (Fig. 7). This intensity may seem weak, but the dimensionless magnetic field strength $B_{dip}$ corresponding to this dipole intensity gives $B_{dip} = 0.0043$, which is higher than the value $B_{dip} \approx 0.0032$ obtained for Paper I. Furthermore, the evolution of the dipole energy suggests that the dipole is tilted toward the equator because the axial dipole is twice lower than the average dipole. With the same method as Paper I, an averaged tilt angle of $\theta_{dip} \approx 100^\circ$ can be computed from the magnetic dipole moment. Overall, this model produces results that are qualitatively consistent with those in Paper I, with a stronger dimensionless dipole.

To study the dipole geometry in more detail, we show in Fig. 8 the time evolution and radial dependence of the modes ($l = 1, m = 0$) and ($l = 1, m = 1$). The dipole is mainly present in the low-density region (outer radii), where the turbulence is the strongest. Both modes are radially coherent and form large-scale structures. The equatorial dipole rotates at a different rotation speed than the simulation frame, which leads to
<figure><img src="image_8.png" /><figcaption>Fig. 7. Time evolution of the axial dipole (dashed orange), total dipole (dashed blue), and total magnetic (solid black) energy densities of the model Standard.</figcaption></figure>

oscillations of the equatorial dipole. The timescale of these oscillations (about 50 ms) is much longer than the rotation period of the frame, which means that the magnetic dipole rotates with a frequency close to that of the simulation frame. The time evolution of the axisymmetric component of magnetic dipole is characterized by slower periodic reversals, which again suggests that a mean-field oscillatory behavior is present in the simulations.

## 3.2. $\alpha\Omega$ dynamo

The previous subsection indicated a mean-field dynamo by pointing out an oscillatory behavior. The general mean-field theory has been developed by Moffatt (1978) and Krause & Raedler (1980) and has been widely used to study dynamos. We used the mean-field concept to understand which processes dominate the generation of the mean magnetic field. The basic idea of a mean-field dynamo is that a large-scale magnetic field is generated by small-scale turbulence. The velocity and magnetic fields are therefore decomposed into a mean and a small-scale component, which we represent using the following notation: $X = \overline{X}^{\phi} + X'$. The definition of mean here is the axisymmetric average operator noted $\overline{\phi}^{\phi}$, which verifies the Reynolds averaging rules. The approach of the mean-field theory is to expand the electromagnetic force (EMF) $\mathcal{E}$ only in terms of the mean quantities ($\overline{u}^{\phi}$ and $\overline{B}^{\phi}$) and the statistical properties of the fluctuating quantities ($u'$ and $B'$). We now focus in more detail on the characterization of this mechanism. The most common realization of a mean-field dynamo with differential rotation is the so-called $\alpha\Omega$ dynamo. The $\Omega$ effect corresponds to the shearing of the magnetic field by differential rotation that generates a toroidal magnetic field from a poloidal field. With our cylindrical differential rotation, the $\Omega$ effect reads

$$ \frac{\partial \overline{B}_{\phi}^{\phi}}{\partial t} = s \overline{B}_{s}^{\phi} \frac{d\Omega}{ds}, \qquad (26) $$

and it should induce an anticorrelation between the radial field $\overline{B}_{\phi}^{\phi}$ in cylindrical coordinates and the azimuthal magnetic field $\overline{B}_{\phi}^{\phi}$.

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Radius-time diagram of the dipolar modes. Top: radius-time diagram of the axial dipolar mode $B_r(l = 1, m = 0)$. The period of the dipole reversal is $P_{dip} \approx 410$ ms. Bottom: radius-time diagram of the equatorial dipolar mode amplitude $B_{\phi}(l = 1, m = 1)$ at a given azimuth $\phi = 0$ in the rotating frame of the simulation. The fast oscillations are due to the rotation of this dipole mode.</figcaption></figure>

The $\alpha$ effect comes instead from the closure relation of the mean EMF,

$$ \mathcal{E} = \overline{u' \times B'}^{\phi}, \qquad (27) $$

where $u' = u - \overline{u}^{\phi}$ and $B' = B - \overline{B}^{\phi}$, expressed as a function of the mean magnetic field $\overline{B}^{\phi}$,

$$ \mathcal{E}_i = \alpha_{ij} \overline{B}_j^{\phi} + \beta_{ij} \left( \overline{\nabla \times B}^{\phi} \right)_j, \qquad (28) $$

where $\alpha_{ij}$ and $\beta_{ij}$ are tensors that do not depend on $\overline{B}^{\phi}$, and i, j refer to spherical r, $\theta$, $\phi$ or cylindrical coordinates s, $\phi$, z. The diagonal components of the $\alpha$ tensor correspond to the component of the EMF in the direction of the mean magnetic field, and their effect is physically described as the twisting of the mean magnetic field lines by the cyclonic turbulence, which forms magnetic field loops that can generate poloidal magnetic field from the toroidal magnetic field and vice versa.

In our case, the generation of the poloidal magnetic field by this effect can be seen as a correlation between the toroidal component of the EMF and the toroidal component of the magnetic
field $B_\phi$ in the form of $\mathcal{E}_\phi = \alpha_{\phi\phi} \overline{B_\phi}^\phi$. The diagonal components of the $\beta$ tensor are in the direction of the mean current $\overline{\mathcal{J}}^\phi = \mu_0^{-1} \nabla \times \overline{\mathcal{B}}^\phi$, and their effect is physically described as a turbulent diffusivity, which adds to the magnetic diffusivity $\eta$. Another effect that has been proposed to complete the dynamo loop in the case of the MRI is the nondiagonal resistivity $\beta_{\phi s}$, which could generate a poloidal field from the toroidal field (Lesur & Ogilvie 2008). In this case, the azimuthal component of the electromotive force $\mathcal{E}_\phi$ would be correlated with the radial current $\overline{\mathcal{J}}_s^\phi$ in cylindrical coordinates. One (or both) of these two effects, in combination with the $\Omega$ effect, may therefore be expected to complete the dynamo cycle, leading to a self-sustained magnetic field.

The usual way to characterize an $\alpha\Omega$ dynamo is to compute space-time diagrams of $B_\phi$, $B_s$, and $\mathcal{E}_\phi$, which are often referred to as butterfly diagrams. Figure 9 shows these azimuthally averaged quantities in the southern hemisphere at $r = 0.86 r_o$ (which lies in the middle of the turbulent region). First, the butterfly diagrams show large coherent structures in the mid-latitudes, which is consistent with a mean-field oscillatory dynamo with a period of $P \approx 410$ ms. This period is estimated by taking the peaks of the mean toroidal magnetic field at different colatitudes and averaging the frequencies for each hemisphere. The visual inspection of the butterfly diagrams suggests that the dynamo can be interpreted as an $\alpha\Omega$ dynamo. The components $\overline{B}_\phi^\phi$ and $\overline{B}_s^\phi$ are anticorrelated, which shows that $\overline{B}_\phi^\phi$ is mainly generated by the $\Omega$ effect. On the other hand, the electromotive force $\mathcal{E}_\phi$ is correlated with $\overline{B}_\phi^\phi$, which suggests that the $\alpha$ effect plays an important role.

To corroborate the visual correlations of the butterfly diagrams, we computed the Pearson correlation coefficient between two quantities $X$ and $Y$ with the following formula:

$$C_P(X, Y) = \frac{\int_t (X - \langle X \rangle_t) dt \int_t (Y - \langle Y \rangle_t) dt}{\sqrt{\left(\int_t (X - \langle X \rangle_t)^2 dt\right) \left(\int_t (Y - \langle Y \rangle_t)^2 dt\right)}}, \quad (29)$$

where $\langle \cdot \rangle_t$ represents a time average. Figure 10 shows the correlation coefficients between $\mathcal{E}_\phi$ and $\overline{B}_\phi^\phi$ and the radial current $\overline{\mathcal{J}}_s^\phi$ in order to test the nondiagonal resistivity hypothesis. We find that $\overline{B}_\phi^\phi$ and $\mathcal{E}_\phi$ are well correlated, while no correlation is found between $\overline{\mathcal{J}}_s$ and $\mathcal{E}_\phi$. This means that the $\alpha$ effect is prominent in our simulations. Moreover, the antisymmetry of the correlation matches the expected symmetry of the components of the tensor $\alpha$.

We can estimate the value of the diagonal components of the $\alpha$ tensor with the formula:

$$\alpha_{ii} = \frac{\langle \overline{B}_i^\phi \overline{\mathcal{E}}_i^\phi \rangle_t}{\langle (\overline{B}_i^\phi)^2 \rangle_t}. \quad (30)$$

This estimation assumes that the EMF is only due to the $\alpha$ effect, which is a good approximation in the case of high correlation values. A theoretical estimation is possible under the second-order correlation approximation, which considers only second-order fluctuating quantities, in addition to several hypotheses (Moffatt 1978; Krause & Raedler 1980). The turbulence is, in fact, assumed to be statistically homogeneous and isotropic, and the mean flow $\overline{u}^\phi$ is usually also neglected. The second-order correlation approximation is valid when one of two dimensionless numbers are small: the magnetic Reynolds number $Rm$, or the Strouhal number $St = V\tau/L$, where $V$, $\tau$ and $L$ are typical values of the velocity, time variation, and length scale of the turbulence. Under these hypotheses, the diagonal components of the $\alpha$ tensor are proportional to the kinetic helicity $h$,

<figure><img src="image_11.png" /><figcaption>Fig. 9. Butterfly diagrams of $\overline{B}_\phi^\phi$ (top), $\overline{B}_s^\phi$ (middle), and $\mathcal{E}_\phi$ (bottom) in the southern hemisphere of model Standard. Each butterfly diagram is computed at the same spherical radius $r = 0.86 r_o$ in the turbulent region.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Time-averaged correlation coefficients (Eq. (29)) between $\mathcal{E}_\phi$ and $\overline{B}_\phi^\phi$ (blue) and between $\mathcal{E}_\phi$ and $\overline{\mathcal{J}}_s^\phi$ (red) taken at $r = 0.86 r_o$ in model Standard.</figcaption></figure>

$$\alpha_{ii} = -\frac{\tau_c}{3} \overline{u' \cdot \nabla \times u'}^\phi = -\frac{\tau_c}{3} h, \quad (31)$$
<figure><img src="image_13.png" /><figcaption>Fig. 11. Time-averaged values of the $\alpha_{\phi\phi}$ component estimated with Eq. (30), (blue) and with the turbulent kinetic helicity (orange), (Eq. (31)) taken at $r = 0.86 r_o$ and with $\tau_c = 2.5$ ms.</figcaption></figure>

where $\tau_c$ is the correlation timescale of the turbulence. Although the conditions that are theoretically necessary for the validity of this formula are not satisfied in our simulation<fnref n="2" />, we find a quantitative agreement between Eqs. (30) and (31) for $\tau_c \approx 2.5$ ms = $0.38 \times 2\pi/\Omega_0$ with a peak value of $\alpha_{\phi\phi} \approx 6 \times 10^5$ cm s⁻¹ (see Fig. 11). This empirical turbulent correlation time is consistent with the timescale of MRI turbulence, which is typically a fraction of the rotation period.

Finally, the value of $\alpha_{\phi\phi}$ can be used to estimate the theoretical frequency of an $\alpha\Omega$ dynamo using the relation (e.g., Busse & Simitev 2006; Gastine et al. 2012; Gressel & Pessah 2015)

$$ \omega_{\alpha\Omega} = \left| \frac{1}{2} \alpha_{\phi\phi} \frac{d\Omega}{d\ln s} k_z \right|^{1/2} = \left| \frac{1}{2} \alpha_{\phi\phi} q \Omega k_z \right|^{1/2}, \quad (32) $$

where $k_z$ is the vertical wavenumber. We computed the shear rate $q = -0.65$ and the rotation rate $\Omega = 688$ rad s⁻¹ in the middle of the turbulent region from a vertically and azimuthally averaged rotation profile, which gives lower values than $q_o$ and $\Omega_0$. We estimated $k_z \approx 2.79 \times 10^{-06}$ cm⁻¹ by taking the longest vertical length in the turbulent region at mid-latitudes and obtained a period of $P_{\alpha\Omega} = 2\pi/\omega_{\alpha\Omega} = 324$ ms, which is roughly agrees with the period inferred from the butterfly diagrams. The small difference between the two values is compatible with the uncertainties in our estimate of the different parameters and may, for example, be due to an overestimation of the $\alpha$ effect, as we take its peak value in the turbulent region. If we consider an average on the angles $\theta \in [30^\circ, 60^\circ]$, we obtain a period $P_{\alpha\Omega} = 393$ ms, which is closer to the measured period $P \approx 410$ ms. Moreover, according to the theoretical Parker-Yoshimura rule (Parker 1955; Yoshimura 1975), the dynamo wave should propagate in the direction of $-\alpha_{\phi\phi}e_\phi \times \nabla\Omega$, that is, from the equator toward the pole in our case. Unfortunately, a difficulty arises here because the full pattern does not have a clear propagation direction, as shown in Fig. 9. Overall, all these results show that the MRI-driven dynamo we obtained in our simulations can be described as an $\alpha\Omega$ dynamo.

## 4. Comparison with other models

After demonstrating with a realistic PNS setup that the MRI can produce a subdominant dipole and a turbulent dynamo with a mean-field behavior, we aim at understanding the influence of the different physical ingredients we added in the anelastic model (shell aspect ratio, buoyancy, thermodynamic background, etc.). We discuss in particular the impact of the density profile (Sect. 4.1), the entropy profile (Sect. 4.2), and the thermal diffusivity (Sect. 4.3). We compare these results with our previous incompressible study (Sect. 4.4).

### 4.1. Impact of the density gradient

We compared the results of four simulations under different approximations: the Boussinesq approximation (i.e., no density gradient) or the anelastic approximation (i.e., with a density gradient), both with or without buoyancy ($N^2 > 0$ or $N^2 = 0$). For Boussinesq models, constant density was taken to be equal to the density $\rho_o = 10^{11}$ g cm⁻³ at the outer boundary of the anelastic model. All the other simulation parameters were kept identical. In order to simplify the interpretation of the results, we start by describing the effect of the density gradient.

The comparison of the snapshots of $B_\phi$ in Fig. 12 gives qualitative insights into the effect of the density gradient. While the maximum intensity of the magnetic field is slightly higher for Boussinesq simulations, the striking difference between Boussinesq and anelastic simulations is their structure: for the former models, MRI-driven turbulence develops through the domain, except at the poles. It is restricted to the outer low-density layers for the latter. It also seems that at mid-latitudes, the magnetic field has structures at slightly larger scales in the anelastic cases. These differences are due to the background density gradient because both anelastic simulations (with or without buoyancy) have a similar magnetic field structure.

The presence of the background density gradient also impacts the angular velocity. The bottom panels of Fig. 12 show that the flow inside the domain rotates more slowly with a density gradient. This can be understood as a consequence of the lower efficiency of the inward transport of the angular moment by the outer boundary forcing in the presence of a density gradient due to the low density at the outer boundary.

Figure 13 compares the time evolution of the turbulent magnetic (green) and kinetic (orange) energy densities and the axisymmetric magnetic energy density (purple). Both Boussinesq simulations have a higher turbulent magnetic energy and kinetic energy. This is at least partly due to the presence of turbulence in only half of the domain for anelastic simulations compared to the full domain for Boussinesq simulations because the maximum field strength difference in the snapshots is too low to explain the difference in magnetic energy. On the other hand, the Boussinesq model with buoyancy has a similar axisymmetric magnetic energy density such that the ratio of the axisymmetric magnetic energy and turbulent magnetic energy is higher for anelastic simulations. The axisymmetric component for the Boussinesq model without buoyancy, called model Incompressible, is quite peculiar because it increases faster than the turbulent magnetic energy and reaches a higher energy than all other models. The ratio of axisymmetric magnetic energy to turbulent magnetic energy also saturates at higher levels. Another important point, however, is that clear oscillations occur only for anelastic simulations, which suggests that there is no mean-field dynamo with a constant density.

In a similar fashion to axisymmetric magnetic energy, the dipole energy density seems to be on the same order for all

<footnote n="2">² The Strouhal number is on the order of unity and the magnetic Reynolds number is on the order of a few hundred.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Meridional cuts of $B_\phi$ (top) and azimuthal average of $\Omega$ (bottom) at a statistically stationary state for incompressible (panels a and e), Boussinesq (panels b and f), anelastic without buoyancy (panels c and g), and anelastic including buoyancy (panels d and h) models.</figcaption></figure>

simulations, except for model Incompressible, which is higher at late times (Fig. 14). Because their total magnetic energy also increases, both Boussinesq models have a similar ratio of dipole field to total magnetic field of ≈3.4%, which is lower than the ratios of the anelastic models ≈4.3%, as we discuss in Sect. 4.4. The dipole energy density is also dominated by its equatorial component in all simulations, which indicates that the inclination of the dipole toward the equator is a robust feature of the MRI.

To confirm the mean-field behavior, we show in Fig. 15 the butterfly diagrams of $\vec{B}_\phi^\phi$. For Boussinesq simulations (panels a and b), the MRI-unstable region is noisier. It is hard to observe clear signs of coherent mean-field patterns with buoyancy (panel a), while without buoyancy (panel b), an increasing quadrupole develops after 400 ms. This quadrupole is rather puzzling because its strength is higher than the incompressible models of Paper I, and it does not match the $\alpha\Omega$ behavior. There is no clear signal in the toroidal component of the EMF $\mathcal{E}_\phi$, and we find that the correlation coefficients between the EMF $\mathcal{E}_\phi$ and $\vec{B}_\phi^\phi$ are low, with an average $C_P(\mathcal{E}_\phi, \vec{B}_\phi^\phi) \approx 0.04$ in both hemispheres.

We find no correlations either between the EMF $\mathcal{E}_\phi$ and $\vec{J}_s^\phi$. The physical origin of this toroidal quadrupole is so far uncertain. We would like to point out that this region is not turbulent, and the patterns found there (which can also be seen in the Boussinesq model with buoyancy in panel a) are probably not caused by a true mean-field mechanism. These results suggest that it is more difficult to have a mean-field behavior without density stratification. The stronger ratio of dipole to total magnetic field in the case of anelastic models might be due to the mean-field dynamo boosting the generation of the dipole.

## 4.2. Impact of buoyancy

To study the influence of buoyancy, we focused first on the comparison of model Incompressible to the Boussinesq model with buoyancy. The comparison between the snapshots of $B_\phi$ shows that the structure of the turbulence is rather similar, while the turbulent magnetic field is stronger without buoyancy (Fig. 12). This result is also found in terms of turbulent energy density for the magnetic and kinetic energies (Fig. 13). It can be expected that the higher magnetic field increases the
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Time evolution of the turbulent kinetic (orange), nonaxisymmetric (green), and axisymmetric toroidal (purple) magnetic energy densities for incompressible (dotted line), Boussinesq (dash-dotted line), anelastic without buoyancy (dashed line), and anelastic (solid line) models.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Time evolution of the axial dipole (green), total dipole (blue), and total magnetic (black) energy densities of the same models as displayed in Fig. 13.</figcaption></figure>

angular momentum transport, thus decreasing the angular frequency near the axis. As mentioned in the previous section, the physical origin of the high axisymmetric energy density of model Incompressible is uncertain, and this feature is not discussed.

For anelastic models, we see a difference in the equatorial plane where MRI-driven turbulence is damped with buoyancy (panels c and d of Fig. 12). Moreover, the butterfly diagrams in Fig. 15 suggest that the model with buoyancy (panel c) persistently suppresses small-scale fluctuations in the equatorial region, while without buoyancy (panel d), turbulence is fully developed. This effect may be expected because in the equatorial plane, the buoyancy quenches the motions in the direction of the differential rotation gradient, which reduces the MRI turbulence (Menou et al. 2004; Guilet & Müller 2015). The weaker turbulence in the equatorial plane with buoyancy may partly explain the lower turbulent magnetic and kinetic energy densities of model Standard compared to the model without buoyancy (Fig. 13).

By contrast, the mid-latitude regions are less impacted by the buoyancy. The magnetic field $B_\phi$ is similar in strength and structure outside of the equatorial plane (Fig. 12). In the anelastic model without buoyancy, a mean-field dynamo also produces a field with a similar amplitude to model Standard at the midlatitudes (see panel c and d in Fig. 15). This leads to a similar axisymmetric toroidal energy density (Fig. 13) and dipole energy density (Fig. 14) for both anelastic models. In terms of dynamo periods, both theoretical dynamo periods are similar, even though the patterns have a longer cycle without buoyancy. The dynamo frequencies agree within 15% with buoyancy and 24% without buoyancy (see Table B.1). For these two models, the dynamo periods are more difficult to measure due to the combination of short and long patches of magnetic fields, and this can lead to differences between the hemispheres (see Table 1 for model standard). Considering the different uncertainties in our estimate, the two models give consistent results with theoretical values, which shows that buoyancy has a weak influence on the dynamo mechanism at $Pr = 0.005$.

### 4.3. Influence of the thermal Prandtl number

The previous section shows that except in the equatorial plane, buoyancy has a rather small impact on the anelastic simulations overall, which might be due to the high thermal diffusivity. To further study the influence of the thermal diffusion, we ran two anelastic simulations with a lower thermal diffusivity corresponding to larger Prandtl numbers, $Pr = 0.02$ and $Pr = 0.1$, respectively. By comparing the timescale for thermal diffusion to compensate for the entropy fluctuations to the timescale of gravity waves, we may expect the thermal diffusion to reduce the effects of buoyancy on scales smaller than the critical length $L_c$ defined by Eq. (17). For our simulation Standard at $Pr = 0.005$, the critical length at mid-latitudes $L_c \approx 0.14D$ is on the same order of magnitude as the turbulent scale of the radial velocity, such that buoyancy effects are expected to be marginal. By contrast, the critical length decreases to $L_c \approx 0.03D$ at $Pr = 0.1$, hence decreasing the range of scales where turbulent motions are suppressed by buoyancy (i.e., for scales larger than $L_c$). As a consequence, the typical scale of the radial velocity is expected to be smaller at higher $Pr$.

To verify this expectation, we compare snapshots of the radial velocity in Fig. 16. For increasing $Pr$, we observe a decrease in the size of the velocity structures and in the maximum amplitude of the radial velocity.

This trend is also confirmed in the radial velocity spectrum of the different simulations (Fig. 17). The peak values of the spectrum in the small scales for $l > 20$ are lower and shifted toward smaller scales with increasing $Pr$. For models with $Pr = 0.02$ and $Pr = 0.1$, they also seem to match the theoretical critical degree $l_c \approx \frac{R}{L_c}$ well, corresponding to the critical length.

The suppression of small scale-turbulence can also be seen to some extent in the corresponding kinetic energy time series in Fig. 18. The kinetic energy of the simulation at $Pr = 0.1$ (dotted lines) is indeed the lowest during most of the simulation time. The turbulent magnetic energy follows the same trend with $Pr$ as the turbulent kinetic energy (green curves).

Despite the differences in the turbulence properties, the axisymmetric toroidal magnetic energy is not strongly dependent on $Pr$ and the characteristic oscillations of a mean-field dynamo are present in all three simulations. This suggests that buoyancy does not strongly impact the mean-field behavior at midlatitudes. The main difference is the weaker turbulence in the equatorial plane for model Standard, while no turbulence can
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Butterfly diagrams of $\vec{B}_{\phi}^\phi$ at $r = 0.92 r_0$ for the models shown in Fig. 12 for Boussinesq with buoyancy (panel a), Boussinesq without buoyancy (panel b), anelastic with buoyancy (panel c), anelastic without buoyancy (panel d), and anelastic with buoyancy with $Pr = 0.02$ (panel e) and $Pr = 0.1$ (panel f). Each butterfly diagram is computed at the same spherical radius $r = 0.92 r_0$ (in the turbulent region). (a) $\nabla \cdot u = 0, N^2 > 0$, $Pr = 0.005$. (b) $\nabla \cdot u = 0, N^2 = 0$, $Pr = 0.005$. (c) $\nabla \cdot (\tilde{\rho}u) = 0, N^2 > 0$, $Pr = 0.005$. (d) $\nabla \cdot (\tilde{\rho}u) = 0, N^2 = 0$, $Pr = 0.005$. (e) $\nabla \cdot (\tilde{\rho}u) = 0, N^2 > 0$, $Pr = 0.02$. (f) $\nabla \cdot (\tilde{\rho}u) = 0, N^2 > 0$, $Pr = 0.1$.</figcaption></figure>

be seen for the other models. The butterfly diagrams in the turbulent region at mid-latitudes are similar for $Pr = 0.02$ (panel e of Fig. 15) and the south hemisphere of the Standard simulation at $Pr = 0.005$ (panel c), with a similar dynamo period (see Table 1). One noticeable difference between panels c and e in Fig. 15 is the poleward propagation of the dynamo wave, which is consistent

with the Parker-Yoshimura rule for the simulation at $Pr = 0.02$, probably due to the lack of turbulence in the equatorial plane.

To determine whether the specific behavior at different $Pr$ is still an $\alpha\Omega$ dynamo, we compare in Fig. 19 the $\alpha_{\phi\phi}$ component for the different Prandtl numbers. At $Pr = 0.02$, we find a mean $\alpha_{\phi\phi} \approx 8.0 \times 10^5 \text{ cm s}^{-1}$ in the turbulent region, which translates
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Radial velocity $u_r$ and $\alpha_{\phi\phi}$ for anelastic models with buoyancy and for increasing Prandtl numbers. *Top*: snapshots of the radial velocity $u_r$. From left to right, the slices correspond to Pr = 0.005 at $t = 773$ ms, Pr = 0.02 at $t = 1128$ ms, and Pr = 0.1 at $t = 1128$ ms. *Bottom*: 2D distribution of $\alpha_{\phi\phi}$ for Pr = 0.005 (**left**), Pr = 0.02 (**middle**), and Pr = 0.1 (**right**). The dashed line represents the spherical radius $r = 0.92 r_o$.</figcaption></figure>

into a period $P_{\alpha\Omega} \approx 300$ ms using the same $k_z$ and Eq. (32). This agrees with the measured dynamo frequency within 17%. Moreover, both values are relatively close to those measured and estimated for the southern hemisphere of model Standard.

The 2D meridional distribution of $\alpha_{\phi\phi}$ (bottom left and bottom center panels of Fig. 16) also shows that there is little difference between these two models for the mean-field dynamo. This can be interpreted as follows: most of the kinetic turbulence is still at smaller scales than the critical wavelength for the buoyancy, as seen in the $u_r$ spectra (Fig. 17), and the MRI is not as much impacted by buoyancy at mid-latitudes. Therefore, the

change in the velocity spectrum is not important enough to modify the $\alpha$ effect strongly, and both models have a similar mean-field dynamo.

At the higher Prandtl number Pr = 0.1, oscillating dynamo cycles are still present, but have different characteristics. In the northern hemisphere, the patterns look like those in the other anelastic simulations, albeit with a lower frequency $P_{\text{North}} \approx 493$ ms. However, in the southern hemisphere, we observe mean-field patterns of similar amplitude, but higher frequency (with a period of $P_{\text{South}} \approx 213$ ms). Moreover, the pattern propagates toward the equator, which is in the opposite direction to the
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Normalized volume-averaged spectrum of the radial velocity squared $u_r$ for increasing Prandtl numbers $Pr = 0.005$ (blue), $Pr = 0.02$ (orange), and $Pr = 0.1$ (green). All the spectra are normalized by the integral over $l$ of the Standard model spectrum. The dashed lines correspond to the critical degree $l_c \approx R/L_c$ above which the scales are impacted by buoyancy.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Time-averaged values of the $\alpha_{\phi\phi}$ component estimated using Eq. (30) at $r = 0.92 r_o$ in the turbulent region for three different Prandtl numbers: $Pr = 0.005$ (blue), $Pr = 0.02$ (orange), and $Pr = 0.1$ (green).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Time evolution of the magnetic and turbulent kinetic energy densities for models with $Pr = 0.1$ (dotted), $Pr = 0.02$ (dashed), and $Pr = 0.005$ (solid). The colors represent the same quantities as in Fig. 13.</figcaption></figure>

<table><caption>**Table 1.** Periods of the different mean-field patterns for $Pr = 0.005$, $Pr = 0.02$, and $Pr = 0.1$.</caption><thead><tr><th>Model</th><th>Pr</th><th>P<sub>North</sub><br>[ms]</th><th>P<sub>South</sub><br>[ms]</th><th>P<sub>αΩ</sub><br>[ms]</th></tr></thead><tbody><tr><td>Standard</td><td>0.005</td><td>530</td><td>401</td><td>342</td></tr><tr><td>Anel Pr0 02</td><td>0.02</td><td>357</td><td>366</td><td>300</td></tr><tr><td>Anel Pr0 1</td><td>0.1</td><td>493</td><td>213</td><td>366</td></tr></tbody></table>

Parker-Yoshimura rule. The meridional distribution of $\alpha_{\phi\phi}$ (bottom left panel of Fig. 16) also highlights the smaller turbulent region and the lower amplitude of $\alpha_{\phi\phi}$ than in the other two models.

At $Pr = 0.1$, we measure $\alpha_{\phi\phi} \approx 4.57 \times 10^5 \text{ cm s}^{-1}$ in the northern and $\alpha_{\phi\phi} \approx 4.47 \times 10^5 \text{ cm s}^{-1}$ in the southern hemisphere (Fig. 19). By assuming that the turbulent region is 20% smaller than the region from simulation Standard, we obtain a theoretical prediction of the period $P_{\alpha\Omega} = 366$ ms in the northern

hemisphere, which agrees within 27% with the numerical value. The similar $\alpha_{\phi\phi}$ value in the southern hemisphere would lead to a similar period and therefore disagrees more strongly with the measured numerical value. This result, combined with the opposite direction of propagation, may suggest a significant deviation from the $\alpha\Omega$ dynamo formalism. This new behavior of the southern hemisphere with $Pr = 0.1$ remains obscure and might result from being close to the dynamo threshold. This idea is also supported by the fact that a more typical $\alpha\Omega$ dynamo occurs in the northern hemisphere, but with a weaker $\alpha$ effect. All in all, the results presented in this section show that a low Prandtl number diminishes the impact of buoyancy on the MRI dynamo, whereas the buoyancy force can limit the MRI-driven dynamo for Prandtl numbers closer to unity.

### 4.4. Comparison to our previous incompressible models

One of the important results of Paper I is the robust linear relation between the dipole intensity and the averaged magnetic field strength (see their Fig. 16). In order to directly compare the ratios measured in simulations that have different densities and parameters, we reproduce the same figure by using the dimensionless dipole strength defined by $B_{\text{dip}} = B_{\text{dip}}/(\sqrt{\mu_0\bar{\rho}_o}D\Omega)$ and the dimensionless total magnetic field strength $B_{\text{tot}} = B_{\text{tot}}/(\sqrt{\mu_0\bar{\rho}_o}D\Omega)$ instead of the magnetic intensities. For anelastic simulations, we find a greater magnetic strength $B_{\text{tot}}$ and dipole strength $B_{\text{dip}}$ than in Paper I. However, the ratio of $B_{\text{dip}}$ to $B_{\text{tot}}$ of anelastic simulations is quite similar, if slightly lower ($\approx 4.3\%$) than the linear relation of Paper I.

For our Boussinesq simulations, the differences with Paper I are stronger. The magnetic strength $B_{\text{tot}}$ is higher by a factor of approximately three due to the stronger turbulence, and the dipole strength $B_{\text{dip}}$ is twice stronger. This gives a different linear relation of $\approx 3.4\%$ between the magnetic dipole and the magnetic field strength. This change in linear relation between the Boussinesq models here and the incompressible models of Paper I is probably due to the difference in the aspect ratio between the models. With a smaller shell gap $D$, the forcing of the differential rotation by the outer boundary is expected to be more efficient, which might lead to stronger turbulence for a similar dipole intensity. This result highlights the importance of developing global models that take the full PNS structure into account.
<figure><img src="image_22.png" /><figcaption>Fig. 20. Dimensionless dipole strength $B_{\text{dipole}}$ as a function of the dimensionless magnetic field strength $B_{\text{tot}}$. The simulations of this work are plotted in blue (anelastic models) and green (no density stratification), and the incompressible simulations from Paper I are shown in red. Dimensionless magnetic field strengths $B$ were used in order to compare results with different densities $\bar{\rho}_o$ and parameters.</figcaption></figure>

# 5. Discussion

## 5.1. Spatial distribution of the turbulence

One puzzling feature of the turbulence in the anelastic models is its concentration in the low-density region. The MRI modes can be damped by the viscosity or the resistivity if the initial magnetic field is lower than a critical value. In order to investigate whether the lack of turbulence in the high-density region is due to the initial magnetic field, we ran a simulation with a higher resolution $(n_r, n_\theta, n_\phi) = (385, 768, 1536)$ and a magnetic field that was twice as high $B_o = 3.3 \times 10^{14}$ G. We find similar results for the magnetic field as model Standard as shown in Table B.1, which confirms that our results converge well and are not an artifact of the initial magnetic field. Our interpretation of the absence of turbulence in the high-density region is that the forcing of the differential rotation is not efficient enough to sustain the turbulence in this region. As shown in the bottom panels of Fig. 12, the angular velocity in the domain is slower than at the outer boundary. This result might change with a different forcing method or when the magnetic Prandtl number is increased with a lower physical resistivity, which would make the turbulence easier to sustain.

## 5.2. $\alpha\Omega$ versus $\alpha^2\Omega$ dynamo

As shown in the previous sections, the dynamos observed in our anelastic models are consistent with an $\alpha\Omega$ mechanism. Although the anticorrelation of $B_s$ and $B_\phi$ in the butterfly diagrams (Fig. 9) suggests that the $\Omega$ effect generates the toroidal field, we would like to assess this more quantitatively. For this purpose, we compared the $\Omega$ effect to the generation of toroidal field by the diagonal components $\alpha_{rr}$ and $\alpha_{\theta\theta}$ (see Appendix D), which would be relevant in the case of an $\alpha^2\Omega$ dynamo. To distinguish between an $\alpha\Omega$ and $\alpha^2\Omega$ dynamo, we computed the ratio of the two dynamo numbers $C_\alpha = \max(\alpha_{rr}, \alpha_{\theta\theta})R/\eta$ and $C_\Omega = \Omega R^2/\eta$ in the turbulent region, which gives

$$ \frac{C_{\Omega}}{C_{\alpha}} = \frac{\Omega R}{\alpha} \approx 100, \qquad (33) $$

where $R$ is taken as the middle radius $(r_i + r_o)/2$. This estimate therefore confirms that the $\Omega$ effect largely dominates the generation of the toroidal field and that the dynamo in our simulations is an $\alpha\Omega$ dynamo.

## 5.3. Properties of the $\alpha$ tensor components

The MRI-driven $\alpha\Omega$ dynamo in our simulations was characterized by estimating the tensor components $\alpha_{ij}$. Our method assumes that the EMF component $E_i$ is only due to the contribution of the mean magnetic field component $B_j$ and neglects the contribution of the other components of the $\alpha$ tensor and the turbulent resistivity tensor $\beta$ components. In the azimuthal EMF $E_\phi$, the contribution from $\alpha_{\phi r}$ and $\alpha_{\phi \theta}$ can be neglected because the Pearson correlation coefficients between $E_\phi$ and $B_\theta$ or $B_r$ are much lower in the turbulent region than they are for $B_\phi$ (panels c and d of Appendix C). The $\alpha_{\phi\phi}$ component is well estimated and is therefore the main contribution to the generation of the poloidal magnetic field. Due to their weaker correlation in the turbulent region, the other $\alpha_{\phi i}$ components (Appendix D) may be incorrectly estimated, and they do not drive the generation of the poloidal field. The other $\alpha$ components do not dominate the generation of toroidal field because the dynamo cannot be described as an $\alpha^2\Omega$ dynamo.

In the case of a statistically symmetric velocity field with respect to the equator, mean-field theory predicts that the components of the $\alpha_{ij}$ tensor (Eq. (28)) are either equatorial-symmetric ($\alpha_{rb}$, $\alpha_{\theta r}$, $\alpha_{\phi \theta}$, $\alpha_{\theta \phi}$) or equatorial-antisymmetric ($\alpha_{rr}$, $\alpha_{\phi r}$, $\alpha_{\phi \phi}$, $\alpha_{r \phi}$, $\alpha_{\theta \theta}$), (Gubbins & Zhang 1993). The estimation of the components $\alpha_{ij}$ (Fig. 11 and Appendix D) shows that these properties are verified in our simulations.

We also compared the signs of the $\alpha$ tensor components to the study of the MRI by Gressel & Pessah (2015), who used the test field method to measure them, taking the contribution of all components into account. For the diagonal components, we also find that the sign of $\alpha_{\phi\phi}$ is opposite to that of $\alpha_{rr}$ (see Appendix D), which corresponds to the radial component $\alpha_{xx}$ in Gressel & Pessah (2015). Interestingly, in our models the off-diagonal elements of the $\alpha$ tensor have opposite signs across the domain with respect to their conjugate, but they do not have the same amplitude, contrary to the expectations of mean-field theory in the case of isotropic turbulence (Krause & Raedler 1980). By contrast, in local simulations of the MRI with stratification, they are found to have the same sign, which has been interpreted as a significant anisotropy of the MRI turbulence (Brandenburg 2008; Gressel & Pessah 2015). This different behavior might be due to the difference in geometry between a local model with a vertical stratification and a spherical global model with a radial stratification.

## 5.4. Diffusive processes

Ideal MHD without any explicit diffusion has sometimes been used to study the MRI, but convergence studies have shown that the strength of MRI turbulence then depends on resolution (Fromang & Papaloizou 2007; Pessah et al. 2007). To avoid this issue, we considered explicit diffusivities and used pseudo-spectral methods (which have low numerical dissipation), making sure that the diffusive scales are well resolved. This constraint forced us to use diffusivities that were sometimes higher than the expected values in a PNS. In particular, the value of the magnetic Prandtl number $Pm$ we adopted in our simulations is much lower than the value expected in a PNS. This is a significant limit to our models because the quantitative
results of the MRI turbulence are heavily dependent on $Pm$, as shown by studies of accretion disks (Lesur & Longaretti 2007; Fromang et al. 2007; Longaretti & Lesur 2010; Meheut et al. 2015) and in Paper I. The magnetic field strength that would be reached at higher $Pm$ is probably higher than the field strengths presented in this work, which should therefore be considered as lower bounds. The turbulent region may also be larger and occupy the entire domain in this regime, if its size is indeed affected by the magnetic diffusivity in the denser region.

Last, we assumed that neutrinos are in the diffusive regime, which is valid when the considered scales are longer than the neutrino mean free path. As the simulation domain is close to the PNS surface, the neutrino mean free path is large ($\approx 5 \times 10^5$ cm at $r = 33$ km) and most of the turbulence would therefore be impacted by a neutrino drag. In this regime, Guilet et al. (2015) showed that the neutrino drag, weaker in the outer layer of the PNS (for a radius $\ge 30$ km), does not impact the linear modes of the MRI. Therefore, the MRI is able to grow freely on small scales in this region. However, the impact of the neutrino drag has been studied only in the linear phase, but never in nonlinear simulations. It would be especially important to assess it if we consider that the turbulence in our simulations is close to the PNS surface, where the diffusive approximation is less valid. The evolution of the MRI in this regime is postponed to future studies.

### 5.5. Validity of the anelastic approximation

The anelastic approximation is used in our models to include density stratification with a reduced computational cost with respect to standard MHD, as it allows filtering out sound waves. The filtering of sound waves is justified when the different velocities of our system, such as fluid velocity $u$, Alfvén velocity $v_A$, and gravity wave $v_g$ velocity, are much lower than the sound speed $c_s$ of the fluid. From the 1D CCSN simulation used to compute our anelastic reference state, we can compute the sound speed, which ranges from $c_s \approx 4 \times 10^9$ cm s$^{-1}$ at the inner boundary to $c_s \approx 2.9 \times 10^9$ cm s$^{-1}$ at the outer boundary. Using the values of Alfvén speed shown in Fig. 20, we can verify a posteriori that $v_A^2/c_s^2 \le 10^{-2}$. Because the turbulent kinetic energy is lower than the magnetic energy, we also have $u^2/c_s^2 \le v_A^2/c_s^2 \le 10^{-2}$. For the gravity waves, we ran some hydrodynamic tests in the anelastic model with $\text{Pr} = 0.1$ and a rotation that was a thousand times slower. The results agree well with the frequency and damping rate expected from large-scale gravity modes ($l \le 6$). To justify the filtering of sound waves, we also verified a posteriori in our anelastic simulations that the oscillation frequency of these modes is lower than the Lamb frequency $w_{\text{Lamb}} = c_s k$, where $k = \sqrt{l(l+1)/r^2}$ is the horizontal wavenumber of the spherical harmonics. Furthermore, the LBR anelastic approximation (Lantz 1992; Braginsky & Roberts 1995) that we used describes the propagation of gravity modes well when compared to compressible equations (Brown et al. 2012). Last, we also verified that relative density perturbations $\delta\rho$ due to entropy perturbations are small: $\delta\rho/\bar{\rho} \sim \bar{\alpha}_T \bar{T} s' / (c_p \bar{\rho}) \ll 1$.

The main limitation in our implementation of the anelastic approximation is that the thermodynamical background was assumed not to evolve over the timescale of the simulation. Our simulations lasted approximately one second, but the 1D CCSN simulation shows that during this time, the width of the outer stably stratified region shrinks from approximately 15 km to 5 km, while the PNS radius contracts from 40 km to 20 km. It is difficult to take this structural evolution in our models into account.

We therefore postpone the study of MRI-driven dynamos at later times to future works.

## 6. Conclusions

We have investigated the effect of the density gradient and stable stratification on the generation of large-scale magnetic fields by the MRI in 3D spherical PNS models. We developed anelastic models using a thermodynamic background describing the stably stratified region of a PNS based on 1D CCSN simulations and compared our results to those obtained with incompressible and Boussinesq models. The main findings of our study are summarized below.

We obtain a quasi-stationary state with a MRI-driven dynamo in presence of density stratification and stable thermal stratification. The averaged magnetic field strength is $\sim 1.4 \times 10^{14}$ G and the (mostly) equatorial magnetic dipole represents $\sim 4.3\%$ of the averaged magnetic field. These results are consistent with our previous incompressible study when they are rescaled for the different densities.

The MRI-driven dynamo shows a qualitatively different behavior in presence of density stratification, with a prominent axisymmetric component of the magnetic field displaying oscillations on a timescale significantly longer than the rotation period. This mean-field dynamo and its oscillation period can be described by an $\alpha\Omega$ mechanism. The dynamo frequencies in our simulations agree relatively well with the theoretical dynamo frequency explained by an $\alpha$ effect, which is consistent with the theoretical calculations based on the local kinetic helicity.

A low thermal Prandtl number prevents the buoyancy from damping the MRI turbulence, except in the equatorial plane. When Pr is increased, the radial velocity is distributed at smaller scales, which can impact the MRI-driven dynamo and its $\alpha\Omega$ mechanism.

These findings support the hypothesis that the MRI is able to generate a strong large-scale magnetic field in presence of realistic physical ingredients. This means that it is an important mechanism that can help explain magnetar formation.

Although the magnetic dipole of $6 \times 10^{12}$ G obtained in our simulations may seem weak to form a magnetar, we note that these results are obtained with a large PNS with a radius of $r_0 = 39.25$ km, that is, before the contraction to the final size of a cold neutron star with a radius close to 12 km. Under the plausible hypothesis that magnetic flux is conserved during this contraction, the magnetic field could be amplified by a factor of $\sim 10$ and the dipole would then be close to the lower end of the magnetar range $10^{14}-10^{15}$ G. In addition to this dipolar magnetic field, the MRI in our anelastic models generates a large-scale toroidal field of $B \approx 8 \times 10^{13}$ G that can be amplified to $\sim 10^{15}$ G by flux conservation. Finally, we note that there are reasons to expect that our simulations may underestimate the intensity of the magnetic field generated by the MRI. Importantly, as highlighted in the discussion, a higher magnetic Prandtl number will quantitatively change the results, leading most likely to an increase in the magnetic field strength. Higher magnetic Prandtl numbers may also enable a fully turbulent state extending to the higher-density region at smaller radii, where we may expect a stronger magnetic field.

An important limitation of our approach concerns the evolution of the PNS structure. The 1D CCSN simulation we used shows that the PNS contracts in about 5 s and becomes almost fully convective with a thin stably stratified outer layer of a few kilometers. We did not take this evolution into account so far
because the required developments are beyond the scope of the present study. Forthcoming improvements will consist of modeling the entire PNS, which includes the stably stratified and convective zones. In the latter, the onset of a convective dynamo with fast rotation can generate magnetar-like magnetic fields (Raynaud et al. 2020). This leads to a strong magnetic field that is buried below the stably stratified zone, which could impact the MRI-driven dynamo. The stably stratified region may also influence the convective dynamo, as shown in some studies of planetary dynamos (Gastine & Wicht 2021). The interplay between the convective dynamo and the MRI-driven dynamo is therefore a key question of magnetorotational explosions because the magnetic field at the PNS surface has the strongest impact on the launch of the explosion (Obergaulinger & Aloy 2017, 2020; Bugli et al. 2020, 2021).

The core-collapse dynamics may be impacted by the large-scale magnetic field that we obtain in our simulations. The dipole field in our models may be slightly low for directly launching jets, especially because the equatorial component is less efficient than the aligned component, as shown by Bugli et al. (2021). Nonethe-less, other large-scale mean-field structures, such as the toroidal mean magnetic field, may still impact the supernova dynamics. The mean-field dynamo in our simulations opens exciting perspectives for modeling the generation of large-scale magnetic fields in core-collapse simulations. Our results can be used to cal-ibrate a subgrid model of the MRI-driven turbulence with an $\alpha\Omega$ dynamo mechanism. This would allow models to describe an in-situ amplification of the PNS magnetic field instead of relying on an unrealistically strong initial magnetic field.

Magnetic fields are also important in the context of binary neutron star mergers because magnetars are invoked to power short gamma-ray bursts and kilonovae, such as the GW170817 event (Metzger et al. 2018b). A stable magnetar could indeed explain the high luminosity of the kilonova associated with the recent short gamma-ray burst GRB200522A (Fong et al. 2021). This interpretation was also invoked to explain an X-ray transient as the aftermath of a binary neutron star merger (Xue et al. 2019). To support this scenario, the MRI has been invoked to amplify the magnetic field, as similar conditions in terms of neutrino radiation and differential rotation can be found in neutron star mergers (Guilet et al. 2017). Realistic simulations have shown that the MRI amplifies the magnetic field from an initial strong axial dipole (Siegel et al. 2013; Kiuchi et al. 2018; Ciolfi et al. 2019; Mösta et al. 2020). It is difficult to study large-scale field generation in realistic models of neutron star mergers because it requires taking many different physical ingredients into account: general relativity, treatment of neutrino physics, equation of state of hot and dense matter, and MHD. The development of idealized models with methods similar to those employed in this study will help to understand the effects of different physical processes and provide a useful reference for comparison and calibration of $\alpha\Omega$ dynamos in merger simulations (Shibata et al. 2021).

Investigating the different scenarios of magnetar formation is a promising avenue of research as more statistics will be available for transients events in the multi-messenger era. For instance, new statistics on FRBs and their host galaxies may give new insights on magnetar formation in further galaxies because magnetars are at least one of the FRB progenitors.

*Acknowledgements.* We thank Hans-Thomas Janka who provided the data of 1D CCSN supernova simulations. 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 A0050410317, A0070410317 and A0090410317). The authors are grateful to Thomas Gastine, Thierry Foglizzo and Sébastien Fromang for their valuable insight and expertise.

## References

Acheson, D. J., & Gibbons, M. P. 1978, Phil. Trans. R. Soc. London Ser. A, 289, 459

Akiyama, S., Wheeler, J. C., Meier, D. L., & Lichtenstadt, I. 2003, ApJ, 584, 954
Balbus, S. A. 1995, ApJ, 453, 380

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Bochenek, C. D., Ravi, V., Belov, K. V., et al. 2020, Nature, 587, 59

Braginsky, S. I., & Roberts, P. H. 1995, Geophys. Astrophys. Fluid Dyn., 79, 1

Brandenburg, A. 2008, Astron. Nachr., 329, 725

Brown, B. P., Vasil, G. M., & Zweibel, E. G. 2012, ApJ, 756, 109

Bucciantini, N., Metzger, B. D., Thompson, T. A., & Quataert, E. 2012, MNRAS, 419, 1537

Bugli, M., Guilet, J., Obergaulinger, M., Cerdá-Durán, P., & Aloy, M. A. 2020, MNRAS, 492, 58

Bugli, M., Guilet, J., & Obergaulinger, M. 2021, MNRAS, 507, 443

Busse, F. H., & Simitev, R. D. 2006, Geophys. Astrophys. Fluid Dyn., 100, 341

CHIME/FRB Collaboration (Andersen, B. C., et al.) 2020, Nature, 587, 54

Christensen, U., & Wicht, J. 2015, in Treatise on Geophysics (Second Edition), ed. G. Schubert (Oxford: Elsevier), 245

Ciolfi, R., Kastaun, W., Kalinani, J. V., & Giacomazzo, B. 2019, Phys. Rev. D, 100, 023005

Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. 2018, MNRAS, 474, 961

Coti Zelati, F., Borghese, A., Israel, G. L., et al. 2021, ApJ, 907, L34

Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52

Deng, H., Mayer, L., Latter, H., Hopkins, P. F., & Bai, X.-N. 2019, ApJS, 241, 26

Dessart, L., Burrows, A., Livne, E., & Ott, C. D. 2008, ApJ, 673, L43

Drout, M. R., Soderberg, A. M., Gal-Yam, A., et al. 2011, ApJ, 741, 97

Duncan, R. C., & Thompson, C. 1992, ApJ, 392, L9

Esposito, P., Rea, N., & Israel, G. L. 2021, in Magnetars: A Short Review and Some Sparse Considerations, eds. T. M. Belloni, M. Méndez, & C. Zhang, 461, 97

Fong, W., Laskar, T., Rastinejad, J., et al. 2021, ApJ, 906, 127

Fromang, S., & Papaloizou, J. 2007, A&A, 476, 1113

Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123

Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661

Gastine, T., & Wicht, J. 2012, Icarus, 219, 428

Gastine, T., & Wicht, J. 2021, Icarus, 368, 114514

Gastine, T., Duarte, L., & Wicht, J. 2012, A&A, 546, A19

Gilman, P. A., & Glatzmaier, G. A. 1981, ApJS, 45, 335

Gompertz, B. P., O'Brien, P. T., & Wynn, G. A. 2014, MNRAS, 438, 240

Götz, D., Mereghetti, S., Molkov, S., et al. 2006, A&A, 445, 313

Gouhier, B., Lignières, F., & Jouve, L. 2021, A&A, 648, A109

Gouhier, B., Jouve, L., & Lignières, F. 2022, A&A, 661, A119

Gressel, O., & Pessah, M. E. 2015, ApJ, 810, 59

Gubbins, D., & Zhang, K. 1993, Phys. Earth Planet. Int., 75, 225

Guilet, J., & Müller, E. 2015, MNRAS, 450, 2153

Guilet, J., Müller, E., & Janka, H.-T. 2015, MNRAS, 447, 3992

Guilet, J., Bauswein, A., Just, O., & Janka, H.-T. 2017, MNRAS, 471, 1879

Hawley, J. F., & Balbus, S. A. 1992, ApJ, 400, 595

Hüdepohl, L. 2014, PhD Thesis, Technische Universität, München, Germany

Hurley, K., Boggs, S. E., Smith, D. M., et al. 2005, Nature, 434, 1098

Inserra, C., Smartt, S. J., Jerkstrand, A., et al. 2013, ApJ, 770, 128

Israel, G. L., Belloni, T., Stella, L., et al. 2005, ApJ, 628, L53

Jones, C. A., Boronski, P., Brun, A. S., et al. 2011, Icarus, 216, 120

Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261

Kiuchi, K., Kyutoku, K., Sekiguchi, Y., & Shibata, M. 2018, Phys. Rev. D, 97, 124039

Kouveliotou, C., Dieters, S., Strohmayer, T., et al. 1998, Nature, 393, 235

Krause, F., & Raedler, K. H. 1980, Mean-field Magnetohydrodynamics and Dynamo Theory (Oxford: Pergamon Press)

Kuroda, T., Arcones, A., Takiwaki, T., & Kotake, K. 2020, ApJ, 896, 102

Lantz, S. R. 1992, PhD Thesis, Cornell University, USA

Lattimer, J. M., & Swesty, D. F. 1991, Nucl. Phys. A, 535, 331

Lesur, G., & Longaretti, P. Y. 2007, MNRAS, 378, 1471

Lesur, G., & Ogilvie, G. I. 2008, A&A, 488, 451

Longaretti, P. Y., & Lesur, G. 2010, A&A, 516, A51

Lorimer, D. R., Bailes, M., McLaughlin, M. A., Narkevic, D. J., & Crawford, F. 2007, Science, 318, 777

Margalit, B., Metzger, B. D., Thompson, T. A., Nicholl, M., & Sukhbold, T. 2018, MNRAS, 475, 2659
Masada, Y., Sano, T., & Shibata, K. 2007, ApJ, 655, 447

Masada, Y., Takiwaki, T., Kotake, K., & Sano, T. 2012, ApJ, 759, 110

Masada, Y., Takiwaki, T., & Kotake, K. 2015, ApJ, 798, L22

Masada, Y., Takiwaki, T., & Kotake, K. 2022, ApJ, 924, 75

Meheut, H., Fromang, S., Lesur, G., Joos, M., & Longaretti, P.-Y. 2015, A&A, 579, A117

Menou, K., Balbus, S. A., & Spruit, H. C. 2004, ApJ, 607, 564

Metzger, B. D., Quataert, E., & Thompson, T. A. 2008, MNRAS, 385, 1455

Metzger, B. D., Giannios, D., Thompson, T. A., Bucciantini, N., & Quataert, E. 2011, MNRAS, 413, 2031

Metzger, B. D., Beniamini, P., & Giannios, D. 2018a, ApJ, 857, 95

Metzger, B. D., Thompson, T. A., & Quataert, E. 2018b, ApJ, 856, 101

Moffatt, H. K. 1978, Cambridge Monographs on Mechanics and Applied Mathematics (Cambridge: Cambridge University Press)

Moiseenko, S. G., Bisnovatyi-Kogan, G. S., & Ardeljan, N. V. 2006, MNRAS, 370, 501

Mösta, P., Richers, S., Ott, C. D., et al. 2014, ApJ, 785, L29

Mösta, P., Ott, C. D., Radice, D., et al. 2015, Nature, 528, 376

Mösta, P., Radice, D., Haas, R., Schnetter, E., & Bernuzzi, S. 2020, ApJ, 901, L37

Nicholl, M., Smartt, S. J., Jerkstrand, A., et al. 2013, Nature, 502, 346

Obergaulinger, M., & Aloy, M. Á. 2017, MNRAS, 469, L43

Obergaulinger, M., & Aloy, M. Á. 2020, MNRAS, 492, 4613

Obergaulinger, M., Cerdá-Durán, P., Müller, E., & Aloy, M. A. 2009, A&A, 498, 241

Obergaulinger, M., Just, O., & Aloy, M. A. 2018, J. Phys. G Nucl. Part. Phys., 45, 084001

Ott, C. D., Burrows, A., Thompson, T. A., Livne, E., & Walder, R. 2006, ApJS, 164, 130

Parker, E. N. 1955, ApJ, 122, 293

Pessah, M. E., & Chan, C.-K. 2008, ApJ, 684, 498

Pessah, M. E., Chan, C.-K., & Psaltis, D. 2007, ApJ, 668, L51

Rampp, M., & Janka, H. T. 2002, A&A, 396, 361

Raynaud, R., Guilet, J., Janka, H. T., & Gastine, T. 2020, Sci. Adv., 6, eaay2732

Raynaud, R., Cerdá-Durán, P., & Guilet, J. 2021, MNRAS, 509, 3410

Reboul-Salze, A., Guilet, J., Raynaud, R., & Bugli, M. 2021, A&A, 645, A109

Rembiasz, T., Guilet, J., Obergaulinger, M., et al. 2016, MNRAS, 460, 3316

Roberts, O. J., Veres, P., Baring, M. G., et al. 2021, Nature, 589, 207

Rodríguez Castillo, G. A., Israel, G. L., Tiengo, A., et al. 2016, MNRAS, 456, 4145

Rowlinson, A., O'Brien, P. T., Metzger, B. D., Tanvir, N. R., & Levan, A. J. 2013, MNRAS, 430, 1061

Schaeffer, N. 2013, Geochem. Geophys. Geosyst., 14, 751

Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211

Shi, J.-M., Stone, J. M., & Huang, C. X. 2016, MNRAS, 456, 2273

Shibata, M., Liu, Y. T., Shapiro, S. L., & Stephens, B. C. 2006, Phys. Rev. D, 74, 104026

Shibata, M., Fujibayashi, S., & Sekiguchi, Y. 2021, Phys. Rev. D, 104, 063026

Siegel, D. M., Ciolfi, R., Harte, A. I., & Rezzolla, L. 2013, Phys. Rev. D, 87, 121302

Simon, J. B., & Hawley, J. F. 2009, ApJ, 707, 833

Svinkin, D., Frederiks, D., Hurley, K., et al. 2021, Nature, 589, 211

Thompson, C., & Duncan, R. C. 1993, ApJ, 408, 194

Tiengo, A., Esposito, P., Mereghetti, S., et al. 2013, Nature, 500, 312

Tilgner, A., & Busse, F. H. 1997, J. Fluid Mech., 332, 359

Turolla, R., Zane, S., & Watts, A. L. 2015, Rep. Progr. Phys., 78, 116901

Wicht, J. 2002, Phys. Earth Planet. Inter., 132, 281

Winteler, C., Käppeli, R., Perego, A., et al. 2012, ApJ, 750, L22

Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015

Xue, Y. Q., Zheng, X. C., Li, Y., et al. 2019, Nature, 568, 198

Yoshimura, H. 1975, ApJ, 201, 740
**Appendix A: Reference state of the anelastic model**

The comparison of the reference model in the simulation to the mixing-length theory model of 1D CCSN simulations is presented in this appendix.

<figure><img src="image_23.png" /><figcaption>**Fig. A.1.** Comparison of the reference state considered in this study (magic model) with the mlt model of the 1D CCSN simulation by Hüdepohl (2014). (a) Background density profile, (b) effective entropy profile, (c) temperature profile, (d) gravity profile, and (e) thermal expansion profile as a function of the normalized radius. Dimensional and dimensionless units are converted by simple multiplication by the reference values at the outer boundary.</figcaption></figure>
## Appendix B: Parameters of our global simulations

**Table B.1.** Overview of the global models. All simulations use $Pm = 16$ and insulating magnetic boundary conditions. The $(n_r, n_\theta, n_\phi) = (257, 512, 1024)$ resolution is called “Medium” and the $(n_r, n_\theta, n_\phi) = (385, 768, 1536)$ resolution is called “High”. $Ra$ corresponds to the Rayleigh number linked to the Brunt-Väisälä frequency (see Eq. 21), $Pr$ to the thermal Prandtl number, $B_o$ to the initial magnetic field strength at the outer boundary. Magnetic field and dipole field strengths are time and space averaged. $P_{sim}$ corresponds to the estimated period of the mean-field dynamo patterns at $r = 0.92 r_o$, $k_z$ to the vertical wavenumber and $q\Omega$ to the value in the turbulent region (at $r = 0.92 r_o$). $k_z$, $q\Omega$ and the alpha component $\alpha_{\phi\phi}$ are used to compute the theoretical $\alpha\Omega$ period $P_{\alpha\Omega}$.

<table><thead><tr><th>Model</th><th>Resolution</th><th>Ra</th><th>Pr</th><th>B<sub>o</sub><br>[10<sup>14</sup>G]</th><th>Magnetic field<br>[10<sup>14</sup>G]</th><th>Dipole<br>[10<sup>14</sup>G]</th><th>P<sub>sim</sub><br>[ms]</th><th>k<sub>z</sub><br>[cm<sup>-1</sup>]</th><th>qΩ<br>[rad.s<sup>-1</sup>]</th><th>α<sub>φφ</sub><br>[cm.s<sup>-1</sup>]</th><th>P<sub>αΩ</sub><br>[ms]</th></tr></thead><tbody><tr><td>Standard</td><td>Medium</td><td>-4.73e5</td><td>0.005</td><td>1.5</td><td>1.40</td><td>0.063</td><td>401</td><td>2.3e-6</td><td>385</td><td>7.7e5</td><td>342</td></tr><tr><td>High Res</td><td>High</td><td>-4.73e5</td><td>0.005</td><td>3.3</td><td>1.36</td><td>0.053</td><td>351</td><td>2.3e-6</td><td>454</td><td>7.0e5</td><td>329</td></tr><tr><td>Incompressible</td><td>High</td><td>0.0</td><td>-</td><td>1.5</td><td>2.36</td><td>0.085</td><td>NO</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Boussinesq</td><td>High</td><td>-4.73e5</td><td>0.005</td><td>1.5</td><td>1.85</td><td>0.0548</td><td>NO</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Anel Ra0</td><td>Medium</td><td>0.0</td><td>0.005</td><td>1.5</td><td>1.41</td><td>0.051</td><td>461</td><td>2.3e-6</td><td>369</td><td>7.8e6</td><td>349</td></tr><tr><td>Anel Pr0 02</td><td>Medium</td><td>-1.892e6</td><td>0.02</td><td>1.5</td><td>1.59</td><td>0.074</td><td>363</td><td>2.3e-6</td><td>479</td><td>8.0e5</td><td>300</td></tr><tr><td>Anel Pr0 1</td><td>Medium</td><td>-9.46e6</td><td>0.1</td><td>1.5</td><td>1.36</td><td>0.062</td><td>493</td><td>2.7e-6</td><td>480</td><td>4.5e5</td><td>366</td></tr></tbody></table>
Appendix C: Correlation coefficients between the EMF and the magnetic field

<figure><img src="image_24.png" /><figcaption>Fig. C.1. Time-averaged values of the correlation coefficient between the EMF and the magnetic field taken at r = 0.86 r_o.</figcaption></figure>
**Appendix D: Estimation of the components of the α tensor**

<figure><img src="image_25.png" /><figcaption>**Fig. D.1.** Time-averaged values of the α tensor components estimated by the correlation between the EMF and the magnetic field taken at $r = 0.86 r_o$.</figcaption></figure>
# Appendix E: Azimuthal average of magnetic and velocity fields at quasi-stationary state

## E.1. Simulation Standard

<figure><img src="image_26.png" /><figcaption>**Fig. E.1.** Snapshots of the azimuthal average of magnetic field components and velocity field components for the simulation Standard. Top: Snapshots of the azimuthal average of the magnetic field components $\bar{B}_{\phi}$, $\bar{B}_{r}$, and $\bar{B}_{\theta}$ at $t = 773$ ms for the model Standard. Bottom: Snapshots of the azimuthal average of the velocity field components $\bar{u}_{r}$, $\bar{u}_{\theta}$, $\bar{u}_{\phi}$, and the turbulent kinetic helicity $h$. The angular velocity component can be found in Fig. 12.</figcaption></figure>
E.2. Simulation Boussinesq

<figure><img src="image_27.png" /><figcaption>**Fig. E.2.** Snapshots of the azimuthal average of magnetic field components and velocity field components for the simulation Boussinesq. Top: Snapshots of the azimuthal average of the magnetic field components $\tilde{B}_{\phi}$, $\tilde{B}_{r}$, and $\tilde{B}_{\theta}$ at $t = 511$ ms for model Boussinesq. Bottom: Snapshots of the azimuthal average of the velocity field components $\tilde{u}_{r}$, $\tilde{u}_{\theta}$, $\tilde{u}_{\phi}$, and the turbulent kinetic helicity $h$. The angular velocity can be found in Fig. 12.</figcaption></figure>
E.3. Simulation Pr=0.1

<figure><img src="image_28.png" /><figcaption>Fig. E.3. Snapshots of the azimuthal average of magnetic field components and velocity field components for the simulation Pr=0.1. Top: Snapshots of the azimuthal average of the magnetic field components $\tilde{B}_{\phi}$, $\tilde{B}_{r}$, and $\tilde{B}_{\theta}$ at $t = 1128$ ms for model Anel Pr0.1. Bottom: Snapshots of the azimuthal average of the velocity field components $\tilde{u}_{r}$, $\tilde{u}_{\theta}$, $\tilde{u}_{\phi}$, and the turbulent kinetic helicity $h$.</figcaption></figure>