© ESO 2018

1. Introduction

The Sun, our nearest late-type star exhibits a magnetic activity cycle with a period of around 11 yr. The cyclic magnetic field is generated by a dynamo operating below the surface, where it converts the energy of rotating convective turbulence into magnetic energy. We are still far from fully understanding the solar dynamo mechanism (e.g., Ossendrijver 2003; Charbonneau 2014). One reason is the limited information about the dynamics in the solar convection zone. Helioseismology has provided us with the profile of temperature and density stratification and the differential rotation (e.g., Schou et al. 1998) in the interior. Further information, such as the meridional circulation profiles, convective velocity strength, or even magnetic field distributions, are currently inconclusive or not even possible (e.g., Basu 2016; Hanasoge et al. 2016). One way to investigate how important differential rotation, meridional circulation, and turbulent convective velocities are for the solar dynamo is to use numerical simulations. Since the early simulations by Gilman (1983), the increase of computing resources has prompted several advances using numerical simulations. Nowadays, global simulations of convective dynamos are able to reproduce cyclic magnetic fields and dynamo solutions resembling many features of the solar magnetic field evolution (Ghizaru et al. 2010; Käpylä et al. 2012; Warnecke et al. 2014; Augustson et al. 2015), even long-time evolution (Augustson et al. 2015; Käpylä et al. 2016; Beaudoin et al. 2016). The cyclic magnetic field in these simulations can be well understood in terms of Parker–Yoshimura rule (Parker 1955; Yoshimura 1975; Warnecke et al. 2014), in which a propagating αΩ dynamo wave is excited; see also Gastine et al. (2012). The α effect (Steenbeck et al. 1966) describes the magnetic field enhancement from helical turbulence and the Ω effect the shearing of magnetic field caused by differential rotation. The propagation direction of the dynamo wave depends on the sign of α and shear: to generate an equatorward propagating wave, the product of α and the radial gradient of Ω must be negative (positive) in the northern (southern) hemisphere. Explaining the solar equatorward propagation of the sunspot appearance by the Parker–Yoshimura rule therefore requires either invoking the near-surface-shear layer (Brandenburg 2005), because the radial gradient is only negative in that layer (Barekat et al. 2014) and α is positive, or changing the sign of α in the bulk of the convection zone (Duarte et al. 2016) where the radial shear is positive. Furthermore, to understand the magnetic field evolution fully in the global numerical simulation one needs suitable analysis tools to extract the important contribution of turbulent dynamo effects. One of these tools is the test-field method (Schrinner et al. 2005, 2007; Warnecke et al. 2018). This method allows the determination of the turbulent transport coefficients directly from the simulations. This includes the measurement of tensorial coefficients such as α, turbulent pumping, and turbulent diffusion. The first application to global convection simulations of solar-like dynamo has already revealed that the turbulent effects can have a significant impact on large-scale magnetic field dynamics (Warnecke et al. 2018; Gent et al. 2017).

Another approach for understanding the solar dynamo makes use of the observation of other stars. Since the Mount Wilson survey, we know that many stars exhibit cyclic magnetic activity (e.g., Noyes et al. 1984a,b; Baliunas et al. 1995). In this survey, the authors observe solar-like stars in the chromospheric Ca II H&K band, which is used as a proxy for magnetic activity. Using this data Brandenburg et al. (1998) and Saar & Brandenburg (1999) found two distinct branches, when they plotted the ratio of rotational period and activity cycle period P_rot/P_cycl over the rotational influence on the stellar convection in terms of the inverse Rossby number. The two branches are called inactive and active branch because of their preferred magnetic activity, divided by the so-called Vaughan–Preston gap (Vaughan & Preston 1980). Their slopes are positive in terms of rotational influence, which means the cycle period decreases faster than linear with increasing rotation. This agrees qualitatively with the findings of Noyes et al. (1984b), who obtained $P_{\text{cycl}}\propto P_{\text{rot}}^{1.25}$. However, in their recent reanalysis of the Mount Wilson data, Oláh et al. (2016) found a relation of $P_{\text{cycl}}\propto P_{\text{rot}}^{0.24}$. The activity branches of Brandenburg et al. (1998) have been recently supported (Brandenburg et al. 2017), but also questioned (Reinhold et al. 2017; Distefano et al. 2017; Boro Saikia et al. 2018; Olspert et al. 2018). One of the shortcomings is clearly the use of the ill-determined convective turnover time τ_c, which is used to the calculated the Rossby number Ro = 4πP_rot/τ_c; for every star τ_c is highly depth-dependent and a different location of a dynamo might invoke a different τ_c. However, these branches can be also obtained using the fractional Ca II H&K emission ${R\prime }_{\text{HK}}$ instead of the Rossby number (see, e.g., Brandenburg et al. 2017; Olspert et al. 2018).

Explaining the observational findings via dynamo models has been challenging. Simple mean-field models of turbulent dynamos produce rotational dependencies of cycle periods similar to the observed cycle periods using overlapping induction layers (Kleeorin et al. 1983). Advective dominated flux transport models have tended to produce an increase of cycle periods with increasing rotation rates (e.g., Dikpati & Charbonneau 1999; Bonanno et al. 2002; Jouve 2010), which is opposite from what is observed. In these models the cycle length is mainly determined by the strength of return flow of the meridional circulation, which is believed to decrease with increasing rotation (e.g., Köhler 1970; Brown et al. 2008; Warnecke et al. 2016; Käpylä et al. 2017; Viviani et al. 2018). However, the models of Kitchatinov & Rüdiger (1999) showed an increase of meridional circulation strength with rotation, leading to a decrease of cycle period with rotation (e.g., Küker et al. 2001), which agrees qualitatively with observations.

Another possibility to explain cycles in dynamo models is via the turbulent (eddy) magnetic diffusivity. In a propagating αΩ dynamo wave, the dynamo drivers, which are responsible for the cycle length, have to balance with the contribution from the turbulent diffusion. Using a turbulent diffusivity of η_t = 2 × 10⁸ m² s⁻¹, a magnetic cycle length of around 23 yr is obtained (e.g., Roberts & Stix 1972), which is pretty close to solar value of 22 yr. A change in the cycle length can be then associated with a change in the turbulent diffusion caused by magnetic or rotational quenching (e.g., Rüdiger et al. 1994). These authors found a cycle dependence of $P_{\text{cycl}}\propto P_{\text{rot}}^{0.1}$.

There have only been a limited number of studies of rotational dependencies of dynamo cycles using global dynamo simulations. Strugarek et al. (2017) found $P_{\text{cycl}}\propto P_{\text{rot}}^{-1.06}$ in a rather limited sample of rotation rates. In the recent work by Viviani et al. (2018), the authors have found no clear dependency of cycle periods with rotation. However, a decrease in cycle period with increasing rotation seems be more likely than an increase, which might be because of the strong oscillatory nonaxisymmetric magnetic fields in these simulations.

