Issue 
A&A
Volume 570, October 2014



Article Number  A43  
Number of page(s)  10  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201423412  
Published online  14 October 2014 
Confirmation of bistable stellar differential rotation profiles
^{1} Department of Physics, Gustaf Hällströmin katu 2a ( PO Box 64), 00014 University of Helsinki, Finland
email: petri.kapyla@helsinki.fi
^{2} ReSoLVE Centre of Excellence, Department of Information and Computer Science, Aalto University, PO Box 15400, 00076 Aalto, Finland
^{3} NORDITA, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
^{4} Department of Astronomy, Stockholm University, 10691 Stockholm, Sweden
Received: 13 January 2014
Accepted: 18 June 2014
Context. Solarlike differential rotation is characterized by a rapidly rotating equator and slower poles. However, theoretical models and numerical simulations can also result in a slower equator and faster poles when the overall rotation is slow.
Aims. We study the critical rotational influence under which differential rotation flips from solarlike (fast equator, slow poles) to an antisolar one (slow equator, fast poles). We also estimate the nondiffusive (Λ effect) and diffusive (turbulent viscosity) contributions to the Reynolds stress.
Methods. We present the results of threedimensional numerical simulations of mildly turbulent convection in spherical wedge geometry. Here we apply a fully compressible setup which would suffer from a prohibitive time step constraint if the real solar luminosity was used. To avoid this problem while still representing the same rotational influence on the flow as in the Sun, we increase the luminosity by a factor of roughly 10^{6} and the rotation rate by a factor of 10^{2}. We regulate the convective velocities by varying the amount of heat transported by thermal conduction, turbulent diffusion, and resolved convection.
Results. Increasing the efficiency of resolved convection leads to a reduction of the rotational influence on the flow and a sharp transition from solarlike to antisolar differential rotation for Coriolis numbers around 1.3. We confirm the recent finding of a largescale flow bistability: contrasted with running the models from an initial condition with unprescribed differential rotation, the initialization of the model with certain kind of rotation profile sustains the solution over a wider parameter range. The antisolar profiles are found to be more stable against perturbations in the level of convective turbulent velocity than the solartype solutions.
Conclusions. Our results may have implications for real stars that start their lives as rapid rotators implying solarlike rotation in the early mainsequence evolution. As they slow down, they might be able to retain solarlike rotation for lower Coriolis numbers, and thus longer in time, before switching to antisolar rotation. This could partially explain the puzzling findings of antisolar rotation profiles for models in the solar parameter regime.
Key words: convection / turbulence / Sun: rotation / stars: rotation
© ESO, 2014
1. Introduction
The solar surface differential rotation is characterized by a fast equator and a monotonic decrease of angular velocity toward the poles. The internal rotation of the Sun is such that only weak radial shear is present in the bulk of the convection zone. Strong shear is present only in the boundary layers at the bottom in the tachocline and in the nearsurface shear layer (e.g. Thompson et al. 2003; Miesch & Toomre 2009, and references therein). The differential rotation is explained by the interaction of anisotropic turbulence and rotation, which is represented by a nondiffusive contribution to the Reynolds stress known as Λ effect in meanfield hydrodynamics (Krause & Rüdiger 1974; Rüdiger 1980, 1989). Furthermore, turbulent latitudinal heat fluxes, or stably stratified layers below (Brun et al. 2011) or above (Warnecke et al. 2013) the convection zone are needed to explain why the TaylorProudman balance is broken in the Sun. Meanfield models with these ingredients or other parameterizations have been successful in reproducing the solar rotation profile (e.g. Brandenburg et al. 1992; Kitchatinov & Rüdiger 1995; Rempel 2005; Hotta & Yokoyama 2011), and recently also the tachocline and the nearsurface shear layer (e.g. Kitchatinov & Olemskoy 2011).
Meanfield models rarely produce antisolar differential rotation unless a strong meridional circulation is imposed (Kitchatinov & Rüdiger 2004). This suggestion is supported by simulations of Dobler et al. (2006), which displayed antisolar differential rotation in their fully convective models as a consequence of strong meridional circulation. An arguably similar pathway was offered by Aurnou et al. (2007), who proposed that strong mixing by turbulent convection would be the primary agent for angular momentum equilibration and thus antisolar differential rotation. This mixing must then be stronger in the meridional plane than in the azimuthal direction to produce negative differential rotation, similar to what is expected to occur in the nearsurface shear layer of the Sun (Kitchatinov & Rüdiger 2005). Differential rotation and meridional circulation are intimately coupled (Kippenhahn 1963), and one implies the other. Antisolar differential rotation implies a counterclockwise circulation in the northern hemisphere, corresponding to poleward motion at the surface. In the context of the solar nearsurface shear layer, this mechanism is referred to as gyroscopic pumping (Miesch & Hindman 2011).
Both meanfield models and simulations have limitations. Furthermore, approximations like firstorder smoothing are used to derive turbulent transport coefficients, such as turbulent viscosity and the socalled Λ effect. However, their range of validity in stellar convection zones is questionable and only limited comparisons between hydrodynamic meanfield models and direct simulations have been performed (see, however Rieutord et al. 1994). Direct numerical simulations of convection, on the other hand, lack smallscale structure and tend to be dominated by structures that span the entire convection zone.
The rotational influence on the flow can be measured by the local Coriolis number Co = 2Ωτ, where Ω is the rotation rate and τ is the turnover time. It is large in the bulk of the solar convection zone, especially if τ is estimated from mixing length theory, which predicts values of Co ranging from 10^{3} near the surface to more than 10 in the deep layers (e.g. Ossendrijver 2003; Brandenburg & Subramanian 2005; Käpylä 2011). This is not captured by present simulations, which might well be due to insufficient density stratification in most numerical simulations so far. Whether the deep layers of stellar convection zones are really like this remains an open question. Observational evidence also points to antisolar differential rotation in some slowly rotating stars (e.g. Strassmeier et al. 2003; Weber et al. 2005) and solarlike rotation in rapidly rotating dwarfs (e.g. Collier Cameron 2002). However, for the star with the best observational evidence for antisolar differential rotation, namely the single K giant star HD 31993 (Strassmeier et al. 2003), the Coriolis number is around unity and thus close to the expected transition.
Recent numerical studies have explored the transition from antisolar to solarlike differential rotation (Brun & Palacios 2009; Chan 2010; Käpylä et al. 2011b,a; Gastine et al. 2013, 2014; Guerrero et al. 2013). Here we concentrate on the effect of superabiabaticity on the rotational influence and resulting largescale flows, and study the bistability of the differential rotation first reported by Gastine et al. (2014) using Boussinesq convection. They found that, near the transition from antisolar to solarlike differential rotation, both kinds of solutions are possible – depending on the initial conditions. This is interesting for stellar applications because most stars rotate rapidly when they are born and slow down due to magnetic braking. Rapid rotation implies solarlike differential rotation, which might then persist despite the angular momentum loss due to stellar winds.
In this paper we set out to study the transition of solarlike rotation profiles into the antisolar regime, extending the work of Gastine et al. (2014) into models of compressible convection. Using the spherical wedge model employed in various previous papers, the essential parts summarized in Sect. 2, we investigate the effect of changing the rotational influence by modifying the radiative conductivity that has an effect on the convective velocities. We perform two types of simulations, presented in Sect. 3. Firstly, we run models from scratch, i.e. without an initially prescribed rotation profile, and locate the solar to antisolar transition in terms of Coriolis number. Secondly, we investigate the dependence of the solutions on initial conditions by running models from solar and antisolar states, with otherwise identical parameters.
2. The Model
Our hydrodynamic model is essentially the same as the one used in Käpylä et al. (2011a). The same model has also been used to model convectiondriven dynamos (Käpylä et al. 2012, 2013; Cole et al. 2014). The computational domain is a wedge in spherical polar coordinates, where (r,θ,φ) are radius, colatitude, and longitude. The radial, latitudinal, and longitudinal extents of the wedge are r_{0} ≤ r ≤ r_{1}, θ_{0} ≤ θ ≤ π − θ_{0}, and 0 ≤ φ ≤ φ_{0}, respectively, where r_{0} = 0.72 R_{⊙} and r_{1} = 0.97 R_{⊙} denote the positions of the bottom and top of the computational domain, and R_{⊙} = 7 × 10^{8} m is the radius of the Sun. Here we consider θ_{0} = π/ 12 and φ_{0} = π/ 2, so we cover a quarter of the azimuthal extent between ± 75° latitude. The dependence on the latitudinal extent of the wedge was studied by Käpylä et al. (2011b), who found that the results are robust as long as the opening angle of the wedge is more than 90 degrees. We solve the compressible hydrodynamic equations,
where D/Dt = ∂/∂t + u·∇ is the advective time derivative, ρ is the density, ν is the constant kinematic viscosity, (4)are radiative and subgridscale (hereafter SGS) heat fluxes, where K is the radiative heat conductivity and χ_{SGS} is the turbulent heat conductivity, which represents the unresolved convective transport of heat (Käpylä et al. 2013) and was referred to as χ_{t} in Käpylä et al. (2011a, 2012). Furthermore, s is the specific entropy, T is the temperature, and p is the pressure. The fluid obeys the ideal gas law with p = (γ − 1)ρe, where γ = c_{P}/c_{V} = 5/3 is the ratio of specific heats at constant pressure and volume, respectively, and e = c_{V}T is the specific internal energy. The rate of strain tensor S is given by (5)where the semicolons denote covariant differentiation (Mitra et al. 2009).
The gravitational acceleration is given by g = −GM_{⊙}r/r^{3}, where G = 6.67 × 10^{11} m^{3} kg^{1} s^{2} is the gravitational constant, and M_{⊙} = 2.0 × 10^{30} kg is the mass of the Sun. We neglect selfgravity of the matter in the convection zone. Furthermore, the rotation vector Ω_{0} is given by Ω_{0} = (cosθ, − sinθ,0)Ω_{0}.
2.1. Initial and boundary conditions
Here we make an effort to connect the model more closely with the parameters of the Sun. Due to the fully compressible formulation of our model, we are faced with a prohibitive time step limitation if we were to use the solar luminosity. As explained in detail in Käpylä et al. (2013), we circumvent this by using a roughly 10^{6} times higher luminosity in the model in comparison to the Sun. As the convective energy flux scales as F_{conv} ~ ρu^{3}, the convective velocity u is roughly 100 times greater in the simulations than in the Sun. To obtain the same rotational influence on the flow as in the Sun, we must therefore increase Ω by the same factor. In general, this can be written as (6)where L_{ratio} = L_{0}/L_{⊙}, with L_{0} and L_{⊙} ≈ 3.84 × 10^{26} W being the luminosities of the model and the Sun, respectively, and Ω_{⊙} ≈ 2.7 × 10^{6} s^{1} is the mean solar rotation rate, corresponding to 430 nHz. In what follows we scale our results back to solar units so that, say for the velocity, we quote . The scaling used here is based on dimensional arguments. It is supported by mixing length theory (Vitense 1953) and simulations (Brandenburg et al. 2005; Miesch et al. 2012), and should be applicable as long as the energy transport is not yet affected by rotation (see e.g. Yadav et al. 2013). Furthermore, we assume that the density and the temperature at the base of the convection zone at r = 0.72 R_{⊙} have the solar values ρ_{0} = 200 kg m^{3} and T_{0} = 2.23×10^{6} K.
The initial state is isentropic and the hydrostatic temperature gradient is given by (7)where n_{ad} = 1.5 is the polytropic index for an adiabatic stratification. We fix the value of ∂T/∂r on the lower boundary. The density profile follows from hydrostatic equilibrium. To speed up the thermal relaxation, the initial condition is chosen not to be in thermodynamic equilibrium, but closer to the final convecting state. We choose the heat conduction profile such that radiative diffusion is responsible for supplying the energy flux in the system and progressively less so further out by choosing a radiative conductivity, K(r) = K_{0} [n(r) + 1], with n(r) = δn(r/r_{0})^{15} + n_{ad} − δn replacing the polytropic index, (8)being a reference conductivity, and ℒ being the nondimensional luminosity, given below. Now n = n_{ad} at the bottom of the convection zone and approaches n_{ad} − δn at the surface. This means that K = (n + 1)K_{0} decreases toward the surface like r^{15} such that the value of δn regulates the flux that is carried by convection (Brandenburg et al. 2005). Initial, final, and hydrostatic profiles of the temperature and density as well as the profiles of Pr_{SGS} = ν/χ_{SGS} and Pr = ν/χ, where χ = K/ρc_{P}, are shown in Fig. 1. We introduce weak smallscale Gaussian noise velocity perturbations in the initial state.
Our simulations are defined by the energy flux imposed at the bottom boundary, F_{b} = −(K∂T/∂r)  _{r = r0} as well as the values of Ω_{0}, ν, and . Furthermore, the radial profile of χ_{SGS} is piecewise constant above r> 0.75 R_{⊙} with at 0.75R_{⊙}<r< 0.95 R_{⊙}, and above r = 0.95 R_{⊙}. Below r = 0.75 R_{⊙}, χ_{SGS} tends smoothly to zero; see Fig. 1 of Käpylä et al. (2011a). We fix the value of χ_{SGS} such that it corresponds to 5 × 10^{8} m^{2} s^{1} in physical units at r = r_{1}.
The radial and latitudinal boundaries are assumed to be impenetrable and stress free, i.e., Density and specific entropy have vanishing first derivatives on the latitudinal boundaries, thus suppressing heat fluxes through them.
On the outer radial boundary we apply a black body condition (11)where σ is the StefanBoltzmann constant. We use a modified value for σ that takes into account that both surface temperature and energy flux through the domain are larger than in the Sun. We choose σ such that the flux at the surface, σT^{4}, carries the total luminosity through the boundary in the initial nonconvecting state.
Fig. 1 Left and middle panels: temperature T in units of 10^{6} K and density ρ in units of kg m^{3}, respectively, for the thermally relaxed state (black solid line), initial condition (red dotted), and hydrostatic solutions (blue dashed). Rightmost panel: Prandtl numbers related to radiative diffusion Pr = ν/χ where χ = K/ρc_{P} (black solid line), and the turbulent heat conductivity Pr_{SGS} = ν/χ_{SGS} (blue dashed) as functions of radius. Data is taken from Run D. 
2.2. Dimensionless parameters
The nondimensional input parameters of our models are the luminosity parameter (12)the normalized pressure scale height at the surface, (13)with T_{1} being the temperature at the surface, the Taylor number (14)where Δr = r_{1} − r_{0} = 0.25 R_{⊙}, as well as the fluid and SGS Prandtl numbers (15)where χ_{m} = K(r_{m}) /c_{P}ρ_{m} is the thermal diffusivity and ρ_{m} is the density, both evaluated at r = r_{m} = 0.845 R_{⊙}. We vary Pr and keep Pr_{SGS} = 0.25 fixed. Finally, the nondimensional viscosity is (16)In addition to ξ, we quote the initial density contrast, . In the current moderately stratified simulations the density contrast changes by less than 10 per cent during the run, see the middle panel of Fig. 1.
All other parameters are used as diagnostics and are not input parameters. These include the fluid Reynolds number (17)where is an estimate of the wavenumber of the largest eddies. The Coriolis number is defined as (18)where is the rms velocity and the subscripts indicate averaging over r, θ, φ, and a time interval during which the run is thermally relaxed. For u_{rms} we omit the contribution from the azimuthal velocity, because it is dominated by the differential rotation (Käpylä et al. 2011b). To have a reasonable estimate of the rms velocity similar to that under isotropic conditions, we compensate for the omission of by the 3/2 factor. The Taylor number can also be written as Ta = Co^{2}Re^{2}(k_{f}R_{⊙})^{4}. Furthermore, we define the Rayleigh number as (19)where s_{hs} is the entropy in the hydrostatic, nonconvecting state. We compute the hydrostatic stratification by evolving a onedimensional model (no convection) with the values and profiles of K and χ_{SGS} given above. We also quote the convective Rossby number (Gilman 1977) (20)We define mean quantities as averages over the φcoordinate and denote them by overbars. We also often average the data in time over the period of the simulations where thermal energy and differential rotation have reached statistically saturated states.
The simulations are performed with the Pencil Code^{1}, which uses a highorder finite difference method for solving the compressible equations of magnetohydrodynamics.
Fig. 2 Timeaveraged energy fluxes from Runs A (top) and E (bottom): radiative (thin solid line), convective (dashed), kinetic energy (dashdotted), SGS (tripledashdotted), and viscous (longdashed) flux. The thick solid line denotes the total flux, whereas the red horizontal dotted lines show the zero and unity line. The red vertical dotted line at r = r_{m} shows the midpoint of the convection zone. 
Fig. 3 Radial dependence of the time averaged fluctuating rmsvelocity where the contributions from the mean flows are omitted for Runs A (solid black line), C (blue dashed), and E (red dotdashed). The black dotted line shows the convective velocity from the mixing length (ML) model of Stix (2002). 
3. Results
Our simulations are summarized in Table 1. We perform three sets of runs. In the first set, we run the model from the initial conditions described in Sect. 2.1 (Runs A–E). Here, Runs A–C turn out to have antisolar differential rotation, while Runs D and E have solarlike differential rotation. Secondly, we study the bistability of the rotation profile by taking either a solarlike (Runs D0–D4) or an antisolar (Runs B0–B10b) solution as initial conditions. Apart from δn, we keep all other parameters fixed. We also estimate Λ effect coefficients from the simulation results.
Summary of the runs.
3.1. Effect of varying radiative flux
We change the radiative conductivity by varying the parameter δn, which regulates the amount of flux that convection has to transport, thus influencing the convective velocities and the Coriolis number. The different contributions to the total energy flux from Runs A and E are shown in Fig. 2. The definitions of the fluxes can be found in Käpylä et al. (2013). We give the fractions of convective and radiative contributions to the total energy flux at the middle of the convection zone in Table 1. For δn = 2.5 (Run A) the convective flux can exceed the total flux, so the radiative flux transports less than 10 per cent of the luminosity in the upper part of the convection zone. In Run E with δn = 1.75 the fractions of radiative diffusion and convection are 37 and 53 per cent, respectively. In the extreme case of δn = 1 (Runs B10 and B10b), convection transports only about 20 per cent of the flux. These cases are comparable to the setups used in earlier works (Käpylä et al. 2010b, 2011b).
The rmsvelocities based on the fluctuating velocity for Runs A, C, and E are shown in Fig. 3. The changes between Runs A and C are rather subtle, the main effect being the decrease in the overall magnitude of u_{rms}. The main difference between Runs C and E is the increased contrast of the values near the boundaries. We also find that all runs show higher velocities than the mixing length (ML) model of Stix (2002) below r = 0.95 R_{⊙}. The increase in u_{rms} near the lower boundary is due to our impenetrable boundary condition. Near the upper boundary our values of u_{rms} are close to those obtained from ML. A likely explanation is the weaker density stratification used in our simulations (Γ ≈ 12) as opposed to what is realized in the ML model (Γ ≈ 60), leading to higher values near the surface in the latter. Computing an average velocity using the data in Fig. 3 and using that to compute a Coriolis number according to Eq. (18), we find Co ≈ 2.13 for the ML model. This is clearly above the values our Runs A–E owing to the lower velocities. However, we note that recent highresolution simulations of nonrotating solar convection suggest significantly higher velocities than anticipated from ML models (Hotta et al. 2014) so the actual Coriolis number of the Sun might also be much smaller than our estimate.
Fig. 4 Timeaveraged rotation profiles from Runs A–E showing in nHz. 
Fig. 5 Timeaveraged meridional circulation from Runs A (left) and D (right). The arrows show the flow , whereas the colour contours show . 
3.2. Differential rotation
The decrease in the rmsvelocity, since δn is lowered, is also reflected in the Coriolis number, which more than doubles between the extreme cases A and B10. The increasing rotational effect on the flow affects the rotation profile realized in the runs. The profiles of are shown in Fig. 4 for Runs A–E, which were run from the initial conditions stated in Sect. 2.1 with δn varying systematically. We characterize the relative radial and latitudinal differential rotation by the quantities (Käpylä et al. 2013) (21)where and are the equatorial rotation rates at the surface and at the base of the convection zone, and is the average rotation rate between the latitudinal boundaries on the outer radius. In the Sun, and , which is what is required for a solution to be classified as “solarlike”.
Fig. 6 Radial (left panel) and latitudinal (right panel) differential rotation, defined by Eqs. (21), from Runs A–E (diamonds), and Set D (blue dotted line with asterisks) and B (red dashed line with triangles). 
For high values of δn, i.e. for low radiative flux (Runs A–C), the rotation profile is antisolar with a negative radial gradient of at all latitudes. The difference between the equator and latitudes ± 75° is larger than 1.5Ω_{⊙} in both cases. In Run D the rotation profile flips to solarlike. Thus, the transition from the antisolar to solarlike regime occurs when 1.85 <δn< 2.0, corresponding to 1.33 < Co < 1.50 in this set of runs. This is compatible with the results of Gastine et al. (2014), who found the transition at a local Rossby number of around Ro_{l} ≈ 1 where Ro_{l} ≈ 2/Co. As noted by Gastine et al. (2014), the transition is abrupt and occurs in a narrow parameter range. In Runs D and E, with the highest Coriolis numbers, the rotation profile is clearly solarlike. However, there are some interesting features in Run E: a strong polar jet appears on the northern hemisphere, and a decrease in is seen near the equator. Such polar vortices have frequently been found in similar simulations, see, e.g., Käpylä et al. (2010b, 2011b,a), and were also found in rapidly rotating convection in relatively thin shells (e.g. Elliott et al. 2000; Gastine & Wicht 2012). In the current setup this tendency is weaker, but we still occasionally observe polar jets (Runs E, B9, and B10; see also the left panel of Fig. 9.
We note that Run A is expected to be closest to the Sun with the highest convective energy flux. Furthermore, our choice of Pr_{SGS} = 0.25 leads to a fairly large contribution of SGSflux within the convection zone, which in the Sun is transported by convection. These results imply that the convective velocities in the Sun would be even higher, leading to lower Co and conditions that are more suitable for antisolar differential rotation. There are, however, hints that the velocities in simulations might be significantly higher than those in the Sun (cf. Hanasoge et al. 2012; Miesch et al. 2012).
3.3. Meridional circulation
Figure 5 shows the meridional circulation from representative Runs A and D. In the relatively slowly rotating antisolar Run A, the flow is concentrated in a single anticlockwise cell mostly outside the tangent cylinder with a peak amplitude of 90 m s^{1}, which is clearly higher than what is observed in the Sun (e.g. Zhao & Kosovichev 2004). In Run D the circulation pattern extends to higher latitudes and consists of several cells at low latitudes. The cell at high latitudes is likely an artefact of the closed θboundary. A similar transition from multiple to single cells has been observed before in different settings (e.g. Käpylä et al. 2011a; Matt et al. 2011; Gastine et al. 2013). The flow amplitude near the surface in Run D is of the order of 30 m s^{1}, which is still somewhat higher than the 20 m s^{1} obtained from helioseismology (Zhao & Kosovichev 2004).
A singlecell poleward circulation with solarlike rotation has been reported from simulations in spherical shells with the ASH code by imposing a latitudinal entropy variation on the bottom boundary (Miesch 2007; Miesch et al. 2011). In our sphericalwedge simulations, such a circulation pattern in combination with a solarlike differential rotation profile has so far occurred only as a transitory phenomenon in runs that have not yet fully relaxed, and they typically end up in the antisolar regime. Recent helioseismic studies suggest that the solar meridional circulation pattern consists of several cells in radius and possibly also in latitude (Zhao et al. 2013; Schad et al. 2013; Kholikov et al. 2014), which is also realized in our more rapidly rotating cases but is at odds with meanfield models of solar rotation (e.g. Rempel 2005; Kitchatinov & Rüdiger 2005).
3.4. Flow bistability
We confirm recent results of Gastine et al. (2014) that near the transition from solarlike to antisolar differential rotation, two stable solutions for the largescale flow exist for the same parameter values, only depending on the initial conditions.
Our results for and are shown in Fig. 6 for three sets of models (cf. Table 1). Firstly, we run models from the initial conditions described in Sect. 2.1. Furthermore, we run two additional sets where we take a snapshot from an antisolar and a solarlike solution as initial conditions. In the last two sets we vary the Rayleigh and Prandtl numbers by changing δn and hence the radiative conductivity K(r), while keeping the other control parameters at fixed values. We find that for these choices of parameters and initial conditions, it is more difficult to switch from antisolar to solarlike differential rotation than vice versa. This is seen by comparing the δn required for solarlike solutions in the different sets of runs: in Set D where we approach from the rapid rotation regime, the switch occurs between 0.68 < Ro_{c}< 0.70 (1.27 < Co < 1.37). In the opposite case of Set B, where we approach from the antisolar branch, the switch occurs between 0.53 < Ro_{c}< 0.55 (2.43 < Co < 2.01). In the case of Runs A–E that were run from scratch, we found 0.63 < Ro_{c}< 0.65 (1.33 < Co < 1.50). Thus, in terms of the Coriolis number the bistability region extends farther into the antisolar regime than the solarlike one. Physically, this might be related to the fact that in this case the strength of the differential rotation is much larger (see the two panels of Fig. 6). We have considered a single value of the Taylor number in our study. We note that according to Gastine et al. (2014), the size of the bistable region is wider with higher Ta.
Fig. 7 Timeaveraged rotation profiles from Runs C and D1. 
Fig. 8 Radial and latitudinal differential rotation as functions of time from Runs D1 and B3 with the same parameters. 
Figure 7 shows the rotation profiles from Runs C and D1 that have the same control parameters but different histories: Run C was run from the initial conditions described in Sect. 2.1, whereas in Run D1 we used the final thermally relaxed state of Run D (=D0) as initial condition. The resulting evolution of and is shown in Fig. 8, where we also show the corresponding results for Run B3 with the same input parameters as Run D1 after restarting from Run B (=B0). Our simulations were typically run for roughly 100 years solar time. By comparison, the viscous and SGS diffusion times, τ_{ν} = (Δr)^{2}/ν and τ_{SGS} = τ_{ν}Pr_{SGS}, in our simulations are 10.5 and 2.6 years, respectively.
Rapidly rotating toroidal jets at high latitudes appear in many numerical simulations of global scale convection (e.g. Miesch et al. 2000; Käpylä et al. 2011a,b). In the Sun the latitudinal gradient of is known to be monotonic in each hemisphere. Models based on solar differential rotation have also been adopted in stellar studies where low latitude bands or polar vortices are ignored. Nonmagnetic meanfield models also tend to produce this type of solution. However, the gas giants in the solar system have alternating bands of slower and faster rotation, although in that case it is not clear whether the convective layer is deep (e.g. Busse 1976; Heimpel & Aurnou 2007) or shallow (Kaspi et al. 2013). There are also theoretical and observational studies suggesting similar profiles for fully convective Mdwarfs and brown dwarfs (e.g. Balbus & Weiss 2010; Crossfield et al. 2014).
Fig. 9 Timeaveraged rotation profiles from Runs B10 and B10b. 
In Fig. 9 we show the rotation profiles obtained in Runs B10 and B10b with the same control parameters. Run B10 was run from a snapshot of Run B0 exhibiting antisolar differential rotation. This run now exhibits polar jets. However, it has been suggested that these jets can be unstable (Wicht et al. 2002). To investigate this possibility further, we take a snapshot from B10 as initial for Run B10b and apply a relaxation term at high latitudes for the azimuthally averaged u_{φ} that yields a monotonic latitudinal gradient of . The relaxation timescale is roughly four days and the term is switched on for two weeks in solar time in the beginning of the simulation. We find that the resulting profile with a more solarlike monotonic behaviour is also stable at least for 50 years, which corresponds to roughly five viscous diffusion times.
3.5. Λ effect and turbulent viscosity
The changes in the differential rotation should somehow be reflected in similar changes in the underlying mechanism responsible for driving it, which is the Λ effect (Rüdiger 1980, 1989). The Λ effect corresponds to a rank three tensor that parameterizes the nondiffusive contributions to the Reynolds stress, in addition to the diffusive contributions that result from turbulent viscosity. The Reynolds stress is given by , where is the fluctuating velocity. The relevant offdiagonal components can be written as where Λ_{V} and Λ_{H} are the vertical and horizontal components of the Λ effect and ν_{t} is the turbulent viscosity. Obviously, we cannot extract both effects selfconsistently from a single Reynolds stress component. Instead, we use a simple mixing length formula to estimate the turbulent viscosity (24)where u_{rms} = u_{rms}(r,θ) varies across the meridional plane, α_{MLT} = 1.7 is taken for the mixing length parameter, and H_{p} = −(∂lnp/∂r)^{1} is the pressure scale height. Similar methods have been applied in earlier studies of isotropically forced turbulence (Snellman et al. 2009) and convection (Pulkkinen et al. 1993; Käpylä et al. 2010a) in Cartesian geometry, where anisotropy is selfconsistently produced by rotation and/or stratification. In Fig. 10 we present the results for Runs A and D, which are representative of antisolar and solarlike rotation regimes.
Fig. 10 From left to right: timeaveraged Reynolds stresses Q_{rφ} and Q_{θφ} normalized by ν_{t}Ω_{⊙}, the turbulent viscosity divided by the molecular viscosity ν_{t}/ν, Λ_{V} and Λ_{H} normalized by ν_{t}, and the anisotropy parameters A_{V} and A_{H}. Top row: Run A; bottom row: Run D1. In the fifth column we only use data some degrees away from the equator so as to avoid the singularity associated with the division by cosθ. The contours in the lower row are oversaturated near the θboundaries in order to highlight the features at lower latitudes. 
In both cases, the Reynolds stress responsible for radial transport of angular momentum, Q_{rφ}, is negative at high latitudes. Somewhat surprisingly, Q_{rφ} is very small near the equator in the relatively slowly rotating Run A, but consistent with early results of Rieutord et al. (1994). In Run D1, in the solarlike regime, a positive contribution to the stress appears at low latitudes. The latter is consistent with Cartesian simulations where the corresponding stress changes sign at low latitudes in the rapid rotation regime (Käpylä et al. 2004). The horizontal stress, Q_{θφ}, leads to angular momentum transport that is directed toward the equator in both runs, except at high latitudes in Run A where the transport is toward the poles.
The mixing length estimate of the turbulent viscosity shows a decrease at high latitudes in comparison to the equatorial regions in Run A. In Run D1 there is a minimum at midlatitudes. In both cases the maximum value of ν_{t}/ν is of the order of 35, which is reasonable given that the Reynolds numbers in both runs are similar. We have here neglected anisotropies that are caused by gravity and rotation that play the roles of preferred directions in the system. Thus, we do not expect the details of the turbulent viscosity to be captured accurately. However, the order of magnitude of ν_{t}/ν ≈ Re seems reasonable, so we proceed to using this estimate to extract the Λ effect.
Solving Eqs. (22) and (23) for Λ_{V} and Λ_{H} yields The profiles of Λ_{V} show some distinct differences between both runs. It is mostly negative for the antisolar Run A and mostly positive for the solarlike Run D1. Interestingly, these differences would not have been so obvious if we had directly compared the vertical stresses Q_{rφ} of both runs. This highlights the usefulness of employing Eq. (25), even though this involves the uncertainty of estimating ν_{t}. Likewise, while the horizontal stresses Q_{θφ} for Runs A and D1 show some similarities, the profiles of Λ_{H} also show differences between both runs, and it is mostly negative for the antisolar Runs A and mostly positive for the solarlike Run D1.
According to meanfield theory (Rüdiger 1980), coefficients Λ_{V} and Λ_{H} are related to the anisotropy parameters (27)via Λ_{V} ≈ 2τA_{V} and Λ_{H} ≈ 2τA_{H}, where τ is the correlation time of the turbulence. Profiles of A_{V} and A_{H} are shown in the last two columns of Fig. 10. In Run A, both Λ_{V} and A_{V} are mostly negative, while both are mostly positive for Run D1, in broad agreement with the radial differential rotation gradient. Also Λ_{H} and A_{H} are mostly negative in Run A and mostly positive in Run D1, again in broad agreement with the latitudinal differential rotation.
According to firstorder smoothing results (e.g. Kitchatinov & Rüdiger 1995, 2005), A_{V} is always negative due to the simplified turbulence model used. This, however, is not always the case in simulations (e.g. Käpylä et al. 2004). Indeed, A_{V} is mostly negative in Run A and attains positive values near the equator in Run D1, in qualitative accordance to our estimates of Λ_{V} for each case. Similarly we find that the sign of A_{H} agrees mostly with that of Λ_{H}. It is remarkable that the results for the Λ effect are consistent with those of the anisotropy parameters given our use of a very simple approximation for ν_{t}. Our results for A_{V} and A_{H} also indicate that the differential rotation has a significant impact on the properties of turbulence.
For the run with solarlike rotation (Run D), Λ_{H} is concentrated near the equator and close to the upper boundary. A similar concentration near the equator has also been observed in earlier studies (e.g. Chan 2001; Käpylä et al. 2004) and can be partly explained by the banana cells near the equator (Käpylä et al. 2011b).
4. Conclusions
Simulations of mildly turbulent threedimensional convection in spherical wedges have allowed us to study aspects of differential rotation relevant to the Sun. In our most solarlike model, Run A, we found antisolar differential rotation. Decreasing the convective velocities by increasing the radiative conductivity, thereby also increasing the Coriolis number, we found a threshold after which the rotation changes to a solarlike profile. In agreement with the recent finding of Gastine et al. (2014) using Boussinesq convection, our stratified and compressible models near the threshold between both states confirm the existence of bistability in that the question of antisolar or solarlike rotation depends on the initial condition or the history of the run. We also confirmed earlier findings showing that highlatitude toroidal jets are possible in runs where solarlike solutions are observed with the same parameters. This could be significant for stellar rotation profiles that are naively assumed to have a monotonic latitudinal gradient of angular velocity. These jets may well be asymmetric and present only at one of the poles (see also Heimpel & Aurnou 2007; Jones & Kuzanyan 2009).
Another interesting but speculative conjecture is that the solar differential rotation is also the result of bistability. Our Runs A and B with energy fluxes closest to what we would expect from the Sun i.e. where resolved convection transports most of the total flux, show antisolar differential rotation. Making the model more realistic with respect to the Sun requires that we (i) decrease the SGSflux, which now contributes more than ten per cent of the total; (ii) increase the density stratification by a factor of five; and (iii) increase the Reynolds number. All of these factors are likely to lead to higher convective velocities and thus lower Coriolis number, in turn leading to even more suitable conditions for antisolar differential rotation. We know, however, that the Sun rotated much more rapidly in its youth, which favours a solarlike rotation profile. Furthermore, the results of Gastine et al. (2014) suggest that the range of parameters in which bistable solutions are possible increases as the Taylor number increases. These two facts lend some credence to the conjecture that the Sun is in a bistable regime. Clearly, further research into this matter is required.
Acknowledgments
We thank the referee for a thorough report and constructive criticism. The simulations were performed using the supercomputers hosted by the CSC – IT Center for Science Ltd. in Espoo, Finland, which is administered by the Finnish Ministry of Education. Financial support from the Academy of Finland grants No. 136189, 140970, 272786 (PJK) and 272157 to the ReSoLVE Centre of Excellence (MJM), the University of Helsinki research project “Active Suns”, as well as the Swedish Research Council grants 62120115076 and 20125797 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.
References
 Aurnou, J., Heimpel, M., & Wicht, J. 2007, Icarus, 190, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Weiss, N. O. 2010, MNRAS, 404, 1263 [NASA ADS] [Google Scholar]
 Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Moss, D., & Tuominen, I. 1992, A&A, 265, 328 [NASA ADS] [Google Scholar]
 Brandenburg, A., Chan, K. L., Nordlund, Å., & Stein, R. F. 2005, Astron. Nachr., 326, 681 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., & Palacios, A. 2009, ApJ, 702, 1078 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 1976, Icarus, 29, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, K. L. 2001, ApJ, 548, 1102 [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]
 Cole, E., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2014, ApJ, 780, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Collier Cameron, A. 2002, Astron. Nachr., 323, 336 [Google Scholar]
 Crossfield, I. J. M., Biller, B., Schlieder, J. E., et al. 2014, Nature, 505, 654 [NASA ADS] [CrossRef] [Google Scholar]
 Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336 [NASA ADS] [CrossRef] [Google Scholar]
 Elliott, J. R., Miesch, M. S., & Toomre, J. 2000, ApJ, 533, 546 [NASA ADS] [CrossRef] [Google Scholar]
 Gastine, T., & Wicht, J. 2012, Icarus, 219, 428 [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]
 Gilman, P. A. 1977, Geophys. Astrophys. Fluid Dynam., 8, 93 [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Nat. Acad. Sci., 109, 11928 [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. 2014, ApJ, 786, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, C. A., & Kuzanyan, K. M. 2009, Icarus, 204, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J. 2011, Astron. Nachr., 332, 43 [NASA ADS] [CrossRef] [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., 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]
 Kaspi, Y., Showman, A. P., Hubbard, W. B., Aharonson, O., & Helled, R. 2013, Nature, 497, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Kholikov, S., Serebryanskiy, A., & Jackiewicz, J. 2014, ApJ, 784, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R. 1963, ApJ, 137, 664 [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]
 Krause, F., & Rüdiger, G. 1974, Astron. Nachr., 295, 93 [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. 2007, Astron. Nachr., 328, 998 [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, Annu. Rev. Fluid Mech., 41, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., Brown, B. P., Browning, M. K., Brun, A. S., & Toomre, J. 2011, in Astrophysical Dynamics: From Stars to Galaxies, eds. N. H. Brummell, A. S. Brun, M. S. Miesch, & Y. Ponty, IAU Symp., 271, 261 [Google Scholar]
 Miesch, M. S., Featherstone, N. A., Rempel, M., & Trampedach, R. 2012, ApJ, 757, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D. 2009, ApJ, 697, 923 [NASA ADS] [CrossRef] [Google Scholar]
 Ossendrijver, M. 2003, A&ARv, 11, 287 [Google Scholar]
 Pulkkinen, P., Tuominen, I., Brandenburg, A., Nordlund, A., & Stein, R. F. 1993, A&A, 267, 265 [NASA ADS] [Google Scholar]
 Rempel, M. 2005, ApJ, 622, 1320 [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 Dynam., 16, 239 [Google Scholar]
 Rüdiger, G. 1989, Differential Rotation and Stellar Convection. Sun and Solartype Stars (Berlin: Akademie Verlag) [Google Scholar]
 Schad, A., Timmer, J., & Roth, M. 2013, ApJ, 778, L38 [NASA ADS] [CrossRef] [Google Scholar]
 Snellman, J. E., Käpylä, P. J., Korpi, M. J., & Liljeström, A. J. 2009, A&A, 505, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stix, M. 2002, The Sun: An Introduction (Berlin: Springer) [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]
 Vitense, E. 1953, Z. Astrophys., 32, 135 [NASA ADS] [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]
 Wicht, J., Jones, C. A., & Zhang, K. 2002, Icarus, 155, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, R. K., Gastine, T., Christensen, U. R., & Duarte, L. D. V. 2013, ApJ, 774, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, J., & Kosovichev, A. G. 2004, ApJ, 603, 776 [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 Left and middle panels: temperature T in units of 10^{6} K and density ρ in units of kg m^{3}, respectively, for the thermally relaxed state (black solid line), initial condition (red dotted), and hydrostatic solutions (blue dashed). Rightmost panel: Prandtl numbers related to radiative diffusion Pr = ν/χ where χ = K/ρc_{P} (black solid line), and the turbulent heat conductivity Pr_{SGS} = ν/χ_{SGS} (blue dashed) as functions of radius. Data is taken from Run D. 

In the text 
Fig. 2 Timeaveraged energy fluxes from Runs A (top) and E (bottom): radiative (thin solid line), convective (dashed), kinetic energy (dashdotted), SGS (tripledashdotted), and viscous (longdashed) flux. The thick solid line denotes the total flux, whereas the red horizontal dotted lines show the zero and unity line. The red vertical dotted line at r = r_{m} shows the midpoint of the convection zone. 

In the text 
Fig. 3 Radial dependence of the time averaged fluctuating rmsvelocity where the contributions from the mean flows are omitted for Runs A (solid black line), C (blue dashed), and E (red dotdashed). The black dotted line shows the convective velocity from the mixing length (ML) model of Stix (2002). 

In the text 
Fig. 4 Timeaveraged rotation profiles from Runs A–E showing in nHz. 

In the text 
Fig. 5 Timeaveraged meridional circulation from Runs A (left) and D (right). The arrows show the flow , whereas the colour contours show . 

In the text 
Fig. 6 Radial (left panel) and latitudinal (right panel) differential rotation, defined by Eqs. (21), from Runs A–E (diamonds), and Set D (blue dotted line with asterisks) and B (red dashed line with triangles). 

In the text 
Fig. 7 Timeaveraged rotation profiles from Runs C and D1. 

In the text 
Fig. 8 Radial and latitudinal differential rotation as functions of time from Runs D1 and B3 with the same parameters. 

In the text 
Fig. 9 Timeaveraged rotation profiles from Runs B10 and B10b. 

In the text 
Fig. 10 From left to right: timeaveraged Reynolds stresses Q_{rφ} and Q_{θφ} normalized by ν_{t}Ω_{⊙}, the turbulent viscosity divided by the molecular viscosity ν_{t}/ν, Λ_{V} and Λ_{H} normalized by ν_{t}, and the anisotropy parameters A_{V} and A_{H}. Top row: Run A; bottom row: Run D1. In the fifth column we only use data some degrees away from the equator so as to avoid the singularity associated with the division by cosθ. The contours in the lower row are oversaturated near the θboundaries in order to highlight the features at lower latitudes. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.