© The Authors 2025

1. Introduction

Beginning with the discovery of sunspot magnetic fields by George Ellery Hale in 1908, the magnetic activity of the Sun and solar-like stars has captivated astronomers for over a century. These magnetic phenomena, which include the 11-year solar cycle, auroras, and solar flares, arise from the dynamo mechanism within stellar convection zones (Moffatt 1978; Dormy & Soward 2007; Charbonneau 2020; Brun & Browning 2017). In such dynamos, the interaction of turbulent convection, differential rotation, and magnetic field generation creates a dynamic and self-sustaining system of magnetic energy with a distribution depending on the astrophysical objects under consideration. Extending our understanding to other stars of similar types reveals both universal principles and fascinating differences in magnetic behaviors (Donati & Landstreet 2009; Jeffers et al. 2023).

Solar-like dynamos are particularly challenging to model due to the vast range of scales and physical effects involved. From the millisecond timescales of Alfvén waves to the billion-year lifespans of main-sequence stars, a comprehensive model must reconcile processes occurring at vastly different temporal and spatial resolutions. Observational data such as helioseismic profiles provide valuable constraints for the dynamo models but simultaneously reveal the intricacies of phenomena such as differential rotation and meridional circulation that need to be reproduced. Numerical simulations have made significant strides in capturing the broad features of solar-like dynamos, yet they remain limited by computational resources and the inherent complexity of the system (Käpylä et al. 2023). In general, dynamo action requires a complex three-dimensional flow structure in order to bypass anti-dynamo theorems (Cowling 1933). Thus, in addition to temporal challenges, modeling stellar magnetism involves bridging length scales ranging from the millimeter-scale eddies of convection to the global-scale flows spanning the entire stellar radius and influenced by rotation and the Coriolis force. The equations governing these magnetohydrodynamic (MHD) systems include nonlinear interactions between these scales that are sensitive to initial and boundary conditions. Dimensionless parameters (such as the Rossby number) quantifying the rotational influence on convective turbulence help categorize different dynamo behaviors with respect to the flow regimes (Christensen & Aubert 2006; Schrinner et al. 2012; Gastine et al. 2012; Raynaud et al. 2015; Menu et al. 2020). However, they emphasize the transitions between flow regimes that are still not well understood in the context of stellar physics, such as from solar-like differential rotation with fast equator and slow poles to inverse anti-solar rotation (Gastine et al. 2014; Karak et al. 2015). These large-scale zonal flows, generated by convective motions in rotating spherical shells (Christensen 2001), are important for the dynamo process, as they shear the poloidal magnetic field into toroidal field lines (the so-called Ω-effect).

Despite these challenges, advances in observational techniques, such as spectropolarimetric Zeeman-Doppler imaging (ZDI) (Donati & Landstreet 2009), asteroseismology, and chromospheric Ca II H and K emission lines (Böhm-Vitense 2007), have revealed active and inactive branches of magnetic activity cycles across different types of stars (Brandenburg et al. 2017). Observations have shown a negative correlation between the ratio of the magnetic cycle to rotation periods and the Rossby number, at least on the inactive branch, and the existence of modulated multi-period cycles, which have been also reported numerically (Strugarek et al. 2018; Brun et al. 2022). However, transitions between these activity branches, as well as the scaling of the cycle period with stellar rotation, and the corresponding dynamo mechanisms are still debated (Olspert et al. 2018; Saikia et al. 2018b). Further discrepancies remain, including the inability of simulations to fully replicate the Sun’s convective amplitudes or equatorward propagation of dynamo waves with a preferred dipolar symmetry and a complex temporal behavior on long timescales (Raynaud & Tobias 2016). The magnetic fields resulting from Boussinesq (Busse & Simitev 2006; Goudard & Dormy 2008; Schrinner et al. 2011, 2012) or anelastic models (Gastine et al. 2012; Schrinner et al. 2014; Strugarek et al. 2017, 2018; Pinçon et al. 2024) exhibit cyclic magnetic behaviors with low magnetic Reynolds numbers and a parity that is often quadrupolar, whereas the Sun’s magnetic field is dominated by its dipolar components that have periodic oscillations with very high magnetic Reynolds numbers.

Surface stellar magnetic fields, reconstructed from the inversion of multiple polarized spectral lines over several stellar rotation periods with ZDI, are of particular interest, as they simultaneously infer activity and magnetic topology of stellar dynamos. Examples of stars exhibiting solar-like magnetic cycles with predominantly dipolar fields during low levels of activity include 61 Cyg, ϵ Eri, and κ Cet, with activity periods of 7 to 13 years (Jeffers et al. 2023), and F-type stars τ Boo and HD 75332, with shorter cycles of 1 to 3 years (Jeffers et al. 2018; Brown et al. 2021). On the other hand, younger K stars such as EK Dra and HN Peg have rapidly evolving magnetic fields without correlation with an S-index. Unexpected potentially acyclic reversals of magnetic topology or minima in activity have also been observed, for example, in LQ Hya, HD 16620 (Willamo et al. 2022; Lehtinen et al. 2022). During magnetic reversals, the dominant magnetic topology switches to being multipolar or quadrupolar (Saikia et al. 2018a). These observations suggest that stellar cycles can become intermittent and aperiodic at certain stages of stellar evolution (Jeffers et al. 2023; Van Saders et al. 2016; Metcalfe et al. 2022). Although surface magnetic fields are influenced by small-scale turbulence, their large-scale components must evolve with the global nature of the dynamo action, and this coupling represents an important motivation for the current study.

In this work, we focus on the dynamical mechanisms of large-scale magnetic behavior in strongly stratified and turbulent models of low-mass stars and show that density stratification can promote dipolar, albeit aperiodic dynamo states in such systems. Such stratification separates the flow into an inner layer of rotationally constrained convective columns and an outer layer of anisotropic turbulent small-scale convection, thus promoting magnetic pumping of small-scale magnetic fluctuations inside the tangent cylinder (i.e., the imaginary cylindrical surface around the inner boundary of our computational domain). This process enhances the dipolar magnetic mode, which competes with a more complex multipolar magnetic state and dynamically interacts with differential rotation, and thus prevents the transition to anti-solar rotation.

This paper is structured as follows. In Sect. 2, we briefly discuss the anelastic model and our choice of parameter regime. In Sect. 3, we describe the dynamo states observed in our simulations. Section 4 focuses on the dipole solutions that manifest themselves in strongly stratified convection, and Sect. 5 presents the argument for magnetic pumping supporting dipolar solutions. In Sect. 6, we discuss of our results in the context of previous work and stellar observations, and we conclude with our findings in Sect. 7.

2. Anelastic simulations

2.1. Setup

We perform direct numerical simulations of the dynamo process in anelastic approximation using MagIC code (Gastine et al. 2014), solving equations for velocity, magnetic fields, and entropy in spherical coordinates (r, θ, ϕ), as described in Appendix A. Cylindrical coordinates (r_c, ϕ, z_c) are also used where necessary. Conducting fluid is contained in a shell between two spheres with radius ratio χ = r_i/r_o = 0.35, co-rotating with angular velocity Ω and representing a thick convective envelope of a solar-like star. Convection is triggered by the entropy contrast ΔS fixed at the boundaries of the spherical shell, under assumption of adiabatic background state. Non-uniform density between the two spheres, modeling density decrease with radius in stellar envelopes, is set by imposing background density profile with polytropic index n = 2. The key dimensionless parameters explored in our work are Rayleigh number Ra and density contrast N_ρ, defined as

$$\begin{array}{c}\hfill Ra=\frac{GMd\mathrm{\Delta }S}{\nu \kappa c_p},N_{\rho }=ln\frac{{\rho }_i}{{\rho }_o},\end{array}$$(1)

where ν is hydrodynamic viscosity, κ thermal diffusivity, c_p heat capacity, d = r_o − r_i the gap between the two spheres, G gravitational constant, and M the central mass. Rayleigh number determines the strength of convection relative to diffusive transport; convection is driven at a certain critical value Ra_cr. This value increases with density contrast (Jones et al. 2009; Raynaud et al. 2015), and so simulations at high N_ρ become challenging. To approximate stellar conditions, it is necessary to set

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathit{Ra}}=Ra/R{a}_{\mathrm{cr}}\gg 1,{\rho }_i/{\rho }_o\gg 1.\end{array}$$(2)

The rest of dimensionless parameters governing the system are Ekman number E = ν/Ωd², Prandtl number Pr = ν/κ, and magnetic Prandtl number Pm = ν/η, where η is magnetic diffusivity. We fixed E = 3 × 10⁻⁴, corresponding to moderate rotation as compared to viscosity, Pr = 1 and Pm = 2. When compared to estimations of these parameters in the stars (e.g., Ra ∼ 10²⁰, E ∼ 10⁻¹⁵, Pr ∼ 10⁻⁷, Pm ∼ 10⁻⁵ for the Sun), these values in simulations can be interpreted as if additional turbulent dissipation was present in the simulations. We note that $\stackrel{\sim }{\mathit{Ra}}$ is calculated based on the entropy difference across the convective zone, which is hard to measure in stellar or planetary context, and includes dissipative effects in its definition, which are considered weak in this context. A modified flux-based Rayleigh number Ra_F^* = RaNuEk³/Pr², eliminating dissipative effects and based on the heat flux that is easier to estimate from observations, has been proposed by Christensen (2002) and Aurnou et al. (2020) as a scaling of convection strength. The Nusselt number Nu is defined as the output luminosity to the basic state luminosity (see Appendix B). Our data show that flux-based Rayleigh number scales linearly with $\stackrel{\sim }{\mathit{Ra}}$ (Table B.1); for simplicity, we keep $\stackrel{\sim }{\mathit{Ra}}$ as a parameter in this work.

A typical simulation in our work was set up as follows. As an initial condition for magnetic, velocity and temperature fields, we used a steady-state flow snapshots from the convective dynamo database of Pinçon et al. (2024). Furthermore, we employ stress-free boundary conditions for velocity field and insulating boundary condition for magnetic field. Typical resolution in our simulations span from N_r × N_ϕ = [65, 576] to [241, 1280] where N_r denotes the number of Chebyshev polynomials in radial direction and N_ϕ the total number of spherical harmonics. The resolution was set to have at least three order of magnitude difference between the energy in the largest and the smallest scales in the energy spectra for each combination of N_ρ and $\stackrel{\sim }{\mathit{Ra}}$, and the simulations were run at high resolution until a steady state was obtained again. After that, the number of spherical harmonics was reduced to speed up the calculation and gather long-term flow statistics and snapshots. For more details on the parameters, flow equations and numerical methods, see Appendices A, B, and the website of MagIC code¹.

