© ESO, 2011

1. Introduction

Many stars show signs of differential rotation. For instance, the solar equator rotates with a shorter period than the solar polar caps. The difference of 132 nHz between the rotation rates found by Ulrich et al. (1988) corresponds to a difference of 0.07 rad/day between the respective angular velocities or a lapping time of 88 days. Helioseismology has found that this pattern persists throughout the whole solar convection zone but not in the radiative zone below (Thompson et al. 2003). Stellar differential rotation can be inferred from the light curves of rotating spotted stars (see, e.g. Henry et al. 1995; Messina & Guinan 2003), by monitoring magnetic activity (Donahue et al. 1996), spectroscopically (Reiners & Schmitt 2003), or by Doppler imaging (Barnes et al. 2000; Donati et al. 2000).

While several studies have found a systematic dependence of the surface differential rotation on the rotation rate, no such dependence has been found when the samples are combined (Hall 1991; Barnes et al. 2005). Moreover, measurements of the surface differential rotation of the rapidly rotating K dwarfs PZ Tel and AB Dor with the Doppler imaging technique show that stars rotating much more rapid than the Sun have surface shear values very similar to solar values – as first predicted by Kitchatinov & Rüdiger (1999). PZ Tel is a young K dwarf with a rotation period of 0.95 days, whose surface differential rotation δΩ = 0.075 rad/day is remarkably close to that of the Sun (Barnes et al. 2000). AB Dor is a rapidly rotating K0 dwarf with a rotation period of 0.51 days, whose surface differential rotation was found to vary with time between 0.09 and 0.05 rad/day (Collier Cameron & Donati 2002).

Barnes et al. (2005) proposed a dependence on the effective temperature (albeit with large scatter) and found the power law $\mathit{\delta }\mathrm{\Omega }\mathrm{=}{\mathrm{\Omega }}_{\mathrm{eq}}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{pole}}\mathrm{\propto }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{8.92}\mathrm{\pm }\mathrm{0.31}}\mathit{,}$(1)where δΩ is the difference between the angular velocities at the equator and the polar caps¹. The power law in Eq. (1) was confirmed by Reiners (2006) using spectroscopic methods. The light curves of the stars ϵ Eri and κ¹ Cet recorded by the MOST satellite (Croll et al. 2006; Walker et al. 2007) both indicate by measuring δΩ ≃ 0.062 and δΩ ≃ 0.064, respectively, a very similar equator-pole difference of the rotation rate as the Sun. The value of 0.11 for the young G star CoRoT-2a (Fröhlich et al. 2009) also agrees with the common picture of a rotation-independent surface shear for G stars.

Equation (1) disagrees, however, with recent observations of the young G dwarfs LQ Lup, R58 and HD 171488. Marsden et al. (2005) report a surface differential rotation of 0.025  ±  0.015 rad/day for R58, while Jeffers & Donati (2009b) find the much larger value of 0.138  ±  0.011 rad/day. Donati et al. (2000) find a surface differential rotation δΩ = 0.12  ±  0.02 rad/day for LQ Lup. Jeffers & Donati (2009a) determined the surface rotation of HD 171488 using the Zeeman-Doppler imaging technique and found a large surface differential rotation of 0.5 rad/day. Marsden et al. (2006) report a smaller but still large value of 0.402  ±  0.044 rad/day. Huber et al. (2009) found no evidence at all of differential rotation but could not exclude it either. Jeffers et al. (2011) confirmed the findings of Marsden et al. (2006).

The values that have been found for LQ Lup and R58 are in line with the Barnes et al. picture but the large values found for HD 171488 are not. The three G dwarfs are similar in their ages, effective temperatures, and radii, yet HD 171488 has a much larger differential rotation than the other two stars. The studies mentioned have focused on rotation rate and effective temperature as the properties determining the surface differential rotation. As the stars are very similar in terms of their stellar parameters such as age and effective temperature, we ask whether there could be a difference in their internal structure that would cause the observed difference in their surface rotation. All three stars have just reached the zero-age main-sequence or are approaching it. Given the rapid retreat of the convection zone in the final part of the pre-main sequence evolution and the uncertainty in the stellar age, we ask whether the strong differential rotation of HD 1714888 is caused by the different depth of its outer convection zone compared to LQ Lup and R58.

In the following, we compute model convection zones of different depth and their large-scale gas motions, i.e. rotation and meridional flow. The models are based on the mean field formulation of fluid dynamics that has been very successful in describing the Sun, reproducing its surface rotation, internal rotation in the convection zone, and its surface meridional flow very well (Kitchatinov & Rüdiger 1999, KR99; Küker & Stix 2001, KS01). Models have been constructed for a variety of stars with spectral types from M to F (Küker & Rüdiger 2008).

The new scheme that we present here allows the computation of stellar rotation laws and meridional flow patterns based on a mean-field model of the large-scale flows in stellar convection zones for high rotation rates when narrow boundary layers exist. It assumes strict spherical symmetry for the basic stratification, ignoring any flattening that might occur for very rapid rotation. However, the impact of rotation on the thermal structure can be taken into account by including a gravity darkening term in the heat transport equation so that the model remains applicable even for moderately flattened stars.

2. Theory of stellar rotation laws

While thermal convection is driven by the stratification of the star, the convective gas motions carry angular momentum as well as heat. Momentum transport by turbulent motions is known as Reynolds stress and can be described as a turbulent viscosity in cases of the simplest shear flows. In a rotating, stratified convection zone, the Reynolds stress is anisotropic and its azimuthal components have a contribution proportional to the angular velocity, Ω, itself rather than its gradient. This form of Reynolds stress (with “Λ effect”) is incompatible with solid body rotation.

2.1. Basic equations

