Free Access
Issue
A&A
Volume 599, March 2017
Article Number A4
Number of page(s) 18
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201628973
Published online 20 February 2017

© ESO, 2017

1. Introduction

Simulations of convection-driven dynamos have recently reached a level of sophistication where they capture effects observed in the Sun such as equatorward migration of activity belts (Schrinner et al. 2011; Käpylä et al. 2012; Augustson et al. 2015; Duarte et al. 2016) and irregular cycle variations such as grand minima and long-term modulations (Passos & Charbonneau 2014; Käpylä et al. 2016a). Most of these simulations are individual numerical experiments, and it is not clear how they are situated in parameter space in relation to each other. Important parameters in this connection concern the relative strengths of different diffusion coefficients, viscosity (ν), and magnetic (η) and thermal (χ) diffusivities present in the system. Their ratios are characterized by the thermal and magnetic Prandtl numbers, Pr = ν/χ and PrM = ν/η, respectively. In the solar convection zone, these Prandtl numbers are Pr ≪ 1 and PrM ≪ 1, while the fluid and magnetic Reynolds numbers, Re = ul/ν and ReM = ul/η, with u and l being the characteristic velocity and length scale, are on the orders of 1012 and 109, respectively. Such parameter regimes are not accessible to current numerical simulations, which are restricted to Pr ≈ 1, PrM ≈ 1, and Reynolds numbers on the order of 102−103. In all simulations by a number of different groups, the dominant contribution to thermal diffusion is made by a subgrid-scale (SGS) coefficient χSGS whose magnitude is much higher than the radiative one. Similar arguments also apply to ν and η, but since the functional form of these diffusion operators is unchanged, we omit the subscript SGS in them. Thus, the relevant thermal Prandtl number in simulations is PrSGS = ν/χSGS. We emphasize that this applies to simulations of all groups, although the nomenclature may be different (see Table A.1). This is also true for groups using realistic luminosities, and thus the correct order of magnitude for the radiative diffusivity (e.g., Brun et al. 2004; Hotta et al. 2016).

When the convection simulations of Käpylä et al. (2012) in wedge geometry showed clear equatorward migration for the first time, it was not obvious that this was related to their choice of PrSGS = 2.5 compared with PrSGS ≲ 1 used in most earlier simulations that showed either quasi-stationary configurations (Brown et al. 2010) or either weak or poleward migration (e.g., Ghizaru et al. 2010; Käpylä et al. 2010b; Brown et al. 2011; Gastine et al. 2012). Recently, Warnecke et al. (2014) showed that the change in the dynamo behavior between PrSGS> 1 and PrSGS< 1 regimes is due to a change in the differential rotation profile, which in the PrSGS ≳ 1 regime leads to a region of negative radial shear that facilitates the equatorward migration. The magnetic Prandtl number, which is proportional to the magnetic Reynolds number, can also strongly affect the results. Increasing ReM by increasing PrM can allow magnetohydrodynamic (MHD) instabilities such as magnetic buoyancy (Parker 1955) and magnetorotational instabilities (e.g., Parfrey & Menou 2007; Masada 2011) to be excited.

Increasing the magnetic Reynolds number can influence the large-scale dynamo by several other avenues. First, the most easily excited dynamo mode can change. Second, a small-scale dynamo is most likely excited after ReM exceeds a threshold value (e.g., Cattaneo 1999), and this may also affect the large-scale dynamo by modifying the velocity field. Third, Boussinesq simulations indicate that differential rotation is strongly quenched as the magnetic Reynolds number increases (Schrinner et al. 2012). This was shown to be associated with a transition from oscillatory multipolar large-scale field configurations to quasi-stationary dipole-dominated dynamos as a function of ReM. One of the main goals of the present paper is therefore to systematically study the effects of varying Prandtl numbers on the differential rotation and dynamo modes excited in the simulations. We note that similar parameter studies have been performed with Boussinesq simulations (e.g., Simitev & Busse 2005; Busse & Simitev 2006). Here we explore the stratified, fully compressible simulations and reach parameter regimes that are significantly more supercritical in terms of both the convection and the dynamo.

Another important aspect is related to the saturation level of the large-scale field in simulations at high ReM. Dynamo theory experienced a crisis in the early 1990s when it was discovered that the energy of the large-scale magnetic field saturates at a level that is inversely proportional to the magnetic Reynolds number (Gruzinov & Diamond 1995; Brandenburg & Dobler 2001). If this were to carry over to the Sun, where the magnetic Reynolds number is on the order of 109 or greater, only very weak large-scale fields would survive. This phenomenon was related to a catastrophic quenching of the α effect (Cattaneo & Vainshtein 1991; Vainshtein & Cattaneo 1992; Cattaneo & Hughes 1996).

Later, this was understood in terms of magnetic helicity: if the system is closed or fully periodic, that is, when magnetic field lines do not cross the boundary of the system, no flux of magnetic helicity in or out can occur, and only the molecular diffusion can change it (Brandenburg 2001). In astrophysical systems this would mean that magnetic helicity would be nearly conserved. However, astrophysical systems are not closed and magnetic helicity can escape, for example, through coronal mass ejections in the Sun (e.g., Blackman & Brandenburg 2003; Warnecke et al. 2011, 2012) or through winds from galaxies (Shukurov et al. 2006; Sur et al. 2007; Del Sordo et al. 2013). In mean-field theory these physical effects are parameterized by fluxes, which lead to alleviation of catastrophic quenching in suitable parameter regimes (e.g., Brandenburg et al. 2009).

Direct numerical simulations of large-scale dynamos have demonstrated that open boundaries lead to alleviation of catastrophic quenching in accordance with the interpretation in terms of magnetic helicity conservation (e.g., Brandenburg & Sandin 2004; Käpylä et al. 2010a). Although the large-scale magnetic field amplitude does not decrease proportional to ReM in the cases when open boundaries are used, there is still a decreasing trend even at the highest currently studied ReM in local simulations of convection-driven dynamos (e.g., Käpylä et al. 2010a). This, however, is compatible with mean-field models, which suggest that the magnetic helicity fluxes become effective only at significantly higher ReM (Brandenburg et al. 2009; Del Sordo et al. 2013).

In convective dynamos in spherical coordinates the computational challenge is even greater and systematic parameter scans have not been performed. Some preliminary attempts have been made, but the results remain inconclusive. An illuminating example is the study of Nelson et al. (2013), where the large-scale axisymmetric field decreases by a factor of two when the magnetic Reynolds number is increased by a factor of four, which is still rather steep. In a recent paper, Hotta et al. (2016) showed that in even higher-ReM simulations the mean magnetic energy recovers, and the authors claim that this is a consequence of an efficient small-scale dynamo that suppresses small-scale flows. Another goal of the present paper is therefore to study the saturation level of the large-scale field in convection-driven dynamos in spherical coordinates with and without a simultaneous small-scale dynamo (hereafter SSD).

In the present study, we employ a spherical wedge geometry by imposing either a perfect conductor or a normal field boundary condition at high latitudes. Earlier mean-field simulations of αΩ dynamos have suggested that solutions with a perfect conductor boundary condition are similar to those in full spherical shells (Jennings et al. 1990). However, more recent work by Cole et al. (2016) has demonstrated that this conclusion is not generally valid and depends on the nature of the solutions. Their work also suggests that the use of a normal field boundary condition at high latitudes might be a better way of obtaining solutions that are applicable to full spherical shells. Owing to this uncertainty, we investigate here cases with both types of boundary conditions.

2. Model

Our model is similar to that of Käpylä et al. (2012) and is described in detail in Käpylä et al. (2013). The model describes the dynamics of magnetized gas in spherical coordinates where only parts of the latitude and longitude ranges are retained. More specifically, the model covers a wedge that spans r0rr1 in radius, θ0θπθ0 in colatitude, and 0 ≤ φφ0 in longitude. Here we use r0 = 0.7 R, r1 = R, and where R = 6.96 × 108 m is the solar radius, θ0 = 15°, and φ0 = 90°.

We solve the following set of compressible hydromagnetics equations A∂t=U×Bημ0J,DlnρDt=·U,DUDt=g2Ω0×U+1ρ[·(2νρS)p+J×B],\begin{eqnarray} &&\frac{\pd {\bm A}}{\pd t} = {\bm U} \times {\bm B} - \eta \mu_0 {\bm J},\\[1mm] &&\frac{D \ln \rho}{Dt} = - \bm\nabla\bm\cdot{\bm U}, \\[1mm] &&\frac{D {\bm U}}{Dt} = {\bm g} - 2\bm\Omega_0\times \bm U + \frac{1}{\rho} [\bm\nabla\!\bm\cdot\!(2\nu\rho \bm{\mathsf{S}}) - \bm\nabla p + {\bm J} \times {\bm B}], \\[1mm] &&T \frac{Ds}{Dt} = \frac{1}{\rho} \left[\eta \mu_0 {\bm J}^2-\bm\nabla\bm\cdot({\bm F}^{\rm rad}+{\bm F}^{\rm SGS})\right]+2\nu \bm{\mathsf{S}}^2, \label{equ:1} \end{eqnarray}where A is the magnetic vector potential, U is the velocity, B = × A is the magnetic field, η is the magnetic diffusivity, μ0 is the permeability of vacuum, J = × B/μ0 is the current density, D/Dt = /∂t + U· is the advective time derivative, ρ is the density, g is the acceleration due to gravity, and Ω0 = (cosθ,−sinθ,0)Ω0 is the angular velocity vector where Ω0 is the rotation rate of the frame of reference, ν is the kinematic viscosity, p is the pressure, and s is the specific entropy with Ds = cVDlnpcPDlnρ. The gas is assumed to obey the ideal gas law, p = (γ−1)ρe, where e = cVT is the specific internal energy and γ = cP/cV is the ratio of specific heats at constant pressure and volume, respectively. The rate of strain tensor is given by Sij=12(Ui;j+Uj;i)13δij·U,\begin{eqnarray} \mathsf{S}_{ij} = \onehalf (U_{i;j} + U_{j;i}) - \onethird \delta_{ij} \bm\nabla\bm\cdot {\bm U}, \end{eqnarray}(5)where the semicolons refer to covariant derivatives; see Mitra et al. (2009). The radiative and SGS fluxes are given by Frad=KT,FSGS=χSGSρTs,\begin{eqnarray} {\bm F}^{\rm rad} = -K\bm\nabla T, \hspace{1cm }{\bm F}^{\rm SGS} = -\chiSGS \rho T \bm\nabla s, \label{equ:Frad} \end{eqnarray}(6)where K = cPρχ is the heat conductivity, and χSGS is the (turbulent) SGS diffusion coefficient for the entropy.

Table 1

Summary of the runs.

thumbnail Fig. 1