In this work, we present the results of spherical convective dynamo models with rotation rates varying by a factor of 30. We determine the cycle dependency on rotational influence in Sect. 1, interpret the data in terms of the Parker–Yoshimura rule with test-field obtained transport coefficients in Sect. 3.2, and compare the finding with observational and other numerical results in Sect. 3.3.

2. Model and setup

The detailed description of the general model can be found in Käpylä et al. (2013) and is not be repeated here. We modeled the convection zone of a solar-like star in spherical polar coordinates (r, θ, φ) using the wedge assumption 0.7 R < r < R, Θ₀ < θ < π− Θ₀ and 0 < φ < π/2 with R being the stellar radius and Θ₀ = 15°. We solved the evolution equations of compressible magnetohydrodynamics for the magnetic vector potential A, which therefore ensures the solenoidality of the magnetic field B = ∇ × A, for the velocity u, specific entropy s and density ρ. The model assumes an ideal gas for the equation of state. The fluid is also influenced by Keplerian gravity and rotation via the Coriolis force. Because of the wedge assumption, we used a periodic boundary condition in the azimuthal (φ) direction. We assume a stress-free condition for the velocity field at other boundaries and perfect conducting latitudinal and bottom boundaries and radial field condition at the top boundary for the magnetic field. The energy is transported into the system via a constant heat flux at bottom boundary and the temperature obeys a black body condition. On the latitudinal boundaries, the energy flux is vanishing using zero derivative for the thermal-dynamical quantities. The detailed setup including the exact equations and expression for the boundary condition can be found in Käpylä et al. (2013, 2017) and Warnecke et al. (2014).

We characterize our runs with the following nondimensional input parameters: the Taylor number, subgrid-scale (SGS), and magnetic Prandtl numbers

$$\begin{array}{ccc}\text{Ta} = {[2{\mathrm{\Omega }}_0{(0.3 R)}^2/v]}^2,& {\text{Pr}}_{\text{SGS}}=\frac{v}{{\chi }_{\mathrm{m}}^{\text{SGS}}},& {\text{Pr}}_{\mathrm{M}}\end{array}=\frac{v}{\eta },$$(1)

where ν and η are the constant kinematic viscosity and magnetic diffusivity, and the SGS heat diffusivity ${\chi }_{\mathrm{m}}^{\text{SGS}}={\chi }_{\text{SGS}}\left(r_{\mathrm{m}}\right)$ is evaluated at r_m = 0.85R. Furthermore, we use the Rayleigh number obtained from the hydrostatic stratification, evolving a onedimensional model, given by

$$\text{Ra=}\frac{GM{\left(0.3R\right)}^4}{v{\chi }_{\mathrm{m}}^{\text{SGS}}{R}^2}{(-\frac{1}{c_{\mathrm{P}}}\frac{\mathrm{d}{s}_{\text{hs}}}{\mathrm{d}r})}_{(r=0.85R)},$$(2)

where s_hs is the hydrostatic entropy. As diagnostic parameters, we quote density contrast as

$${\mathrm{\Gamma }}_{\rho }\equiv \rho (r=0.7R)/\rho \left(R\right),$$(3)

the fluid and magnetic Reynolds numbers and the Péclet number,

$$\begin{array}{ccc}\text{Re=}\frac{{\mu }_{\text{rms}}}{v{k}_f},& {\text{Re}}_{\mathrm{M}}=\frac{{\mu }_{\text{rms}}}{\eta k_f},& \text{Pe=}\frac{{\mu }_{\text{rms}}}{{\chi }_{\mathrm{m}}^{\text{SGS}}{k}_f}\end{array},$$(4)

where k_f = 2π/0.3R ≈ π 21/R is an estimate of the wavenumber of the largest eddies. We define the Coriolis number as

$$\text{Co} = 2{\mathrm{\Omega }}_0{\tau }_{\mathrm{c}}\equiv {\text{Ro}}^{-1},$$(5)

where τ_c = 1/u_rmsk_f is the convective turnover time and $u_{\text{rms}}=\sqrt{(3/2){\langle u_r^2+u_{\theta }^2\rangle }_{r\theta \phi t}}$ is the rms velocity and the subscripts indicate averaging over r, θ, φ and a time interval covering the saturated state. The duration of the saturated state is indicated by τ_sat and covers several magnetic diffusion times. The kinetic and magnetic energy density are given by

$$\begin{array}{cc}{E}_{\text{kin}}=\frac{1}{2}{\langle \rho {\mathit{u}}^2\rangle }_{r\theta \phi t,}& E_{\text{mag}}=\frac{1}{2{\mu }_0}{\langle {\mathit{B}}^2\rangle }_{r\theta \phi t\mathrm{.}}\end{array}$$(6)

All values for these nondimensional input and diagnostic parameters are shown in Table 1 for all runs.

The wedge assumption in the azimuthal (φ) direction allows us to suppress the nonaxisymmetric dynamo mode with azimuthal degree m = 1, 2, 3 and therefore use the mean-field decomposition to describe the large-scale velocity and magnetic field. We use an overbar to refer to the mean, azimuthal averaged quantity and a prime for the fluctuating quantity; for example, B = B̄ + B^′.

To determine some of the turbulent transport coefficients in these simulations, we make use of the test-field method (Schrinner et al. 2005, 2007; Warnecke et al. 2018). This method uses linear independent test-fields, which do not backreact on the flow, to determine the electromotive forces of these test fields using the mean and fluctuating flow fields of the simulations. This allows us to obtain all components of the turbulent transport tensors. In this work, we only use the φφ component of the α tensor.

We present some of the results in physical units with a normalization based on the solar rotation rate Ω_⊙ = 2.7 × 10⁻⁶ s⁻¹, solar radius R = 7 × 10⁸ m, density at the bottom of the convection zone ρ(0.7R) = 200 kg/ m³, and μ₀ = 4π × 10⁻⁷ H m⁻¹. Furthermore, the rotation of the simulations is given in terms of solar rotation rate with Ω̃ ≡ Ω₀/ Ω_Ω. However, the rotational influence on the convection is much better described by the use of the Coriolis number Co. The simulations were performed using the PENCIL CODE¹.

3. Results

For all the simulations we keep all input parameters constant except that we increase the rotation rate for 0.5 to 15 solar rotation rate, corresponding to Co = 0.7 to 40.3. We only lower the diffusivities (ν, η, ${\chi }_{\mathrm{m}}^{\text{SGS}}$) for the run with the highest rotation rate (M15) to keep the Reynolds and Péclet numbers at a similar level, however the Prandtl numbers are kept fixed. We name the runs ‘M’ because of their magnetic nature followed by their solar rotation rate. Run M5 has been discussed as Run I in Warnecke et al. (2014), Run A1 in Warnecke et al. (2016), Run D3 in Käpylä et al. (2017) and Warnecke et al. (2018), and Run G^W in Viviani et al. (2018). We analyzed Run M3 as Run B1 in Warnecke et al. (2016). Runs M10 and M15 are similar to Runs I^W and J^W of Viviani et al. (2018). As calculated in Warnecke et al. (2018), the Rayleigh number for Run M5 is around 100 times the critical value. We expect that this factor increases for lower rotation and decreases for higher rotation.