Our model consists of a 1D background model, which assumes hydrostatic equilibrium and spherical symmetry and a system of partial differential equations that describe the convective heat flux, the transport of angular momentum, and the meridional flow, assuming axisymmetry. The equation for the convective heat transport is then given by $\mathrm{\nabla }\mathrm{·}\mathrm{(}{F}\mathrm{conv}\mathrm{+}{F}\mathrm{rad}\mathrm{)}\mathrm{-}\mathit{\rho T}{u̅}\mathrm{·}\mathrm{\nabla }\mathit{S}\mathrm{=}\mathrm{0}\mathit{,}$(2)where F^conv and F^rad are the convective and radiative heat flux, respectively, ρ the mass density,  the mean gas velocity, T the gas temperature, and S is the specific entropy $\mathit{S}\mathrm{=}{\mathit{C}}_{\mathrm{v}}\log \left(\frac{\mathit{p}}{{\mathit{\rho }}^{\mathit{\gamma }}}\right)\mathit{,}$(3)where p is the gas pressure, C_v the specific heat at constant volume, and γ the adiabatic index. In spherical polar coordinates, the Reynolds equation can be rewritten as a system of two partial differential equations. The azimuthal component of the Reynolds equation expresses the conservation of angular momentum $\mathrm{\nabla }\mathrm{·}{t}\mathrm{=}\mathrm{0}\mathit{,}$(4)where t is the angular momentum flux vector ${\mathit{t}}_{\mathit{i}}\mathrm{=}\mathit{r}\sin \mathit{\theta }\mathrm{(}\mathit{\rho r}\sin \mathit{\theta }\mathrm{\Omega }{\mathit{u}}_{\mathit{i}}^{\mathrm{m}}\mathrm{+}\mathit{\rho }{\mathit{Q}}_{\mathit{i\phi }}\mathrm{)}\mathit{,}$(5)which has two transporters of angular momentum: i) the meridional flow (u^m) and ii) the zonal flux of the turbulent angular momentum (Q_iφ).

The equation for the meridional flow is derived from the Reynolds equation by taking the azimuthal component of its curl ${\left[\mathrm{\nabla }\mathrm{\times }\frac{\mathrm{1}}{\mathit{\rho }}\mathrm{\nabla }\mathrm{·}\mathit{R}\right]}_{\mathit{\phi }}\mathrm{+}\mathit{r}\sin \mathit{\theta }\frac{\mathit{\partial }{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{\partial z}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\rho }}^{\mathrm{2}}}\mathrm{(}\mathrm{\nabla }\mathit{\rho }\mathrm{\times }\mathrm{\nabla }\mathit{p}{\mathrm{)}}_{\mathit{\phi }}\mathrm{+}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{=}\mathrm{0}\mathit{,}$(6)where R_ij =  −ρQ_ij is the Reynolds stress and ∂/∂z  =  cosθ∂/∂r − sinθ∂/∂θ is the derivative along the axis of rotation. For very rapid rotation, the second term on the left hand side (LHS) of Eq. (6) dominates and the rotation rate is constant along cylindrical surfaces, as stated by the Taylor-Proudman theorem.

2.2. Transport coefficients

If spherical polar coordinates are chosen, the Λ effect appears in two components of the Reynolds stress $\begin{array}{ccc}{\mathit{Q}}_{\mathit{r\phi }}& \mathrm{=}& {\mathit{Q}}_{\mathit{r\phi }}^{\mathrm{visc}}\mathrm{+}\mathit{V}{\mathit{\nu }}_{\mathrm{t}}\sin \mathit{\theta }\mathrm{\Omega }\mathit{,}\\ {\mathit{Q}}_{\mathit{\theta \phi }}& \mathrm{=}& {\mathit{Q}}_{\mathit{\theta \phi }}^{\mathrm{visc}}\mathrm{+}\mathit{H}{\mathit{\nu }}_{\mathrm{t}}\cos \mathit{\theta }\mathrm{\Omega }\mathit{,}\end{array}$\arraycolsep1.75ptwhere ${\mathit{Q}}_{\mathit{r\phi }}^{\mathrm{visc}}$ and ${\mathit{Q}}_{\mathit{\theta \phi }}^{\mathrm{visc}}$ contain only first order derivatives of Ω with respect to r and θ and therefore vanish for uniform rotation. The coefficients V and H refer to the vertical and horizontal parts of the Λ effect, while ν_t represents the eddy viscosity in the convection zone. As the corresponding terms contain the angular velocity Ω itself, they do not vanish for rigid rotation. Thus, rigid rotation is not stress-free.

