EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 570 (October 2014) A&A, 570 (2014) A43 Full HTML
Free Access
Issue
A&A
Volume 570, October 2014
Article Number A43
Number of page(s) 10
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201423412
Published online 14 October 2014

© 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 near-surface 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 non-diffusive contribution to the Reynolds stress known as Λ effect in mean-field 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 Taylor-Proudman balance is broken in the Sun. Mean-field 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 near-surface shear layer (e.g. Kitchatinov & Olemskoy 2011).

Mean-field models rarely produce anti-solar 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 anti-solar 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 anti-solar 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 near-surface shear layer of the Sun (Kitchatinov & Rüdiger 2005). Differential rotation and meridional circulation are intimately coupled (Kippenhahn 1963), and one implies the other. Anti-solar differential rotation implies a counter-clockwise circulation in the northern hemisphere, corresponding to poleward motion at the surface. In the context of the solar near-surface shear layer, this mechanism is referred to as gyroscopic pumping (Miesch & Hindman 2011).

Both mean-field models and simulations have limitations. Furthermore, approximations like first-order smoothing are used to derive turbulent transport coefficients, such as turbulent viscosity and the so-called Λ effect. However, their range of validity in stellar convection zones is questionable and only limited comparisons between hydrodynamic mean-field models and direct simulations have been performed (see, however Rieutord et al. 1994). Direct numerical simulations of convection, on the other hand, lack small-scale 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 anti-solar differential rotation in some slowly rotating stars (e.g. Strassmeier et al. 2003; Weber et al. 2005) and solar-like rotation in rapidly rotating dwarfs (e.g. Collier Cameron 2002). However, for the star with the best observational evidence for anti-solar 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 anti-solar to solar-like 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 large-scale 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 anti-solar to solar-like 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 solar-like 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 solar-like rotation profiles into the anti-solar 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 anti-solar transition in terms of Coriolis number. Secondly, we investigate the dependence of the solutions on initial conditions by running models from solar and anti-solar 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 convection-driven 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 r0 ≤ r ≤ r1, θ0 ≤ θ ≤ π − θ0, and 0 ≤ φ ≤ φ0, respectively, where r0 = 0.72 R⊙ and r1 = 0.97 R⊙ denote the positions of the bottom and top of the computational domain, and R⊙ = 7 × 108 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 subgrid-scale (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 γ = cP/cV = 5/3 is the ratio of specific heats at constant pressure and volume, respectively, and e = cVT 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/r3, where G = 6.67 × 10-11 m3 kg-1 s-2 is the gravitational constant, and M⊙ = 2.0 × 1030 kg is the mass of the Sun. We neglect self-gravity 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 106 times higher luminosity in the model in comparison to the Sun. As the convective energy flux scales as Fconv ~ ρu3, 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 Lratio = L0/L⊙, with L0 and L⊙ ≈ 3.84 × 1026 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 T0 = 2.23×106 K.

The initial state is isentropic and the hydrostatic temperature gradient is given by (7)where nad = 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) = K0 [n(r) + 1], with n(r) = δn(r/r0)-15 + nad − δn replacing the polytropic index, (8)being a reference conductivity, and ℒ being the non-dimensional luminosity, given below. Now n = nad at the bottom of the convection zone and approaches nad − δn at the surface. This means that K = (n + 1)K0 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 PrSGS = ν/χSGS and Pr = ν/χ, where χ = K/ρcP, are shown in Fig. 1. We introduce weak small-scale Gaussian noise velocity perturbations in the initial state.