2.2. Choice of parameter regime

In this section, we explore N_ρ = (2, 4, 6), corresponding to density ratios of ρ_i/ρ_o = (7.4, 54.6, 403.4) between the inner and the outer spheres. We also provide a motivation for the choice of parameter regimes for our dynamo study.

First, we compared the flow structures for these three values of N_ρ at the same level of turbulence, $\stackrel{\sim }{\mathit{Ra}}=32$. In rotating convection, the Coriolis force tends to remove vertical gradients and align flow structures along the rotation axis. When convective turbulence sets in, inertia forces due to nonlinear interactions of flow fluctuations are in competition with the Coriolis force and tend to break this alignment. The result is a rather three-dimensional structure of the flow, as is visible from the surface distribution of the vertical component of flow vorticity in the bulk of the domain, r/r_o = 0.5 for N_ρ = 2 (Figure 1a). Since the density contrast is low in this case, the structure of convection does not vary considerably along the radius, and thus the distribution of vorticity across different length scales is similar both in the bulk of the domain, r/r_o = 0.5, and near the surface, r/r_o = 0.9 (Figure 1a,d). However, simultaneously increasing convection and stratification has a profound effect on turbulent flow properties: For the same level of convection, its length scale decreases toward the surface for N_ρ = 4 (Figure 1b,e). Further increase of stratification to N_ρ = 6 results in extremely strong dynamical contrast between the surface and the bulk, where large-scale, rotationally constrained convective columns develop at r/r_o = 0.5 (Figure 1c), and small-scale, three-dimensional vortical structures at the surface (Figure 1f). The flow in the domain is essentially separated in a Coriolis-affected interior and a strongly turbulent surface layer. The influence of surface layers on dynamics was often overlooked in the previous works, although such layers possess strong kinetic helicity and are thus natural candidates to generate magnetic fields.

Fig. 1. Instantaneous surfaces of the vertical vorticity component of the flow at r/r0 = 0.5 (top) and r/r0 = 0.9 (bottom) illustrating scale separation developing with Nρ in strongly turbulent convection, $\stackrel{\sim }{\mathit{Ra}}=32$. (a,d) Nρ = 2; (b,e) Nρ = 4; (c,f) Nρ = 6.

Quantitatively, the separation of scales in stratified convection can be described by the local convective Rossby number, representing the balance between inertia and Coriolis forces at each radii,

$$\begin{array}{c}\hfill R{o}_l(r)=\frac{u_{\mathrm{rms}}(r)}{l_{\mathrm{conv}}(r)\mathrm{\Omega }},l_{\mathrm{conv}}=\frac{\pi E^{\mathrm{kin}}(r)}{{\displaystyle {\sum }_l}l{E}_l^{\mathrm{kin}}(r)},m\ne 0,\end{array}$$(3)

where u_rms is the root-mean-square convective velocity and l_conv is the typical length scale of convection calculated from the distribution of kinetic energy E^kin with spherical harmonic degree l (Christensen & Aubert 2006), given in the units of d. The energy contained in the axisymmetric, azimuthal zonal flow ⟨u_ϕ⟩_ϕ, averaged over longitude ϕ, is subtracted to separate proper convective motions from differential rotation driven by them. For stars, it is not straightforward to estimate this local parameter, and so global Rossby number, the ratio of stellar rotation period to empirically inferred convective turnover time, is used instead (Noyes et al. 1984).

Figure 2a presents the dependence of the local Ro number with radius for varying N_ρ ($\stackrel{\sim }{\mathit{Ra}}=32$, solid). For N_ρ = 2, it is relatively large, Ro_l ∼ 0.5, and is nearly constant across the spherical shell. It increases proportionally in the entire domain with increasing $\stackrel{\sim }{\mathit{Ra}}$ and enhancement of convection. Increasing density stratification consistently decreases Ro_l in the deep interiors, developing there a rotationally constrained zone Ro_l ≪ 1; for N_ρ = 6, it occupies about a half of the entire simulation domain. There, the dominance of the Coriolis force enforce Taylor-Proudman constraint to produce elongated convective columns, as seen in Figure 1c; effectively, the flow in the bulk becomes less dependent on the vertical coordinate z. On the other hand, the radial profile of Ro_l becomes much steeper throughout the domain if compared to N_ρ = 2. The area with Ro_l > 1, developing for N_ρ = 4 and 6, corresponds to the strongly turbulent surface layer with inertia-dominated small-scale vortices (Figure 1f). The increase in Ro_l at the surface with Ra (Figure 2a) is accompanied by a decrease of the length scale of convection (not shown here), meaning that local l_conv(r) at the surface becomes about 5 times smaller than in the bulk. This essentially creates a scale separation between convection in the bulk and the surface layer, which can be observed in Figure 1(c,f), with large scales much more energetic in the bulk, and small scales at the surface. To develop solutions with such pronounced scale separation, high values $\stackrel{\sim }{\mathit{Ra}}$ are needed. Figure 2a shows that Ro_l gradually increases with increase of $\stackrel{\sim }{\mathit{Ra}}$ at the outer shell (dashed lines, fixed N_ρ = 6), developing a surface area where Ro > 1. Thus, it is the simultaneous increase of density stratification and enhancement of convection that results in Coriolis force dominating in the bulk, and turbulent inertia dominating the outer layer.

Fig. 2. (a) Effect of varying Nρ (in red, Ra/Racr = 32) and Ra (in blue, Nρ = 6) on the local Rossby number calculated based on the convective length scale lconv (3). Dotted line denotes Rol = 1, e.g., a balance between inertia and the Coriolis force. (b) Density contrasts Nρ and Rayleigh numbers $\stackrel{\sim }{\mathit{Ra}}$ spanned by simulations in this work. Different symbols denote three dynamical regimes of the flow, discussed in more details in Sect. 3. Filled (empty) triangles, predominantly multipolar fields with solar (anti-solar) differential rotation (Sect. 3.2). Circles, quadrupolar oscillatory dynamos with solar-like differential rotation (Sect. 3.3). Diamonds, reversing dipole-multipole solutions with transitions between three-layered and anti-solar rotation (Sect. 3.4). Solar- and anti-solar differential rotation is shown in the inset plots depicting zonal flows ⟨uϕ⟩ϕ for Nρ = 2, $\stackrel{\sim }{\mathit{Ra}}=32$ (left) and $\stackrel{\sim }{\mathit{Ra}}=96$ (right).

Previous studies of non-stratified convection showed that convective turbulence with strong inertia and Ro_l > 1 develops multipolar, disorganized dynamo solutions (Christensen & Aubert 2006). Similar behavior was found for weakly stratified anelastic convection Gastine et al. (2012), Raynaud et al. (2015), who found that the stability range of dipolar dynamos shrinks in $\stackrel{\sim }{\mathit{Ra}}$ with increase of N_ρ. In this work, we are interested in the joint effect of stratification and turbulence, and scale separation, on the dynamo properties. We thus explore parameter space around solutions with N_ρ = 4, 6 and $\stackrel{\sim }{\mathit{Ra}}=32$, more representative of stellar flows. The entire set of our simulations, all in fully turbulent convective regimes, is schematically presented in Figure 2b. Due to numerical constraints, we probe a wider range of turbulent states for N_ρ = 2, with $\stackrel{\sim }{\mathit{Ra}}\in [32,96]$, than for N_ρ = 6 ($\stackrel{\sim }{\mathit{Ra}}\in [8,32]$).

An important feature of this parameter space is transition from solar-like differential rotation, with faster equator and slower poles to anti-solar differential rotation with slower equator and faster poles (see inset plots in Figure 2b). For N_ρ = 2, this transition takes place between $\stackrel{\sim }{\mathit{Ra}}=72$ and $\stackrel{\sim }{\mathit{Ra}}=80$; our most stratified flows with $\stackrel{\sim }{\mathit{Ra}}$ appear to be near the threshold of such a transition. This affects the properties of dynamo and convection; for example, the structure of convection is modified by differential rotation for N_ρ = 6 and $\stackrel{\sim }{\mathit{Ra}}$, which creates a local increase in convective length scale near r_o and the corresponding local deficiency in Ro_l (Figure 2a), comparable to the thickness of weak retrograde zonal jet at the surface (see Sect. 5.4 and Figure 5c). We note that all simulations discussed below are magnetohydrodynamic dynamos; in the absence of feedback of magnetic stresses on the flow, transition to anti-solar rotation takes place at lower $\stackrel{\sim }{\mathit{Ra}}$.

3. Types of dynamo states in simulations

In this section, we examine the three different types of dynamo solutions observed in the flow for the set of simulations in Figure 2b. While all of our simulations are fully turbulent dynamos with energy spread across a range of scales, their dynamics is different if the large scales are considered. To illustrate these differences, we use the following characteristics: butterfly diagrams of the radial component of magnetic field B_r, filtered snapshots of the large-scale magnetic field, and zonal flows. In the following, we discuss in details these diagnostic quantities and resulting dynamo states.

3.1. Diagnostic quantities

Butterfly diagrams in Figure 3 present the dependence of the azimuthally averaged radial component of the magnetic field ⟨B_r⟩_ϕ at the surface of the outer sphere, r = r_o, as a function of rotational time Ωt and latitude θ. They serve as a proxy for both magnetic topology – dipolar or multipolar, and dynamics. A more quantitative representation of the dynamo topology is the surface dipolarity,

Fig. 3. Butterfly diagrams of the axisymmetric, azimuthally averaged magnetic component ⟨Br⟩ϕ at the outer boundary. (a) Nρ = 2, Ra/Racr = 32; (b) Nρ = 6, Ra/Racr = 16; and (c) Nρ = 6, Ra/Racr = 32. To highlight large-scale components, the color maps were saturated with cutoffs of 0.2, 0.25, and 0.1, respectively; the magnitude of magnetic fluctuations is larger by a corresponding factor. Black curves correspond to the surface dipolarity with cutoff parameter of lcut = 11 (4).

$$\begin{array}{c}\hfill f_{\mathrm{dip}}^{{\ell }_{\mathrm{cut}}}={\left(\frac{{\displaystyle {\sum }_{m=-1}^1}{E}_{\mathrm{mag}}^{l=1,m}(r_o)}{{\displaystyle {\sum }_{\ell =1}^{{\ell }_{\mathrm{cut}}}{\sum }_{m=-\ell }^{m=+\ell }}{E}_{\mathrm{mag}}^{l,m}(r_o)}\right)}^{1/2},\end{array}$$(4)