Top row: temporally averaged rotation profiles from Runs A3, B3, C3, and D3 with PrM = 1 and PrSGS varying from 0.25 (left) to 2 (right). Lower row: same as above, but for Runs C2, C3, C4, and C5 with PrSGS = 1 and PrM varying from 0.5 (left) to 5 (right).

2.1. Initial and boundary conditions

We follow the description given in Käpylä et al. (2014) to transform our results into physical units. As our equations are fully compressible, we cannot afford to use the solar luminosity, which would lead to the sound speed dominating the time step. Thus we increase the luminosity substantially, and to compensate for this, we increase the rotation rate to achieve the same rotational influence as in the Sun. Assuming a scaling of the luminosity with the convective energy flux as Lρu3, we find that the convective velocities increase to one-third power of the luminosity (e.g., Brandenburg et al. 2005; Karak et al. 2015). The ratio between model and solar luminosities is L0/L ≈ 6.4 × 105 in the current simulations. We correspondingly increase the rotation rate by a factor of (L0/L)1386\hbox{$(L_0/L_\odot)^{\frac{1}{3}}\approx 86$} for the solar case. We reiterate that the main effect of increasing the luminosity is to increase the Mach number and bring the acoustic and dynamical timescales closer to each other, which facilitates the computations with a fully compressible method (cf. Fig. 1 of Käpylä et al. 2013). Another possibility would be to apply the so-called reduced sound speed method where the sound speed is artificially changed so that the Mach number does not become too small at the base of the convection zone (Hotta et al. 2012; Käpylä et al. 2016b).

The higher luminosity also helps to reach a thermal equilibrium within reasonable simulation running time. Furthermore, we assume that the density and the temperature at the base of the convection zone are the same as in the Sun, that is, ρbot = 200 kg m-3 and T = 2.23 × 106 K.

As the initial condition for the hydrodynamics we use an isentropic atmosphere. The temperature gradient is given by ∂T∂r=GM/r2cV(γ1)(nad1),\begin{equation} \frac{\pd T}{\pd r} = -\frac{GM_\odot/r^2}{\cV(\gamma-1)(n_{\rm ad}-1)}, \end{equation}(7)where G = 6.67 × 10-11 N m2 kg-2 is Newton’s gravitational constant, M = 1.989 × 1030 kg is the solar mass, and nad = 1.5 is the polytropic index for an adiabatic stratification.

The initial state is not in thermal equilibrium but closer to the final convecting state to reduce the time needed to reach a statistically stationary state. The heat conductivity has a profile given by K(r) = K0 [n(r) + 1], where n(r) = δn(r/r0)15 + nadδn, with K0=(/4π)cV(γ1)(nad+1)ρbotGMR\hbox{$K_0=(\mathcal{L}/4\pi)\cV(\gamma-1)(n_{\rm ad}+1)\rho_{\rm bot} \sqrt{GM_\odot \,R_\odot}$}, and where is a dimensionless luminosity defined below. We keep δn = 1.9 fixed in our simulations, resulting a situation where radiation transports all of the flux into the domain, but its contribution diminishes rapidly toward the surface (see, e.g., Fig. 2 of Käpylä et al. 2011a).

The SGS diffusivity χSGS for the entropy has a piecewise constant profile, such that the value χSGS = 6.1 × 108 m s-2 is fixed above r = 0.98 R in all runs, and the value χSGSm=χSGS(r=rm=0.85R)\hbox{$\chiSGSm=\chiSGS(r=r_m=0.85\,R_\odot)$} in the bulk of the convection zone is varied by the corresponding Prandtl number (see below). The value below r = 0.75 R is set equal to 0.1χSGSm\hbox{$0.1\chiSGSm$}. The constant values in the different layers connect smoothly over a transition depth of d = 0.01 R.

The boundary conditions for the flow are assumed to be impenetrable and stress free, that is, Ur=0,Uθ∂r=Uθr,Uφ∂r=Uφr(r=r0,r1),Ur∂θ=Uθ=0,Uφ∂θ=Uφcotθ(θ=θ0θ0).\begin{eqnarray} && U_r=0,\quad \frac{\pd U_\theta}{\pd r}=\frac{U_\theta}{r},\quad \frac{\pd U_\phi}{\pd r}=\frac{U_\phi}{r} \quad (r=r_0, r_1),\\[2mm] && \frac{\pd U_r}{\pd \theta}=U_\theta=0,\quad \frac{\pd U_\phi}{\pd \theta}=U_\phi \cot \theta \quad (\theta=\theta_0,\pi-\theta_0). \quad \end{eqnarray}The lower radial boundary is assumed to be perfectly conducting, and on the outer radial boundary the field is purely radial. The latitudinal boundaries are either perfectly conducting (PC) or a normal field (NF) condition is assumed. In terms of the magnetic vector potential, these are given by Ar∂r=Aθ=Aφ=0(r=r0),Ar=0,Aθ∂r=Aθr,Aφ∂r=Aφr(r=r1),Ar=Aθ∂θ=Aφ=0(θ=θ0θ0)(PC),Ar∂θ=Aθ=0,Aφ∂θ=cotθAφr(θ=θ0θ0)(NF).\begin{eqnarray} &&\frac{\pd A_r}{\pd r}= A_\theta=A_\phi =0 \,\quad (r=r_0),\\ &&A_r=0, \;\; \frac{\pd A_{\theta}}{\pd r}=-\frac{A_{\theta}}{r},\;\; \frac{\pd A_{\phi}}{\pd r}=-\frac{A_{\phi}}{r} \quad (r=r_1),\quad\\ &&A_r=\frac{\pd A_\theta}{\pd\theta}=A_\phi=0 \quad (\theta=\theta_0,\pi-\theta_0) \quad \mbox{(PC)},\\ &&\frac{\pd A_r}{\pd\theta}\!=\!A_\theta\!=\!0, \frac{\pd A_\phi}{\pd\theta}\!=\!-\frac{\cot\theta A_\phi}{r}\ (\theta=\theta_0,\pi-\theta_0) \ \mbox{(NF)}. \end{eqnarray}For the density and specific entropy we assume vanishing first derivatives on the latitudinal boundaries.

At the lower boundary we specify Frrad+FrSGS=L04πr02,\begin{equation} F_r^{\rm rad}+F_r^{\rm SGS}=\frac{L_0}{4\pi r_0^2}, \end{equation}(14)which leads to constant input luminosity into the system. On the outer radial boundary we apply a radiative boundary condition σT4=Frrad+FrSGS,\begin{equation} \sigma T^4 = F_r^{\rm rad}+F_r^{\rm SGS}, \label{eq:bbb} \end{equation}(15)where σ is the Stefan-Boltzmann constant. We use a modified value for σ that takes into account that the luminosity and the temperature at the surface are higher than in the Sun. The value of σ is chosen so that the surface flux, σT4, carries the total luminosity through the boundary in the initial non-convecting state.

2.2. System parameters and diagnostics quantities

The parameters governing our models are the dimensionless luminosity =L0ρ0(GM)3/2R1/2,\begin{equation} \mathcal{L} = \frac{L_0}{\rho_0 (GM_\odot)^{3/2} \,R_\odot^{1/2}}, \end{equation}(16)the normalized pressure scale height at the surface, ξ=(γ1)cVT1GM/R,\begin{equation} \xi = \frac{(\gamma-1) c_{\rm V}T_1}{GM_\odot/R_\odot}, \end{equation}(17)with T1 being the temperature at the surface in the initial state, the Taylor number Ta=(2Ω0Δr2/ν)2,\begin{equation} \Ta=(2\Omega_0 \Delta r^2/\nu)^2, \end{equation}(18)where Δr = r1r0 = 0.3 R, as well as the fluid, SGS, and magnetic Prandtl numbers Pr=νχm,PrSGS=νχSGSm,PrM=νη,\begin{equation} \Pra=\frac{\nu}{\chi_{\rm m}},\quad \Pra_{\rm SGS}=\frac{\nu}{\chiSGSm},\quad \PrM=\frac{\nu}{\eta}, \end{equation}(19)where χm = K(rm) /cPρm is the thermal diffusivity and ρm is the density, both evaluated at r = rm. We keep Pr = 71.7 fixed and vary PrSGS and PrM in most of our models, with the exception of Set G, where PrSGS and PrM are set to unity and Pr is varied by changing the value of ν (see Table 1). The Rayleigh number is defined as Ra=GM(Δr)4νχSGSmR2(1cPdshsdr)rm,\begin{equation} \Ra=\frac{GM_\odot(\Delta r)^4}{\nu \chiSGSm R_\odot^2} \bigg(-\frac{1}{c_{\rm P}}\frac{{\rm d} s_{\rm hs}}{{\rm d}r} \bigg)_{r_{\rm m}}, \label{equ:Co} \end{equation}(20)where shs is the entropy in the hydrostatic (hs), non-convecting state obtained from a one-dimensional model where no convection can develop with the prescriptions of K and χSGS given above.

The remaining parameters are used only as diagnostics. These include the fluid and magnetic Reynolds numbers, and the Péclet number Re=urmsνkf,ReM=urmsηkf,Pe=urmsχSGSmkf,\begin{equation} \Rey=\frac{\urms}{\nu \kf},\quad \ReM=\frac{\urms}{\eta \kf},\quad \Pe=\frac{\urms}{\chiSGSm \kf}, \end{equation}(21)where kf=2π/Δr21R-1\hbox{$\kf=2\pi/\Delta r\approx21 \,R_\odot^{-1}$} is an estimate of the wavenumber of the largest eddies. Rotational influence on the flow is given by the Coriolis number Co=2Ω0urmskf,\begin{equation} \Co=\frac{2\Omega_0}{\urms \kf}, \label{eq:Coriolis} \end{equation}(22)where urms=(3/2)Ur2+Uθ2rθφt\hbox{$\urms=\sqrt{(3/2)\brac{U_r^2+U_\theta^2}_{r\theta\phi t}}$} is the rms velocity and the subscripts indicate averaging over r, θ, φ, and a time interval during which the run is thermally relaxed. We omit the contribution from the azimuthal velocity in urms because it is dominated by differential rotation (Käpylä et al. 2011b).

We define mean quantities as averages over the φ-coordinate and denote them by overbars. We also often average the data in time over the period of the simulations where thermal energy, differential rotation, and large-scale magnetic fields have reached statistically saturated states.

The simulations are performed with the Pencil Code1, which uses a high-order finite difference method for solving the compressible equations of magnetohydrodynamics.

3. Data analysis: D2 statistic

