Magnetically controlled stellar differential rotation near the transition from solar to antisolar profiles
^{1} NORDITA, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
email: bidyakarak@gmail.com
^{2} Department of Physics, Gustaf Hällströmin katu 2a ( PO Box 64), 00014 University of Helsinki, Finland
^{3} ReSoLVE Centre of Excellence, Department of Information and Computer Science, Aalto University, PO Box 15400, 00076 Aalto, Finland
^{4} Department of Astronomy, Stockholm University, 10691 Stockholm, Sweden
^{5} Tartu Observatory, 61602 Tõravere, Estonia
Received: 3 July 2014
Accepted: 21 January 2015
Context. Latetype stars rotate differentially owing to anisotropic turbulence in their outer convection zones. The rotation is called solarlike (SL) when the equator rotates fastest and antisolar (AS) otherwise. Hydrodynamic simulations show a transition from SL to AS rotation as the influence of rotation on convection is reduced, but the opposite transition occurs at a different point in the parameter space. The system is bistable, i.e., SL and AS rotation profiles can both be stable.
Aims. We study the effect of a dynamogenerated magnetic field on the largescale flows, particularly on the possibility of bistable behaviour of differential rotation.
Methods. We solve the hydromagnetic equations numerically in a rotating spherical shell that typically covers ± 75° latitude (wedge geometry) for a set of different radiative conductivities controlling the relative importance of convection. We analyse the resulting differential rotation, meridional circulation, and magnetic field and compare the corresponding modifications of the Reynolds and Maxwell stresses.
Results. In agreement with earlier findings, our models display SL rotation profiles when the rotational influence on convection is strong and a transition to AS when the rotational influence decreases. We find that dynamogenerated magnetic fields help to produce SL differential rotation compared to the hydrodynamic simulations. We do not observe any bistable states of differential rotation. In the AS cases we find coherent singlecell meridional circulation, whereas in SL cases we find multicellular patterns. In both cases, we obtain poleward circulation near the surface with a magnitude close to that observed in the Sun. In the slowly rotating cases, we find activity cycles, but no clear polarity reversals, whereas in the more rapidly rotating cases irregular variations are obtained. Moreover, both differential rotation and meridional circulation have significant temporal variations that are similar in strength to those of the Sun.
Conclusions. Purely hydrodynamic simulations of differential rotation and meridional circulation are shown to be of limited relevance as magnetic fields, selfconsistently generated by dynamo action, significantly affect the flows.
Key words: convection / turbulence / Sun: magnetic fields / Sun: rotation / stars: rotation
© 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 solarlike (SL) differential rotation. The opposite case where the equator rotates slower than the poles, is referred to as antisolar (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 mainsequence 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 – meanfield models and global convection simulations. In the meanfield approach, smallscale 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 firstorder smoothing, involve approximations that are illsuited for stellar conditions and may yield inaccurate results. A numerical method for determining the turbulent transport coefficients relevant for the electromotive force is the testfield 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 meanfield models are often based on educated guesses or they are even used as free parameters. Despite these shortcomings, hydrodynamical meanfield 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 meanfield models (e.g., Kitchatinov & Rüdiger 2004). At the same time, meanfield 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 BabcockLeighton process in the socalled 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 largescale 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 largescale contribution to the magnetic field (Brun et al. 2004).
According to meanfield hydrodynamics, differential rotation is generated from the anisotropy of the Reynolds stress, which is parameterised in terms of the socalled Λ 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 bistability 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 largescale effect (MalkusProctor effect; see Malkus & Proctor 1975; Brandenburg et al. 1992) and smallscale 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 largescale 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 largescale and smallscale contributions) to understand the temporal variations of the largescale 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_{0} ≤ r ≤ r_{1}, θ_{0} ≤ θ ≤ π − θ_{0}, and 0 ≤ φ ≤ φ_{0}, respectively. Here, r_{0} = 0.72 R_{⊙} and r_{1} = 0.97 R_{⊙} are the positions of the bottom and top of the computational domain, R_{⊙} is the radius of the Sun, θ_{0} = π/ 12 is the colatitude of the polar cap, and φ_{0} = π/ 2 is the longitudinal extent. The following hydromagnetic equations are solved: where A is the magnetic vector potential, B = ∇ × A is the magnetic field, J = ∇ × B/μ_{0} is the current density with μ_{0} 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^{3} is the gravitational acceleration with M_{⊙} being the mass of the Sun, Ω_{0} = (cosθ, − sinθ,0)Ω_{0} is the angular velocity vector, is the rate of strain tensor, where the semicolons denote covariant differentiation. The radiative and subgrid scale (SGS) heat fluxes are given by (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_{V}T is the specific internal energy.
2.2. Initial and boundary conditions
The initial hydrostatic state is isentropic, so the temperature is given by (6)where n_{ad} = 1.5 is the polytropic index and the value of ∂T/∂r at r = r_{0} 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 depthdependent polytropic index n(r) = δn(r/r_{0})^{15} + n_{ad} − δn for the radiative conductivity K(r) = K_{0} [n(r) + 1], where the reference conductivity is , with ℒ being the nondimensional 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 Ω_{0}, ν, η, and 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 in 0.75 R_{⊙}<r< 0.95 R_{⊙}, and 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_{1} it corresponds to 5 × 10^{8} m^{2} s^{1} in physical units. We also assume the density and temperature at r = r_{0} to have solar values, ρ_{0} = 200 kg m^{3} and T_{0} = 2.23 × 10^{6} 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^{4} = −K∇_{r}T − χ_{SGS}ρT∇_{r}s, where σ is related to the StefanBoltzmann 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 nonconvecting state the flux at the surface carries the total luminosity through the boundary. We use smallscale, lowamplitude 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^{6} 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^{3}. 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_{0}/L_{⊙})^{1/3}Ω_{⊙} and u_{sim} = (L_{0}/L_{⊙})^{1/3}u_{⊙}, where L_{0} is the luminosity in simulation, L_{⊙} ≈ 3.84 × 10^{26} 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 nondimensional and therefore not affected. All computations are performed with the Pencil code^{1}.
2.3. Dimensionless parameters and diagnostics
First, we define the nondimensional input parameters. The luminosity parameter ℒ and the normalized pressure scale height at the surface ξ are given by (7)where T_{1} is the temperature at the surface. The influence of rotation is measured by the Taylor number, (8)where Δr = r_{1} − r_{0} is the thickness of the convecting shell. The fluid, magnetic, and SGS Prandtl numbers are defined as (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 nondimensional viscosity, (10)and the Rayleigh and convective Rossby numbers (Gilman 1977), (11)where the entropy gradient of the nonconvecting hydrostatic solution is evaluated in the middle of the convection zone, r = r_{m}. We also quote the initial density contrast .
As diagnostic quantities we define the fluid and magnetic Reynolds numbers, and the Coriolis number as (12)where 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^{2}Re^{2}(k_{f}R_{⊙})^{4}.
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 bistability 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 Ω_{0} 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, dashdotted: kinetic energy, dashtripledotted: 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 = r_{m}. 

Open with DEXTER 
Fig. 2 Variation of averaged over the whole convection zone as a function of u_{rms} 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 u_{rms} and . 

Open with DEXTER 
In the statistically stationary state, the total luminosity L_{tot}(r) = 4πr^{2}ℱ_{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 rmsvelocity is compatible with a onethird 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).
Summary of the runs.
Fig. 3 Distribution of angular velocity in the meridional plane from Runs A–E. The 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). 

