---
abstract: 'We present results from simulations of rotating magnetized turbulent convection in spherical wedge geometry representing parts of the latitudinal and longitudinal extents of a star. Here we consider a set of runs for which the density stratification is varied, keeping the Reynolds and Coriolis numbers at similar values. In the case of weak stratification, we find quasi-steady dynamo solutions for moderate rotation and oscillatory ones with poleward migration of activity belts for more rapid rotation. For stronger stratification, the growth rate tends to become smaller. Furthermore, a transition from quasi-steady to oscillatory dynamos is found as the Coriolis number is increased, but now there is an equatorward migrating branch near the equator. The breakpoint where this happens corresponds to a rotation rate that is about 3–7 times the solar value. The phase relation of the magnetic field is such that the toroidal field lags behind the radial field by about $\pi/2$, which can be explained by an oscillatory $\alpha^2$ dynamo caused by the sign change of the $\alpha$-effect about the equator. We test the domain size dependence of our results for a rapidly rotating run with equatorward migration by varying the longitudinal extent of our wedge. The energy of the axisymmetric mean magnetic field decreases as the domain size increases and we find that an $m=1$ mode is excited for a full $2\pi$ azimuthal extent, reminiscent of the field configurations deduced from observations of rapidly rotating late-type stars.'
author:
- 'Petri J. Käpylä$^{1,2}$, Maarit J. Mantere$^{3,1}$, Elizabeth Cole$^{1}$, Jörn Warnecke$^{2,4}$ and Axel Brandenburg$^{2,4}$'
bibliography:
- 'paper.bib'
title: Effects of enhanced stratification on equatorward dynamo wave propagation
---

Introduction
============

The large-scale magnetic field of the Sun, manifested by the 11 year sunspot cycle, is generally believed to be generated within or just below the turbulent convection zone [e.g., @O03 and references therein]. The latter concept is based on the idea that strong shear in the tachocline near the bottom of the convection zone amplifies the toroidal magnetic field which then becomes buoyantly unstable and erupts to the surface [e.g., @Pa55a]. This process has been adopted in many mean-field models of the solar cycle in the form of a non-local $\alpha$-effect [e.g., @KO11], which is based on early ideas of [@Bab61] and [@Lei69] that the source term for poloidal field can be explained through the tilt of active regions. Such models assume a reduced turbulent diffusivity within the convection zone and a single cell anti-clockwise meridional circulation which acts as a conveyor belt for the magnetic field. These so-called flux transport models [e.g., @DC99] are now widely used to study the solar cycle and to predict its future course [@DG06; @CCJ07].

The flux transport paradigm is, however, facing several theoretical challenges: $10^5$gauss magnetic fields are expected to reside in the tachocline [@DSC93], but such fields are difficult to explain with dynamo theory [@GK11] and may have become unstable at much lower field strengths [@ASR05]. Furthermore, flux transport dynamos require a rather low value of the turbulent diffusivity within the convection zone [several $10^{11}\,{\rm cm}^2\,{\rm s}^{-1}$; see @BERB02], which is much less than the standard estimate of several $10^{12}\,{\rm cm}^2\,{\rm s}^{-1}$ based on mixing length theory, which, in turn, is also verified numerically [e.g., @KKB09a]. Several other issues have already been addressed within this paradigm, for example, the parity of the dynamo [@BERB02; @CNC04; @DdTGAW04] and the possibility of a multicellular structure of the meridional circulation [@JB07], which may be more complicated than that required in the flux transport models [@Ha11; @MFRT12; @ZBKDH13]. These difficulties have led to a revival of the distributed dynamo [e.g., @Br05; @Pi13] in which magnetic fields are generated throughout the convection zone due to turbulent effects [e.g., @KR80; @KKT06; @PS09].

Early studies of self-consistent three-dimensional magnetohydrodynamic (MHD) simulations of convection in spherical coordinates produced oscillatory large-scale dynamos [@Gi83; @Gl85], but the dynamo wave was found to propagate toward the poles rather than the equator—as in the Sun. These models are referred to as direct numerical simulations (DNS), i.e., all operators of viscous and diffusive terms are just the original ones, but with vastly increased viscosity and diffusivity coefficients. More recent anelastic large-eddy simulations (LES) with rotation rates somewhat higher than that of the Sun have produced non-oscillatory [@BBBMT10] and oscillatory [@BMBBT11; @NBBMT13] large-scale magnetic fields, depending essentially on the rotation rate and the vigor of the turbulence. However, similar models with the solar rotation rate have either failed to produce an appreciable large-scale component [@BMT04] or, more recently, oscillatory solutions with almost no latitudinal propagation of the activity belts [@GCS10; @RCGBS11]. These simulations covered a full spherical shell and used realistic values for solar luminosity and rotation rate, necessitating the use of anelastic solvers and spherical harmonics [e.g., @BMT04] or implicit methods [e.g. @GCS10]. Here we exploit an alternative approach by modeling fully compressible convection in wedge geometry [see also @RC01] with a finite-difference method. We omit the polar regions and usually cover only a part of the longitudinal extent, e.g., $90\degr$ instead of the full $360\degr$. At the cost of omitting connecting flows across the poles and introducing artificial boundaries there, the gain is that higher spatial resolution can be achieved. Furthermore, retaining the sound waves can be beneficial when considering possible helio- or asteroseismic applications. Our model is a hybrid between DNS and LES in that we supplement the thermal energy flux by an additional subgrid scale (SGS) term to stabilize the scheme and to further reduce the radiative background flux. Recent hydrodynamic [@KMB11; @KMGBC11] and MHD [@KKBMT10] studies have shown that this approach produces results that are in accordance with fully spherical models. Moreover, the first turbulent dynamo solution with solar-like migration properties of the magnetic field was recently obtained using this type of setup [@KMB12a]. Extended setups that include a coronal layer as a more realistic upper radial boundary have been successful in producing dynamo-driven coronal ejections [@WKMB12]. As we show in a companion paper [@WKMB13], a solar-like differential rotation pattern might be another consequence of including an outer coronal layer.

Here we concentrate on exploring further the recent discovery of equatorward migration in spherical wedge simulations [@KMB12a]. In particular, we examine a set of runs for which the rotational influence on the fluid, measured by the Coriolis number, which is also called the inverse Rossby number, is kept approximately constant while the density stratification of the simulations is gradually increased.

