© ESO, 2015

1. Introduction

Differential rotation is an important ingredient for the generation of stellar magnetic fields. The internal rotation rate of the Sun has been mapped by helioseismology, revealing that the angular velocity within the convection zone mildly increases (decreases) as a function of radius at low (high) latitudes and that the radial shear is concentrated in shallow layers at the base of the convection zone and near the surface (e.g., Brown et al. 1989; Schou et al. 1998; Thompson et al. 2003). Thus the solar equator rotates faster than its poles. This kind of rotation profile is called solar-like (SL) differential rotation. The opposite case where the equator rotates slower than the poles, is referred to as anti-solar (AS) differential rotation. Because of the difficulties in observing slowly rotating stars that might possess AS differential rotation, it is not clear how common it is in main-sequence stars. However, it has been observed in some K giants (e.g., Strassmeier et al. 2003; Weber et al. 2005; Kővári et al. 2015).

Historically, the differential rotation and magnetic fields of the Sun and other stars have been modelled following two approaches – mean-field models and global convection simulations. In the mean-field approach, small-scale turbulence is parameterised by expressing the Reynolds stress in the momentum equation in terms of the mean velocity, the turbulent electromotive force in the induction equation in terms of the mean magnetic field, and the turbulent heat flux in the entropy equation in terms of the mean entropy. These parameterisations involve turbulent transport coefficients that need to be calculated for highly turbulent flows of stellar interiors. Analytical approaches, such as first-order smoothing, involve approximations that are ill-suited for stellar conditions and may yield inaccurate results. A numerical method for determining the turbulent transport coefficients relevant for the electromotive force is the test-field method (Schrinner et al. 2005, 2007), but for angular momentum or heat transport no similar methods have been developed yet. This means that the turbulent transport coefficients used in mean-field models are often based on educated guesses or they are even used as free parameters. Despite these shortcomings, hydrodynamical mean-field models are capable of producing SL differential rotation (Brandenburg et al. 1992; Kitchatinov & Rüdiger 1995; Rempel 2005; Kitchatinov & Olemskoy 2011) as well as basic properties of the rotation in some other stars (Küker & Rüdiger 2011; Kitchatinov & Olemskoy 2011; Hotta & Yokoyama 2011). However, obtaining AS differential rotation is less straightforward for mean-field models (e.g., Kitchatinov & Rüdiger 2004). At the same time, mean-field dynamo models also reproduce some features of solar and stellar magnetic cycles either by including turbulent inductive effects (e.g., Käpylä et al. 2006; Pipin & Kosovichev 2011), or by applying the Babcock-Leighton process in the so-called flux transport dynamo models (e.g., Choudhuri et al. 1995; Dikpati & Charbonneau 1999; Karak 2010; Karak et al. 2014; Miesch & Dikpati 2014).

On the other hand, there have been some successes in modelling the differential rotation and magnetic fields using global convection simulations, mainly in recent years (Miesch et al. 2006; Ghizaru et al. 2010; Racine et al. 2011; Käpylä et al. 2012, 2013; Warnecke et al. 2013; Augustson et al. 2014). However, owing to the extreme parameter regimes of the Sun, realistic simulations are not possible at present. Nevertheless, the simulations are able to reproduce solar values of the Coriolis number Co, which measures the relative importance of rotation and turbulent convection. For certain values of Co, but with different values for other parameters, such as the fluid and magnetic Reynolds and Prandtl numbers, simulations occasionally produce AS differential rotation (e.g., Matt et al. 2011; Käpylä et al. 2014), poleward migration of the large-scale magnetic fields (Gilman 1983; Käpylä et al. 2010b; Nelson et al. 2013), no clear magnetic cycles (Brown et al. 2010), or sometimes even no appreciable large-scale contribution to the magnetic field (Brun et al. 2004).

According to mean-field hydrodynamics, differential rotation is generated from the anisotropy of the Reynolds stress, which is parameterised in terms of the so-called Λ effect (Rüdiger 1980, 1989). A radially increasing (SL) angular velocity results if horizontal turbulent velocities dominate over vertical velocities, while AS rotation follows if radial motions (even laminar ones) are dominant. The importance of the Λ effect depends on the rotational influence on the turbulence, i.e., the value of Co, which is the ratio of the convective turnover time to the rotation period. At large Co, the SL rotation is more favourable and the transition from SL to AS rotation depends on the Coriolis number (Brun & Palacios 2009; Chan 2010; Käpylä et al. 2011a,b; Guerrero et al. 2013; Gastine et al. 2013, 2014). Using Boussinesq convection, Gastine et al. (2014) discovered that near the transition from AS to SL rotation, both states are possible, depending on the initial conditions of the simulations. This has been independently verified by Käpylä et al. (2014) in fully compressible hydrodynamic convection simulations. If this discovery were to apply to the Sun, this might have important consequences because young rapidly rotating stars, which preferably possess SL rotation, slowly spin down due to loss of the angular momentum and can persist in the SL rotation state even when their rotation is slow. Some doubts have already been expressed by Fan & Fang (2014), who found that the bi-stability disappears when magnetic fields are present.

In the magnetohydrodynamic (MHD) case, the situation is more complicated than in the case of pure hydrodynamics. A dynamically significant magnetic field, which possibly varies cyclically, introduces extra time dependent effects into the system that are capable of influencing the fluid flow both through large-scale effect (Malkus-Proctor effect; see Malkus & Proctor 1975; Brandenburg et al. 1992) and small-scale effects (Reynolds and Maxwell stresses and therefore the Λ effect; see Kitchatinov et al. 1994). Therefore, the magnetic field tries to destabilize the equilibrium states of the rotation. To explore to what extent the presence of a dynamically significant large-scale magnetic field affects the bistable nature of the differential rotation, we perform several simulations with the same setup as in Käpylä et al. (2014), but our simulations include magnetic fields. Similar to their work, we perform two types of simulations. In one of them we run the simulations from scratch, i.e., with an initially rigid rotation profile. Then we take either a SL state or an AS state and vary the rotational influence by varying the radiative heat conductivity to identify the transition. We analyse the activity cycles using diagnostic tools of stellar activity in Sect. 3.7. Next we measure the temporal variations of the Lorentz forces (both from large-scale and small-scale contributions) to understand the temporal variations of the large-scale flows observed in the simulations (Sect. 3.8). Finally, we compute the contributions of Reynolds stress, Maxwell stress, and the stresses from the azimuthally averaged mean flow and mean magnetic field to the angular momentum balance (Sect. 3.9). Then we study the influence of the magnetic field on the angular momentum transport by comparing the results with the hydrodynamic simulations.

2. The model

2.1. Basic equations

Our model is similar to many earlier studies (Käpylä et al. 2012, 2013; Cole et al. 2014). The hydrodynamic part of this model has been used in Käpylä et al. (2014). We model a spherical wedge with radial, latitudinal, and longitudinal extents r₀ ≤ r ≤ r₁, θ₀ ≤ θ ≤ π − θ₀, and 0 ≤ φ ≤ φ₀, respectively. Here, r₀ = 0.72 R_⊙ and r₁ = 0.97 R_⊙ are the positions of the bottom and top of the computational domain, R_⊙ is the radius of the Sun, θ₀ = π/ 12 is the colatitude of the polar cap, and φ₀ = π/ 2 is the longitudinal extent. The following hydromagnetic equations are solved: $\begin{array}{ccc}& & \frac{\mathit{\partial }{A}}{\mathit{\partial t}}\mathrm{=}{u}\mathrm{\times }{B}\mathrm{-}{\mathit{\mu }}_{\mathrm{0}}\mathit{\eta }{J}\mathit{,}\\ & & \frac{\mathit{D}\ln \mathit{\rho }}{\mathit{Dt}}\mathrm{=}\mathrm{-}{\nabla }\mathrm{·}{u}\mathit{,}\\ & & \frac{\mathit{D}{u}}{\mathit{Dt}}\mathrm{=}{g}\mathrm{-}\mathrm{2}{{\Omega }}_{\mathrm{0}}\mathrm{\times }{u}\mathrm{+}\frac{\mathrm{1}}{\mathit{\rho }}\left({J}\mathrm{\times }{B}\mathrm{+}{\nabla }\mathrm{·}\mathrm{2}\mathit{\nu \rho }{S}\mathrm{-}{\nabla }\mathit{p}\right)\mathit{,}\\ & & \end{array}$where A is the magnetic vector potential, B = ∇ × A is the magnetic field, J = ∇ × B/μ₀ is the current density with μ₀ being the vacuum permeability, u is the velocity, D/Dt = ∂/∂t + u·∇ is the advective derivative, ρ is the density, s is the specific entropy, T is the temperature, p is the pressure, ν is the constant kinematic viscosity, g = −GM_⊙r/r³ is the gravitational acceleration with M_⊙ being the mass of the Sun, Ω₀ = (cosθ, − sinθ,0)Ω₀ is the angular velocity vector, ${{S}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{u}}_{\mathit{i}\mathrm{;}\mathit{j}}\mathrm{+}{\mathit{u}}_{\mathit{j}\mathrm{;}\mathit{i}}\mathrm{)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}{\nabla }\mathrm{·}{u}$ is the rate of strain tensor, where the semicolons denote covariant differentiation. The radiative and subgrid scale (SGS) heat fluxes are given by ${{F}}^{\mathrm{rad}}\mathrm{=}\mathrm{-}\mathit{K}{\nabla }\mathit{T} \mathrm{and} {{F}}^{\mathrm{SGS}}\mathrm{=}\mathrm{-}{\mathit{\chi }}_{\mathrm{SGS}}\mathit{\rho T}{\nabla }\mathit{s,}$(5)respectively. Here K is the radiative heat conductivity and χ_SGS is the turbulent heat diffusivity, which represents the unresolved convective transport of heat (Käpylä et al. 2013). The fluid obeys the ideal gas law p = (γ − 1)ρe, where γ = c_P/c_V = 5/3 is the ratio of specific heats at constant pressure and volume, and e = c_VT is the specific internal energy.

2.2. Initial and boundary conditions