Open with DEXTER 
3.2. Differential rotation
In Fig. 3 we show the rotation profile 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 midlatitudes, 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 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 highlatitude 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 poleequator 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 nonmagnetic 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 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 nonmagnetic solarlike cases of Käpylä et al. (2014), we never obtain polar vortices or jetlike structures with magnetic fields. It is also interesting to note that the poleequator 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°. 

Open with DEXTER 
Fig. 5 Variations of angular velocities at r = 0.96 R_{⊙} from Runs B (solid line), B′ (red dashed), and B′′ (blue dashdotted). 

Open with DEXTER 
Fig. 6 Difference of between the HD and MHD runs. The left and right panels are for Runs A and E, respectively. 

Open with DEXTER 
3.3. Identifying the SL to AS transition
We measure the relative radial and latitudinal differential rotation by the quantities (13)where and are the equatorial rotation rates at the surface and at the base of the convection zone, and is the rotation rate at latitudes ± 55° computed as an average of 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 and , listed in Table 1, help us to identify AS and SL differential rotation. The SL differential rotation implies and . 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 bistability 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 Ro_{c}. 

Open with DEXTER 
In Fig. 7 we show differential rotation parameters and 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 as the difference of rotation between equator and the endpoints of the domain at θ = θ_{0} and π − θ_{0} as in Käpylä et al. (2013, 2014), instead of Ω_{55}, as we do here, we would have obtained smaller values of because of the slowly rotating highlatitude 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, . The sign of w_{3} 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 bistability
Next we study the flow bistability 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 and the background colour shows u_{θ}. Upper (lower) panels are from magnetohydrodynamic (hydrodynamic) simulations. 

