Issue |
A&A
Volume 667, November 2022
|
|
---|---|---|
Article Number | A94 | |
Number of page(s) | 26 | |
Section | Astrophysical processes | |
DOI | https://doi.org/10.1051/0004-6361/202142368 | |
Published online | 11 November 2022 |
MRI-driven αΩ dynamos in protoneutron stars
1
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
2
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), 14476 Potsdam, Germany
3
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
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 1014–1015 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 hypernovae, 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 αΩ 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 αΩ 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
© A. Reboul-Salze et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe-to-Open model.
This Open access funding provided by Max Planck Society.
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 Ṗ, two important timing parameters from which a magnetic field
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) associated with long gamma-ray bursts (Drout et al. 2011) and superluminous 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 millisecond 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 superluminous 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 transients 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 cannot explain the formation of millisecond magnetars. To have both fast rotation and a strong magnetic field, an in-situ magnetic 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 magnetohydrodynamics (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), numerical 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 hundred milliseconds; Obergaulinger et al. 2009; Guilet & Müller 2015; Rembiasz et al. 2016). In protoneutron stars, the physical 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 viscosity (Masada et al. 2007; Guilet et al. 2015). Finally, the MRI turbulence 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 intensities 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 nonlinear phase of the MRI and in Sect. 4 for the comparison between different global models of the MRI. Finally, we discuss the validity of our assumptions in Sect. 5 and draw our conclusions in Sect. 6.
2. Numerical setup
2.1. 1D protoneutron star model
The simulations performed in this article are designed to represent the outer region of a fast-rotating PNS, which is stable to the convection and unstable to the MRI. To model this stably 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⊙ 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⊙ progenitor s27.0 by Woosley et al. (2002). The calculations were performed with the code Prometheus-Vertex, which combines the 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 ri = 25.5 km to ro = 39.25 km, corresponding to the PNS surface defined by the density ρo = 1011 g cm−3. All background profiles of density, temperature, entropy gradient, and gravitational acceleration 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
was computed using the following thermodynamic relation:
![]() |
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. |
where cp = 3.6 × 108 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, κ and ν, respectively. The magnetic diffusivity η was also explicitly included in our simulations. Finally, we assumed all the diffusivities κ, ν, and η 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 Ω0 = 958 rad s−1 read (Jones et al. 2011)
where u is the flow velocity, B is the magnetic field, P is the pressure, s′ is the entropy perturbation, and μ0 is the vacuum permeability. In this system, viscous force and viscous heating are given by and Φν = ∂juiσij, where
is the rate of strain tensor and eij = (∂jui + ∂iuj)/2 is the deformation tensor. Tensors are expressed with the Einstein summation convention and the Kronecker symbol δij.
2.3. Numerical methods
In order to integrate the system of Eqs. (3)–(7) in time, we used the benchmarked pseudo-spectral code MagIC1 (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,
where W and Z are the poloidal and toroidal kinetic scalar potentials, respectively, while b and aj are the magnetic potentials. The scalar potentials and the pressure P are decomposed on spherical harmonics for the colatitude θ and the longitude ϕ 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 (nr, nθ, nϕ) = (257, 512, 1024) or a higher resolution of (nr, nθ, nϕ) = (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ν ≃ 70 and lη ≃ 100, respectively.
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,
where s is the cylindrical radius, qo = 1.5 corresponds to the shear rate
at the outer boundary, and Ωi = 3822 rad s−1 is the rotation rate of the inner core. This rotation rate Ωi was computed so that the ratio of total angular momentum over the moment of inertia was equal to the frame rotation rate Ω0 defined in Sect. 2.2. This was computed with the following formula:
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 qo = 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
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,
![]() |
Fig. 2. Ratio of the squared Brunt–Väisälä frequency N2 (Eq. (16)) to the squared initial rotation frequency Ω2 (Eq. (10)). The buoyancy influence is strongest in the red regions. |
The initial root mean-square magnetic field strength at the outer boundary ranged from Bo = 1.5 × 1014 G to Bo = 3.3 × 1014 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 Bo = 1.5 × 1014 G reads (Balbus & Hawley 1991)
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 (ur = 0). At the inner boundary, we used a stress-free condition, where the viscous stress vanished. For the outer boundary, we forced uϕ 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 N2 over the rotation frequency squared Ω2, with
where Ye is the electron fraction and 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 N2/Ω2 ∼ 103 − 104 in the radiative zone of intermediate and massive stars (Fuller et al. 2019). The Brunt–Väisälä frequency at the outer boundary is No = 4136 s−1. Accordingly, we might expect the buoyancy to dampen MRI-driven turbulence (Guilet & Müller 2015), but the thermal diffusivity κ also plays an important role as it reduces the impact of the stable stratification on radial scales smaller than a critical length Lc. One way to estimate Lc is to compare the timescale for thermal diffusion to the timescale of gravity waves, from which we obtain
![]() |
Fig. 3. Meridional slices at ϕ = 0 of the initial poloidal magnetic field Br (left) and of the toroidal magnetic field Bϕ at t = 773 ms (right) for model Standard. |
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 ≃ 3.3 × 106 km, which are κ = 1.61 × 1014 cm s−2 and ν = 8.03 × 1011 cm s−2, respectively. While all our models share the same constant kinematic viscosity, we considered three different values for the thermal diffusivity, κ = {1.61 × 1014, 4.02 × 1013, 8.03 × 1012} cm s−2, which means that the thermal Prandtl number is
Our standard value of the thermal Prandtl number 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 Lc = 0.14L = 1.97 × 105 cm to Lc = 0.031L = 4.4 × 104 cm. In our simulations, the ratio of the Brunt–Väisälä frequency No at the outer boundary to frame rotation Ω0 was fixed by the following formula:
where the Ekman number E (characterizing the importance of viscosity over the Coriolis force) reads
and the Rayleigh number Ra is defined by
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 × 105 to Ra = −9.46 × 106.
With these dimensionless numbers, we can estimate the regime of boundary forcing in our simulations. The viscous timescale τν = L2/ν can be compared to the Ekman spin-up timescale and the Eddington-Sweet circulation timescale
, (Gouhier et al. 2021, 2022). This gives
for the simulation Standard. This order (i.e., τE ≪ τED ≪ τν) corresponds to the Eddington-Sweet regime, but as Pr increases for some simulations, it approaches the viscous regime (i.e., τE ≪ τν ≪ τED).
We used a constant resistivity of η = 5.0 × 1010 cm s−2, which gives a magnetic Prandtl number
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 ≈ 1013 (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
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.
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 Bo = 1.5 × 1014 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 ur. 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 and saturated magnetic field shown in Fig. 3 suggests that the magnetic field has lost the memory of its initial configuration.
![]() |
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. |
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 ≃ 1.4 × 1014 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 ≃ 2.25 × 1014 G.
![]() |
Fig. 5. Time evolution of the magnetic and turbulent kinetic energy densities for model Standard with Bo = 1.5 × 1014 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. |
In order to compare the intensity in the full volume to our previous incompressible study (Paper I), where the density ρo = 4 × 1013 g cm−3 was higher, we used the dimensionless magnetic field strength defined by the Lorentz number,
We have ℬ ≃ 0.096 for model Standard and ℬ ≃ 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 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.
![]() |
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. |
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 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 Bdip ≃ 6.3 × 1012 G, which represents ≃4.5% of the total magnetic field (Fig. 7). This intensity may seem weak, but the dimensionless magnetic field strength ℬdip corresponding to this dipole intensity gives ℬdip = 0.0043, which is higher than the value ℬdip ≃ 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 θdip ≈ 100° 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.
![]() |
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. |
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 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.
![]() |
Fig. 8. Radius-time diagram of the dipolar modes. Top: radius-time diagram of the axial dipolar mode Br(l = 1, m = 0). The period of the dipole reversal is Pdip ≃ 410 ms. Bottom: radius-time diagram of the equatorial dipolar mode amplitude Br(l = 1, m = 1) at a given azimuth ϕ = 0 in the rotating frame of the simulation. The fast oscillations are due to the rotation of this dipole mode. |
3.2. αΩ 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: . The definition of mean here is the axisymmetric average operator noted
, which verifies the Reynolds averaging rules. The approach of the mean-field theory is to expand the electromotive force (EMF) ℰ only in terms of the mean quantities (
and
) 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 αΩ dynamo. The Ω 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 Ω effect reads
and it should induce an anticorrelation between the radial field in cylindrical coordinates and the azimuthal magnetic field
.
The α effect comes instead from the closure relation of the mean EMF,
where and
, expressed as a function of the mean magnetic field
,
where αij and βij are tensors that do not depend on , and i, j refer to spherical r, θ, ϕ or cylindrical coordinates s, ϕ, z. The diagonal components of the α 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ϕ in the form of . The diagonal components of the β tensor are in the direction of the mean current
, and their effect is physically described as a turbulent diffusivity, which adds to the magnetic diffusivity η. Another effect that has been proposed to complete the dynamo loop in the case of the MRI is the nondiagonal resistivity βϕs, which could generate a poloidal field from the toroidal field (Lesur & Ogilvie 2008). In this case, the azimuthal component of the electromotive force ℰϕ would be correlated with the radial current
in cylindrical coordinates. One (or both) of these two effects, in combination with the Ω effect, may therefore be expected to complete the dynamo cycle, leading to a self-sustained magnetic field.
The usual way to characterize an αΩ dynamo is to compute space-time diagrams of Bϕ, Bs, and ℰϕ, which are often referred to as butterfly diagrams. Figure 9 shows these azimuthally averaged quantities in the southern hemisphere at r = 0.86 ro (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 ≃ 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 αΩ dynamo. The components and
are anticorrelated, which shows that
is mainly generated by the Ω effect. On the other hand, the electromotive force ℰϕ is correlated with
, which suggests that the α effect plays an important role.
![]() |
Fig. 9. Butterfly diagrams of |
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:
where ⟨ ⋅ ⟩t represents a time average. Figure 10 shows the correlation coefficients between ℰϕ and and the radial current
in order to test the nondiagonal resistivity hypothesis. We find that
and ℰϕ are well correlated, while no correlation is found between
and ℰϕ. This means that the α effect is prominent in our simulations. Moreover, the antisymmetry of the correlation matches the expected symmetry of the components of the tensor α.
![]() |
Fig. 10. Time-averaged correlation coefficients (Eq. (29)) between ℰϕ and |
We can estimate the value of the diagonal components of the α tensor with the formula
This estimation assumes that the EMF is only due to the α 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 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 = 𝒱τ/ℒ, where 𝒱, τ and ℒ are typical values of the velocity, time variation, and length scale of the turbulence. Under these hypotheses, the diagonal components of the α tensor are proportional to the kinetic helicity h,
where τ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 simulation2, we find a quantitative agreement between Eqs. (30) and (31) for τc ≃ 2.5 ms = 0.38 × 2π/Ω0 with a peak value of αϕϕ ≃ 6 × 105 cm s−1 (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.
![]() |
Fig. 11. Time-averaged values of the αϕϕ component estimated with Eq. (30), (blue) and with the turbulent kinetic helicity (orange), (Eq. (31)) taken at r = 0.86 ro and with τc = 2.5 ms. |
Finally, the value of αϕϕ can be used to estimate the theoretical frequency of an αΩ dynamo using the relation (e.g., Busse & Simitev 2006; Gastine et al. 2012; Gressel & Pessah 2015)
where kz is the vertical wavenumber. We computed the shear rate q = −0.65 and the rotation rate Ω = 688 rad s−1 in the middle of the turbulent region from a vertically and azimuthally averaged rotation profile, which gives lower values than qo and Ω0. We estimated kz ≃ 2.79 × 10−06 cm−1 by taking the longest vertical length in the turbulent region at mid-latitudes and obtained a period of PαΩ = 2π/ωαΩ = 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 α effect, as we take its peak value in the turbulent region. If we consider an average on the angles θ ∈ [30° ,60° ], we obtain a period PαΩ = 393 ms, which is closer to the measured period P ≃ 410 ms. Moreover, according to the theoretical Parker-Yoshimura rule (Parker 1955; Yoshimura 1975), the dynamo wave should propagate in the direction of −αϕϕeϕ × ∇Ω, 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 αΩ 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 (N2 > 0 or N2 = 0). For Boussinesq models, constant density was taken to be equal to the density ρo = 1011 g cm−3 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ϕ 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 throught 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.
![]() |
Fig. 12. Meridional cuts of Bϕ (top) and azimuthal average of Ω (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. |
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.
![]() |
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. |
In a similar fashion to axisymmetric magnetic energy, the dipole energy density seems to be on the same order for all 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.
![]() |
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. |
To confirm the mean-field behavior, we show in Fig. 15 the butterfly diagrams of . 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 αΩ behavior. There is no clear signal in the toroidal component of the EMF ℰϕ, and we find that the correlation coefficients between the EMF ℰϕ and
are low, with an average
in both hemispheres. We find no correlations either between the EMF ℰϕ and
. 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.
![]() |
Fig. 15. Butterfly diagrams of |
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ϕ 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 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ϕ 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 mid-latitudes (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.
Periods of the different mean-field patterns for Pr = 0.005, Pr = 0.02, and Pr = 0.1.
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 Lc defined by Eq. (17). For our simulation Standard at Pr = 0.005, the critical length at mid-latitudes Lc ≃ 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 Lc ≃ 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 Lc). 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.
![]() |
Fig. 16. Radial velocity ur and αϕϕ for anelastic models with buoyancy and for increasing Prandtl numbers. Top: snapshots of the radial velocity ur. 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 αϕϕ for Pr = 0.005 (left), Pr = 0.02 (middle), and Pr = 0.1 (right). The dashed line represents the spherical radius r = 0.92 ro. |
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 well, corresponding to the critical length.
![]() |
Fig. 17. Normalized volume-averaged spectrum of the radial velocity squared ur 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 lc ≈ R/Lc above which the scales are impacted by buoyancy. |
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).
![]() |
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. |
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 mid-latitudes. The main difference is the weaker turbulence in the equatorial plane for model Standard, while no turbulence can 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 αΩ dynamo, we compare in Fig. 19 the αϕϕ component for the different Prandtl numbers. At Pr = 0.02, we find a mean αϕϕ ≃ 8.0 × 105 cm s−1 in the turbulent region, which translates into a period PαΩ ≃ 300 ms using the same kz 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.
![]() |
Fig. 19. Time-averaged values of the αϕϕ component estimated using Eq. (30) at r = 0.92 ro in the turbulent region for three different Prandtl numbers: Pr = 0.005 (blue), Pr = 0.02 (orange), and Pr = 0.1 (green). |
The 2D meridional distribution of αϕϕ (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 ur 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 α 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 PNorth ≃ 493 ms. However, in the southern hemisphere, we observe mean-field patterns of similar amplitude, but higher frequency (with a period of PSouth ≃ 213 ms). Moreover, the pattern propagates toward the equator, which is in the opposite direction to the Parker-Yoshimura rule. The meridional distribution of αϕϕ (bottom left panel of Fig. 16) also highlights the smaller turbulent region and the lower amplitude of αϕϕ than in the other two models.
At Pr = 0.1, we measure αϕϕ ≃ 4.57 × 105 cm s−1 in the northern and αϕϕ ≃ 4.47 × 105 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αΩ = 366 ms in the northern hemisphere, which agrees within 27% with the numerical value. The similar αϕϕ 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 αΩ 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 αΩ dynamo occurs in the northern hemisphere, but with a weaker α 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 and the dimensionless total magnetic field strength
instead of the magnetic intensities. For anelastic simulations, we find a greater magnetic strength ℬtot and dipole strength ℬdip than in Paper I. However, the ratio of ℬdip to ℬtot of anelastic simulations is quite similar, if slightly lower (≃4.3%) than the linear relation of Paper I.
For our Boussinesq simulations, the differences with Paper I are stronger. The magnetic strength ℬtot is higher by a factor of approximately three due to the stronger turbulence, and the dipole strength ℬdip is twice stronger. This gives a different linear relation of ≃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.
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 (nr, nθ, nϕ) = (385, 768, 1536) and a magnetic field that was twice as high Bo = 3.3 × 1014 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. αΩ versus α2Ω dynamo
As shown in the previous sections, the dynamos observed in our anelastic models are consistent with an αΩ mechanism. Although the anticorrelation of Bs and Bϕ in the butterfly diagrams (Fig. 9) suggests that the Ω effect generates the toroidal field, we would like to assess this more quantitatively. For this purpose, we compared the Ω effect to the generation of toroidal field by the diagonal components αrr and αθθ (see Appendix D), which would be relevant in the case of an α2Ω dynamo. To distinguish between an αΩ and α2Ω dynamo, we computed the ratio of the two dynamo numbers Cα = max(αrr, αθθ)R/η and CΩ = ΩR2/η in the turbulent region, which gives
where R is taken as the middle radius (ri + ro)/2. This estimate therefore confirms that the Ω effect largely dominates the generation of the toroidal field and that the dynamo in our simulations is an αΩ dynamo.
5.3. Properties of the α tensor components
The MRI-driven αΩ dynamo in our simulations was characterized by estimating the tensor components αij. Our method assumes that the EMF component ℰi is only due to the contribution of the mean magnetic field component Bj and neglects the contribution of the other components of the α tensor and the turbulent resistivity tensor β components. In the azimuthal EMF ℰϕ, the contribution from αϕr and αϕθ can be neglected because the Pearson correlation coefficients between ℰϕ and Bθ or Br are much lower in the turbulent region than they are for Bϕ (panels c and d of Appendix C). The αϕϕ 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 αϕi components (Appendix D) may be incorrectly estimated, and they do not drive the generation of the poloidal field. The other α components do not dominate the generation of toroidal field because the dynamo cannot be described as an α2Ω dynamo.
In the case of a statistically symmetric velocity field with respect to the equator, mean-field theory predicts that the components of the αij tensor (Eq. (28)) are either equatorial-symmetric (αrθ, αθr, αϕθ, αθϕ) or equatorial-antisymmetric (αrr, αϕr, αϕϕ, αrϕ, αθθ), (Gubbins & Zhang 1993). The estimation of the components αij (Fig. 11 and Appendix D) shows that these properties are verified in our simulations.
We also compared the signs of the α 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 αϕϕ is opposite to that of αrr (see Appendix D), which corresponds to the radial component αxx in Gressel & Pessah (2015). Interestingly, in our models the off-diagonal elements of the α 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 (≃5 × 105 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 ≥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 vA, and gravity wave vg velocity, are much lower than the sound speed cs of the fluid. From the 1D CCSN simulation used to compute our anelastic reference state, we can compute the sound speed, which ranges from cs ≃ 4 × 109 cm s−1 at the inner boundary to cs ≃ 2.9 × 109 cm s−1 at the outer boundary. Using the values of Alfvén speed shown in Fig. 20, we can verify a posteriori that . Because the turbulent kinetic energy is lower than the magnetic energy, we also have
. For the gravity waves, we ran some hydrodynamic tests in the anelastic model with 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 ≤ 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 wLamb = csk, where
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 δρ due to entropy perturbations are small:
.
![]() |
Fig. 20. Dimensionless dipole strength ℬdipole as a function of the dimensionless magnetic field strength ℬ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 ℬ were used in order to compare results with different densities |
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 ∼1.4 × 1014 G and the (mostly) equatorial magnetic dipole represents ∼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 αΩ mechanism. The dynamo frequencies in our simulations agree relatively well with the theoretical dynamo frequency explained by an α 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 αΩ 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 × 1012 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 ro = 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 ∼10 and the dipole would then be close to the lower end of the magnetar range 1014 − 1015 G. In addition to this dipolar magnetic field, the MRI in our anelastic models generates a large-scale toroidal field of B ≃ 8 × 1013 G that can be amplified to ∼1015 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). Nonetheless, 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 calibrate a subgrid model of the MRI-driven turbulence with an αΩ 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 αΩ 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.
Acknowledgments
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 [CrossRef] [Google Scholar]
- Akiyama, S., Wheeler, J. C., Meier, D. L., & Lichtenstadt, I. 2003, ApJ, 584, 954 [NASA ADS] [CrossRef] [Google Scholar]
- Balbus, S. A. 1995, ApJ, 453, 380 [NASA ADS] [CrossRef] [Google Scholar]
- Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
- Bochenek, C. D., Ravi, V., Belov, K. V., et al. 2020, Nature, 587, 59 [NASA ADS] [CrossRef] [Google Scholar]
- Braginsky, S. I., & Roberts, P. H. 1995, Geophys. Astrophys. Fluid Dyn., 79, 1 [Google Scholar]
- Brandenburg, A. 2008, Astron. Nachr., 329, 725 [NASA ADS] [CrossRef] [Google Scholar]
- Brown, B. P., Vasil, G. M., & Zweibel, E. G. 2012, ApJ, 756, 109 [Google Scholar]
- Bucciantini, N., Metzger, B. D., Thompson, T. A., & Quataert, E. 2012, MNRAS, 419, 1537 [NASA ADS] [CrossRef] [Google Scholar]
- Bugli, M., Guilet, J., Obergaulinger, M., Cerdá-Durán, P., & Aloy, M. A. 2020, MNRAS, 492, 58 [NASA ADS] [CrossRef] [Google Scholar]
- Bugli, M., Guilet, J., & Obergaulinger, M. 2021, MNRAS, 507, 443 [NASA ADS] [CrossRef] [Google Scholar]
- Busse, F. H., & Simitev, R. D. 2006, Geophys. Astrophys. Fluid Dyn., 100, 341 [NASA ADS] [CrossRef] [Google Scholar]
- CHIME/FRB Collaboration (Andersen, B. C., et al.) 2020, Nature, 587, 54 [Google Scholar]
- Christensen, U., & Wicht, J. 2015, in Treatise on Geophysics (Second Edition), ed. G. Schubert (Oxford: Elsevier), 245 [CrossRef] [Google Scholar]
- Ciolfi, R., Kastaun, W., Kalinani, J. V., & Giacomazzo, B. 2019, Phys. Rev. D, 100, 023005 [NASA ADS] [CrossRef] [Google Scholar]
- Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. 2018, MNRAS, 474, 961 [NASA ADS] [CrossRef] [Google Scholar]
- Coti Zelati, F., Borghese, A., Israel, G. L., et al. 2021, ApJ, 907, L34 [NASA ADS] [CrossRef] [Google Scholar]
- Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52 [NASA ADS] [CrossRef] [Google Scholar]
- Deng, H., Mayer, L., Latter, H., Hopkins, P. F., & Bai, X.-N. 2019, ApJS, 241, 26 [NASA ADS] [CrossRef] [Google Scholar]
- Dessart, L., Burrows, A., Livne, E., & Ott, C. D. 2008, ApJ, 673, L43 [NASA ADS] [CrossRef] [Google Scholar]
- Drout, M. R., Soderberg, A. M., Gal-Yam, A., et al. 2011, ApJ, 741, 97 [NASA ADS] [CrossRef] [Google Scholar]
- Duncan, R. C., & Thompson, C. 1992, ApJ, 392, L9 [Google Scholar]
- 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 [NASA ADS] [Google Scholar]
- Fong, W., Laskar, T., Rastinejad, J., et al. 2021, ApJ, 906, 127 [NASA ADS] [CrossRef] [Google Scholar]
- Fromang, S., & Papaloizou, J. 2007, A&A, 476, 1113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661 [NASA ADS] [Google Scholar]
- Gastine, T., & Wicht, J. 2012, Icarus, 219, 428 [Google Scholar]
- Gastine, T., & Wicht, J. 2021, Icarus, 368, 114514 [NASA ADS] [CrossRef] [Google Scholar]
- Gastine, T., Duarte, L., & Wicht, J. 2012, A&A, 546, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Gilman, P. A., & Glatzmaier, G. A. 1981, ApJS, 45, 335 [CrossRef] [Google Scholar]
- Gompertz, B. P., O’Brien, P. T., & Wynn, G. A. 2014, MNRAS, 438, 240 [NASA ADS] [CrossRef] [Google Scholar]
- Götz, D., Mereghetti, S., Molkov, S., et al. 2006, A&A, 445, 313 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Gouhier, B., Lignières, F., & Jouve, L. 2021, A&A, 648, A109 [EDP Sciences] [Google Scholar]
- Gouhier, B., Jouve, L., & Lignières, F. 2022, A&A, 661, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Gressel, O., & Pessah, M. E. 2015, ApJ, 810, 59 [NASA ADS] [CrossRef] [Google Scholar]
- Gubbins, D., & Zhang, K. 1993, Phys. Earth Planet. Int., 75, 225 [CrossRef] [Google Scholar]
- Guilet, J., & Müller, E. 2015, MNRAS, 450, 2153 [NASA ADS] [CrossRef] [Google Scholar]
- Guilet, J., Müller, E., & Janka, H.-T. 2015, MNRAS, 447, 3992 [NASA ADS] [CrossRef] [Google Scholar]
- Guilet, J., Bauswein, A., Just, O., & Janka, H.-T. 2017, MNRAS, 471, 1879 [NASA ADS] [CrossRef] [Google Scholar]
- Hawley, J. F., & Balbus, S. A. 1992, ApJ, 400, 595 [NASA ADS] [CrossRef] [Google Scholar]
- Hüdepohl, L. 2014, PhD Thesis, Technische Universität, München, Germany [Google Scholar]
- Hurley, K., Boggs, S. E., Smith, D. M., et al. 2005, Nature, 434, 1098 [NASA ADS] [CrossRef] [Google Scholar]
- Inserra, C., Smartt, S. J., Jerkstrand, A., et al. 2013, ApJ, 770, 128 [NASA ADS] [CrossRef] [Google Scholar]
- Israel, G. L., Belloni, T., Stella, L., et al. 2005, ApJ, 628, L53 [NASA ADS] [CrossRef] [Google Scholar]
- Jones, C. A., Boronski, P., Brun, A. S., et al. 2011, Icarus, 216, 120 [NASA ADS] [CrossRef] [Google Scholar]
- Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261 [Google Scholar]
- Kiuchi, K., Kyutoku, K., Sekiguchi, Y., & Shibata, M. 2018, Phys. Rev. D, 97, 124039 [NASA ADS] [CrossRef] [Google Scholar]
- Kouveliotou, C., Dieters, S., Strohmayer, T., et al. 1998, Nature, 393, 235 [NASA ADS] [CrossRef] [Google Scholar]
- Krause, F., & Raedler, K. H. 1980, Mean-field Magnetohydrodynamics and Dynamo Theory (Oxford: Pergamon Press) [Google Scholar]
- Kuroda, T., Arcones, A., Takiwaki, T., & Kotake, K. 2020, ApJ, 896, 102 [NASA ADS] [CrossRef] [Google Scholar]
- Lantz, S. R. 1992, PhD Thesis, Cornell University, USA [Google Scholar]
- Lattimer, J. M., & Swesty, D. F. 1991, Nucl. Phys. A, 535, 331 [CrossRef] [Google Scholar]
- Lesur, G., & Longaretti, P. Y. 2007, MNRAS, 378, 1471 [NASA ADS] [CrossRef] [Google Scholar]
- Lesur, G., & Ogilvie, G. I. 2008, A&A, 488, 451 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Longaretti, P. Y., & Lesur, G. 2010, A&A, 516, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Lorimer, D. R., Bailes, M., McLaughlin, M. A., Narkevic, D. J., & Crawford, F. 2007, Science, 318, 777 [Google Scholar]
- Margalit, B., Metzger, B. D., Thompson, T. A., Nicholl, M., & Sukhbold, T. 2018, MNRAS, 475, 2659 [NASA ADS] [CrossRef] [Google Scholar]
- Masada, Y., Sano, T., & Shibata, K. 2007, ApJ, 655, 447 [NASA ADS] [CrossRef] [Google Scholar]
- Masada, Y., Takiwaki, T., Kotake, K., & Sano, T. 2012, ApJ, 759, 110 [NASA ADS] [CrossRef] [Google Scholar]
- Masada, Y., Takiwaki, T., & Kotake, K. 2015, ApJ, 798, L22 [Google Scholar]
- Masada, Y., Takiwaki, T., & Kotake, K. 2022, ApJ, 924, 75 [NASA ADS] [CrossRef] [Google Scholar]
- Meheut, H., Fromang, S., Lesur, G., Joos, M., & Longaretti, P.-Y. 2015, A&A, 579, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Menou, K., Balbus, S. A., & Spruit, H. C. 2004, ApJ, 607, 564 [NASA ADS] [CrossRef] [Google Scholar]
- Metzger, B. D., Quataert, E., & Thompson, T. A. 2008, MNRAS, 385, 1455 [CrossRef] [Google Scholar]
- Metzger, B. D., Giannios, D., Thompson, T. A., Bucciantini, N., & Quataert, E. 2011, MNRAS, 413, 2031 [Google Scholar]
- Metzger, B. D., Beniamini, P., & Giannios, D. 2018a, ApJ, 857, 95 [NASA ADS] [CrossRef] [Google Scholar]
- Metzger, B. D., Thompson, T. A., & Quataert, E. 2018b, ApJ, 856, 101 [NASA ADS] [CrossRef] [Google Scholar]
- Moffatt, H. K. 1978, Cambridge Monographs on Mechanics and Applied Mathematics (Cambridge: Cambridge University Press) [Google Scholar]
- Moiseenko, S. G., Bisnovatyi-Kogan, G. S., & Ardeljan, N. V. 2006, MNRAS, 370, 501 [NASA ADS] [CrossRef] [Google Scholar]
- Mösta, P., Richers, S., Ott, C. D., et al. 2014, ApJ, 785, L29 [CrossRef] [Google Scholar]
- Mösta, P., Ott, C. D., Radice, D., et al. 2015, Nature, 528, 376 [CrossRef] [Google Scholar]
- Mösta, P., Radice, D., Haas, R., Schnetter, E., & Bernuzzi, S. 2020, ApJ, 901, L37 [CrossRef] [Google Scholar]
- Nicholl, M., Smartt, S. J., Jerkstrand, A., et al. 2013, Nature, 502, 346 [NASA ADS] [CrossRef] [Google Scholar]
- Obergaulinger, M., & Aloy, M. Á. 2017, MNRAS, 469, L43 [Google Scholar]
- Obergaulinger, M., & Aloy, M. Á. 2020, MNRAS, 492, 4613 [Google Scholar]
- Obergaulinger, M., Cerdá-Durán, P., Müller, E., & Aloy, M. A. 2009, A&A, 498, 241 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Obergaulinger, M., Just, O., & Aloy, M. A. 2018, J. Phys. G Nucl. Part. Phys., 45, 084001 [NASA ADS] [CrossRef] [Google Scholar]
- Ott, C. D., Burrows, A., Thompson, T. A., Livne, E., & Walder, R. 2006, ApJS, 164, 130 [NASA ADS] [CrossRef] [Google Scholar]
- Parker, E. N. 1955, ApJ, 122, 293 [Google Scholar]
- Pessah, M. E., & Chan, C.-K. 2008, ApJ, 684, 498 [NASA ADS] [CrossRef] [Google Scholar]
- Pessah, M. E., Chan, C.-K., & Psaltis, D. 2007, ApJ, 668, L51 [Google Scholar]
- Rampp, M., & Janka, H. T. 2002, A&A, 396, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Raynaud, R., Guilet, J., Janka, H. T., & Gastine, T. 2020, Sci. Adv., 6, eaay2732 [NASA ADS] [CrossRef] [Google Scholar]
- Raynaud, R., Cerdá-Durán, P., & Guilet, J. 2021, MNRAS, 509, 3410 [NASA ADS] [CrossRef] [Google Scholar]
- Reboul-Salze, A., Guilet, J., Raynaud, R., & Bugli, M. 2021, A&A, 645, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Rembiasz, T., Guilet, J., Obergaulinger, M., et al. 2016, MNRAS, 460, 3316 [NASA ADS] [CrossRef] [Google Scholar]
- Roberts, O. J., Veres, P., Baring, M. G., et al. 2021, Nature, 589, 207 [NASA ADS] [CrossRef] [Google Scholar]
- Rodríguez Castillo, G. A., Israel, G. L., Tiengo, A., et al. 2016, MNRAS, 456, 4145 [CrossRef] [Google Scholar]
- Rowlinson, A., O’Brien, P. T., Metzger, B. D., Tanvir, N. R., & Levan, A. J. 2013, MNRAS, 430, 1061 [NASA ADS] [CrossRef] [Google Scholar]
- Schaeffer, N. 2013, Geochem. Geophys. Geosyst., 14, 751 [NASA ADS] [CrossRef] [Google Scholar]
- Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211 [Google Scholar]
- Shi, J.-M., Stone, J. M., & Huang, C. X. 2016, MNRAS, 456, 2273 [NASA ADS] [CrossRef] [Google Scholar]
- Shibata, M., Liu, Y. T., Shapiro, S. L., & Stephens, B. C. 2006, Phys. Rev. D, 74, 104026 [NASA ADS] [CrossRef] [Google Scholar]
- Shibata, M., Fujibayashi, S., & Sekiguchi, Y. 2021, Phys. Rev. D, 104, 063026 [NASA ADS] [CrossRef] [Google Scholar]
- Siegel, D. M., Ciolfi, R., Harte, A. I., & Rezzolla, L. 2013, Phys. Rev. D, 87, 121302 [NASA ADS] [CrossRef] [Google Scholar]
- Simon, J. B., & Hawley, J. F. 2009, ApJ, 707, 833 [NASA ADS] [CrossRef] [Google Scholar]
- Svinkin, D., Frederiks, D., Hurley, K., et al. 2021, Nature, 589, 211 [NASA ADS] [CrossRef] [Google Scholar]
- Thompson, C., & Duncan, R. C. 1993, ApJ, 408, 194 [NASA ADS] [CrossRef] [Google Scholar]
- Tiengo, A., Esposito, P., Mereghetti, S., et al. 2013, Nature, 500, 312 [CrossRef] [Google Scholar]
- Tilgner, A., & Busse, F. H. 1997, J. Fluid Mech., 332, 359 [NASA ADS] [CrossRef] [Google Scholar]
- Turolla, R., Zane, S., & Watts, A. L. 2015, Rep. Progr. Phys., 78, 116901 [Google Scholar]
- Wicht, J. 2002, Phys. Earth Planet. Inter., 132, 281 [Google Scholar]
- Winteler, C., Käppeli, R., Perego, A., et al. 2012, ApJ, 750, L22 [Google Scholar]
- Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015 [NASA ADS] [CrossRef] [Google Scholar]
- Xue, Y. Q., Zheng, X. C., Li, Y., et al. 2019, Nature, 568, 198 [NASA ADS] [CrossRef] [Google Scholar]
- Yoshimura, H. 1975, ApJ, 201, 740 [Google Scholar]
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.
![]() |
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. |
Appendix B: Parameters of our global simulations
Overview of the global models. All simulations use Pm = 16 and insulating magnetic boundary conditions. The (nr, nθ, nϕ) = (257, 512, 1024) resolution is called “Medium” and the (nr, nθ, nϕ) = (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, Bo to the initial magnetic field strength at the outer boundary. Magnetic field and dipole field strengths are time and space averaged. Psim corresponds to the estimated period of the mean-field dynamo patterns at r = 0.92 ro, kz to the vertical wavenumber and qΩ to the value in the turbulent region (at r = 0.92 ro). kz, qΩ and the alpha component αϕϕ are used to compute the theoretical αΩ period PαΩ.
Appendix C: Correlation coefficients between the EMF and the magnetic field
![]() |
Fig. C.1. Time-averaged values of the correlation coefficient between the EMF and the magnetic field taken at r = 0.86 ro. |
Appendix D: Estimation of the components of the α tensor
![]() |
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 ro. |
Appendix E: Azimuthal average of magnetic and velocity fields at quasi-stationary state
E.1. Simulation Standard
![]() |
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 |
E.2. Simulation Boussinesq
![]() |
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 |
E.3. Simulation Pr=0.1
![]() |
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 |
All Tables
Periods of the different mean-field patterns for Pr = 0.005, Pr = 0.02, and Pr = 0.1.
Overview of the global models. All simulations use Pm = 16 and insulating magnetic boundary conditions. The (nr, nθ, nϕ) = (257, 512, 1024) resolution is called “Medium” and the (nr, nθ, nϕ) = (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, Bo to the initial magnetic field strength at the outer boundary. Magnetic field and dipole field strengths are time and space averaged. Psim corresponds to the estimated period of the mean-field dynamo patterns at r = 0.92 ro, kz to the vertical wavenumber and qΩ to the value in the turbulent region (at r = 0.92 ro). kz, qΩ and the alpha component αϕϕ are used to compute the theoretical αΩ period PαΩ.
All Figures
![]() |
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. |
In the text |
![]() |
Fig. 2. Ratio of the squared Brunt–Väisälä frequency N2 (Eq. (16)) to the squared initial rotation frequency Ω2 (Eq. (10)). The buoyancy influence is strongest in the red regions. |
In the text |
![]() |
Fig. 3. Meridional slices at ϕ = 0 of the initial poloidal magnetic field Br (left) and of the toroidal magnetic field Bϕ at t = 773 ms (right) for model Standard. |
In the text |
![]() |
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. |
In the text |
![]() |
Fig. 5. Time evolution of the magnetic and turbulent kinetic energy densities for model Standard with Bo = 1.5 × 1014 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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
Fig. 8. Radius-time diagram of the dipolar modes. Top: radius-time diagram of the axial dipolar mode Br(l = 1, m = 0). The period of the dipole reversal is Pdip ≃ 410 ms. Bottom: radius-time diagram of the equatorial dipolar mode amplitude Br(l = 1, m = 1) at a given azimuth ϕ = 0 in the rotating frame of the simulation. The fast oscillations are due to the rotation of this dipole mode. |
In the text |
![]() |
Fig. 9. Butterfly diagrams of |
In the text |
![]() |
Fig. 10. Time-averaged correlation coefficients (Eq. (29)) between ℰϕ and |
In the text |
![]() |
Fig. 11. Time-averaged values of the αϕϕ component estimated with Eq. (30), (blue) and with the turbulent kinetic helicity (orange), (Eq. (31)) taken at r = 0.86 ro and with τc = 2.5 ms. |
In the text |
![]() |
Fig. 12. Meridional cuts of Bϕ (top) and azimuthal average of Ω (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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
Fig. 15. Butterfly diagrams of |
In the text |
![]() |
Fig. 16. Radial velocity ur and αϕϕ for anelastic models with buoyancy and for increasing Prandtl numbers. Top: snapshots of the radial velocity ur. 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 αϕϕ for Pr = 0.005 (left), Pr = 0.02 (middle), and Pr = 0.1 (right). The dashed line represents the spherical radius r = 0.92 ro. |
In the text |
![]() |
Fig. 17. Normalized volume-averaged spectrum of the radial velocity squared ur 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 lc ≈ R/Lc above which the scales are impacted by buoyancy. |
In the text |
![]() |
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. |
In the text |
![]() |
Fig. 19. Time-averaged values of the αϕϕ component estimated using Eq. (30) at r = 0.92 ro in the turbulent region for three different Prandtl numbers: Pr = 0.005 (blue), Pr = 0.02 (orange), and Pr = 0.1 (green). |
In the text |
![]() |
Fig. 20. Dimensionless dipole strength ℬdipole as a function of the dimensionless magnetic field strength ℬ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 ℬ were used in order to compare results with different densities |
In the text |
![]() |
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. |
In the text |
![]() |
Fig. C.1. Time-averaged values of the correlation coefficient between the EMF and the magnetic field taken at r = 0.86 ro. |
In the text |
![]() |
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 ro. |
In the text |
![]() |
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 |
In the text |
![]() |
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 |
In the text |
![]() |
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 |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.