which is the square root of the ratio between the energy contained in the dipolar spherical harmonic with m ∈ [ − 1, 0, 1], l = 1, and the energy contained in the first l_cut spherical harmonics. The classic definition of a dipolar dynamo is f_dip > 0.5; if f_dip < 0.5, then the dynamo is considered multipolar. We use here l_cut = 11 as a default output of MagIC code; we found that decreasing the cutoff length scale in (4) to l_cut = 7, perhaps more realistic for spectropolarimetric observations and ZDI, does not strongly affect the estimation of f_dip^l_cut. Together with ⟨B_r⟩_ϕ at the surface, Figure 3 shows this quantity for different solution types. The color maps were saturated to highlight large scale fields; the maximum values of magnetic fluctuations close to the surface are actually much stronger in all our simulations.

Spectropolarimetric observations typically allow one to reproduce magnetic field topology up to 5 − 7 spherical harmonics reliably due to mutual cancellation of small-scale magnetic features (Kochukhov 2021). To provide a qualitative comparison of our models to these observations, we performed spectral filtering of the surface radial magnetic field, B_r, with l_cut = 5 as a cutoff parameter. Figure 4 compares these large-scale filtered fields with the corresponding instantaneous maps of the radial magnetic field near the surface, r/r_o = 0.95.

Fig. 4. Instantaneous snapshots of the surface radial magnetic field Br at r/ro = 0.95 (top) with the same fields filtered to l ≤ 5 (bottom) for different parameters: (a,d) Nρ = 6, Ra/Racr = 32, multipolar fields; (b,e) Nρ = 6, Ra/Racr = 16, quadupolar waves; and (c,f) Nρ = 6, Ra/Racr = 32, during dipole period. The color maps in panels (a-c) were saturated to 0.5|Br|max to highlight the spatial structure of magnetic features.

Finally, we identify the structure of differential rotation by the axisymmetric component of the toroidal velocity field ⟨u_ϕ⟩_ϕ, with brackets ⟨⋯⟩_ϕ denoting a spatial average of ϕ direction. In the following, we refer to this quantity as zonal flows. The local sign of the zonal flow is directly related to the mean deviation ⟨u_ϕ⟩_ϕ/r_c of the flow from the solid-body rotation Ω. The zonal flows for different dynamo solutions in our simulations are presented in Figure 5, where the corresponding fields were also averaged in time.

Fig. 5. Zonal flows (in color) for different types of dynamo solutions corresponding to Figure 3. (a) Nρ = 2, $\stackrel{\sim }{\mathit{Ra}}=32$; (b) Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=16$; and (c) Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$, during the low magnetic energy interval with predominance of the dipole on large scales; and (d) Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$, during the high-energy interval without the dipole. Flux lines corresponds to the poloidal component of Poynting flux (see Sect. 5.4). Coordinates are denoted in the cylindrical coordinate system (rc, zc). Fluxes and zonal flows were averaged in time over multiple snapshots.

3.2. Oscillatory dynamos in weakly stratified systems

The first type of dynamo behavior develops at low stratification, N_ρ = 2, and corresponds to the less rotationally constrained runs with higher values of Ro_l throughout the domain (green triangles in Figure 2b). It is a predominantly multipolar solution with a signature of distorted waves propagating toward the poles (Figure 3a). These waves, however, are very weak compared to the overall turbulent background (see Appendix B and Table B.1 therein). In this regime, the surface dipolarity parameter remains always smaller than 0.5, so this flow state lacks any prolonged dipole periods. Consistent with this, the instantaneous snapshots of B_r reveal very small-scale fields from which it is hard to predict the distribution of their large scales (Figure 4a). When filtered, it results in a multipolar configuration (Figure 4d) on large scales, corresponding to low stratification. Such periodic dynamo solutions have been extensively studied in weakly stratified dynamo models, examples including Schrinner et al. (2011, 2012), Karak et al. (2015), Strugarek et al. (2017), Viviani et al. (2019). They are usually attributed to Parker-Yoshimura dynamo waves driven by the mean-field α²Ω or αΩ dynamo mechanisms; the direction of their propagation is defined by the gradient of angular velocity. For this reason, here we do not focus on the particular properties of these waves. See the excellent reviews by Brun & Browning (2017) and Charbonneau (2020), and the references therein for more details.

The differential rotation is solar-like at low $\stackrel{\sim }{\mathit{Ra}}$ in this weakly stratified case. Figure 5a shows the resulting profile of ⟨u_ϕ⟩_ϕ, averaged in time, for $\stackrel{\sim }{\mathit{Ra}}=32$ (in color). Two opposite mean-flow directions develop, a prograde jet at the equator and a retrograde jet in the bulk. These flow structures are aligned vertically, and their relative strengths are comparable. Even though the flow transitions to anti-solar rotation for $\stackrel{\sim }{\mathit{Ra}}\ge 80$, with inversed meridional distribution (see the inset on the right in Figure 2b), the average dipolarity of the flow does not increase with $\stackrel{\sim }{\mathit{Ra}}$ (Figure 8a) and the dynamo configuration remains non-dipolar.

3.3. Oscillatory dynamos in rapidly rotating and strongly stratified systems (large N_ρ and Ro_l < 1 for all radii)

The second type of the dynamo solutions develops when stratification is increased to N_ρ = 4 or 6, but the levels of turbulence are weaker, e.g., $\stackrel{\sim }{\mathit{Ra}}\le 16$ (see the corresponding values of Reynolds number in Table B.1, increasing with the level of supercriticality). These are rather strongly rotationally constrained simulations with Ro_l < 1 throughout the domain (Figure 2a), corresponding to blue circles in Figure 2b. The dynamo is dominated in this case by large-scale, clearly manifested dynamo waves (Figure 3b). These waves are predominantly symmetric with respect to equator; in the following, we refer to such symmetry of magnetic field as quadrupolar. The timescale of these waves is much slower than in previous, weakly stratified case (N_ρ = 2), with about T/T_rot = 500 for wave period as opposed to at most T/T_rot = 100 for weakly stratified flows. Compared to the dynamo configuration in Figure 3a, the waves are much stronger and so their spatiotemporal coherency is visible from magnetic fluctuations at the surface 4b. When filtered from small-scales, large-scale magnetic field in this case represent a quadrupolar, axisymmetric field (Figure 4e), with surface distribution reflecting the moment in time with respect to the total period of the wave (Figure 3e). Contrary to predominantly multipolar solution in Figure 3a, in the quadrupolar waves state all the large-scale modes exhibit nearly periodic behavior, as is reflected in the oscillation of dipolar component. The dipolarity in this case is also quite small and does not exceed $f_{\mathrm{dip}}^{11}\le 0.3$ during its maxima.

Such oscillatory dynamos also exhibit solar-like differential rotation (Figure 5b), although the topology of the zonal flow is different. In this case, the prograde jet has a conical shape and is confined predominantly to the area close to the equator; it exhibits periodic oscillations, i.e., expands and shrinks during magnetic cycle in Figure 3b. The opposing, retrograde flow in the bulk is much weaker than in the previous, low-stratified case; the overall magnitude of the zonal flow also remains weak.

3.4. Temporal variations of the magnetic field symmetry in strongly stratified and turbulent systems

Finally, we describe the third type of dynamo solutions that develops at high levels of turbulence, $\stackrel{\sim }{\mathit{Ra}}=32$ and high stratification, N_ρ = 4 or 6, featuring surface layers with Ro_l > 1 (Figure 2a). These simulations are thus the most representative of stellar flow regimes, especially for slowly rotating stars, and corresponds to the diamond symbols in Figure 2b. It represents an aperiodic dynamo with time intervals of predominantly dipolar configuration (i.e., anti-symmetric with respect to the equator) and intervals of multipolar magnetic solutions without visible domination of any particular spherical harmonic (Figure 3c). During strong dipolar events, the surface dipolarity is large and reaches the values of $f_{\mathrm{dip}}^{11}=0.5-0.6$, formally corresponding to dipolar magnetic fields. The dipolar component is the most intense at high latitudes, and less visible near the equator, where magnetic behavior is less coherent.

The large-scale, dipole component exhibits active dynamics, with two subsequent magnetic reversals around Ωt = 3600, then loss of coherent magnetic field and development of the multipolar magnetic state between Ωt ∈ [4500, 5500], a subsequent dipole re-appearance at about Ωt = 5500, reversal, and loss at Ωt = 6000. We ran these simulations for at least one-two magnetic diffusion timescale, and we observed that the dynamics persist in time. A similar behavior takes place at the lower density stratification, N_ρ = 4, and strong convection $\stackrel{\sim }{\mathit{Ra}}=32$.

When filtered, magnetic fields show predominantly dipolar configuration at N_ρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$ during a dipolar interval (Figure 4f). During non-coherent, multipolar intervals, or dipole reversals, the field topology appears similar to the one for N_ρ = 2 in panel (d). Compared to the total field magnitude in the upper panels, unfiltered fields are order of magnitude stronger for all the three cases.

In this case, two different types of differential rotation are observed at different times in the same simulations. During dipolar periods, a thin retrograde jet develops at the surface, and the prograde jet shifts to the mid-gap between the spheres. The inner flow area, including the area inside the tangent cylinder, develops a very weak retrograde rotation and remains quiescent overall (Figure 5c). In the following, we refer to this configuration as three-layer differential rotation. During multipolar periods, the dynamo entirely transitions to stronger anti-solar differential rotation, with a retrograde jet expanding at the surface and prograde jet shifting into the bulk of the domain, toward the tangent cylinder (Figure 5d). The peak of the shear at the boundary between the surface retrograde jet and inner prograde jet thus also shifts in the interior of the domain.

We note here that these aperiodic dipolar solutions exhibit bi-stability with respect to periodic waves from Sect. 3.3 at N_ρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$, as shown in Figure 3(b), depending on initial conditions. Slightly increasing Ra would potentially eliminate this bi-stability.

4. Dynamical properties of the flow and the field in strongly stratified and turbulent systems

4.1. Predator-prey interaction of zonal flows and magnetic fields