Open with DEXTER 
Summary of the diagnostic quantities of Runs A–E.
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 multicell meridional circulation for SL differential rotation and single or twocell 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 nearsurface shear layer, which helps to produce a poleward flow in the upper layers through inward angular momentum transport possibly by the downflow 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 and 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 largescale, spatiotemporal organization of the magnetic field can be seen from a timelatitude or butterfly diagram of , 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 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 (in Gauss) at r = 0.74 R_{⊙} from Runs A, B, BC, C, D, and E (top to bottom). 

Open with DEXTER 
Fig. 10 Same as Fig. 9 but here is shown at r = 0.96 R_{⊙}. 

Open with DEXTER 
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 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 Ω_{0}, 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}/ Ω_{0} 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 (14)where f(t_{i}),i = 1,...,N is the input time series, σ^{2} 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 In the latter condition, Δt is the socalled correlation length. In our case, the time series is evenly sampled, in which case a modification to the general tophat 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. 

Open with DEXTER 
For the particular case when Δt is longer than the full data span, the D^{2}(P) statistics is essentially a slight reformulation of the wellknown 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 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^{2} 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 resamples 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, multiperiodicity 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}/ Ω_{0} 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 largescale magnetic field over the whole convection zone measured by vs. log Co, and c) ω_{cyc}/ Ω_{0} vs. . The red colour shows the values computed for the southern hemisphere. 

Open with DEXTER 
In Fig. 12b we show the mean magnetic field. Here the errors of are computed as the largest departure of the mean from any onethird 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}/ Ω_{0} increases with , 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 meanfield dynamo, so the increasing trend in the cycle frequency ratio suggests that the α effect increases with magnetic field strength, which is referred to as antiquenching. This interpretation hinges on some illknown 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 antiquenching itself. Antiquenching 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 largescale magnetic field () 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 , at r = 0.73 R_{⊙} and r = 0.95 R_{⊙}, as well as the latitudinal and radial differential rotation and , 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 and u_{θ}(0.95 R_{⊙}, ± 32°) ≈ −0.36, − 0.38.) This kind of weak anticorrelation 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 dynamogenerated 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 and u_{θ}(0.73 R_{⊙}, ± 32°) ≈ −0.22, − 0.47). We see that (Fig. 13c) also shows a weak anticorrelation with the magnetic variations (having correlation coefficient ≈ − 0.25). The strong magnetic fields during maxima change by a few per cent (≈ 6%). However (Fig. 13d) shows positive correlation (correlation coefficient ≈ 0.36) and the overall variation is larger (≈ 12%). Because of this variation of at the equator, the values of and (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 largescale magnetic field over the whole convection zone normalized by B_{eq}; 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 ;d) 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 , and f) latitudinal shear , defined in Eq. (13), as functions of time. 

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

Open with DEXTER 
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 largescale 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 anticorrelation 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 , , , , and , 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 largescale 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 largescale and smallscale magnetic fields. Firstly, it can act directly on the largescale flow, which is known as “macrofeedback” (caused by the mean Lorentz force) and has been applied in several meanfield models (e.g., Schüssler 1979; Brandenburg et al. 1992). Secondly, it can affect the largescale flow by affecting the convective motions and the best example of this is the magnetic quenching of the Λ effect, which is known as “microfeedback” (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 largescale magnetic field and the smallscale contribution , 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 smallscale Lorentz force, which enters into the total stress and thus the Λ effect, is typically stronger than the Lorentz force of the largescale 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 smallscale and largescale contributions to the Lorentz force are responsible for producing temporal variations in the largescale 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 largescale Lorentz force and smallscale Lorentz force during a magnetic maximum (left two panels) and minimum (right two panels). Values are given in units of 10^{9} N m^{3}. 

Open with DEXTER 
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 Q_{rφ}, Q_{θφ}, ν_{t}, Λ_{V}, and Λ_{H} for Run A. The top and bottom panels show data averaged over four maxima and minima, respectively. 

Open with DEXTER 
To make contact with meanfield 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 nondiffusive contributions. Simulations have confirmed many aspects of Λ effect in meanfield 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 , and the Maxwell stress 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., . The radial and latitudinal angular momentum transports are determined by the offdiagonal components of Q_{ij} and M_{ij}, namely Q_{rφ}, Q_{θφ}, M_{rφ}, and M_{θφ}, respectively. In the meanfield 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) The coefficients Λ_{V} and Λ_{H} are the vertical and horizontal Λ effects, which are nondiffusive 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}(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}Ω_{0}, and the third panel shows ν_{t} normalized by the microphysical 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 timeaveraged over the last few maxima and minima. 