In this work, we do not discuss all properties of the rotational influence of the hydrodynamical dynamics, i.e., the angular momentum evolution, we instead focus on the discussion and analysis of the dynamo cycles and their possible origin.

Before we do this in detail, we look at the differential rotation generated in these simulation by the interplay of rotation and turbulent convection. In Fig. 1, we show the time averaged differential rotation Ω = Ω₀ + ū/r sin θ for all runs. In agreement with earlier findings (Gastine et al. 2014; Käpylä et al. 2014; Fan & Fang 2014; Karak et al. 2015; Viviani et al. 2018), for slow rotation of Co = 0.7–2.9 the equator is rotating slower than the poles, that is, so-called antisolar differential rotation; for rapid rotation the Co = 3.7–40 the poles are rotating slower than the equator, that is, so-called solar-like differential rotation. The overall relative latitudinal and radial shear is the strongest for the slowest rotation and decreases in the solar-like differential regime for higher rotation rates. This agrees quantitatively with results of observational studies (Reinhold T. et al. 2013; Lehtinen et al. 2016), results of mean-field models (Kitchatinov & Rüdiger 1999) and results of previous global simulations (e.g., Käpylä et al. 2011, 2014; Gastine et al. 2014; Viviani et al. 2018). The difference between the northern and southern hemisphere, for example in Runs M0.5 and M2, is caused by a hemispheric dynamo producing a stronger magnetic field in one hemisphere than in the other; see Sect. 3.1. An interesting feature is the occurrence of a region of a minimum of Ω at mid-latitude in the solar-like differential rotation regime. This region corresponds to area of large negative shear and has been shown to be responsible for the equatorward migrating dynamo wave in Run M5 (Warnecke et al. 2014, 2016, 2018; Käpylä et al. 2017). This region seems to become less pronounced for more rapidly rotating simulations.

Fig. 1. Normalized differential rotation Ω/ Ω0 with Ω = Ω0 + ū/r sin θ for all runs. The value Ω has been calculated as a time average over the saturated state.

3.1. Magnetic cycles

All simulations discussed in this work show large-scale dynamo action. Lowering the rotation rate below Ω̃ = 0.5, produces only a weak large-scale magnetic field (Ω̃ = 0.4) or no dynamo action (Ω̃ < 0.3) for the same parameters. A small-scale dynamo is not present in any of these simulations (Warnecke et al. 2018). In Fig. 2, we show the near-surface mean azimuthal magnetic field B̅_ϕ as a function of time and latitude, the so-called butterfly diagram for all runs. For the slow rotating Run M0.5, we find the dynamo produces a large-scale magnetic field that is most pronounced in one hemisphere and has no polarity reversals; only the amplitude shows weak cyclic variations. For Runs M1 to M2.5 the magnetic field is of chaotic nature and has polarity reversal, which does not seem to be cyclic. This is similar to previous findings for slowly rotating convective dynamos (Fan & Fang 2014; Karak et al. 2015; Hotta et al. 2016; Käpylä et al. 2017). Interestingly, Run M1 shows some indication of an oscillating magnetic field with antisolar differential rotation, however, from the current running time, we cannot draw any certain conclusions. Indications of cyclic solutions in the antisolar regime have been also found by Karak et al. (2015), but only recently Viviani et al. (2018) have obtained clear cyclic solutions with many polarity cycles using the full azimuthal extent. Run M3 shows indication of a cyclic magnetic field that is most pronounced at the poles. The magnetic cycle is even more pronounced in the middle of convection zone, as shown in Warnecke et al. (2016). The Runs M4 to M15 show a clear cyclic magnetic field with regular polarity reversals. The cycle length for these runs seems to be very similar. For Runs M4 to M10 the magnetic field shows a clear equatorward propagation similar to that observed for the solar activity belt. Furthermore, a shorter, much weaker poleward migrating cycle is present in addition to the equatorward migrating mode. The short poleward migrating cycle seems to become stronger for increasing rotation. It has been associated with a local α² dynamo mode in addition to the αΩ dynamo causing the equatorward migration (Käpylä et al. 2016; Warnecke et al. 2018). For Run M15 the field shows indication of both equatorward and poleward migration.

Fig. 2. Mean azimuthal magnetic field B̅φ as a function of time in years and latitude near the surface (r = 0.98R) for all runs. The time interval shows the full duration of the saturated state for Runs M0.5 to M3 and an interval of 32 yr for Runs M4 to M15 to illustrate the similarity in cycle length. The black and white dashed horizontal line indicates the equator.

To quantify the cycle period of the magnetic field of all runs, we calculate the power spectrum of the magnetic field and use the strongest peak as the cycle frequency. In the following we distinguish between magnetic cycle and activity cycle. The magnetic cycle includes a full polarity reversal, corresponding to the 22 yr on the Sun, the activity cycle uses the maximum and minima of the magnetic energy, which corresponds to 11 yr for the Sun. We use two ways to calculate our cycle period as follows: first, we determine the magnetic cycle and take the half and, second, we determine the activity cycle. Both results are shown in Table 1. For the magnetic cycle we take the radial and azimuthal mean magnetic field component and calculate the power spectrum for each latitude and at three radii (r = 0.98, 0.85, 0.72). Then we average the spectra over latitude to obtain six spectra. Now, we determine the frequency of the largest peak of each spectrum with a corresponding error. The error is estimated by taking the full width half maximum unless this value is narrower than the local grid spacing in the frequency space, in which case we take local grid spacing as an error estimate. From the six frequencies, we calculate the weighted average of the corresponding periods. We take the half of the averaged magnetic cycle period and show it as P_cycl with its error E_Pcycl in Table 1.

For the activity cycle, we use the root-mean-squared value of the magnetic field ${\overline{B}}^{rms}=\sqrt{{\overline{B}}_{r^2}+{\overline{B}}_{{\theta }^2}+{\overline{B}}_{{\phi }^2}}$ near the surface and as a function of radius. As above we calculate the spectrum for each latitude and radius and averaged over them. The cycle period is determined in the same way as for magnetic cycle without dividing by two. The results are shown as P̃_cycl with its error Ẽ_Pcycl in Table 1. In Fig. 3, we show as an example the power spectra of Run M7. In addition to the cycle with an activity period of around 2.7 yr, we notice the weak short cycle, which is also visible in Fig. 2.