The initial hydrostatic state is isentropic, so the temperature is given by $\frac{\mathit{\partial T}}{\mathit{\partial r}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}{\mathit{r}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{V}}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}{\mathit{n}}_{\mathrm{ad}}\mathrm{+}\mathrm{1}\mathrm{)}}\mathit{,}$(6)where n_ad = 1.5 is the polytropic index and the value of ∂T/∂r at r = r₀ is fixed. The density stratification follows from hydrostatic equilibrium. The initial state chosen is not in thermodynamic equilibrium but closer to the final convecting state to reduce the needed computational time to reach a thermally relaxed state. The heat conductivity profile is chosen such that radiative diffusion is responsible for supplying the energy flux into the system. Radiative diffusion becomes progressively less efficient towards the surface (Käpylä et al. 2011a). As in Käpylä et al. (2013, 2014), this is achieved by taking a depth-dependent polytropic index n(r) = δn(r/r₀)^-15 + n_ad − δn for the radiative conductivity K(r) = K₀ [n(r) + 1], where the reference conductivity is ${\mathit{K}}_{\mathrm{0}}\mathrm{=}\mathrm{(}\mathrm{ℒ}\mathit{/}\mathrm{4}\mathit{\pi }\mathrm{)}{\mathit{c}}_{\mathrm{V}}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}{\mathit{n}}_{\mathrm{ad}}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\rho }}_{\mathrm{0}}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}{\mathit{R}}_{\mathrm{\odot }}}$, with ℒ being the non-dimensional luminosity. Note that n = n_ad at the bottom and n →n_ad − δn towards the surface. Hence, K decreases towards the surface like r^-15 such that the value of δn regulates the flux that is carried by convection (Brandenburg et al. 2005; Käpylä et al. 2014).

Along with the imposed energy flux at the bottom boundary F_b = −(K∂T/∂r)_{r = r0}, the values of Ω₀, ν, η, and ${\overline{\mathit{\chi }}}_{\mathrm{SGS}}\mathrm{=}{\mathit{\chi }}_{\mathrm{SGS}}$ at the middle of the convection zone r = r_m = 0.845 R_⊙ are defined. The turbulent heat conductivity χ_SGS is piecewise constant above r> 0.75 R_⊙ with ${\mathit{\chi }}_{\mathrm{SGS}}\mathrm{=}{\overline{\mathit{\chi }}}_{\mathrm{SGS}}$ in 0.75 R_⊙<r< 0.95 R_⊙, and ${\mathit{\chi }}_{\mathrm{SGS}}\mathrm{=}\mathrm{1.35}{\overline{\mathit{\chi }}}_{\mathrm{SGS}}$ at r ≥ 0.95 R_⊙. At r< 0.75 R_⊙, χ_SGS tends smoothly to zero (see Fig. 1 of Käpylä et al. 2011a). We fix χ_SGS in such a way that at r = r₁ it corresponds to 5 × 10⁸ m² s^-1 in physical units. We also assume the density and temperature at r = r₀ to have solar values, ρ₀ = 200 kg m^-3 and T₀ = 2.23 × 10⁶ K.

Radial and latitudinal boundaries are assumed to be impenetrable and stress free for the flow, whereas for the magnetic field we assume a radial field condition on the outer radial boundary and perfect conductor conditions on the lower radial and latitudinal boundaries; see Käpylä et al. (2013) for details. Density and entropy have vanishing first derivatives on the latitudinal boundaries. A black body condition with σT⁴ = −K∇_rT − χ_SGSρT∇_rs, where σ is related to the Stefan-Boltzmann constant, is applied on the upper radial boundary. However the value of σ is modified to attain the desired values of surface temperature and energy flux. Moreover, we choose σ in such a way that in the initial non-convecting state the flux at the surface carries the total luminosity through the boundary. We use small-scale, low-amplitude Gaussian noise as initial condition for the velocity and magnetic fields.

As discussed by Käpylä et al. (2013, 2014), in our fully compressible simulation the time step is severely limited if we were to use the solar luminosity. This would imply both huge Rayleigh and small Mach numbers. This problem is avoided by taking roughly 10⁶ times higher luminosity in the simulation than in the Sun. However, the convective velocity u becomes then in our simulation 100 times larger than in the Sun because the convective energy flux F_conv scales as ρu³. Therefore to achieve the same rotational influence on the flow as in the Sun, we need to increase Ω by the same factor. Consequently, we have the relations: Ω_sim = (L₀/L_⊙)^1/3Ω_⊙   and   u_sim = (L₀/L_⊙)^1/3u_⊙, where L₀ is the luminosity in simulation, L_⊙ ≈ 3.84 × 10²⁶ W is the solar luminosity, and Ω_⊙ is the average solar rotation rate ≈ 2.7 × 10^-6 s^-1. This allows us to quote values for angular velocity, meridional circulation, and magnetic field in physical units that can be compared with solar values and with results of other groups. However, we often quote ratios between different quantities that are obviously non-dimensional and therefore not affected. All computations are performed with the Pencil code¹.

2.3. Dimensionless parameters and diagnostics

First, we define the non-dimensional input parameters. The luminosity parameter ℒ and the normalized pressure scale height at the surface ξ are given by $\mathrm{ℒ}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{(}\mathit{G}{\mathit{M}}_{\mathrm{\odot }}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{R}}_{\mathrm{\odot }}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{,} \mathit{\xi }\mathrm{=}\frac{\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{c}}_{\mathrm{V}}{\mathit{T}}_{\mathrm{1}}}{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}{\mathit{R}}_{\mathrm{\odot }}}\mathit{,}$(7)where T₁ is the temperature at the surface. The influence of rotation is measured by the Taylor number, $\mathrm{Ta}\mathrm{=}\mathrm{(}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{0}}\mathrm{\Delta }{\mathit{r}}^{\mathrm{2}}\mathit{/}\mathit{\nu }{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(8)where Δr = r₁ − r₀ is the thickness of the convecting shell. The fluid, magnetic, and SGS Prandtl numbers are defined as $\mathrm{Pr}\mathrm{=}\frac{\mathit{\nu }}{{\mathit{\chi }}_{\mathrm{m}}}\mathit{,} \mathrm{Pm}\mathrm{=}\frac{\mathit{\nu }}{\mathit{\eta }}\mathit{,} {\mathrm{Pr}}_{\mathrm{SGS}}\mathrm{=}\frac{\mathit{\nu }}{{\overline{\mathit{\chi }}}_{\mathrm{SGS}}}\mathit{,}$(9)respectively, where χ_m = K(r_m) /c_Pρ_m is the thermal diffusivity and ρ_m is the density, both evaluated r = r_m. Furthermore, we define the non-dimensional viscosity, $\begin{array}{c}\mathrm{˜}\\ \mathit{\nu }\end{array}\mathrm{=}\frac{\mathit{\nu }}{\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}{\mathit{R}}_{\mathrm{\odot }}}}\mathit{,}$(10)and the Rayleigh and convective Rossby numbers (Gilman 1977), $\begin{array}{ccc}\mathrm{Ra}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}\mathrm{(}\mathrm{\Delta }\mathit{r}{\mathrm{)}}^{\mathrm{4}}}{\mathit{\nu }{\overline{\mathit{\chi }}}_{\mathrm{SGS}}{\mathit{R}}_{\mathrm{\odot }}^{\mathrm{2}}}(\mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}_{\mathrm{P}}}\frac{\mathrm{d}\mathit{s}}{\mathrm{d}\mathit{r}}{)}_{{\mathit{r}}_{\mathrm{m}}}\mathit{,}{\mathrm{Ro}}_{\mathrm{c}}\mathrm{=}{\left(\frac{\mathrm{Ra}}{{\mathrm{Pr}}_{\mathrm{SGS}}\mathrm{Ta}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(11)where the entropy gradient of the non-convecting hydrostatic solution is evaluated in the middle of the convection zone, r = r_m. We also quote the initial density contrast ${\mathrm{\Gamma }}_{\mathit{\rho }}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\equiv }\mathit{\rho }\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}\mathit{/}\mathit{\rho }\mathrm{\left(}{\mathit{r}}_{\mathrm{1}}\mathrm{\right)}$.

As diagnostic quantities we define the fluid and magnetic Reynolds numbers, and the Coriolis number as $\mathrm{Re}\mathrm{=}\frac{{\mathit{u}}_{\mathrm{rms}}}{\mathit{\nu }{\mathit{k}}_{\mathrm{f}}}\mathit{,} \mathrm{Rm}\mathrm{=}\frac{{\mathit{u}}_{\mathrm{rms}}}{\mathit{\eta }{\mathit{k}}_{\mathrm{f}}}\mathrm{=}\mathrm{Pm}\mathrm{Re}\mathit{,} \mathrm{Co}\mathrm{=}\frac{\mathrm{2}{\mathrm{\Omega }}_{\mathrm{0}}}{{\mathit{u}}_{\mathrm{rms}}{\mathit{k}}_{\mathrm{f}}}\mathit{,}$(12)where ${\mathit{u}}_{\mathrm{rms}}\mathrm{=}\sqrt{\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{⟨}{\mathit{u}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{u}}_{\mathit{\theta }}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{r\theta \phi t}}}$ is the volume and time averaged rms velocity during the time when the simulation is thermally relaxed. We exclude u_φ from u_rms as it is dominated by the differential rotation, and use k_f = 2π/ Δr as an estimate of the wavenumber of the largest eddies. The Taylor number can also be written as Ta = Co²Re²(k_fR_⊙)⁴.

We define mean values as averages over longitude and time and denote these by an overbar. Sometimes we also perform additional averaging over latitude and/or radius, which we always mention explicitly.

3. Results

First, we perform a set of simulations for different values of the radiative conductivity starting from the initial condition described in Sect. 2.2. These are referred to as Runs A–E in Table 1. Except for the allowance of a magnetic field, all other input parameters of our Runs A–E are identical to those of Käpylä et al. (2014). However, we have an additional Run BC in this set, whose radiative conductivity lies between those of Runs B and C. It turns out that Runs A–BC produce AS differential rotation, while Runs C–E produce SL differential rotation. Next, we perform two further sets of simulations where we use Runs A and D with AS and SL differential rotation, respectively, as progenitors to study the possibility of bi-stability of the rotation profile.

3.1. Energy fluxes in our dynamo runs