In this section, we focus our attention on the dynamics of strongly stratified and turbulent aperiodic dynamo states, described in the previous Sect. 3.4. To monitor transitions between the two types of differential rotation in this run, three-layer and anti-solar, we follow the radial position of u_ϕ^max at the equator, θ = π/2, i.e., the maximum in the prograde jet. Since both configurations of zonal flows are aligned with the rotation axis (Figure 5c,d), u_ϕ^max is a good proxy for the type of the zonal flow: r/r_o ∼ 0.4 for anti-solar and r/r_o ∼ 0.65 for three-layer configuration. We plot this quantity in Figure 6a, together with the total magnetic energy of the flow, $E_{\mathrm{mag}}^{\mathrm{tot}}$. It appears that the total magnetic energy correlates well with the type of the zonal flow, developing strong peaks when rotation is anti-solar, and local minima otherwise. Therefore, these dynamo regimes exhibit a bi-stability between the flow state with low magnetic energy and three-layer differential rotation, and the flow state with high magnetic energy and two-layer, anti-solar zonal flow. The dynamo transitions between the two states aperiodically, spending in each state about 500 − 1000 rotation times. However, these runs are computationally too expensive to gather accurate temporal statistics for these stochastically driven transitions.

Fig. 6. (a) Total magnetic energy as a function of time (in black) versus the radius of ⟨uϕ⟩ϕmax in the equatorial plane (θ = π/2), in blue. The large radii correspond to a three-layer rotation, while low values denote anti-solar rotation (compare Figures 5c,d). (b) The magnitude of the zonal flow, Δuϕeq = |⟨uϕ⟩max − ⟨uϕ⟩min|, in the equatorial plane (in green). The rms of the magnetic Lorentz force is shown in red. (c) The two energy states of the flow – high and low – in the phase space of the total magnetic energy and the energy of the axisymmetric (zonal) flow, with the points colored by the rms of the Lorentz force, for Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$.

To explain the physical mechanism behind these reversals, i.e., whether magnetic fields are driving re-configuration of the zonal flows, we clarify the causality between the growth and decay of the magnetic field and the amplitude of the zonal flows. To this end, we plot a proxy for the strength of zonal flows, Δu_ϕ^eq = |⟨u_ϕ⟩^max − ⟨u_ϕ⟩^min| at the equator, together with the rms of the Lorentz force F^Lorentz ∝ (∇×B)×B, as a function of time in Figure 6b. The strength of the shear increases by a factor of two during the transitions from three-layer to anti-solar rotation, well before the peaks in $E_{\mathrm{mag}}^{\mathrm{tot}}$ develop. This indicates that the development of anti-solar rotation that temporally enhances magnetic energy. On the other hand, it is when magnetic field is the strongest that the magnitude of differential rotation begins to decrease; these peaks coincide with the maxima in the rms of the Lorentz force, thus acting as magnetic damping on the zonal flow. Subsequently, the system returns to low-energy, three-layer rotation state. Figure 6c summarizes these dynamical transitions in the phase space of the total magnetic energy, the axisymmetric kinetic energy of the zonal flow, and the Lorentz force.

4.2. Correlation between the flow state and strength of dipolar modes

Finally, we link the flow cycles in Sect. 4.1 and the appearance of large-scale axisymmetric magnetic dipoles in Sect. 3.4. At first, we compare the radial distribution of magnetic energy contained in the first 6 spherical harmonics with azimuthal wavenumber m = 0, corresponding to the large-scale axisymmetric field structures, during high- and low-energy states for the run with N_ρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$.

Figure 7a shows the radial profile of these modal energies at a moment in time (Ωt ≈ 6900) during the high-energy state with strong anti-solar rotation. The energy in all large-scale harmonics are comparable, with relatively small-scale component of l = 6 peaking in the middle of the computational domain. On the other hand, the energy of axisymmetric dipolar mode with l = 1, m = 0 is much larger than the rest during a later moment with low magnetic energy and suppressed, three-layer differential rotation, even if the energy content decays by a factor of ten in all modes (Figure 7b). The energy of the dipolar mode peaks at small radii of r/r_o ≈ 0.4 (Figure 7b), which receive their energy contribution predominantly from the quiescent areas within the tangent cylinder (Figure 5c). This indicates that axisymmetric magnetic flux accumulates within the tangent cylinder, consistently with stronger surface fields at high latitudes from the butterfly diagrams (Figure 3c). The first minima of the dipolar mode correlates with the location of the prograde jet in Figure 5c. When the axisymmetric modes in Figure 7b are normalized by the radial distribution of the cumulative large-scale energy, the dipolar mode is still consistently larger then the rest across the domain.

Fig. 7. (a,b) Magnetic energy in the first six axisymmetric spherical harmonics (m = 0) as a function of spherical radius r at Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$. (a) High-energy state, anti-solar rotation, trot ≈ 6900. (b) Low-energy state, three-layer rotation, trot ≈ 7300. The relative buildup in dipolarity with respect to large scales takes place only during low-energy states with weakened three-layer rotation. (c) The energy in axisymmetric dipole normalized with the total large-scale magnetic energy up to 7th spherical harmonic, ∑El, m ≤ 7, at $\stackrel{\sim }{\mathit{Ra}}=32$. The modal energy was averaged over dipolar intervals for Nρ = 4, 6, and over the whole time evolution for Nρ = 2.

There is also a systematic increase in axisymmetric dipole energy with increase of density stratification. To see it, we average the axisymmetric dipolar mode with l = 1, m = 0 over several snapshots during dipolar periods and normalize it by the total large-scale magnetic energy up to the 7th spherical harmonic, ∑E_{l, m ≤ 7}. The averages are taken over the low-energy dipole intervals for N_ρ = 4 and 6 with three-layer differential rotation, and over the whole time integration of N_ρ = 2, since it does not exhibit intervals of predominant dipolarity. In Figure 7(c), we compare this relative radial distribution of axisymmetric dipoles for the three values of N_ρ. The strength of the dipolar component indeed consistently increases with density stratification across the whole domain, especially at low radii. The buildup of dipolarity at the surface is correlated with the increase of the dipolar energy in the bulk in the low-energy states, although it is partially conditioned by the magnetic boundary conditions – the energy of all magnetic modes decays at the surface (see Figure 7a,b). For this reason, we find that relative volume-distributed characteristics of modal energy, such as those in Figure 7a,b, are better indicators of low- and high-energy dynamo states in Figure 6 than, for example, local in time, weaker peaks in surface dipolarity in Figure 3.

5. Mechanisms of dipole field generation

5.1. Buildup of dipolarity with an increase of N_ρ and $\stackrel{\sim }{\mathit{Ra}}$

Even though only flow states such as those in Figure 3c have prolonged intervals of dipolar magnetic topology in our simulations, the average values of the dipolarity parameter $f_{\mathrm{dip}}^{11}$ systematically increases with $\stackrel{\sim }{\mathit{Ra}}$ and N_ρ. To illustrate this, we plot in Figure 8 the values of $f_{\mathrm{dip}}^{11}$, averaged in time, for the three different paths in the parameter space (Figure 2b). In the first one, with low stratification, N_ρ = 2 and variable convective strength $\stackrel{\sim }{\mathit{Ra}}$, the relative strength of dipolar mode remains approximately the same even when the convection is three times stronger (Figure 8a). In the second one, with high stratification N_ρ = 6 (same figure), dipolarity linearly increases with Rayleigh number with a scaling of $f_{\mathrm{dip}}^{11}\sim 0.01\stackrel{\sim }{\mathit{Ra}}$. The third path, where convection intensity is fixed, $\stackrel{\sim }{\mathit{Ra}}=32$, and stratification is varied, is shown in Figure 8b. It also exhibit a linear scaling of $f_{\mathrm{dip}}^{11}\sim 0.03N_{\rho }$, with a somewhat stronger coefficient than the scaling with $\stackrel{\sim }{\mathit{Ra}}$. We note that these are the mean values of dipolarity; its instantaneous maxima during dipolar intervals are higher (Figure 3). With $f_{\mathrm{dip}}^{11}$ at N_ρ = 6 and $\stackrel{\sim }{\mathit{Ra}}=32$ as a reference, these results suggest that the strength of the dipolar mode would systematically increase further with a simultaneous increase in both N_ρ and $\stackrel{\sim }{\mathit{Ra}}$ and thus for realistic stellar densities further away from the convective onset. In the rest of this section, we aim to explain the mechanisms of this gradual buildup in dipolarity.

Fig. 8. (a) Average surface dipolarity of the dynamo solutions, given by equation (4), as a function of Ra at Nρ = 6 (circles) and Nρ = 2 (triangles). (b) The same quantity as a function of Nρ at Ra = 32. The averages were performed over dipolarity time series such as those in Figure 3.

5.2. Surface compression of small-scale magnetic fields

For relatively low stratification, the area occupied by the positive upflows and negative downflows is approximately the same at the surface. The color map in Figure 9a illustrate this through the distribution of radial velocity component u_r, perpendicular to the surface plane, for N_ρ = 2 and $\stackrel{\sim }{\mathit{Ra}}=32$. As density gradient between the two spherical shell becomes more pronounced, convective cells also become more anisotropic. For N_ρ = 6, $\stackrel{\sim }{\mathit{Ra}}$, although the upward motions occupy larger surface, they are compensated by much stronger and localized downward sweeps (Figure 9b). Such asymmetry results from pressure fluctuations enhancing buoyancy driving in downflows and braking buoyancy in upflows (Hurlburt et al. 1984). Hence, the probability of having extreme negative velocity fluctuations at the surface for N_ρ = 6 is much larger than for N_ρ = 2, and so the whole distribution of u_r is skewed toward negative values (Figure 9c, in red). This does not take place at N_ρ = 2 where the distribution of u_r is much more symmetric (Figure 9c, in green), and probabilities of up- and down-flows are nearly equal.

Fig. 9. (a) Instantaneous snapshots of the radial magnetic field Br (contours) and radial velocity component vr (color) for Nρ = 2, Ra/Racr = 32 at r/ro = 0.9. The contour levels are [0.1Brmin, 0.05Brmin, 0.05Brmax, 0.1Brmax]. (b) The same for the stronger stratification, Ra/Racr = 32, Nρ = 6. (c) Two-dimensional probability density of the joint distribution of b2 and vr near the surface, r/ro = 0.9. Left, Nρ = 6; right, Nρ = 2. $\stackrel{\sim }{\mathit{Ra}}=32$. Colors are spaced logarithmically. (d) Unbiased sample skewness coefficient G1 (5) as a function of $\stackrel{\sim }{\mathit{Ra}}$. Red, Nρ = 6, blue, Nρ = 4, green, Nρ = 2.