To detect possible cycles and to estimate their average lengths, we have chosen to use D2 phase dispersion statistic (Pelt 1983). It has recently been applied to irregularly spaced long-term photometry of solar-like stars (Lindborg et al. 2013; Olspert et al. 2015) as well as to more regularly sampled magnetoconvection simulation data (Karak et al. 2015; Käpylä et al. 2016a). In the previous applications the statistic has been used exclusively on one-dimensional time series (e.g., by fixing a certain latitude and radius in the azimuthally averaged data). In the current study we use a generalized form of the statistic given by D2(P,Δtcoh)=i=1N1j=i+1Ng(ti,tj,P,Δtcoh)||f(ti)f(tj)||22σ2i=1N1j=i+1Ng(ti,tj,P,Δtcoh),\begin{equation} \label{eq:D2} D^2(P,\Delta t_{\rm coh})\!=\!\frac{\sum\limits_{i=1}^{N-1} {\sum\limits_{j=i+1}^N {g(t_i}} ,t_j ,P,\Delta t_{\rm coh})||{\bm f}(t_i ) - {\bm f}(t_j )||^2}{2\sigma^2\sum\limits_{i = 1}^{N - 1} {\sum\limits_{j = i + 1}^N {g(t_i } } ,t_j ,P,\Delta t_{\rm coh})}, \end{equation}(23)where f(ti) is the vector of observed variables at time moment ti, σ2 = N-2i,j>i | | f(ti)−f(tj) | | 2 is the variance of the full time series, g(ti,tj,P,Δtcoh) is the selection function, which is significantly greater than zero only when tjtikP,k=±1,±2,...,and|tjti|Δtcoh=lcohP,\begin{eqnarray} t_j - t_i &\approx& kP,k = \pm 1, \pm 2, \ldots, \quad \mbox{and}\\ \left| {t_j - t_i } \right| &\lessapprox& \Delta t_{\rm coh} = l_{\rm coh}P, \end{eqnarray}where P is the trial period and Δtcoh is the so-called coherence time, which is the measure of the width of the sliding time window wherein the data points are taken into account by the statistic. The number of trial periods fitting into this interval, lcoh = Δtcoh/P, is called the coherence length. With this definition there is no restriction on the dimensionality of the data, but Eq. (23) leaves open the choice of the vector norm. In most cases it is natural to use the Euclidean norm, which we also assume in our analysis. In the following we use this statistic to analyze the radial and azimuthal components of the magnetic field at regions near the surface over latitude intervals, where the cycles are the most pronounced (see Sect. 4.2.2).

thumbnail Fig. 2

Estimates of radial and latitudinal differential rotation ΔΩ(r)\hbox{$\Delta_\Omega^{(r)}$} and ΔΩ(θ)\hbox{$\Delta_\Omega^{(\theta)}$}, respectively, according to Eq. (26) for Sets A–F as indicated by the legend.

4. Results

We performed six sets of simulations (Sets A–F), each with a constant value of PrSGS but changing PrM in the range 0.25 ≤ PrM ≤ 5; see Table 1. Furthermore, in an additional Set G we fixed PrSGS = PrM = 1 and varied the Reynolds and Péclet numbers. The rotation rate was varied such that Sets A–D and G have Ω0 = 5Ω, whereas in Sets E and F we used and Ω, respectively. Sets E and F were included to study the reliability of our findings at slower rotation. In some sets (A, D, E, and F) numerical problems prevented us from using PrM = 5. In these cases we used a lower value that produces numerically stable solutions. The latitudinal PC boundary conditions led to numerical problems as a result of an unidentified instability at high magnetic Reynolds numbers near the latitudinal boundaries. Cases where this occurred were rerun with the NF conditions. We did not find major qualitative differences in the behavior of the large-scale field between PC and NF runs with otherwise identical parameters. In the following, we omit the PrM = 0.25 runs that did not lead to dynamos in Sets E and F. We note that Run B2 has been presented as Run II in Warnecke et al. (2014), Run D3 as Run A1 in Warnecke et al. (2016a,b), and Run C3 covers the first 120 yr of the run presented in Käpylä et al. (2016a). Furthermore, runs similar to Run D3 (but with PrM = 2.5 instead of 2.0) are presented as Runs B4m and C1 in Käpylä et al. (2012, 2013).

4.1. Large-scale flows and their generators

4.1.1. Differential rotation and meridional circulation

The sign of the radial gradient of Ω plays a crucial role in determining the propagation direction of dynamo waves in αΩ dynamos: with a positive α effect in the northern hemisphere, a negative radial gradient of Ω is required to obtain solar-like equatorward migration and vice versa (Parker 1955; Yoshimura 1975). It is remarkable that this rule also seems to apply to the fully nonlinear convective dynamo simulations (Warnecke et al. 2014, 2016b). The current simulations can produce equatorward migration only in cases where a region with a negative radial gradient of Ω occurs at mid-latitudes (e.g., Käpylä et al. 2012, 2013; Augustson et al. 2015) or if the sign of the kinetic helicity, which is a proxy of the α effect, is inverted in the bulk of the convection zone (Duarte et al. 2016).

The simulations of Käpylä et al. (2012), Augustson et al. (2015), Käpylä et al. (2016a), and Warnecke et al. (2016a) showing equatorward migration and a region of negative radial shear at mid-latitudes have PrSGS ≳ 1. This is in contrast to earlier simulations with PrSGS< 1, which did not show equatorward migration (e.g., Brun et al. 2004; Brown et al. 2011; Nelson et al. 2013) and had consistently positive gradients of Ω. Given that the dynamo wave propagation is apparently heavily influenced by this, it is important to study the effect that PrSGS has on the rotation profiles.

We show representative results of the temporally averaged rotation profiles as a function of PrSGS from Runs A3, B3, C3, and D3 in the top row of Fig. 1. In the lowest SGS Prandtl number case (Run A3), the angular velocity Ω=Ω0+Uφ/rsinθ\hbox{$\Omega=\Omega_0 + \mean{U}_\phi/r\sin\theta$} decreases monotonically from the equator toward the poles. Much of the latitudinal variation occurs at high latitudes near the latitudinal boundaries. However, the rotation profile is qualitatively similar to those obtained from low PrSGS models in fully spherical shells (e.g., Brun et al. 2004; Brown et al. 2010). In Run B3, a dip at mid-latitudes is developing, which is seen to deepen in the higher PrSGS Runs C3 and D3. We note that a similar transition occurs also when the density stratification is increased with PrSGS = 2...5 (Käpylä et al. 2011a, 2013). The overall magnitude of the differential rotation also decreases from low to high SGS Prandtl numbers. However, much of this variation occurs already between Runs A3 and B3, whereas the differences between Runs B3, C3, and D3 are much smaller.

The time-averaged rotation profiles from runs with PrSGS = 1 with varying PrM are shown in the lower row of Fig. 1 from Runs C2–C5. The absolute shear decreases steeply as PrM and ReM increase so that in the highest ReM case the differential rotation is appreciable only near the latitudinal boundaries. There are also qualitative changes such that the negative shear layer at mid-latitudes is almost absent in Run C5 and a near-surface shear layer is developing in Runs C4 and C5.

To quantify the radial and latitudinal differential rotation, we use the quantities (Käpylä et al. 2013) ΔΩ(r)=ΩeqΩbotΩeq,ΔΩ(θ)=ΩeqΩpoleΩeq,\begin{eqnarray} \Delta_\Omega^{(r)}=\frac{\Omega_{\rm eq}-\Omega_{\rm bot}}{\Omega_{\rm eq}}, \quad \Delta_\Omega^{(\theta)}=\frac{\Omega_{\rm eq}-\Omega_{\rm pole}}{\Omega_{\rm eq}},\label{equ:pDRt} \end{eqnarray}(26)where Ωeq = Ω(r1/ 2) and Ωbot = Ω(r0/ 2) are the rotation rates at the surface and at the bottom of the convection zone at the equator. Furthermore, Ωpole=12[Ω(r1,θ0)+Ω(r1θ0)]\hbox{$\Omega_{\rm pole}=\onehalf[\Omega(r_1,\theta_0)+\Omega(r_1,\pi-\theta_0)]$} is the average rotation rate between the latitudinal boundaries at the outer boundary. We show ΔΩ(r)\hbox{$\Delta_\Omega^{(r)}$} and ΔΩ(θ)\hbox{$\Delta_\Omega^{(\theta)}$} for the runs of Sets A–F in Fig. 2. We find that both radial and latitudinal differential rotation are modestly quenched for ReM ≲ 30 in Sets A–E. For higher values of ReM, both ΔΩ(r)\hbox{$\Delta_\Omega^{(r)}$} and ΔΩ(θ)\hbox{$\Delta_\Omega^{(\theta)}$} decrease steeply, and for the highest values of ReM the differential rotation is almost completely quenched. Set F with the lowest rotation rate stands apart from the other runs. The main difference in this set is that the differential rotation is anti-solar, that is, with a slow equator and faster poles. There the radial differential rotation decreases monotonically, but the decrease is only roughly 20 per cent in the range ReM = 20...151, whereas ΔΩ(θ)\hbox{$\Delta_\Omega^{(\theta)}$} remains roughly constant above ReM = 39. A possible explanation of the difference between Set F and the other sets is that the mean magnetic fields are stronger in the latter sets (see Table 2 and Sect. 4.2.4), which leads to a stronger back-reaction to the flow.

Table 2

Volume and time-averaged kinetic and magnetic energy densities realized in the simulations in units of 105 J m-3.

We list in Table 2 the kinetic energy densities of the total flow Ekin=12ρU2V\hbox{$E_{\rm kin}=\brac{\onehalf \mean{\rho \UUU^2}}_V$}, differential rotation EkinDR=12ρU2φV\hbox{$E_{\rm kin}^{\rm DR}=\brac{\onehalf \mean{\rho} \mean{U}_\phi^2}_V$}, meridional circulation EkinMC=12ρ(U2r+U2θ)V\hbox{$E_{\rm kin}^{\rm MC}=\brac{\onehalf \mean{\rho} (\mean{U}_r^2+\mean{U}_\theta^2)}_V$}, the fluctuating velocity Ekinfluct=12ρu2V\hbox{$E_{\rm kin}^{\rm fluct}=\brac{\onehalf \mean{\rho \uuu^2}}_V$}, and the magnetic energy densities related to the total field Emag=B2/2μ0V\hbox{$E_{\rm mag}=\brac{\mean{\BBB^2}/2\mu_0}_V$}, azimuthally averaged toroidal Emagtor=B2φ/2μ0V\hbox{$E_{\rm mag}^{\rm tor}=\brac{\mean{B}_\phi^2/2\mu_0}_V$} and poloidal fields Emagpol=(B2r+B2θ)/2μ0V\hbox{$E_{\rm mag}^{\rm pol}=\brac{(\mean{B}_r^2+\mean{B}_\theta^2)/2\mu_0}_V$}, and the fluctuating magnetic energy Emagfluct=b2/2μ0V\hbox{$E_{\rm mag}^{\rm fluct}=\brac{\mean{\bbb^2}/2\mu_0}_V$} with u=UU\hbox{$\uuu=\UUU-\meanv{U}$}, b=BB\hbox{$\bbb=\BBB-\meanv{B}$} and where ⟨ ⟩ V indicates a volume average in our simulations. We find that the kinetic energy decreases monotonically as the magnetic Reynolds number is increased, regardless of the SGS Prandtl number. This is mostly due to quenching of the differential rotation, whereas the fluctuating kinetic energy is much less affected. Differences between the sets of runs are large, however. In Set F, the total kinetic energy drops by less than a factor of two and the energy of the differential rotation by a factor of four, whereas in Sets A, B, C, and E, the Ekin reduces by roughly an order of magnitude and Ekin(DR)\hbox{$E_{\rm kin}^{\rm (DR)}$} by two orders of magnitude. Set D falls between the two extreme cases. The energy of the meridional flow is negligible in comparison to both differential rotation and fluctuating (non-axisymmetric) contributions.