For slow rotation, the convective heat flux is proportional to the gradient of the specific entropy ${\mathit{F}}_{\mathit{i}}^{\mathrm{conv}}\mathrm{=}\mathrm{-}\mathit{\rho T}{\mathit{\chi }}_{\mathrm{t}}\frac{\mathit{\partial S}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathit{,}$(9)where χ_t the turbulent heat conductivity. The turbulent heat conductivity is determined by the convection velocity u_c and the convective turnover time τ_c${\mathit{\chi }}_{\mathrm{t}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\tau }}_{\mathrm{c}}{\mathit{u}}_{\mathrm{c}}^{\mathrm{2}}\mathit{.}$(10)Using standard mixing length theory, we find that ${\mathit{u}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{{\mathit{l}}_{\mathrm{m}}^{\mathrm{2}}\mathit{g}}{\mathrm{4}{\mathit{C}}_{\mathrm{p}}}\frac{\mathit{\partial S}}{\mathit{\partial r}}$(11)\newpage\noindentand ${\mathit{\tau }}_{\mathrm{c}}\mathrm{=}\frac{{\mathit{l}}_{\mathrm{m}}}{{\mathit{u}}_{\mathrm{c}}}\mathit{,}$(12)where l_m is the mixing length. Equation (9) describes a strictly radial, diffusive heat flux. In a rotating convection zone, the heat flux is not strictly aligned with the entropy gradient, i.e. ${\mathit{F}}_{\mathit{i}}^{\mathrm{conv}}\mathrm{=}\mathrm{-}\mathit{\rho T}{\mathit{\chi }}_{\mathrm{t}}{\mathrm{\Phi }}_{\mathit{ij}}\frac{\mathit{\partial S}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}$(13)(Rüdiger 1989; Kitchatinov et al. 1994), where the dimensionless coefficients Φ_ij are functions of the Coriolis number Ω ^∗  = 2τ_cΩ fulfilling  | Φ_ij |  ≤ 1. The flux vector is then tilted towards the rotation axis.

The radiative heat flux is given by ${\mathit{F}}_{\mathit{i}}^{\mathrm{rad}}\mathrm{=}\mathrm{-}\frac{\mathrm{16}\mathit{\sigma }{\mathit{T}}^{\mathrm{3}}{\mathrm{\nabla }}_{\mathit{i}}\mathit{T}}{\mathrm{3}\mathit{\kappa \rho }}\mathit{,}$(14)where σ is the Stefan-Boltzmann constant and κ the opacity.

2.3. Background model

We assume that the star is essentially in hydrostatic equilibrium and that all gas motions constitute small perturbations of that equilibrium that do not change the structure of the star. We particularly assume that no net gas transport occurs $\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{u̅}\mathrm{\right)}\mathrm{=}\mathrm{0.}$(15)Our model convection zone is a spherical layer with adiabatic stratification that is heated from below and cooled at the surface. We assume an ideal gas with the equations of state $\mathit{p}\mathrm{=}\frac{\mathrm{ℛ}}{\mathit{\mu }}\mathit{\rho T}\mathrm{=}{\mathit{C}}_{\mathrm{p}}\frac{\mathit{\gamma }\mathrm{-}\mathrm{1}}{\mathit{\gamma }}\mathit{\rho T,}$(16)where p is the gas pressure, ℛ the gas constant, and C_p the heat capacity for constant pressure, and hydrostatic equilibrium $\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\mathit{g\rho }\mathit{.}$(17)From the law of gravity $\mathit{g}\mathrm{=}\mathrm{-}\frac{\mathit{GM}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathit{,}$(18)where G is the gravitational constant and M(r) the mass contained in the sphere of radius r, we derive $\frac{\mathrm{d}\mathit{g}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\mathrm{2}\frac{\mathit{g}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{r}}\mathrm{+}\mathrm{4}\mathit{\pi G\rho }\mathit{.}$(19)For adiabatic stratification where pρ ^− γ = const., the density can be expressed in terms of the temperature $\mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{e}}{\left(\frac{\mathit{T}}{{\mathit{T}}_{\mathrm{e}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{,}$(20)where ρ_e and T_e are the values of density and temperature at a reference radius r_e. The condition of hydrostatic equilibrium can then be rewritten in terms of the temperature, i.e. $\frac{\mathrm{d}\mathit{T}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathit{g}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{C}}_{\mathrm{p}}}\mathrm{·}$(21)Equations (19) − (21) form a system that can be integrated from the reference radius, where the values of density, gravity, and temperature are known from the full stellar evolution model.

2.4. Boundary conditions and numerical method

As boundary conditions, we assume stress-free and impenetrable boundaries for the gas motion and an imposed radial heat flux. At the lower boundary, the heat flux is constant corresponding to the stellar luminosity for the heat transport ${\mathit{F}}_{\mathit{r}}\mathrm{=}\frac{\mathit{L}}{\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{b}}^{\mathrm{2}}}\mathit{,}$(22)where L is the stellar luminosity and r_b the radius at which the boundary is located. Since we assume a particular form of radiative heat flux and define the bottom of the convection zone as the point where the radiative heat flux equals the total heat flux, this condition is equivalent to imposing zero convective heat flux.

At the upper boundary, we assume that the convective heat flux equals the radiative flux and is converted into radiation (Gilman & Glatzmaier 1981), according to the linearized condition ${\mathit{F}}_{\mathit{r}}\mathrm{=}\frac{\mathit{L}}{\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{t}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{4}\frac{\mathit{S}}{{\mathit{C}}_{\mathrm{p}}}\mathrm{)}\mathit{,}$(23)where r_t is the radius of that boundary. The boundary conditions on the rotation axis apply because of the axisymmetry.

To solve the above system, we introduce the stream function ψ for the meridional flow ${\mathit{u}}_{\mathit{r}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\rho }{\mathit{r}}^{\mathrm{2}}\sin \mathit{\theta }}\frac{\mathit{\partial \psi }}{\mathit{\partial \theta }}\mathit{,}{\mathit{u}}_{\mathit{\theta }}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho r}\sin \mathit{\theta }}\frac{\mathit{\partial \psi }}{\mathit{\partial r}}\mathit{,}$(24)and expand the specific entropy, the angular velocity, and the stream function in terms of spherical harmonics $\begin{array}{ccc}& & \mathit{S}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{S}}_{\mathit{n}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\overline{\mathit{P}}}_{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}\mathit{,}\\ & & \mathrm{\Omega }\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{\Omega }}_{\mathit{n}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{\overline{\mathit{P}}\begin{array}{}\mathrm{1}\\ \mathrm{2}\mathit{n}\mathrm{-}\mathrm{1}\end{array}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}}{\sin \mathit{\theta }}\mathit{,}\\ & & \end{array}$(25)In the above equation, ${\overline{\mathit{P}}}_{\mathit{l}}$ and $\overline{\mathit{P}}\begin{array}{}\mathrm{1}\\ \mathit{l}\end{array}$ are the normalized standard and adjoint Legendre polynomials. The expansions are convenient because the functions of colatitude θ on the right-hand side are the eigenfunctions of the angular parts of the corresponding diffusion operators. In using only Legendre polynomials of only even or odd degree in a certain expansion, we assume symmetry with respect to the equator. The expansion shown in Eq. (25) usually converges rapidly such that N = 5 is sufficient for cases such as the solar rotation. More rapidly rotating stars require values up to 20.

Applying the expansions in Eq. (25) to the system of partial differential equations in Eqs. (2), (4), and (6) transforms the equations into a system of ordinary nonlinear differential equations with the independent variable r. We convert the equations into the system of first-order differential equations by introducing new dependent variables. In particular, the Reynolds stress components R_rθ and R_rφ are adopted as new dependent variables. The resulting system of first-order differential equations is solved by the standard relaxation method described by Press et al. (1992).