Open with DEXTER 
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 meanfield 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: , (left two panels), , (middle two), and (right two), normalized by ν_{t}Ω_{⊙}. Upper and lower panels are during magnetic maximum and minimum, respectively. 

Open with DEXTER 
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 firstorder 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 cyclerelated 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 largescale flow.
Fig. 19 Same as Fig. 18 but from Run E and averaged over last few maxima and minima. 

Open with DEXTER 
Next we compute the radial and latitudinal components to the angular momentum transport by the meridional circulation, and . 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 are not wellcorrelated 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, and , 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 (nonaxisymmetric) 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 and . 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 largescale 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 largescale 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 bistability 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 dynamogenerated 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 monostable 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 jetlike 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 largescale 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 largescale 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 largescale flows have corresponding variations. Run A, which is more solarlike in terms of its highest convective flux and lowest Co value, shows ≈ 6% variation in about its mean. However the variation in meridional circulation is as large as 60% and in largescale shear, and , 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 largescale flows, which is primarily caused by the variations of the largescale (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 is about 4%, whereas in both meridional circulation and largescale shear the variation is about 50%.
Many authors simulate the largescale 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 dynamogenerated 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 largescale 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 62120115076 and 20125797, 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 HPCEuropa2 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
 Antia, H. M., Basu, S., & Chitre, S. M. 2008, ApJ, 681, 680 [NASA ADS] [CrossRef] [Google Scholar]
 Augustson, K., Brun, A. S., Miesch, M. S., & Toomre, J. 2014, ApJ, submitted [arXiv:1410.6547] [Google Scholar]
 Ballot, J., Brun, A. S., & TurckChièze, S. 2007, ApJ, 669, 1190 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Beaudoin, P., Charbonneau, P., Racine, E., & Smolarkiewicz, P. K. 2013, Sol. Phys., 282, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Moss, D., & Tuominen, I. 1992, A&A, 265, 328 [NASA ADS] [Google Scholar]
 Brandenburg, A., Jennings, R. L., Nordlund, Å., et al. 1996, J. Fluid Mech., 306, 325 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Brandenburg, A., Saar, S. H., & Turpin, C. R. 1998, ApJ, 498, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Chan, K. L., Nordlund, Å., & Stein, R. F. 2005, Astron. Nachr., 326, 681 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, T. M., ChristensenDalsgaard, J., Dziembowski, W. A., et al. 1989, ApJ, 343, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 2011, ApJ, 731, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., & Palacios, A. 2009, ApJ, 702, 1078 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., & Toomre, J. 2002, ApJ, 570, 865 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, K. L. 2010, in IAU Symp. 264, eds. A. G. Kosovichev, A. H. Andrei, & J.P. Rozelot, 219 [Google Scholar]
 Chatterjee, P., Mitra, D., Rheinhardt, M., & Brandenburg, A. 2011, A&A, 534, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chou, D.Y., & Dai, D.C. 2001, ApJ, 559, L175 [NASA ADS] [CrossRef] [Google Scholar]
 Choudhuri, A. R., Schüssler, M., & Dikpati, M. 1995, A&A, 303, L29 [NASA ADS] [Google Scholar]
 Cole, E., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2014, ApJ, 780, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Dikpati, M., & Charbonneau, P. 1999, ApJ, 518, 508 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., & Fang, F. 2014, ApJ, 789, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Gastine, T., Wicht, J., & Aurnou, J. M. 2013, Icarus, 225, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76 [NASA ADS] [CrossRef] [Google Scholar]
 Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A. 1977, Geophys. Astrophys. Fluid Dynam., 8, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A. 1983, ApJS, 53, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., & Rightmire, L. 2010, Science, 327, 1350 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., Upton, L., & Colegrove, O. 2013, Science, 342, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Hazra, G., Karak, B. B., & Choudhuri, A. R. 2014, ApJ, 782, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Heimpel, M., & Aurnou, J. 2007, Icarus, 187, 540 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., & Yokoyama, T. 2011, ApJ, 740, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 798, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., & Brandenburg, A. 2008, A&A, 488, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2004, A&A, 422, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2006, Astron. Nachr., 327, 884 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Brandenburg, A., Korpi, M. J., Snellman, J. E., & Narayan, R. 2010a, ApJ, 719, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., Brandenburg, A., Mitra, D., & Tavakol, R. 2010b, Astron. Nachr., 331, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011a, Astron. Nachr., 332, 883 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJ, 755, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karak, B. B. 2010, ApJ, 724, 1021 [NASA ADS] [CrossRef] [Google Scholar]
 Karak, B. B., & Choudhuri, A. R. 2011, MNRAS, 410, 1503 [NASA ADS] [Google Scholar]
 Karak, B. B., & Choudhuri, A. R. 2012, Sol. Phys., 278, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Karak, B. B., & Choudhuri, A. R. 2013, Res. Astron. Astrophys., 13, 1339 [NASA ADS] [CrossRef] [Google Scholar]
 Karak, B. B., Kitchatinov, L. L., & Choudhuri, A. R. 2014, ApJ, 791, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Kholikov, S., Serebryanskiy, A., & Jackiewicz, J. 2014, ApJ, 784, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Olemskoy, S. V. 2011, MNRAS, 411, 1059 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 1995, A&A, 299, 446 [NASA ADS] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 2004, Astron. Nachr., 325, 496 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 2005, Astron. Nachr., 326, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., Rüdiger, G., & Küker, M. 1994, A&A, 292, 125 [NASA ADS] [Google Scholar]
 Kitchatinov, L. L., Pipin, V. V., Makarov, V. I., & Tlatov, A. G. 1999, Sol. Phys., 189, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Kővári, Z., Kriskovics, L., Künstler, A., et al. 2015, A&A, 573, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Küker, M., & Rüdiger, G. 2011, Astron. Nachr., 332, 933 [NASA ADS] [CrossRef] [Google Scholar]
 Küker, M., Arlt, R., & Rüdiger, G. 1999, A&A, 343, 977 [NASA ADS] [Google Scholar]
 Lindborg, M., Mantere, M. J., Olspert, N., et al. 2013, A&A, 559, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malkus, W. V. R., & Proctor, M. R. E. 1975, J. Fluid Mech., 67, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Matt, S. P., Do Cao, O., Brown, B. P., & Brun, A. S. 2011, Astron. Nachr., 332, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., & Dikpati, M. 2014, ApJ, 785, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., & Hindman, B. W. 2011, ApJ, 743, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., & Toomre, J. 2009, Ann. Rev. Fluid Mech., 41, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 641, 618 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2013, ApJ, 762, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, A., Brandenburg, A., Jennings, R. L., et al. 1992, ApJ, 392, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Passos, D., Charbonneau, P., & Beaudoin, P. 2012, Sol. Phys., 279, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Pelt, J. 1983, in Statistical Methods in Astronomy, ed. E. J. Rolfe, ESA SP, 201, 37 [Google Scholar]
 Pipin, V. V., & Kosovichev, A. G. 2011, ApJ, 727, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Pulkkinen, P., Tuominen, I., Brandenburg, A., Nordlund, A., & Stein, R. F. 1993, A&A, 267, 265 [NASA ADS] [Google Scholar]
 Racine, É., Charbonneau, P., Ghizaru, M., Bouchat, A., & Smolarkiewicz, P. K. 2011, ApJ, 735, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2006, ApJ, 647, 662 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., Brandenburg, A., Mangeney, A., & Drossart, P. 1994, A&A, 286, 471 [NASA ADS] [Google Scholar]
 Rüdiger, G. 1980, Geophys. Astrophys. Fluid Dyn., 16, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G. 1989, Differential Rotation and Stellar Convection. Sun and Solartype Stars (Berlin: Akademie Verlag) [Google Scholar]
 Rüdiger, G., Egorov, P., & Ziegler, U. 2005, Astron. Nachr., 326, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Saar, S. H., & Brandenburg, A. 1999, ApJ, 524, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Schad, A., Timmer, J., & Roth, M. 2013, ApJ, 778, L38 [NASA ADS] [CrossRef] [Google Scholar]
 Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Schrinner, M., Rädler, K.H., Schmitt, D., Rheinhardt, M., & Christensen, U. 2005, Astron. Nachr., 326, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Schrinner, M., Rädler, K.H., Schmitt, D., Rheinhardt, M., & Christensen, U. R. 2007, Geophys. Astrophys. Fluid Dyn., 101, 81 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Schüssler, M. 1979, A&A, 72, 348 [NASA ADS] [Google Scholar]
 Stellingwerf, R. F. 1978, ApJ, 224, 953 [NASA ADS] [CrossRef] [Google Scholar]
 Strassmeier, K. G., Kratzwald, L., & Weber, M. 2003, A&A, 408, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thompson, M. J., ChristensenDalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Tobias, S. M., Brummell, N. H., Clune, T. L., & Toomre, J. 2001, ApJ, 549, 1183 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2013, ApJ, 778, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Weber, M., Strassmeier, K. G., & Washuettl, A. 2005, Astron. Nachr., 326, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, J., Bogart, R. S., Kosovichev, A. G., Duvall, Jr., T. L., & Hartlep, T. 2013, ApJ, 774, L29 [NASA ADS] [CrossRef] [Google Scholar]
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, dashdotted: kinetic energy, dashtripledotted: 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 = r_{m}. 

Open with DEXTER  
In the text 
Fig. 2 Variation of averaged over the whole convection zone as a function of u_{rms} 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 u_{rms} and . 

Open with DEXTER  
In the text 
Fig. 3 Distribution of angular velocity in the meridional plane from Runs A–E. The 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). 