The radiative conductivity in our model is controlled by the parameter δn, regulating the fractional flux that convection has to transport. Increasing δn reduces the radiative flux and increases the convective flux and thus u_rms. Having Ω₀ fixed, changing δn affects the rotational influence on the convection via the convective velocities (for further details see Käpylä et al. 2014). Hence, different values of δn in Runs A–E imply different values of Co. Thus, Run E is more rotationally dominated than Run A. For Run A with δn = 2.5, the convective flux dominates over the other fluxes and the radiative flux transports a very small fraction of the luminosity. For the definitions of the fluxes, we refer to Eqs. (26)–(31) of Käpylä et al. (2013).

Fig. 1 Contributions of different energy fluxes of Runs A (top) and E (bottom). The black (red) lines correspond to magnetohydrodynamic (hydrodynamic) case. Thin solid: radiative, dashed: convective, long dashed: viscous, dash-dotted: kinetic energy, dash-triple-dotted: SGS, and thick solid: total. Dotted horizontal lines indicate the zero and unity values. The vertical dotted line indicates the position of the middle of the convection zone r = rm.

Fig. 2 Variation of ${\mathrm{(}{\mathit{F}}_{\mathrm{conv}}\mathit{/}\mathit{\rho }{\mathit{u}}_{\mathrm{rms}}^{\mathrm{3}}\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$ averaged over the whole convection zone as a function of urms from different runs. Black asterisks: Runs A−E, red diamonds: Runs E1–E4, and blue triangles: Runs A1−A8. The dotted line shows the linear dependence between urms and ${\mathrm{(}{\mathit{F}}_{\mathrm{conv}}\mathit{/}\mathit{\rho }{\mathit{u}}_{\mathrm{rms}}^{\mathrm{3}}\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$.

In the statistically stationary state, the total luminosity L_tot(r) = 4πr²ℱ_tot(r) is constant, where ℱ_tot is the time averaged total energy flux. In Fig. 1 we show the radial dependence of the contributions from radiation (L_rad) and convection (L_conv), as well as kinetic (L_kin), viscous (L_visc), and subgrid scale (L_SGS) energy fluxes in the convection zone for Runs A and E. For comparison we also show the fluxes from the corresponding hydrodynamic simulations of Käpylä et al. (2014) with red lines. We see that for Run A, the convective, kinetic, and SGS energy fluxes have decreased in the lower part of the convection zone in the magnetic case. In Run E there is very little change in the SGS flux and only a small reduction of the convective and kinetic energy fluxes is visible. In Table 1 we show the fractions of the radiative and the convective fluxes at the middle of the convection zone for all the runs. In Runs A, B, and BC, the convective flux is above 75% and the radiative flux is less than 25%. By contrast, Runs C, D, and E have a convective flux of less than 70% and the radiative flux is larger than 25%.

We find that the rms-velocity is compatible with a one-third power proportionality to the convective energy flux. This is the case at least in the narrow range of parameters studied here; see Fig. 2. This is in agreement with the scaling used in connection with the artificially high luminosities used in our simulations (see also Brandenburg et al. 2005).

Fig. 3 Distribution of angular velocity $\overline{\mathrm{\Omega }}$ in the meridional plane from Runs A–E. The $\overline{\mathrm{\Omega }}$ is computed from Ω first by the longitudinal average and then the time average over the last few cycles. The arrows in the leftmost panel show the colatitudes at which the latitudinal differential rotation is computed in Eq. (13).

3.2. Differential rotation

In Fig. 3 we show the rotation profile $\overline{\mathrm{\Omega }}\mathrm{=}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathit{/}\mathit{r}\sin \mathit{\theta }\mathrm{+}{\mathrm{\Omega }}_{\mathrm{0}}$ for Runs A–E. We see that in Runs A, B, and BC the equator rotates slower than the mid- and high latitudes, which is opposite to the rotation profile observed in the Sun. However, we find that the regions near the latitudinal boundaries have slower rotation than mid-latitudes, which was not observed in the hydrodynamical simulations of Käpylä et al. (2014). We have repeated Run B by increasing and decreasing the latitudinal extent. Runs B′ and B′′ in Table 1 correspond to these two cases where the latitudinal end points are at ± 84° and ± 66°, respectively. In Fig. 4, we show $\overline{\mathrm{\Omega }}$ profiles for these two runs, whereas in Fig. 5 we show their latitudinal variations at r = 0.96 R_⊙. From these two plots we see that the rotation profiles are very similar in all these runs up to about ± 50° latitudes and the significant departures appear only near the boundaries. However, the slowly rotating high-latitude branch still exists in all runs. Therefore this is probably not due to our restricted latitudinal extent, but may be a real feature that might be caused by the magnetic fields.

By comparing the rotation profiles from hydrodynamic simulations of Käpylä et al. (2014) for Runs A–E, we see that the pole-equator differential rotation is generally weaker in the magnetic runs, the reduction being strongest in the runs with a slowly rotating equator. This will be discussed in Sect. 3.5, where we give the ratio of the rotational energy of the hydrodynamic and magnetic simulations. The overall reduction of the pole–equator differential rotation by magnetic fields agrees with what has been reported before (e.g., Brun et al. 2004; Beaudoin et al. 2013; Fan & Fang 2014). In particular, for slow rotation Fan & Fang (2014) found a switch from AS to SL differential rotation, which agrees with our results reported below. By contrast, Brun et al. (2004) and Beaudoin et al. (2013) only consider runs with SL rotation, but in the models by Brun et al. (2004) the angular velocity at high latitudes is not reduced much by magnetic fields. In the models of Beaudoin et al. (2013), on the other hand, there is a lower overshoot layer, where angular velocity is constant and equal to that at high latitudes when there is a magnetic field. This is, however, not the case in their non-magnetic runs, so the magnetic field leads to a strong reduction in their case, which is similar to ours, but different from what Brun et al. (2004) found. We caution that comparing runs with and without overshoot layers can be misleading because of different physical effects involved.

To illustrate the reduction of differential rotation in our magnetic runs, we show in Fig. 6 the differences in $\overline{\mathrm{\Omega }}$ between the hydrodynamic and magnetic versions of Runs A and E. We see that for Run A the difference is even larger than the average rotation rate of the Sun. Possible reasons for the reduction of differential rotation will be discussed in Sect. 3.9.

In Runs C–E, the equator rotates faster than the polar regions, as is also the case in the Sun. However, unlike some of the non-magnetic solar-like cases of Käpylä et al. (2014), we never obtain polar vortices or jet-like structures with magnetic fields. It is also interesting to note that the pole-equator difference in the rotational velocity is comparable to that of the Sun. The solar value (~430 nHz) is marked in each colourbar of Fig. 3.

Fig. 4 Same as Run B in Fig. 3, except here the latitudinal extent is different. For the left panel, the latitudinal boundaries are at ± 84°, whereas for the right panel they are at ± 66°.

Fig. 5 Variations of angular velocities $\overline{\mathrm{\Omega }}$ at r = 0.96 R⊙ from Runs B (solid line), B′ (red dashed), and B′′ (blue dash-dotted).

Fig. 6 Difference of $\overline{\mathrm{\Omega }}$ between the HD and MHD runs. The left and right panels are for Runs A and E, respectively.

3.3. Identifying the SL to AS transition

We measure the relative radial and latitudinal differential rotation by the quantities $\begin{array}{ccc}{\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{eq}}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{bot}}}{{\mathrm{\Omega }}_{\mathrm{eq}}}\mathit{,} {\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{eq}}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{55}}}{{\mathrm{\Omega }}_{\mathrm{eq}}}\mathit{,}& & \end{array}$(13)where ${\mathrm{\Omega }}_{\mathrm{eq}}\mathrm{=}\overline{\mathrm{\Omega }}\mathrm{(}{\mathit{r}}_{\mathrm{1}}\mathit{,\pi }\mathit{/}\mathrm{2}\mathrm{)}$ and ${\mathrm{\Omega }}_{\mathrm{bot}}\mathrm{=}\overline{\mathrm{\Omega }}\mathrm{(}{\mathit{r}}_{\mathrm{0}}\mathit{,\pi }\mathit{/}\mathrm{2}\mathrm{)}$ are the equatorial rotation rates at the surface and at the base of the convection zone, and ${\mathrm{\Omega }}_{\mathrm{55}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\left[}\overline{\mathrm{\Omega }}\mathrm{\right(}{\mathit{r}}_{\mathrm{1}}\mathit{,}{\mathrm{35}}^{\mathrm{◦}}\mathrm{)}\mathrm{+}\overline{\mathrm{\Omega }}\mathrm{(}{\mathit{r}}_{\mathrm{1}}\mathit{,}{\mathrm{145}}^{\mathrm{◦}}\mathrm{\left)}\mathrm{\right]}$ is the rotation rate at latitudes ± 55° computed as an average of $\overline{\mathrm{\Omega }}$ at 35° and 145° colatitudes on the outer radius. The arrows in the leftmost panel of Fig. 3 show the positions of these points in the r − θ plane. The values of ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}$ and ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$, listed in Table 1, help us to identify AS and SL differential rotation. The SL differential rotation implies ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{>}\mathrm{0}$ and ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\mathit{>}\mathrm{0}$. Following this definition, we see that Runs A, B, and BC are classified as AS and Runs C–E as SL. Hence we see that there is a transition from AS to SL rotation around Co ~ 1.4. Although this transition has been reported extensively in the literature (e.g., Gilman 1977; Rieutord et al. 1994; Käpylä et al. 2011b; Gastine et al. 2013; Guerrero et al. 2013), the possibility of a bi-stability has only recently been discovered (Gastine et al. 2014; Käpylä et al. 2014), i.e., AS and SL rotation profiles can be obtained for the same input parameters. We compare the current results with the hydrodynamical results and we see that the transition happens at a slightly larger value of Ro_c, also manifested by the change of the AS hydrodynamical counterpart of Run C to SL in the MHD regime. The transition observed in the current dynamo cases is less abrupt than that of earlier hydrodynamic studies. Therefore we conclude that the magnetic field helps to produce SL differential rotation, which has also been found in the recent anelastic simulations of Fan & Fang (2014).

Fig. 7 Radial (left panel) and latitudinal (right panel) differential rotation, defined by Eq. (13), for Runs A–E (green crosses), A0–A8 (red dashed line with triangles), and D0–D4 (black asterisks with dotted line). The other points are taken from the hydrodynamical simulations (black diamonds: Runs A–E, blue dotted line with asterisks: Runs D0–D4, and red dashed line with squares: Runs B0–B10) of Käpylä et al. (2014). The horizontal black dotted lines show the zero value. A red (black) arrow shows the transition point when differential rotation in the MHD (HD) simulations changes from AS to SL rotation with the decrease of Roc.