The rapidly rotating stars have very thin boundary layers near the top and bottom of their convection zones. To resolve the layers, a nonuniform numerical grid (zeros of Chebyshev polynomials) is applied $\begin{array}{ccc}{\mathit{r}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left({\mathit{r}}_{\mathrm{t}}\mathrm{+}{\mathit{r}}_{\mathrm{b}}\mathrm{-}\mathrm{(}{\mathit{r}}_{\mathrm{t}}\mathrm{-}{\mathit{r}}_{\mathrm{b}}\mathrm{)}\cos \left(\mathit{\pi }\frac{\mathit{i}\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{n}\mathrm{-}\mathrm{2}}\right)\right)\mathit{,}& & \\ \mathrm{2}\mathrm{\le }\mathit{i}\mathrm{\le }\mathit{n}\mathrm{-}\mathrm{1}\mathit{,}{\mathit{r}}_{\mathrm{1}}\mathrm{=}{\mathit{r}}_{\mathrm{b}}\mathit{,}{\mathit{r}}_{\mathit{n}}\mathrm{=}{\mathit{r}}_{\mathrm{t}}\mathit{,}& & \end{array}$(26)where n is the number of radial grid points. The grid has fine spacing near the boundaries. To solve for the solar rotation law, a value of n as small as 20 is sufficient, but we usually use higher resolutions (e.g. n = 200), especially for rapid rotation.

3. The Sun

We first simulate the solar rotation law and compare the results to those from earlier models by Kitchatinov & Rüdiger (1999) and Küker & Stix (2001) by solving the same set of equations with different numerical methods. KR99 used a simple Kramers’ law for the opacity, while KS01 used the opacity from Ahrens et al. (1992). We use a subroutine from Hansen & Kawaler (1994) that is based on the analytical opacity law by Stellingwerf (1975). A more elaborate treatment is possible but not likely to improve the model as we do not treat the uppermost layers of the convection zone and the atmosphere. As the solar convection zone is not fully ionized up to the photosphere, the adiabatic gradient, ∇_ad = (∂log T/∂log P)_S, and hence the specific heat capacity are not constant. While ∇_ad = 0.4 holds throughout most of the convection zone, it drops to values as low as 0.1 in the upper 35 000 km. As we consider only constant values, we either have to exclude the uppermost layer or adjust ∇_ad to reproduce the depth of the convection zone from the stellar evolution model.

Fig. 1 The solar differential rotation and meridional flow as computed with the new scheme. Top left: surface rotation rate versus colatitude. Top right: rotation rate as function of radius at 0, 15, 30, 45, 60, 75, and 90° latitude, respectively, from top to bottom. Bottom left: streamlines of the meridional flow. The flow is counter-clockwise, i.e. directed towards the pole at the surface and towards the equator at the bottom. Bottom right: flow speed as a function of latitude at the top (solid blue) and bottom (dashed red) of the convection zone. Positive values indicate a flow away from the north pole.

Fig. 2 Temperature deviation as function of the latitude at the bottom (left) and top (right) of the solar convection zone. The polar axis is always warmer than the equator. The differences, however, are small.

The top panels in Fig. 1 show the rotation pattern in the solar convection zone. It agrees well with those computed with the KR99 codes and the KS01 scheme with the modifications described in Bonanno et al. (2007). The surface rotation is most rapid at the equator with a difference $\mathit{\delta }\mathrm{\Omega }\mathrm{=}\mathrm{0.07}\mathrm{rad}\mathit{/}\mathrm{day}$(27)between the rotation rates of the equator and the polar caps, corresponding to a relative shear, δΩ/Ω_eq, of 29 percent. The variation with radius at mid-latitudes is weak. The angular velocity decreases with radius at high latitudes, while at the equator it slightly increases with radius (“superrotation”) in the deep convection zone. In the surface layer, there is a negative gradient of Ω (“subrotation”) at all latitudes. As for the rotation profiles computed with the KR99 and KS01 schemes, this rotation profile is in excellent agreement with the findings of helioseismology.

The bottom panels of Fig. 1 show the meridional flow. There is one cell per hemisphere with the surface flow directed towards the poles and the return flow at the bottom of the convection zone. The amplitude of this “counter-clockwise” flow is 14 m/s at the (model) surface and 6 m/s at the bottom. The difference between the flow speeds at the top and bottom, respectively, is smaller than might be expected from the density stratification and the requirement of mass conservation but the concentration of the return flow to a thin layer allows a relatively fast flow despite the larger mass density.

Figure 2 shows the quantity $\mathit{\delta T}\mathrm{=}\mathrm{(}\mathit{S}\mathrm{-}{\mathit{S}}_{\mathrm{0}}\mathrm{)}\frac{\mathit{T}}{{\mathit{C}}_{\mathrm{p}}}$(28)at the top and the bottom of the solar convection zone, where S₀ is the specific entropy at the bottom of the convection zone at the equator. This choice of S₀ implies negative values of δT at the top of the convection zone. In general, the polar axis is warmer than the more equatorward parts of the fluid. Unfortunately, the small size of the temperature difference does not allow any empirical confirmation. Though the deviation from the adiabatic stratification is only small, the resulting baroclinic term has a profound impact on the large-scale meridional flow and the differential rotation (see below).

3.1. Rapidly rotating Sun

To study the impact of the basic rotation we compute the rotation law of a hypothetical rapidly rotating Sun with the P = 1.33 day rotation period of HD 171488. Figure 3 shows the resulting rotation pattern. The isocontours are cylindrically shaped according to the Taylor-Proudman theorem. The radial profiles in the right-hand diagram show an increase in the rotation rate with increasing radius at low latitudes. Unlike the case of the real Sun, there is a pronounced increase in the rotation rate in the upper part of the convection zone at the equator. At both radial boundaries, there are pronounced boundary layers with huge gradients in the rotation rate caused by the stress-free boundary conditions. The surface rotation is most rapid at the equator and decreases monotonically towards the polar caps but the slope of the decline has a minimum at mid-latitudes. With a value of 0.08 rad/day instead of 0.07 rad/day, the horizontal shear is only slightly stronger than for the Sun despite there being a factor of 20 between the average rotation rates.

The meridional flow has only one flow cell per hemisphere with the surface flow directed towards the poles. The amplitudes are 34 m/s at the surface and 17 m/s at the bottom of the convection zone. As for the real Sun with its much slower rotation, the return flow is half as fast as the surface flow but even more concentrated at the bottom of the convection zone. The flow is fastest at mid-latitudes.

Fig. 3 Rotation of a hypothetical rapidly rotating Sun with P = 1.33 d.

3.2. MOST stars: ϵ Eri and κ¹ Cet

ϵ Eri and κ¹ Cet are active dwarf stars for which surface differential rotation has been derived from light curves recorded by the MOST satellite. Croll et al. (2006) derived a rotation law of the form $\mathit{P}\mathrm{\left(}\mathit{B}\mathrm{\right)}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{eq}}}{\mathrm{1}\mathrm{-}\mathit{k}{\sin }^{\mathrm{2}}\mathit{B}}$(29)for the K2 V star ϵ Eri, where P(B) is the rotation period at the latitude B, and is P_eq the rotation period at the equator (B = 0). The parameter k gives the shear of the rotation law. The general rule derived from the theory of the stellar rotation laws of Kitchatinov & Rüdiger (1999) $\mathit{k}\mathrm{\simeq }\frac{{\mathit{P}}_{\mathrm{eq}}}{\mathrm{100} \mathrm{days}}$(30)is perfectly fulfilled by these stars.