Fig. 3. Magnetic power spectra for Run M7. Panel a: spectra of the mean radial magnetic field B̄r (solid) and the mean azimuthal magnetic field B̄φ (dashed) near the surface (r = 0.98; black), in the middle of the convection zone (r = 0.72; blue). Panel b: spectra of ${\overline{B}}^{rms}=\sqrt{{\overline{B}}_{r^2}+{\overline{B}}_{{\theta }^2}+{\overline{B}}_{{\phi }^2}}$ near the surface (r = 0.85; red) and averaged over radius (black). All spectra are obtained for each latitude and then averaged. The peaks in a correspond to magnetic cycle periods and in b to activity cycles periods. The solid vertical lines indicate the cycle periods determined from the weighted average of the spectra in a (Pcycl), the dashed lines indicate the cycle periods determined from the weighted average of the spectra in b ( P̃cycl).

For the Runs M4 to M15, the cycle periods can be determined well with a small error in P_cycl. Furthermore, for these runs P̃_cycl agrees very well with P_cycl, even though their errors are higher. The larger errors of P̃_cycl are caused by the summation over phaseshifted magnetic field components and this can result in a less pronounced peak in the spectrum of B̄^rms. This can be also seen in the fact that Ẽ_Pcycl is for all runs significant larger than E_Pcycl. For the slowly rotating runs, the periods are not well determined; there are with significant differences between the two methods and large errors as well. In the following, we therefore focus on the analysis of the runs with a clear cycle period.

In Fig. 4 we show the cycle periods as a function of Coriolis number. We find two clearly separated group of runs: the slowly rotating simulations with long ill-determined cycles (Runs M1 to M3) and the moderately to rapidly rotating runs with short cycles (Runs M4 to M15). In the latter group the cycle periods increase slightly with increasing rotation. We perform a power law fit for theses runs and obtain P_cycl ∝ Co^{0.25 ± 0.04}, or in terms of rotation period $P_{\text{cycl}}\propto P_{\text{rot}}^{-0.33\pm 0.05}$. This value is in disagreement with the observation of Noyes et al. (1984b) and Oláh et al. (2016), but agree qualitatively with scaling of advective dominated fluxtransport dynamo models (e.g., Dikpati & Charbonneau 1999; Bonanno et al. 2002; Jouve 2010). Strugarek et al. (2017) also found an increase of cycle period with increasing rotation, however their power law fit reveals a much steeper increase with rotation $P_{\text{cycl}}\propto P_{\text{rot}}^{-1.06}$ . The cycle period calculated for Run M5 agrees with the cycle periods obtain for similar runs using the D2 phase dispersion statistics and ensemble empirical mode decomposition (Käpylä et al. 2016, 2017).

Fig. 4. Cycle periods as a function of Coriolis number Co showing Pcycl in black and P̃cycl orange. The green dashed line indicates a power law fit of the Runs M4 to M15.

3.2. Cause of magnetic cycles

Earlier studies of similar simulations as Run M5 show that the equatorward migrating mean magnetic field can be well explained with a Parker–Yoshimura (Parker 1955; Yoshimura 1975) αΩ-dynamo wave propagating equatorward (Warnecke et al. 2014, 2016, 2018). Following the calculation of Parker (1955) and Yoshimura (1975), we can compute the cycle frequency of the dynamo wave using (see also Stix 1976)

$${\omega }_{\text{PY}}={|\frac{{\alpha }_{\phi \phi }{k}_{\theta }}{2}r \cos \theta \frac{\partial \mathrm{\Omega }}{\partial r}|}^{1/2},$$(7)

where k_θ is the latitudinal wave number. The corresponding activity cycle period is then given by

$$P_{\text{PY}}=\frac{2\pi }{2{\omega }_{\text{PY}}}\mathrm{.}$$(8)

As pointed out by Warnecke et al. (2014), the location at which we measure the shear and the α_φφ is crucial to get a meaningful result for the direction and therefore the period of the dynamo wave. Following this work, we calculate P_PY in the region where (i) ${\overline{B}}_{\phi }^{\text{rms}}={\left({\langle {\overline{B}}_{\phi }^2\rangle }_t\right)}^{1/2}$ is large, in our case at least larger than the half of the maximum value; (ii) the radial shear ∂Ω/ ∂r is negative; and (iii) α_φφ is positive. The last two criteria are needed to excite an equatorward migrating dynamo wave, following the Parker–Yoshimura sign rule. To make sure that these drivers are really responsible for exciting a dynamo wave at this location, the production of azimuthal magnetic field must be large at this location, leading to criterion (i). The lower limit of half of the maximum value is a reasonable choice; a slightly different value has only little effect on the cycle period determination and on its dependence on rotation rate. The criteria have been also successfully used to confirm Paker-Yoshimura dynamo waves in similar simulations (Warnecke et al. 2014, 2016, 2018; Käpylä et al. 2016). Using these criteria, we then average over these regions. In Fig. 5, we show these regions for Run M1, M3, and M7. To calculate P_PY, we choose

$$k_{\theta }=\frac{1}{R(1-2{\mathrm{\Theta }}_0/\pi ),}=1.2/R,$$(9)

Fig. 5. Mean magnetic field, effect and radial shear profiles for Runs M1, M3, and M7. We show the rms mean azimuthal magnetic field averaged over the saturate state ${\overline{B}}_{\phi }^{\text{rms}}$ (top row), αφφ determined with the test-field method (middle row) and the radial shear r sin θ∂Ω/∂r (bottom row). The dashed lines indicate the region for which we calculate PPY.

where the factor 1 – 2Θ₀/π takes into account the absence of the poles in our simulations. However, the actually value of k_θ only affects the value P_PY with a −1/2 dependency, but not the scaling with rotation.

In Table 1, we list all computed values for P_PY in Col. 11 and these agree well with the values of P_cycl and P̃_cycl for the runs with well-determined cycles (M4 to M15). For oscillatory solution of planetary dynamos, Gastine et al. (2012) also found good agreement between rotational dependency of measured and dynamo wave predicted cycle length. In Fig. 6a, we show for these runs the cycle periods P_cycl and P̃_cycl together with the predicted period P_PY. A power law fit results in P_PY ∝ Co^0.51±0.05, which is close to P_cycl ∝ Co^0.25±0.04. Therefore, the Parker–Yoshimura dynamo wave explains well the weakly dependency of cycle frequency with rotation that we find for the moderately and rapidly rotating simulations. We now go a step further and check which driver of the dynamo wave causes this rotational dependency. For this we plot in Fig. 6b the rotational dependency of the radial shear and α effect in terms of |r cos θ|∂Ω/∂r and α_φφ; as for P_PY both quantities are averaged over the region of interest. The strength of the shear strongly weakens for larger rotation with an estimated scaling of Co^{−1.33±0.18}. For Run M15 the shear in the region is just below zero explaining the mixture of equatorward and poleward migration pattern shown in Fig. 2. For α_φφ, we find an increase with rotation corresponding to a scaling of Co⁰.70±0.25, so much less than linear. The strong decrease in shear causes the cycles to become larger with rotation; assuming a constant α_φφ, shear alone would leading a scaling of P_cycl ∝ Co^0.67±0.09. The α effect, on the other hand, leads to a decrease of cycle length with rotation; that is, P_cycl ∝ Co^{−0.35±0.12}. Because the increase of cycle length due to shear is stronger than the decrease due to the α effect, the resulting cycle length shows only a weak increase with rotation.