In Fig. 7 we show differential rotation parameters ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}$ and ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$ computed from Eq. (13) for all the runs as functions of Ro_c and Ra. To compare with the hydrodynamic simulations of Käpylä et al. (2014), we have analysed the latitudinal differential rotation in their data with our new definition (13). The hydrodynamic values shown in Fig. 7 are now considerably smaller for the AS branch than the ones reported by Käpylä et al. (2014). Note that, had we defined ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$ as the difference of rotation between equator and the endpoints of the domain at θ = θ₀ and π − θ₀ as in Käpylä et al. (2013, 2014), instead of Ω₅₅, as we do here, we would have obtained smaller values of ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$ because of the slowly rotating high-latitude regions in our magnetic Runs A–C. Another way of characterizing the SL or AS differential rotation can be done following Käpylä et al. (2011b), who approximated the surface rotation profile in terms of Gegenbauer polynomials, which is, $\mathrm{\Omega }\mathrm{=}{\mathrm{\Omega }}_{\mathrm{0}}{\sum }_{\mathit{\ell }\mathrm{=}\mathrm{1}\mathit{,}\mathrm{3}\mathit{,}\mathrm{5}}{\mathit{\omega }}_{\mathit{\ell }}{\mathit{P}}_{\mathit{\ell }}^{\mathrm{1}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}\mathit{/}\sin \mathit{\theta }$. The sign of w₃ indicates whether a rotation is SL or AS. Following this procedure, we obtain the same conclusion for the classification of the SL and AS differential rotation.

3.4. Checking for flow bi-stability

Next we study the flow bi-stability and take AS and SL cases as initial conditions. Firstly, we have performed a set of simulations by starting from the saturated state of Run A with AS differential rotation and decreasing δn slowly, which corresponds to Runs A1–A8 in Table 1. The differential rotation parameters of these runs are shown as red triangles in Fig. 7. Secondly, we start from Run D with SL differential rotation and increase δn slowly to produce Runs D1–D4. These are shown as black asterisks in Fig. 7. We see that both sets of simulations produce similar results and there is no evidence for the existence of multiple solutions at the same parameters. Therefore, we conclude that the bistable nature of the differential rotation, recently discovered by Gastine et al. (2014) and Käpylä et al. (2014), disappears when dynamically important magnetic fields are allowed to be generated. This conclusion is supported by the recent study of Fan & Fang (2014), who find a stable SL differential rotation independent of the history of their convective dynamo simulations.

3.5. Meridional circulation

Fig. 8 Meridional circulation from Runs A and E. The arrows show the direction of flow ${\overline{{u}}}_{\mathrm{m}}\mathrm{\equiv }\mathrm{\left(}\mathit{u̅}{}_{\mathit{r}}\mathit{,}\mathit{u̅}{}_{\mathit{\theta }}\mathrm{\right)}$ and the background colour shows uθ. Upper (lower) panels are from magnetohydrodynamic (hydrodynamic) simulations.

For all the AS cases (Runs A, B, and BC), we find single cell meridional circulation with poleward flow near the surface and equatorward flow near the bottom of the convection zone. This is also the usual assumption in flux transport dynamo models (e.g., Dikpati & Charbonneau 1999), although in such models only the equatorward motion at the bottom of the convection zone matters (Hazra et al. 2014). However, as we go to the SL differential rotation cases, i.e., from Run C to Runs D and E, the meridional circulation becomes weaker and shows multiple cells in radius and latitude, which has been detected in recent observations (Zhao et al. 2013; Schad et al. 2013; Kholikov et al. 2014). Guerrero et al. (2013) also find multi-cell meridional circulation for SL differential rotation and single or two-cell circulation for AS differential rotation in their hydrodynamic simulations. Again in our SL cases, we observe some equatorward flow near the bottom of the convection zone, at least in mid to high latitudes.

Table 1 shows that E_mer/E_kin decreases rapidly from Runs A to E (with increasing Co the energy in the azimuthal component increases). The upper two panels of Fig. 8 show the meridional circulation for an AS (Run A) and a SL case (Run E). Note that, irrespective of the differential rotation profile, we obtain a poleward flow near the surface and its amplitude is in agreement with solar surface observations (see e.g., Hathaway & Rightmire 2010; Zhao et al. 2013). The fact that it is poleward both for AS and SL rotation suggests that the meridional circulation is not just a consequence of differential rotation.This has been discussed in detail by Rüdiger (1989), who has shown that baroclinic forcing can also be an important driver for the origin of meridional circulation; see Miesch & Toomre (2009) (Sect. 3), who have demonstrated how in simulations the convective (and magnetic) angular momentum flux maintain meridional circulation in the solar convection zone through the gyroscopic pumping. We note that in our simulations we do not have the near-surface shear layer, which helps to produce a poleward flow in the upper layers through inward angular momentum transport possibly by the down-flow plumes (Kitchatinov & Rüdiger 2005; Miesch & Hindman 2011; Hotta et al. 2015).

It is important to compare these results with the hydrodynamical counterparts of the same models shown in the two lower panels. We see that the hydrodynamic flow is much stronger, although the overall pattern is not very different. In Table 2, we compare the energy ratios ${\mathit{E}}_{\mathrm{mer}}^{\mathrm{HD}}\mathit{/}{\mathit{E}}_{\mathrm{mer}}^{\mathrm{MHD}}$ and ${\mathit{E}}_{\mathrm{rot}}^{\mathrm{HD}}\mathit{/}{\mathit{E}}_{\mathrm{rot}}^{\mathrm{MHD}}$ of meridional circulation and rotation respectively with their hydrodynamic counterparts. The magnetic field clearly suppresses the circulation in the AS cases, Runs A–BC, but also in Run C, whereas in the other SL cases the effect is small. Moreover, the flow shows significant temporal variation, which will be explored later in Sect. 3.8.

3.6. Magnetic variability and butterfly diagrams

The large-scale, spatio-temporal organization of the magnetic field can be seen from a time-latitude or butterfly diagram of ${\overline{\mathit{B}}}_{\mathit{\phi }}$, for example. The degree of radial coherence can be judged by looking at different depths. We show such butterfly diagrams for Runs A–E both at r = 0.74 R_⊙ (Fig. 9) and at r = 0.96 R_⊙ (Fig. 10), where ${\overline{\mathit{B}}}_{\mathit{\phi }}$ is given in Gauss. Here we only show the first 70 years of each simulation, although in some cases we ran for longer times (see below). We see that Runs A, B, and BC, which produce AS differential rotation, show prominent activity cycles, but no clear polarity reversals. Therefore, these activity cycles are different from the Hale polarity cycle of the Sun.However, we want to stress here that polarity reversals are verified only for the Sun, whereas for stellar cycles, most commonly detected either from photometry or from Ca H & K lines with spectroscopy, such information is not retrievable, and therefore the reported variability might equally well be related to variations in the magnetic field strength as to polarity reversals.

Fig. 9 Butterfly diagrams: contours of the toroidal field ${\overline{\mathit{B}}}_{\mathit{\phi }}$ (in Gauss) at r = 0.74 R⊙ from Runs A, B, BC, C, D, and E (top to bottom).

Fig. 10 Same as Fig. 9 but here ${\overline{\mathit{B}}}_{\mathit{\phi }}$ is shown at r = 0.96 R⊙.

For Runs A and B, the activity cycles are associated with a weak poleward propagation at high latitudes. Another aspect we notice is that, as we go from Run A to Run BC, the activity cycles appear at a later time. This is surprising given that the rotational influence on the turbulence actually increases. For Run C, with SL rotation, the variations are rather irregular and show periods of very weak mean fields. The variations in Runs D and E are irregular and a different dynamo mode appears to be excited in these simulations, which is also characterized by a weaker influence on the differential rotation than in the AS cases. By comparing Figs. 9 and 10, we also see that ${\overline{\mathit{B}}}_{\mathit{\phi }}$ is stronger near the bottom of the convection zone, which could be caused by downward pumping of the mean magnetic field (Nordlund et al. 1992; Brandenburg et al. 1996; Tobias et al. 2001).

In comparison to earlier work, we note that in several recent global dynamo simulations (Ghizaru et al. 2010; Racine et al. 2011; Brown et al. 2011; Käpylä et al. 2012, 2013; Cole et al. 2014), regular cycles develop, resulting in butterfly diagrams with polarity reversals and sometimes even equatorward migration of toroidal field at low latitudes. However, all these simulations are for larger Coriolis numbers (for example, Co ≥ 5 in Cole et al. 2014), in which regime more coherent fields have been found to be favoured (Brown et al. 2010, 2011).

3.7. Diagnostic stellar activity diagrams

The nature of stellar cycles can be characterized by the ratio of cycle frequency ω_cyc = 2π/P to the rotation rate Ω₀, where P is an estimate of the cycle period. There is a tendency for stars of different characteristics to group at different positions in a “diagnostic” diagram of ω_cyc/ Ω₀ versus Co (Brandenburg et al. 1998; Saar & Brandenburg 1999). The more rapidly rotating stars of Käpylä et al. (2013) were found to be located in three groups with increasing slope in two of them and decreasing slope in one.