We note, however, that the temporal variations of differential rotation increase as the magnetic Reynolds number is increased, see Table 2. This is demonstrated in Figs. 3a–c, where the energies of the differential rotation and mean toroidal magnetic field are shown for Runs G1–G3. In the lowest-ReM case, both are fairly stable with Emagtor\hbox{$E_{\rm mag}^{\rm tor}$} being typically an order of magnitude smaller than EkinDR\hbox{$E_{\rm kin}^{\rm DR}$} with a few excursions with strong magnetic field and weak differential rotation (e.g., around t ≈ 22 yr and t ≈ 50–80 yr). In an extended version of this run, the quiescent magnetic field leads to stronger-than-average differential rotation for the last 70 yr of that run (see Figs. 4a and 4b of Käpylä et al. 2016a). In Run G2 with ReM = 66, two more distinct states appear to be present: either the differential rotation is strong and the magnetic fields weak (t ≈ 5–20 and t ≳ 42 yr), or the two are comparable (t ≈ 25–42 yr). These events are associated with a change of the large-scale dynamo mode from an oscillatory equatorward migrating mode (strong DR, weak magnetic field) to a quasi-stationary one (DR and magnetic field energies comparable), see Fig. 3d. Panels e and f of Fig. 3 show that the latitudinal differential rotation decreases by roughly 30 per cent from the high to the low state. Similar but apparently more violent variations are seen in the highest-ReM case (Run G3), but there the time series is too short to draw reliable conclusions.

thumbnail Fig. 3

Energies of differential rotation (black solid lines) and mean toroidal magnetic field (red dashed) as functions of time from Runs G1, G2, and G3 (panels a)c)). Panel d) shows the azimuthally averaged azimuthal magnetic field at r = 0.98 R from Run G2. The vertical dashed lines indicate the high and low states of differential rotation. Panels e) and f) show the time-averaged rotation profiles in Run G2 from the high and low states indicated as blue dotted and dashed lines in panel b).

Our results appear to stand apart from similar studies in full spherical shells (e.g., Nelson et al. 2013; Hotta et al. 2016) in that the differential rotation is strongly quenched as a function of the magnetic Reynolds number. However, in Nelson et al. (2013) the values of Rm′ (=2πReM) correspond to a range of 8...32 in ReM in our units where the radial and latitudinal differential rotation decrease by about 30 per cent. This is roughly consistent with our results. On the other hand, Hotta et al. (2016) reached higher values of ReM than in the present study, but no strong quenching was reported. The reason might be that their models are rotating substantially slower than ours, leading to weaker magnetic fields and a weaker back-reaction to the flow. Furthermore, in these models, the differential rotation is strongly influenced by their SGS heat flux, which transports one third of the luminosity. Another obvious candidate for explaining the difference is the wedge geometry used in the current simulations. However, we note that earlier simulations with a similar setup did not show a marked trend in the energy of the differential rotation as the azimuthal extent of the domain was varied (see Table 1 of Käpylä et al. 2013). However, results of Boussinesq simulations of convective dynamos have shown a similar change as a function of the magnetic Prandtl number (Schrinner et al. 2012). The drop in the amplitude of the differential rotation was associated with a change in the dynamo mode from an oscillatory multipolar solution to a quasi-steady dipolar configuration (cf. Fig. 15 of Schrinner et al. 2012) that prevents strong differential rotation from developing. We do not find a strong dipole component in our simulations (see Sect. 4.2.3). However, the strong suppression of the differential rotation often coincides with the appearance of a small-scale dynamo (see Table 1 and the discussion in the Sect. 4.1.2) or a change in the large-scale dynamo mode as discussed above.

4.1.2. Angular momentum transport