An equatorial rotation period P_eq = 11.2 days and a value $\mathit{k}\mathrm{=}\mathrm{0.1}{\mathrm{1}}_{-0.02}^{\mathrm{+}\mathrm{0.03}}$ for the differential rotation has been reported for ϵ Eri corresponding to δΩ = 0.062 rad/day. For the G5 V star κ¹ Cet, Walker et al. (2007) found P_eq = 8.77 days and $\mathit{k}\mathrm{=}\mathrm{0.0}{\mathrm{9}}_{-0.005}^{\mathrm{+}\mathrm{0.006}}$ corresponding to δΩ = 0.064 rad/day.

Kitchatinov & Olemskoy (2011) computed rotation laws for model stars with 0.8 and 1.0 solar masses for ϵ Eri and κ¹ Cet. The results were k = 0.127 for ϵ Eri and k = 0.13 for κ¹ Cet.

Here the calculations for these stars are based on improved stellar models computed with the MESA/STARS code (Paxton et al. 2011). A star with M _⋆  = 0.85   M_⊙, Z = 0.02 and an age of 0.44 Gyr serves as a model for ϵ Eri. It has an effective temperature of 5076 K, a radius of 0.76 R_⊙, and a luminosity of 0.34 L_⊙. The bottom of its outer convection zone is located at 0.69 R _⋆ . For this model star, we find a surface differential rotation δΩ = 0.51 rad/day or k = 0.10. This is in agreement with the observed value.

For κ¹ Cet, we assume a solar mass and metallicity and an age of 1 Gyr. The model star has an effective temperature of 5677 K, radius 0.915 R_⊙, and luminosity 0.78 L_⊙. The resulting rotation law shows a surface shear δΩ = 0.077 rad/day, or k = 0.12. This is slightly larger than the observed value. An earlier model from the same evolutionary track with a stellar age of 167 Myr, an effective temperature of 5649 K, a radius of 0.9 R_⊙, and an luminosity of 0.74 L_⊙ yields k = 0.11 or δΩ = 0.072 rad/day. This is still not in perfect agreement with the observed value but reasonably close.

Possible reasons for the remaining discrepancy are the stellar model, which might not be a sufficiently accurate representation of the real star, and an underestimate of the total surface shear by the sin²θ derived from spot rotation as the observed spots do not cover the whole range of latitudes from pole to equator.

4. Numerical experiments

We carry out several numerical experiments to compare the roles played by the two drivers of differential rotation, i.e. the Λ effect and the baroclinic flow.

4.1. Sun without Λ effect

As mentioned, our model includes two effects that are capable to maintain differential rotation, the Λ effect and the baroclinic flow that is caused by a horizontal temperature gradient. The latter effect can be seen when the baroclinic term in Eq. (6) is rewritten in terms of the temperature $\frac{\mathrm{1}}{{\mathit{\rho }}^{\mathrm{2}}}\mathrm{(}\mathrm{\nabla }\mathit{\rho }\mathrm{\times }\mathrm{\nabla }\mathit{p}{\mathrm{)}}_{\mathit{\phi }}\mathrm{\approx }\mathrm{-}\frac{\mathit{g}}{\mathit{rT}}\frac{\mathit{\partial \delta T}}{\mathit{\partial \theta }}\mathrm{·}$(31)The baroclinic term can be a powerful driver of meridional flows, which in turn can drive differential rotation. In our theory, the latitudinal temperature profile is due to the anisotropic heat-conductivity tensor in the presence of rotation. The polar axis is always slightly warmer than the equator but with an unobservable temperature excess. A positive pole-equator temperature difference will drive a clockwise flow. Its angular momentum transport leads to accelerated rotation at low latitudes and slower rotation at the polar caps. A baroclinic flow (also “thermal wind”) can therefore maintain differential rotation with a solar-type surface rotation. As an illustration, we repeat our computation for the Sun where the Λ terms have been canceled within the Reynolds stress.

The resulting rotation and flow patterns are shown in Fig. 4. The rotation is indeed solar-type but where δΩ = 0.04 is weaker than with the Λ effect and the isocontours are distinctly disk-shaped at the poles even in the mid-latitudes. The equator region shows basically no structure.

The typical superrotation of the deep convection zone beneath the equator known as a result of the helioseismology does not occur. The reason is simple: if only a (fast) meridional circulation transports the angular momentum with a pattern symmetric with respect to the equator then the angular momentum becomes uniform in the equatorial part of the convection zone (except, of course, the boundary layers). Hence, there is r²Ω ≃ const. independent of the flow direction but in contradiction to observations².

With amplitudes of 4.7 m/s at the top and 2.1 m/s at the bottom of the convection zone, the barocline-induced clockwise meridional flow is substantially weaker than in the complete model – and it has the “wrong” direction.

Fig. 4 Solar rotation law and meridional flow without Λ effect. Top: the rotation law. Note the steep negative gradient of the angular velocity at the poles and the almost rigid rotation beneath the equator. Bottom: the meridional flow is anti-solar: it is positive (equatorward) at the top of the convection zone and negative (poleward) at the bottom of the convection zone.

4.2. Sun without baroclinic flow

The barocline term is then neglected while the Λ effect is considered, which maintains the differential rotation together with the meridional flow caused by the former through the second term on the LHS of Eq. (6) (“Biermann-Kippenhahn flow”, see Kippenhahn 1963; Köhler 1970).

Figure 5 shows the resulting rotation pattern. The differential rotation is solar-type, i.e. the rotation period is shortest at the equator and longest at the polar caps. With δΩ = 0.014 rad/day, the surface rotation is much more rigid than observed. The isocontour plot shows a distinctly cylindrically shaped pattern in the bulk of the convection zone, while at the top and bottom boundaries the pattern deviates from the cylinder geometry and the rotation rate falls off with increasing radius at all latitudes. The counterclockwise meridional flow has a very similar geometry as in the full model but is faster with amplitudes of 17 and 8 m/s at the top and bottom boundaries.

The subrotation along the polar axis, which is typical for the solar rotation law and produced by baroclinic flows (cf. Fig. 4), does not exist here. Köhler (1970) and Rüdiger et al. (1998) have given several numerical examples of this direct consequence of the Taylor-Proudman theorem.