In the present work, the values of Co are much smaller, so it is important to repeat such an analysis for the more slowly rotating stars of the present paper.However, even by visual inspection of the butterfly diagrams, it is evident that these variations are not strictly harmonic, and therefore Fourier transform is not useful for the analysis. Instead, we use the phase dispersion method (Pelt 1983; Lindborg et al. 2013). It is based on the statistic ${\mathit{D}}^{\mathrm{2}}\mathrm{\left(}\mathit{P}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\frac{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}\mathrm{-}\mathrm{1}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{N}}\mathit{g}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathit{,}{\mathit{t}}_{\mathit{j}}\mathit{,P,}\mathrm{\Delta }\mathit{t}\mathrm{\left)}\mathrm{\right[}\mathit{f}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathrm{-}\mathit{f}\mathrm{\left(}{\mathit{t}}_{\mathit{j}}\mathrm{\right)}{\mathrm{]}}^{\mathrm{2}}}{\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}\mathrm{-}\mathrm{1}}\sum _{\mathit{j}\mathrm{=}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{N}}\mathit{g}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathit{,}{\mathit{t}}_{\mathit{j}}\mathit{,P,}\mathrm{\Delta }\mathit{t}\mathrm{)}}\mathit{,}$(14)where f(t_i),i = 1,...,N is the input time series, σ² is its variance, g(t_i,t_j,P,Δt) is the selection function, which in the general case is significantly different from zero only when $\begin{array}{ccc}& & {\mathit{t}}_{\mathit{j}}\mathrm{-}{\mathit{t}}_{\mathit{i}}\mathrm{\approx }\mathit{kP,k}\mathrm{=}\mathrm{\pm }\mathrm{1}\mathit{,}\mathrm{\pm }\mathrm{2}\mathit{,}\mathit{...}\mathrm{and}\\ & & \left|{\mathit{t}}_{\mathit{j}}\mathrm{-}{\mathit{t}}_{\mathit{i}}\right|\mathrm{\le }\mathrm{\Delta }\mathit{t}\mathit{.}\end{array}$In the latter condition, Δt is the so-called correlation length. In our case, the time series is evenly sampled, in which case a modification to the general top-hat selection function is needed to prevent artefacts: in this study we choose g as the product of two gaussians: one with a half width at half maximum (HWHM) of Δt and the other with an HWHM of a preselected phase separation limit, for which we adopt 0.1.

Fig. 11 Phase dispersion analysis results; from top to bottom, Runs A–E, in the middle column for the northern hemisphere, on the right for the southern hemisphere. In the left column, we show the original time series analysed (blue/red for north/south). The time series are normalized by their means, so values are not shown.

For the particular case when Δt is longer than the full data span, the D²(P) statistics is essentially a slight reformulation of the well-known Stellingwerf statistics (Stellingwerf 1978). As the correlation length is made shorter, we match nearby cycles in a progressively narrower region, and consequently estimate a certain mean period, which does not need to be coherent for the full time span. After having continued Runs C, D and E for a longer time, we apply this statistic to the time series of $\overline{\mathit{B}}\begin{array}{c}\mathrm{2}\\ \mathit{\phi }\end{array}$ at r = 0.96 R_⊙, averaged over 10° to 50° latitude, separately for north and south. This is shown in the left column of Fig. 11, where we now depict the full length of the simulation, excluding the initial exponential growth phase of the dynamo.

To determine the possible average period of the cycle in the time series we proceed as follows: first we calculate the D² statistics for preselected period and correlation length ranges. Then, by inspecting the phase dispersion diagrams (middle and rightmost columns of Fig. 11), we detect the correlation lengths for which there exist only distinct flat minima around certain periods (this must be the case for the shortest correlation length if the lower bound has been chosen correctly). After eliminating doubles or halves of the actual periods, we obtain a set of candidate periods. The period with the strongest minimum is selected as the best candidate, others being possible modulating periods.

In Fig. 11, the results for the phase dispersion analysis are depicted; the periods with error estimates for all the runs can be found from Table 2. The error estimates were calculated by building bootstrap re-samples for the differences between pairs of data points with approximately the same time separation. For Run A, a stable period of slightly less than three years is prominent especially for the southern hemisphere, which is manifested by the fact that the minimum does not split even if the correlation length is increased. For the northern hemisphere, a period with roughly the same value is obtained, but it is less coherent and splits at correlation lengths of roughly 25 years. A secondary period of roughly four years is also present in both hemispheres. The hemispheres appear less synchronized for Run B: again, two significant periods are detected for both hemispheres, the most prominent periods being eight and seven years, respectively (albeit with large error estimates). Again, the dominant periods persist even when the correlation length is increased. For Run BC, a stable period of about 11 years is found in both hemispheres. In the southern hemisphere, it is more persistent as the correlation length is increased. For Run C, only the northern hemisphere shows a prominent period around seven years, while the splitting of the minima starts already for the correlation length of ten years in the south. Run D shows prominent periods around seven and nine years for the northern and southern hemispheres, respectively. Although splitting occurs, especially in the north, we regard these periods as significant. Finally in the case of Run E, multi-periodicity was detected for the northern hemisphere with a stronger period around 4.5 years and a weaker one around 2.5 years, while the latter period was detected also for the southern hemisphere.

The plot of ω_cyc/ Ω₀ versus Co is shown in Fig. 12a, where we see that stars with AS and SL rotation are located in two different positions in this diagram. We know that the Sun lies on the upper left branch in such a diagram for real stars (Brandenburg et al. 1998), but in Fig. 12a this branch corresponds to AS rotation, which obviously disagrees with the observations. However, it is plausible that there are stars with AS rotation that have simply not yet been observed. Our simulations therefore suggest a possible prediction in that these stars might occupy a possibly separate branch further to the left of the solar branch, which already corresponds to the domain of inactive stars. We thus expect there to be either a separate branch or a part of the branch with inactive stars like the Sun with AS rotation however.

Fig. 12 Diagnostic diagrams showing a) ωcyc/ Ω0, b) the strength of the large-scale magnetic field over the whole convection zone measured by ${\overline{\mathit{B}}}_{\mathrm{rms}}\mathrm{=}\mathrm{⟨}\mathrm{⟨}{\mathit{B}}_{\mathit{r}}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}{\mathit{B}}_{\mathit{\theta }}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}{\mathit{B}}_{\mathit{\phi }}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{r\theta t}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ vs. log Co, and c) ωcyc/ Ω0 vs. ${\overline{\mathit{B}}}_{\mathrm{rms}}$. The red colour shows the values computed for the southern hemisphere.

In Fig. 12b we show the mean magnetic field. Here the errors of ${\overline{\mathit{B}}}_{\mathrm{rms}}$ are computed as the largest departure of the mean from any one-third of the full time series. We see that for the AS differential rotation runs (Runs A, B, and BC) the magnetic field decreases with rotation rate, whereas for the SL rotation cases (Runs C, D, and E) it increases.

In Fig. 12c we plot the cycle frequency ratio against the mean magnetic field. Note that ω_cyc/ Ω₀ increases with ${\overline{\mathit{B}}}_{\mathrm{rms}}$, which agrees with the observed trend found for active and inactive stars. Following the interpretation of Brandenburg et al. (1998), the cycle frequency ratio is essentially a measure of the α effect in a mean-field dynamo, so the increasing trend in the cycle frequency ratio suggests that the α effect increases with magnetic field strength, which is referred to as anti-quenching. This interpretation hinges on some ill-known assumptions, for example, the turbulent transport in this model is assumed to depend only on the largest scale of the mean magnetic field and of course the rather unconventional assumption of anti-quenching itself. Anti-quenching of both α and turbulent diffusivity has actually been detected in simulations (Chatterjee et al. 2011), and may be possible more easily in a sphere than in a Cartesian layer, but we have at present no further indication that this interpretation is applicable to our model.

3.8. Magnetic modulation of the flow

We have seen that some of our runs show clear activity cycles. Therefore we expect to see a corresponding modulation of the flow. In Fig. 13, we show for Run A the temporal variation of the mean large-scale magnetic field ($\overline{\mathit{B}}$) normalized by B_eq, the latitudinal component of the meridional circulation u_θ(r, ± 32°) at r ≈ 0.95 R_⊙ and r ≈ 0.73 R_⊙, the mean rotation rate $\overline{\mathrm{\Omega }}\mathrm{(}\mathrm{0.95}{\mathit{R}}_{\mathrm{\odot }}\mathit{,}\mathrm{\pm }{\mathrm{32}}^{\mathrm{◦}}\mathrm{)}$, $\overline{\mathrm{\Omega }}\mathrm{\left(}\mathit{r,}{\mathrm{0}}^{\mathrm{◦}}\mathrm{\right)}$ at r = 0.73 R_⊙ and r = 0.95 R_⊙, as well as the latitudinal and radial differential rotation ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}$ and ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$, defined in Eq. (13). We see that the meridional circulation varies with the magnetic field, becoming weaker during maximum and stronger during minimum, the overall temporal variation being about 50% in this case. (The linear correlation coefficient between $\overline{\mathit{B}}$ and u_θ(0.95 R_⊙, ± 32°) ≈ −0.36, − 0.38.) This kind of weak anti-correlation between the activity cycle and the meridional flow has been found in solar observations (Chou & Dai 2001; Hathaway & Rightmire 2010) and is believed to arise at least in part from the Lorentz force of the dynamo-generated magnetic fields (see, e.g., Rempel 2006, Karak & Choudhuri 2012, Passos et al. 2012).The meridional circulation at the bottom is also weakly correlated with the activity cycle (correlation coefficients between $\overline{\mathit{B}}$ and u_θ(0.73 R_⊙, ± 32°) ≈ −0.22, − 0.47). We see that $\overline{\mathrm{\Omega }}\mathrm{(}\mathrm{0.95}{\mathit{R}}_{\mathrm{\odot }}\mathit{,}\mathrm{\pm }{\mathrm{32}}^{\mathrm{◦}}\mathrm{)}$ (Fig. 13c) also shows a weak anti-correlation with the magnetic variations (having correlation coefficient ≈ − 0.25). The strong magnetic fields during maxima change $\overline{\mathrm{\Omega }}$ by a few per cent (≈ 6%). However $\overline{\mathrm{\Omega }}\mathrm{(}\mathrm{0.95}{\mathit{R}}_{\mathrm{\odot }}\mathit{,}{\mathrm{0}}^{\mathrm{◦}}\mathrm{)}$ (Fig. 13d) shows positive correlation (correlation coefficient ≈ 0.36) and the overall variation is larger (≈ 12%). Because of this variation of $\overline{\mathrm{\Omega }}$ at the equator, the values of ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}$ and ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$ (Figs. 13e, f) show a positive correlation with the magnetic field (correlation coefficients 0.36, 0.21) with the overall variation being ~ 75% and 166%, respectively.