Fig. 6. Panel a: comparison of the cycle periods Pcycl (black) and P̃cycl (orange) with predicted cycle periods using a Parker–Yoshimura dynamo wave PPY (blue) for Runs M4 to M15. The dashed blue line indicates a power law fit to PPY. Panel b: contributions to the Parker– Yoshimura dynamo wave containing the radial shear (black line; left y-axis) and αφφ (red; right y-axis) for Runs M4 to M15. The dashed lines indicate the corresponding power law fits.

The surprising issue with the interpretation of the magnetic field evolution as a Parker–Yoshimura dynamo wave is that for runs rotating slower than the M4 (Co = 6.5) it fails. Equation (7) for these runs predict cycle periods of similar length as for the more rapidly rotating runs, but the actual magnetic field shows no clear cyclic evolution. For example, Run M3 shows similar condition for a dynamo wave as in Run M7; there exists a localized region, in which the mean toroidal field is strong, α_φφ is positive, and the shear is negative; see Fig. 5. In the simulations with antisolar differential rotation (Run M0.5 to M2), we find instead a more extended region of strong mean toroidal field, positive α_φφ and negative shear, however strength of the shear and the α effect should be sufficient to excite an αΩ dynamo wave. One of the reasons for the absence of an αΩ dynamo wave can be the larger turbulent magnetic diffusion due to higher convective velocities as shown in Fig. 9. This is in agreement with previous studies of rotating convection in Cartesian boxes (Käpylä et al. 2009). To make a reliable statement, whether a αΩ dynamo is actually operating in these simulations and the reason for not exciting dynamo wave need to be studied in more detail using all the turbulent transport coefficients. We will postpone such a study for the future.

We can now also interpret the scaling of the shear and the α effect in terms of mean-field models (e.g., Krause & Rädler 1980; Rüdiger 1989). From models of differential rotation, one typically finds that the absolute radial and latitudinal differential rotation stays nearly constant for increasing rotation (e.g., Kitchatinov & Rüdiger 1999); this disagrees with our findings. However, we stress that these models take the latitudinal averaged values at the bottom and surface to compute the radial differential rotation, whereas we compute the local radial shear in the region of interest. In mean-field dynamo models α_φφ is related to the mean kinetic helicity $\overline{\mathit{u}\mathit{\prime }\mathbf{\nabla }\mathit{\times }\mathit{u}\mathit{\prime }}$ and therefore can be linear related to the Ω. Taking the convective turnover time τ_c into account as well leads to scaling of α_φφ ∝ Co for weak rotation. (e.g., Krause & Rädler 1980). As shown in Warnecke et al. (2018), approximating the diagonal α components with ${\alpha }_K=-1/3{\tau }_{\mathrm{c}}\overline{\mathit{u}\mathit{\prime }\mathbf{\nabla }\mathit{\times }\mathit{u}\mathit{\prime }}$ is not correct and can lead to the overestimation of α_φφ. Indeed, we find α_φφ in the region of interest depends weaker on rotation as predicted from mean-field models. The overall scaling of α_φφ averaged of the simulations might be different, but the import value of α_φφ determining the cycle period comes from this region. Warnecke et al. (2018) also took into account the nonlinear quenching of the α effect due to magnetic helicity conservation (see Brandenburg & Subramanian 2005, for details) and use the form introduced by Pouquet et al. (1976)$\alpha =-1/3{\tau }_{\mathrm{c}}\overline{\mathit{u}\mathit{\prime }\cdot \mathbf{\nabla }\mathit{\times }\mathit{u}\mathit{\prime }}+1/3{\tau }_{\mathrm{c}}/\rho \overline{\mathbf{\nabla }\mathit{\times }\mathit{b}\mathit{\prime }\mathit{\cdot }\mathit{b}\mathit{\prime }}$, but still could not find an agreement with the actual measured α_φφ; see Figs. 1 and 2 of Warnecke et al. (2018).

Furthermore, we plot the ratio of rotation period and cycle period P_rot/P_cycl over Coriolis number (see Fig. 7). We find a scaling of P_rot/P_cycl ∝ Co^{−0.98±0.04} for the runs with well-determined cycles. This scaling fits well with the cycle period predicted by a dynamo wave. Interestingly, Run M0.5 fit well to this relation, even though we find no polarity reversals there. In the interpretation of stellar observation, P_rot/P_cycl is often used to determine the quenching of the α effect. If one assumes a linear dependency of α and ∂Ω/Ωr on Ω together with Eq. (7), P_rot/P_cycl over Co give an estimate over the rotational quenching of α(e.g., Brandenburg et al. 1998). Moreover, one can go a step further and plot P_rot/P_cycl over magnetic activity, which is related to the surface magnetic field strengths (e.g., Schrijver et al. 1989). In case of dynamo simulations we can instead use the ratio of magnetic and kinetic energy, so-called dynamo efficiency, to mimic the magnetic activity as in Fig. 8. Then the P_rot/P_cycl dependence on magnetic activity can be interpreted as the magnetic quenching of the α effect. However, stellar observations indicate an increase of the α effect with magnetic activity in the inactive and active branch (Brandenburg et al. 1998, 2017; Saar & Brandenburg 1999).

Fig. 7. Ratio of rotation period and cycle period Prot/Pcycl over the Coriolis number Co. The black asterisks indicate Pcycl and the orange asterisks P̃cycl. The blue solid line shows the predicted cycle length from Eqs. (7) and (8), the green dashed line the fit of Fig. 4 and the red dashed line a power law fit of Runs M4 to M15.

Fig. 8. Ratio of rotation period and cycle period Prot/Pcycl over magnetic and kinetic energy ratio. The black asterisks indicate Pcycl and the orange asterisks P̃cycl.The red dashed line indicates a power law fit of Run s M4 to M15.

Fig. 9. Turbulent (eddy) magnetic diffusivity over Coriolis number Co. The diffusivity is determined using an estimate of the turbulence in the simulations ηt = 1/3urms/kf (black) and using the cycle periods as in Eq. (12) (red) and Eq. (13) following Roberts & Stix (1972) (blue). The green dashed line indicates a power law fit of the black asterisk of Runs M4 to M15.