Fig. 5 Solar rotation law without baroclinic terms. The structure disappears now at the polar axis but appears at the equatorial region. The equatorial plane shows superrotation but the polar axis rotates almost rigidly.

5. Young G dwarfs

5.1. CoRoT-2a

This G7 V star is a young Sun, i.e. it has solar mass but is much younger with an age of 0.5 Gyr. The star was observed by the CoRoT satellite and found to have a planet with an orbital period of 1.743 days (Alonso et al. 2008). In addition to planetary transits, the light curve displays a periodic variation that is most easily explained by a rotating, spotted surface with basic rotation period of 4.5 days. Spot modeling using circular spots provides an excellent fit assuming three spots and a solar-type surface differential rotation with a difference of 0.11 rad/day (Fröhlich et al. 2009).

Fig. 6 Rotation pattern of CoRoT-2a, a young Sun rotating with P = 4.5 days.

Our model star is based on a model from an evolutionary track for the Sun. The age is 0.5 Gyr, the radius 0.9 solar radii, and the luminosity 75 percent of the solar value. The bottom of the convection zone is at a fractional radius x = 0.73.

Figure 6 shows the resulting rotation pattern. The surface differential rotation of 0.09 rad/day reproduces the value observed by Fröhlich et al. (2009). The rotation pattern is more cylindrical than that of the Sun, but not as much as that of our rapidly rotating Sun. Similarly, the boundary layers are more pronounced than in the Sun but less pronounced than for the rapidly rotating Sun. The flow pattern is similar to that in the solar convection zone. The amplitude is slightly larger with a maximum value of 18.6 m/s at the top and 10.6 m/s at the bottom of the convection zone.

5.2. R58, LQ Lup and HD 171488

We now address the differential rotation measurements of the rapidly rotating G dwarfs R58, LQ Lup, and HD 171488. R58 (HD 307938) is a young active G dwarf in the open cluster IC 2602. It has a photospheric temperature of 5800 K and rotates with a period of 0.56 days (Marsden et al. 2005). LQ Lup (RX J1508.6-4423) is a post-T Tauri star with an effective temperature of 5750  ±  50 K, a rotation period of 0.31 days, and a radius of 1.22  ±  0.12 solar radii (Donati et al. 2000), corresponding to a mass of 1.16  ±  0.04 solar masses and an age of 25  ±  10 Myr. The equator-pole Ω-difference is 0.12 rad/day. HD 171488 (V889 Her) is a young active G dwarf. Strassmeier et al. (2003) found a photospheric temperature of 5830 K, a rotation period of 1.337 days, a mass of 1.06  ±  0.05 solar masses, a radius of 1.09  ±  0.05 solar radii, and an age of 30−50 Myr. Marsden et al. (2006) find a rotation period of 1.313 days, a photospheric temperature of 5800 K, and a radius of 1.15  ±  0.08 solar radii.

With 0.4–0.5 rad/day the equator-pole Ω-difference of HD 171488 is exceptionally large. The k-value of 0.10 disagrees slightly with the k = 0.013 predicted by Eq. (30) by a factor of 7.5. As the corresponding factors for R58 and LQ Lup are much smaller (factor 2 − 3), the following calculations take the case of HD 171488 as the main example. As the stars have similar radii and effective temperatures, we investigate the impact differences in the internal structure on the surface rotation patterns. We study two models that differ mainly in terms of the metallicity and, as a consequence, the depth of the convection zone. The models were computed with the MESA/STARS code. We use the metallicity as a control parameter for the depth of the convection zone and vary the mass and the mixing length parameter to adjust the radius and effective temperature to values close to those observed for the three young G dwarfs. The rotation period is the 1.33 day period of HD 171488. Both model stars are 30 Myr old. We also compute one model with an artificially truncated convection zone and another with a very short rotation period of 1.33 days. The results are summarized in Table 1.

5.2.1. Deep convection zone

The first model star has the high metallicity of Z = 0.03. A mass of 1.11 solar masses and a value of 1.6 for the mixing length parameter lead to an effective temperature of 5685 K and a radius of 1.14 solar radii. The bottom of the convection zone is located at x = 0.765 and has a temperature of 1.46  ×  10⁶ K. The rotation law is shown in the top part of Fig. 7. The isocontour plot shows cylindrically shaped contours similar to those for the rapidly rotating Sun. The surface rotation law is solar-type with a rapidly rotating equator. The surface shear is much stronger. We find a value of 0.11 rad/day, which is about 1.5 times the solar value. The radial profiles show superrotation beneath the equator and subrotation along the rotation axis. The boundary layers are much less pronounced and the surface layer resembles the real Sun, i.e. the rotation rate decreases with increasing radius at all latitudes. The meridional flow is of solar-type (counterclockwise) with the surface flow towards the poles and the return flow located at the bottom of the convection zone. The amplitudes are 23 and 14 m/s, respectively.

Fig. 7 Internal rotation of young G dwarfs with rotation period of 1.33 days. Top: deep convection zone, δΩ = 0.11 rad/day at the surface. Bottom: shallow convection zone. δΩ = 0.5 rad/day at the surface.

5.2.2. Shallow convection zone

Our second model star has a metallicity of 0.01. In addition, a mass of 1.08 solar masses and a mixing length parameter of 1.0 correspond to a radius of 1.14 solar radii and an effective temperature of 5750 K. The bottom of the convection zone lies at a fractional radius of 0.89 and a temperature of 6  ×  10⁵ K. The bottom part of Fig. 7 shows the resulting rotation law. The surface rotation has a total shear of δΩ = 0.5 rad/day. The contour plot shows a disk-shaped pattern with flat profiles at low latitudes and a decrease in the rotation rate with increasing radius at the polar caps. The radial profiles show pronounced boundary layers at the radial boundaries. The flow pattern is mostly solar-type. There is one large flow cell with poleward flow at the top and the return flow at the bottom of the convection zone. In addition, there is a small cell of clockwise flow at low latitudes. The flow speeds at the top reach 18 m/s and 2.2 m/s, respectively. The return flow at the bottom of the convection zone reaches 3.5 m/s.

5.2.3. Truncated convection zone