This property of the probability distribution can be characterized by the Fisher-Pearson standardized moment coefficient,

$$\begin{array}{c}\hfill G_1=\frac{\sqrt{k(k-1)}}{k-2}\frac{(1/k){\displaystyle {\sum }_{i=1}^k}{(u_{r,i}(r)-\overline{u_r}(r))}^3}{{\left[(1/k){\displaystyle {\sum }_{i=1}^k}{(u_{r,i}(r)-\overline{u_r(r)})}^2\right]}^{(3/2)}},\end{array}$$(5)

which serves as a measure of the asymmetry of the velocity distribution at given r (here r/r_o = 0.9). We calculated this coefficient from a series of k flow snapshots for different values of $\stackrel{\sim }{\mathit{Ra}}$ and N_ρ, and plotted it in Figure 9d. For N_ρ = 2, the skewness coefficient is relatively small and nearly saturates at the largest values of $\stackrel{\sim }{\mathit{Ra}}$. With increase of stratification to N_ρ = 4 or 6, the left tail of the probability density distribution for u_r becomes longer and strong negative downflows progressively more plausible. Here, the distribution skewness does not show signs of saturation with $\stackrel{\sim }{\mathit{Ra}}$, indicating that surface convection will become even more anisotropic if enhanced. We note that the mean values of u_r become more positive as the skewness of its distribution increases, since the surface area of upward motions expands (Figure 9b).

This process is related mathematically to development of convergent flows in stratified flows. In anelastic approximation, mass conservation equation takes the form

$$\begin{array}{c}\hfill \mathrm{\nabla }·(\rho (r)\mathbf{u})=0,\rho (r)\mathrm{\nabla }·\mathbf{u}+u_r\frac{\partial \rho (r)}{\partial r}=0.\end{array}$$(6)

In compressible fluids, ∇ ⋅ u ≠ 0 and is proportional to u_r according to (6). Since density is a positive-definite function that decreases along spherical shell, ∂ρ/∂r < 0, convective downflows with u_r < 0 are linked with the regions of convergent flows, ∇ ⋅ u < 0. These flows are involved in the equation for magnetic energy:

$$\begin{array}{c}\hfill \frac{1}{2}\frac{\partial B^2}{\partial t}=-B^2(\mathrm{\nabla }·\mathbf{u})+\mathbf{B}(\mathbf{B}·\mathrm{\nabla })\mathbf{u}-\mathbf{B}(\mathbf{u}·\mathrm{\nabla }\mathbf{B})+\cdots ,\end{array}$$(7)

where the first term reflects compression of magnetic field lines and amplification of magnetic energy due to convergent flows. We illustrate this process with 2D joint probability density distributions for u_r and the radial magnetic field B_r (Figure 9c). The probability densities were computed by binning the data for u_r and B_r into 100 and 50 bins of equal size, respectively, in the same interval of u_r ∈ [ − 2000, 2000] and B_r² ∈ [0, 2000].

The turbulent layer, forming at the surface for N_ρ = 6, is associated with strong magnetic induction process near the surface, generating small-scale magnetic structures. Similarly to Ro_l, magnetic Reynolds number Rm = u_rmsd/η (see Table B.1) also becomes progressively higher with increasing N_ρ and $\stackrel{\sim }{\mathit{Ra}}$, enhancing induction. As a result, the characteristic value of magnetic energy fluctuations near the surface also increases with N_ρ, up to 2 times in Figure 9c. For N_ρ = 6, strong downward motion from the negative, long tail of the velocity distribution are associated with the enhanced probability of stronger magnetic energy, unlike the upflows of the “shorter", positive tail (Figure 9c, in red). This effect is much less present for N_ρ = 2 (in green), as the energy of magnetic fluctuations is distributed together with the more symmetric distribution of velocity fluctuations. In other words, negative fluctuations of velocity field appear more correlated with the more energetic field lines in stratified turbulence. We attribute this effect to strong anisotropic sweeps being able to compress magnetic field fluctuations induced by small-scale convection at the surface through the term proportional to B²∇⋅u in (7). This process is illustrated by the contours of B_r in Figure 9, which are much more aligned with the localized downflows for N_ρ = 6 (panel b) than for N_ρ = 2, where they are more or less evenly spread across the domain surface (panel a). This compression provides a mechanism of inward transport of magnetic energy generated by strongly turbulent, rotationally unconstrained small-scale motions in outer layer with Ro_l > 1 (Figure 2a), suggesting importance of this region for the dynamo. To see how it is results in the transfer of magnetic energy to the large scales, we relate it to the global, mean-field advection of magnetic field by the flow – magnetic pumping.

5.3. Magnetic pumping

Magnetic pumping is a mechanism of large-scale magnetic field generation in the mean-field dynamo theory (Moffatt 1978). Using a set of approximations, the electromotive force, ε, generating the mean-field from correlated fluctuations of velocity and magnetic field can be parametrized as

$$\begin{array}{c}\hfill \epsilon =⟨\mathit{u}\times \mathit{B}⟩\approx \alpha ⟨\mathbf{B}⟩+\mathit{\gamma }\times ⟨\mathit{B}⟩+\cdots ,\end{array}$$(8)

where ⟨B⟩ is large-scale, axisymmetric magnetic field. In the expression above, higher-order contributions, proportional to derivatives of ⟨B⟩ and corresponding to the turbulent diffusion, were omitted for simplicity. In (8), α-effect generates large-scale field through perturbing magnetic field lines by cyclonic motions, and γ × ⟨B⟩ describes the pumping of magnetic field with a speed γ. The mean-field induction equation shows that γ is an addition to the mean velocity that also acts on the mean magnetic field:

$$\begin{array}{c}\hfill {⟨\mathit{B}⟩}_t=\mathrm{\nabla }\times (⟨\mathit{U}⟩\times ⟨\mathit{B}⟩+\mathit{\gamma }\times ⟨\mathit{B}⟩+\alpha ⟨\mathbf{B}⟩-{\eta }_T\mathrm{\nabla }\times ⟨\mathit{B}⟩).\end{array}$$(9)

In essence, magnetic pumping is a turbulence-generated component of the electromotive force (Drobyshevski 1977) that can help redistribute magnetic energy to large scales. In rotating convection, it can be viewed as a systematic transport of magnetic energy by anisotropic convective columns toward the bulk of the domain: convective structures are stretched vertically if their location is perturbed inward, and this increases the length scales of magnetic flux tubes that are carried by them (Schrinner et al. 2012). Equation (9) shows that this effect is capable of generating large-scale magnetic fields. However, evaluating γ is not straightforward. The most rigorous way is the test-field method, but it requires solving a set of additional equations (Schrinner et al. 2012) for auxiliary variables. Another frequently used method is the method of snapshots, i.e., point-wise reconstruction of the tensor coefficients from time series of ε and mean magnetic fields at each point in the computational domain (Simard et al. 2016). This method is more straightforward to use, yet it suffers from noise in strongly turbulent conditions such as those explored here and may not always give correct results (Warnecke et al. 2018). To unambiguously quantify the transport of magnetic energy across the domain for different parameter regimes in our simulations, we used an indirect measure of this process, the Poynting flux, averaged over longitude and in time:

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{Poynting}}\propto {⟨(\mathbf{v}\times \mathbf{B})\times \mathbf{B}⟩}_{\varphi ,t}={⟨B^2{\mathbf{u}}_{\perp }⟩}_{\varphi ,t}.\end{array}$$(10)

This vector describes the direction in which magnetic energy is transported by the velocity u_⊥ perpendicular to magnetic field lines.

In Figure 10 we plot the distribution of the radial component of Poynting flux with radius, integrated in latitude θ, for three series of simulations in Figure 2b. The influence of stratification on the flux at a given $\stackrel{\sim }{\mathit{Ra}}=32$ is relatively clear: while for N_ρ = 2 the flux is nearly zero throughout the domain, it becomes negative for N_ρ = 4, reflecting the inward transport of magnetic energy toward the inner sphere (Figure 10a). For N_ρ = 6, this flux becomes more localized, with a strong peak developing near the surface, at r/r_o = 0.9. This flux results from consistent enhancement of turbulence near the outer shell, and becomes more and more pronounced as $\stackrel{\sim }{\mathit{Ra}}$ grows (Figure 10b). This indicates that the energy generated and compressed by anisotropic convection in the turbulent surface layer is systematically transported into the deep interiors. On the other hand, in weakly stratified flows (N_ρ = 2), the enhancement of convection results in magnetic flux becoming predominantly positive throughout the domain, and the outward transport of magnetic energy (Figure 10c). This is potentially related to inversion of magnetic pumping effect: energetic turbulent motions, developing near the inner sphere, expel outward generated by them magnetic fluctuations. Although a region of negative flux develops around r/r_o ≈ 0.9 at the outer surface for $\stackrel{\sim }{\mathit{Ra}}=96$, it is entirely compensated by even stronger flux expulsion region at r/r_o ∈ [0.35, 0.8]. Correspondingly, the average energy of large-scale dipole fields remains low with increase of $\stackrel{\sim }{\mathit{Ra}}$ in weakly stratified flow, while for strongly stratified cases dipolarity increases (Figure 8a).

Fig. 10. (a) Poynting flux (10), integrated in latitude θ, as a function of spherical radius r and Nρ, $\stackrel{\sim }{\mathit{Ra}}=32$. (b) Poynting flux as a function of $\stackrel{\sim }{\mathit{Ra}}$, Nρ = 6. (c) Same but for Nρ = 2.

We note that the dependence of Poynting flux on the radius is non-monotonic. For example, the flow with N_ρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$ case has two regions of negative flux buildup, separated by a minimum at r/r_o ≈ 0.7. It also varies with parameter regime: unlike for $\stackrel{\sim }{\mathit{Ra}}=32$, $\stackrel{\sim }{\mathit{Ra}}\le 16$, the flux is enhanced nearly throughout the whole domain. This suggests that the flux distribution is a 2D function of both r and θ. To clarify this point, we show the streamlines of the poloidal component of the Poynting flux, ${\mathbf{F}}_{\mathrm{Poynting}}^P=F_r{\mathbf{e}}_{\mathbf{r}}+F_{\theta }{\mathbf{e}}_{\theta }$ in Figure 5. The corresponding magnetic regimes are described in Sect. 3. For N_ρ = 2, where the integrated in θ flux is weak, the flow is separated in three different regions – inward-flux zone near the equator, outward around the mid-gap, and another inward-flux zone near the inner sphere (Figure 5a). For N_ρ = 6 and $\stackrel{\sim }{\mathit{Ra}}=16$, where the large-scale coherent waves develop, Poynting flux is consistently directed inward across the entire domain, except for a narrow zone at the equator. Unstable dipole dynamos (N_ρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$) develop three interchanging positive and negative regions of magnetic flux (Figure 5c), similarly to weakly stratified case in panel (a). But unlike N_ρ = 2, in this case the width of the inner region of inward-directed flux is considerably larger, and occupies more than half of the domain; there is a clear negative (positive) vertical flux of magnetic energy in upper (lower) hemisphere. This suggests more efficient transport of magnetic flux in the polar regions, similarly to flux distribution of axisymmetric large-scale dynamo waves (Figure 5b). We note that dipolar field component is most prominent in the polar regions inside the tangent cylinder (Figure 3c and also meridional slices of ⟨B_r⟩_ϕ, not shown).