The model {#sec:model}
=========

Our model is the same as that in [@KMB12a]. We consider a wedge in spherical polar coordinates, where $(r,\theta,\phi)$ denote radius, colatitude, and longitude. The radial, latitudinal, and longitudinal extents of the wedge are $r_0
\leq r \leq R$, $\theta_0 \leq \theta \leq \pi-\theta_0$, and $0 \leq
\phi \leq \phi_0$, respectively, where $R$ is the radius of the star and $r_0=0.7\,R$ denotes the position of the bottom of the convection zone. Here we take $\theta_0=\pi/12$ and in most of our models we use $\phi_0=\pi/2$, so we cover a quarter of the azimuthal extent between $\pm75^\circ$ latitude. We solve the compressible hydromagnetic equations[^1], $$\frac{{\partial}\bm A}{{\partial}t} = {\bm u}\times{\bm B} - \mu_0\eta {\bm J},$$ $$\frac{D \ln \rho}{Dt} = -\bm\nabla\cdot\bm{u},$$ $$\frac{D\bm{u}}{Dt} = \bm{g} -2\bm\Omega_0\times\bm{u}+\frac{1}{\rho}
\left(\bm{J}\times\bm{B}-\bm\nabla p
+\bm\nabla \cdot 2\nu\rho\bm{\mathsf{S}}\right),$$ $$T\frac{D s}{Dt} = \frac{1}{\rho}\left[-\bm\nabla \cdot
\left({\bm F^{\rm rad}}+ {\bm F^{\rm SGS}}\right) +
\mu_0 \eta {\bm J}^2\right] +2\nu \bm{\mathsf{S}}^2,
\label{equ:ss}$$ where ${\bm A}$ is the magnetic vector potential, $\bm{u}$ is the velocity, ${\bm B} =\bm\nabla\times{\bm A}$ is the magnetic field, ${\bm J} =\mu_0^{-1}\bm\nabla\times{\bm B}$ is the current density, $\mu_0$ is the vacuum permeability, $D/Dt = {\partial}/{\partial}t + \bm{u} \cdot \bm\nabla$ is the advective time derivative, $\rho$ is the density, $\nu$ is the kinematic viscosity, $\eta$ is the magnetic diffusivity, both assumed constant, $${\bm F^{\rm rad}}=-K{\boldsymbol{\nabla}} T\quad\mbox{and}\quad
{\bm F^{\rm SGS}} =-{\chi_{\rm SGS}}\rho  T{\boldsymbol{\nabla}} s$$ are radiative and SGS heat fluxes, where $K$ is the radiative heat conductivity and ${\chi_{\rm SGS}}$ is the turbulent heat conductivity, which represents the unresolved convective transport of heat and was referred to as $\chi_{\rm t}$ in [@KMB12a], $s$ is the specific entropy, $T$ is the temperature, and $p$ is the pressure. The fluid obeys the ideal gas law with $p=(\gamma-1)\rho e$, where $\gamma=c_{\rm P}/c_{\rm V}=5/3$ is the ratio of specific heats at constant pressure and volume, respectively, and $e=c_{\rm V} T$ is the specific internal energy. The rate of strain tensor $\bm{\mathsf{S}}$ is given by $$\mathsf{S}_{ij}=\onehalf(u_{i;j}+u_{j;i})
-\onethird \delta_{ij}\bm\nabla\cdot\bm{u},$$ where the semicolons denote covariant differentiation [@MTBM09].

The gravitational acceleration is given by $\bm{g}=-GM\hat{\bm{r}}/r^2$, where $G$ is the gravitational constant, $M$ is the mass of the star (without the convection zone), and $\hat{\bm{r}}$ is the unit vector in the radial direction. Furthermore, the rotation vector $\bm\Omega_0$ is given by $\bm\Omega_0=(\cos\theta,-\sin\theta,0)\Omega_0$.

Initial and boundary conditions {#sec:initcond}
-------------------------------

The initial state is isentropic and the hydrostatic temperature gradient is given by $$\frac{{\partial}T}{{\partial}r}=-\frac{GM/r^2}{c_{\rm V}(\gamma-1)(n_{\rm ad}+1)},$$ where $n_{\rm ad}=1.5$ is the polytropic index for an adiabatic stratification. We fix the value of ${\partial}T/{\partial}r$ on the lower boundary. The density profile follows from hydrostatic equilibrium. The heat conduction profile is chosen so that radiative diffusion is responsible for supplying the energy flux in the system, with $K$ decreasing more than two orders of magnitude from bottom to top [@KMB11]. We do this by choosing a variable polytropic index $n=2.5\,(r/r_0)^{-15}-1$, which equals 1.5 at the bottom of the convection zone and approaches $-1$ closer to the surface. This means that $K=(n+1)K_0$ decreases toward the surface like $r^{-15}$ such that most of the flux is carried by convection [@BCNS05]. Here, $K_0$ is a constant that will be defined below.

Our simulations are defined by the energy flux imposed at the bottom boundary, $F_{\rm b}=-(K {\partial}T/{\partial}r)|_{r=r_0}$ as well as the values of $\Omega_0$, $\nu$, $\eta$, and ${\overline{\chi}_{\rm SGS}}={\chi_{\rm SGS}}(r_{\rm m}=0.85\, R)$. Furthermore, the radial profile of ${\chi_{\rm SGS}}$ is piecewise constant above $r>0.75R$ with ${\chi_{\rm SGS}}={\overline{\chi}_{\rm SGS}}$ at $0.75 < r <0.98$, and ${\chi_{\rm SGS}}=12.5{\overline{\chi}_{\rm SGS}}$ above $r=0.98R$. Below $r=0.75R$, ${\chi_{\rm SGS}}$ tends smoothly to zero; see Fig. 1 of [@KMB11].

The radial and latitudinal boundaries are assumed to be impenetrable and stress free, i.e., $$\begin{aligned}
&&\!\!\!
u_r=0,\quad \frac{{\partial}u_\theta}{{\partial}r}=\frac{u_\theta}{r},\quad \frac{{\partial}u_\phi}{{\partial}r}=\frac{u_\phi}{r} \quad (r=r_0, R),\\
&&\!\!\!
\frac{{\partial}u_r}{{\partial}\theta}=u_\theta=0,\quad \frac{{\partial}u_\phi}{{\partial}\theta}=u_\phi \cot \theta \quad (\theta=\theta_0,\pi-\theta_0).
\quad\end{aligned}$$ For the magnetic field we assume perfect conductors on the latitudinal and lower radial boundaries, and radial field on the outer radial boundary. In terms of the magnetic vector potential these translate to $$\begin{aligned}
&&\!\!\!
\frac{{\partial}A_r}{{\partial}r}= A_\theta=A_\phi =0 \,\quad
(r=r_0),\\
&&\!\!\!
A_r=0, \;\; \frac{{\partial}A_{\theta}}{{\partial}r}=-\frac{A_{\theta}}{r},\;\; \frac{{\partial}A_{\phi}}{{\partial}r}=-\frac{A_{\phi}}{r} \quad (r=R),\quad\\
&&\!\!\!
A_r=\frac{{\partial}A_\theta}{{\partial}\theta}=A_\phi=0 \quad
(\theta=\theta_0,\pi-\theta_0).\end{aligned}$$ We use small-scale low amplitude Gaussian noise as initial condition for velocity and magnetic field. On the latitudinal boundaries we assume that the density and entropy have vanishing first derivatives, thus suppressing heat fluxes through the boundaries.

On the upper radial boundary we apply a black body condition $$\sigma T^4  = -K\nabla_r T - {\chi_{\rm SGS}}\rho T \nabla_r s,
\label{eq:bbb}$$ where $\sigma$ is the Stefan–Boltzmann constant. We use a modified value for $\sigma$ that takes into account that both surface temperature and energy flux through the domain are larger than in the Sun. The value of $\sigma$ can be chosen so that the flux at the surface carries the total luminosity through the boundary in the initial non-convecting state. However, in many cases we have changed the value of $\sigma$ during runtime to speed up thermal relaxation.

Dimensionless parameters {#sec:dimless}
------------------------

To facilitate comparison with other work using different normalizations, we present our results by normalizing with physically meaningful quantities. We note, however, that in the code we used non-dimensional quantities by choosing $$R = GM = \rho_0 = c_{\rm P} = \mu_0 = 1,$$ where $\rho_0$ is the initial density at $r=r_0$. The units of length, time, velocity, density, entropy, and magnetic field are therefore $$\begin{aligned}
&[x] = R,\,\,\, [t]=\sqrt{R^3/GM},\,\,\, [u] = \sqrt{GM/R},& \nonumber
\\ &[\rho] =\rho_0,\,\,\, [s] = c_{\rm P},\,\,\, [B]= \sqrt{\rho_0
  \mu_0 GM/R}.&\end{aligned}$$ The radiative conductivity is proportional to $K_0=(\mathcal{L}/4\pi)c_{\rm V}(\gamma-1)(n_{\rm ad}+1)\rho_0\sqrt{GMR}$, where $\mathcal{L}$ is the non-dimensional luminosity, given below. The corresponding nondimensional input parameters are the luminosity parameter $$\mathcal{L} = \frac{L_0}{\rho_0 (GM)^{3/2} R^{1/2}},$$ the normalized pressure scale height at the surface, $$\xi = \frac{(\gamma-1) c_{\rm V}T_1}{GM/R},$$ with $T_1$ being the temperature at the surface, the Taylor number $${{\rm Ta}}=(2\Omega_0 R^2/\nu)^2,$$ the fluid and magnetic Prandtl numbers $${{\rm Pr}}=\frac{\nu}{\chi_{\rm m}},\quad {{\rm Pr}}_{\rm SGS}=\frac{\nu}{{\overline{\chi}_{\rm SGS}}},\quad {{\rm Pm}}=\frac{\nu}{\eta},$$ where $\chi_{\rm m}=K/c_{\rm P} \rho_{\rm m}$ and $\rho_{\rm m}$ are the thermal diffusivity and density at $r=r_{\rm m}$, respectively. Finally, we have the non-dimensional viscosity $$\tilde{\nu}=\frac{\nu}{\sqrt{GMR}}.$$ Instead of $\xi$, we often quote the initial density contrast, $\Gamma_\rho^{(0)}\equiv\rho(r_0)/\rho(R)$. The density contrast can change during the run. We list the final values of $\Gamma_\rho$ from the thermally saturated stage in Table \[tab:runs\].

Other useful diagnostic parameters are the fluid and magnetic Reynolds numbers $${{\rm Re}}=\frac{{u_{\rm rms}}}{\nu {k_{\rm f}}},\quad {{\rm Rm}}=\frac{{u_{\rm rms}}}{\eta {k_{\rm f}}},$$ where ${k_{\rm f}}=2\pi/\Delta r\approx21 R^{-1}$ is an estimate of the wavenumber of the largest eddies, and $\Delta r=R-r_0=0.3\,R$ is the thickness of the layer. The Coriolis number is defined as $${{\rm Co}}=\frac{2\Omega_0}{{u_{\rm rms}}{k_{\rm f}}},$$ where ${u_{\rm 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$, $\theta$, $\phi$, and a time interval during which the run is thermally relaxed and which covers several magnetic diffusion times. The averaging procedures employ the correct volume or surface elements of spherical polar coordinates. Note that for ${u_{\rm rms}}$ we omit the contribution from the azimuthal velocity, because its value is dominated by effects from the differential rotation [@KMGBC11] and compensate for this with the $3/2$ factor. The Taylor number can also be written as ${{\rm Ta}}={{\rm Co}}^2{{\rm Re}}^2({k_{\rm f}}R)^4$. Due to the fact that the initial stratification is isentropic, we quote the turbulent Rayleigh number ${{\rm Ra}_{\rm t}}$ from the thermally relaxed state of the run, $$\begin{aligned}
{{\rm Ra}_{\rm t}}\!=\!\frac{GM(\Delta r)^4}{\nu {\overline{\chi}_{\rm SGS}}R^2} \bigg(-\frac{1}{c_{\rm P}}\frac{{\rm d}{\langle s \rangle}_{\theta \phi t}}{{\rm d}r} \bigg)_{r_{\rm m}}.
\label{equ:Co}\end{aligned}$$ We also quote the value of $k_\omega=\omega_{\rm rms}/{u_{\rm rms}}$, where $\bm\omega=\bm\nabla \times {\bm u}$, and $\omega_{\rm rms}$ is the volume averaged rms value of $\bm\omega$. The magnetic field is expressed in equipartition field strengths, ${B_{\rm eq}}(r)=\langle \mu_0 \rho \bm{u}^2 \rangle^{1/2}_{\theta\phi t}$, where all three components of $\bm{u}$ are included. We define mean quantities as averages over the $\phi$-coordinate and denote them by overbars. However, as we will see, there can also be significant power in non-axisymmetric spherical harmonic modes with low azimuthal degree $m=1$ and 2, which will be discussed at the end of the paper.

The simulations were performed with the [Pencil Code]{}[^2], which uses a high-order finite difference method for solving the compressible equations of magnetohydrodynamics.

[ccccccccccccccccc]{}

A1 & $128\times256\times128$ & 71 & $1.5$ & $1.0$ & $1.0$ & 0.29 & 2.0 & 2.1 & $1.7$ & $3.8$ & $0.92$ & $0.83$ & $26$ & $26$ & $8.6$\
A2 & $128\times256\times128$ & 71 & $1.5$ & $1.0$ & $1.8$ & 0.29 & 2.0 & 2.1 & $1.7$ & $3.8$ & $0.92$ & $0.11$ & $24$ & $24$ & $12.8$\
B1 & $128\times256\times128$ & 82 & $2.5$ & $1.0$ & $0.64$ & 0.09 & 5.0 & 5.3 & $2.9$ & $3.8$ & $10.9$ & $1.1$ & $22$ & $22$ & $8.1$\
B2 & $128\times256\times128$ & 82 & $2.5$ & $1.0$ & $1.4$ & 0.09 & 5.0 & 5.2 & $2.9$ & $3.8$ & $10.9$ & $1.1$ & $20$ & $20$ & $13.7$\
C1 & $128\times256\times128$ & 56 & $2.5$ & $1.0$ & $1.4$ & 0.02 & 30 & 22 & $2.9$ & $3.8$ & $1.4\cdot10^{3}$ & $2.1$ & $35$ & $35$ & $7.8$\
C2 & $128\times256\times128$ & 56 & $2.5$ & $1.0$ & $4.0$ & 0.02 & 30 & 21 & $2.9$ & $3.8$ & $1.4\cdot10^{3}$ & $2.7$ & $31$ & $31$ & $14.8$\
D1 & $128\times256\times128$ & 503 & $7.5$ & $3.0$ & $0.16$ & 0.008 & 100 & 85 & $4.7$ & $0.63$ & $3.9\cdot10^{4}$ & $1.2$ & $11$ & $34$ & $8.0$\
D2 & $256\times512\times256$ & 269 & $4.0$ & $2.0$ & $1.0$ & 0.008 & 100 & 74 & $2.5$ & $0.63$ & $3.9\cdot10^{4}$ & $2.4$ & $25$ & $50$ & $9.1$\
E1 & $128\times256\times64$ & 56 & $2.5$ & $1.0$ & $1.4$ & 0.02 & 30 & 22 & $2.9$ & $3.8$ & $1.4\cdot10^{3}$ & $2.1$ & $34$ & $34$ & $7.9$\
E2 & $128\times256\times128$ & 56 & $2.5$ & $1.0$ & $1.4$ & 0.02 & 30 & 22 & $2.9$ & $3.8$ & $1.4\cdot10^{3}$ & $2.1$ & $35$ & $35$ & $7.8$\
E3 & $128\times256\times256$ & 56 & $2.5$ & $1.0$ & $1.4$ & 0.02 & 30 & 22 & $2.9$ & $3.8$ & $1.4\cdot10^{3}$ & $2.4$ & $35$ & $35$ & $7.9$\
E4 & $128\times256\times512$ & 67 & $3.0$ & $1.0$ & $1.0$ & 0.02 & 30 & 23 & $3.5$ & $3.8$ & $1.4\cdot10^{3}$ & $2.2$ & $28$ & $28$ & $8.1$ \[tab:runs\]

[cccccccccccc]{} A1 & 0.084 & 0.010 & 0.000 & 0.580 & 0.418 & 0.045 & 0.396 & 0.013 & 0.089 & 62\
A2 & 0.095 & 0.009 & 0.000 & 0.490 & 0.553 & 0.068 & 0.338 & 0.009 & 0.050 & 62\
B1 & 0.028 & 0.013 & 0.000 & 0.705 & 0.345 & 0.038 & 0.487 & 0.034 & 0.142 & 68\
B2 & 0.098 & 0.012 & 0.000 & 0.757 & 0.222 & 0.056 & 0.427 & 0.023 & 0.072 & 72\
C1 & 0.006 & 0.021 & 0.001 & 0.440 & 0.346 & 0.138 & 0.203 & 0.047 & 0.068 & 93\
C2 & 0.105 & 0.019 & 0.001 & 0.326 & 0.706 & 0.198 & 0.238 & 0.016 & 0.030 & 94\
D1 & 0.003 & 0.011 & 0.002 & 0.222 & 0.472 & 0.166 & 0.135 & 0.011 &-0.000 & 89\
D2 & 0.003 & 0.013 & 0.000 & 0.617 & 0.222 & 0.133 & 0.190 & 0.045 & 0.058 & 116\
E1 & 0.007 & 0.021 & 0.001 & 0.478 & 0.393 & 0.133 & 0.328 & 0.048 & 0.069 & 92\
E2 & 0.006 & 0.021 & 0.001 & 0.440 & 0.346 & 0.138 & 0.203 & 0.047 & 0.068 & 93\
E3 & 0.005 & 0.021 & 0.001 & 0.375 & 0.380 & 0.120 & 0.172 & 0.037 & 0.055 & 92\
E4 & 0.024 & 0.020 & 0.001 & 0.410 & 0.477 & 0.016 & 0.080 & 0.028 & 0.054 & 89 \[tab:runs2\]

Relation to reality {#sec:reality}
-------------------

In simulations, the maximum possible Rayleigh number is much smaller than in real stars due to the higher diffusivities. This implies higher energy fluxes and thus larger Mach numbers [@BCNS05]. To have realistic Coriolis numbers, the angular velocity in the Coriolis force has to be increased in proportion to one third power of the increase of the energy flux, but the centrifugal acceleration is omitted, as it would otherwise be unrealistically large [cf. @KMGBC11]. In the present models this would mean that the centrifugal acceleration is of the same order of magnitude as gravity, thus significantly altering the hydrostatic balance.

We note that we intend to use low values of $\mathcal{L}$ so that the Mach number is sufficiently below unity. This is particularly important when the stratification is strong. In our current formulation the unresolved turbulent heat conductivity, ${\chi_{\rm SGS}}$, acts on the total entropy and thus contributes to the radial heat flux. In the current models with ${{\rm Pr}}_{\rm SGS}$ greater than unity, the SGS-flux accounts for a few per cent of the total flux within the convection zone. Using smaller values of ${{\rm Pr}}_{\rm SGS}$ at the same Reynolds number would lead to a greater contribution due to the SGS-flux. To minimize the effects of the SGS-flux within the convection zone, we use the smallest possible value of ${\chi_{\rm SGS}}$ that is still compatible with numerical stability.

The span of time scales in our model is strongly compressed so as to comprise the full range all the way to the viscous, thermal, and resistive time scales. In the following we define acoustic, convective, thermal, resistive, and viscous time scales as follows: $$\tau_{\rm ac}=\sqrt{R^3/GM},\quad
\tau_{\rm conv}=H_{\rm P0}/{u_{\rm rms}}^{\rm(ref)},$$ $$\tau_{\rm th}=H_{\rm P0}^2/\chi_0,\quad
\tau_{\rm res}=H_{\rm P0}^2/\eta,\quad
\tau_{\rm visc}=H_{\rm P0}^2/\nu,$$ where $H_{\rm P0}$ is the pressure scale height at $r_0$, ${u_{\rm rms}}^{\rm(ref)}=(F_0/\rho_0)^{1/3}$ is a convective reference velocity based on the luminosity of the model, $L=4\pi r_0^2 F_0$, and $F_0$ is the total flux at $r_0$.

A visual comparison of these different time scales for the Sun and Run C1 is given in Figure \[fig:ptime\_scales\]. In order to allow for slow thermal and resistive relaxation processes, we require that their respective time scales are shorter than the run time $T$ of the simulation. As stated in Section \[sec:dimless\], the acoustic time scale of our model is equal to that of the Sun. This implies that all the other time scales must be significantly reduced: $\tau_{\rm conv}$ by a factor 70–100, $\tau_{\rm th}$ and $\tau_{\rm res}$ by a factor $10^7$, and $\tau_{\rm visc}$ by a factor $10^{14}$. This is accomplished by taking values of ${\cal L}$ that are not as small as in the Sun (where ${\cal L}\approx5\times10^{-11}$), but typically $3.8\times10^{-5}$ for Run C1. This just corresponds to taking values of the Rayleigh number that are on the order of $10^6$ rather than solar values (in excess of $10^{24}$). Likewise, shorter thermal, resistive, and viscous time scales are obtained by choosing values of the magnetic and fluid Reynolds number that are not as large as in the Sun and by choosing magnetic and fluid Prandtl numbers that are not as small as in the Sun.

![ Visual comparison of acoustic, convective, thermal, resistive, and viscous time scales both in the Sun (upper part, in red) and in Run C1 (lower part). In our models, resistive and viscous time scales are often equal, and the thermal time scale is 1.5–7.5 times longer; see Columns 4 and 5 of Table \[tab:runs\]. The simulation time scales are confined between the length of the time step, $\delta t$, and the maximum run time, $T$. []{data-label="fig:ptime_scales"}](ptime_scales.ps){width="46.00000%"}

For the purpose of comparing dynamo time scales of the model with the Sun, it is useful to rescale them such that $\tau_{\rm conv}$ coincides with that of the Sun. We can then compare the rotation rates of our models in Table \[tab:runs\] with that of the Sun: Runs A1 and A2 are 2 and 3 times solar, B1 and B2 are 3 and 4.4 times solar, C1 (including all of Set E) and C2 are 4.4 and 7 times solar, and D1 and D2 are 4.3 and 6 times solar.

Results {#sec:results}
=======

We perform runs for four values of $\xi$, corresponding to initial density contrasts $\Gamma_\rho^{(0)}=2, 5, 30$, and $100$. These runs are referred to as Sets A–D. In Set E we use $\Gamma_\rho^{(0)}=30$ and vary $\phi_0$ with all other parameters being kept the same as in Run C1, except in Run E4 where we use 20% higher viscosity and magnetic diffusivity than in the other runs in Set E. For each series, we consider different values of ${{\rm Ta}}$ and, as a consequence, of Co and Re. The hydrodynamic progenitors of the Runs B1, C1, and D1 correspond to Runs A4, B4, and C4, respectively, from [@KMB11]. The rest of the simulations were run from the initial conditions described in Section \[sec:initcond\].

Earlier studies applying fully spherical simulations have shown that organized large-scale magnetic fields appear provided the rotation of the star is rapid enough [@BBBMT10] and that at even higher rotation rates, cyclic solutions with poleward migration of the activity belts are obtained [@BMBBT11]. A similar transition has been observed in the spherical wedge models of [@KKBMT10] and [@KMB12a]. However, in the former case the oscillatory mode showed poleward migration, whereas in the latter an equatorward branch appears near the equator. Furthermore, in these runs the dynamo mode changes from one showing a high frequency cycle with poleward migration near the equator to another mode with lower frequency and equatorward migration when the magnetic field becomes dynamically important.

There are several differences between the models of [@KKBMT10] and [@KMB12a]: the amount of density stratification (a density contrast of 3 in comparison to 30), the efficiency of convective energy transport (20% versus close to 100% in the majority of the domain achieved by the use of ${\chi_{\rm SGS}}$; see also Figure \[fig:pflu\]), and the top boundary condition for entropy (constant temperature versus black body radiation). Here we concentrate on studying the influence of the density stratification on models similar to those presented in [@KMB12a].

![Luminosity of the energy fluxes from Run E4: radiative conduction (thin solid line), enthalpy (dashed), kinetic energy (dot-dashed), and unresolved subgrid scale (dotted) fluxes. The thick solid line is the sum of all contributions. The two dashed red lines indicate the zero and unity lines. []{data-label="fig:pflu"}](pflu_E4.eps){width="46.00000%"}

Thermal boundary effects and energy balance
-------------------------------------------

In [@KMB11] we started to apply the black–body boundary condition, Equation (\[eq:bbb\]), that has previously been used in mean-field models with thermodynamics [@R89; @BMT92; @KM00]. Instead of using the physical value for the Stefan–Boltzmann constant, we estimate the value of $\sigma$ so that the flux at the upper boundary is approximately that needed to transport the total luminosity of the star through the surface; see Table \[tab:runs\]. However, the final thermally relaxed state of the simulation can significantly deviate from the initial state. In combination with the nonlinearity of Equation (\[eq:bbb\]), the final stratification is usually somewhat different from the initial one; see Figure \[fig:pstrat\] for an illustrative example from Run C1. The final density stratification in this case is around 22, down from 30 in the initial state.

The main advantage of the black–body condition is that it allows the temperature at the surface more freedom than in our previous models where a constant temperature was imposed [@KKBMT10; @KMGBC11]. In particular, as the temperature is allowed to vary at the surface, this can be used as a diagnostic for possible irradiance variations. These issues are discussed further in Section \[sec:irra\].

Considering the energy balance, we show the averaged radial energy fluxes for Run E4 in Figure \[fig:pflu\]. We find that the simulation is thermally relaxed and that the total luminosity is close to the input luminosity, i.e., $L_{\rm tot}-L_0\approx0$. The fluxes are defined as: $$\begin{aligned}
\mathcal{F}_{\rm rad} &=& -K  \langle \nabla_r T \rangle,\\
\mathcal{F}_{\rm conv} &=& c_{\rm P}\langle(\rho u_r)^{\prime}T^{\prime} \rangle,\\
\mathcal{F}_{\rm kin} &=& \onehalf \left\langle\rho u_r{\bm u}^2\right\rangle,\\
\mathcal{F}_{\rm visc} &=& -2\nu \left\langle \rho u_i S_{ir} \right\rangle,\\
\mathcal{F}_{\rm turb} &=& -{\chi_{\rm SGS}}\langle\rho T
\nabla_r s\rangle,\\
\mathcal{F}_{\rm Poyn} &=& \langle E_\theta B_\phi - E_\phi B_\theta
\rangle/\mu_0,\end{aligned}$$ where ${\bm E}=\eta \mu_0 {\bm J}-{\bm u}\times{\bm B}$, the primes denote fluctuations, and angle brackets abbreviate $\langle\cdot\rangle_{\theta,\phi,t}$. The radiative flux carries energy into the convection zone and drops steeply as a function of radius so that it contributes only a few per cent in the middle of the convection zone. The resolved convection is responsible for transporting the energy through the majority of the layer, whereas the unresolved turbulent transport carries energy through the outer surface. The viscous and Poynting fluxes are much smaller, and are thus omitted in this figure. The flux of kinetic energy is also very small in the rapid rotation regime considered here [see also @ABBMT12].

![Initial (solid lines) and saturated (dashed) radial profiles of temperature $T$, density $\rho$, and pressure $p$, normalized by their respective values at the bottom of the domain (indicated by the subscript zero) from Run C1. The inset shows the specific entropy $s/c_{\rm P}$ from the same run.[]{data-label="fig:pstrat"}](pstrat_C1.eps){width="46.00000%"}

Dynamo excitation and large-scale magnetic fields
-------------------------------------------------

The azimuthally averaged toroidal magnetic fields from Sets A–D listed in Tables \[tab:runs\] and \[tab:runs2\] are shown in Figures \[fig:butterflyA\]–\[fig:butterflyD\]. The full time evolution from the introduction of the seed magnetic field to the final saturated state is shown for each run. Note that the magnitude of the seed field in terms of the equipartition strength is different in each set so direct comparisons between different sets are not possible. We measure the average growth rates during the kinematic stage, $$\lambda=\langle d\ln B_{\rm rms}/dt \rangle_t,$$ and find that $\lambda$ is greater for smaller stratification; see Column 2 of Table \[tab:runs2\] for $\tilde\lambda=\lambda/({u_{\rm rms}}{k_{\rm f}})$. Comparing Runs A1, B1, C1, and D2 with roughly comparable Reynolds and Coriolis numbers shows that the normalized growth rate decreases monotonically from 0.084 in Run A1 to just 0.003 in Run D2. Another striking feature is that $\tilde\lambda$ increases by a factor of nearly 20 from Run C1 to C2 whose only difference is that the latter has a roughly two times higher Coriolis number. It turns out that in all of the cases (Runs A1, A2, B1, B2, C2, and E4) with the highest growth rates, a dynamo mode with poleward migration at low latitudes, is excited first. In some of the runs this mode is later overcome by another one that can be quasi-stationary (Runs A1 and B1) or oscillatory with equatorward migration and a much longer cycle period (Runs C2 and E4).

Table \[tab:runs2\] shows that, even though the growth rates decrease dramatically with increasing stratification, many properties of the saturated stages are similar. In particular, the ratio of magnetic to kinetic energies does not seem to systematically depend on stratification, but rather on the Coriolis number, which varies only little between different runs.

In Figure \[fig:butterflyA\] we show the azimuthally averaged toroidal magnetic field ${\overline{B}}_\phi$ near the surface of the computational domain ($r=0.98\,R$) for two runs (A1 and A2; see Table \[tab:runs\]) with $\Gamma_\rho^{(0)}=2$. We find that in Run A1 with ${{\rm Co}}\approx8.7$ the mean magnetic field is initially oscillatory with poleward propagation of the activity belts. At $t{u_{\rm rms}}{k_{\rm f}}\approx400$ the dynamo mode changes to a quasi-steady configuration. In Run A2 a poleward mode persists throughout the simulation, although the oscillation period is irregular and significant hemispherical asymmetry exists. This behavior is similar to Run A4 presented in [@KKBMT10] with comparable stratification ($\Gamma\approx 3$) and Reynolds ($\approx 20$) numbers, but a somewhat lower Coriolis number[^3] ($\approx 4.7$). The transition to oscillatory solutions thus occurs at a lower ${{\rm Co}}$ in the models of [@KKBMT10]. A possible explanation is that in the present models we lack a lower overshoot layer which could affect the dominant dynamo mode.

In Set B with $\Gamma_\rho^{(0)}=5$ the situation is similar: in Run B1 with ${{\rm Co}}\approx8.1$ there is a poleward mode near the equator with a short cycle period which is visible from early times; see Figure \[fig:butterflyB\]. However, after around $t {u_{\rm rms}}{k_{\rm f}}=1200$ there is a dominating non-oscillatory mode that is especially clear at high latitudes. There are still hints of the poleward mode near the equator. In Run B2 with ${{\rm Co}}\approx13.7$, however, the poleward mode also prevails at late times. As in Run A2, the cycles show significant variability and hemispheric asymmetry. The runs in Sets A and B also show signs of non-axisymmetric ‘nests’ of convection [cf. @Bu02; @BBBMT08] in the hydrodynamical and kinematic stages. Once the magnetic field becomes dynamically important, these modes either vanish or they are significantly damped.

Increasing the stratification further to $\Gamma_\rho^{(0)}=30$ (Set C) the dynamo solutions at lower rotation rates, ${{\rm Co}}\lesssim 5$, are still quasi-steady; see Figure 2 of [@KMB12a]. However, a watershed regarding the oscillatory modes at higher ${{\rm Co}}$ seems to have been reached so that the irregular poleward migration seen in Sets A and B is replaced by more regular equatorward patterns. In Run C1 with ${{\rm Co}}\approx8.7$ the poleward migration near the equator is also visible in the kinematic stage where the equatorward mode is not yet excited; see Figure \[fig:butterflyC\]. The poleward mode near the equator is more prominent in the early stages of Run C2 with ${{\rm Co}}\approx14.7$, but subdominant at late times.

![${\overline{B}}_\phi$ near the surface of the star at $r=0.98\,R$ as a function of latitude ($=90\degr-\theta$) and time for Runs A1 (top) and A2 (bottom). The white dotted line denotes the equator $90\degr-\theta=0$. []{data-label="fig:butterflyA"}](m128x256x128d2_Bphi_small.ps "fig:"){width="48.00000%"} ![${\overline{B}}_\phi$ near the surface of the star at $r=0.98\,R$ as a function of latitude ($=90\degr-\theta$) and time for Runs A1 (top) and A2 (bottom). The white dotted line denotes the equator $90\degr-\theta=0$. []{data-label="fig:butterflyA"}](m128x256x128d1_Bphi_small.ps "fig:"){width="48.00000%"}

![Same as Figure \[fig:butterflyA\], but for Runs B1 (top) and B2 (bottom).[]{data-label="fig:butterflyB"}](m128x256x128b4_Bphi_small.ps "fig:"){width="48.00000%"} ![Same as Figure \[fig:butterflyA\], but for Runs B1 (top) and B2 (bottom).[]{data-label="fig:butterflyB"}](m128x256x128b5_Bphi_small.ps "fig:"){width="48.00000%"}

![Same as Figure \[fig:butterflyA\], but for Runs C1 (top) and C2 (bottom). Note the difference in cycle frequency between the early times when the frequency is similar to that of Run B2 (Figure \[fig:butterflyB\]) and late times. []{data-label="fig:butterflyC"}](m128x256x128a4_Bphi_small.ps "fig:"){width="48.00000%"} ![Same as Figure \[fig:butterflyA\], but for Runs C1 (top) and C2 (bottom). Note the difference in cycle frequency between the early times when the frequency is similar to that of Run B2 (Figure \[fig:butterflyB\]) and late times. []{data-label="fig:butterflyC"}](m128x256x128a5_Bphi_small.ps "fig:"){width="48.00000%"}

![Same as Figure \[fig:butterflyA\], but for Runs D1 (top) and D2 (bottom).[]{data-label="fig:butterflyD"}](m128x256x128c4_Bphi_small.ps "fig:"){width="48.00000%"} ![Same as Figure \[fig:butterflyA\], but for Runs D1 (top) and D2 (bottom).[]{data-label="fig:butterflyD"}](m256x512x256c5_Bphi_small.ps "fig:"){width="48.00000%"}

For $\Gamma_\rho^{(0)}=100$ (Set D) the general picture is similar to that in Set C. Quasi-steady configurations at lower rotation rates change into equatorward migrating solutions at sufficiently high values of ${{\rm Co}}$. We find that this transition occurs between ${{\rm Co}}=5$ and $8$, similar to Set C; see Fig. 2 of [@KMB12a]. For Set D the equatorward mode is visible for both of its runs; see Figure \[fig:butterflyD\]. In Run D1 no poleward migration at low latitudes is seen in the kinematic stage. Also, the poleward migrating branch at high latitudes is missing in the non-linear stage. Both of these features are present in Run D2. The apparently slower growth of the magnetic field in Run D1 is due to a two orders of magnitude lower seed magnetic field than in Run D2.

Diagnostic stellar activity diagrams
------------------------------------

To identify the possibility of different types of dynamos, it is useful to classify them in diagrams relating their characteristic properties. In the geodynamo literature it has become customary to consider the Elsasser number as a measure of the magnetic energy. It correlates well with ${{\rm Rm}}$ [@CA06], but this is partially explained by the fact that ${{\rm Rm}}$ itself enters in the definition of the Elsasser number. Geodynamo models are mostly dominated by a strong dipolar component. [@GDW12] have shown that such solutions fall on a branch that is distinct from the cyclic solutions studied here, and that the latter solutions become favored once density stratification is large and rotation is sufficiently rapid so that large-scale non-axisymmetric fields become dominant [see also @NBBMT13]. However, this type of analysis is not well suited for the present work, where ${{\rm Rm}}$ and ${{\rm Co}}$ vary only little. Furthermore, these tools do not characterize the nature of magnetic cycles, which is the focus of this section.

To connect our results with observations of magnetically active stars we compute the ratio of cycle to rotation frequency $\omega_{\rm cyc}/\Omega_0$, where $\omega_{\rm cyc}=2\pi/T_{\rm cyc}$ is the cycle frequency of magnetic energy of the mean field and $T_{\rm cyc}$ its period. Plotting this ratio as a function of the Coriolis number for stars exhibiting chromospheric activity has shown that stars tend to group along inactive and active branches [@BST98], and for higher Coriolis numbers along a super-active branch [@SB99]. Six of our simulations (Runs A2, B2, C1, C2, D1, and D2), excluding the runs in Set E (which are very similar to each other and to Run C1), show cycles and can thus be used in this analysis. We compute the cycle frequency from the highest peak of a temporal Fourier transformation of the time series for ${\overline{B}}_\phi$ averaged over a latitudinal strip of $\pm10\degr\ldots30\degr$ near the surface. The results are shown in Figure \[fig:pocycom\](a).

Three of the models, Runs C1, D1, and D2, fall on a branch labeled ‘${\cal A}$?’ for active stars, while Run C2 might be suggestive of the superactive stars of [@SB99], labeled here ‘${\cal S}$??’. Runs A2 and B2 show irregular cycles and group along the branch labeled ‘${\cal I}$?’ for inactive stars. The question marks on these labels in Figure \[fig:pocycom\](a) indicate that the association with real branches is quite uncertain and somewhat premature, because there are too few models. We cannot be sure that there are no models connecting the group of Runs C1, D1, D2 with that of A2 and B2 through a single line with a steeper slope. Nevertheless, this plot allows us to see that, while the separation in the ratio $\omega_{\rm cyc}/\Omega_0$ is slightly less for the two groups of runs compared with active and inactive stars, their relative ordering in the value of ${{\rm Co}}$ is actually the other way around. One would therefore not have referred to Runs A2 and B2 as inactive just because their $\omega_{\rm cyc}/\Omega_0$ ratio agrees with that of inactive stars. In fact, their $E_{\rm mag}/E_{\rm kin}$ ratios (a measure of stellar activity) in Table \[tab:runs2\] are typically larger than for Runs C1, D1, and D2.

As is visible from Figure \[fig:pocycom\](c), there is no clear relation between ${{\rm Co}}$ and $E_{\rm mag}/E_{\rm kin}$, which is different from stars for which there is a clear relation between ${{\rm Co}}$ (referred to as the inverse Rossby number in that context) and stellar activity; see [@BST98] for details and references. Furthermore, there are also no indications of branches in the graph of $\omega_{\rm cyc}/\Omega_0$ versus $E_{\rm mag}/E_{\rm kin}$; see Figure \[fig:pocycom\](b). Instead, there might just be one group in it, possibly with a positive correlation, i.e., $\omega_{\rm cyc}/\Omega_0$ might increase with $E_{\rm mag}/E_{\rm kin}$. Such a possibility does indeed arise when considering the frequency ratio versus the dimensional rotation rate [@OKS00]. However, as discussed by [@BST98], a positive slope is not easily explained in the framework of standard mean-field dynamo theory, where the frequency ratio is usually a decreasing function of normalized rotation rate and activity parameter [@To98; @SB99].

In conclusion, we reiterate that the quantity $\omega_{\rm cyc}/\Omega_0$ is an important and robust property of cyclic dynamo models and its dependence on other properties of the model should therefore be a useful characteristics that can be compared with other models and ultimately with actual stars. Here we have made a first attempt in classifying model results in this way.

![Diagnostic diagrams from six runs that show cyclic activity. (a) Ratio of cycle and rotation frequencies vs. $\log{{\rm Co}}$. The dotted and dashed lines are given by $c_i {{\rm Co}}^{\sigma_i}$, where $\sigma_i$ correspond to those in [@BST98] for active (labeled ‘${\cal A}$?’) and inactive (‘${\cal I}$?’) stars, while $c_i$ are used as fit parameters. The label ‘${\cal S}$??’ indicates the possibility of the superactive branch in [@SB99]. (b) Ratio of cycle and rotation frequencies vs. magnetic to kinetic energy $E_{\rm mag}/E_{\rm kin}$. (c) Time averaged $E_{\rm mag}/E_{\rm kin}$ vs. $\log{{\rm Co}}$.[]{data-label="fig:pocycom"}](pocycom.eps "fig:"){width="46.00000%"} ![Diagnostic diagrams from six runs that show cyclic activity. (a) Ratio of cycle and rotation frequencies vs. $\log{{\rm Co}}$. The dotted and dashed lines are given by $c_i {{\rm Co}}^{\sigma_i}$, where $\sigma_i$ correspond to those in [@BST98] for active (labeled ‘${\cal A}$?’) and inactive (‘${\cal I}$?’) stars, while $c_i$ are used as fit parameters. The label ‘${\cal S}$??’ indicates the possibility of the superactive branch in [@SB99]. (b) Ratio of cycle and rotation frequencies vs. magnetic to kinetic energy $E_{\rm mag}/E_{\rm kin}$. (c) Time averaged $E_{\rm mag}/E_{\rm kin}$ vs. $\log{{\rm Co}}$.[]{data-label="fig:pocycom"}](pocycom_E.eps "fig:"){width="46.00000%"} ![Diagnostic diagrams from six runs that show cyclic activity. (a) Ratio of cycle and rotation frequencies vs. $\log{{\rm Co}}$. The dotted and dashed lines are given by $c_i {{\rm Co}}^{\sigma_i}$, where $\sigma_i$ correspond to those in [@BST98] for active (labeled ‘${\cal A}$?’) and inactive (‘${\cal I}$?’) stars, while $c_i$ are used as fit parameters. The label ‘${\cal S}$??’ indicates the possibility of the superactive branch in [@SB99]. (b) Ratio of cycle and rotation frequencies vs. magnetic to kinetic energy $E_{\rm mag}/E_{\rm kin}$. (c) Time averaged $E_{\rm mag}/E_{\rm kin}$ vs. $\log{{\rm Co}}$.[]{data-label="fig:pocycom"}](pocycom_Co.eps "fig:"){width="46.00000%"}

![Time averaged mean rotation profiles ${\overline{\Omega}}/\Omega_0$ (gray/color scale and line contours) from Sets A, B, C, and D.[]{data-label="fig:pOm"}](pOm2_small.ps){width="0.9\columnwidth"}

Differential rotation and meridional circulation
------------------------------------------------

![ Meridional circulation in the northern hemisphere of the convection zone of Run C1 (left panel) and Run D1 (right) shown as vectors of the mass flux ${\overline{\rho}}({\overline{u}}_r,{\overline{u}}_\theta,0)$, which is also averaged over a time span of around 250 turnover times in the saturated state. The black solid lines indicate the surface ($r=R$) and the bottom of the convection zone ($r=0.7\,R$), and the red solid line indicates the position of the inner tangent cylinder. Note, that for Run D1 (right), the mass flux have been multiplied by a factor of 5 to emphasize the structure. []{data-label="fig:meri"}](merion_conv.ps "fig:"){width="0.49\columnwidth"} ![ Meridional circulation in the northern hemisphere of the convection zone of Run C1 (left panel) and Run D1 (right) shown as vectors of the mass flux ${\overline{\rho}}({\overline{u}}_r,{\overline{u}}_\theta,0)$, which is also averaged over a time span of around 250 turnover times in the saturated state. The black solid lines indicate the surface ($r=R$) and the bottom of the convection zone ($r=0.7\,R$), and the red solid line indicates the position of the inner tangent cylinder. Note, that for Run D1 (right), the mass flux have been multiplied by a factor of 5 to emphasize the structure. []{data-label="fig:meri"}](merion_conv_II.ps "fig:"){width="0.49\columnwidth"}

Non-uniform rotation of the convection zone of the Sun is an important ingredient in maintaining the large-scale magnetic field. Furthermore, the sign of the radial gradient of the mean angular velocity plays a crucial role in deciding whether the dynamo wave propagates toward the pole or the equator in $\alpha$–$\Omega$ mean-field models [e.g., @Pa55b; @Pa87a]. In the following we use the local angular velocity defined as ${\overline{\Omega}}=\Omega_0+{\overline{u}}_\phi/r\sin\theta$. Azimuthally averaged rotation profiles from the runs in Sets A to D are shown in Figure \[fig:pOm\]. The rotation profiles of Runs E1, E3, and E4 are very similar to that of Run C1. We quantify the radial and latitudinal differential rotation by $$\begin{aligned}
\Delta_\Omega^{(r)} &=& \frac{{\overline{\Omega}}_{\rm eq}-{\overline{\Omega}}_{\rm bot}}{{\overline{\Omega}}_{\rm eq}},\\
\Delta_\Omega^{(\theta)} &=& \frac{{\overline{\Omega}}_{\rm eq}-{\overline{\Omega}}_{\rm pole}}{{\overline{\Omega}}_{\rm eq}},\end{aligned}$$ where ${\overline{\Omega}}_{\rm eq}={\overline{\Omega}}(R,\pi/2)$ and ${\overline{\Omega}}_{\rm
  bot}={\overline{\Omega}}(r_0,\pi/2)$ are the angular velocities at the top and bottom at the equator, respectively, and ${\overline{\Omega}}_{\rm
pole}=[{\overline{\Omega}}(R,\theta_0) + {\overline{\Omega}}(R,\pi-\theta_0)]/2$. It has long been recognized that dynamo-generated magnetic fields can have an important effect on the angular velocity [@Gi83; @Gl85; @Gl87]. Indeed, magnetic fields affect the turbulence that gives rise to Reynolds stress and turbulent convective heat flux [e.g., @KPR94; @KKT04]. Furthermore, the large-scale flows are directly influenced by the Lorentz force when the magnetic field is strong enough [e.g., @MP75]. A magnetically caused decrease of $\Delta_\Omega^{(\theta)}$ has also been observed in LES models [e.g., @BMT04]. Comparing the latitudinal differential rotation in Run B1 with that of the otherwise identical hydrodynamic Run A4 of [@KMB11], we find that $\Delta_\Omega^{(\theta)}$ decreases only slightly from 0.15 to 0.14. For $\Delta_\Omega^{(r)}$ the change is more dramatic—from 0.079 to 0.034. The fraction of kinetic energy contained in the differential rotation, $E_{\rm
  rot}/E_{\rm kin}$, drops from 0.91 to 0.71. A similar decrease is observed in Run C1 in comparison to its hydrodynamical parent Run B4 of [@KMB11] with $\Delta_\Omega^{(\theta)}$ changing from 0.08 to 0.07, $\Delta_\Omega^{(r)}$ from 0.066 to 0.047, and $E_{\rm rot}/E_{\rm
  kin}$ dropping from 0.58 to 0.44. Similar changes have also been seen in dynamos from forced turbulence in Cartesian domains [@B01], in addition to those from convective turbulence in spherical shells [@BMT04].

In all cases in Figure \[fig:pOm\], we see a rapidly spinning equator with a positive radial gradient of ${\overline{\Omega}}$. The latitudinal variation of angular velocity is, however, not always monotonic and there can be local minima at mid-latitudes, as is seen, for example, in Run C1. Similar features have previously been seen [see, e.g., @METCGG00; @KMGBC11] and might be related to the lack of small-scale turbulence. Especially at larger stratification one would expect smaller-scale turbulent structures to emerge, but this means large Reynolds numbers and thus requires sufficient resolution, which is not currently possible.

The amount of latitudinal differential rotation (here 0–0.09; see Table \[tab:runs2\]) is clearly less than in the Sun where $\Delta_\Omega^{(\theta)}\approx0.2$ between the equator and latitude $60\degr$ [e.g., @Schouea98]. Furthermore, $\Delta_\Omega^{(\theta)}$ generally decreases within each set of runs as ${{\rm Co}}$ increases, except for Runs D1 and D2 where the value increases; see Table \[tab:runs2\]. However, in Run D1 the lower Reynolds number possibly contributes to the weak differential rotation in comparison to Run D2 with comparable ${{\rm Co}}$. The rotation profiles appear to be dominated by the Taylor-Proudman balance, except at very low latitudes where the baroclinic term is significant; see Figure 9 of [@WKMB13]. In this companion paper, we show that an outer coronal layer seems to favor a solar-like rotation, which shows even radially orientated contours of constant rotation. Such ‘spoke-like’ rotation profiles have thus far only been obtained in mean-field models involving anisotropic heat transport [e.g., @BMT92; @KR95] or a subadiabatic tachocline [@Re05], and in purely hydrodynamic LES models where a latitudinal entropy gradient is enforced at the lower boundary [@MBT06], or where a stably stratified layer is included below the convection zone [@BMT11].

The meridional circulation is weak in all cases and typically shows multiple cells in the radial direction. In Figure \[fig:meri\], we plot the mean mass flux, ${\overline{\rho}}({\overline{u}}_r,{\overline{u}}_\theta,0)$, of the meridional circulation for Runs C1 and D1. In Run C1 the circulation pattern is mostly concentrated in the equatorial region outside the inner tangent cylinder, where we find a solar-like anti-clockwise cell at low latitudes ($<30^\circ$) in the upper third of the convection zone. There are additional cells deeper down and also at higher latitudes. Only the cell near the surface seems to have the same curvature as the surface, while the others, in particular the strong one above the inner tangent cylinder, seem to be parallel to the rotation axis. This is similar to earlier results by [@KMB12a] where the meridional circulation pattern was shown in terms of the velocity. The circulation pattern in Run D1 is qualitatively quite similar, but the velocity is smaller by roughly a factor of five. Similar patterns of multi-cellular meridional circulation have also been seen in anelastic simulations using spherical harmonics [see, e.g., @NBBMT13] and in models with an outer coronal layer [@WKMB13]. In addition, as we will show in the next section, the importance of meridional circulation relative to the turbulent magnetic diffusivity is rather low, which is another reason why it cannot play an important role in our models.

![Local dynamo parameter $c_\alpha$ from Sets A, B, C, and D.[]{data-label="fig:pCalp_sat"}](pCalp_sat2_small.ps){width="0.9\columnwidth"}

![Local dynamo parameter $c_\Omega$ from Sets A, B, C, and D. We omit regions closer to $2.5\degr$ from the latitudinal boundaries.[]{data-label="fig:pCOm"}](pCOm2_small.ps){width="0.895\columnwidth"}

![Local dynamo parameter $c_U$ from Sets A, B, C, and D.[]{data-label="fig:pCU"}](pCU2_small.ps){width="0.9\columnwidth"}

Estimates of local dynamo parameters {#DynamoParameters}
------------------------------------

To estimate the dynamo parameters related to $\alpha$-effect, radial differential rotation, and meridional circulation, we consider local ($r$- and $\theta$-dependent) versions of dynamo numbers, referred to as local dynamo parameters that are defined by $$c_\alpha=\frac{\alpha \Delta r}{{\eta_{\rm t0}}},\quad
c_\Omega=\frac{{\partial}{\overline{\Omega}}/{\partial}r (\Delta r)^3}{{\eta_{\rm t0}}},\quad
c_U=\frac{{\overline{u}}_{\rm mer}^{\rm rms} \Delta r}{{\eta_{\rm t0}}},$$ where ${\partial}{\overline{\Omega}}/{\partial}r$ is the $r$- and $\theta$-dependent radial gradient of ${\overline{\Omega}}$, $\Delta r=R-r_0$ is the thickness of the layer, and $\alpha$ is a proxy of the $\alpha$-effect [@PFL76], $$\alpha = -\onethird \tau ({\overline{\mbox{\boldmath $\omega$}\cdot{\bm u}}}
-{\overline{{\bm j}\cdot{\bm b}}}/{\overline{\rho}}),$$ with $\tau=\alpha_{\rm MLT}H_{\rm P}/{u_{\rm rms}}(r,\theta)$ being the local convective turnover time and $\alpha_{\rm MLT}$ the mixing length parameter. We use $\alpha_{\rm MLT}=5/3$ in this work. We estimate the turbulent diffusivity by ${\eta_{\rm t0}}=\tau{u_{\rm rms}}^2(r,\theta)/3$. Furthermore, ${\overline{u}}_{\rm mer}^{\rm rms}=\sqrt{{\overline{u}}_r^2 +
  {\overline{u}}_\theta^2}$ is the rms value of the meridional circulation.

The results for the local dynamo parameters are shown in Figures \[fig:pCalp\_sat\]–\[fig:pCU\]. Generally, the values of $c_\alpha$ are fairly large, and those of $c_\Omega$ surprisingly small, suggesting that the dynamos might mainly be of $\alpha^2$ type. In the following, however, we focus on relative changes between different runs. It turns out that there is a weak tendency for $c_\alpha$ to increase as a function of $\Gamma_\rho$ (from Sets A to D) and ${{\rm Co}}$ (from subsets 1 to 2). In Set B, however, $c_\alpha$ decreases by a third from Run B1 to B2. The spatial distribution of $c_\alpha$ becomes more concentrated near the radial boundaries as $\Gamma_\rho$ increases.

We find that differential rotation is strongest near the equator in all cases. Sets A and B have extended regions outside the inner tangent cylinder and at low latitudes where $c_\Omega$ is large, but in all cases $c_\Omega$ is clearly smaller than $c_\alpha$. This is surprising given the fact that the energy of the mean toroidal field is greater than that of the mean poloidal field by a significant factor (see $E_{\rm pol}$ and $E_{\rm tor}$ in Table \[tab:runs2\]) which would be expected if differential rotation dominates over the $\alpha$-effect in maintaining the field. In Runs C1, C2, and D1, $c_\alpha$ and $c_\Omega$ have comparable magnitudes whereas in Run D2 the maximum of $c_\Omega$ is roughly twice that of $c_\alpha$. However, in these cases the toroidal and poloidal field energies are roughly comparable (see Table \[tab:runs2\]). For Set C (and especially for Run C2) there are broad regions where $c_\Omega$ is negative. In this connection we recall that in the diagnostic diagrams (Figure \[fig:pocycom\]), C2 appears as an outlier and far away from the ${\cal A?}$ and ${\cal I?}$ branches. Furthermore, in the more strongly stratified models, $c_\alpha$ shows enhanced values at low latitudes. However, for the most strongly stratified models this is only true of Run D2, which is rotating slightly faster than Run D1. This is interesting in view of the fact that many mean-field dynamos produce too strong fields at high latitudes, which is then ‘artificially’ reduced by an ad-hoc factor proportional to $\sin^2\theta$ [@RB95] or other such variants [@DdTGAW04; @PK11] for $\alpha$. We note that in local convection simulations, the $\alpha$-effect has been found to peak at mid-latitudes for rapid rotation [@KKOS06].

We find that $c_U$ is always small in comparison to both $c_\alpha$ and $c_\Omega$. Note, however, that the range of $c_U$ does increase as we go from Set A to Set D. Figure \[fig:pCU\] also shows the concentration of coherent meridional circulation cells in the equatorial regions with a multi-cell structure.

In flux transport dynamos, $c_U$ has values of several hundreds [@KRS01]. This is a consequence of choosing a small value of the turbulent magnetic diffusivity. In our simulations, on the other hand, $c_U$ is much smaller. This is a consequence of faster turbulent motions, making the turbulent diffusivity large and therefore $c_U$ small. Whether or not this also applies to more realistic models remains to be seen.

Phase relation and nature of the dynamo
---------------------------------------

The relative magnitudes of the estimated values of $c_\alpha$ and $c_\Omega$, and also the comparable amplitudes of ${\overline{B}}_r$ and ${\overline{B}}_\phi$, shown in Figure 4(a) of [@KMB12a], strongly suggest that the dynamos of this study are not of $\alpha^2\Omega$ type, as is usually expected for the Sun. This can be motivated further through direct inspection of the $\Omega$ term in the equation for the mean toroidal field. Following [@SPD12], we compare the $\Omega$-effect, $r\sin\theta {\overline{\bm B}}_{\rm pol}\cdot\bm\nabla {\overline{\Omega}}$, with the mean toroidal field. The results for Runs C1 and D1 are shown in Figure \[fig:BpgOm\], where we have scaled the $\Omega$ term by the magnetic cycle period, $T_{\rm cyc}$. A fraction of this would be responsible for the production of mean toroidal field for the next cycle. For Run C1, the magnitude of this term is actually large compared with ${\overline{B}}_\phi$, and the two are clearly correlated at latitudes below $\pm35\degr$, which is also where equatorward migration is seen. For Run D1, however, no clear correlation is seen even at low latitudes. The possibility of $\alpha^2\Omega$ type dynamo action therefore remains unclear, and especially for Run D1 it may not be the dominant mechanism. To explore the possibility that our dynamo is of $\alpha^2$ type, we now consider the phase relation between ${\overline{B}}_r$ and ${\overline{B}}_\phi$; see Figure \[fig:pphase\].

![$\Omega$-effect, as quantified by $T_{\rm cyc} r \sin \theta {\overline{\bm B}}_{\rm pol} \cdot
  \bm\nabla {\overline{\Omega}}$, where ${\overline{\bm B}}_{\rm
    pol}=({\overline{B}}_r,{\overline{B}}_\theta)$ (left panels), and the mean toroidal magnetic field ${\overline{B}}_\phi$ (right panels) normalized by ${B_{\rm eq}}$ from the saturated states of Runs C1 (upper panels) and D1 (lower panel). The data is averaged over the longitude and approximately 60 convective turnover times in both cases.[]{data-label="fig:BpgOm"}](m128x256x128a4_BpgOm_small.ps "fig:"){width="\columnwidth"} ![$\Omega$-effect, as quantified by $T_{\rm cyc} r \sin \theta {\overline{\bm B}}_{\rm pol} \cdot
  \bm\nabla {\overline{\Omega}}$, where ${\overline{\bm B}}_{\rm
    pol}=({\overline{B}}_r,{\overline{B}}_\theta)$ (left panels), and the mean toroidal magnetic field ${\overline{B}}_\phi$ (right panels) normalized by ${B_{\rm eq}}$ from the saturated states of Runs C1 (upper panels) and D1 (lower panel). The data is averaged over the longitude and approximately 60 convective turnover times in both cases.[]{data-label="fig:BpgOm"}](m128x256x128c4_BpgOm_small.ps "fig:"){width="\columnwidth"}

![ Phase relations of ${\overline{B}}_r$ (thick red lines) and ${\overline{B}}_\phi$ (black dashed lines) for Run C1 at (a) $70^\circ$ and (b) $30^\circ$ latitude, compared with results of mean-field dynamos of (c) $\alpha^2$ type and (d), (e) $\alpha^2\Omega$ type, with positive and negative shear, respectively. The amplitudes have been rescaled to unity. Note that only the $\alpha^2$ dynamo has approximately the phase relation seen in the simulations. []{data-label="fig:pphase"}](pphase.ps){width="\columnwidth"}

![image](pmoll_uu_E1_small.ps){height="5cm"}![image](pmoll_uu_E2_small.ps){height="5cm"}![image](pmoll_uu_E3_small.ps){height="5cm"}![image](pmoll_uu_E4_small.ps){height="5cm"} ![image](pmoll_bb_E1_small.ps){height="5cm"}![image](pmoll_bb_E2_small.ps){height="5cm"}![image](pmoll_bb_E3_small.ps){height="5cm"}![image](pmoll_bb_E4_small.ps){height="5cm"}

![Same as Figure \[fig:butterflyA\] but for Run E4.[]{data-label="fig:butterflyE4"}](m128x256x512a4_Bphi_small.ps){width="48.00000%"}

For $\alpha^2\Omega$ dynamos, the phase relation between ${\overline{B}}_r$ and ${\overline{B}}_\phi$ is commonly used to determine the sign of the radial differential rotation [@St76; @Yo76]. By contrast, the sign of $\alpha$ is determined by the sense of migration of the dynamo wave. For negative radial shear, ${\overline{B}}_r$ and ${\overline{B}}_\phi$ are approximately in antiphase with ${\overline{B}}_\phi$ preceding ${\overline{B}}_r$ by $\approx3\pi/4$. For positive radial shear, ${\overline{B}}_r$ and ${\overline{B}}_\phi$ are approximately in phase with ${\overline{B}}_\phi$ lagging ${\overline{B}}_r$ by $\approx\pi/4$. In our simulations, radial shear is indeed positive, but ${\overline{B}}_\phi$ precedes ${\overline{B}}_r$ by a certain amount; see Figures \[fig:pphase\](a) and (b). This cannot be explained by an $\alpha^2\Omega$ dynamo where (for positive radial shear) ${\overline{B}}_\phi$ [*lags*]{} ${\overline{B}}_r$ by $\pi/4$.

Another possibility are oscillatory $\alpha^2$ dynamos of the type recently found by [@MTKB10] using direct numerical simulations of forced turbulence in a spherical wedge. Those models have also been used to study the effects of an outer coronal layer to shed magnetic helicity [@WBM11]. Oscillatory $\alpha^2$ dynamos were first studied by [@BS87] and [@Ra87]; see also the monograph of [@RH04]. Such solutions have also been studied in connection with the geodynamo, where the $\alpha$-effect might change sign in the middle of the outer liquid iron core [@SG03]. By contrast, in the simulations of [@MTKB10] and [@WBM11], $\alpha$ changes sign about the equator. They used a perfect conductor boundary condition at high latitudes and found equatorward migrating dynamo waves. With a vacuum condition, on the other hand, mean-field simulations have predicted poleward migration [@BCC09]. Those simulations were done in Cartesian geometry, where $(x,y,z)$ can be identified with $(r,\phi,-\theta)$. Looking at their Figure 2, it is clear that ${\overline{B}}_y$ lags ${\overline{B}}_x$ by $\pi/2$.

We have verified the phase relations of the Cartesian model of [@BCC09] with a one-dimensional spherical model[^4], where $\alpha=\alpha_0\cos\theta$ has been assumed, which changes sign about the equator at $\theta=\pi/2$. The dynamo number for the marginally excited case is $\alpha_0 R/{\eta_{\rm t0}}\approx23.63$ and, as expected, ${\overline{B}}_\phi$ lags ${\overline{B}}_r$ by $\pi/2$; see Figure \[fig:pphase\](c). The amplitudes have been rescaled to unity. The corresponding behavior for an $\alpha^2\Omega$ dynamo is shown in Figures \[fig:pphase\](d) and (e), where ${\overline{B}}_\phi$ either precedes ${\overline{B}}_r$ by $\pi/4$ or lags ${\overline{B}}_r$ by $3\pi/4$. In this case, we have used a Cartesian model with constant $\alpha$, constant shear, $S=du_y/dx={\rm const}$, and periodic boundaries in a domain $0<z<L$, where the critical dynamo number is $\alpha S L^3/{\eta_{\rm t0}}^2\approx8\pi^2$. In this model, the Cartesian coordinates $(x,y,z)$ correspond to $(-r,\phi,\theta)$, so positive (negative) values of $S$ correspond to negative (positive) radial angular velocity gradients. Neither of the phase relations of these two models agrees with those of the DNS.

Another hint pointing toward an $\alpha^2$ dynamo in Run C1 is, that the magnetic field is particularly strong in the middle of the convection zone ($0.8\,R<r<0.9\,R$), from where dynamo waves seem to propagate toward the surface and the bottom of the convection zone; see Figure 3(a) of [@KMB12a]. Even though there exists no tachocline at the bottom or a near-surface shear layer at the top of our convection zone, the $\Omega$-effect appears to be larger toward the bottom and top of the convection zone; see Figure \[fig:pCOm\]. Therefore, an $\Omega$-effect would produce the magnetic field mainly at the bottom and the top of the convection zone, which is not the case in our simulation. The case of Run D1 is more clear because the $\Omega$-effect is weak except near the boundaries (Figure \[fig:pCOm\]) and the toroidal field shows no correlation with it; see Figure \[fig:BpgOm\]. We therefore suggest that oscillatory $\alpha^2$ dynamos of the type found by [@MTKB10] might explain the origin of equatorward migrating dynamo waves in the spherical wedge simulations of [@KMB12a]. It is also possible that this mechanism explains the poleward migration at high latitudes, but detailed comparisons must await a proper determination of $\alpha$-effect and turbulent diffusivity tensors. A first step toward this has recently been attempted by [@RCGBS11] who estimated the tensor components of $\alpha$ by correlating the electromotive force with the mean magnetic field using singular value decomposition. These results were applied in mean-field models of [@SCB13], in an effort to explain the dynamos seen by [@GCS10]. However, this analysis is flawed in the sense that the diffusive part of the electromotive force cannot be separated from the one related to the $\alpha$-effect. This has been shown to lead to erroneous estimates of $\alpha$ [@KKB10]. The only reliable way to compute the turbulent diffusion tensor is currently possible with the test-field method [@SRSRC05; @SRSRC07]. We postpone such analysis to a future publication.

Effect of domain size
---------------------

We recently reported equatorward migration of activity belts in a spherical wedge simulation [@KMB12a]. There we gave results from simulations with a $\phi$-extent of $\pi/2$. However, at large values of the Coriolis number, the $\alpha$-effect becomes sufficiently anisotropic and differential rotation weak so that non-axisymmetric solutions become possible; see [@MB95] for corresponding mean-field models with dominant $m=1$ modes in the limit of rapid rotation. To allow for such modes, we now choose a $\phi$-extent of up to $2\pi$ for the same model as in [@KMB12a]. In the present case, we find that for ${{\rm Co}}\approx7.8$ it is possible that non-axisymmetric dynamo modes of low azimuthal order ($m=1$ or 2) can be dominant. This was not possible in the simulations of [@KMB12a]. The same applies to non-axisymmetric modes excited in hydrodynamic convection [e.g., @Bu02; @BBBMT08; @KMGBC11; @ABBMT12].

We test the robustness of the equatorward migration by performing runs with $\phi_0=\pi/4$, $\pi/2$, $\pi$, and $2\pi$ with otherwise similar parameters; see Table \[tab:runs\]. We find that the same dynamo mode producing equatorward migration is ultimately excited in all of these runs. The only qualitatively different run is that with $\phi_0=2\pi$ where the poleward mode near the equator grows much faster than in the other cases. However, after $t {u_{\rm rms}}{k_{\rm f}}\approx1500$ the equatorward mode takes over similarly as in the runs with a smaller $\phi_0$.

The velocity field shows no marked evidence of low-degree non-axisymmetric constituents, but there are indications of $m=1$ structures in the instantaneous magnetic field (Figure \[fig:pmoll\_uu1bb3\]); see also http://youtu.be/u55sAtN2Fqs for an animation of the toroidal magnetic field. This is also reflected by the fraction of the axisymmetric part of the total magnetic energy; see Columns 5 and 6 of Table \[tab:runs2\]. We find that the energy of the mean toroidal field decreases monotonically when $\phi_0$ is increased so that there is a factor of three in $E_{\rm tor}/E_{\rm mag}$ between the extreme cases of Runs E1 and E4. The axisymmetric part still exhibits an oscillatory mode with equatorward migration in all runs in Set E. The most prominent exception is visible in Figure \[fig:butterflyE4\], where we show the butterfly diagram of the $m=0$ contribution for Run E4. Clearly, equatorward migratory events are now rare and superimposed on a background of small-scale, high-frequency poleward migratory field.

![image](pps.eps){width="\textwidth"}

We compute power spectra of the azimuthal component of the magnetic field from the Run E4 over three $10\degr$ latitude strips from each hemisphere, centered around latitudes of $\pm25\degr$, $\pm45\degr$, and $\pm65\degr$. The results for the three lowest degrees $m=0, 1, 2$ are shown in Figure \[fig:pps\]. We find that at low ($\pm25\degr$) and high ($\pm65\degr$) latitudes the axisymmetric ($m=0$) mode begins to dominate after around $1000$ turnover times and shows a cyclic pattern consistent with that seen in the time-latitude diagram of the azimuthally averaged field. After $t {u_{\rm rms}}{k_{\rm f}}\approx 1600$, however, the $m=1$ mode becomes stronger in the southern hemisphere, coinciding with the growth of the $m=1$ mode at mid-latitudes ($\pm45\degr$) where it dominates earlier in both hemispheres. This is in rough agreement with some observational results of rapid rotators, which show the most prominent non-axisymmetric temperature [e.g., @Thomas01; @Heidi07; @Marjaana11] and magnetic structures [@Oleg11] at the latitudinal range around 60$\degr$–80$\degr$, while the equatorial and polar regions are more axisymmetric; some temperature inversions even show almost completely axisymmetric distributions in the polar regions and rings of azimuthal field at low latitudes [e.g., @Donati03]. The strength of the axisymmetric versus the non-axisymmetric part in such objects has also been reported to vary over time with a time scale of a few years [@Oleg11].

Irradiance variations {#sec:irra}
---------------------

In contrast to the constant temperature condition used earlier, the black body boundary condition (\[eq:bbb\]) allows the temperature to vary at the surface of the star and thus enables the study of irradiance variations due to the magnetic cycle [@Spr00]. Such variations might even be responsible for driving torsional oscillations in the Sun [@Spr03; @Re06]. In Figure \[fig:pTTmxy\] we compare time–latitude surface representations ($r=R$) of azimuthally averaged temperature variations relative to its temporal average, $\Delta{\overline{T}}(\theta,t)={\overline{T}}(\theta,t)-{\langle {\overline{T}} \rangle}_t(\theta)$, with those of the azimuthally averaged radial magnetic field, ${\overline{B}}_r(\theta,t)$, for Run C1 in the saturated state of the dynamo. We also show scatter plots of $\Delta{\overline{T}}/{\overline{T}}$ versus ${\overline{B}}_r/B_{\rm eq}$ at $\pm70^\circ$ and $\pm30^\circ$ latitude to demonstrate that there are many instances where enhanced surface magnetic activity leads to a local decrease in surface temperature. We see that $$\Delta{\overline{T}}/{\overline{T}} \approx -Q_T {\overline{B}}_r^2/B_{\rm eq}^2$$ with ‘quenching’ coefficients $Q_T$ of $\approx0.14$ at high latitudes and $\approx0.33$ at low latitudes. However, there is also considerable scatter, even though our data is already longitudinally averaged. Without such averaging, the correlation between individual structures on the surface would be rather poor. The temperature modulation is best seen near the poles; see Figure \[fig:pTTmxy\]. This could be a consequence of a strong radial magnetic field that builds up some 50–100 turnover times earlier and thus precedes the temperature signal. A weaker modulation is also seen near the equator. The peak values of $\Delta{\overline{T}}/{\overline{T}}$ at high latitudes are 15%–20% of the surface temperature; see the last two panels of Figure \[fig:pTTmxy\]. This is relatively large compared with earlier work using mean-field models [@BMT92], which showed remarkably little relative variation of the order of $10^{-3}$ in the bulk of the convection zone and even less at the surface. This difference in the modulation amplitude is probably related to the importance of latitudinal variations that were also present in the mean-field model of [@BMT92] and referred to as thermal shadows [@Pa87b].

Conclusions {#sec:conclusions}
===========

We have studied the effects of density stratification on the dynamo solutions found in simulations of rotating turbulent convection in spherical wedge geometry for four values of $\Gamma_\rho$, which is the ratio of the densities at the bottom and at the surface of the convection zone. In addition, we vary the rotation rate for each value of $\Gamma_\rho$. For all stratifications we find quasi-steady large-scale dynamos for lower rotation and oscillatory solutions when rotation is rapid enough. The transition from quasi-steady to oscillatory modes seems to occur at a lower ${{\rm Co}}$ for higher stratification. Furthermore, for low values of $\Gamma_\rho$ the oscillatory solutions show only poleward propagation of the activity belts whereas at higher $\Gamma_\rho$ an equatorward branch appears at low latitudes.

![ Top panel: azimuthally averaged temperature fluctuations, normalized with their temporal average, $\Delta{\overline{T}}/{\overline{T}}$, at the surface as a function of time. Second panel: azimuthally averaged radial field at the surface, ${\overline{B}}_r/B_{\rm eq}$, as function of time. Lower panels: scatter plots of $\Delta{\overline{T}}/{\overline{T}}$ vs.${\overline{B}}_r/B_{\rm eq}$ at $\pm70^\circ$ and $\pm30^\circ$ latitude. Here, solid (dashed) red lines and filled (open) symbols refer to northern (southern) latitudes. The slopes are $Q_T\approx0.14$ and 0.33 at $\pm70$ and $\pm30$ latitude. All plots show quantities from the saturated stage of Run C1. []{data-label="fig:pTTmxy"}](BTcorr_small.ps){width="\columnwidth"}

The equatorward branch was first noted by [@KMB12a] using a wedge with $\phi_0=90\degr$ longitude extent. Here we test the robustness of this result by varying $\phi_0$ from $45\degr$ to full $360\degr$. We find a very similar pattern of the axisymmetric part of the field in all cases. However, the energy of the axisymmetric magnetic field decreases with increasing $\phi_0$. In the simulation with the full $\phi$-extent of $2\pi$ we observe an $m=1$ mode which is visible even by visual inspection (see Figure \[fig:pmoll\_uu1bb3\]). Such field configurations have been observed in rapidly rotating late-type stars [see e.g., @Oleg11] and our simulation is one of the first to reproduce such features [see also @GD08; @GDW12]. We are currently investigating the rapid rotation regime with more targeted runs which will be reported in a separate publication [@CKMB13].

The ratio between cycle to rotation frequency, $\omega_{\rm cyc}/\Omega_0$, is argued to be an important non-dimensional output parameter of a cyclic dynamo. For the Sun and other relatively inactive stars, this ratio is around 0.01, while for the more active stars it is around 0.002. For our models we find values in the range 0.002–0.01, but for most of the runs it is around 0.004. Although it is premature to make detailed comparisons with other stars and even the Sun, it is important to emphasize that kinematic mean-field dynamos produce the correct cycle frequency only for values of the turbulent magnetic diffusivity that are at least 10 times smaller than what is suggested by standard estimates [@Ch90]. In our case, these longer cycle periods (or smaller cycle frequencies) might be a result of nonlinearity as they are only obtained in the saturated regime of the dynamo. The detailed reason for this is unclear, but it has been speculated that it is connected with a slow magnetic helicity evolution [@Br05]. On the other hand, magnetic helicity effects are expected to become important only at values of ${{\rm Rm}}$ between 100 and 1000 [@DSGB13], which is much larger than what has been reached in the present work. Equally unclear is the reason for equatorward migration, which, as we have seen, might be a consequence of nonlinearity, as well. It will therefore be important to provide an accurate determination of all the relevant turbulent transport coefficients. The explanation favored in the present paper is that the dynamo wave is that expected for an oscillatory $\alpha^2$ dynamo caused by the change of sign of $\alpha$ about the equator. This is evidenced by our finding that ${\overline{B}}_\phi$ lags ${\overline{B}}_r$ by about $\pi/2$, which cannot be explained by an $\alpha^2\Omega$ dynamo.

We thank Matthew Browning, Martin Schrinner, L$\acute{\rm u}$cia Duarte, and Matthias Rheinhardt, as well as the referee for making a number of useful suggestions. The simulations were performed using the supercomputers hosted by CSC – IT Center for Science Ltd. in Espoo, Finland, who are administered by the Finnish Ministry of Education. Financial support from the Academy of Finland grants No. 136189, 140970 (PJK) and 218159, 141017 (MJM), as well as the Swedish Research Council grant 621-2007-4064, and the European Research Council under the AstroDyn Research Project 227952 are acknowledged as well as the HPC-Europa2 project, funded by the European Commission - DG Research in the Seventh Framework Programme under grant agreement No. 228398. The authors thank NORDITA for hospitality during their visits.

[^1]: Note that in Equation (4) of [@KMB12a] the Ohmic heating term $\mu_0 \eta {\bm J}^2$ and a factor $\rho$ in the viscous dissipation term $2\nu \bm{\mathsf{S}}^2$ were missing, but they were actually included in the calculations.

[^2]: http://pencil-code.googlecode.com/

[^3]: Note that the values of ${{\rm Re}}$ and ${{\rm Co}}$ have been recalculated with the same definition of ${u_{\rm rms}}$ as in the current paper.

[^4]: [http://www.nordita.org/\$\^\\ensuremath{\\sim}\$brandenb/PencilCode/MeanFieldSpherical.html](http://www.nordita.org/$^\ensuremath{\sim}$brandenb/PencilCode/MeanFieldSpherical.html)