The models above have (roughly) the same effective temperature but differ in terms of mass and structure. To isolate the effect of the depth of the convection zone, we compute a model of a young solar mass star and artificially reduce the depth of the convection zone by imposing the lower boundary at a fractional radius x = 0.88 instead of x = 0.78 and compare the rotation pattern with that of the full model. In both cases, the top boundary is at x = 0.98. The truncated convection zone thus has half the depth of the full model. In both cases, the rotation law is solar-type but while the full convection zone produces a surface shear of 0.12 rad/day, the truncated convection zone has 0.32 rad/day. The full model has a temperature difference of 24 K at the top and 124 K at the bottom. For the truncated model, the temperature difference between the polar cap and the equator is 47 K at the top of the convection zone and 254 K at its bottom.

5.2.4. Very rapid rotation

Fig. 8 Differential rotation of the G dwarf with the shallow convection zone and a rotation period of 0.33 days. Left: isocontours of the angular velocity. Right: the rotation profiles at various latitudes.

To illustrate the impact of the rotation rate, we repeat the computation for the shallow convection zone model with rotation period of 0.33 days. Figure 8 shows the resulting rotation pattern. The left plot shows solar-type rotation in both the surface differential rotation and the predominantly radial isocontours. This is still not cylindrically-shaped but closer to the Taylor-Proudman state than the disk-shaped pattern we found for a rotation period of 1.33 days. The surface shear is 0.38 rad/day and the temperature difference between polar caps and equator is 1785 K at the bottom and 184 K at the top. The maximum flow speeds are 14 m/s at the top and 8 m/s at the bottom of the convection zone.

6. Discussion

6.1. Shallow versus deep convection zone

Our “shallow” and “deep” convection zone models have similar radii and effective temperatures but the stars differ in terms of mass and metallicity. We have also chosen different values for the mixing length parameter. As a result, the depth of the convection zone differs between the models, and we find very different results for the differential rotation of these convection zones. The model with a shallow convection zone has much stronger differential rotation than the ones with deeper convection zones. Our “shallow convection zone” model does indeed reproduce the very large value of 0.5 rad/day found by Jeffers & Donati (2009a,b), while the “deep convection zone” model with its smaller value of 0.11 is in rough agreement with the values observed for LQ Lup and R58.

The result of the shallow convection zone model is in line with previous findings for F stars. Extreme values of surface shear were observed for F stars by Reiners (2006) and found in theoretical models for stars with very shallow convection zones by Küker & Rüdiger (2007). In the latter study, the largest values of surface shear were found for the most massive stars studied. These stars are not only the hottest, their convection zones are also extremely shallow.

To produce a thin convection zone, we had to lower the mixing length parameter to a value of 1.0. This is not only different from the value used for the “deep” convection zone model but much lower than usually used in stellar evolution models. Hence, the model is probably unrealistic and does not directly apply to LQ Lup. Not being specialists in the field, we leave the question of whether a shallow convection zone is possible open but we speculate whether strong magnetic fields could inhibit the convective heat transport in a way that might be mimicked by this choice of parameter.

The “truncated convection zone” model avoids the uncertainties in the “shallow convection zone” model and shows even more clearly the impact that the depth of the convection zone has on the surface rotation. Both the vertical and horizontal temperature gradients are steepest for the model with the shallow convection zone and flattest for the model with the deep convection zone. At the bottom of the convective zone, the differences are 750 K and 100 K for a convection zone of small or large depth, respectively. The corresponding values at the top are 55 K and 19 K.

6.2. Thermal wind equilibrium

We always find a higher temperature at high latitudes than at the equator. This is a consequence of the tilt of the convective heat transport vector towards the axis of rotation caused by the Coriolis force. In spherical polar coordinates, the components of the heat flux read $\begin{array}{ccc}& & {\mathit{F}}_{\mathit{r}}\mathrm{=}\mathrm{-}\mathit{\rho T}{\mathit{\chi }}_{\mathit{rr}}\frac{\mathit{\partial S}}{\mathit{\partial r}}\mathrm{-}\mathit{\rho T}{\mathit{\chi }}_{\mathit{r\theta }}\frac{\mathit{\partial S}}{\mathit{\partial \theta }}\mathit{,}\\ & & {\mathit{F}}_{\mathit{\theta }}\mathrm{=}\mathrm{-}\mathit{\rho T}{\mathit{\chi }}_{\mathit{\theta r}}\frac{\mathit{\partial S}}{\mathit{\partial r}}\mathrm{-}\mathit{\rho T}{\mathit{\chi }}_{\mathit{\theta \theta }}\frac{\mathit{\partial S}}{\mathit{\partial \theta }}\mathrm{·}\end{array}$\arraycolsep1.75ptThe first term in the horizontal component precludes a purely radial stratification. Any variation in the specific entropy with radius will cause a horizontal heat flux thus produce a horizontal gradient. This case is profoundly different from that of a latitude-dependent but still purely radial heat flux, which would have a much smaller impact and cause a much weaker differential rotation (Rüdiger et al. 2005).

The impact on the maintenance of differential rotation can be seen from the equation for the meridional flow. For rapid rotation (i.e. large Taylor number), Reynolds stress and nonlinear terms can be neglected and Eq. (6) is reduced to the thermal wind equation $\mathrm{2}\mathit{r}\sin \mathit{\theta }\mathrm{\Omega }\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial z}}\mathrm{-}\frac{\mathit{g}}{\mathit{rT}}\frac{\mathit{\partial \delta T}}{\mathit{\partial \theta }}\mathrm{=}\mathrm{0.}$(34)It follows that a gradient in the angular velocity along the axis of rotation is needed to balance the horizontal temperature gradient. Hence, a disk-shaped rotation pattern is generated in the polar area (only) because of the action of the baroclinic flow. For warmer poles (δT sinks equatorward) the z-gradient of Ω is negative as observed. The meridional flow without the baroclinic component only yields dΩ/dz ≃ 0 (Köhler 1970, see also Fig. 5). Hence, the empirical finding of negative dΩ/dz at the polar axis by helioseismology proves the existence of baroclinic flows in the solar convection zone.

As we have demonstrated by means of Fig. 4, the superrotation beneath the solar equator is a direct consequence of the Reynolds stress in the convection zone. All models without the Λ effect but with a clockwise (equatorward at the surface) circulation lead to dΩ/dr ≲ 0 in the bulk of the convection zone beneath, close to the equator. These clockwise flows are able to accelerate the equator relative to the poles but are unable to reproduce the superrotation beneath the equator. The same relation between baroclinic flows and differential rotation has been found for stellar radiative zones by Rieutord (2006) and Espinosa Lara & Rieutord (2007).