Fig. 13 From Run A: a) the large-scale magnetic field over the whole convection zone $\overline{\mathit{B}}\mathrm{=}\mathrm{⟨}\mathrm{⟨}{\mathit{B}}_{\mathit{r}}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}{\mathit{B}}_{\mathit{\theta }}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}{\mathit{B}}_{\mathit{\phi }}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{r\theta }}^{\mathrm{1}\mathit{/}\mathrm{2}}$ normalized by Beq; b) the latitudinal component of meridional circulation uθ(r, ± 32°) (smoothed over 5 months) at r ≈ 0.95 R⊙ (black and red) and r ≈ 0.73 R⊙ (blue and green); c) azimuthally averaged angular velocity $\overline{\mathrm{\Omega }}\mathrm{(}\mathrm{0.95}{\mathit{R}}_{\mathrm{\odot }}\mathit{,}\mathrm{\pm }{\mathrm{32}}^{\mathrm{◦}}\mathrm{)}$;d) $\overline{\mathrm{\Omega }}\mathrm{\left(}\mathit{r,}{\mathrm{0}}^{\mathrm{◦}}\mathrm{\right)}$ at r = 0.73 R⊙ (red dashed) and r = 0.95 R⊙ (black).The dashed (solid) line corresponds to the southern (northern) hemisphere; e) radial shear ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}$, and f) latitudinal shear ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$, defined in Eq. (13), as functions of time.

Fig. 14 Same as Fig. 13, but from Run E, which produces SL differential rotation.

We have seen that Runs B, BC, and C produce clear activity cycles similar to Run A and in all these runs we do see a corresponding variation in the flow. However, for Runs C, D, and E, which produce SL differential rotation, the magnetic variations are not so regular. In Fig. 14, we show the temporal variations for the SL differential rotation case Run E. We see that the large-scale magnetic field does not have a regular cycle. The meridional circulation does not appear to show a close correlation with the magnetic field (with correlation coefficients being about −0.1 and −0.3 for surface and bottom meridional circulation, respectively). However significant variations (up to 45%) exist. An irregular variation in meridional circulation is found to be crucial in modelling many aspects of the solar cycle in flux transport dynamo model (Karak & Choudhuri 2011, 2013). The situation is similar in the case of differential rotation and rotational shear: the early part of the time series (t = 30...50 yr) show an anti-correlation with the magnetic field strength, but at later times this correlation is not so obvious (the overall linear correlation coefficients are −0.32, −0.22, −0.43, 0.28, −0.55, and 0.19 for $\overline{\mathrm{\Omega }}\mathrm{(}\mathrm{0.95}{\mathit{R}}_{\mathrm{\odot }}\mathit{,}\mathrm{\pm }{\mathrm{32}}^{\mathrm{◦}}\mathrm{)}$, $\overline{\mathrm{\Omega }}\mathrm{\left(}\mathrm{0.95}{\mathit{R}}_{\mathrm{\odot }}\mathit{,}{\mathrm{0}}^{\mathrm{◦}}\mathrm{\right)}$, $\overline{\mathrm{\Omega }}\mathrm{\left(}\mathrm{0.73}{\mathit{R}}_{\mathrm{\odot }}\mathit{,}{\mathrm{0}}^{\mathrm{◦}}\mathrm{\right)}$, ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}$, and ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$, respectively). The variation in differential rotation is about 4%, whereas for the radial and latitudinal shear it is about 50% and 60%, respectively.

Significant variations observed in the large-scale flows in all the simulations motivate us to measure the Lorentz force. The Lorentz force can change the flow by acting in two ways, through large-scale and small-scale magnetic fields. Firstly, it can act directly on the large-scale flow, which is known as “macro-feedback” (caused by the mean Lorentz force) and has been applied in several mean-field models (e.g., Schüssler 1979; Brandenburg et al. 1992). Secondly, it can affect the large-scale flow by affecting the convective motions and the best example of this is the magnetic quenching of the Λ effect, which is known as “micro-feedback” (for an application, see Küker et al. 1999).

To get an idea of these effects we measure the φ component of the Lorentz force (which appears in the zonal momentum equation) from the large-scale magnetic field $\mathrm{(}\overline{{J}}\mathrm{\times }\overline{{B}}{\mathrm{)}}_{\mathit{\phi }}$ and the small-scale contribution $\mathrm{(}\overline{{j}\mathrm{\times }{b}}{\mathrm{)}}_{\mathit{\phi }}$, which are shown in Fig. 15 both during magnetic maximum (left two panels) and minimum (right two panels) from Run A. We see that the small-scale Lorentz force, which enters into the total stress and thus the Λ effect, is typically stronger than the Lorentz force of the large-scale field and both have significant temporal variations, becoming weaker during magnetic minimum. Clearly both contributions are important, which was already emphasized by Beaudoin et al. (2013), who discussed in detail the differences in the integrated stresses for hydrodynamic and magnetic models in the case of SL rotation. Our results also show that both small-scale and large-scale contributions to the Lorentz force are responsible for producing temporal variations in the large-scale flows. Beaudoin et al. (2013) also emphasized the importance of magnetic fields in providing coupling to the radiative interior, which is absent in their hydrodynamic models. This is also the case in our magnetic models, which lack the presence of a lower overshoot layer. Finally, we note that while the Lorentz force from the mean field shows, during magnetic maximum, a systematic variation with distance from the axis, there are variations on similarly small scales both from the mean and fluctuating fields. This property is related to poor scale separation and may also apply to real stars.

Fig. 15 From Run A: contributions of the φ component of the large-scale Lorentz force $\overline{{J}}\mathrm{\times }\overline{{B}}$ and small-scale Lorentz force $\overline{{j}\mathrm{\times }{b}}$ during a magnetic maximum (left two panels) and minimum (right two panels). Values are given in units of 10-9 N m-3.

3.9. Turbulent angular momentum transport

Numerical simulations have been used on various occasions to study angular momentum transport in both hydrodynamic and magnetic cases, but they usually focus on the regime of SL rotation (Brun & Toomre 2002; Brun et al. 2004; Beaudoin et al. 2013). Furthermore, the contributions to the stress are often integrated over latitude or radius. In such representations, the stresses from opposite signs tend to cancel.

Fig. 16 Normalized profiles of Qrφ, Qθφ, νt, ΛV, and ΛH for Run A. The top and bottom panels show data averaged over four maxima and minima, respectively.

