# What Sets the Magnetic Field Strength and Cycle Period in Solar-type Stars?

G. Guerrero, B. Zaire, P. K. Smolarkiewicz, E. M. de Gouveia Dal Pino, A. G. Kosovichev, N. N. Mansour

► To cite this version:

G. Guerrero, B. Zaire, P. K. Smolarkiewicz, E. M. de Gouveia Dal Pino, A. G. Kosovichev, et al. What Sets the Magnetic Field Strength and Cycle Period in Solar-type Stars?. The Astrophysical Journal, 2019, 880, 10.3847/1538-4357/ab224a. insu-03674431

HAL Id: insu-03674431

https://insu.hal.science/insu-03674431v1

Submitted on 13 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# What Sets the Magnetic Field Strength and Cycle Period in Solar-type Stars?

G. Guerrero<fnref n="1" />, B. Zaire<fnref n="2" />,<fnref n="1" />, P. K. Smolarkiewicz<fnref n="3" />, E. M. de Gouveia Dal Pino<fnref n="4" />, A. G. Kosovichev<fnref n="5" />, and N. N. Mansour<fnref n="6" />

Received 2018 October 17; revised 2019 May 2; accepted 2019 May 9; published 2019 July 17

## Abstract

Two fundamental properties of stellar magnetic fields have been determined by observations for solar-like stars with different Rossby numbers (Ro), namely, the magnetic field strength and the magnetic cycle period. The field strength exhibits two regimes: (1) for fast rotation, it is independent of Ro, and (2) for slow rotation, it decays with Ro following a power law. For the magnetic cycle period, two regimes of activity, the active and inactive branches, have also been identified. For both of them, the longer the rotation period, the longer the activity cycle. Using global dynamo simulations of solar-like stars with Rossby numbers between ~0.4 and ~2, this paper explores the relevance of rotational shear layers in determining these observational properties. Our results, consistent with nonlinear $\alpha^2\Omega$ dynamos, show that the total magnetic field strength is independent of the rotation period. Yet at surface levels, the origin of the magnetic field is determined by Ro. While for Ro $\lesssim$ 1, it is generated in the convection zone, for Ro $\gtrsim$ 1, strong toroidal fields are generated at the tachocline and rapidly emerge toward the surface. In agreement with the observations, the magnetic cycle period increases with the rotational period. However, a bifurcation is observed for Ro ~ 1, separating a regime where oscillatory dynamos operate mainly in the convection zone from the regime where the tachocline has a predominant role. In the latter, the cycles are believed to result from the periodic energy exchange between the dynamo and the magneto-shear instabilities developing in the tachocline and the radiative interior.

**Key words:** dynamo – stars: magnetic field – stars: rotation – stars: solar-type

**Supporting material:** animations

## 1. Introduction

Modern observations have revealed the existence of large-scale magnetic fields in most types of stars across the H-R diagram. Among these are solar-type stars, with convective envelopes, large-scale magnetic fields, as well as fully convective stars, which are either in a pre-main-sequence phase or represent main-sequence M-type dwarfs. There is no doubt that the magnetic field is relevant in every phase of the life of stars. It also plays a critical role in the evolution of planetary disks and, ultimately, may define criteria for habitability (do Nascimento et al. 2016).

For late-type and solar-like stars that have convective envelopes, large-scale magnetic fields, as well as different field topologies observed on their surface (Petit et al. 2008; Gregory et al. 2012), are convincing evidence of a dynamo mechanism operating in the stellar interiors. The dynamo is the result of a complex system of electric currents induced by the differential rotation in a processes known as the $\Omega$-effect and the helical turbulent convective motions and fields producing the so-called $\alpha$-effect (Parker 1955; Steenbeck et al. 1966). Furthermore, as will be detailed below, observations show clear correlations of the magnetic field strength and the activity cycle period with the stellar Rossby number, $Ro = P_{\text{rot}}/\tau_c$, where $P_{\text{rot}}$ is the period of rotation, and $\tau_c$ is the convective turnover time. These correlations provide information about the dynamo process that might help to decode its elusive details.

The relationship between the magnetic field strength and Ro has two well-defined regimes. They are evident in the stellar X-ray luminosity data, $L_X$ (Pizzolato et al. 2003; Wright et al. 2011), as

well as in direct measurements of the mean magnetic field, $\langle B \rangle$ (Vidotto et al. 2014). For $Ro \gtrsim 0.1$, the magnetic activity shows a power-law behavior, $\langle B \rangle \propto Ro^{-1.38}$ (Vidotto et al. 2014). For $Ro \lesssim 0.1$, observations indicate a regime of activity independent of Ro, which is often called the saturated phase. Recent observational results by Wright & Drake (2016) point out that the two regimes described above occur in both fully and partially convective stars. These results question the canonical theory in which a rotational shear layer at the interface between the radiative and the convection zones (CZs) is fundamental. This layer is called the tachocline.

A number of stars, especially of types F, G, and K, exhibit chromospheric variations consistent with cyclic magnetic activity (Baliunas et al. 1995; Saar & Brandenburg 1999; Brandenburg et al. 2017). The seminal studies of Noyes et al. (1984b) and Brandenburg et al. (1998) identified correlations between the magnetic cycles and the stellar rotation. Their results suggested the existence of two main branches, dividing active (A) from inactive (I) stars (Saar & Brandenburg 1999). The A and I branches show positive dependences between the ratio $P_{\text{rot}}/P_{\text{cyc}}$ (where $P_{\text{cyc}}$ is the magnetic cycle period) and either $Ro^{-1}$ or $\langle R'_{HK} \rangle$ (where $\langle R'_{HK} \rangle = F'_{HK}/F_{\text{bol}}$ is the mean fractional CaII H and K flux relative to the stellar bolometric flux, $F_{\text{bol}}$). The same branches, A and I, are confirmed by Brandenburg et al. (2017) using recalibrated measurements and new data from the *Kepler* satellite.

Böhm-Vitense (2007) compared the period of rotation with the magnetic cycle period ($P_{\text{rot}}$ versus $P_{\text{cyc}}$) of the Mount Wilson sample of stars and found the same two branches, A and I,

<footnote n="1">¹ Physics Department, Universidade Federal de Minas Gerais, Av. Antonio Carlos, 6627, Belo Horizonte, MG 31270-901, Brazil; guerrero@fisica.ufmg.br</footnote>

<footnote n="2">² IRAP, Université de Toulouse, CNRS/UMR 5277, CNES, UPS, 14 avenue E. Belin, Toulouse, F-31400 France</footnote>

<footnote n="3">³ European Centre for Medium-Range Weather Forecasts, Reading RG2 9AX, UK</footnote>

<footnote n="4">⁴ Astronomy Department, Universidade de São Paulo, IAG-USP, Rua do Matão, 1226, São Paulo, SP, 05508-090, Brazil</footnote>

<footnote n="5">⁵ New Jersey Institute of Technology, Newark, NJ 07103, USA</footnote>

<footnote n="6">⁶ NASA, Ames Research Center, Moffett Field, Mountain View, CA 94040, USA</footnote>
having different positive slopes. Brandenburg et al. (2017) claimed that each of the trends found under this representation describes a family of lines for different values of the convective turnover time rather than a universal behavior. The comparison of $P_{\text{rot}}/P_{\text{cyc}}$ versus $\langle R'_{HK} \rangle$ provides an universal trend given that all of the quantities are observables and do not depend on the unknown convective turnover time, $\tau_c$.

The correlations between the magnetic field amplitude and the magnetic cycle period with the Rossby number have challenged theoreticians and modelers over the last decades. For the scaling of the field strength with Ro, explanations rely on the so-called mean-field dynamo number $D = C_\alpha C_\Omega$, where $C_\alpha$ and $C_\Omega$ are nondimensional quantities that compare the inductive effects of the turbulent $\alpha$-effect and the shear against the dissipative effects of turbulence (Noyes et al. 1984b). Nevertheless, the scaling of the dynamo coefficients with the Rossby number is unknown, and the hypothesis based on the linear mean-field theory remains unproven (Noyes et al. 1984b; Saar & Brandenburg 1999; Blackman & Thomas 2015).

Under the same linear dynamo theory, the activity cycle period is proportional to $D^{-1/2}$ (Stix 1976; Noyes et al. 1984b). Yet this is an incomplete approach, since it does not consider the back reaction of the magnetic field on the flow. The mean-field simulations of Pipin & Kosovichev (2016) that consider the dynamic evolution of the $\alpha$-effect (i.e., a form of including the magnetic field back reaction on the flow) produce a magnetic cycle period that increases with the period of rotation, while the magnetic field amplitude decreases with the increase of the rotational period. These results are in agreement with the observations, except for the saturated phase, which was not considered in their model. When the nonlinearity is considered through a simple algebraic quenching, the opposite relation is obtained; i.e., $P_{\text{cyc}}$ decreases with the increase of $P_{\text{rot}}$. Compared with other mean-field results, Pipin & Kosovichev (2016) clearly demonstrated the importance of the nonlinear processes occurring in the dynamo. For instance, flux transport mean-field simulations, in which the cycle period is mainly determined by the meridional circulation and the buoyant rise of magnetic flux tubes, result in correlations that are at odds with the observations (Jouve et al. 2010; Karak et al. 2014). Similarly, the global numerical simulations of solar-like stars performed by Strugarek et al. (2017) and Warnecke (2018) obtained a $P_{\text{cyc}}$ decreasing with the increase of $P_{\text{rot}}$. More recently, Viviani et al. (2018) reported high-resolution simulations of stars with rotation between 1 and 30 times the solar rotation rate. In the $P_{\text{rot}}/P_{\text{cyc}}$ representation, they found that the slow-rotating cases, displaying antisolar differential rotation, fall close to the I branch but have a negative slope. Interestingly, the dynamo solutions in the fast-rotating cases are all nonaxisymmetric and fall in a different branch of activity for superactive stars. One common aspect of the global simulations above (Strugarek et al. 2017; Warnecke 2018; Viviani et al. 2018) is the absence of the radial shear layers that are well observed in the Sun and should also exist in stars with radiative zones. Thus, these models do not generate strong toroidal fields, neglecting the most important source of nonlinearity. It is then worth exploring the influence of these regions in the dynamo mechanism and the determination of stellar magnetic cycles.

In this paper, we study the scaling of the magnetic field strength and the magnetic cycle period with the Rossby number in global convective dynamo simulations including rotational shear layers. The numerical model employed here is the same as described in Guerrero et al. (2016a). A detailed analysis of

the angular momentum budget and the generation of torsional oscillations was presented in Guerrero et al. (2016b). In this paper, we present an extensive series of simulations where the only varying parameter is the rotation rate of the reference frame, and therefore Ro. Our previous results have demonstrated that the presence of tachoclines results in dynamos where the evolution of the plasma is governed in large extent by deep-seated magnetic fields. Here we show how the scaling laws obtained in these dynamos exhibit similarities with the observations. The goal of this paper is to provide a theoretical analysis explaining the physics behind the resulting scaling laws.