Figure 9 shows the results of a detailed computation of the terms on the LHS of Eq. (34) and their sum versus the fractional radius for a latitude of 45°. The left panel shows the plot for the Sun and the right panel that for the deep convection zone G dwarf model. The same quantities were computed for the rapid Sun and the shallow convection zone G dwarf models.

For the real Sun, the two terms cancel in the bulk of the solar convection zone where the meridional flow is slow and smooth but not close to the boundaries. For the deep convection zone, the boundary layers are much thinner but also much more pronounced. These findings indicate that for rapid rotation the meridional flow is mainly driven by the boundary layers, while the bulk of the convection zone is in thermal wind equilibrium. In the boundary layers, the meridional circulation is strongly sheared so that the right hand side of Eq. (34), which consists of third derivatives of the meridional velocity components multiplied by the eddy viscosity, becomes large.

Fig. 9 The drivers of angular momentum versus fractional radius at 45° latitude for the Sun (left, slow rotation) and the deep convection zone G dwarf (right, rapid rotation). Dashed line: baroclinic term. Dash-dotted line: centrifugal term. Solid line: both.

The numerical experiments show that the baroclinic term is an important but not the sole driver of the solar differential rotation. Neglecting the Λ effect leaves a surface shear of 0.04 rad/day, while switching off the baroclinic term reduces it to 0.014 rad/day. This is not because the Λ effect is an inefficient transporter of angular momentum but because it is balanced by the meridional flow. Disregarding the meridional flow altogether increases the surface shear to 0.18 rad/day. This is a well-known effect: meridional flows driven by differential rotation then act to reduce it (Rüdiger et al. 1998). The baroclinic flow has the opposite direction and builds up the differential rotation along the polar axis. The observed superrotation beneath the equator, however, cannot be produced by meridional circulation but should be a direct consequence of the Λ effect.

7. Conclusions

Stellar differential rotation is driven by Reynolds stress, by the centrifugal-induced (“Biermann-Kippenhahn”) flow, and (for stratified density) by the baroclinic flow. In our model we have found that both meridional circulations have opposite directions. The baroclinic flow becomes important for more rapid rotation as the convective heat flux deviates from the radial direction. This deviation can be interpreted as a tilt of the heat flux vector towards the rotation axis and causes the poles to be slightly warmer than the regions of lower latitudes. The baroclinic flow is responsible for the strong negative gradient dΩ/dz along the rotation axis, and the Reynolds stress produces the typical positive dΩ/dr along the equatorial mid plane, which are results that are both known from helioseismology. As a result of both impacts, the isolines of Ω known for the Sun are almost radial at mid-latitudes.

Studying real stellar models involves varying a range of parameters such as mass, radius, effective temperature, and the depth of the convection zone, which, however, are interdependent. Kitchatinov & Olemskoy (2011) studied differential rotation along the lower ZAMS for fixed rotation rate and found that the surface shear is a function of the effective temperature alone. Using some artificial models, we find that for the same effective temperature the depth of the convection zone has a significant impact on stellar rotation. Shallow convection zones produce a stronger surface shear.

The presented mean-field theory for G stars naturally explains the rotation laws of the Sun, of the MOST stars ϵ Eri and κ¹ Cet, and of rapidly rotating stars such as R58 and LQ Lup. The discrepancy between these stars and the much stronger differential rotation observed for HD 171488 – if real – is hard to explain for stars of similar age and spectral type. If, however, HD 171488 had a shallower convection zone for some reason, its strong surface differential rotation would be easily produced.

During the pre-MS evolution, the convection zone retreats from the central region and the originally fully convective star forms a radiative zone around the core. This change in the depth of the convection zone should be reflected by the surface differential rotation.

Acknowledgments

This work has been funded by the Deutsche Forschungsgemeinschaft (RU 488/21). L.L.K. is thankful to the support by the Russian Foundation for Basic Research (projects 10-02-00148, 10-02-00391).

References

All Tables

All Figures

Fig. 1 The solar differential rotation and meridional flow as computed with the new scheme. Top left: surface rotation rate versus colatitude. Top right: rotation rate as function of radius at 0, 15, 30, 45, 60, 75, and 90° latitude, respectively, from top to bottom. Bottom left: streamlines of the meridional flow. The flow is counter-clockwise, i.e. directed towards the pole at the surface and towards the equator at the bottom. Bottom right: flow speed as a function of latitude at the top (solid blue) and bottom (dashed red) of the convection zone. Positive values indicate a flow away from the north pole.

In the text

	Fig. 2 Temperature deviation as function of the latitude at the bottom (left) and top (right) of the solar convection zone. The polar axis is always warmer than the equator. The differences, however, are small.
In the text

	Fig. 3 Rotation of a hypothetical rapidly rotating Sun with P = 1.33 d.
In the text

	Fig. 4 Solar rotation law and meridional flow without Λ effect. Top: the rotation law. Note the steep negative gradient of the angular velocity at the poles and the almost rigid rotation beneath the equator. Bottom: the meridional flow is anti-solar: it is positive (equatorward) at the top of the convection zone and negative (poleward) at the bottom of the convection zone.
In the text

	Fig. 5 Solar rotation law without baroclinic terms. The structure disappears now at the polar axis but appears at the equatorial region. The equatorial plane shows superrotation but the polar axis rotates almost rigidly.
In the text

	Fig. 6 Rotation pattern of CoRoT-2a, a young Sun rotating with P = 4.5 days.
In the text

	Fig. 7 Internal rotation of young G dwarfs with rotation period of 1.33 days. Top: deep convection zone, δΩ = 0.11 rad/day at the surface. Bottom: shallow convection zone. δΩ = 0.5 rad/day at the surface.
In the text

	Fig. 8 Differential rotation of the G dwarf with the shallow convection zone and a rotation period of 0.33 days. Left: isocontours of the angular velocity. Right: the rotation profiles at various latitudes.
In the text

	Fig. 9 The drivers of angular momentum versus fractional radius at 45° latitude for the Sun (left, slow rotation) and the deep convection zone G dwarf (right, rapid rotation). Dashed line: baroclinic term. Dash-dotted line: centrifugal term. Solid line: both.
In the text