5.4. Poynting flux and zonal flows

Figure 5 also demonstrates that the local direction of the Poynting flux and the effective surface from which magnetic flux is advected inside the domain is linked to the distribution of zonal flows. The distribution of magnetic flux lines for weakly stratified turbulence at N_ρ = 2 and $\stackrel{\sim }{\mathit{Ra}}=32$ shows that the convergent flux zone at the equator is the bulk is correlated with the area of strong shear between the prograde and the retrograde jets, where part of the Poynting flux is transformed from poloidal to toroidal (Figure 5a). The maximum of the retrograde jet in the bulk coincides with divergent flux zone, and seemingly isolates the area of flux accumulation inside the tangent cylinder. Since the jet is located close to the tangent cylinder, the surface from which the flux accumulates is effectively small. The flows with strong stratification and weaker turbulence (e.g., N_ρ = 6, $\stackrel{\sim }{\mathit{Ra}}=16$), featuring large-scale oscillatory dynamos, also exhibit anti-solar rotation (Figure 5b). However, inner retrograde flow in the bulk is very weak in this case, and so the flux is directed predominantly radially everywhere and collected from nearly all latitudes.

The flows with N_ρ = 4, 6, and $\stackrel{\sim }{\mathit{Ra}}=32$, develop two different types of differential rotation: weaker three-layer and stronger anti-solar (Figure 5c,d). Although a fraction of Poynting flux is sheared near the surface by the interface between the equatorial retrograde jet and its neighboring weak prograde flow during the intervals with three-layer rotation, the prograde jet is located further away from the inner sphere (Figure 5c). The corresponding quiescent area in the bulk, where radially directed flux accumulates, occupies the area twice as large than that of N_ρ = 2 (Figure 5a). In the second configuration of anti-solar rotation, Poynting flux is again directed downward in the equatorial region of the surface retrograde jet. However, this flux is sheared at the interface between the two layers of the zonal flow and becomes less coherent within the prograde jet area. This results in strengthening of axisymmetric toroidal component of magnetic field in the equatorial region during multipolar, high-energy intervals (not shown here). In the meanwhile, vertical advection of magnetic energy inside the tangent cylinder is disrupted, and the large-scale dipole component is lost. For weakly stratified and strongly turbulent cases such as N_ρ = 2, $\stackrel{\sim }{\mathit{Ra}}=96$, similar configuration of zonal flows, accompanied by a stronger prograde jet near the tangent cylinder (Figure 2b, inset), results in positive sign of Poynting flux and expulsion of magnetic energy from the bulk (Figure 10c).

6. Discussion

Recent observations have shown magnetic dynamics either in the form of periodic oscillations or reversals are encountered in low-mass stars (Jeffers et al. 2023). While many of them have a dipolar magnetic topology, others show a mix of dipolar and quadrupolar magnetic fields (Willamo et al. 2022). The magnetic topology and dynamics in our simulations, as given by Figures 3 and 4, are representative of these observations. In our most strongly stratified and turbulent runs, such as those in Figure 3c, which are closer to stellar conditions, the magnetic field behaves aperiodically. The reversals and switches between the dipolar and multipolar topologies through dipole decay and buildup take place on a timescale of 50 − 100 rotation units. With the average solar rotation rate of 28 days as a reference, the reversals of the magnetic field in our model would take place on a timescale of 4-8 years, which is compatible with spectropolarimetric observations. We note that to make a one-to-one comparison with ZDI magnetic maps, one needs to take into account the inclination of the rotation axis and the number of observations per rotation period of a star (Hackman et al. 2024).

Our analysis of correlation between radial velocity and magnetic field at the surface, together with enhancement of inward-directed Poynting flux, indicates that magnetic pumping mechanism is responsible for development of such large-scale axisymmetric magnetic states in strongly stratified turbulence. Development of turbulent, rotationally unconstrained thin surface layer at N_ρ = 4 and N_ρ = 6 (Figure 2b), accounting for 10 − 20% of the entire domain, is a natural outcome of such convection. Such a layer can contribute to both induction of magnetic fluctuations, and to transport of this energy inward (Figure 9c). We find that the increase in dipolarity with N_ρ and $\stackrel{\sim }{\mathit{Ra}}$ in our simulations (Figure 8) is correlated with enhancement of inwardly directed magnetic Poynting flux (Figure 10a,b), as opposed to low-stratified turbulent cases where magnetic energy is expelled outward (Figure 10c). These results are in agreement with previous studies of magnetic pumping effect in Cartesian geometry, for example with Tobias et al. (1998), who found that strong downflows in stratified convection can advect magnetic field sheets initially located at the top of convection zone downward, thereby concentrating magnetic energy at radiative-convective interface. This effect is also considered one of the main ingredients in solar dynamo theory, as it is able to counteract magnetic buoyancy bringing magnetic flux tubes to the surface (Brandenburg 2005). More recent studies in global spherical geometry also found compression to be an important mechanism for generation of strong magnetic fields, which enforce solar-like prograde rotation at the equator through the magnetic transport of momentum (Hotta et al. 2022).

The magnetic pumping coefficient γ is a velocity-like vector, and therefore magnetic pumping is in fact a three-dimensional effect (Ossendrijver et al. 2002). To confirm its presence in our simulations, in Figure 11a,b we show the radial and the latitudinal components of γ, as derived in equation (8), calculated using the method of snapshots (Simard et al. 2016) on the data from dipolar solution at N_ρ = 4, $\stackrel{\sim }{\mathit{Ra}}=32$. The radial component of γ-effect is negative in the areas inside the tangent cylinder at high latitudes, indicating overall inward advection of magnetic field in these areas (Figure 11a). On the other hand, in the outer, equatorial region of strong shear, γ_r is positive, indicating outward pumping. The latitudinal component, γ_θ, is anti-symmetric with respect to the equator, with reflects equivalent converging flows in the equatorial plane and divergent flows in deep interiors (Figure 11b). These results are consistent with manifestation of the large-scale field predominantly at high latitudes (Figure 3). Moreover, the spatial distribution of γ_r and γ_θ is compatible with previous results for stratified convection leading to dipolar solutions (Schrinner et al. 2012), which suggests that axisymmetric magnetic pumping mechanism, induced by the combination of stratification and strong turbulence at the surface, is indeed responsible for the encountered in simulations dipoles. Such mechanism could be also responsible for the dipolar stellar-like dynamo with co-existing small- and large-scale fields at the surface reported in Yadav et al. (2015), although with much smaller values of Ro_l ≤ 0.04 and thus more rotationally constrained. As N_ρ is increased, γ obtained with the method of snapshots becomes considerably more noisy yet its structure is similar; we leave its detailed analysis for the future work.

Fig. 11. (a,b) Components of the mean-field γ-coefficient calculated for Nρ = 4, $\stackrel{\sim }{\mathit{Ra}}=32$ using the method of snapshots (Simard et al. 2016). (a) Radial component, γr. (b) Latitudinal component, γθ. (c) Total magnetic energy near the surface (r/ro = 0.9) for run Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$ as a function of spherical harmonic degree l for high- (tH) and low-energy (tL) dynamo states.

Another measure of spatially distributed transport of electromagntic energy, meridional distributions of Poynting flux, is vertically aligned due to the system rotation and depends on the distribution of zonal flows (Figure 5), suggesting that differential rotation plays an important role in re-distributing the flux across the domain. Stronger zonal flows in the area of tangent cylinder prevent magnetic flux from accumulating there and result in non-axisymmetric solutions with smaller scales (Figure 5a,d). Weakened zonal flows in this area are accompanied by more axisymmetric magnetic fields, either quadrupolar or dipolar (Figure 5b,c). In agreement with Figure 11, more quiescent polar zones seems to be the locations where magnetic energy transport is most prominent (Figure 5b,c). In all these cases, the distribution of Poynting flux is strongly constrained by zonal flows, with the regions of strong differential rotation acting as radial flux attractors and transforming it into poloidal flux and fields.

The transition between periodic large-scale waves and dipolar solutions has been attributed previously to increase of the period of oscillatory dynamo solutions due to fluctuations in α-effect in simplified one-dimensional mean-field dynamo models (Richardson & Proctor 2010). We suppose that similar mechanism can be responsible for destabilization of the periodic large-scale waves at N_ρ = 4 and 6 with increase of $\stackrel{\sim }{\mathit{Ra}}$, although a rigorous analysis of temporal fluctuation of α-effect and the corresponding mean-field equation is needed to confirm this in the future. In our simulations, such non-oscillatory dipolar solutions are unstable and exhibit either reversals of polarity or stochastic transitions to a multipolar state with higher energy, as described in Sect. 4. We interpret this unstable behavior as a predator-prey interaction of zonal flow and magnetic field: the system attempts to escape to a state with strong anti-solar rotation and strong turbulent magnetic fields, while increasing magnetic tension from these fields pushes it back to the original state with weaker, dipolar magnetic field and weaker three-layer rotation. This is particularly visible from the phase space map in Figure 6c, with flow trajectories colored by the Lorentz force. Such predator-prey transitions have been observed before both in hydrodynamic convection (Busse & Simitev 2006) and in dynamo simulations as Type 2 modulation producing minima in magnetic activity (Raynaud & Tobias 2016). We suggest that such quasi-stable flow states result from the flow being at the verge of transition between anti-solar and solar rotation. While magnetic fields tend to postpone such a transition compared to pure hydrodynamic flows (Fan & Fang 2014; Karak et al. 2015), such a transition nevertheless takes place in the range of about $\stackrel{\sim }{\mathit{Ra}}\in [70-80]$ in our magnetohydrodynamic simulations for the least numerically demanding regime of N_ρ = 2 (Figure 2b). Development of the strongly turbulent shear layer with Ro_l ≤ 1 at the surface is correlated with the transition between solar- and anti-solar rotation for all explored N_ρ, which suggests that the presence of this layer can induce this transition. We thus expect that with further enhancement of convective forcing the three-layer differential rotation (Figure 5c) would be lost and anti-solar state would become the only solution. Yet the aperiodic switches between the two different shear topologies in our aperiodic dipole solutions suggest our parameter regime can be close to bi-stability of magnetically controlled transition between solar- and anti-solar rotation, essentially for a stronger density contrast than ρ_i/ρ_o = 12 reported by Karak et al. (2015), who did not find bi-stable behavior. The three-layer structure of differential rotation, obtained thanks to the high N_ρ, is an interesting new result that resembles the near-surface shear layer in the Sun (Vasil et al. 2024); in the future, we will extend our study to even higher N_ρ to study the stability of this structure with $\stackrel{\sim }{\mathit{Ra}}$ and magnetic field. The influence of the latter can be particularly important, since magnetic field impacts both alignment or non-alignment of the flow with the axis of rotation (Pinçon et al. 2024), and the sign of differential rotation (Brun et al. 2022).