Open with DEXTER  
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°. 

Open with DEXTER  
In the text 
Fig. 5 Variations of angular velocities at r = 0.96 R_{⊙} from Runs B (solid line), B′ (red dashed), and B′′ (blue dashdotted). 

Open with DEXTER  
In the text 
Fig. 6 Difference of between the HD and MHD runs. The left and right panels are for Runs A and E, respectively. 

Open with DEXTER  
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 Ro_{c}. 

Open with DEXTER  
In the text 
Fig. 8 Meridional circulation from Runs A and E. The arrows show the direction of flow and the background colour shows u_{θ}. Upper (lower) panels are from magnetohydrodynamic (hydrodynamic) simulations. 

Open with DEXTER  
In the text 
Fig. 9 Butterfly diagrams: contours of the toroidal field (in Gauss) at r = 0.74 R_{⊙} from Runs A, B, BC, C, D, and E (top to bottom). 

Open with DEXTER  
In the text 
Fig. 10 Same as Fig. 9 but here is shown at r = 0.96 R_{⊙}. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 12 Diagnostic diagrams showing a) ω_{cyc}/ Ω_{0}, b) the strength of the largescale magnetic field over the whole convection zone measured by vs. log Co, and c) ω_{cyc}/ Ω_{0} vs. . The red colour shows the values computed for the southern hemisphere. 

Open with DEXTER  
In the text 
Fig. 13 From Run A: a) the largescale magnetic field over the whole convection zone normalized by B_{eq}; 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 ;d) 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 , and f) latitudinal shear , defined in Eq. (13), as functions of time. 

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

Open with DEXTER  
In the text 
Fig. 15 From Run A: contributions of the φ component of the largescale Lorentz force and smallscale Lorentz force during a magnetic maximum (left two panels) and minimum (right two panels). Values are given in units of 10^{9} N m^{3}. 

Open with DEXTER  
In the text 
Fig. 16 Normalized profiles of Q_{rφ}, Q_{θφ}, ν_{t}, Λ_{V}, and Λ_{H} for Run A. The top and bottom panels show data averaged over four maxima and minima, respectively. 

Open with DEXTER  
In the text 
Fig. 17 Similar to Fig. 16, but for Run E and timeaveraged over the last few maxima and minima. 

Open with DEXTER  
In the text 
Fig. 18 From Run A: , (left two panels), , (middle two), and (right two), normalized by ν_{t}Ω_{⊙}. Upper and lower panels are during magnetic maximum and minimum, respectively. 

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

Open with DEXTER  
In the text 