In the next section, we describe the numerical model, and the results are described in Section 3. We discuss the implications of our results for solar and stellar dynamos in Section 4. The technical details of the analysis are presented in Appendices A–C.

## 2. The Model

We consider a full spherical shell domain, $0 \le \phi \le 2\pi$, $0 \le \theta \le \pi$, with the bottom boundary at $r_b = 0.61R_\odot$ and the top boundary at $r_t = 0.96R_\odot$. The simulations have a grid resolution of $128 \times 64 \times 64$ points in longitude ($\phi$), latitude ($\theta$), and radius ($r$), respectively.

We solve a set of anelastic MHD equations in the following form:

$$ \nabla \cdot (\rho_s \mathbf{u}) = 0, \qquad (1) $$

$$ \frac{D\mathbf{u}}{Dt} + 2\Omega \times \mathbf{u} = -\nabla\left(\frac{p'}{\rho_s}\right) + g\frac{\Theta'}{\Theta_s} + \frac{1}{\mu_0\rho_s}(\mathbf{B} \cdot \nabla)\mathbf{B}, \quad (2) $$

$$ \frac{D\Theta'}{Dt} = -\mathbf{u} \cdot \nabla\Theta_e - \frac{\Theta'}{\tau}, \qquad (3) $$

$$ \frac{DB}{Dt} = (\mathbf{B} \cdot \nabla)\mathbf{u} - \mathbf{B}(\nabla \cdot \mathbf{u}), \qquad (4) $$

where $D/Dt = \partial/\partial t + \mathbf{u} \cdot \nabla$ is the total time derivative, $\mathbf{u}$ is the velocity field in a rotating frame with angular velocity $\Omega = (\Omega_r, \Omega_\theta, \Omega_\phi) = \Omega_0(\cos\theta, -\sin\theta, 0)$, $p'$ is the pressure perturbation variable that accounts for both the gas and magnetic pressure, $\mathbf{B}$ is the magnetic field, and $\Theta'$ is the potential temperature perturbation with respect to an ambient state $\Theta_e$ (see Guerrero et al. 2013; Cossette et al. 2017, for comprehensive discussions). Furthermore, $\rho_s$ and $\Theta_s$ are the density and potential temperature of the reference state, which is chosen to be isentropic (i.e., $\Theta_s = \text{const}$) and in hydrostatic equilibrium; $g = GM/r^2\dot{e}_r$ is the gravity acceleration, where G and M are the gravitational constant and stellar mass, respectively, and $\mu_0$ is the magnetic permeability. The potential temperature, $\Theta$, is related to the specific entropy: $s = c_p \ln \Theta + \text{const}$.

The simulations were performed using the EULAG-MHD code,<fnref n="7" /> a spin-off of the hydrodynamical model EULAG predominantly used in atmospheric and climate research (Prusa et al. 2008). The time evolution is calculated using a bespoke semi-implicit approach derivable from the trapezoidal-rule path integration of the prognostic Equations (2)–(4). At the heart of the approach, there is a nonoscillatory (viz. high-resolution) forward-in-time Multidimensional Positive Definite Advection

<footnote n="7">⁷ The code is available at the dedicated website: http://www.astro.umontreal.ca/~paulchar/grps/eulag-mhd.html.</footnote>
<table><caption>**Table 1**
Simulation Parameters and Results</caption><thead><tr><th>Model</th><th>P<sub>rot</sub><br>(days)</th><th>⟨u<sub>rms</sub>⟩<sub>CZ</sub><br>(m s<sup>-1</sup>)</th><th>τ<sub>c</sub><br>(s)</th><th>Ro</th><th>Ro<sub>1</sub></th><th>D'<sub>10<sup>3</sup></sub></th><th>P<sub>cyc</sub><br>(yr)</th><th>⟨B<sub>φ</sub>⟩<sub>TAC</sub><br>(T)</th><th>⟨B<sub>φ</sub>⟩<sub>CZ</sub><br>(T)</th><th>⟨B<sub>φ</sub>⟩<sub>NSL</sub><br>(T)</th><th>⟨B<sub>p</sub>⟩<sub>TAC</sub><br>(T)</th><th>⟨B<sub>p</sub>⟩<sub>CZ</sub><br>(T)</th><th>⟨B<sub>p</sub>⟩<sub>NSL</sub><br>(T)</th></tr></thead><tbody><tr><td>RC07</td><td>7.0</td><td>35.18</td><td>2.71</td><td>0.36</td><td>0.26</td><td>2.93</td><td>3.0</td><td>0.095</td><td>0.063</td><td>0.072</td><td>0.241</td><td>0.035</td><td>0.024</td></tr><tr><td>RC14</td><td>14.0</td><td>37.69</td><td>3.41</td><td>0.56</td><td>0.60</td><td>3.74</td><td>9.6</td><td>0.167</td><td>0.128</td><td>0.152</td><td>0.483</td><td>0.045</td><td>0.039</td></tr><tr><td>RC18</td><td>18.0</td><td>38.99</td><td>2.83</td><td>0.87</td><td>0.80</td><td>1.84</td><td>...</td><td>0.215</td><td>0.044</td><td>0.062</td><td>0.302</td><td>0.024</td><td>0.023</td></tr><tr><td>RC21</td><td>21.0</td><td>38.85</td><td>2.87</td><td>1.01</td><td>0.91</td><td>0.99</td><td>30.1</td><td>0.246</td><td>0.037</td><td>0.062</td><td>0.240</td><td>0.022</td><td>0.024</td></tr><tr><td>RC24</td><td>24.0</td><td>38.83</td><td>2.96</td><td>1.12</td><td>1.02</td><td>0.35</td><td>16.3</td><td>0.286</td><td>0.058</td><td>0.088</td><td>0.110</td><td>0.040</td><td>0.041</td></tr><tr><td>RC28</td><td>28.0</td><td>38.90</td><td>3.11</td><td>1.24</td><td>1.19</td><td>-0.27</td><td>16.3</td><td>0.318</td><td>0.064</td><td>0.099</td><td>0.106</td><td>0.040</td><td>0.044</td></tr><tr><td>RC35</td><td>35.0</td><td>39.84</td><td>3.37</td><td>1.43</td><td>1.51</td><td>-1.06</td><td>19.1</td><td>0.409</td><td>0.066</td><td>0.087</td><td>0.111</td><td>0.034</td><td>0.037</td></tr><tr><td>RC42</td><td>42.0</td><td>40.70</td><td>3.64</td><td>1.59</td><td>1.87</td><td>-1.54</td><td>22.9</td><td>0.353</td><td>0.068</td><td>0.088</td><td>0.106</td><td>0.029</td><td>0.034</td></tr><tr><td>RC49</td><td>49.0</td><td>42.19</td><td>3.83</td><td>1.76</td><td>2.24</td><td>-1.92</td><td>26.2</td><td>0.456</td><td>0.075</td><td>0.090</td><td>0.106</td><td>0.025</td><td>0.028</td></tr><tr><td>RC56</td><td>56.0</td><td>43.83</td><td>3.95</td><td>1.95</td><td>2.69</td><td>-2.18</td><td>...</td><td>0.406</td><td>0.079</td><td>0.101</td><td>0.169</td><td>0.022</td><td>0.026</td></tr><tr><td>RC63</td><td>63.0</td><td>44.07</td><td>4.08</td><td>2.12</td><td>3.03</td><td>-2.28</td><td>...</td><td>0.340</td><td>0.079</td><td>0.089</td><td>0.149</td><td>0.021</td><td>0.023</td></tr></tbody></table>

**Note.** The convective turnover time, $\tau_c$, is computed from the spectra of the nonaxisymmetric velocity and magnetic fields, as explained in Appendix A. Here $R_0 = P_{rot}/2\pi\tau_c$ is the Rossby number evaluated with $\tau_c$ in this table; $R_0 = P_{rot}/\tau_c^*$ is computed with $\tau_c^*$ estimated at one pressure height scale above the bottom of the convection zone, as in Noyes et al. (1984a). The dynamo number, $D_c' = C_a C_{\phi}^{\prime r}$, corresponds to the average over the TAC in the polar region. The primes in this definition mean that it is computed from polynomial fits to the dynamo coefficients. The period is computed by using the Fourier transform, as explained in Appendix C. The quantities in angular brackets correspond to the time and volume average over the regions TAC, CZ, and NSL (see the text).

Transport Algorithm (MPDATA) broadly documented in the literature (see Smolarkiewicz 2006 for an overview and Waruszewski et al. 2018 for recent advancements). A comprehensive description of the MHD implementation is presented in Smolarkiewicz & Charbonneau (2013).

The truncation terms in MPDATA evince viscosity comparable to the explicit subgrid scale (SGS) viscosity used in large-eddy simulation (LES) models (Elliott & Smolarkiewicz 2002; Domaradzki et al. 2003; Margolin et al. 2006). Thus, the results of MPDATA are often interpreted as implicit LES, or ILES (Smolarkiewicz & Margolin 2007). This implicit SGS approach has been fundamental to successfully reproducing the solar tachocline and deep-seated magnetic dynamos with timescales compatible with the solar cycle (Ghizaru et al. 2010; Racine et al. 2011; Smolarkiewicz & Charbonneau 2013; Guerrero et al. 2013, 2016a, 2016b).

For the velocity field, we use impermeable, stress-free conditions at the top and bottom surfaces of the shell, whereas the magnetic field is assumed to be radial at these boundaries. Finally, for the thermal boundary condition, we consider zero radial derivative of the radial convective flux of potential temperature perturbations at the bottom and zero convective radial flux of the potential temperature perturbations at the top surface. All simulations start from a random noise, centered about zero and the same for each experiment, in the potential temperature perturbations, velocity, and magnetic field. For the vector fields the noise is divergence-free. All simulations are run until reaching a statistically steady state, using a constant time step $\Delta t = 1800$ s. A list of the simulation runs used in the current paper is presented in Table 1.

## 3. Results

### 3.1. Large-scale Flows and Magnetic Field

Figure 1 shows the differential rotation (left panel of each model) and meridional circulation (right panels) of some representative models from (a) RC07 to (i) RC63 (the profiles of models RC28 and RC56 are not presented here, since they appear in Guerrero et al. 2016a, 2016b). In the differential rotation profiles, colored contours depict the variations of the mean angular velocity, $\bar{\Omega}$, calculated as a temporal and azimuthal

average, with respect to the rotating frame. The results qualitatively show that the gradients of angular velocity become prominent with the increase of the rotation period. Observational results also indicate that the latitudinal differential rotation increases with the rotation period (see Section 6.1 of Lehtinen et al. 2016, and references therein). The radial differential rotation for stars other than the Sun is evasive to observations. Nevertheless, it is commonly assumed that the shear is stronger for rapidly rotating stars (e.g., Noyes et al. 1984b; Blackman & Thomas 2015). However, in our fast-rotating simulations (RC07–RC21), there is almost no radial shear at the tachocline and in the near-surface layer. As the rotation diminishes progressively from model RC21 to model RC63, the clearest gradients are observed at the tachocline, as well as in the near-surface layers (a quantitative analysis is presented in Section 3.2).

In the meridional circulation panels, the colored contours show the mean latitudinal velocity ($\bar{u}_{\theta}$). The contour lines show the stream function. For the fast-rotating models, we observe a pattern of multiple slow convection cells circulating over thin cylinders. The number of cells decreases with the increase of the rotation period. For instance, from model RC28 to model RC63, only two cells are developed in each meridional quadrant. One is a broad counterclockwise cell going from $\sim 0.72$ to $\sim 0.86R_\odot$ and in all latitudes. The second one is clockwise. It is located above $\sim 0.86R_\odot$ and close to the equator. The amplitude of $\bar{u}_{\theta}$ is monotonically increasing from the faster- (RC07) to the slower- (RC63) rotating models, as can be noticed in the corresponding color bars.

Figure 2 shows the time–latitude (at $r = 0.95R_\odot$) and time–radius (at $\theta = 24^\circ$ latitude) evolution of the mean magnetic field associated with the mean flows described above ($\bar{B}_r$ with colored contours and $\bar{B}_\phi$ with contour lines). Several dynamo modes can be distinguished for different values of $R_0$. For the fastest-rotating model (smallest $R_0$), RC07, the rotation profile is almost homogeneous in the whole domain, and the shear is small (Figure 1(a)). Even though there is a steady magnetic field in the radiative zone, an oscillatory dynamo with a short period of 3 yr develops predominantly in the upper convection zone. In model RC14 (panel (b)), the radial shear at the tachocline starts to develop, leading to the formation of an antisymmetric steady dynamo. In this model, the magnetic field is
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Differential rotation (left) and meridional circulation (right) of models (a) RC07 to (i) RC63. The colors in the differential rotation show isorotation contours in the rotating frame. The solid lines divide the domain in the six different analysis regions. In the meridional circulation panels, the colored contours show the latitudinal velocity, $\bar{u}_{\theta}$. The solid (dashed) lines represent clockwise (counterclockwise) circulation.</figcaption></figure>

oscillatory in the convection zone; however, it does not show hemispheric polarity reversals. In model RC18 (panel (c)), the dynamo also develops mainly at the tachocline. The field is steady, and no cyclic variations of the field are observed in the convection zone. In model RC21 (panel (d)), the solution exhibits bimodality; i.e., two dynamo modes are simultaneously excited with the magnetic field periodic in the equatorial region but steady at the poles. In the radius–time diagram, periodic reversals of the toroidal field (solid and dashed lines) can be observed while the poloidal field remains steady. Models RC24–RC49 (panels (e)–(g)) are periodic with well-defined magnetic field polarity reversals. The magnetic field is generated mainly at the tachocline, but the dynamo action occurs in the entire convection zone. The radial magnetic field is also observed in the convectively stable layer and reverses polarity during the toroidal field maxima. The parity of these models varies in time and is not well defined. The slowly rotating models RC56 (the butterfly diagram of which is not shown in Figure 2 but is presented in Figure 6(c) of Guerrero et al. 2016a) and RC63 (panel (h)) are all antisymmetric steady dynamos. We have also performed simulations (not shown here) for longer rotational periods, i.e.,

112 and 140 days. In these cases, there is still dynamo action, but the magnetic field is weak, with no back reaction on the convection zone dynamics. For models with a rotation period longer than 224 days, the dynamo instability does not develop.

Complementary to Figure 2, Figure 3 depicts the latitudinal distribution of the fields in the meridional plane (r–θ). It highlights the series of animations available with the figure. On the left of each panel, the line integral convolution (LIC) representation depicts the distribution of the poloidal field lines, with the color indicating the magnitude and direction of the mean latitudinal field, $\bar{B}_{\theta}$. On the right of each panel, the colored contours show the distribution of the mean toroidal magnetic field, $\bar{B}_{\phi}$. They make clear that a layer of strong toroidal field is formed at and below the tachocline, especially for models with $P_{\text{rot}} \gtrsim 21$ days.

Figures 2 and 3 evidence the complexity of large-scale dynamos. Different dynamo modes can be excited, depending on local conditions. Examples of this are the cases RC07–RC21, where the steady and oscillatory modes are mixed. In the time–latitude diagrams of models RC24–RC49, it can be seen that the field generated at the base of the convection zone is superposed with
<figure><img src="image_4.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_11.png" /><figcaption>**Figure 2.** Time–latitude (at $r = 0.95R_⊙$) and time–radius (at $\theta = 24°$) butterfly diagrams showing the evolution of models between (a) RC07 and (h) RC63. The colored contours show the radial magnetic field with its amplitude in teslas depicted in the color tables. The solid (dashed) contour lines show the positive (negative) toroidal field.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 3.** Snapshots indicating the distribution of the mean magnetic field components for some characteristic simulations. On the left, the LIC representation depicts the distribution of the meridional field lines in the convection zone, with the color indicating the magnitude of the mean latitudinal field, $\bar{B}_{\phi}(t, \theta, r)$. The colored contours on the right correspond to the azimuthal mean magnetic field, $\bar{B}_{\phi}(t, \theta, r)$. Because the fields at the tachocline are about 1 order of magnitude larger than the fields at the convection zone, the contours are highly saturated. For a better comprehension of these results and behavior, an animation showing all panels is available. The video duration is 13 s. (Full animations of each panel are also available at http://lilith.fisica.ufmg.br/~guerrero/cycle_global.html.) In the movie, especially for the upper regions of simulations with periods larger than 28 days, we note that the LIC part seems to evolve faster. This is due to the way in which the LIC is generated from the convolution between the vector field and a random white-noise background, which causes an artificial advance with respect to the rapid timescale of the magnetic flux emergence in the frames.</figcaption></figure>

(An animation of this figure is available.)

the field generated near the surface, forming irregular branches that end up mixing the parity of the dynamo solutions. Because the simulated magnetic field is generated by sources at different locations, in the next section, we explore the magnetic field amplitude and dynamo coefficients by performing volume averaging over three radial regions: the tachocline region (TAC), in the range $0.63R_{\odot} \le r < 0.74R_{\odot}$; the bulk of the convection zone (CZ), in the range $0.74R_{\odot} \le r < 0.89R_{\odot}$; and the near-surface layer region (NSL), in the range $0.89R_{\odot} \le r < 0.96R_{\odot}$. In latitude, we separate the domain into polar (POL) and equatorial (EQU) regions. The averages over these regions are denoted by angle brackets, $\langle \rangle$.

## 3.2. Mean-field Analysis: Magnetic Field and Dynamo Coefficients

For a better understanding of the dynamo solutions as a function of the Rossby number, we perform a systematic analysis of the simulation results in terms of the mean-field dynamo framework. Upon the condition that the results are axisymmetric, the magnetic field, as well as the velocity field, may be decomposed into their large-scale and turbulent components. This analysis leads to the governing equation for a mean-field $\alpha^2\Omega$ dynamo (Moffatt 1978),

$$ \frac{\partial \bar{B}}{\partial t} = \nabla \times (\bar{u} \times \bar{B}) + \nabla \times (\alpha \bar{B}) - \nabla \times (\eta \nabla \times \bar{B}), \quad (5) $$

where $\bar{B} = (\bar{B}_r, \bar{B}_{\theta}, \bar{B}_{\phi})$ is the magnetic field averaged over longitude and $\bar{u} = (\bar{u}_r, \bar{u}_{\theta}, \bar{u}_{\phi})$ is the velocity field averaged over

longitude and time in an interval during the dynamo saturated phase (see Figure 1). In the second term on the right-hand side (rhs), the $\alpha$ term stems for the $\alpha$-effect, which has kinetic and magnetic contributions, namely, $\alpha = \alpha_k + \alpha_m$. These terms generate the large-scale magnetic field from small-scale helical motions and currents, respectively. In the third term, $\eta = \eta_m + \eta_t$ is the sum of the molecular and turbulent magnetic diffusivities.

If the mean-velocity field is expressed as $\bar{u} = r \sin\theta\Omega\hat{e}_{\phi} + \bar{u}_p$, with $\bar{u}_p = (\bar{u}_r, \bar{u}_{\theta}, 0)$, and the mean magnetic field as $\bar{B} = \bar{B}_{\phi} + \bar{B}_p$, with $\bar{B}_p = (\bar{B}_r, \bar{B}_{\theta}, 0)$, then Equation (5) can be written as

$$ \frac{\partial \bar{B}}{\partial t} = [r \sin \theta \mathbf{B}_p \cdot \nabla \Omega] + \nabla \times (\bar{u}_p \times \bar{B}) + \nabla \times (\alpha \bar{B}) - \nabla \times (\eta \nabla \times \bar{B}). \quad (6) $$

Here the first term on the rhs represents the rotational shear that generates the toroidal field from the poloidal one. The second term corresponds to the advective transport by the meridional circulation, $\bar{u}_p$. The $\alpha$-effect, the third term in the above equation, represented by a pseudoscalar, $\alpha$, operates on both components of the field. Rigorously, $\alpha$ is a second-order tensor that acts as a source of the toroidal and poloidal fields and also advects them (Moffatt 1978). The turbulent magnetic diffusivity, $\eta_t$, comes from a third-order tensor, $\beta$, which also might have source terms (Brandenburg et al. 2008). For our analysis, we estimate $\alpha$ and $\eta_t$ using the first-order smoothing approximation (FOSA), as detailed in Appendix B. This approximation assumes isotropy such that both $\alpha$ and $\beta$ become scalars. We notice that this approximation is valid only
<figure><img src="image_13.png" /><figcaption>**Figure 4.** Top panels: toroidal (red) and poloidal (blue) magnetic field energy densities as a function of Ro. The left, middle, and right columns show averages over the TAC, CZ, and NSL regions, respectively. In all panels, the solid (dashed) lines depict the values at the EQU (POL) latitudes. The magnetic energy densities are normalized to the kinetic energy density of each model averaged over the entire volume, $e_k$. The black lines in the top panels correspond to the magnetic energy density averaged over the entire northern hemisphere for the TAC, CZ, and NSL regions. The second, third, and fourth rows correspond to the dynamo coefficients, $C_\alpha$, $C_\Omega$, and $C_u$ respectively, as a function of Ro. The thin solid and dashed black lines in the second and third rows show the polynomial fits to these coefficients. The bottom panels show the dynamo numbers, $D'_r$ (black) and $D'_\theta$ (red). The shaded area in the left panel exhibits the regions in the Ro-$D'_r$ space where the dynamos are oscillatory at the tachocline. In the right panel, the shaded area shows the range of Ro for which the magnetic energy densities decay. Finally, the black dotted lines indicate the zero level of each quantity.</figcaption></figure>

for low magnetic Reynolds numbers.<fnref n="8" /> However, it has been shown that the coefficient profiles obtained with FOSA are qualitatively compatible with those obtained by directly inverting the electromotive force (Racine et al. 2011) and those obtained by the test-field method (Warnecke et al. 2018). In spite of possible differences with their actual values, the systematic use of the same technique for all of the simulations

provides a reliable picture of the change of the helicities with the Rossby number. Details of the computation of $\alpha_k$ and $\alpha_m$, as well as of the turbulent diffusion coefficient, are presented in Appendix B. Meridional profiles of the kinetic, magnetic, and total $\alpha$-effect for representative models between RC07 and RC63 are shown in Figures 13(a)–(h). The radial profile of $\eta_t$ for the same models is presented in Figure 12(d).

Writing Equation (6) in a nondimensional form and using a characteristic dynamo timescale, $\tau_{dyn} = R_\odot^2/\eta_0$, where $\eta_0$ is a suitable value of the diffusivity coefficient, we can define

<footnote n="8">⁸ A complete determination of the dynamo coefficients can be performed via the so-called test-field method (Warnecke et al. 2018).</footnote>
nondimensional dynamo coefficients that compare the inductive, $C_α = ατ_{dyn}/R_⊙$ and $C_Ω = ΔΩτ_{dyn}$, and advective, $C_u = \bar{u}_pτ_{dyn}/R_⊙$, effects with diffusion. Here $\bar{u}_p = \sqrt{\bar{u}_r^2 + \bar{u}_θ^2}$ is the amplitude of the meridional motions. Because both the radial and latitudinal derivatives of $\bar{Ω}$ contribute to the generation of $\bar{B}_φ$ (by stretching the poloidal field lines in the azimuthal direction), the parameter $C_Ω$ has two components, $C_{Ωr} = R_⊙∂_rΩτ_{dyn}$ and $C_{Ωθ} = R_⊙(∂_θΩ/r)τ_{dyn}$.

The energy density of the magnetic field components, as well as the dynamo coefficients as a function of the Rossby number, Ro, are presented in Figure 4. The left, middle, and right columns correspond to the TAC, CZ, and NSL regions, respectively. In each panel, the solid and dashed lines depict the volume averages over the POL and EQU regions and the corresponding radial extent. The first row shows the poloidal ($e_{B_p} = \bar{B}_p^2/2μ_0$; blue lines and symbols) and toroidal ($e_{B_θ} = \bar{B}_θ^2/2μ_0$; red) magnetic energy densities. They are normalized to the kinetic energy density, $e_k = ⟨ρ_k⟩u_{rms}^2/2$, where $⟨ρ_k⟩$ is the mean value of the isentropic density, and $u_{rms}$ is the time-volume averaged rms velocity. The second, third, and fourth rows show, respectively, $C_α$, $C_Ω$, and $C_u$. The fifth row depicts the radial and latitudinal dynamo numbers, $D'_r = C'_α C^r_Ω$ and $D'_θ = C'_α C^θ_Ω$. The primes in these definitions come from the fact that the dynamo numbers are computed from polynomial fits to the dynamo coefficients. We next describe the main characteristics of these quantities from the TAC to the NSL regions.

### 3.2.1. TAC Region

In the TAC region, at both POL and EQU altitudes, the magnetic energy shows two separate branches. For $Ro \lesssim 1$, the poloidal field energy is larger than the toroidal one. Both trends cross at $Ro \sim 1$, and the toroidal field energy becomes larger for $Ro \ge 1$. The total magnetic energy averaged in the entire hemisphere (solid black line) is, thus, roughly constant for all of the Rossby numbers. It reaches roughly 50% of the kinetic energy. Note that most of the magnetic energy is in the tachocline magnetic fields; therefore, the normalized magnetic energy density, averaged over the entire domain, is independent of $Ro$. For $Ro \lesssim 1$, the $α$-effect is negative, and its amplitude decreases as $Ro$ increases. The meridional profiles in Figures 13(a)–(h) indicate that this term has a magnetic origin. In the EQU region, $C_α$ reaches zero at $Ro \sim 0.5$, while in the POL region, it does this at $Ro \sim 1$, roughly the same point where the toroidal and magnetic energy densities have similar values. Note that in Figures 2(a)–(c), corresponding to the same range of $Ro$, the radial field in the stable layer is positive. For $Ro > 1$, $C_α$ increases in POL and is roughly null in EQU. Regarding $C^r_Ω$, it is negative in POL and positive in EQU. It means that the rotation goes from faster to slower at POL latitudes and from slower to faster at EQU latitudes. As is observed in the Sun, in all simulations the largest shear occurs in the POL region. In both latitudinal zones, the radial shear increases with $Ro$. The latitudinal shear, $C^θ_Ω$, is rather similar in both latitudinal zones. It is positive and also increases slowly with $Ro$. The meridional flow coefficient, $C_u$, has small values at all latitudes and is roughly independent of $Ro$.

Although the trends are clear, the curves of the dynamo coefficients fluctuate. We found it useful and cleaner to plot the dynamo numbers, a multiplication of $C_α$ and $C_Ω$, by using polynomial fits to these coefficients. The radial (black) and

latitudinal (red) dynamo numbers are presented in the bottom row of Figure 4. In the TAC region, it is evident that the most relevant dynamo number is $D'_r$ at the POL latitudes (dashed black line). It is positive for the fast-rotating cases and decreases with the Rossby number until $Ro \sim 1$, where it becomes negative. Its amplitude increases with $Ro$ for the slow-rotating cases. In the EQU zone (solid black line), the radial dynamo number is positive (except for $Ro \sim 0.5$). It has small values because of small $C^{EQU}_α$ and is roughly independent of $Ro$. The latitudinal dynamo number, $D'_r$, at POL is negative for the fast-rotating cases and positive for the models with $Ro \gtrsim 1$. In the EQU region, it is roughly null. The shaded region shows the interval of Rossby numbers for oscillatory dynamos, $1 \lesssim Ro \lesssim 1.7$ (cases RC21, which exhibits bimodality, to RC49).

In Figure 5, we explore how the dynamo sources contribute to the spatiotemporal evolution of the mean magnetic fields. We compare the source terms in Equation (6) with the dynamo-generated mean magnetic fields in simulations (a) RC21, (b) RC35, and (c) RC49. All of the quantities are averaged in longitude and over the radial extent of the TAC region, and only the northern hemisphere is presented for clarity. The first row shows a time-latitude butterfly diagram with $\bar{B}_r$ presented in colored contours and $\bar{B}_φ$ with solid and dashed contour lines, as in Figure 2. The second and third row compare, respectively, the shear term, $r \sin θ(\bar{B}_p \cdot \nabla)Ω$, and the azimuthal component of the $α$ source term, $(\nabla \times α\bar{B})|_φ$, with the toroidal field, $\bar{B}_φ$. The fourth row depicts the colored contours of the radial component of the $α$ source term, $(\nabla \times α\bar{B})|_r$, and the contour lines of $\bar{B}_r$. The fifth row shows the time-latitude evolution of $\bar{B}_φ$ plotted as contour lines over the colored contours of $α = α_k + α_m$.

From the second row of Figure 5, it can be noticed that positive (negative) values of the shear source terms correlate well with positive (negative) values of the toroidal field at low latitudes in the three presented cases. The amplitude of the shear source increases with the rotational period (from the left to the right panel). As a matter of fact, it is the amplitude of the equatorial shear that makes the dynamo of simulation RC21 oscillatory. In simulations RC35 and RC49, the shear also correlates with the toroidal field near the pole, above 85°. As will be discussed in Section 3.4, for these cases, the reversal of the toroidal field starts at these latitudes. At intermediate latitudes, the signs of the shear source and the toroidal field are opposite. The panels in the third row reveal that the $α$ source term is principally responsible for the generation of the toroidal field, near the poles for simulation RC21 and in latitudes between 30° and 80° for RC35 and RC49. In the latter cases, representative of the slow-rotating simulations, its contribution increases in amplitude and spatial extent from the faster to the slower rotational rate. The fourth row of panels shows a clear correlation between the radial component of the $α$ source term and the radial field. (We have verified that in simulation RC21 and other simulations with faster rotation rates, the $α$ source term correlates better with the latitudinal field, $\bar{B}_θ$. However, since $\bar{B}_θ$ changes sign within the TAC region, this correlation is observed by averaging over smaller radial extents. For consistency, we have decided to present only the correlations between the $α$ source term and $\bar{B}_r$.) In simulations RC35 and RC49, the quantity $(\nabla \times α\bar{B})|_r$ is concentrated at the POL latitudes, from where it migrates equatorward. It induces a change of polarity of the radial field that follows the same pattern of migration until $(\nabla \times α\bar{B})|_r$ reverses sign.
<figure><img src="image_14.png" /><figcaption>**Figure 5.** Mean-field dynamo sources in the TAC region compared with the toroidal and radial mean magnetic fields for models (a) RC21, (b) RC35, and (c) RC49. All quantities are averaged in longitude and over the radial TAC extent. From top to bottom, the panels display the radial, $\bar{B}_r$ (colored contours), and toroidal, $\bar{B}_\phi$ (contour lines), mean magnetic fields; the sources of the toroidal field, $r \sin\theta(\bar{B}_p \cdot \nabla)\Omega$ and $(\nabla \times \alpha\bar{B})|_\phi$ (colored contours), both compared with $\bar{B}_\phi$; the source of the radial field, $(\nabla \times \alpha\bar{B})|_r$ (color), compared with $\bar{B}_r$ (contour lines); and $\alpha = \alpha_k + \alpha_m$ (colored contours), compared with $\bar{B}_\phi$. In the color map, the dimensions are [T] for the magnetic field, $10^{-8}$ [T s⁻¹] for the source terms, and [m s⁻¹] for the α-effect. The red vertical lines in the bottom right panel indicate the time interval described in Section 3.4.</figcaption></figure>

According to this analysis, we can infer that in this region of the domain, where the strongest magnetic field is generated, the mean-field coefficients capture well the physics of the dynamo mechanism. Thus, we can conclude that the dynamos operating in the simulations are of $\alpha^2\Omega$ type, with the α-effect generated in the stable layer. The amplitude of this quantity is not constant in time but varies dynamically with the cycle evolution, as can be seen in the bottom panels of Figure 5.

### 3.2.2. CZ Region

In the CZ region (middle column of Figure 4), the toroidal magnetic energy is larger than the poloidal for all values of Ro and seems to be independent of it. On the other hand, for Ro $\gtrsim$ 1.2, the poloidal field energy decays in all latitudes. The total magnetic energy (black dotted line) in this part of the domain is about 3% of the kinetic energy. Here $C_\alpha$ is positive at

EQU and consistent with zero at POL. In this case, $\alpha$ is mainly due to the kinetic helicity of the convective motions.

The radial shear, $C_\Omega^r$, is negative at POL with a minimum for Ro $\sim$ 1 (the same as in the TAC). At lower latitudes, it is positive (negative) for the fast- (slow-) rotating cases. The latitudinal shear, $C_\Omega^\theta$, is positive for all Ro with larger amplitudes in the POL latitudes. The meridional flow term, $C_\mu$, for Ro $\gtrsim$ 1 has values comparable to (or even larger than) $C_\alpha$, especially at EQU latitudes. Therefore, it is likely that this motion plays a significant role in transporting the magnetic field inside the CZ. Finally, the dynamo numbers, $D_r'$ and $D_\theta'$, are roughly zero at POL due to the small values of $C_\alpha$. At EQU, $D_\theta'$ follows the profile of $C_\alpha$, while $D_r'$ follows $C_\Omega^r$.

In contrast to the analysis of the TAC region, the interpretation of Figure 4 in the CZ region (top middle panel) in terms of mean-field theory is less conclusive. For Ro $\lesssim$ 1, an analysis similar to that of Figure 5 (not presented here for the sake of brevity) indicates local dynamo action. However, for
<figure><img src="image_15.png" /><figcaption>**Figure 6.** Radial profiles of $\alpha$ (blue lines), $\partial_r\Omega$ (red lines), and $\eta_t$ (black lines) averaged over time and over POL (dashed lines) and EQU (solid lines) for simulations (a) RC28, (b) RC42, (c) RC49, and (d) RC63. The black dotted line depicts zero. For clarity, the profiles of $\alpha$, $\partial_r\Omega$, and $\eta_t$ from all simulations were normalized to $10 \text{ m s}^{-1}$, $2 \times 10^{14} \text{ 1/(ms)}$, and $2 \times 10^9 \text{ m}^2 \text{ s}^{-1}$, respectively.</figcaption></figure>

Ro $\gtrsim$ 1, no clear correlation between the source terms and the magnetic field is observed. For instance, Figure 4 shows that the poloidal field decays in the POL and EQU regions, while concomitantly, its source term, $C_\alpha$, is either zero or $\sim$3, respectively. As this appears counterintuitive, we note that the local magnetic diffusivity $\eta_t$, as specified in Equation (11), generally varies between the pole and the equator. Because the coefficients in Figure 4 were computed using $\tau_{dyn} = R_\odot/\eta_0$, where $\eta_0$ is the average of the $\eta_t$ over the entire CZ, this variability is not reflected in the dynamo coefficients. For substantiation, Figure 6 displays normalized radial profiles of $\alpha$ (blue lines), $\partial_r\Omega$ (red lines), and $\eta_t$ (black lines) at POL (dashed lines) and EQU (solid lines) regions for models (a) RC28, (b) RC42, (c) RC49, and (d) RC63. The figure shows an orderly difference between the EQU and POL magnetic diffusivity increasing with Ro, which, in principle, could be responsible for different decay rates of the poloidal field energy at POL and EQU. Moreover, the time evolution of the mean magnetic field (see the movies in the supplementary material) reveals nonlocal contributions to the local magnetic field, i.e., magnetic buoyancy. The energy transported by nonlocal processes is not captured by Equation (6); therefore, it cannot be quantified by the dynamo coefficients. The intricacy of the magnetic fields in the CZ can also undergo the advective action of the $\alpha$-effect, the so-called turbulent pumping (Guerrero & de Gouveia Dal Pino 2008), or the meridional circulation.

### 3.2.3. NSL Region

The right column of Figure 4 presents the magnetic energy, dynamo coefficients, and numbers corresponding to the NSL. The top panel shows that while the toroidal field energy density at EQU increases with Ro, despite the fluctuations, at POL, it decreases. Since it is larger than the poloidal energy density, the total energy in the entire hemisphere is independent of Ro and about 5% of the kinetic energy. For Ro $\lesssim$ 1, the poloidal energy density, $e_{B_\phi}/e_k$, has large fluctuations at POL (especially odd is case RC21) and increases slowly at EQU. Starting from case RC28, Ro $\sim$ 1.2, the poloidal energy shows a clear decay. It is fast at EQU, $e_{B_\phi}/e_k \propto$ Ro$^{-6.4}$ (black dashed lines), and slow at POL, $e_{B_\phi}/e_k \propto$ Ro$^{-2.3}$ (a similar trend is observed in $e_{B_\phi}/e_k$ at POL).

In the NSL, the $\alpha$-effect coefficient, $C_\alpha$, is positive in both POL and EQU latitudes for model RC07 and negative for all other cases. The shear coefficient, $C_\Omega^\theta$, is positive at both POL and EQU, with amplitudes similar to those in the CZ. On the other hand, the radial shear coefficient, $C_\Omega^r$, as in the solar near-surface shear layer, is negative at POL and EQU (except for

case RC07 at EQU). At POL, it reaches its minimum value for case RC35, Ro $\sim$ 1.25, followed by a step decrease for larger values of Ro.

The dynamo coefficients show trends that partially explain the behavior of the magnetic field energy density in this region. Nevertheless, similar to the CZ region, nonlocal processes seem to be relevant for defining the amplitude of the magnetic fields.

In Figure 7, we perform the same analysis as in Figure 5 but for models (a) RC07, (b) RC21, and (c) RC49 and with all quantities averaged in longitude and over the radial extent of the NSL region. Together with Figure 4, this figure illustrates the dynamo behavior for different Rossby numbers at the outer layers. Of particular interest are the following points.

1. In the POL region, the toroidal field energy increases with Ro (Figure 4). This trend does not follow the scaling of $C_\Omega^r$, which, at similar latitudes, has a minimum at Ro $\sim$ 1 and decreases sharply for the largest Ro. We note that the scaling of $e_{B_\phi}/e_k$ with Ro in this region reproduces the behavior of the toroidal energy at the TAC, evidently with less energy. The second and third rows of Figure 7 show that for model RC07 (column (a)), which is an example of simulations with Ro < 1, there are no dynamo sources at high latitudes. Correspondingly, the toroidal field is rather small. For the simulation RC21 (column (b)), where Ro $\sim$ 1, the local shear and $\alpha$-effect are responsible for the generation of a toroidal field. For simulation RC49 (column (c)), representing Ro > 1, the values of both the shear and the $\alpha$-effect are significant; nevertheless, the toroidal field shows poor correlation with these quantities. The movies presented in the supplementary material clearly show that in simulations RC24–63, the toroidal field at the TAC is transported from the TAC to the NSL.

2. In the EQU region, for Ro < 1, the toroidal magnetic energy, $e_{B_\phi}/e_k$, seems to be independent of its values at the TAC. The second and third panels of Figure 7(a) suggest that $C_\alpha$ and $C_\Omega^r$ contribute locally to the generation of the toroidal field. The source terms seem to be out of phase with the magnetic field because the generation occurs slightly below in the CZ (see also panel (a) of Figure 4 and the movie corresponding to simulation RC07). Local toroidal field generation is also observed in simulation RC21, where the field correlates well with the shear term (second panel of Figure 7(b)). In model RC49 (representative of Ro > 1), there is only a marginal correlation between the toroidal field and its local source
<figure><img src="image_16.png" /><figcaption>Figure 7. Same as Figure 5 but for simulations (a) RC07, (b) RC21, and (c) RC49. The average in this case is done over the NSL radial extent.</figcaption></figure>

terms. Thus, the decay of the toroidal energy for $Ro \gtrsim 1.2$ is likely a consequence of the decaying toroidal field at the TAC and the enhanced magnetic diffusivity at lower latitudes.

3. Unlike $\bar{B}_{\phi}$, which is generated by the shear and the $\alpha$-effect, the poloidal magnetic field is solely generated by $\alpha$, specifically by the term $\nabla \times \alpha \bar{B}_{\phi}$. Since the poloidal source term, $C_{\alpha}$, is roughly constant with $Ro$ in the NSL region, one should expect that the poloidal magnetic energy roughly follows the trend of the toroidal energy. This proportionality is observed at the POL region for $Ro \lesssim 1$. However, for $Ro > 1$, while $e_{B_{\phi}}/e_k$ increases, $e_{B_{p}}/e_k$ decays. Relevant hints to understand this behavior can be found by watching the supplementary movies for simulations RC24–RC49. We find that when the toroidal field (on the right side quadrants) quickly rises from the TAC to the NSL, the existent poloidal flux is rapidly redistributed in the bulk of the convection zone, and a new poloidal field of opposite polarity is generated. This is clear evidence of the nonlocal effects present in the simulations. Although it is not easy to make a quantitative analysis, we suggest that the decay of $e_{B_{p}}/e_k$ with $Ro$ is due to the fact that the total magnetic energy is independent of the rotation rate. Thus, the more

toroidal field is rapidly deposited into the upper layers of the domain, the less poloidal field might reside there.

4. At the EQU region, $e_{B_{p}}/e_k$ roughly follows $e_{B_{\phi}}/e_k$ for all values of $Ro$. Note in Figure 7 (fourth row) how the field is locally generated for simulation RC07 (column (a)), has a minimum for simulation RC21 (column (b)), and shows diffusive values in model RC49 (column (c)), in agreement with what is observed in Figure 4. In simulations RC56 and RC63, the toroidal field at the TAC is steady at POL latitudes; therefore, a weak poloidal field develops at EQU.

### 3.3. Comparison with the Observations

The top boundary of our model is placed at $r_t = 0.96$ of the stellar radius. However, this does not preclude the relevance of simulated NSL properties to observations. Unfortunately, it is not yet clear how magnetic fields erupt to the surface to form star spots and how this emergence process depends on fluid properties such as rotation and convective motions. For instance, it is not clear what is the correspondence between the $\langle R'_{HK} \rangle$ flux and the magnetic field in stellar interiors. Also, there is no complete interpretation of the magnetic fields
<figure><img src="image_17.png" /><figcaption>**Figure 8.** (a) Scaling of the NSL magnetic field components with Ro; the conventions are the same as in Figure 4. (b) Correlation between the toroidal and poloidal magnetic energies at the NSL. The light blue (green) dots correspond to the EQU (POL) region. The dashed lines show the power law.</figcaption></figure>

inferred by the ZDI technique, and a few shortcomings of this method have been recently identified (Lehmann et al. 2019).

To establish some connection between the results presented above and the observations, we remind the reader that the results of Vidotto et al. (2014) correspond to solar-like stars with Rossby numbers spanning from 0.3 to 3. All of the stars in this sample follow the relation $\bar{B} \propto \text{Ro}^{-1.38}$. In addition, See et al. (2015) reported a power-law relation between the toroidal and poloidal magnetic energies with $\bar{B}_{\phi}^2 \propto (\bar{B}_p^2)^{1.25}$ for stars with masses between 0.5 and 1.5 $M_{\oplus}$ in the same range of Ro.

As depicted in Figure 8(a), our simulations correspond to $0.36 < \text{Ro} < 2.12$ ($0.26 < \text{Ro} < 3.03$ according to the definition of Noyes et al. 1984b). Within this interval, we found different behaviors at the EQU and POL regions and, unlike the observations, two different scalings with Ro. In the EQU region, for $\text{Ro} \lesssim 1.2$, both $\bar{B}_{\phi}$ and $\bar{B}_p$ seem independent of Ro. For $\text{Ro} \gtrsim 1.2$, our simulations predict $\bar{B}_{\phi}^{\text{NSL}} \propto \text{Ro}^{-1.2}$ and $\bar{B}_p^{\text{NSL}} \propto \text{Ro}^{-2.9}$. In the POL region, the toroidal field increases with the Rossby number as $\bar{B}_{\phi}^{\text{NSL}} \propto \text{Ro}^{1.1}$. Yet the poloidal field increases similarly to $\bar{B}_{\phi}$ for $\text{Ro} \lesssim 1.2$ and then decreases as $\bar{B}_p^{\text{NSL}} \propto \text{Ro}^{-0.9}$ for large values of Ro.

As discussed in the previous section, there is no straightforward interpretation of these scaling relations; however, we can summarize our findings as follows. For $\text{Ro} \lesssim 1$, there is local dynamo action that occurs mostly at the EQU region. On the other hand, for $\text{Ro} \gtrsim 1$, there is a rapid transport of $\bar{B}_{\phi}$ from the bottom to the top of the convection zone. These effects, however, are more pronounced and effectively change the scaling relations for $\text{Ro} \gtrsim 1.2$. Thus, the toroidal field at surface levels scales with Ro in the same way as it does at the TAC; i.e., it increases and is located closer and closer to the poles with increasing Ro. Consequently, the toroidal field at the equator diminishes as Ro increases. Furthermore, large fractions of poloidal magnetic flux are quickly removed from the NSL and redistributed in the convection zone. We believe that the faster decay of the poloidal field at EQU compared to POL is due to the latitudinal variation of the turbulent magnetic diffusivity. Also, it is worth noticing that the contribution of other transport mechanisms like meridional circulation or turbulent pumping cannot be ruled out.

The correlation between the toroidal and poloidal magnetic field energies is presented in Figure 8(b). The light blue and

green points correspond to the EQU and POL regions, respectively. It can be seen that the poloidal field energy is an increasing power-law function of the toroidal energy with coefficients 0.7 at EQU and 1.3 at POL. However, the correlation is better defined at the equator than at the poles, confirming that at lower latitudes, the field strength depends mostly on local dynamo contributions. At the pole, the relation shows more dispersion, especially for stronger magnetic fields. This is understandable if nonlocal sources are contributing to the toroidal field but not to the poloidal one, which, moreover, is expelled from the places where the toroidal flux increases.

### 3.4. The Dynamo Period

The magnetic cycles in Figure 2 are clear and well defined for most of the models. We compute the period by using a Fourier transform of the data. The technical details are presented in Appendix C. In Figure 9, we show the (a) $P_{\text{cyc}}$ versus $P_{\text{rot}}$ and (b) $\log(P_{\text{rot}}/P_{\text{cyc}})$ versus $\log(1/\text{Ro})$ representations of the dynamo cycle period for simulations RC07–RC49 (black stars). In both panels, we have plotted the observational data as reported by Brandenburg et al. (2017). The blue and red stars correspond to the active and inactive branches, respectively. The 11 yr cycle of the Sun is represented with a yellow star. In panel (a), we notice that, similar to the observations, the magnetic cycle period increases with the period of rotation. If simulation RC21 (i.e., the third black star from left to right in panel (a)) is discarded, the trend appears linear and nearly separates the active (blue) from the inactive (red) branch.

In Figure 9(b), all values of $P_{\text{rot}}/P_{\text{cyc}}$ fall closer to the active branch. However, two different trends can be identified: one with a positive inclination for cases RC07–RC21 (black stars for $\log(1/\text{Ro}) \gtrsim 0$) and one with a negative inclination for cases RC21–RC49 ($\log(1/\text{Ro}) \lesssim -0.1$). Simulation RC18, which shows no signals of periodicity in the convection zone, indicates that simulation RC21 actually is close to the transition between bifurcated regimes, one with the oscillatory dynamo operating mainly in the convection zone and the other with the tachocline having a predominant role. We recall that in model RC21, a relevant oscillatory toroidal field developed at the tachocline coexists with a steady poloidal field in the stable layer. Thus, model RC21 effectively belongs to the branch with negative inclination, together with models RC24–RC49.
<figure><img src="image_18.png" /><figcaption>**Figure 9.** (a) $P_{\text{cyc}}$ vs. $P_{\text{rot}}$ and (b) $\log(P_{\text{rot}}/P_{\text{cyc}})$ vs. $\log(1/R_o)$ representations of the magnetic cycle period against the rotation period. The simulation results are shown with black stars. The blue and red stars correspond, respectively, to active and inactive branches of activity, with data taken from Brandenburg et al. (2017). Stars correspond to F and G stars, and asterisks correspond to K-type stars. The yellow stars show the 11 yr solar activity cycle. In the left panel, the black dashed line compares the rotation period with $B_{\phi}^{\text{TAC}}$. These quantities are normalized to the cycle period of model RC28.</figcaption></figure>

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

In Guerrero et al. (2016a), we discussed the instabilities that may occur in the tachoclines and radiative zones. We showed that the so-called magneto-shear instability is the most likely source of the magnetic α-effect. This instability belongs to the Tayler instability family modified by the presence of shear. The Tayler instability is related to the decay of a large-scale toroidal field in a stable stratified layer (Tayler 1973). In the nonlinear phase, it results in a saturated state with nonzero helicity (Bonanno & Urpin 2012), which, in turn, might develop large-scale magnetic fields. The growth rate of this instability is inversely related to the ratio $\Omega_0/B_\phi$ (Bonanno & Urpin 2013). This means that it is inhibited by fast rotation or enhanced by strong toroidal magnetic fields. On the other hand, the shear contributes to this complex process by replenishing the toroidal field (Miesch 2007; Szklarski & Arlt 2013). However, its relevance still needs to be quantified.

The analysis of dynamo results presented in Lawson et al. (2015) and Guerrero et al. (2016a) suggests that there is an exchange of energy between the development of a large-scale toroidal field and turbulent motions and currents, i.e., the interaction between the dynamo and magneto-shear. The timescale of this exchange can set the activity cycle period, at least in dynamos where reversals of both magnetic field components take place in the stable layer (RC24–RC49). This process also allows for the existence of the bimodality observed in RC21 as reported by Rogers (2011) for axisymmetric dynamo simulations.

The sequence of images in Figure 10 presents in more detail the reversal process for simulation RC49. The figure shows snapshots of the magnetic field lines around $r = 0.66R_\odot$ during a half-cycle period; the colors represent the direction of the toroidal field, and the thickness of the lines represents the magnitude of the magnetic field. For the reader's convenience, the top panels of the figure repeat the butterfly diagrams of Figure 5 presenting the evolution of $\bar{B}_\phi$, $\bar{B}_r$, and $\alpha$ in the TAC region. The vertical dashed lines indicate the times of the corresponding snapshots.

The leftmost vertical dashed lines in the top panels of Figure 10 correspond to a time when a negative toroidal field covers almost the entire northern hemisphere. At this moment, the radial field is negative and in a growing phase. Panel (a) of

Figure 10 corresponds to this stage. A few years later, while a new positive toroidal field is originated at POL latitudes, a negative $\alpha$ develops at the same locations. This can be seen in panels (a) and (b) of Figure 10, where the magnetic field, mostly toroidal and negative, has a positive inclination near the pole with respect to the equator. This tilt, characteristic of the so-called clamshell instability (Cally 2003), generates poloidal field components that are wound up by the differential rotation at the poles first and then at the equator, where the shear is stronger; see panels (a)–(d) of Figure 10. This generates a new positive toroidal field that migrates poleward. In panel (d), the radial field reaches its maximum, and the remaining negative toroidal field at POL latitudes continues to be unstable. However, this time it shows a negative tilt with respect to the horizontal direction, giving rise to positive values of $\alpha$; see the top right panel of Figure 10. This $\alpha$-effect seems to be responsible for the generation of both a positive radial field and a positive toroidal field at intermediate latitudes. Note in panels (c)–(f) that although the field lines are erratic, some of them are oriented in the latitudinal direction first and the azimuthal direction later. The positive toroidal field finally covers almost the entire hemisphere. A new positive tilt is observed near the poles, together with the remnants of the negative field. This is the initial configuration for the second part of the full cycle.

We have described the reversal process for model RC49. For other oscillatory simulations with fast rotation, the process is fundamentally the same (see the movie of the field reversal for simulation RC28 in the supplementary material). However, according to Figure 4, increasing the rotation leads to larger poloidal fields and smaller toroidal fields. Thus, for a field configuration analogous to that of Figure 10(a), the toroidal field is less coherent and decays faster, giving rise to well-organized field lines in the latitudinal direction. These, in turn, are rapidly stretched by the equatorial shear, resulting in a faster cycle.

In Figure 9(a), the black dotted line depicts the amplitude of the toroidal field as a function of the rotational period. Its trend agrees with that of the cycle period, suggesting proportionality between the two quantities. As discussed in the third paragraph of this subsection, this is at odds with the linear theory of the Tayler instability, which predicts fast
<figure><img src="image_20.png" /><figcaption>**Figure 10.** Snapshots of the magnetic field lines around $r = 0.66R_⊙$ covering one polarity reversal (half-cycle period) of model RC49. The blue (red) colors correspond to toroidal fields pointing eastward (westward), and the thickness of the lines is proportional to the magnitude of the magnetic field. The top panels of the figure repeat the butterfly diagrams of Figure 5 presenting the evolution of $\bar{B}_\phi$, $\bar{B}_r$, and $\alpha$ (from left to right) in the TAC region. The dashed vertical lines show the times of the corresponding snapshots. An animation of the field line evolution of simulations RC28 and RC49 is available. The video duration is 7 s.</figcaption></figure>

(An animation of this figure is available.)

development of the unstable modes for the strongest fields. Because rotation also stabilizes toroidal fields, the fact that the toroidal field decays faster for simulations with rapid rotation is also against the linear theory. Nevertheless, our results correspond to a steady nonlinear state of the simulations, which is hardly comparable to the linear phase of the instability. Miesch et al. (2007) studied the development of shear-current instabilities in the tachocline through nonlinear MHD simulations. Even though the unstable modes reported by them have similarities to those presented in Figure 10, in their simulations, a latitudinal shear is imposed as an initial condition. In our case, the shear is mostly radial, and it is continuously replenished by the rotation and the convective motions in the convective zone. Therefore, the mechanisms driving the instabilities are different, and the results cannot be compared. The stability of the magnetic fields in stable stratified atmospheres under these circumstances has not yet been explored in detail. A dynamo operating in these radiative zones has been envisaged by, e.g., Spruit (2002), Zahn et al. (2007), and Bonanno (2013); nevertheless, the results presented here (and perhaps also the simulations of Ghizaru et al. 2010, and subsequent papers) are the first to capture the entire process from first principles. Idealized simulations, where these processes are studied separately, are still necessary to clarify the contribution of each one (G. Guerrero et al. 2019, in preparation).

# 4. Summary and Conclusions

We analyzed the results of global dynamo simulations in models that have a solar-like stratification and span a wide range of rotation rates. In all of the simulated cases, the resulting mean flows exhibit a fast equator and slower poles, and the formation of radial shear layers in two locations, the interface between a convectively stable layer at the bottom of the domain and the convection zone (the tachocline) and the uppermost layer of the domain (the near-surface shear layer). The magnetic fields evolve in a variety of dynamo modes, from oscillatory dynamos with short periods occurring mostly in the convection zone to dynamos mostly happening at the tachocline with periodic, steady, and mixed modes. The MHD properties of some of the simulations have been studied in detail in previous works (Guerrero et al. 2016a, 2016b). In this work, we focused our analysis on the magnetic activity, i.e., the magnetic field strength and cycle period, and compared our results with recent observational findings.

The butterfly diagrams of the simulated stellar dynamos demonstrate the complexity of the systems (Figure 2). One of the characteristics is the existence of dynamo sources in different parts of the domain that add features to the spatiotemporal evolution. For instance, oscillatory dynamos for rotational periods between 24 and 49 days have their magnetic timescale governed by the tachocline dynamics;
however, due to the dynamo action in the near-surface layer, new branches of activity appear, modifying the solution.

In view of these intricacies, we have performed an analysis considering three different shells within the domain (TAC, covering the shear region and the stable stratified layer; CZ, enclosing the bulk of the convection zone; and NSL, covering the uppermost fraction of the model) and two latitudinal zones (EQU, from the equator to 45°; and POL, from 45° to the north pole). The volume rms values of the normalized magnetic field energy density; the nondimensional dynamo coefficients $C_\alpha$, $C_\Omega$, and $C_u$; and the dynamo numbers $D'_r$ and $D'_\theta$ were computed for these regions. These quantities appear in Figure 4 as a function of the Rossby number. This figure summarizes and quantifies what can also be observed in Figures 1, 2, and 13. Our most relevant findings are as follows.

1. The total magnetic energy density is independent of the rotational period, yet the energy density in the toroidal and poloidal components of the magnetic field is determined by Ro and reaches different amplitudes at different depths. Although the dynamo numbers in the three regions have similar values, the magnetic energy in the TAC is 1 order of magnitude larger than in the CZ and NSL. This is a consequence of the stable stratified layer, which allows longer storage and further amplification of the magnetic fields.

2. In the TAC region, for $Ro \lesssim 1$, the poloidal field is larger than the toroidal one. They have similar energies at $Ro \sim 1$, and for $Ro > 1$, the toroidal field energy dominates. This is an outcome of the scaling of $C_\alpha$ (at higher latitudes) and $C_\Omega^r$ (in the entire hemisphere) with the Rossby number. For $Ro \lesssim 1$, $C_\alpha$ has larger amplitude and negative values, while $C_\Omega^r$ is smaller. For $Ro \sim 1$, $C_\alpha$ changes sign, while $C_\Omega^r$ reaches larger values that enhance the generation of the toroidal field.

3. At $Ro \sim 1$, when $C_\alpha$ becomes positive and the radial shear at the TAC is sufficiently strong to generate large toroidal fields, the dynamo simulations start to develop deep-seated cyclic dynamos (models RC21–RC49). For $Ro \gtrsim 1.7$ (models RC56–RC63), $C_\alpha$ and $C_\Omega^r$ are still significant; however, Figures 1(h)–(i) and 13(h) reveal that these coefficients are prominent in latitudes closer to the poles. The magnetic field in these slow-rotating cases is steady and concentrated at higher latitudes (Figure 2(h)).

4. Case RC21 exhibits bimodality at the TAC; i.e., a steady radial field is mixed with an oscillatory toroidal field. This transitional case is confirmed by simulations RC18, with a steady dynamo in the TAC and no cyclic behavior in the CZ, and RC24, where an oscillatory dynamo is observed in the TAC.

5. In the NSL region, for $Ro \gtrsim 1.2$, the normalized magnetic energy densities, $e_{B_\phi}$ and $e_{B_\rho}$, decay at EQU latitudes as $Ro^{-2.8}$ and $Ro^{-6.4}$, respectively. At higher latitudes, $e_{B_\phi}$ increases as $Ro^{1.8}$ and $e_{B_\rho}$ decays as $Ro^{-2.3}$; see the scaling laws for the magnetic field components in Figure 8. In Section 3.2.3 we have argued that finding an explanation for these scaling laws is not straightforward. In some parts of the domain, the local generation is relevant; in others, there are nonlocal sources of magnetic flux that cannot be easily quantified. Furthermore, we demonstrated that the turbulent magnetic diffusivity varies in the latitudinal direction, differently affecting

the decay of the field in the POL and EQU regions. It is important to bear in mind that we have computed the dynamo turbulent coefficients by considering the FOSA approximation. The results of Warnecke et al. (2018) indicate that the turbulent diffusivity may even be anisotropic; i.e., it may have different values for the different components of the field.

6. The scaling of the magnetic field energy with Ro and the relation between the toroidal and poloidal field energy densities in the most external layers of the domain are reminiscent of the observational results (Pizzolato et al. 2003; Wright et al. 2011; Vidotto et al. 2014; See et al. 2015; Wright & Drake 2016). Because the magnetic field may have different local and nonlocal sources for different regimes and regions within the domain, these observations cannot be explained by a simple model. Before attaining a satisfactory explanation, it is necessary to better understand the correlation between the different observational signatures and the magnetic field in the interior of the stars.

In Figure 9, we show two different representations, often used in the literature, correlating the rotation period, $P_{\text{rot}}$, with the magnetic cycle period, $P_{\text{cyc}}$. The simulation results clearly show that $P_{\text{cyc}}$ increases with $P_{\text{rot}}$, in agreement with the observations of stars of types G–K (Böhm-Vitense 2007; Brandenburg et al. 2017). Most of the cycle periods obtained in this work are consistent with the A branch described in the literature. However, the bimodal dynamo observed in case RC21 appears out of the curve. We conclude that this case is in the transition between cyclic dynamos operating in the convection zone and cyclic dynamos in the tachocline. As mentioned above, this transition is characterized by the enhanced generation of a deep-seated toroidal field due to an EQU radial shear that increases monotonically with the Rossby number. The nonlinear effects resulting from this strong field (1 order of magnitude larger than in the rest of the domain) are substantial; they can even modify the character of the convective motions (e.g., Passos et al. 2017; Beaudoin et al. 2018) and/or give rise to other instabilities different from the dynamo (Lawson et al. 2015; Guerrero et al. 2016a).

We suggest that the resulting cycle period may be explained by the energy exchange between the dynamo and shear-current instabilities occurring below the tachocline. In the dynamical phase of evolution of the velocity and magnetic fields, both instabilities reach an equilibrium state of energy exchange that behaves like an $\alpha^2\Omega$ dynamo and determines the cycle period. A similar oscillatory exchange of energy in a shear layer between the magnetic field and the nonaxisymmetric kinetic energy was reported by Miesch (2007). However, it is still necessary to quantify how the period of these oscillations depends on the strength of the toroidal field and the amount of shear. Since the entire process is hard to disentangle in convection simulations, a step-by-step analysis of these instabilities is left to a forthcoming work (G. Guerrero et al. 2019, in preparation).

The resulting magnetic fields are a direct consequence of the self-consistent development of tachoclines in our simulations. This characteristic increase of $P_{\text{cyc}}$ with $P_{\text{rot}}$ is precisely what distinguishes them from other convective dynamo simulations, where $P_{\text{cyc}}$ decreases with the increase of $P_{\text{rot}}$ (Strugarek et al. 2017; Warnecke 2018). It is worth mentioning, however, that the observations are still inconclusive in this regard.
Nevertheless, the estimation of stellar magnetic fields is currently an exiting and quite active field of research. Luckily, in the near future, observations will provide a clear picture of the relation between stellar cycle and rotational period. From the numerical point of view, current dynamo models are able to reproduce both scenarios.

We thank an anonymous referee for insightful comments that helped to improve the paper. This work was partly funded by FAPEMIG grant APQ-01168/14 (G.G.), FAPESP grant 2013/ 10559-5 (E.M.G.D.P.), CNPq grant 306598/2009-4 (E.M.G.D. P.), and NASA grants NNX09AJ85G, NNX17AE76A, and NNX14AB70G. The simulations were performed in the NASA cluster Pleiades and the Brazilian supercomputer SDumont of the National Laboratory of Scientific Computation (LNCC). This work was partly funded by FAPEMIG grant APQ-01168/14 and supported by Universidade Federal de Minas Gerais (G.G.).

## Appendix A

### The Convective Turnover Time and the Rossby Number

One of the bottlenecks in the interpretation of stellar activity is the determination of the convective turnover time, $\tau_c$. It is useful for computing the Rossby number, as well as the turbulent dynamo coefficients. Observationally, there is a correlation between $\tau_c$ and the chromospheric emission $\langle R'_{HK} \rangle$ (Noyes et al. 1984a). In simulations, this quantity can be computed through the energy spectrum. Since we are interested in the scales associated with the turbulent quantities, the spectrum is computed from the transformation of, for example, the nonaxisymmetric velocity, $\mathbf{u}' = \mathbf{u} - \bar{\mathbf{u}}$, into spherical harmonics, $Y_l^m(\theta, \phi) = P_l^m(\cos\theta)\exp(im\phi)$, by

$$ \mathbf{u}'(l, m) = \sum_{l,m} Q_l^m Y_l^m(\theta, \phi) + S_l^m r \nabla Y_l^m(\theta, \phi) - T_l^m \mathbf{r} \times \nabla Y_l^m(\theta, \phi), \quad (7) $$

<figure><img src="image_21.png" /><figcaption>Figure 11. Energy spectra of the kinetic (solid lines) and magnetic (dashed lines) energies computed for $r = 0.7R_\odot$ (left) and $r = 0.9R_\odot$ (right) for some representative models. The black dashed line shows the Kolmogorov $k^{-5/3}$ energy scaling.</figcaption></figure>

where the expansion coefficients, $Q_l^m$, $S_l^m$, and $T_l^m$, are computed with the optimized library SHTNs (Schaeffer 2013).

The kinetic energy spectra is computed by

$$ \tilde{E}_k = 4\pi \sum_{m=-l}^{l} [ |Q_l^m|^2 + l(l+1)(|S_l^m|^2 + |T_l^m|^2) ]. \quad (8) $$

With a similar decomposition, we obtain the spectral energy of the nonaxisymmetric magnetic field, $\tilde{E}_m$. In Figure 11, we present the kinetic (solid lines) and magnetic (dashed lines) energy spectra for the radial levels $r = 0.7R_\odot$ (left) and $r = 0.9R_\odot$ (right). Note that in the convection zone, the kinetic energy dominates over the magnetic energy. For the faster-rotating simulations, the spectra peak at larger values of $k$. With the decrease of the rotation (large period), it moves toward smaller values of $k$ and reaches large values for the energy. This behavior is expected from simulations dominated by rotation toward simulations dominated by convection. In the radiative zone, the magnetic energy is dominant, and the spectra peak at small wavenumbers. Since we are interested in the scales that carry most of the energy as a function of depth, we compute a similar spectrum for each radial level of the simulation and compute such length scales as

$$ \ell_{k,m}(r) = \frac{r \int_k \frac{\tilde{E}_{k,m}(k,r)}{k} dk}{\int_k \tilde{E}_{k,m}(k,r) dk}, \quad (9) $$

where $k = l + 1/2$, according to the Jeans rule. The convective turnover time is computed by $\tau_c(r) = \ell_k(r)/u_{rms}(r)$. The radial profiles of $\ell$, $\tau_c$, and $u_{rms}$ are depicted in Figures 12(a)–(c), respectively. Finally, the Rossby number is computed as $Ro = P_{\text{rot}}/(2\pi\tau_c^{\text{NSL}})$ (see Table 1).

<figure><img src="image_22.png" /></figure>
<figure><img src="image_23.png" /><figcaption>Figure 12. (a) Length scale of the energy-carrying eddies, $\ell(r)$; (b) convective turnover times, $\tau_c$; (c) radial profile of the $u_{rms}$ velocity; and (d) turbulent magnetic diffusivity. Different colors are assigned to some representative simulations between RC07 and RC63 (see annotations in panel (d)). The dashed lines in panel (b) depict the Alfvén timescale, $\tau_A$. See the text for definitions.</figcaption></figure>

## Appendix B: Mean-field Turbulent Coefficients

We use the FOSA to compute the dynamo turbulent coefficients. Under this formalism, the $\alpha$-effect is given by

$$ \alpha = \alpha_k + \alpha_m = -\frac{\tau_c}{3} \langle \omega' \cdot u' \rangle + \frac{\tau_c}{3} \langle j' \cdot B' \rangle / \rho_e, \quad (10) $$

where $\omega' = \nabla \times u'$ and $j' = \nabla \times B'$ are the small-scale vorticity and current, respectively, and $\tau_c$ is the convective turnover time defined in Appendix A. Note that the convective motions spread along the unstable region; therefore, $\tau_c$ in Equation (10) is valid only down to $r \sim 0.74R_\odot$. This is not relevant for $\alpha_k$, since the velocity drops to small values in the radiative zone. However, it presents a problem for $\alpha_m$, as a significant amount of current helicity, $\langle j' \cdot b' \rangle$, develops below this radius.

As can be noticed in the left panel of Figure 11, at $r = 0.7R_\odot$, the nonaxisymmetric magnetic energy dominates over the kinetic energy. We believe that the magnetic field develops at and below the tachocline because of shear-current instabilities. Therefore, we associate the timescale of the magnetic $\alpha$-effect below the tachocline with one of the timescales associated with these instabilities, namely, the Alfvén time, $\tau_A = \ell_m/v_A$ (dashed lines in the top right panel of Figure 12), where $v_A = B'_{rms}/(\mu_0\rho_e)^{1/2}$. Here $B'_{rms} = (\langle B_r'^2 + B_\theta'^2 + B_\phi'^2 \rangle_{\phi,\theta,t})^{1/2}$ is the nonaxisymmetric magnetic field averaged over $\phi, \theta$, and time.

The profiles of the kinetic (left), magnetic (middle), and total (right) $\alpha$-effect for simulations RC07–RC63 are presented in panels (a)–(h) of Figure 13. A quantitative analysis of these profiles is presented in Section 3.2. Finally, the turbulent diffusivity coefficient presented in Figure 12(d) is evaluated as

$$ \eta_t = \frac{1}{3} \tau_c u'^2. \quad (11) $$
<figure><img src="image_24.png" /><figcaption>**Figure 13.** Kinetic, magnetic, and total α-effect, from left to right, for representative models between RC07 and RC63, in panels (a)–(h).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 14.** Spectral density computed from the volume-averaged time series of the poloidal field ($B_p$). The dotted lines show the frequency where the spectral density peaks. For models RC24–RC49, the period was computed in the regions NSL (solid lines) and TAC (dashed lines). For clear visualization, the spectral density in the region TAC was multiplied by a factor of $10^2$.</figcaption></figure>

## Appendix C
### Cycle Period Evaluation

The magnetic cycle period, $P_{\text{cyc}}$, is evaluated through the Fourier transformation of the rms vertical magnetic field (see Figure 14). Since the oscillatory behavior observed in models RC07–RC21 (left panel) is prominent close to the equator, to obtain a better estimation of the period, we consider a latitude range of $\pm 10^\circ$ latitude. In radius, these models are oscillatory only in the CZ; therefore, the radial average is made in the NSL. While for models RC14 and RC21, the frequency with maximal spectral density is clearly defined, model RC07 seems to have multiple periodicities. We have chosen the peak in the spectral density that matches better with the periodicity observed in the butterfly diagram ($P_{\text{cyc}} = 3.0$ yr). For models RC24–RC49 (right panel of Figure 14), we consider latitudes between 0° and 90°. We compute the period in the shells NSL (solid lines) and TAC (dashed lines). For all of these models, the periods agree in both regions; they are presented in Table 1.

## ORCID iDs

G. Guerrero https://orcid.org/0000-0002-2671-8796

B. Zaire https://orcid.org/0000-0002-9328-9530

P. K. Smolarkiewicz https://orcid.org/0000-0001-7077-3285

E. M. de Gouveia Dal Pino https://orcid.org/0000-0001-8058-4752

A. G. Kosovichev https://orcid.org/0000-0003-0364-4883

N. N. Mansour https://orcid.org/0000-0002-3927-3917

## References

Baliunas, S. L., Donahue, R. A., Soon, W. H., et al. 1995, ApJ, 438, 269

Beaudoin, P., Strugarek, A., & Charbonneau, P. 2018, ApJ, 859, 61

Blackman, E. G., & Thomas, J. H. 2015, MNRAS, 446, L51

Böhm-Vitense, E. 2007, ApJ, 657, 486

Bonanno, A. 2013, SoPh, 287, 185

Bonanno, A., & Urpin, V. 2012, ApJ, 747, 137

Bonanno, A., & Urpin, V. 2013, MNRAS, 431, 3663

Brandenburg, A., Mathur, S., & Metcalfe, T. S. 2017, ApJ, 845, 79

Brandenburg, A., Rädler, K.-H., Rheinhardt, M., & Käpylä, P. J. 2008, ApJ, 676, 740

Brandenburg, A., Saar, S. H., & Turpin, C. R. 1998, ApJL, 498, L51

Cally, P. S. 2003, MNRAS, 339, 957

Cossette, J.-F., Charbonneau, P., Smolarkiewicz, P. K., & Rast, M. P. 2017, ApJ, 841, 65

do Nascimento, J.-D., Jr., Vidotto, A. A., Petit, P., et al. 2016, ApJL, 820, L15

Domaradzki, J. A., Xiao, Z., & Smolarkiewicz, P. K. 2003, PhFl, 15, 3890

Elliott, J. R., & Smolarkiewicz, P. K. 2002, IJNMF, 39, 855

Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJL, 715, L133

Gregory, S. G., Donati, J.-F., Morin, J., et al. 2012, ApJ, 755, 97

Guerrero, G., & de Gouveia Dal Pino, E. M. 2008, A&A, 485, 267

Guerrero, G., Smolarkiewicz, P. K., de Gouveia Dal Pino, E. M., Kosovichev, A. G., & Mansour, N. N. 2016a, ApJ, 819, 104

Guerrero, G., Smolarkiewicz, P. K., de Gouveia Dal Pino, E. M., Kosovichev, A. G., & Mansour, N. N. 2016b, ApJL, 828, L3

Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176

Jouve, L., Brown, B. P., & Brun, A. S. 2010, A&A, 509, A32

Karak, B. B., Kitchatinov, L. L., & Choudhuri, A. R. 2014, ApJ, 791, 59

Lawson, N., Strugarek, A., & Charbonneau, P. 2015, ApJ, 813, 95

Lehmann, L. T., Hussain, G. A. J., Jardine, M. M., Mackay, D. H., & Vidotto, A. A. 2019, MNRAS, 483, 5246

Lehtinen, J., Jetsu, L., Hackman, T., Kajatkari, P., & Henry, G. W. 2016, A&A, 588, A38

Margolin, L. G., Smolarkiewicz, P. K., & Wyszogrodzki, A. A. 2006, JAM, 73, 469

Miesch, M. S. 2007, ApJL, 658, L131

Miesch, M. S., Gilman, P. A., & Dikpati, M. 2007, ApJS, 168, 337

Moffatt, H. K. 1978, Magnetic Field Generation in Electrically Conducting Fluids (Cambridge: Cambridge Univ. Press)

Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984a, ApJ, 279, 763

Noyes, R. W., Weiss, N. O., & Vaughan, A. H. 1984b, ApJ, 287, 769

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

Passos, D., Miesch, M., Guerrero, G., & Charbonneau, P. 2017, A&A, 607, A120

Petit, P., Dintrans, B., Solanki, S. K., et al. 2008, MNRAS, 388, 80

Pipin, V. V., & Kosovichev, A. G. 2016, ApJ, 823, 133

Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147

Prusa, J. M., Smolarkiewicz, P. K., & Wyszogrodzki, A. A. 2008, CF, 37, 1193

Racine, É., Charbonneau, P., Ghizaru, M., Bouchat, A., & Smolarkiewicz, P. K. 2011, ApJ, 735, 46

Rogers, T. M. 2011, ApJ, 735, 100

Saar, S. H., & Brandenburg, A. 1999, ApJ, 524, 295

Schaeffer, N. 2013, GGG, 14, 751

See, V., Jardine, M., Vidotto, A. A., et al. 2015, MNRAS, 453, 4301

Smolarkiewicz, P. K. 2006, IJNMF, 50, 1123

Smolarkiewicz, P. K., & Charbonneau, P. 2013, JCoPh, 236, 608
Smolarkiewicz, P. K., & Margolin, L. G. 2007, in Studies in Geophysics, ed. F. F. Grinstein & L. G. Margolin (Cambridge: Cambridge Univ. Press), 413

Spruit, H. C. 2002, A&A, 381, 923

Steenbeck, M., Krause, F., & Rädler, K.-H. 1966, ZNatA, 21, 369

Stix, M. 1976, in IAU Symp. 71, Basic Mechanisms of Solar Activity, ed. V. Bumba & J. Kleczek (Dordrecht: Reidel Publishing Co.), 367

Strugarek, A., Beaudoin, P., Charbonneau, P., Brun, A. S., & do Nascimento, J.-D. 2017, Sci, 357, 185

Szklarski, J., & Arlt, R. 2013, A&A, 550, A94

Tayler, R. J. 1973, MNRAS, 161, 365

Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014, MNRAS, 441, 2361

Viviani, M., Warnecke, J., Käpylä, M. J., et al. 2018, A&A, 616, A160

Warnecke, J. 2018, A&A, 616, A72

Warnecke, J., Rheinhardt, M., Tuomisto, S., et al. 2018, A&A, 609, A51

Waruszewski, M., Kühnlein, C., Pawlowska, H., & Smolarkiewicz, P. K. 2018, JCoPh, 359, 361

Wright, N. J., & Drake, J. J. 2016, Natur, 535, 526

Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48

Zahn, J.-P., Brun, A. S., & Mathis, S. 2007, A&A, 474, 145