In our simulations, the situation is different. As described above, the radial shear decreases and the α effect increases with higher rotation rate. Therefore, we cannot use P_rot/P_cycl to estimate the quenching of the α effect. The scaling of P_rot/P_cycl ∝ Co^{−0.99±0.04} as shown in Fig. 7 might seem to be expected because we plot rotation rate over rotation rate, but the Coriolis number also includes the strength of convection in terms of τ_c, which is influenced by rotation as well. We do not find any indication of branches with positive slopes similar to the inactive or active branch as postulated by Brandenburg et al. (1998). Furthermore, the slope is different from what is found from the super-active branch, which has a slope of Co^−0.43 (Saar & Brandenburg 1999).

In Fig. 8, we plot P_rot/P_cycl over ratio of magnetic and kinetic energy, which can be interpreted as the dynamo efficiency, we find a scaling of P_rot/P_cycl ∝ (E_mag/E_kin)^{−1.29±0.05}. The scaling seems to agree qualitatively with what is found in Viviani et al. (2018). Also, we do not find any indication of a positive slope and therefore a similarity to the activity branches. It seems clear that the runs with well-determined cycles cannot be interpreted in terms of activity branches with positive slopes. However, our scaling agrees qualitatively with the suggested transitional branch by Distefano et al. (2017).

Another way to analyze the scaling of the cycle period with rotation is via the turbulent eddy magnetic diffusivity. In the saturated state the dynamo drivers has to balance to contribution of the magnetic diffusion. In a αΩ dynamo wave, the balance is written

$${\omega }_{\text{PY}}-k^2{\eta }_{\mathrm{t}}=0,$$(10)

where k is a wavenumber and ω_PY is given by Eq. (7). Therefore, we can use this equation to calculate η_t based on the cycle frequency and compare with estimated values using the turbulent flow of the simulations. Using the first-order smoothing approximation (FOSA; see, e.g., Krause & Rädler 1980) and isotropic and homogeneous turbulence, the turbulent eddy diffusivity can be estimated as

$${\eta }_{\mathrm{t}}=\frac{1}{3}\frac{{\mu }_{\text{rms}}}{k_{\mathrm{f}}}\mathrm{.}$$(11)

If we assume k = k_θ in Eq. (10), we can calculate η_t based on the cycle period

$${\eta }_{\mathrm{t}}=\frac{2\pi }{2 P_{\text{cycl}}{k}_{\theta }^2}\mathrm{.}$$(12)

Another way is to use the radius R and the depth of the convection zone 0.3R of the star in relation η_t with the cycle period (Roberts & Stix 1972), i.e.,

$${\eta }_t=\frac{0.3R^2}{2P_{\text{cycl}}}\mathrm{.}$$(13)

We show the values for both expressions for the runs with welldetermined cycles (Runs M4 to M15) together with values of Eq. (11) as a function of Coriolis number in Fig. 9. The value η_t of Eq. (13) fits remarkable well with Eq. (11). Even the scaling of η_t = Co^{−0.27±0.02} determined from Eq. (11) is the same as expected from the cycle periods η_t ∝ 1/P_cycl ∝ Co^{−0.25±0.04} (see Fig. 4). This good agreement is indeed interesting and not fully expected because the estimation of η_t in Eq. (11) is based on strong assumptions, which are most likely not fulfilled in these simulations. The values obtained through Eq. (12) have obviously the same scaling as Eq. (13), however the values are around a factor of ten higher. This is because of the different values of scales and wave numbers included in this calculation. If we use k_θ instead of k_f in Eq. (11), the curves would lie closer together. Therefore, the scaling of cycle period with rotation rate of P_cycl ∝ Co^0.25±0.04 can also be explained well by the rotational quenching of the turbulent (eddy) magnetic diffusivity.

3.3. Comparison with observational and other numerical studies

There is only limited amount of numerical studies of cycles of solar and stellar dynamos. In the following, we compare our results with the recent work of Strugarek et al. (2017) and Viviani et al. (2018). The study by Strugarek et al. (2017) include seven models, all showing cyclic dynamo solutions. Even though the authors interpret their models in the vicinity of the Sun, their Rossby numbers indicate a rapid rotational regime. If we convert their numbers to our definition of Coriolis numbers or inverse Rossby number as defined in Brandenburg et al. (1998), respectively, we obtain values of Co = 18 to 62. These high numbers are due to their low convective velocities compared to other models (Käpylä et al. 2017). For comparison the estimated Coriolis number of the Sun is Co = 6.2 (e.g., Brandenburg et al. 2017). In Fig. 10, we plot the models of Strugarek et al. (2017) together with our models. There seems to be no overlap between their and our simulations. However, their simulations also show a decrease of P_rot/P_cycl with rotation following an even steeper slope. The similar slope might be caused by a similar dynamo mechanism and the small difference and the offset might be because of different system parameters, as Prandtl numbers and/or Reynolds numbers. As shown in Käpylä et al. (2017), the Reynolds numbers in typical simulations with the EULAC code can be lower compared to models of other, similar codes. This might also explain the dominantly axisymmetric large-scale magnetic field solution in Strugarek et al. (2017). The simulations of Viviani et al. (2018) show clearly that there is a transition from an axisymmetric to nonaxisymmetric large-scale magnetic field solution at Co ≥ 3. However, if the resolution, therefore the Reynolds and Rayleigh numbers, are not high enough, the nonaxisymmetric magnetic solution cannot be obtained. This is most important for large rotation rates, where the convection is rotationally quenched. Furthermore, Strugarek et al. (2017) have claimed that their dynamo cycle is caused by a nonlinear feedback of the torsional oscillation on the magnetic field. Such an interpretation is not very likely to be correct because a cyclic torsional oscillation requires a cyclic dynamo in the first place. In the light of the results of this work we are inclined to think that their cyclic magnetic field is also caused by a Parker–Yoshimura dynamo wave. Indeed, their differential rotation profiles show localized regions of strong negative shear, where the mean magnetic field propagates equatorward as well.

Fig. 10. Ratio of rotation period and cycle period Prot/Pcycl over the Coriolis number Co. The black asterisks indicate Pcycl. The dashed red line indicates the fit Co−0.98 of Fig. 7 and the dashed black line the transition from antisolar to solar like differential rotation. We include the simulations of Viviani et al. (2018; blue squares, with crosses for wedge runs) and Strugarek et al. (2017; green diamond), the observational studies of Lehtinen et al. (2016; purple crosses) and of Brandenburg et al. (2017; light gray triangles for K dwarfs, dark gray circles for F, G dwarfs, including the Sun: yellow cross).