The azimuthally averaged z-component of the angular momentum is governed by the equation ∂t(ρϖ2Ω)+·{ϖ[ϖρUΩ+ρuφu2νρS·φ̂\begin{eqnarray} % && \frac{\pd}{\pd t} ({\rho} \varpi^2 \Omega) + {\nabla\cdot}\{\varpi [\varpi{\rho {U}}\Omega + {\rho}\ {u_\phi {u}} - 2\nu{\rho} {{\mathsf{S}}}{{\cdot}}{{\hat\phi}} \nonumber \\ && \hspace{4cm} -\mu_0^{-1}({B}_\phi{{B}} + {b_\phi {b}})] \}=0, \label{equ:angmom} \end{eqnarray}(27)where ϖ = rsinθ is the lever arm, and velocity and magnetic field have been decomposed into mean and fluctuating parts according to Ui=Ui+ui\hbox{$U_i = \mean{U}_i + u_i$} and Bi=Bi+bi\hbox{$B_i = \mean{B}_i + b_i$}. We denote the Reynolds and Maxwell stresses as 𝒬ij=uiuj\hbox{$\mathcal{Q}_{ij}=\mean{u_i u_j}$} and ij=bibj/μ0ρ\hbox{$\mathcal{M}_{ij}=-\mean{b_i b_j}/\mu_0 \mean{\rho}$}, respectively. We note that Eq. (27) differs from the formulation of Brun et al. (2004), for instance, in that we have retained terms containing the mass flux ρU in the contributions corresponding to the meridional circulation. This is because in our fully compressible setup, ·ρU is non-zero. This is particularly important for the average radial mass flux ρUr\hbox{$\mean{\rho U_r}$}. However, it turns out that the effects of compressibility in the Reynolds stress and the viscous terms are negligible. This means that terms of the form (ρui)uj\hbox{$\mean{(\rho u_i)'u_j}$}, where the prime denotes fluctuations, ρuiuj\hbox{$\mean{\rho u_i u_j}$}, and νρui;j\hbox{$\nu\mean{\rho u_{i;j}}$}, are neglected.

thumbnail Fig. 4

From left to right: Reynolds, Maxwell, and total stress components \hbox{$\qqrp$}, , T, \hbox{$\qqtp$}, θφ, and Tθφ for Runs C2 (top), C3, C4, and C5 (bottom). Data nearer than 2.5° from the latitudinal boundaries are not shown, so as to emphasize the structures at lower latitudes. The yellow contours denote the zero levels in each panel.

The main generator of differential rotation is commonly thought to be the Reynolds stress and in particular its non-diffusive contribution that is due to the Λ effect (Rüdiger 1980, 1989). Distinguishing the contributions from the Λ effect and the turbulent viscosity is currently only possible using assumptions regarding either of the two coefficients and computing the other from the Reynolds stress (e.g., Käpylä et al. 2014; Karak et al. 2015; Warnecke et al. 2016a). We do not attempt this here, but study the total turbulent stress \hbox{$T_{ij}=\mathcal{Q}_{ij}+\mathcal{M}_{ij}$} and the anisotropy parameters realized in simulations with varying magnetic Reynolds number. Figure 4 shows the off-diagonal Reynolds stress components, \hbox{$\qqrp$} and \hbox{$\qqtp$}, from Runs C2 to C5 with ReM varying between 14 and 146. The spatial structures of the two stress components remain relatively similar in Runs C2, C3, and C4. Neither of the stresses shows a monotonous behavior as functions of ReM, both being generally smaller in Run C3 than in Runs C2 and C4. However, in Run C5 the overall magnitude is decreased by a factor of two for \hbox{$\qqrp$} and about 15 per cent for \hbox{$\qqtp$}. No corresponding decrease is observed in the diagonal components of the stress, whose overall magnitude is given by the fluid Reynolds number; see Table 1.

We find that for magnetic Reynolds numbers up to roughly 30, the off-diagonal Maxwell stresses are smaller than the corresponding Reynolds stresses. The magnitude of the vertical component is between a quarter and a half of \hbox{$\mathcal{Q}_{r\phi}$}, whereas for the horizontal component the difference is greater. With higher ReM, the Maxwell stresses attain similar profiles as the Reynolds stresses but with opposite signs. The magnitudes of the Maxwell stresses increase with ReM such that they become comparable to and locally even larger than the Reynolds stresses at the highest magnetic Reynolds numbers. The total vertical stress T decreases for Runs C2–C4 such that the effect is most clear near the equator. For Run C5 the vertical stress is dominated by the Maxwell stress near the equator. The horizontal component Tθφ decreases monotonically as ReM increases. The near cancellation of the total stress at high ReM is likely to contribute significantly to the quenching of differential rotation. These results are apparently at odds with those of Karak et al. (2015), who reported that the Maxwell stress is an order of magnitude smaller than the Reynolds stress. However, their simulations were in a regime near the transition between anti-solar and solar-like rotation profiles with relatively weak and intermittent large-scale dynamos.

In Set F with anti-solar differential rotation the Reynolds stresses show only a weak decreasing trend as a function of ReM with a corresponding increase in the Maxwell stress (not shown). The vertical Maxwell stress has an opposite sign in comparison to the Reynolds stress in all cases, while a similar tendency for the horizontal stress is not so clear. The magnitude of the vertical Maxwell stress is roughly half of the corresponding Reynolds stress, while the amplitude of the horizontal Maxwell stress is significantly weaker than the horizontal Reynolds stress. The much weaker quenching of the total turbulent stress in Set F is consistent with a clearly milder decrease of the differential rotation than in the other sets.

Assuming that the turbulent viscosity follows a naive mixing length estimate, the quenching of the differential rotation can indicate that the Λ effect is more severely quenched by the large-scale magnetic field at high magnetic Reynolds numbers. Another possibility is that small-scale magnetic fields generated by an efficient small-scale dynamo contribute to enhancing the turbulent viscosity. In a recent paper, Hotta et al. (2016) suggested that the small-scale dynamo at high ReM suppresses velocity at small scales and facilitates the growth of the large-scale magnetic fields. First-order smoothing estimates for isotropic and homogeneous turbulence (M. Rheinhardt, private communication) suggest that turbulent viscosity acquires a contribution νt(SSD)=13brms(0)kf-1\hbox{$\nut^{\rm (SSD)} = \onethird b_{\rm rms}^{(0)} \kf^{-1}$}, where brms(0)\hbox{$b_{\rm rms}^{(0)}$} is the rms-value of the saturated fluctuating magnetic field due to small-scale dynamo. We find that the fluctuating magnetic field energy grows monotonically as a function of ReM; see Table 2. A corresponding increase of turbulent viscosity would be compatible with a strong decrease of differential rotation. However, we find two counter-examples: in Set A, no SSD is present but the differential rotation still experiences strong quenching, and in Set F, an SSD is present but differential rotation remains strong. Furthermore, the dependence of the Λ effect on small-scale magnetic fields is currently unknown. Thus the question of the effect of small-scale dynamo on the turbulent transport of angular momentum at large scales remains open.

thumbnail Fig. 5

Anisotropy parameters AV (left) and AH (right) for the same runs as in Fig. 4. The yellow contours denote the zero levels in each panel.

We also show the anisotropy parameters (Fig. 5) AV=𝒬φφ𝒬rr𝒬φφ+𝒬rr,AH=𝒬φφ𝒬θθ𝒬φφ+𝒬θθ,\begin{eqnarray} A_{\rm V} &=& \frac{\qqpp-\qqrr}{\qqpp+\qqrr},\\ A_{\rm H} &=& \frac{\qqpp-\qqtt}{\qqpp+\qqtt}, \end{eqnarray}which are proportional to the vertical and horizontal Λ effects in mean-field hydrodynamics (Rüdiger 1980) under the assumption of slow rotation (see also the numerical results of Käpylä & Brandenburg 2008). Again the differences between Runs C2 and C3 are relatively small, whereas in Runs C4 and C5 the mostly positive AV at mid-latitudes in lower-ReM runs gives way to negative values. According to mean-field theory, this corresponds to a sign change of the vertical Λ effect that is responsible for generating radial differential rotation (Rüdiger 1980). The magnitude and the spatial distribution of the horizontal anisotropy parameter AH are significantly different in Run C5 in comparison to the other runs, whereas the differences between the other runs are less significant. The negative values of AH at mid-latitudes in Runs C2–C4 coincide with the minimum in the angular velocity, suggesting that the horizontal Λ effect is negative there. In Run C5 the overall magnitude of AH is diminished near the surface at low latitudes. Along with the change of sign of AV, this is likely to contribute to the reduced differential rotation as a function of ReM. However, more theoretical work is needed to distinguish the effects of the small- and large-scale magnetic fields on the angular momentum transport.

thumbnail Fig. 6

Radial (divFr, left panel) and latitudinal (divFθ, middle) parts of the divergence of the angular momentum fluxes. The right panel shows the total divergence. The units are given in the legend. Data taken from Run C3.

The discussion above is valid when the angular momentum is in a statistically steady state demanding that the sum of two terms, divFr + divFθ, vanishes, that is, 1r2(r2r)∂r+1rsinθ(sinθθ)∂θdivFr+divFθ=0,\begin{eqnarray} \frac{1}{r^2} \frac{\pd (r^2\mathcal{F}_r)}{\pd r} + \frac{1}{r\sin\theta}\frac{\pd (\sin\theta \mathcal{F}_\theta)}{\pd \theta}\equiv \mbox{div}F_r + \mbox{div}F_\theta = 0, \label{equ:divF} \end{eqnarray}(30)where r=ϖ[ϖρUrΩ+ρ𝒬νρϖΩ∂rμ0-1(BφBr+brbφ)],\begin{equation} \mathcal{F}_r\!=\!\varpi\bigg[\varpi \mean{\rho U_r}\Omega\!+\!\mean\rho\mathcal{Q}_{r\phi}\!-\!\nu \rho \varpi \frac{\pd \Omega}{\pd r}\!-\!\mu_0^{-1}(\mean{B}_\phi\mean{B}_r\!+\!\mean{b_r b_\phi})\bigg], \end{equation}(31)and θ=ϖ[ϖρUθΩ+ρ𝒬θφνρϖrΩ∂θμ0-1(BφBθ+bθbφ)].\begin{equation} \mathcal{F}_\theta\!=\!\varpi\bigg[\varpi \mean{\rho U_\theta}\Omega\!+\!\mean\rho\mathcal{Q}_{\theta\phi}\!-\!\nu \rho \frac{\varpi}{r} \frac{\pd \Omega}{\pd \theta}\!-\!\mu_0^{-1}(\mean{B}_\phi\mean{B}_\theta\!+\!\mean{b_\theta b_\phi})\bigg]. \end{equation}(32)Figure 6 shows representative results from Run C3 for the terms divFr and divFθ, as well as their sum. We find that the radial and latitudinal contributions have similar structures but opposite signs and that their sum is at most roughly eight per cent of the individual components. The elongated structures at low latitudes within the tangent cylinder are due to the meridional circulation that yields the dominant contribution to the divergence. We note that neither the radial nor the latitudinal fluxes need to individually cancel for the divergence to vanish.

4.2. Large-scale magnetic fields and dynamo cycles

4.2.1. General considerations

We find no dynamos in the lowest magnetic Reynolds number cases Runs A1, A2, B1, C1, and D1, which are in the range ReM = 5...9. On the other hand, the highest magnetic Reynolds numbers clearly exceed critical values for small-scale dynamo action to occur in simpler setups (e.g., Schekochihin et al. 2005). In the present simulations the small- and large-scale dynamos can be excited at the same time, and separating the two is not directly possible. We therefore resort to runs where we artificially suppress the large-scale (φ-averaged) magnetic fields at each time step. This eliminates the large-scale dynamo and growing magnetic fields can be associated with a small-scale dynamo. We have performed such runs for each of the cases where a dynamo is observed and find that a small-scale dynamo is excited in Runs B5, C5, D4, D5, E3, E4, F4, G2, and G3.

Table 3

Cycle lengths detected using D2 statistic.

4.2.2. Cycle detection using D2 analysis

Using the D2 statistic discussed in Sect. 3, we separately analyzed radial and azimuthal components of the mean magnetic field. In the former case we included all the data near the polar region (55° ≤ | 90°−θ | ≤ 75°) and in the latter case data around mid-latitude region (10° ≤ | 90°−θ | ≤ 45°). We considered data above 0.94 R and analyzed the northern and southern hemispheres separately.

thumbnail Fig. 7

Azimuthally averaged azimuthal magnetic field near the surface from Runs B2 (top left), D3 (top right), B3 (bottom left), and F3 (bottom right) with poleward (PW), equatorward (EW), irregular (IR), and quasi-stationary (QS) dynamo solutions, respectively.

The results of the analysis are listed in Table 3. We required that at least five full cycles were covered for each trial period in the period search range. To meet this criterion, we adjusted the upper limit of the period search range according to the length of the dataset, while the lower limit was always fixed at one year. We set the lower limit for the coherence length range to two cycles per given period and the upper value was determined by the dataset length. However, in some cases the longest possible trial period was still shorter than two years, in which cases the analysis was deemed infeasible. This situation was encountered for Run G3, which has therefore not been included in Table 3. Furthermore, Runs A1, A2, B1, C1, and D1 without dynamos were not analyzed.

All the cycle lengths given in Table 3 are significant with the p-values lower than 1%. More precisely, in our case the p-value represents a probability that a cycle with a given period would appear by chance out of white-noise data with the same distribution as the original data. In those cases where the obtained D2 spectrum contained no significant minima, no cycle was detected. However, here we must also note that as a result of narrowing the period search range in some cases, it is possible that the real cycle length is located outside the range and we did not detect it. This is particularly relevant for the runs at high magnetic Reynolds numbers where the time series are short.

In all the sets, especially when the highest Reynolds number cases are excluded, the cycle lengths are found to increase as functions of the magnetic Prandtl number. This shows that the dynamo period is sensitive to the strength of the magnetic diffusion such that when the diffusion is decreased, and correspondingly the diffusion timescale increases, the dynamo period increases.

4.2.3. Dynamo modes

Given that recent simulations reproduce solar-like magnetic cycles with equatorward migration (hereafter EW) of activity belts, it is of interest to probe the parameter space to determine when such solutions are excited. After identifying cyclic solutions with the D2 statistics, we classified them as EW or PW (poleward) based on the migration direction of the dynamo wave at low latitudes between ± 10° and ± 45°. We also find quasi-stationary (hereafter QS) solutions in Set F. When a clear classification by visual inspection could not be made, we classify the solution as irregular (IR). In some cases features from more than one class can be present, for instance, an equatorward cyclic variation in addition to a quasi-stationary background. This solution is classified as EW/QS, but typically such cases are quite uncertain and are indicated by brackets in the last column of Table 3. No large-scale dynamo action is denoted by ND. In Fig. 7 we show time-latitude plots of the azimuthally averaged azimuthal magnetic field Bφ\hbox{$\mean{B}_\phi$} from four runs exemplifying each of the main dynamo modes discussed above.

We first consider the Sets A–D with Ω0 = 5Ω. For PrSGS = 0.25 (Set A) and moderate ReM we obtain solutions with clear poleward migration. This type of dynamos was first obtained in the pioneering studies of Gilman (1983) and Glatzmaier (1985). Solutions showing poleward migration have been reported more recently by many groups (e.g., Käpylä et al. 2010b; Brown et al. 2011; Schrinner et al. 2012). In one case (Run A3) we see a hemispheric dynamo with poleward migration, similar to those reported by Busse (2002) and Gastine et al. (2012), for example. In the largest ReM case (Run A5) the clear oscillatory solutions of the runs with smaller Reynolds numbers give way to possibly irregularly reversing or quasi-stationary large-scale fields. The time series is too short, however, for possible cycles to be detected. For PrSGS = 0.5 (Set B) we find a PW mode excited in the case ReM = 12 (Run B2). At intermediate ReM (Runs B3 and B4) the solutions mostly show irregularly reversing fields, although features from both PW and EW modes can be discerned at times; see the lower left panel of Fig. 7. The highest Reynolds number case (Run B5) appears to return to a poleward oscillatory mode with a longer cycle period. However, the large-scale field shows only four reversals in the 24-yr duration of the simulation and the D2 statistic captures this only in one of the four analyzed cases. Therefore the classification of this run as PW is deemed less reliable than those of the other runs in this set.

thumbnail Fig. 8

Azimuthally averaged azimuthal magnetic field Bφ\hbox{$\mean{B}_\phi$} (color contours) in units of kG and the fields lines of the poloidal field (continuous and broken lines for clockwise and anticlockwise loops, respectively). Left: data from Run B2 averaged over three months near a cycle maximum at t = 70 yr. Right: data from Run A5 averaged over the last ten years of the run.

For PrSGS ≥ 1 (Sets C and D) we often find solutions with EW migration (Runs C2, C3, D2, D3, and D4). In Run D4 the solution is clearly EW, but the time series covers only three full cycles and is therefore not detected by D2. The transition from poleward migrating to equatorward migrating solutions at intermediate magnetic Reynolds numbers coincides with the change of the rotation profile from runs with consistently positive radial gradient of Ω to ones with a minimum of Ω at mid-latitudes where Ω /∂r< 0 as PrSGS increases. The change in the dynamo mode fits with the interpretation in terms of a classical dynamo wave obeying the Parker-Yoshimura rule (Warnecke et al. 2014, 2016b). For the highest magnetic Reynolds number cases in Sets C and D we again observe possibly quasi-stationary or irregular configurations. This classification is based on significantly shorter time series than in the other cases, and it is possible that there are cycles that are much longer than in the low-ReM cases or that a prolonged transient is still in progress. These results appear to be in agreement with those of Schrinner et al. (2012), who found a transition from oscillatory dynamos to non-oscillatory ones at a roughly comparable ReM for Boussinesq convection in spherical shells. They also showed that in the high-PrM regime the two bistable branches of dynamo solutions merge and that only strongly dipolar dynamos with weak differential rotation survive. However, Gastine et al. (2012) have shown that for sufficiently density stratified cases the dipolar branch does not exist at least for moderate values of ReM. Our simulations with moderate density stratification of Hρ ≈ 3, where Hρ = −(lnρ/∂r)-1 is the density scale height, do not indicate that a dipolar mode takes over at high ReM, see Fig. 8 for a representative result from Run A5. This is consistent with the results of Gastine et al. (2012), who found no dipole-dominated dynamos for Hρ ≳ 2.

The realized solutions from the runs in Sets A–D are plotted in the (PrSGS,PrM)-plane in Fig. 9. We find that the low and intermediate ReM runs produce cyclic solutions with PW for low PrSGS and EW for high PrSGS. The case PrSGS = 0.5 works as a watershed between the two cyclic regimes. The dynamos at high ReM are fundamentally different from their more diffusive counterparts in that the differential rotation is almost absent. Although the ratio of the energies in the mean poloidal to toroidal components is not significantly different in the high-ReM runs in comparison to lower-ReM runs in each set, the dynamos at the high-ReM regime can be of α2 type. Substantiating this claim, however, requires that the turbulent transport coefficients relevant for the maintenance of large-scale magnetic fields are extracted from the simulations and applied in corresponding mean-field models, which is not within the scope of the present study.

thumbnail Fig. 9

Color contours: radial derivative of Ω at r = 0.85 R averaged from latitudes + 25° and −25°. The dynamo modes realized in Sets A–D are overplotted with abbreviations ND, PW, EW, IR, and QS. The thick yellow curve indicates the zero level of Ω /∂r. Overlapping symbols denote solutions that show characteristics from several modes, while question marks indicate that the classification is uncertain because of the insufficient length of the data set.

thumbnail Fig. 10

Energy densities of the total magnetic field (left panel), and the azimuthally averaged fields (right) from Sets C (black dashed line) and G (black solid), from two sets of runs by Nelson et al. (2013; solid and dashed red lines) and from Hotta et al. (2016; blue solid lines). The red solid line consists of data from cases D3, D3a, and D3b of Nelson et al. (2013) with PrM = const. = 0.5. Correspondingly, the dashed red line shows data from cases D3, D3-pm1, and D3-pm2 of Nelson et al. (2013) with PrM varying from 0.5 to 2. From Hotta et al. (2016) we show data for cases “Low”, “Medium”, and “High”.

In Set E we find a PW solution for the lowest ReM (Run E1), whereas at higher ReM we find IR (E2) and irregular/quasi-stationary (E3 and E4) configurations. Generally, Set E behaves similarly as Sets A–C, that is, a transition from oscillatory dynamos at intermediate ReM to quasi-stationary or irregularly varying solutions at high ReM. However, Run B1 of Warnecke et al. (2016a) is similar to Run E2, but with PrSGS = 2 instead of 1 shows EW. Set F is qualitatively different from the other sets similarly as for the differential rotation. The large-scale magnetic fields show a quasi-stationary configuration in all of the runs in Set F. In Run G2 the modulation of the differential rotation was shown to be related to a changing large-scale dynamo mode that is either EW or QS; see panels b and d of Fig. 3. In the highest-ReM case, Run G3, the two competing modes appear to be present again, but the short data set length renders such classifications preliminary at best.

4.2.4. Saturation level of large-scale magnetic fields

For the total magnetic energy we find a monotonically increasing trend as a function of ReM in all sets, see representative results in the left panel of Fig. 10 and Table 2. The absolute value of the axisymmetric parts of the poloidal and toroidal magnetic fields shows monotonically decreasing trends only in Set F with anti-solar differential rotation. In the other sets the energies of the poloidal and toroidal mean fields typically do not behave monotonically as functions of ReM. However, given the large temporal variations, our results are compatible with mean magnetic field energies converging to constant values at high magnetic Reynolds numbers, see the right panel of Fig. 10 for the results of Sets C and G. The same conclusion also applies to Set F. The findings for the mean magnetic fields appear to contradict the results of Nelson et al. (2013), who found a monotonically decreasing trend as a function of ReM; see Fig. 10.

The main difference to the simulations of Nelson et al. (2013) is that their models were done with a full spherical shell as opposed to the wedge geometry used here, and that magnetic field boundary condition at the outer radial boundary uses extrapolation to a potential field rather than a radial field condition as in the present study. Furthermore, Nelson et al. (2013) considered anelastic models, whereas in our case the gas is fully compressible. Simulations with forced turbulence in spherical shells with coronal envelopes have shown that a higher field strength can be achieved when the magnetic field at the boundary is not restricted to being purely radial (Warnecke & Brandenburg 2014). Furthermore, Nelson et al. (2013) used profiles for ν and η, which are absent in our study. The constant diffusion coefficients used here may lead to a steeper increase of the magnetic field energy than in the runs of Nelson et al. (2013), where the local value of ReM is significantly higher in the deeper layers. Detailed comparisons of diffusion schemes and parameter values are presented in Tables A.1 and A.2.

It is difficult to assess which of these differences is the most important. However, a plausible candidate is the change in topology of the field in the full-sphere simulations of Nelson et al. (2013). As described in Käpylä et al. (2013), the large-scale magnetic field can become dominated by low-order non-axisymmetric modes, and therefore applying an axisymmetric mean will average out such field contributions (see their Fig. 16 and Table 2). Evidence of such non-axisymmetric modes can be seen in the instantaneous magnetic fields in Figs. 4 and 6 of Nelson et al. (2013). However, it is not possible to assess how the degree of non-axisymmetry behaves as a function of ReM in the results of Nelson et al. (2013), and a direct comparison is therefore not possible.

In another recent study, Hotta et al. (2016) presented results from less rapidly rotating (Ω = Ω) simulations at high values of ReM. Such simulations are less likely to contain significant non-axisymmetric modes. Their main claim is that whereas the mean magnetic energy decreases at intermediate ReM, it recovers as ReM is increased further. We have included their results in Fig. 10 for comparison. However, only a rough comparison is possible as only three data points are available and a jump occurs between the two lowest-ReM runs, which is caused by a switch from explicit diffusion to a numerical slope-limited scheme. The absolute values for the mean fields are also lower than ours, probably because their large-scale dynamo is less efficient than in the present study because of their slower rotation. Only the trend as a function of ReM can therefore be compared with the current results or with those of Nelson et al. (2013). The results of Hotta et al. (2016) are more in line with ours, but the degree of temporal variations is indicated only in passing. Hotta et al. (2016) stated that during the last 200 days the mean magnetic field energy in the simulation high is almost 2.5 times higher than in the time average over the full duration (50 yr) of the simulation. It is not possible to accurately assess the steepness of the increasing trend of the energy of the mean field as a function of ReM.

5. Conclusions

We find that, as the SGS Prandtl number (responsible for SGS turbulent transport) is increased, the rotation profiles realized in the simulations develop a region of negative shear at mid-latitudes. At moderate ReM this latitude coincides with a transition from poleward-migrating to equatorward-migrating dynamo modes. This can be explained by interpreting the solutions as dynamo waves propagating along the isocontours of constant shear (Parker 1955; Yoshimura 1975). However, it appears that, as PrM is sufficiently high (corresponding to high ReM), the regular cycles give way to quasi-stationary or irregularly varying solutions. However, as a result of computational constraints, the time series of these runs are typically significantly shorter than of the lower-ReM runs, so that cyclic solutions with long cycle periods cannot be ruled out.

We also find that the cycle length of the dynamo solution undergoes a systematic increase when the magnetic Prandtl number is increased: seemingly independent of the other parameters, the decreasing magnetic diffusivity leads to longer cycles. The highest Reynolds number runs are, however, too short for our period analysis to work conclusively.

We find a strong dependence of the differential rotation on the magnetic Reynolds number so that for the highest values of ReM both radial and latitudinal shear are almost absent. The strongest quenching tends to appear in cases where a small-scale dynamo is present. However, there are exceptions. The physical reason for the quenching is therefore unclear, but several mechanisms appear plausible. First, the small-scale magnetic field can enhance turbulent viscosity or quench the Λ-effect responsible for maintaining differential rotation. Suppression of small-scale flows that are due to the small-scale magnetic fields has been suggested by Hotta et al. (2016). We find that at intermediate ReM, the Maxwell stress is on the same order of magnitude as the Reynolds stress and becomes comparable at high ReM. Because of their opposite signs, the total stress is diminished. We expect that for intermediate and in particular for high ReM, the Maxwell stress plays an important role in the angular momentum transport. Second, the dependence of turbulent transport on the large-scale magnetic field is ReM-dependent and can lead to enhanced quenching in the parameter regime studied here in comparison to earlier more laminar simulations. However, for a limited parameter range we find that the system vacillates between two states where either the differential rotation is strong and mean magnetic field relatively weak, or vice versa. Furthermore, the temporal variations appear to increase as the Rayleigh and Reynolds numbers are increased. Finally, we find that in cases where the differential rotation is anti-solar, the quenching is much less prominent, possibly due to significantly weaker mean fields generated in those cases.

The total magnetic energy grows monotonically as a function of ReM in all of our runs. The energy of the azimuthally averaged fields is in all cases consistent with increasing or constant mean fields at high ReM. This is apparently at odds with the anelastic full-sphere simulations of Nelson et al. (2013), who found a steeply declining trend with magnetic Reynolds number. However, their relatively rapidly rotating (Ω = 3Ω) simulations also allow significant low-order non-axisymmetric modes to develop that will not show up in azimuthal averaging. A fair comparison is therefore not possible without a proper assessment of the non-axisymmetric contributions. Our results are more in line with those of Hotta et al. (2016), who used a model rotating at the solar rate where the large-scale fields are more clearly axisymmetric.

A possible source of discrepancies between the current and previous studies is the use of wedge geometry. Rigorous comparisons between wedges and fully spherical simulations have not been done so far. However, runs with similar parameters (Rayleigh, Taylor, and Prandtl numbers) produce results that are in qualitative agreement; compare the results of Käpylä et al. (2010b) regarding cyclic dynamos with those of Brown et al. (2011), for instance. Similarly, the more recent simulations showing equatorward migration (Käpylä et al. 2012; Augustson et al. 2015) appear to support the validity of the wedge approach. Furthermore, wedges with full 2π extent in longitude produce dynamos dominated by non-axisymmetric large-scale fields with azimuthal dynamo waves (Käpylä et al. 2013; Cole et al. 2014) that are also routinely seen in anelastic full-sphere simulations (e.g., Yadav et al. 2015). Last, the transition from anti-solar to solar-like differential rotation occurs at very similar Coriolis numbers in a wide range of simulations, including wedges with artificially increased luminosity and rotation rate (Gastine et al. 2014). The various modeling approaches also use a range of different SGS models and parameter values (see Tables A.1 and A.2) and are still able to reproduce similar large-scale phenomena. These comparisons suggest that the current results obtained in wedges are likely to be realized in fully spherical simulations in the same parameter regime.

Our final conclusion is that the current simulations are not near an asymptotic regime where the large-scale results would be independent of the microphysical diffusion coefficients. This is most strikingly demonstrated by the steep quenching of the differential rotation and the disappearance of regularly oscillating large-scale magnetic fields at high values of ReM.


Acknowledgments

The authors thank an anonymous referee for useful comments and Paul Charbonneau and Hideyuki Hotta for providing information regarding their models. 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 and in the HLRS supercomputing center in Stuttgart, Germany through the PRACE allocation “SOLDYN”. Financial support from the Academy of Finland grants No. 136189, 140970, 272786 (P.J.K.), and 272157 to the ReSoLVE Centre of Excellence (P.J.K., M.J.K.), as well as the Swedish Research Council grants 621-2011-5076 and 2012-5797, and the Research Council of Norway under the FRINATEK grant 231444 are acknowledged. J.W. acknowledges funding by the Max-Planck/Princeton Center for Plasma Physics and funding 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

  1. Augustson, K., Brun, A. S., Miesch, M., & Toomre, J. 2015, ApJ, 809, 149 [Google Scholar]
  2. Blackman, E. G., & Brandenburg, A. 2003, ApJ, 584, L99 [NASA ADS] [CrossRef] [Google Scholar]
  3. Brandenburg, A. 2001, ApJ, 550, 824 [NASA ADS] [CrossRef] [Google Scholar]
  4. Brandenburg, A., & Dobler, W. 2001, A&A, 369, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Brandenburg, A., & Sandin, C. 2004, A&A, 427, 13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Brandenburg, A., Chan, K. L., Nordlund, Å., & Stein, R. F. 2005, Astron. Nachr., 326, 681 [NASA ADS] [CrossRef] [Google Scholar]
  7. Brandenburg, A., Candelaresi, S., & Chatterjee, P. 2009, MNRAS, 398, 1414 [NASA ADS] [CrossRef] [Google Scholar]
  8. Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424 [Google Scholar]
  9. Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 2011, ApJ, 731, 69 [NASA ADS] [CrossRef] [Google Scholar]
  10. Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [NASA ADS] [CrossRef] [Google Scholar]
  11. Busse, F. H. 2002, Phys. Fluids, 14, 1301 [NASA ADS] [CrossRef] [Google Scholar]
  12. Busse, F. H., & Simitev, R. D. 2006, Geophys. Astrophys. Fluid Dynam., 100, 341 [CrossRef] [Google Scholar]
  13. Cattaneo, F. 1999, ApJ, 515, L39 [NASA ADS] [CrossRef] [Google Scholar]
  14. Cattaneo, F., & Hughes, D. W. 1996, Phys. Rev. E, 54, 4532 [NASA ADS] [CrossRef] [Google Scholar]
  15. Cattaneo, F., & Vainshtein, S. I. 1991, ApJ, 376, L21 [NASA ADS] [CrossRef] [Google Scholar]
  16. Cole, E., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2014, ApJ, 780, L22 [NASA ADS] [CrossRef] [Google Scholar]
  17. Cole, E., Brandenburg, A., Käpylä, P. J., & Käpylä, M. J. 2016, A&A, 593, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Del Sordo, F., Guerrero, G., & Brandenburg, A. 2013, MNRAS, 429, 1686 [NASA ADS] [CrossRef] [Google Scholar]
  19. Duarte, L. D. V., Wicht, J., Browning, M. K., & Gastine, T. 2016, MNRAS, 456, 1708 [NASA ADS] [CrossRef] [Google Scholar]
  20. Gastine, T., Duarte, L., & Wicht, J. 2012, A&A, 546, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76 [NASA ADS] [CrossRef] [Google Scholar]
  22. Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gilman, P. A. 1983, ApJS, 53, 243 [NASA ADS] [CrossRef] [Google Scholar]
  24. Glatzmaier, G. A. 1985, ApJ, 291, 300 [NASA ADS] [CrossRef] [Google Scholar]
  25. Gruzinov, A. V., & Diamond, P. H. 1995, Phys. Plasmas, 2, 1941 [NASA ADS] [CrossRef] [Google Scholar]
  26. Hotta, H., Rempel, M., Yokoyama, T., Iida, Y., & Fan, Y. 2012, A&A, 539, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Hotta, H., Rempel, M., & Yokoyama, T. 2016, Science, 351, 1427 [NASA ADS] [CrossRef] [Google Scholar]
  28. Jennings, R., Brandenburg, A., Tuominen, I., & Moss, D. 1990, A&A, 230, 463 [NASA ADS] [Google Scholar]
  29. Käpylä, P. J., & Brandenburg, A. 2008, A&A, 488, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2010a, A&A, 518, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Käpylä, P. J., Korpi, M. J., Brandenburg, A., Mitra, D., & Tavakol, R. 2010b, Astron. Nachr., 331, 73 [NASA ADS] [CrossRef] [Google Scholar]
  32. Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011a, Astron. Nachr., 332, 883 [NASA ADS] [CrossRef] [Google Scholar]
  33. Käpylä, P. J., Mantere, M. J., Guerrero, G., Brandenburg, A., & Chatterjee, P. 2011b, A&A, 531, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJ, 755, L22 [NASA ADS] [CrossRef] [Google Scholar]
  35. Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [NASA ADS] [CrossRef] [Google Scholar]
  36. Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Käpylä, M. J., Käpylä, P. J., Olspert, N., et al. 2016a, A&A, 589, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Käpylä, P. J., Brandenburg, A., Kleeorin, N., Käpylä, M. J., & Rogachevskii, I. 2016b, A&A, 588, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Karak, B. B., Käpylä, P. J., Käpylä, M. J., et al. 2015, A&A, 576, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Lindborg, M., Mantere, M. J., Olspert, N., et al. 2013, A&A, 559, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Masada, Y. 2011, MNRAS, 411, L26 [NASA ADS] [CrossRef] [Google Scholar]
  42. Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D. 2009, ApJ, 697, 923 [NASA ADS] [CrossRef] [Google Scholar]
  43. Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2013, ApJ, 762, 73 [NASA ADS] [CrossRef] [Google Scholar]
  44. Olspert, N., Käpylä, M. J., Pelt, J., et al. 2015, A&A, 577, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Parfrey, K. P., & Menou, K. 2007, ApJ, 667, L207 [NASA ADS] [CrossRef] [Google Scholar]
  46. Parker, E. N. 1955, ApJ, 121, 491 [NASA ADS] [CrossRef] [Google Scholar]
  47. Passos, D., & Charbonneau, P. 2014, A&A, 568, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Pelt, J. 1983, in Statistical Methods in Astronomy, ed. E. J. Rolfe, ESA SP, 201, 37 [Google Scholar]
  49. Rüdiger, G. 1980, Geophys. Astrophys. Fluid Dynam., 16, 239 [Google Scholar]
  50. Rüdiger, G. 1989, Differential Rotation and Stellar Convection. Sun and Solar-type Stars (Berlin: Akademie Verlag) [Google Scholar]
  51. Schekochihin, A. A., Haugen, N. E. L., Brandenburg, A., et al. 2005, ApJ, 625, L115 [NASA ADS] [CrossRef] [Google Scholar]
  52. Schrinner, M., Petitdemange, L., & Dormy, E. 2011, A&A, 530, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Schrinner, M., Petitdemange, L., & Dormy, E. 2012, ApJ, 752, 121 [NASA ADS] [CrossRef] [Google Scholar]
  54. Shukurov, A., Sokoloff, D., Subramanian, K., & Brandenburg, A. 2006, A&A, 448, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Simitev, R., & Busse, F. H. 2005, J. Fluid Mech., 532, 365 [NASA ADS] [CrossRef] [Google Scholar]
  56. Smolarkiewicz, P. K., & Charbonneau, P. 2013, J. Comp. Phys., 236, 608 [CrossRef] [Google Scholar]
  57. Stix, M. 2002, The Sun: An Introduction (Berlin: Springer) [Google Scholar]
  58. Strugarek, A., Beaudoin, P., Brun, A. S., et al. 2016, Adv. Space Res., 58, 1538 [NASA ADS] [CrossRef] [Google Scholar]
  59. Sur, S., Shukurov, A., & Subramanian, K. 2007, MNRAS, 377, 874 [NASA ADS] [CrossRef] [Google Scholar]
  60. Vainshtein, S. I., & Cattaneo, F. 1992, ApJ, 393, 165 [NASA ADS] [CrossRef] [Google Scholar]
  61. Warnecke, J., & Brandenburg, A. 2014, in IAU Symp. 302, eds. P. Petit, M. Jardine, & H. C. Spruit, 134 [Google Scholar]
  62. Warnecke, J., Brandenburg, A., & Mitra, D. 2011, A&A, 534, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Warnecke, J., Brandenburg, A., & Mitra, D. 2012, J. Spa. Weather Spa. Clim., 2, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, ApJ, 796, L12 [NASA ADS] [CrossRef] [Google Scholar]
  65. Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2016a, A&A, 596, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Warnecke, J., Rheinhardt, M., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2016b, A&A, submitted [arXiv:1601.03730] [Google Scholar]
  67. Yadav, R. K., Gastine, T., Christensen, U. R., & Reiners, A. 2015, A&A, 573, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Yoshimura, H. 1975, ApJ, 201, 740 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Comparison to other simulation methods

Table A.1

Diffusion schemes applied in a few comparable studies.

Table A.2

Parameter values with conversion factors to the definition in the present study from the same studies as in Table A.1.

The purpose of this appendix is to compare the diffusion schemes and estimates of the Prandtl numbers in the present study with several methods presented in the literature. We consider the papers by Nelson et al. (2013) and Hotta et al. (2016) that are discussed in the main text as well as the study of Ghizaru et al. (2010), which represents another established simulation method; see Tables A.1 and A.2 for details.

Estimating the Prandtl number from simulations using the solar luminosity requires that the value of the radiative diffusivity in the solar convection zone is known. For the following we use an estimate of this quantity at r = 0.85 R. Standard solar models (e.g., Stix 2002) indicate that radiation carries roughly ten per cent of the flux at this radius, such that Lrad ≈ 0.1 L. Using Eq. (6), we can write Lrad=4πr2Frad=4πr2KrT.\appendix \setcounter{section}{1} \begin{equation} L^{\rm rad} = 4\pi r^2 F^{\rm rad} = - 4\pi r^2 K \pd_r T. \end{equation}(A.1)We then use the definition of the heat conductivity K = cPρχ, assume that the temperature gradient is close to the adiabatic one, rT ≈ −g/cP, and insert g = GM/r2 to obtain χ=Lrad4πGMρ·\appendix \setcounter{section}{1} \begin{equation} \chi = \frac{L^{\rm rad}}{4\pi G M_\odot \rho}\cdot \label{equ:radchi} \end{equation}(A.2)Using the values quoted above and ρ ≈ 50 kg m-3 (Stix 2002), we find that χ ≈ 5 × 102 m2 s-1 for the Sun at r = 0.85 R.

The solar luminosity is adopted in the study of Nelson et al. (2013), yielding a radiative diffusion coefficient on the order of the estimate from Eq. (A.2). The authors did not present a detailed description of their model, however, but referred to Brown et al. (2010), where the corresponding quantity is denoted as κr. Furthermore, the values of ν at mid-convection zone depth are in the range 6...13 × 107 m2 s-1 for their cases D3[a, b, pm1, pm2], and a value of 2 × 106 m2 s-1 is estimated for case S3 with a Smagorinsky SGS model. These yield Prandtl numbers on the order of 105 and 104, respectively. For the model “Low” of Hotta et al. (2016), only the surface value ν = 108 m2 s-1 is given. As ν is proportional to ρ− 1 / 2, the value at mid-convection zone depth is a few times 107 m2 s-1 and the Prandtl number is of the order of 105. For their cases “Medium”, “High”, and “High-S”, a slope-limited diffusion scheme is used, so that estimating the diffusion coefficients is non-trivial. However, they computed the Reynolds number (hereafter ReSL) based on the Taylor microscale that was obtained from kinetic energy spectra. They found that the Reynolds numbers for “High” and “High-S” are roughly an order of magnitude larger than for “Low”, whereas for “Low” and “Medium” they are comparable. Assuming that the turbulent velocities in all cases are of similar strengths, we can estimate “slope-limited diffusion viscosities” from νSL = νLowReLow/ ReSL. We therefore infer that the value of the Prandtl number in the model “Medium” is 105 and an order of magnitude smaller in the models “High” and “High-S”. In these runs, the surface value of the turbulent heat conductivity is κ = 2 × 109 m2 s-1 (H. Hotta, private communication), yielding PrSGS in the range 0.05...2 × 10-3 for the runs of Hotta et al. (2016).

In principle, even a third Prandtl number can be defined based on the diffusion coefficient applied only to the mean (spherically symmetric) entropy profile. In the present study and also in that of Hotta et al. (2016), this coincides with the coefficient relevant for diffusing entropy fluctuations. However, the coefficients can also be different, which is the formulation often used in simulations performed with the anelastic spherical harmonic (ASH) code including those presented in Nelson et al. (2013). However, the value of the coefficient for mean entropy diffusion (κ0) is not provided in the reference with a more detailed model description (Brown et al. 2010).

The Eulerian-Lagrangian (EULAG) code of Smolarkiewicz & Charbonneau (2013) employed by Ghizaru et al. (2010) uses quite a different approach and replaces radiative conduction by Newtonian cooling toward a prescribed thermodynamic state. The timescale of the cooling τc is typically on the order of 20 months (P. Charbonneau, private communication; see also Strugarek et al. 2016). This timescale would correspond to a radiative/SGS diffusion coefficient of χ = Δr2/τc ≈ 6 × 108 m2 s-1, where Δr = 0.25 R. Although this comparison yields some insight about the entropy evolution, the cooling and diffusion processes cannot be directly equated. Therefore the concept of a Prandtl number does not appear suitable in that case. The diffusion of velocity, magnetic fields, and entropy fluctuations is due to the numerical scheme making the estimates of the other Prandtl numbers also problematic. Presumably the diffusivities of all variables at a given resolution are roughly similar such that PrSGS and PrM are on the order of unity. Hydrodynamic EULAG models with an otherwise similar setup to that of Ghizaru et al. (2010) yield estimates of viscosity and entropy diffusion in ranges νeff = 0.6...1.2 × 108 m2 s-1 and κeff = 1...8 × 107 m2 s-1, with PrSGS = νeff/κeff ≈ 1...8 (Strugarek et al. 2016).

We show the conversion factors for the definitions used in the current study for Péclet, fluid and magnetic Reynolds, and Coriolis numbers in Table A.2. The conversion for the results of Nelson et al. (2013) is straightforward, except that for the Coriolis number we have used the relation Ta = Co′2Re′2 to compute Co′ = 2πCo. The definition of the Rossby number in Nelson et al. (2013), Ro = ω/ 2Ω, is based on the vorticity ω = × u and corresponds to Co = Ro-1kω/kf, where kω = ωrms/urms. The values of Co for the models of Nelson et al. (2013) correspond well to our simulations with Ω = 3Ω, see Table 1. In the case of the EULAG simulations we use the values quoted by Passos & Charbonneau (2014) for the Reynolds numbers (30...60 in their notation) divided by 2π and the estimate of PrSGS quoted above to compute Pe. The estimate of the Coriolis number is based on simulation data provided by P. Charbonneau. Finally, for the study of Hotta et al. (2016), the conversion of the fluid and magnetic Reynolds numbers is based on values given in the paper, whereas the estimates of Pe and Co are based on simulation data provided by H. Hotta. The values of Co are again in excellent agreement with our corresponding simulations with Ω = Ω; see Table 1.

We note that the simulations of Nelson et al. (2013) and Hotta et al. (2016) operate in a low-PrSGS regime that is also realized in the Sun, although in a much more extreme fashion. Such a regime is required especially in simulations with solar luminosity and rotation rate to lower the convective velocities and to achieve solar-like differential rotation (Käpylä et al. 2014; Hotta et al. 2016). However, the tradeoff is that the Péclet numbers are low and the evolution of entropy is significantly influenced by the SGS diffusion. It is not clear whether such an approach is more realistic than having PrSGS and PrM on the order of unity and Re, ReM, and Pe≫ 1. We also note that the run times of the highest resolution runs are short: four years for “S3” in Nelson et al. (2013) and 500 days for “High-S” in Hotta et al. (2016), whereas transients and/or secular evolution in the high-ReM regime can have a significantly longer timescale; see Figs 3b–d.

All Tables

Table 1

Summary of the runs.

Table 2

Volume and time-averaged kinetic and magnetic energy densities realized in the simulations in units of 105 J m-3.

Table 3

Cycle lengths detected using D2 statistic.

Table A.1

Diffusion schemes applied in a few comparable studies.

Table A.2

Parameter values with conversion factors to the definition in the present study from the same studies as in Table A.1.

All Figures

thumbnail Fig. 1

Top row: temporally averaged rotation profiles from Runs A3, B3, C3, and D3 with PrM = 1 and PrSGS varying from 0.25 (left) to 2 (right). Lower row: same as above, but for Runs C2, C3, C4, and C5 with PrSGS = 1 and PrM varying from 0.5 (left) to 5 (right).

In the text
thumbnail Fig. 2

Estimates of radial and latitudinal differential rotation ΔΩ(r)\hbox{$\Delta_\Omega^{(r)}$} and ΔΩ(θ)\hbox{$\Delta_\Omega^{(\theta)}$}, respectively, according to Eq. (26) for Sets A–F as indicated by the legend.

In the text
thumbnail Fig. 3

Energies of differential rotation (black solid lines) and mean toroidal magnetic field (red dashed) as functions of time from Runs G1, G2, and G3 (panels a)c)). Panel d) shows the azimuthally averaged azimuthal magnetic field at r = 0.98 R from Run G2. The vertical dashed lines indicate the high and low states of differential rotation. Panels e) and f) show the time-averaged rotation profiles in Run G2 from the high and low states indicated as blue dotted and dashed lines in panel b).