To make contact with mean-field theory invoking the Λ effect (Rüdiger 1980, 1989), it is useful to look at the profiles without integration over latitude or radius and to separate between diffusive and non-diffusive contributions. Simulations have confirmed many aspects of Λ effect in mean-field theory (Pulkkinen et al. 1993; Rieutord et al. 1994). The aim here is to interpret the differences between hydrodynamic and magnetic cases in the AS and SL regimes in terms of corresponding changes in the underlying Λ effect. For this purpose we first compute the contributions of the Reynolds stress ${\mathit{Q}}_{\mathit{ij}}\mathrm{=}\overline{{\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{u}}_{\mathit{j}}^{\mathrm{\prime }}}$, and the Maxwell stress ${\mathit{M}}_{\mathit{ij}}\mathrm{=}\mathrm{(}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}{\mathrm{)}}^{-1}\overline{{\mathit{B}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{B}}_{\mathit{j}}^{\mathrm{\prime }}}$ to the angular momentum balance in the convection zone. Here, primes denote fluctuating quantities, which are calculated by subtracting the longitudinal mean from the original quantity, e.g., ${\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{=}{\mathit{u}}_{\mathit{i}}\mathrm{-}{\overline{\mathit{u}}}_{\mathit{i}}$. The radial and latitudinal angular momentum transports are determined by the off-diagonal components of Q_ij and M_ij, namely Q_rφ, Q_θφ, M_rφ, and M_θφ, respectively. In the mean-field theory of hydrodynamics, the Reynolds stress contributions to angular momentum transport are approximated in terms of the turbulent viscosity ν_t and the Λ effect (Rüdiger 1980, 1989) $\begin{array}{ccc}{\mathit{Q}}_{\mathit{r\phi }}& \mathrm{=}& \\ {\mathit{Q}}_{\mathit{\theta \phi }}& \mathrm{=}& \end{array}$The coefficients Λ_V and Λ_H are the vertical and horizontal Λ effects, which are non-diffusive contributions to the Reynolds stress that arise from the interaction of anisotropic turbulence and rotation (Kitchatinov & Rüdiger 1995). As in earlier work (Käpylä et al. 2014), we adopt the mixing length formula to estimate ν_t${\mathit{\nu }}_{\mathrm{t}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{u}}_{\mathrm{rms}}{\mathit{\alpha }}_{\mathrm{MLT}}{\mathit{H}}_{\mathit{p}}\mathit{,}$(19)where α_MLT = 1.7 and H_p(r) = −(∂lnp/∂r)^-1. By computing u_rms = u_rms(r,θ) using φ averages, we get the profile of ν_t in the meridional plane.

The two leftmost panels of Fig. 16 show Q_rφ and Q_θφ normalized by ν_tΩ₀, and the third panel shows ν_t normalized by the micro-physical viscosity ν from Run A. We see that Q_rφ is negative in most of the convection zone, which implies inward transport of angular momentum. This is expected given the fact that the differential rotation in this case is AS. The negative Q_rφ is in agreement with Rieutord et al. (1994) and Käpylä et al. (2014). However for Run E, which produces SL differential rotation, Q_rφ is positive at low latitudes (Fig. 17). On the other hand, the latitudinal stress Q_θφ is positive (negative) in the northern (southern) hemisphere, which implies equatorward angular momentum transport. This is true for both Runs A and E (Figs. 16 and 17). The recent observation of Q_θφ (Hathaway et al. 2013) finds the same sign.

Fig. 17 Similar to Fig. 16, but for Run E and time-averaged over the last few maxima and minima.

The magnetic field is expected to change the total (Reynolds and Maxwell) stress over the activity cycle. In fact, magnetic quenching of the Reynolds stress is known to have important consequences in mean-field models, particularly in connection with explaining the origin of the torsional oscillation of the Sun or grand minima (e.g., Kitchatinov et al. 1999; Küker et al. 1999). Therefore, to see the variations over the activity cycle, we show Q_rφ and Q_θφ both during maximum and minimum phases. Note that these are not computed from one snapshot, but from averages over four maxima or minima phases. Comparing top and bottom panels of Fig. 16, we find that Q_rφ is slightly stronger during magnetic maximum, although Q_θφ is weaker.Even if we ignore the fact that our model is very far from that of Beaudoin et al. (2013), a quantitative comparison of our results is not straightforward. Moreover, Beaudoin et al. (2013) show the temporal variations of fluxes integrated across spherical shells; see their definition of I_r in Eq. (22) and Fig. 7. During magnetic maximum, our Q_rφ becomes stronger, whereas I_r in Beaudoin et al. (2013) becomes weaker and its overall variation is higher than ours. The temporal variations of Reynolds stresses show spatial coherence over all depths as seen in Beaudoin et al. (2013).

Next, we see in the middle panels of Fig. 16 that ν_t is significantly weaker during maximum because of magnetic quenching. The normalized value is around 40, which is similar to the value of Re. We note that ν_t decreases towards high latitudes, but increases towards deeper regions of the convection zone.

Fig. 18 From Run A: ${\mathit{R}}_{\mathit{r\phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{u}}}_{\mathit{r}}\mathrm{(}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{r}\sin \mathit{\theta }\mathrm{)}$, ${\mathit{R}}_{\mathit{\theta \phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{u}}}_{\mathit{\theta }}\mathrm{(}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{r}\sin \mathit{\theta }\mathrm{)}$ (left two panels), ${\mathit{M}}_{\mathit{r\phi }}\mathrm{=}\overline{{\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}{\mathit{B}}_{\mathit{\phi }}^{\mathrm{\prime }}}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$, ${\mathit{M}}_{\mathit{\theta \phi }}\mathrm{=}\overline{{\mathit{B}}_{\mathit{\theta }}^{\mathrm{\prime }}{\mathit{B}}_{\mathit{\phi }}^{\mathrm{\prime }}}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$ (middle two), ${\mathit{M}}_{\mathit{r\phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{B}}}_{\mathit{r}}{\overline{\mathit{B}}}_{\mathit{\phi }}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$ and ${\mathit{M}}_{\mathit{\theta \phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{B}}}_{\mathit{\theta }}{\overline{\mathit{B}}}_{\mathit{\phi }}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$ (right two), normalized by νtΩ⊙. Upper and lower panels are during magnetic maximum and minimum, respectively.

After solving Eqs. (17) and (18) for Λ_V and Λ_H, using the computed values of Q_rφ, Q_θφ, and ν_t, we find Λ_V and Λ_H again during magnetic maximum and minimum. These are shown in the last two panels of Fig. 16. We see that Λ_V is negative in most of the convection zone, which is in agreement with the findings from the first-order smoothing approximation (Kitchatinov & Rüdiger 1995, 2005), and with earlier numerical studies (e.g., Pulkkinen et al. 1993; Käpylä et al. 2004, 2014; Rüdiger et al. 2005; Käpylä & Brandenburg 2008). However, for Run E with SL rotation, Λ_V is positive at low latitudes and does not change much from maximum to minimum, which is why we show in Fig. 17 averages over all time. Another important property is that Λ_V is of the same order as ν_t, which is consistent with earlier findings (Käpylä et al. 2010a). We find that Λ_H is positive for both Runs A and E, which is again in agreement with the analytical results. Comparing the present results with the hydrodynamic counterparts of Käpylä et al. (2014, see their Fig. 10), we find that the Λ effect, is weaker. This could be a reason for getting strongly suppressed differential rotation, compared to the hydrodynamic simulations, particularly in AS cases (Runs A–BC). Furthermore we see a significant cycle-related variation in the Λ effect. In AS cases, both Λ_V and Λ_H are smaller during maximum (see last two columns in Fig. 16). This could be one of the reasons for significant temporal variations in the large-scale flow.

Fig. 19 Same as Fig. 18 but from Run E and averaged over last few maxima and minima.

Next we compute the radial and latitudinal components to the angular momentum transport by the meridional circulation, ${\mathit{R}}_{\mathit{r\phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{u}}}_{\mathit{r}}\mathrm{(}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{r}\sin \mathit{\theta }\mathrm{)}$ and ${\mathit{R}}_{\mathit{\theta \phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{u}}}_{\mathit{\theta }}\mathrm{(}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{r}\sin \mathit{\theta }\mathrm{)}$. We show these in Fig. 18 for Run A both, during magnetic maximum and minimum. Again, they are normalized by ν_tΩ_⊙ using the ν_t computed earlier. We see that these stresses, particularly the latitudinal component, are the most dominating stresses for transporting angular momentum. This is because in this run we have strong meridional circulation, mostly poleward near surface and equatorward near the bottom. We also observe temporal modulations of transport due to meridional circulation, which become stronger during magnetic minimum. Therefore this temporal variation could lead to a variation in the differential rotation. In fact, Beaudoin et al. (2013) found significant temporal variation in the angular momentum flux from meridional circulation and suggested that it could be the primary source of the torsional oscillations.The temporal variations of ${\mathit{R}}_{\mathit{r\phi }}^{\mathrm{LS}}$ are not well-correlated spatially, which is also seen in Beaudoin et al. (2013, see Fig. 7B).For Run E, which is shown in Fig. 19, we find smaller values of the stresses from the mean flows than in Run A, which is expected because of a much weaker meridional circulation in this case.

In the middle four panels of Fig. 18, we show the radial and latitudinal components of the Maxwell stress, ${\mathit{M}}_{\mathit{r\phi }}\mathrm{=}\overline{{\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}{\mathit{B}}_{\mathit{\phi }}^{\mathrm{\prime }}}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$ and ${\mathit{M}}_{\mathit{\theta \phi }}\mathrm{=}\overline{{\mathit{B}}_{\mathit{\theta }}^{\mathrm{\prime }}{\mathit{B}}_{\mathit{\phi }}^{\mathrm{\prime }}}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$, respectively, from Run A, both during magnetic maximum and minimum. We see that these are an order of magnitude smaller than Q_rφ and Q_θφ and they have the same signs as Q_rφ and Q_θφ, respectively. During magnetic maximum, both M_rφ and M_θφ (top middle panels of Fig. 18) are significantly larger than during magnetic minimum (lower middle panels). Therefore the Maxwell stresses may be another source of the temporal variation in the differential rotation.We recall that Beaudoin et al. (2013) also find large variations of the Maxwell stress in their simulation with highly correlated spatial variations.In Run E, shown in middle panels of Fig. 19, we again find similar values of M_rφ and M_θφ and acting in the opposite to Q_rφ and Q_θφ, respectively.

The Maxwell stress discussed above contributes to the redistribution of angular momentum by fluctuating (non-axisymmetric) magnetic fields. The mean axisymmetric field also contributes to the angular momentum balance (see, e.g., Brun et al. 2004; Beaudoin et al. 2013). This mean contribution is the result of correlations of the mean toroidal field with the radial and latitudinal components of the mean magnetic fields, which are ${\mathit{M}}_{\mathit{r\phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{B}}}_{\mathit{r}}{\overline{\mathit{B}}}_{\mathit{\phi }}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$ and ${\mathit{M}}_{\mathit{\theta \phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{B}}}_{\mathit{\theta }}{\overline{\mathit{B}}}_{\mathit{\phi }}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$. The two rightmost columns of panels in Fig. 18 show these two terms from Run A, both during magnetic maximum and minimum. For Run E, they are shown in the two rightmost panels of Fig. 19. Again we see that both components of the stress from the mean field are much smaller than the Reynolds stresses Q_rφ and Q_θφ and comparable to the fluctuating contributions, M_rφ and M_θφ. Importantly, they become much weaker during magnetic minimum.The variation in the large-scale magnetic tension is much larger than the variations in the other stresses and also much larger than what Beaudoin et al. (2013) found. They concluded that the variation in the large-scale magnetic tension is not responsible for the torsional oscillations. Comparing the top and bottom rows of the four rightmost panels of Fig. 18, we see that the temporal coherence of the spatial structures is not preserved during the activity cycle, which was also observed by Beaudoin et al. (2013). However, we must bear in mind that here we are comparing our Run A, which produces AS differential rotation, with a SL case in Beaudoin et al. (2013). In our SL cases, there are no prominent cycles and the temporal variations in the stresses are less pronounced than in Run A.

4. Conclusion

Motivated by the recently discovered bi-stability of stellar differential rotation in hydrodynamic spherical convection simulations (Gastine et al. 2014; Käpylä et al. 2014), we have performed several sets of magnetohydrodynamic simulations. Except for the allowance of dynamo-generated magnetic fields, our models are essentially the same as those of Käpylä et al. (2014). Taking different radiative conductivities, the convective velocities and hence the rotational influence are varied in the simulations. Runs A, B, and BC (Co = 1.34, 1.35, and 1.38) produce AS differential rotation, whereas Runs C, D, and E (Co = 1.44, 1.67, and 1.75) produce SL differential rotation. When we take an AS (SL) rotation profile as initial condition and perform a set of simulations by slowly increasing (decreasing) the radiative conductivity, we find similar states of differential rotation as in simulations that were started from scratch, i.e., with an initially rigid rotation profile. Therefore the bistable states of differential rotation seem to disappear in the MHD simulations, and we only find mono-stable solutions.

Besides the disappearance of the bistable differential rotation in MHD simulations, we find several other new results: (i) the abrupt transition from AS to SL rotation in hydrodynamic simulations now seems to become more gradual and this transition happens at slightly larger values of Ro_c (≈ 0.66) and thus smaller values of Co (≈ 1.4). This means that the magnetic field helps to produce SL differential rotation, which is also in agreement with the recent study of Fan & Fang (2014). (ii) The polar vortices or jet-like structures that were observed in previous hydrodynamic simulations (e.g., Heimpel & Aurnou 2007; Käpylä et al. 2014) are now absent. (iii) Both differential rotation and meridional circulation are now strongly suppressed in comparison to the hydrodynamic values and lie within the observed range. The suppression of these large-scale flows in the magnetic runs, particularly in AS differential rotation cases, could be a consequence of a reduction of the Λ effect and the presence of the Lorentz force of the magnetic field acting on the flow (Malkus & Proctor 1975). (iv) The large-scale flows show significant time variation as a consequence of the magnetic variations. All cases with AS differential rotation (Runs A–BC) show clear activity cycles and their large-scale flows have corresponding variations. Run A, which is more solar-like in terms of its highest convective flux and lowest Co value, shows ≈ 6% variation in $\overline{\mathrm{\Omega }}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}$ about its mean. However the variation in meridional circulation is as large as 60% and in large-scale shear, ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}$ and ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$, the variations are about 75% and 160%, respectively. All runs that produce SL differential rotation (Runs C–E) also show some magnetic variations, although they are not as regular and prominent as for AS differential rotation. In these runs we also see a detectable temporal variation in the large-scale flows, which is primarily caused by the variations of the large-scale (axisymmetric) magnetic tension, the Maxwell stresses, and the stress from the mean flow. For Run E, which has the smallest convective flux and SL differential rotation, the variation in $\overline{\mathrm{\Omega }}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}$ is about 4%, whereas in both meridional circulation and large-scale shear the variation is about 50%.