Viviani et al. (2018) has probed a large range of rotation rates in particular in the rapid rotational regime. Their simulations using a similar setup as in this work, but for most of their simulations they use a full 2π extend in the azimuthal direction (2π runs) and obtain nonaxisymmetric large-scale field solutions for moderately to rapidly rotating runs. In Fig. 10, we also overplot their simulations. Their π/2 wedge runs agrees well with our runs and our obtained scaling of Co^{−0.98±0.04}. This also means that their estimates of cycle periods based on the temporal variation in the large-scale magnetic energy seems to describe the magnetic cycle well. Their 2π runs with rapid rotation (Co ≥ 3) clearly show a different scaling that is similar to Co^−0.43 of the super-active branch (Saar & Brandenburg 1999; Viviani et al. 2018). This might mean that the nonaxisymmetric large-scale magnetic field solutions have a different scaling than the pure axisymmetric magnetic field solutions. The argument of different scaling is also supported by the fact that the slowly rotating simulations with axisymmetric large-scale magnetic field solutions of Viviani et al. (2018) can be well described by the cycle period scaling estimated in this work (see Fig. 10). However, surprising is the fact that for nearly the same rotation rate the 2π runs show much shorter and much clearer cycles than their corresponding π/2 wedge simulations.

Finally, we compare our results with observational obtained stellar cycles. One sample comes from Lehtinen et al. (2016), where the authors use photometry to measure cyclic variation in solar-like stars. From this sample, we only plot the cycles that are identified as better than poor and plot these cycles in Fig. 10. The cycles of rapid rotating stars are surprisingly well placed on our scaling relation even though their magnetic field is nonaxisymmetric. We note that Lehtinen et al. (2016) did not do any direct measurements of the magnetic field, but inferred the degree of nonaxisymmetry from the spot distributions. The cycles of slowly rotating stars seems to fall on a parallel line with a similar scaling. Lehtinen et al. (2016) found that the crossover from the transitional branch to the super-active branch happens at around Co = 26.3 and chromospheric activity value of ${R\prime }_{\text{HK}}=-4.4$. We further include the subsample of stars from the Mount Wilson sample analyzed by Brandenburg et al. (2017). Our scaling relation falls through these stars; however, a lot of stars are not captured by our scaling. As shown by Brandenburg et al. (2017), the stars around Co = 10 form the inactive branch with a positive slope and the stars with higher rotation form the active branch with a positive slope as well. The active branch is not as confined as the inactive branch. In recent works the existence of these branches have been questioned (Reinhold et al. 2017; Distefano et al. 2017). Two new studies of the Mount Wilson sample data (Boro Saikia et al. 2018; Olspert et al. 2018) only find an indication of an inactive branch and otherwise a distribution similar to our scaling showing indication of an transitional branch. There, the authors reanalyze the full Mount Wilson sample without relying of the cycle determination of Baliunas et al. (1995). The problem with determining cycles from chromospheric and photospheric activity time series is due to the method used for the period search. For stars in the inactive branch the cycles are clean and can be well determined, but more active stars with higher rotation rates exhibit a complex behavior with multiple cycles (e.g., Oláh et al. 2016). There, the calculated cycle period can differ depending on which method is used (e.g., Olspert et al. 2018). We also note that a one-to-one comparison with observation is not always possible. Even though using the Coriolis number is more meaningful than then the rotational period, because it indicates whether the star or simulation is in the slow or rapidly rotating regime, the convective turnover time is usually ill determined for observed stars. Furthermore, as we know the turnover time, and therefore the corresponding Coriolis number in the Sun changes several orders of magnitude from the solar surface to the bottom of convection zone (e.g., Stix 2002), it is difficult even for the Sun to estimate a single number as a meaningful Coriolis number. Therefore, Brandenburg et al. (2017) and Olspert et al. (2018) have used the chromospheric activity instead of the Coriolis number for their analysis. Using this, Olspert et al. (2018) have found that the cycle periods of this work fit very well with their observed stellar cycles, see their Fig. 6.

Interestingly the cycle data taken from Brandenburg et al. (2017) have indicated that the Sun lies close to our scaling relation, actually close to Run M4. If we go a step further and assume the 22 yr magnetic cycle of the Sun is caused by a Parker–Yoshimura dynamo wave with a cycle-rotation scaling similar to our simulations, we can estimate the corresponding Coriolis number and then with the definition of Eq. (5) we also calculate the corresponding value of u_rms. We calculate a Coriolis number of Co = 8.5 corresponding to u_rms = 21.5 m s⁻¹. This velocity would be located at around 160 Mm depth (r = 0.78 R) according to the mixing length model of Spruit (1974) or r = 0.72 R for the model of Stix (2002). Therefore, this kind of dynamo wave cannot be driven by the near-surface shear layer, it must instead be driven by a positive radial shear and an inversion of sign of α_φφ in the deeper part of the convection zone to get an equatorward migrating magnetic field (Duarte et al. 2016).

4. Conclusions

We use three-dimensional MHD global dynamo simulations to investigate the rotational dependency of magnetic activity cycles. For moderately and rapidly rotating runs (Co = 6.5), we find well-defined cycles in range between 2 yr and 4 yr. For slowly rotating runs, we find irregular cycles with mostly longer periods. There the cycle periods can only be ill-determined. Using the φφ component of the α tensor measured with the test-field method along with the radial shear, we find a good agreement of the cycle period predicted by a Parker–Yoshimura dynamo wave for moderately to rapid rotating runs. There we find that the cycle period only weakly depends on rotation (P_cycl ∝ Co^0.25±0.04 and $P_{\text{cycl}}\propto P_{\text{rot}}^{-0.33\pm 0.05}$). Also this scaling is well reproduced by a Parker–Yoshimura dynamo wave. The value α_φφ increases only weakly with rotation ( α_φφ ∝ Co^0.70±0.25) and the strength of negative radial shear decreases larger than linear with rotation ∂Ω/∂r ∝ Co^{−1.33±0.18}. This is not in agreement with the firstorder smoothing approximation (FOSA) of mean-field theory, where α depends linear on the rotation for small rotation rates (e.g., Krause & Rädler 1980). Also models of differential rotation show that the radial shear does not have a strong dependency on rotation (Kitchatinov & Rüdiger 1999). However, these models usually look at the global quantities and we determine our scaling from the localized region responsible for driving the dynamo wave.

Looking at the ratio of rotation and cycle period over Coriolis number, we do not find any indication of activity branches with positive slopes as found from observation of stellar cycles by Brandenburg et al. (1998) and Saar & Brandenburg (1999). The negative slope of our simulations with well-determined cycles seems to be more in agreement with the transitional branch postulated by Distefano et al. (2017) and confirmed by Boro Saikia et al. (2018) and Olspert et al. (2018). Furthermore, our results suggest that the cyclic magnetic fields found in the work by Strugarek et al. (2017) are also caused by a Parker–Yoshimura dynamo wave because indeed their simulations produce strong negative shear in the location where the magnetic field is oscillating. By assuming the solar magnetic cycle is caused by a Parker–Yoshimura dynamo wave following a similar scaling as in our simulations, we can conclude that the dynamo in the Sun operates near the bottom of convection zone, where turbulent velocities are around 20 m s⁻¹.