In the text
thumbnail Fig. 4

From left to right: Reynolds, Maxwell, and total stress components \hbox{$\qqrp$}, , T, \hbox{$\qqtp$}, θφ, and Tθφ for Runs C2 (top), C3, C4, and C5 (bottom). Data nearer than 2.5° from the latitudinal boundaries are not shown, so as to emphasize the structures at lower latitudes. The yellow contours denote the zero levels in each panel.

In the text
thumbnail Fig. 5

Anisotropy parameters AV (left) and AH (right) for the same runs as in Fig. 4. The yellow contours denote the zero levels in each panel.

In the text
thumbnail Fig. 6

Radial (divFr, left panel) and latitudinal (divFθ, middle) parts of the divergence of the angular momentum fluxes. The right panel shows the total divergence. The units are given in the legend. Data taken from Run C3.

In the text
thumbnail Fig. 7

Azimuthally averaged azimuthal magnetic field near the surface from Runs B2 (top left), D3 (top right), B3 (bottom left), and F3 (bottom right) with poleward (PW), equatorward (EW), irregular (IR), and quasi-stationary (QS) dynamo solutions, respectively.

In the text
thumbnail Fig. 8

Azimuthally averaged azimuthal magnetic field Bφ\hbox{$\mean{B}_\phi$} (color contours) in units of kG and the fields lines of the poloidal field (continuous and broken lines for clockwise and anticlockwise loops, respectively). Left: data from Run B2 averaged over three months near a cycle maximum at t = 70 yr. Right: data from Run A5 averaged over the last ten years of the run.

In the text
thumbnail Fig. 9

Color contours: radial derivative of Ω at r = 0.85 R averaged from latitudes + 25° and −25°. The dynamo modes realized in Sets A–D are overplotted with abbreviations ND, PW, EW, IR, and QS. The thick yellow curve indicates the zero level of Ω /∂r. Overlapping symbols denote solutions that show characteristics from several modes, while question marks indicate that the classification is uncertain because of the insufficient length of the data set.

In the text
thumbnail Fig. 10

Energy densities of the total magnetic field (left panel), and the azimuthally averaged fields (right) from Sets C (black dashed line) and G (black solid), from two sets of runs by Nelson et al. (2013; solid and dashed red lines) and from Hotta et al. (2016; blue solid lines). The red solid line consists of data from cases D3, D3a, and D3b of Nelson et al. (2013) with PrM = const. = 0.5. Correspondingly, the dashed red line shows data from cases D3, D3-pm1, and D3-pm2 of Nelson et al. (2013) with PrM varying from 0.5 to 2. From Hotta et al. (2016) we show data for cases “Low”, “Medium”, and “High”.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.