Our simulations are defined by the energy flux imposed at the bottom boundary, Fb = −(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 × 108 m2 s-1 in physical units at r = r1.

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 Stefan-Boltzmann 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, σT4, carries the total luminosity through the boundary in the initial non-convecting state.

thumbnail Fig. 1

Left and middle panels: temperature T in units of 106 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/ρcP (black solid line), and the turbulent heat conductivity PrSGS = ν/χSGS (blue dashed) as functions of radius. Data is taken from Run D.

2.2. Dimensionless parameters

The non-dimensional input parameters of our models are the luminosity parameter (12)the normalized pressure scale height at the surface, (13)with T1 being the temperature at the surface, the Taylor number (14)where Δr = r1 − r0 = 0.25 R⊙, as well as the fluid and SGS Prandtl numbers (15)where χm = K(rm) /cPρm is the thermal diffusivity and ρm is the density, both evaluated at r = rm = 0.845 R⊙. We vary Pr and keep PrSGS = 0.25 fixed. Finally, the non-dimensional 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 urms 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 = Co2Re2(kfR⊙)4. Furthermore, we define the Rayleigh number as (19)where shs is the entropy in the hydrostatic, non-convecting state. We compute the hydrostatic stratification by evolving a one-dimensional 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 Code1, which uses a high-order finite difference method for solving the compressible equations of magnetohydrodynamics.

thumbnail Fig. 2

Time-averaged energy fluxes from Runs A (top) and E (bottom): radiative (thin solid line), convective (dashed), kinetic energy (dash-dotted), SGS (triple-dash-dotted), and viscous (long-dashed) 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 = rm shows the midpoint of the convection zone.

thumbnail Fig. 3

Radial dependence of the time averaged fluctuating rms-velocity where the contributions from the mean flows are omitted for Runs A (solid black line), C (blue dashed), and E (red dot-dashed). 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 anti-solar differential rotation, while Runs D and E have solar-like differential rotation. Secondly, we study the bistability of the rotation profile by taking either a solar-like (Runs D0–D4) or an anti-solar (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.

Table 1

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 rms-velocities 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 urms. 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 urms near the lower boundary is due to our impenetrable boundary condition. Near the upper boundary our values of urms 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 high-resolution simulations of non-rotating 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.

thumbnail Fig. 4

Time-averaged rotation profiles from Runs A–E showing in nHz.

thumbnail Fig. 5

Time-averaged 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 rms-velocity, 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 “solar-like”.

thumbnail 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 anti-solar 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 solar-like. Thus, the transition from the anti-solar to solar-like 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 Rol ≈ 1 where Rol ≈ 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 solar-like. 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 PrSGS = 0.25 leads to a fairly large contribution of SGS-flux 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 anti-solar 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 anti-solar Run A, the flow is concentrated in a single anti-clockwise 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 single-cell poleward circulation with solar-like 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 spherical-wedge simulations, such a circulation pattern in combination with a solar-like 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 anti-solar 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 mean-field 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 solar-like to anti-solar differential rotation, two stable solutions for the large-scale 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 anti-solar and a solar-like 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 anti-solar to solar-like differential rotation than vice versa. This is seen by comparing the δn required for solar-like solutions in the different sets of runs: in Set D where we approach from the rapid rotation regime, the switch occurs between 0.68 < Roc< 0.70 (1.27 < Co < 1.37). In the opposite case of Set B, where we approach from the anti-solar branch, the switch occurs between 0.53 < Roc< 0.55 (2.43 < Co < 2.01). In the case of Runs A–E that were run from scratch, we found 0.63 < Roc< 0.65 (1.33 < Co < 1.50). Thus, in terms of the Coriolis number the bistability region extends farther into the anti-solar regime than the solar-like 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.

thumbnail Fig. 7

Time-averaged rotation profiles from Runs C and D1.

thumbnail 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 = τνPrSGS, 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. Non-magnetic mean-field 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 M-dwarfs and brown dwarfs (e.g. Balbus & Weiss 2010; Crossfield et al. 2014).

thumbnail Fig. 9

Time-averaged 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 anti-solar 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 solar-like 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 non-diffusive 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 off-diagonal 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 self-consistently from a single Reynolds stress component. Instead, we use a simple mixing length formula to estimate the turbulent viscosity (24)where urms = urms(r,θ) varies across the meridional plane, αMLT = 1.7 is taken for the mixing length parameter, and Hp = −(∂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 self-consistently produced by rotation and/or stratification. In Fig. 10 we present the results for Runs A and D, which are representative of anti-solar and solar-like rotation regimes.

thumbnail Fig. 10

From left to right: time-averaged Reynolds stresses Qrφ and Qθφ normalized by νtΩ⊙, the turbulent viscosity divided by the molecular viscosity νt/ν, ΛV and ΛH normalized by νt, and the anisotropy parameters AV and AH. 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, Qrφ, is negative at high latitudes. Somewhat surprisingly, Qrφ 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 solar-like 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 mid-latitudes. 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 anti-solar Run A and mostly positive for the solar-like Run D1. Interestingly, these differences would not have been so obvious if we had directly compared the vertical stresses Qrφ 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 anti-solar Runs A and mostly positive for the solar-like Run D1.

According to mean-field theory (Rüdiger 1980), coefficients ΛV and ΛH are related to the anisotropy parameters (27)via ΛV ≈ 2τAV and ΛH ≈ 2τAH, where τ is the correlation time of the turbulence. Profiles of AV and AH are shown in the last two columns of Fig. 10. In Run A, both ΛV and AV are mostly negative, while both are mostly positive for Run D1, in broad agreement with the radial differential rotation gradient. Also ΛH and AH are mostly negative in Run A and mostly positive in Run D1, again in broad agreement with the latitudinal differential rotation.

According to first-order smoothing results (e.g. Kitchatinov & Rüdiger 1995, 2005), AV 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, AV 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 AH 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 AV and AH also indicate that the differential rotation has a significant impact on the properties of turbulence.

For the run with solar-like 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 three-dimensional convection in spherical wedges have allowed us to study aspects of differential rotation relevant to the Sun. In our most solar-like model, Run A, we found anti-solar 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 solar-like 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 anti-solar or solar-like rotation depends on the initial condition or the history of the run. We also confirmed earlier findings showing that high-latitude toroidal jets are possible in runs where solar-like 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 anti-solar differential rotation. Making the model more realistic with respect to the Sun requires that we (i) decrease the SGS-flux, 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 anti-solar differential rotation. We know, however, that the Sun rotated much more rapidly in its youth, which favours a solar-like 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 621-2011-5076 and 2012-5797 and the European Research Council under the AstroDyn Research Project 227952, are acknowledged, as well as the HPC-Europa2 project, funded by the European Commission – DG Research in the Seventh Framework Programme under grant agreement No. 228398.

References

  1. Aurnou, J., Heimpel, M., & Wicht, J. 2007, Icarus, 190, 110 [NASA ADS] [CrossRef] [Google Scholar]
  2. Balbus, S. A., & Weiss, N. O. 2010, MNRAS, 404, 1263 [NASA ADS] [Google Scholar]
  3. Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 [NASA ADS] [CrossRef] [Google Scholar]
  4. Brandenburg, A., Moss, D., & Tuominen, I. 1992, A&A, 265, 328 [NASA ADS] [Google Scholar]
  5. Brandenburg, A., Chan, K. L., Nordlund, Å., & Stein, R. F. 2005, Astron. Nachr., 326, 681 [NASA ADS] [CrossRef] [Google Scholar]
  6. Brun, A. S., & Palacios, A. 2009, ApJ, 702, 1078 [NASA ADS] [CrossRef] [Google Scholar]
  7. Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
  8. Busse, F. H. 1976, Icarus, 29, 255 [NASA ADS] [CrossRef] [Google Scholar]
  9. Chan, K. L. 2001, ApJ, 548, 1102 [NASA ADS] [CrossRef] [Google Scholar]
  10. Chan, K. L. 2010, in IAU Symp. 264, eds. A. G. Kosovichev, A. H. Andrei, & J.-P. Rozelot, 219 [Google Scholar]
  11. Cole, E., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2014, ApJ, 780, L22 [NASA ADS] [CrossRef] [Google Scholar]
  12. Collier Cameron, A. 2002, Astron. Nachr., 323, 336 [Google Scholar]
  13. Crossfield, I. J. M., Biller, B., Schlieder, J. E., et al. 2014, Nature, 505, 654 [NASA ADS] [CrossRef] [Google Scholar]
  14. Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336 [NASA ADS] [CrossRef] [Google Scholar]
  15. Elliott, J. R., Miesch, M. S., & Toomre, J. 2000, ApJ, 533, 546 [NASA ADS] [CrossRef] [Google Scholar]
  16. Gastine, T., & Wicht, J. 2012, Icarus, 219, 428 [NASA ADS] [CrossRef] [Google Scholar]
  17. Gastine, T., Wicht, J., & Aurnou, J. M. 2013, Icarus, 225, 156 [NASA ADS] [CrossRef] [Google Scholar]
  18. Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gilman, P. A. 1977, Geophys. Astrophys. Fluid Dynam., 8, 93 [Google Scholar]
  20. Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
  21. Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Nat. Acad. Sci., 109, 11928 [NASA ADS] [CrossRef] [Google Scholar]
  22. Heimpel, M., & Aurnou, J. 2007, Icarus, 187, 540 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hotta, H., & Yokoyama, T. 2011, ApJ, 740, 12 [NASA ADS] [CrossRef] [Google Scholar]
  24. Hotta, H., Rempel, M., & Yokoyama, T. 2014, ApJ, 786, 24 [NASA ADS] [CrossRef] [Google Scholar]
  25. Jones, C. A., & Kuzanyan, K. M. 2009, Icarus, 204, 227 [NASA ADS] [CrossRef] [Google Scholar]
  26. Käpylä, P. J. 2011, Astron. Nachr., 332, 43 [NASA ADS] [CrossRef] [Google Scholar]
  27. Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2004, A&A, 422, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Käpylä, P. J., Brandenburg, A., Korpi, M. J., Snellman, J. E., & Narayan, R. 2010a, ApJ, 719, 67 [NASA ADS] [CrossRef] [Google Scholar]
  29. Käpylä, P. J., Korpi, M. J., Brandenburg, A., Mitra, D., & Tavakol, R. 2010b, Astron. Nachr., 331, 73 [NASA ADS] [CrossRef] [Google Scholar]
  30. Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011a, Astron. Nachr., 332, 883 [NASA ADS] [CrossRef] [Google Scholar]
  31. 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]
  32. Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJ, 755, L22 [NASA ADS] [CrossRef] [Google Scholar]
  33. Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [NASA ADS] [CrossRef] [Google Scholar]
  34. Kaspi, Y., Showman, A. P., Hubbard, W. B., Aharonson, O., & Helled, R. 2013, Nature, 497, 344 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kholikov, S., Serebryanskiy, A., & Jackiewicz, J. 2014, ApJ, 784, 145 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kippenhahn, R. 1963, ApJ, 137, 664 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kitchatinov, L. L., & Olemskoy, S. V. 2011, MNRAS, 411, 1059 [NASA ADS] [CrossRef] [Google Scholar]
  38. Kitchatinov, L. L., & Rüdiger, G. 1995, A&A, 299, 446 [NASA ADS] [Google Scholar]
  39. Kitchatinov, L. L., & Rüdiger, G. 2004, Astron. Nachr., 325, 496 [NASA ADS] [CrossRef] [Google Scholar]
  40. Kitchatinov, L. L., & Rüdiger, G. 2005, Astron. Nachr., 326, 379 [NASA ADS] [CrossRef] [Google Scholar]
  41. Krause, F., & Rüdiger, G. 1974, Astron. Nachr., 295, 93 [NASA ADS] [CrossRef] [Google Scholar]
  42. Matt, S. P., Do Cao, O., Brown, B. P., & Brun, A. S. 2011, Astron. Nachr., 332, 897 [Google Scholar]
  43. Miesch, M. S. 2007, Astron. Nachr., 328, 998 [NASA ADS] [CrossRef] [Google Scholar]
  44. Miesch, M. S., & Hindman, B. W. 2011, ApJ, 743, 79 [NASA ADS] [CrossRef] [Google Scholar]
  45. Miesch, M. S., & Toomre, J. 2009, Annu. Rev. Fluid Mech., 41, 317 [NASA ADS] [CrossRef] [Google Scholar]
  46. Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [NASA ADS] [CrossRef] [Google Scholar]
  47. 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]
  48. Miesch, M. S., Featherstone, N. A., Rempel, M., & Trampedach, R. 2012, ApJ, 757, 128 [NASA ADS] [CrossRef] [Google Scholar]
  49. Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D. 2009, ApJ, 697, 923 [NASA ADS] [CrossRef] [Google Scholar]
  50. Ossendrijver, M. 2003, A&ARv, 11, 287 [Google Scholar]
  51. Pulkkinen, P., Tuominen, I., Brandenburg, A., Nordlund, A., & Stein, R. F. 1993, A&A, 267, 265 [NASA ADS] [Google Scholar]
  52. Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
  53. Rieutord, M., Brandenburg, A., Mangeney, A., & Drossart, P. 1994, A&A, 286, 471 [NASA ADS] [Google Scholar]
  54. Rüdiger, G. 1980, Geophys. Astrophys. Fluid Dynam., 16, 239 [Google Scholar]
  55. Rüdiger, G. 1989, Differential Rotation and Stellar Convection. Sun and Solar-type Stars (Berlin: Akademie Verlag) [Google Scholar]
  56. Schad, A., Timmer, J., & Roth, M. 2013, ApJ, 778, L38 [NASA ADS] [CrossRef] [Google Scholar]
  57. 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]
  58. Stix, M. 2002, The Sun: An Introduction (Berlin: Springer) [Google Scholar]
  59. Strassmeier, K. G., Kratzwald, L., & Weber, M. 2003, A&A, 408, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Thompson, M. J., Christensen-Dalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599 [NASA ADS] [CrossRef] [Google Scholar]
  61. Vitense, E. 1953, Z. Astrophys., 32, 135 [NASA ADS] [Google Scholar]
  62. Warnecke, J., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2013, ApJ, 778, 141 [NASA ADS] [CrossRef] [Google Scholar]
  63. Weber, M., Strassmeier, K. G., & Washuettl, A. 2005, Astron. Nachr., 326, 287 [NASA ADS] [CrossRef] [Google Scholar]
  64. Wicht, J., Jones, C. A., & Zhang, K. 2002, Icarus, 155, 425 [NASA ADS] [CrossRef] [Google Scholar]
  65. Yadav, R. K., Gastine, T., Christensen, U. R., & Duarte, L. D. V. 2013, ApJ, 774, 6 [NASA ADS] [CrossRef] [Google Scholar]
  66. Zhao, J., & Kosovichev, A. G. 2004, ApJ, 603, 776 [NASA ADS] [CrossRef] [Google Scholar]
  67. 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