Finally, we show the energy spectra during low- and high-energy states in the outer region (r/r_o = 0.9) for N_ρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$, as a function of spherical harmonic degree l (Figure 11c). The spectra are in agreement with the time series of magnetic energy for the two flow regimes in Figure 6, with energy decreased across all spherical harmonics during quasi-dipolar state. The l-spectra in this regime shows a certain degree of separation between the large and the small scales, with a peak on small scales around l = 20 and buildup of energy at the dipolar mode with l = 1. It is this buildup that is observed as enhanced dipolarity, when magnetic field is filtered at the surface. During the high-energy state, all spherical harmonics but the smallest ones receive additional amount of energy, so that the spectrum fills up and does not show anymore the signs of scale separation. It is possible that this scale separation appears due to spatial anisotropy of stratified convection, separating generation of small-scale magnetic fluctuations at the surface from large-scale fields generated through magnetic pumping. However, even stronger density stratification is necessary to see whether the trend of scale separation becomes more pronounced; we leave exploring this idea in depth for the future work.

7. Conclusions

Previous studies of stratified convective dynamos by Gastine et al. (2012) and Raynaud et al. (2015) reported dipole collapse with an increase in density stratification and explored parameter regimes with lower values of $\stackrel{\sim }{\mathit{Ra}}$ and N_ρ ≤ 3 than those explored in this work. Essentially, these dipolar dynamos correspond to classic steady low-Ro dynamos and are quickly destabilized by inertia when Ro ∼ 0.1 (Christensen & Aubert 2006). In this work, we explored more stratified and turbulent models, up to N_ρ = 4, 6 and $\stackrel{\sim }{\mathit{Ra}}=32$, which resulted in larger integral values of Ro > 0.1 (Table B.1).

In all of our models, magnetic fields feature a broad spread of energy across spherical harmonics with many small-scale magnetic fluctuations at the surface (Figure 11c). According to the classic definition of dipolarity as the ratio between the energy of the dipole and the total magnetic energy (used by Gastine et al. (2012)), such fields would be classified as multipolar and small scale. Here, we compared the strength of the dipolar mode with the energy of up to 7 − 11 large-scale spherical harmonics (equation 4). The strength of the dipolar component consistently increases with both stratification and turbulence vigor, suggesting similar magnetic regimes are possible in low-mass stars, which frequently feature dipolar magnetic fields. On the large scales, our models reproduce magnetic topologies commonly observed in these stars: multipolar, oscillatory, and dipolar. They also reproduce the main dynamical ingredients of such fields, such as dipolar reversals, and transitions between dipolar and multipolar regimes. With the averaged solar rotation rate as a reference, reversals of the magnetic field in our model would take place on the timescale of 4 to 8 years, which is compatible with spectropolarimetric observations of stellar cycles (Jeffers et al. 2023). We argue that this large-scale definition of dipolarity (4), defined up to a suitable low-order truncation threshold l_cut, better characterizes magnetic topology in stellar dynamo models given that large-scale fields detected by spectropolarimetry are an order of magnitude smaller than total fields from stellar photometric observations.

We identified that one of the direct consequences of increasing both stratification and turbulence is the development of the turbulent surface layer with surface Ro_l > 1 and weakened influence of the Coriolis force, while the deep interiors of the flow become increasingly rotationally constrained (Figure 2a). Strong downflows promoted by the anisotropic convection of the turbulent surface layer, together with stratified and more quiescent convective columns in the flow interiors, induce an axisymmetric magnetic pumping mechanism that scales with N_ρ and $\stackrel{\sim }{\mathit{Ra}}$. This physical mechanism could be relevant for dominant magnetic mode selection in other irregular magnetic cycles previously reported in the slowly rotating stellar models of Yadav et al. (2016), Strugarek et al. (2018), Viviani et al. (2018) at lower stratification. Our results suggest that magnetic pumping consistently enhances the dipolar component in our simulations to the point where prolonged dipolar periods on the order of 500 − 1000 rotation times are observed for N_ρ = 4 and 6, $\stackrel{\sim }{\mathit{Ra}}=32$. These runs are on the bounds between solar- and anti-solar transition, and they exhibit dynamical interaction of the magnetic field and the zonal flows, transitioning between the low-energy dipolar state with weakened three-layer rotation and the high-energy multipolar states with a stronger anti-solar rotation. As in Brun et al. (2022), we observed that the Lorentz force can change the configuration of the differential rotation in simulations where Ro_l > 1; this phenomenon is favored when N_ρ is high. These dynamical interactions suggest additional observational constraints based on the relation between stellar fields and internal differential rotation. Such constraints could be obtained from joint spectropolarimetry and asteroseismology studies.

It is yet to be shown how magnetic pumping effects extrapolate to higher values of Ra and N_ρ and whether dipolar solutions are constrained by the transition between solar-antisolar rotation. However, the enhancement of large scales, together with the emergence of the surface dipole and quadrupole fields in strongly turbulent convection suggests a new framework for modeling dynamical magnetic effects in the DNS models of stellar convection. This framework supposes seeking models and parameter regimes that (i) capture dynamics of large magnetic field scales and (ii) capture the energy balance between large- and small scales that depends on the depth. Such models will be the aim of our future work. Although magnetic dipoles are not steady states in our simulations, their existence at the interface of oscillatory quadrupolar and non-coherent multipolar solutions indicates a parameter path for further exploration of such states in stellar dynamos, which is extremely expensive numerically in these parameter regimes.

Acknowledgments

The authors acknowledge financial support from the French program ‘PROMETHEE’ (Protostellar Magnetism: Heritage vs Evolution) managed by Agence Nationale de la Recherche (ANR), and would like to thank E. Alecian and the PROMETHEE team for inspiring discussions that helped to shape this work. The authors are also grateful to F. Daniel and S. Tobias for their insightful comments, and the anonymous reviewer for their constructive feedback. This work was supported by the “Action Thématique de Physique Stellaire” (ATPS) of CNRS/INSU PN Astro co-funded by CEA and CNES, and the HPC facility MesoPSL for the computational resources.

References

Appendix A: Anelastic model

In this work, we use LBR anelastic approximation, as described by Braginsky & Roberts (1995), Lantz & Fan (1999), Jones et al. (2011). The reference thermodynamic state of the system is assumed to be close to adiabatic hydrostatic equilibrium, and follows the polytropic structure for pressure, density and temperature,

$$\begin{array}{c}\hfill P=P_{\mathrm{c}}{w}^{n+1},\rho ={\rho }_{\mathrm{c}}{w}^n,T=T_{\mathrm{c}}w,\end{array}$$(A.1)

with

$$\begin{array}{cc}& w=c_0+\frac{c_{1}d}{r},c_0=\frac{2w_{\mathrm{o}}-\chi -1}{1-\chi },\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& c_1=\frac{(1+\chi )(1-w_{\mathrm{o}})}{{(1-\chi )}^2},w_{\mathrm{o}}=\frac{\chi +1}{\chi exp(N_{\rho }/n)+1},\hfill \end{array}$$(A.3)

with reference values T_c, ρ_c and P_c at the mid-gap between the sphere boundaries, and ξ = r_i/r_o.

In this approximation, we solve three equations for velocity field u, magnetic field B and the entropy of the flow S:

$$\begin{array}{cc}\hfill \frac{\partial \mathit{u}}{\partial t}+\mathit{u}·\mathbf{\nabla }\mathit{u}=& -\frac{\mathit{Pm}}{E}\mathbf{\nabla }\left(\frac{P^{\prime }}{w^n}\right)+\frac{RaP{m}^2}{\mathit{Pr}}\frac{S}{r^2}{\mathit{e}}_r-\frac{2Pm}{E}{\mathit{e}}_z\times \mathit{u}\hfill \\ \hfill & +\frac{\mathit{Pm}}{E{w}^n}(\mathbf{\nabla }\times \mathit{B})\times \mathit{B}+Pm{\mathit{F}}_{\nu },\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill \frac{\partial S}{\partial t}+\mathit{u}·\mathbf{\nabla }S=& \frac{\mathit{Pm}}{w^{n+1}Pr}\mathbf{\nabla }·\left(w^{n+1}\mathbf{\nabla }S\right)\hfill \\ \hfill & +\frac{\mathit{Di}}{w}[\frac{{(\mathbf{\nabla }\times \mathit{B})}^2}{E{w}^n}+Q_{\nu }]\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}=& \mathbf{\nabla }\times (\mathit{u}\times \mathit{B})+{\mathbf{\nabla }}^2\mathit{B}\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill \mathbf{\nabla }·\mathit{B}=& 0\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill \mathbf{\nabla }·\left(w^n\mathit{u}\right)=& 0.\hfill \end{array}$$(A.8)

Here P′ denotes the pressure perturbation of the reference state. The viscous force F_ν can be expressed F_ν = w⁻ⁿ∇S, where S is the rate of strain tensor,

$$\begin{array}{c}\hfill {\mathsf{S}}_{\mathit{ij}}=2w^n(e_{\mathit{ij}}-\frac{1}{3}{\delta }_{\mathit{ij}}\mathbf{\nabla }·\mathbf{u}),e_{\mathit{ij}}=\frac{1}{2}(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}).\end{array}$$(A.9)