Many authors simulate the large-scale flows in stellar convection zones using hydrodynamical models assuming that the magnetic field does not have a significant effect (Brun & Toomre 2002; Ballot et al. 2007; Gastine et al. 2013, 2014; Guerrero et al. 2013; Käpylä et al. 2014, to mention just a few). It is difficult to quantify the effects of dynamo-generated magnetic fields on flows in the Sun from simulations as all the existing models are still far from the real Sun. However, from observations (Chou & Dai 2001; Hathaway & Rightmire 2010; Antia et al. 2008) we do see significant variations in both meridional circulation and rotational shear as well as a small variation in differential rotation in the form of torsional oscillations, which are believed to be (at least partially) coming from cyclic variations of the magnetic fields. The present study suggests that magnetic fields cannot be neglected in simulating the large-scale flows in the solar convection zone.

Acknowledgments

We thank the anonymous referee and Jörn Warnecke for a careful reading the paper and for suggestions, which improved the presentation. B.B.K. wishes to thank the University of Helsinki for hospitality during the initiation of this work. Financial support from the Academy of Finland grants No. 136189, 140970, 272786 (P.J.K.) and 272157 to the ReSoLVE Centre of Excellence (MJK), as well as the Swedish Research Council grants 621-2011-5076 and 2012-5797, the Research Council of Norway under the FRINATEK grant 231444, 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 computations have been carried out at the National Supercomputer Centres in Linköping and Umeå and the Center for Parallel Computers at the Royal Institute of Technology in Sweden, the Nordic High Performance Computing Center in Iceland, and the supercomputers hosted by CSC – IT Center for Science in Espoo, Finland.

References

All Tables

All Figures

Fig. 1 Contributions of different energy fluxes of Runs A (top) and E (bottom). The black (red) lines correspond to magnetohydrodynamic (hydrodynamic) case. Thin solid: radiative, dashed: convective, long dashed: viscous, dash-dotted: kinetic energy, dash-triple-dotted: SGS, and thick solid: total. Dotted horizontal lines indicate the zero and unity values. The vertical dotted line indicates the position of the middle of the convection zone r = rm.

In the text

Fig. 2 Variation of ${\mathrm{(}{\mathit{F}}_{\mathrm{conv}}\mathit{/}\mathit{\rho }{\mathit{u}}_{\mathrm{rms}}^{\mathrm{3}}\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$ averaged over the whole convection zone as a function of urms from different runs. Black asterisks: Runs A−E, red diamonds: Runs E1–E4, and blue triangles: Runs A1−A8. The dotted line shows the linear dependence between urms and ${\mathrm{(}{\mathit{F}}_{\mathrm{conv}}\mathit{/}\mathit{\rho }{\mathit{u}}_{\mathrm{rms}}^{\mathrm{3}}\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$.

In the text

	Fig. 3 Distribution of angular velocity $\overline{\mathrm{\Omega }}$ in the meridional plane from Runs A–E. The $\overline{\mathrm{\Omega }}$ is computed from Ω first by the longitudinal average and then the time average over the last few cycles. The arrows in the leftmost panel show the colatitudes at which the latitudinal differential rotation is computed in Eq. (13).
In the text

	Fig. 4 Same as Run B in Fig. 3, except here the latitudinal extent is different. For the left panel, the latitudinal boundaries are at ± 84°, whereas for the right panel they are at ± 66°.
In the text

	Fig. 5 Variations of angular velocities $\overline{\mathrm{\Omega }}$ at r = 0.96 R⊙ from Runs B (solid line), B′ (red dashed), and B′′ (blue dash-dotted).
In the text

	Fig. 6 Difference of $\overline{\mathrm{\Omega }}$ between the HD and MHD runs. The left and right panels are for Runs A and E, respectively.
In the text

Fig. 7 Radial (left panel) and latitudinal (right panel) differential rotation, defined by Eq. (13), for Runs A–E (green crosses), A0–A8 (red dashed line with triangles), and D0–D4 (black asterisks with dotted line). The other points are taken from the hydrodynamical simulations (black diamonds: Runs A–E, blue dotted line with asterisks: Runs D0–D4, and red dashed line with squares: Runs B0–B10) of Käpylä et al. (2014). The horizontal black dotted lines show the zero value. A red (black) arrow shows the transition point when differential rotation in the MHD (HD) simulations changes from AS to SL rotation with the decrease of Roc.

In the text

	Fig. 8 Meridional circulation from Runs A and E. The arrows show the direction of flow ${\overline{{u}}}_{\mathrm{m}}\mathrm{\equiv }\mathrm{\left(}\mathit{u̅}{}_{\mathit{r}}\mathit{,}\mathit{u̅}{}_{\mathit{\theta }}\mathrm{\right)}$ and the background colour shows uθ. Upper (lower) panels are from magnetohydrodynamic (hydrodynamic) simulations.
In the text

	Fig. 9 Butterfly diagrams: contours of the toroidal field ${\overline{\mathit{B}}}_{\mathit{\phi }}$ (in Gauss) at r = 0.74 R⊙ from Runs A, B, BC, C, D, and E (top to bottom).
In the text

	Fig. 10 Same as Fig. 9 but here ${\overline{\mathit{B}}}_{\mathit{\phi }}$ is shown at r = 0.96 R⊙.
In the text

	Fig. 11 Phase dispersion analysis results; from top to bottom, Runs A–E, in the middle column for the northern hemisphere, on the right for the southern hemisphere. In the left column, we show the original time series analysed (blue/red for north/south). The time series are normalized by their means, so values are not shown.
In the text

Fig. 12 Diagnostic diagrams showing a) ωcyc/ Ω0, b) the strength of the large-scale magnetic field over the whole convection zone measured by ${\overline{\mathit{B}}}_{\mathrm{rms}}\mathrm{=}\mathrm{⟨}\mathrm{⟨}{\mathit{B}}_{\mathit{r}}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}{\mathit{B}}_{\mathit{\theta }}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}{\mathit{B}}_{\mathit{\phi }}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{r\theta t}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ vs. log Co, and c) ωcyc/ Ω0 vs. ${\overline{\mathit{B}}}_{\mathrm{rms}}$. The red colour shows the values computed for the southern hemisphere.

In the text

Fig. 13 From Run A: a) the large-scale magnetic field over the whole convection zone $\overline{\mathit{B}}\mathrm{=}\mathrm{⟨}\mathrm{⟨}{\mathit{B}}_{\mathit{r}}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}{\mathit{B}}_{\mathit{\theta }}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}{\mathit{B}}_{\mathit{\phi }}{\mathrm{⟩}}_{\mathit{\phi }}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{r\theta }}^{\mathrm{1}\mathit{/}\mathrm{2}}$ normalized by Beq; b) the latitudinal component of meridional circulation uθ(r, ± 32°) (smoothed over 5 months) at r ≈ 0.95 R⊙ (black and red) and r ≈ 0.73 R⊙ (blue and green); c) azimuthally averaged angular velocity $\overline{\mathrm{\Omega }}\mathrm{(}\mathrm{0.95}{\mathit{R}}_{\mathrm{\odot }}\mathit{,}\mathrm{\pm }{\mathrm{32}}^{\mathrm{◦}}\mathrm{)}$;d) $\overline{\mathrm{\Omega }}\mathrm{\left(}\mathit{r,}{\mathrm{0}}^{\mathrm{◦}}\mathrm{\right)}$ at r = 0.73 R⊙ (red dashed) and r = 0.95 R⊙ (black).The dashed (solid) line corresponds to the southern (northern) hemisphere; e) radial shear ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{r}\mathrm{\right)}}$, and f) latitudinal shear ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$, defined in Eq. (13), as functions of time.

In the text

	Fig. 14 Same as Fig. 13, but from Run E, which produces SL differential rotation.
In the text

	Fig. 15 From Run A: contributions of the φ component of the large-scale Lorentz force $\overline{{J}}\mathrm{\times }\overline{{B}}$ and small-scale Lorentz force $\overline{{j}\mathrm{\times }{b}}$ during a magnetic maximum (left two panels) and minimum (right two panels). Values are given in units of 10-9 N m-3.
In the text

	Fig. 16 Normalized profiles of Qrφ, Qθφ, νt, ΛV, and ΛH for Run A. The top and bottom panels show data averaged over four maxima and minima, respectively.
In the text

	Fig. 17 Similar to Fig. 16, but for Run E and time-averaged over the last few maxima and minima.
In the text

Fig. 18 From Run A: ${\mathit{R}}_{\mathit{r\phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{u}}}_{\mathit{r}}\mathrm{(}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{r}\sin \mathit{\theta }\mathrm{)}$, ${\mathit{R}}_{\mathit{\theta \phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{u}}}_{\mathit{\theta }}\mathrm{(}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{0}}\mathit{r}\sin \mathit{\theta }\mathrm{)}$ (left two panels), ${\mathit{M}}_{\mathit{r\phi }}\mathrm{=}\overline{{\mathit{B}}_{\mathit{r}}^{\mathrm{\prime }}{\mathit{B}}_{\mathit{\phi }}^{\mathrm{\prime }}}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$, ${\mathit{M}}_{\mathit{\theta \phi }}\mathrm{=}\overline{{\mathit{B}}_{\mathit{\theta }}^{\mathrm{\prime }}{\mathit{B}}_{\mathit{\phi }}^{\mathrm{\prime }}}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$ (middle two), ${\mathit{M}}_{\mathit{r\phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{B}}}_{\mathit{r}}{\overline{\mathit{B}}}_{\mathit{\phi }}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$ and ${\mathit{M}}_{\mathit{\theta \phi }}^{\mathrm{LS}}\mathrm{=}{\overline{\mathit{B}}}_{\mathit{\theta }}{\overline{\mathit{B}}}_{\mathit{\phi }}\mathit{/}\mathit{\rho }{\mathit{\mu }}_{\mathrm{0}}$ (right two), normalized by νtΩ⊙. Upper and lower panels are during magnetic maximum and minimum, respectively.

In the text

	Fig. 19 Same as Fig. 18 but from Run E and averaged over last few maxima and minima.
In the text