Table 1

Summary of the runs.

In the text

All Figures

thumbnail Fig. 1

Left and middle panels: temperature T in units of 106 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/ρcP (black solid line), and the turbulent heat conductivity PrSGS = ν/χSGS (blue dashed) as functions of radius. Data is taken from Run D.
In the text
thumbnail Fig. 2

Time-averaged energy fluxes from Runs A (top) and E (bottom): radiative (thin solid line), convective (dashed), kinetic energy (dash-dotted), SGS (triple-dash-dotted), and viscous (long-dashed) 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 = rm shows the midpoint of the convection zone.
In the text
thumbnail Fig. 3

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

Time-averaged rotation profiles from Runs A–E showing in nHz.
In the text
thumbnail Fig. 5

Time-averaged meridional circulation from Runs A (left) and D (right). The arrows show the flow , whereas the colour contours show .
In the text
thumbnail 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
thumbnail Fig. 7

Time-averaged rotation profiles from Runs C and D1.
In the text
thumbnail Fig. 8

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

Time-averaged rotation profiles from Runs B10 and B10b.
In the text
thumbnail Fig. 10

From left to right: time-averaged Reynolds stresses Qrφ and Qθφ normalized by νtΩ⊙, the turbulent viscosity divided by the molecular viscosity νt/ν, ΛV and ΛH normalized by νt, and the anisotropy parameters AV and AH. 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 (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

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

Initial download of the metrics may take a while.