The cycle period dependence on rotation of our simulations can be also well explained via the rotational quenching of the turbulent diffusivity, whose contribution has to balance with the dynamo driver in the saturated stage. We find the nearly the same scaling for the turbulent eddy diffusivity with Coriolis number as expected from a direct inverse proportionality of the diffusivity with cycle period (e.g., Roberts & Stix 1972).

For the slowly rotating runs a Parker–Yoshimura dynamo wave does not seem to be excited, as the predicted periods do not fit with the measured periods. This might be due to a higher turbulent diffusion in this rotation regime as found in Käpylä et al. (2009). In particular, these simulations need to be further investigated using the full set of turbulent transport coefficients. Moreover, we note that in the rapidly rotating regime large-scale magnetic fields of stars can become highly nonaxisymmetric as found recently in observations (e.g., Lehtinen et al. 2016) and simulations (Viviani et al. 2018). The work of Viviani et al. (2018) have indicated a weaker scaling of P _rot/P _cycl with Coriolis number than we find in our work and this might be due to the nonaxisymmetric large-scale magnetic field solution, which are suppressed in our work because of the wedge assumption.

We stress that the cycles determined from our simulation are directly linked to the magnetic field. Observational stellar cycles are mostly measured from chromospheric activity variations (Ca II H&K) or photometry. These are only proxies of the magnetic field strength and might not capture all the features of cyclic variations. Therefore, it would be useful to determine the integrated variability caused by the simulated magnetic cycles as well. Furthermore, we will in future also investigate how coronal heating and therefore the X-ray luminosity will depend on the cyclic magnetic field. For this it is crucial to combine convective dynamo models with a coronal envelope as in Warnecke et al. (2011, 2012, 2013, 2016). This is in particular important to study the role of helicity connecting the dynamo active stellar convection zones with stellar coronae, as the magnetic helicity might play an important role in the heating of coronae (Warnecke et al. 2017).

Acknowledgments

We thank the referee Günther Rüdiger and our colleagues Maarit J. Käpylä, Mariangela Viviani, and Jyri J. Lehtinen for comments on the manuscript and discussion leading to this work. The simulations have been carried out on supercomputers at GWDG, on the Max Planck supercomputer at RZG in Garching, in the facilities hosted by the CSC—IT Center for Science in Espoo, Finland, which are financed by the Finnish ministry of education. J.W. acknowledges funding by the Max-Planck/Princeton Center for Plasma Physics and from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement No. 623609.

References

All Tables

All Figures

	Fig. 1. Normalized differential rotation Ω/ Ω0 with Ω = Ω0 + ū/r sin θ for all runs. The value Ω has been calculated as a time average over the saturated state.
In the text

	Fig. 2. Mean azimuthal magnetic field B̅φ as a function of time in years and latitude near the surface (r = 0.98R) for all runs. The time interval shows the full duration of the saturated state for Runs M0.5 to M3 and an interval of 32 yr for Runs M4 to M15 to illustrate the similarity in cycle length. The black and white dashed horizontal line indicates the equator.
In the text

Fig. 3. Magnetic power spectra for Run M7. Panel a: spectra of the mean radial magnetic field B̄r (solid) and the mean azimuthal magnetic field B̄φ (dashed) near the surface (r = 0.98; black), in the middle of the convection zone (r = 0.72; blue). Panel b: spectra of ${\overline{B}}^{rms}=\sqrt{{\overline{B}}_{r^2}+{\overline{B}}_{{\theta }^2}+{\overline{B}}_{{\phi }^2}}$ near the surface (r = 0.85; red) and averaged over radius (black). All spectra are obtained for each latitude and then averaged. The peaks in a correspond to magnetic cycle periods and in b to activity cycles periods. The solid vertical lines indicate the cycle periods determined from the weighted average of the spectra in a (Pcycl), the dashed lines indicate the cycle periods determined from the weighted average of the spectra in b ( P̃cycl).

In the text

	Fig. 4. Cycle periods as a function of Coriolis number Co showing Pcycl in black and P̃cycl orange. The green dashed line indicates a power law fit of the Runs M4 to M15.
In the text

	Fig. 5. Mean magnetic field, effect and radial shear profiles for Runs M1, M3, and M7. We show the rms mean azimuthal magnetic field averaged over the saturate state ${\overline{B}}_{\phi }^{\text{rms}}$ (top row), αφφ determined with the test-field method (middle row) and the radial shear r sin θ∂Ω/∂r (bottom row). The dashed lines indicate the region for which we calculate PPY.
In the text

Fig. 6. Panel a: comparison of the cycle periods Pcycl (black) and P̃cycl (orange) with predicted cycle periods using a Parker–Yoshimura dynamo wave PPY (blue) for Runs M4 to M15. The dashed blue line indicates a power law fit to PPY. Panel b: contributions to the Parker– Yoshimura dynamo wave containing the radial shear (black line; left y-axis) and αφφ (red; right y-axis) for Runs M4 to M15. The dashed lines indicate the corresponding power law fits.

In the text

	Fig. 7. Ratio of rotation period and cycle period Prot/Pcycl over the Coriolis number Co. The black asterisks indicate Pcycl and the orange asterisks P̃cycl. The blue solid line shows the predicted cycle length from Eqs. (7) and (8), the green dashed line the fit of Fig. 4 and the red dashed line a power law fit of Runs M4 to M15.
In the text

	Fig. 8. Ratio of rotation period and cycle period Prot/Pcycl over magnetic and kinetic energy ratio. The black asterisks indicate Pcycl and the orange asterisks P̃cycl.The red dashed line indicates a power law fit of Run s M4 to M15.
In the text

	Fig. 9. Turbulent (eddy) magnetic diffusivity over Coriolis number Co. The diffusivity is determined using an estimate of the turbulence in the simulations ηt = 1/3urms/kf (black) and using the cycle periods as in Eq. (12) (red) and Eq. (13) following Roberts & Stix (1972) (blue). The green dashed line indicates a power law fit of the black asterisk of Runs M4 to M15.
In the text

Fig. 10. Ratio of rotation period and cycle period Prot/Pcycl over the Coriolis number Co. The black asterisks indicate Pcycl. The dashed red line indicates the fit Co−0.98 of Fig. 7 and the dashed black line the transition from antisolar to solar like differential rotation. We include the simulations of Viviani et al. (2018; blue squares, with crosses for wedge runs) and Strugarek et al. (2017; green diamond), the observational studies of Lehtinen et al. (2016; purple crosses) and of Brandenburg et al. (2017; light gray triangles for K dwarfs, dark gray circles for F, G dwarfs, including the Sun: yellow cross).

In the text