The dissipation parameter Di and the viscous heating Q_νare given by

$$\begin{array}{c}\hfill Di=\frac{c_{1}Pr}{\mathit{PmRa}},Q_{\nu }=2[e_{\mathit{ij}}{e}_{\mathit{ij}}-\frac{1}{3}{(\mathbf{\nabla }·\mathit{u})}^2].\end{array}$$(A.10)

Equations above were scaled with the shell width d = r_o − r_i, magnetic diffusion timescale d²/η, by entropy difference ΔS = S(r_i)−S(r₀), pressure scale Ωρ_cη, density scale ρ_c, temperature by T_c, and magnetic field scale $\sqrt{\mathrm{\Omega }{\rho }_{\mathrm{c}}\mu \eta }$, where μ is the magnetic permeability. Resulting dimensionless parameters Pm, Pr, E, and Ra are described in the text in section 2.

Appendix B: Simulation parameters

In stratified convection, both the critical Rayleigh number Ra_cr, and the azimuthal wavenumber of the dominant convective mode, increase with density contrast (see, e.g., Raynaud et al. 2018). The higher is density stratification, the higher is the numerical resolutions needed to resolve all the flow scales. In table B.1 we give the values of Ra_cr employed in our simulations that allowed us, together with section 2 and figure 2b, to reconstruct the whole parameter space spanned by simulations in this work as well as several diagnostic quantities of each run. The term Ra_cr denotes the critical value of Ra for the onset of convection, $\stackrel{\sim }{\mathit{Ra}}=Ra/R{a}_{\mathit{cr}}$ denotes the distance from the onset, and Ra_F^* = RaNuEk³/Pr² denotes modified flux-based Rayleigh number. The term Nu is defined as the ratio between the luminosity of the convective model, where the integral is taken over the top spherical boundary, and the basic state luminosity of the conductive state:

$$\begin{array}{c}\hfill Nu=\frac{{\displaystyle {\int }_S}-\kappa \rho T\frac{\partial S}{\partial r}{r}^{2}sin\theta d\theta d\varphi }{4\pi n{c}_1{\xi }_i^n{(exp{N}_{\rho }-1)}^{-1}}.\end{array}$$(B.1)

In addition, we give the Rossby number Ro = u_rms^I/(l_conv^IΩ), based on the integrated over r convective length scale and rms velocity in equation (3), as introduced by Christensen & Aubert (2006). The Reynolds number and magnetic Reynolds numbers are defined as Re = u_rms^Id/ν and Rm = PmRe, respectively, and represent kinetic energy characteristics of each run. Finally, we present the energetics for magnetic field of the runs in the last two columns of table B.1, measured as the ratio of large-scale axisymmetric components $E_{\mathit{mag}}^{l\le 7,m=0}$ and the energy contained in all large scale harmonics $E_{\mathit{mag}}^{l\le 7}$, similarly to section 4.2. In the last column of table 2 we compare the ratio between the axisymmetric and the total magnetic energy. These are typically order of magnitude smaller compared to the large-scale estimates, except for the runs n41, n62, n63 which feature energetic quadrupolar dynamo waves described in section 3.3.

The difference between weakly stratified (section 3.2) and strongly stratified oscillatory solution (section 3.3) is that the latter contains much more energy in the axisymmetric modes, and shows approximately the same proportion of axisymmetry with respect to both large-scale and total magnetic energy (runs 41, n62, n63). On the other hand, the weakly stratified oscillatory solutions show (runs n21-24) have much weaker axisymmetric component, especially in comparison to the total magnetic energy. We note that axisymmetric energy decreases in proportion to the total magnetic energy for more stratified runs (N_ρ = 4 and 6). This is due to the development of surface layer generating small-scale, non-axisymmetric magnetic fluctuations, which are then factored in total energy estimates. The time-averaged characteristics in table B.1 also do not take into account change in energy as the dynamo undergoes dipolar-multipolar transitions.

All Tables

All Figures

	Fig. 1. Instantaneous surfaces of the vertical vorticity component of the flow at r/r0 = 0.5 (top) and r/r0 = 0.9 (bottom) illustrating scale separation developing with Nρ in strongly turbulent convection, $\stackrel{\sim }{\mathit{Ra}}=32$. (a,d) Nρ = 2; (b,e) Nρ = 4; (c,f) Nρ = 6.
In the text

Fig. 2. (a) Effect of varying Nρ (in red, Ra/Racr = 32) and Ra (in blue, Nρ = 6) on the local Rossby number calculated based on the convective length scale lconv (3). Dotted line denotes Rol = 1, e.g., a balance between inertia and the Coriolis force. (b) Density contrasts Nρ and Rayleigh numbers $\stackrel{\sim }{\mathit{Ra}}$ spanned by simulations in this work. Different symbols denote three dynamical regimes of the flow, discussed in more details in Sect. 3. Filled (empty) triangles, predominantly multipolar fields with solar (anti-solar) differential rotation (Sect. 3.2). Circles, quadrupolar oscillatory dynamos with solar-like differential rotation (Sect. 3.3). Diamonds, reversing dipole-multipole solutions with transitions between three-layered and anti-solar rotation (Sect. 3.4). Solar- and anti-solar differential rotation is shown in the inset plots depicting zonal flows ⟨uϕ⟩ϕ for Nρ = 2, $\stackrel{\sim }{\mathit{Ra}}=32$ (left) and $\stackrel{\sim }{\mathit{Ra}}=96$ (right).

In the text

Fig. 3. Butterfly diagrams of the axisymmetric, azimuthally averaged magnetic component ⟨Br⟩ϕ at the outer boundary. (a) Nρ = 2, Ra/Racr = 32; (b) Nρ = 6, Ra/Racr = 16; and (c) Nρ = 6, Ra/Racr = 32. To highlight large-scale components, the color maps were saturated with cutoffs of 0.2, 0.25, and 0.1, respectively; the magnitude of magnetic fluctuations is larger by a corresponding factor. Black curves correspond to the surface dipolarity with cutoff parameter of lcut = 11 (4).

In the text

Fig. 4. Instantaneous snapshots of the surface radial magnetic field Br at r/ro = 0.95 (top) with the same fields filtered to l ≤ 5 (bottom) for different parameters: (a,d) Nρ = 6, Ra/Racr = 32, multipolar fields; (b,e) Nρ = 6, Ra/Racr = 16, quadupolar waves; and (c,f) Nρ = 6, Ra/Racr = 32, during dipole period. The color maps in panels (a-c) were saturated to 0.5|Br|max to highlight the spatial structure of magnetic features.

In the text

Fig. 5. Zonal flows (in color) for different types of dynamo solutions corresponding to Figure 3. (a) Nρ = 2, $\stackrel{\sim }{\mathit{Ra}}=32$; (b) Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=16$; and (c) Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$, during the low magnetic energy interval with predominance of the dipole on large scales; and (d) Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$, during the high-energy interval without the dipole. Flux lines corresponds to the poloidal component of Poynting flux (see Sect. 5.4). Coordinates are denoted in the cylindrical coordinate system (rc, zc). Fluxes and zonal flows were averaged in time over multiple snapshots.

In the text

Fig. 6. (a) Total magnetic energy as a function of time (in black) versus the radius of ⟨uϕ⟩ϕmax in the equatorial plane (θ = π/2), in blue. The large radii correspond to a three-layer rotation, while low values denote anti-solar rotation (compare Figures 5c,d). (b) The magnitude of the zonal flow, Δuϕeq = |⟨uϕ⟩max − ⟨uϕ⟩min|, in the equatorial plane (in green). The rms of the magnetic Lorentz force is shown in red. (c) The two energy states of the flow – high and low – in the phase space of the total magnetic energy and the energy of the axisymmetric (zonal) flow, with the points colored by the rms of the Lorentz force, for Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$.

In the text

Fig. 7. (a,b) Magnetic energy in the first six axisymmetric spherical harmonics (m = 0) as a function of spherical radius r at Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$. (a) High-energy state, anti-solar rotation, trot ≈ 6900. (b) Low-energy state, three-layer rotation, trot ≈ 7300. The relative buildup in dipolarity with respect to large scales takes place only during low-energy states with weakened three-layer rotation. (c) The energy in axisymmetric dipole normalized with the total large-scale magnetic energy up to 7th spherical harmonic, ∑El, m ≤ 7, at $\stackrel{\sim }{\mathit{Ra}}=32$. The modal energy was averaged over dipolar intervals for Nρ = 4, 6, and over the whole time evolution for Nρ = 2.

In the text

	Fig. 8. (a) Average surface dipolarity of the dynamo solutions, given by equation (4), as a function of Ra at Nρ = 6 (circles) and Nρ = 2 (triangles). (b) The same quantity as a function of Nρ at Ra = 32. The averages were performed over dipolarity time series such as those in Figure 3.
In the text

Fig. 9. (a) Instantaneous snapshots of the radial magnetic field Br (contours) and radial velocity component vr (color) for Nρ = 2, Ra/Racr = 32 at r/ro = 0.9. The contour levels are [0.1Brmin, 0.05Brmin, 0.05Brmax, 0.1Brmax]. (b) The same for the stronger stratification, Ra/Racr = 32, Nρ = 6. (c) Two-dimensional probability density of the joint distribution of b2 and vr near the surface, r/ro = 0.9. Left, Nρ = 6; right, Nρ = 2. $\stackrel{\sim }{\mathit{Ra}}=32$. Colors are spaced logarithmically. (d) Unbiased sample skewness coefficient G1 (5) as a function of $\stackrel{\sim }{\mathit{Ra}}$. Red, Nρ = 6, blue, Nρ = 4, green, Nρ = 2.

In the text

	Fig. 10. (a) Poynting flux (10), integrated in latitude θ, as a function of spherical radius r and Nρ, $\stackrel{\sim }{\mathit{Ra}}=32$. (b) Poynting flux as a function of $\stackrel{\sim }{\mathit{Ra}}$, Nρ = 6. (c) Same but for Nρ = 2.
In the text

Fig. 11. (a,b) Components of the mean-field γ-coefficient calculated for Nρ = 4, $\stackrel{\sim }{\mathit{Ra}}=32$ using the method of snapshots (Simard et al. 2016). (a) Radial component, γr. (b) Latitudinal component, γθ. (c) Total magnetic energy near the surface (r/ro = 0.9) for run Nρ = 6, $\stackrel{\sim }{\mathit{Ra}}=32$ as a function of spherical harmonic degree l for high- (tH) and low-energy (tL) dynamo states.

In the text